---- Gauss-Jordan reduction
-- By Mike Gage  gage@cc.rochester.edu  Apr 1995

--ееееееееееееееееееееееееееееееееееееееееееееееееееееееееееееееее
--- DEFINITIONS
-- start:--initializes the matrices from the matrix B 
-- PB=LU factorization appears in global variables L and U.
-- use table U to see partial results

-- these row commands are defined in the include file above
-- switch(i, j)  -- interchange the rows i and j -- also constructs permutation matrix P
-- add(i,m,j)    -- add m times row j to row i -- constructs matrices L and U
-- mult(i,m)     -- multiply row i by m        -- constructs matrices L and U


--backSubU(U,b)  --solves Ux=b for x where U is an upper triangular matrix (back substitution)
--backSubL(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=U:=B,L:=identityMatrix(count(B)),P:=identityMatrix(count(B)),
              RowOps:= 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 and the matrix RowOps.

pivot(i,j,k)=add(i,-U[i,k]/U[j,k],j),L[i,k]

------end of definitions

-- At every stage mutlMatrix(P,B)=multMatrix(L,U) and RowOps is the inverse of L
--ееееееееееееееееееееееееееееееееееееееееееееееееееееееееееееееее

-- Here is a worked out example: Use pivots to make the matrix upper triangular
B={
{2,	-1,	 0,1,	0, 0},
{-1, 2, -1, 0, 1, 0},
{0,  -1, 2, 0, 0, 1}}
start:

table multMatrix(L,U)
      2.000     -1.000      0.000      1.000      0.000      0.000
     -1.000      2.000     -1.000      0.000      1.000      0.000
      0.000     -1.000      2.000      0.000      0.000      1.000
pivot(2,1,1):-0.500; table U
      2.000     -1.000      0.000      1.000      0.000      0.000
      0.000      1.500     -1.000      0.500      1.000      0.000
      0.000     -1.000      2.000      0.000      0.000      1.000
table L
      1.000      0.000      0.000
     -0.500      1.000      0.000
      0.000      0.000      1.000
table multMatrix(L,U)
      2.000     -1.000      0.000      1.000      0.000      0.000
     -1.000      2.000     -1.000      0.000      1.000      0.000
      0.000     -1.000      2.000      0.000      0.000      1.000
pivot(3,2,2):-0.667; table U
      2.000     -1.000      0.000      1.000      0.000      0.000
      0.000      1.500     -1.000      0.500      1.000      0.000
      0.000      0.000      1.333      0.333      0.667      1.000
pivot(1,2,2):-0.667; table U
      2.000      0.000     -0.667      1.333      0.667      0.000
      0.000      1.500     -1.000      0.500      1.000      0.000
      0.000      0.000      1.333      0.333      0.667      1.000
table multMatrix(L,U)
      2.000     -1.000      0.000      1.000      0.000      0.000
     -1.000      2.000     -1.000      0.000      1.000      0.000
      0.000     -1.000      2.000      0.000      0.000      1.000
pivot(1,3,3):-0.500;pivot(2,3,3):-0.750; table U
      2.000      0.000      0.000      1.500      1.000      0.500
      0.000      1.500      0.000      0.750      1.500      0.750
      0.000      0.000      1.333      0.333      0.667      1.000
table multMatrix(L,U)
      2.000     -1.000      0.000      1.000      0.000      0.000
     -1.000      2.000     -1.000      0.000      1.000      0.000
      0.000     -1.000      2.000      0.000      0.000      1.000
table L
      1.000     -0.667      0.000
     -0.500      1.333     -0.750
      0.000     -0.667      1.500
mult(1,1/2):; mult(2,1/1.5):; mult(3,1/1.33333333):
table U
      1.000      0.000      0.000      0.750      0.500      0.250
      0.000      1.000      0.000      0.500      1.000      0.500
      0.000      0.000      1.000      0.250      0.500      0.750
table L
      2.000     -1.000      0.000
     -1.000      2.000     -1.000
      0.000     -1.000      2.000
table multMatrix(L,U)
      2.000     -1.000      0.000      1.000      0.000      0.000
     -1.000      2.000     -1.000      0.000      1.000      0.000
      0.000     -1.000      2.000      0.000      0.000      1.000
table RowOps
      0.750      0.500      0.250
      0.500      1.000      0.500
      0.250      0.500      0.750

-------------------------------------------------------------------------
--- commands for doing Gauss elimination on rows
--- this version is slower, but more flexible.

-- commands for row switching 
switch_row(X,i,j)[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_row(X,i,m,j)[k,l]=X[i,l]+m*X[j,l] when k=i,
                         X[k,l] otherwise dim[count(X), count(X[1])]
add_col(X,i,m,j)[l,k]=X[l,i]+m*X[l,j] when k=i,
                         X[l,k] otherwise dim[count(X), count(X[1])]
mult_row(X,i,m)[k,l]=m*X[i,l] when i=k, X[k,l] otherwise dim[count(X),count(X[1])]
mult_col(X,i,m)[l,k]=m*X[l,i] when i=k, X[l,k] otherwise dim[count(X),count(X[1])]

--switch rows i and j
switch(i,j)=U:=switch_row(U,i,j),
            RowOps:= switch_row(L,i,j),
            L:=transpose(switch_row(i,j,transpose(switch_row(i,j,L)))),
                                      P:=switch_rowl(i,j,P)  

add(i,m,j)=U:=add_row(U,i,m,j),    --add m times row j to row i and store in U
           L:= add_col(L,j,-m,i),    -- add -m * times col i to row j and store in L
											RowOps:=add_row(RowOps,i,m,j)  -- record row operations
mult(i,m)=U:=mult_row(U,i,m),      --multiply row i by m and store result in U
          L:=mult_col(L,i,1/m),				-- multiply col i by 1/m and store result in L
										RowOps:=mult_row(RowOps,i,m)   -- record row operations   
           
-- At all times B = multMatrix(L,U) should hold so L and U represent a factorization, even though L is in general no longer lower triangular.
-- The matrix RowOps records the row operations performed.


-- the commands for substitution and back substitution:
-- solves Ux=vec for the vector x where U is an upper triangular matrix
backSubU(upper,vec)[i]=(vec[i]-sum(upper[i,j]*backSubU(upper,vec)[j],j,i+1,count(vec)))/upper[i,i]
-- solves Lx=vec for the vector x where L is a lower triangular matrix
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]