########################
# Question 5           #
# Gaussian Elimination #
########################
maxInCol := proc(M,c,r)
  local i, locMax:
  locMax := 1:
  for i from r to RowDimension(M) do
    if abs(M[i,c]) >= abs(M[locMax,c]) then
      locMax := i:
   fi:
  od:
  RETURN(locMax):
end proc:
GaussElimPP := proc(A) #input augmented matrix A|b of coefficients of eqs and solutions, such that Ax=b
local n, m, x, p, k, lk, i, s, f, j:  #local variables (n, m dimensions, lk multipluer, s sum, others loop vars)
  #Dimensions of A
  n := RowDimension(A):
  m := ColumnDimension(A):

  #set up a list x for back substitution
  x := []:
  for i from 1 to n do
    x := [op(x),0]:
  od:

  for p from 1 to (m-2) do
    #find the largest magnitude in column, swap that row with row p
    A[p],A[maxInCol(A,p,p)] := A[maxInCol(A,p,p)],A[p]:

    #row ops
    for k from p+1 to n do
      lk := A[k,p]/A[p,p]:   #multiplier
      for i from p to m do
        A[k,i] := A[k,i] - lk*A[p,i]:  #row operation
      od: #for i p to m
    od: #for k p+1 to n
  od: #for p 1 to (m-2)

  #check for solveability
  s := 0:
  for f from 1 to m-1 do
    s := s + A[n,f]:
  od:
  if s = 0 then
    if A[n,m] = 0 then
      printf("There are no unique solutions"):
      RETURN(0):
    else
      printf("There are no solutions"):
      RETURN(0):
    fi:
  else
    #back substitution
    for j from n to 1 by -1 do
      x[j] := (A[j,m] - add(A[j,h]*x[h],h=j+1..n))/A[j,j]:
    od:
    RETURN(x):
  fi:
end proc:
GaussElimPP(Matrix([[1,3,2,3],[2,-1,-3,-8],[5,-2,1,9]]));
GaussElimPP(Matrix([[1,1,1,2],[1,-1.1,1,-1],[1,1,1.1,3]]));