# procedure XRef,  Maple 6, packable
# application of Gaussian elimination: factorization

RefPkg['XRef'] := proc(A::Matrix)
# returns the row echelon form of the matrix A
# an optional second argument stops reduction at the given column
# an optional third argument, which must be a quoted variable,
#    assigns the transforming matrix to that variable
# an optional fourth argument, which must be a quoted variable,
#    assigns the inverse transforming matrix to that variable
# an optional fifth argument, which must be a quoted variable,
#    assigns a (column) permutation matrix P to that variable such that
#
#          A = (trafo)^-1 &* modified-ref(A) &* P,
#
#    where submatrix(modified-ref(A),,1..n) is a row echelon form provided
#    all zero rows are removed.  For non-singular matrices with
#    n = coldim(A) we have
#
#          A = L &* U &* P.
#
#    P enforces trafo to be lower triangular, while retaining
#    back-substitution for solving: modified-ref(A) &* x = b
  local m,n, i,j,k, AA,L,Linv, alpha, P;
  option  `Copyright (c) 1997;2002 by Erich Kaltofen`;
  description `optional arguments m, T, L, P`;

  AA := copy(A); # so that argument remains unchanged
  m := LinearAlgebra[RowDimension](AA);
  if nargs > 1 # number of columns considered for elim
     then n := args[2];
     else n := LinearAlgebra[ColumnDimension](AA);
  fi;
  if nargs > 2 # name of transforming matrix
     then L := LinearAlgebra[IdentityMatrix](m,compact=false);
  fi;
  if nargs > 3 # name of inv trafo matrix
     then Linv := LinearAlgebra[IdentityMatrix](m,compact=false);
  fi;
  if nargs > 4 # name of permutation matrix
     then P := LinearAlgebra[IdentityMatrix](LinearAlgebra[ColumnDimension](A),
                                             compact=false);
  fi;
  i := 1; # current row
  j := 1; # current column
  while i<m and j<=n do
      # permute columns if requested
      if nargs > 4 and AA[i,j] = 0
         then while i<m do
                  for k from j to n do
                      if AA[i,k] <> 0 then break; fi;
                  od;
                  if k <= n
                     then LinearAlgebra[ColumnOperation](AA,[j,k],inplace=true);
                          LinearAlgebra[RowOperation](P,[j,k], inplace=true);
                          userinfo(3, {'RefPkg', 'XRef'}, `NoName`,
                             print(`i=`,i,`AA=`,AA,`P=`,P));
                          break;
                     else i := i+1; # zero row, search next row
                  fi;
              od;
      fi;
      # find a pivot element in column j
      for k from i to m do
         if AA[k,j] <> 0 then break; # out of for loop
         fi;
      od; # end pivot search
      if k <= m then
         # found a pivot element
         # place pivot element in i-th row
         LinearAlgebra[RowOperation](AA, [i, k], inplace=true);
         if nargs > 2 then LinearAlgebra[RowOperation](L, [i, k],inplace=true); fi;
         if nargs > 3 then LinearAlgebra[ColumnOperation](Linv, [i, k], inplace=true); fi;
         # now AA[i,j] <> 0, so eliminate rest of column
         for k from i+1 to m do
             # set AA[k,j] to 0
             alpha := -AA[k,j]/AA[i,j];
             LinearAlgebra[RowOperation](AA, [k,i],alpha, inplace=true);
             if nargs > 2 then LinearAlgebra[RowOperation](L, [k,i],alpha,inplace=true); fi;
             if nargs > 3 then LinearAlgebra[ColumnOperation](Linv, [i,k], -alpha,inplace=true); fi;
         od;
         userinfo(3, {'RefPkg', 'XRef'}, `NoName`, print(`i=`,i,`j=`,j,AA));
         if nargs > 2 then userinfo(3, {'RefPkg', 'XRef'}, `NoName`, print(`L=`,L)); fi;
         if nargs > 3 then userinfo(3, {'RefPkg', 'XRef'}, `NoName`, print(`L^(-1)=`,Linv)); fi;
         i := i+1; j := j+1;
      else
         # all elements in j-th col below i-1-st row are 0
         j := j+1; # skip to next col, but stay in row
      fi;
  od;
  if nargs > 2 then assign(args[3], L); fi;
  if nargs > 3 then assign(args[4], Linv); fi;
  if nargs > 4 then assign(args[5], P); fi;
  RETURN(AA);
end: