---- Gauss Elimination
-- By Mike Gage  gage@cc.rochester.edu  Apr 1995
--ееееееееееееееееееееееееееееееееееееееееееееееееееееееееееееееее
--- DEFINITIONS
-- these row commands are defined below
-- switch(i, j)  -- interchange the rows i and j
-- add(i,m,j)    -- add m times row j to row i
-- mult(i,m)     -- multiply row i by m


--backSubU(U,b)  --solves Ux=b for x where U is an upper triangular matrix (back substitution)

--columnSelect(A,i,j) -- selects columns i through j of matrix A. i must be less than or equal to j

-- define start 
start(X)=U:=X,L:=identityMatrix(count(B)),P:=identityMatrix(count(B)) 

-- define pivot(i,j,k) -- adds a multiple of row j to row i so as to make U[i,k] 0.  Make sure that U[j,k] is not zero!  Also constructs the matrix L at the same time.

pivot(i,j,k)=L[i,k]:=U[i,k]/U[j,k],add(i,-U[i,k]/U[j,k],j),L[i,k]
-- IMPORTANT NOTE:  This method of constructing the lower triangular matrix L works only if we do the  row operations in the standard order for constructing an upper triangular matrix.  If you plan to use non-standard order of row reduction, or to create zeros above the diagonal use the method in Gauss-Jordan elimination.  It is computationally slower because it doesn't use any shortcuts, but it creates matrices L and U such that B=multMatrix(L, U) for any sequence of row reduction operations.

------end of definitions
--еееееееееееееееееееееееееееееееееееееееееееееееееееееееееееееее

-- Example
B={
{2,	-1,	 0,	0,	0},
{-1, 2, -1, 0, 0},
{0,  -1, 2, -1,0},
{0,0,-1,2,5}}
start(B):--initializes the matrices 
-- B=LU factorization appears in global variables L and U.
-- use table U to see partial results
table U
      2.000     -1.000      0.000      0.000      0.000
     -1.000      2.000     -1.000      0.000      0.000
      0.000     -1.000      2.000     -1.000      0.000
      0.000      0.000     -1.000      2.000      5.000


-----------------------------------------------------------------
--- commands for doing Gauss elimination on rows
-- define permutation matrices
sigma(i,j,n)[k,l] = 1 when k=j and i=l,1 when i=k and j=l, 1 when k=l and i!=k and j!=l, 0 dim[n,n]

-- commands for row switching 
switch_internal(i,j,X)[k,l]= X[i,l] when j=k,
                          X[j,l] when i=k,
                         X[k,l] otherwise  dim[count(X),count(X[1])]
-- commands for adding and multiplying
add_internal(i,m,j)[k,l]=U[i,l]+m*U[j,l] when k=i,
                         U[k,l] otherwise dim[count(U), count(U[1])]
mult_internal(i,m)[k,l]=m*U[i,l] when i=k, U[k,l] otherwise dim[count(U),count(U[1])]

--switch rows i and j
switch(i,j)=U:=switch_internal(i,j,U),L:=transpose(switch_internal(i,j,
                                                       transpose(switch_internal(i,j,L)))),
                                      P:=switch_internal(i,j,P)  

add(i,m,j)=U:=add_internal(i,m,j)    --add m times row j to row i
mult(i,m)=U:=mult_internal(i,m)      --multiply row i by m
           
-- modify this to calculate L as well

-- back substitution command

backSubU(upper,vec)[i]=(vec[i]-sum(upper[i,j]*backSubU(upper,vec)[j],j,i+1,count(vec)))/upper[i,i]
backSubL(lower,vec)[i]=(vec[i]-sum(lower[i,j]*backSubL(lower,vec)[j],j,1,i-1))/lower[i,i]

--columns selector

columnSelect(A,i,j)[k]=A[k][i:j] when i<j,
                      A[k,i] when  i=j,0/0 otherwise  dim[count(A)]

identityMatrix(n)[i,j]=1 when i=j,0 dim[n,n]

--------------------------------------------------------------------------------
---------- matrix operations ------------
-- three multiply operations for different dimension matrices and vectors
-- it's too complicated to write one function which does all three operations in the mathpad language.  The essential problem is following the *mathematical* convention that a list with a single element can be written as the element without including the surrounding brackets. 

-- multiplying a matrix times another matrix to obtain a new matrix
multMatrix(A,B)[i,j] = sum(A[i,k]*B[k,j],k,
         1,count(B))  dim[count(A),count(B[1])]

-- multiplying a matrix by a vector the result is a vector, not a 1 by n matrix
multVector(A,vec)[i] = sum(A[i,j]*vec[j],j,1,count(vec)) dim[count(A)]

-- the dot product of two vectors produces a number, not a one by one matrix
dot(v,w) =sum(v[j]*w[j],j,1,count(w))

--outer product of vectors (also known as the tensor product)
outer(v,w)[i,j] =v[i]*w[j]

normVector(vec) = dot(vec,vec)

normMatrix(A) = trace(multMatrix(transpose(A),A))

transpose(A)[i,j] = A[j,i] dim
                        [count(A[1]),count(A)]

-- this is not the best way to invert a matrix - it is not numerically stable
invert(A) = adjoint(A)/det(A)

adjoint(A) = transpose(cofactor(A))

cofactor(A)[i,j] = (-1)^(i+j)*
                   det(submatrix(A,i,j))

-- trace (sum the diagonal)
trace(A) = sum(A[j,j],j,1,count(A))

-- remove the ith row and jth column
submatrix(A,k,l)[i,j] = A[i,j] when i<k and j<l,
            A[i+1,j] when i>=k and j<l,
            A[i+1,j+1] when i>=k and j>=l,
            A[i,  j+1] when i<k and j>=l dim
                        [count(A)-1,count(A)-1]