# Maple V.3, no intended for a package

with(linalg):
# example from class
A := matrix([ [1, 2, 1, 3], [3, -1, -3, -1], [2, 3, 1, 4] ]):
# example from class
B := matrix([ [1, 1, 1, 1, 1, 1], [-1, -1, 0, 0, 1, -1],
  [-2, -2, 0, 0, 3, 1], [0,0,1,1,3,-1], [1, 1, 2, 2,4,1] ]):
# example p. 26, (l)
C := matrix([ [1, -3, 1, 1], [2, 1, -1, 2],
              [1, 4, -2, 1], [5, -8, 2, 5] ]):


ref := proc(A)
# returns the row echelon form of the matrix A
  local m,n, i,j,k, AA;
  # check the type of the argument to be sure
  if not(type(A, matrix))
     then ERROR(`arg not a matrix`, A);
  fi;
  AA := copy(A); # so that argument remains unchanged
  m := rowdim(AA);
  if nargs > 1 # optional second argument
     then n := args[2];
     else n := coldim(AA);
  fi;
  i := 1; # current row
  j := 1; # current column
  while i<m and j<=n do
      # 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
         AA := swaprow(AA, i, k);
         # now AA[i,j] <> 0, so eliminate rest of column
         for k from i+1 to m do
             # set AA[k,j] to 0
             AA := addrow(AA, i, k, -AA[k,j]/AA[i,j]);
         od;
         i := i+1; j := j+1;
         if printlevel >= 2 then print(AA); fi;
      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;
  RETURN(evalm(AA));
end:

# this code creates the online documentation for ref
`help/text/ref` := TEXT(
`FUNCTION: ref - compute the row echelon form of a matrix`,
`   `,
`CALLING SEQUENCES: ref(A);  or`,
`                   ref(A,m);`,
`   `,
`PARAMETERS: A - a matrix`,
`            m - non-negative integer`,
`   `,
`SYNOPSIS:   `,
`- The call ref(A); returns the row echelon form of the`,
`  matrix A, where elimination stops after the last row/column,`,
`  whichever comes first`,
`- The call ref(A,m) returns the row echelon form of the`,
`  matrix A, where elimination stops in column m, unless it`,
`  reaches the bottom row first`,
`   `,
`EXAMPLE:   `,
`   `,
`> B := matrix([ [1, 1, 1, 1, 1, 1], [-1, -1, 0, 0, 1, -1],`,
`     [-2, -2, 0, 0, 3, 1], [0,0,1,1,3,-1], [1, 1, 2, 2,4,1] ]);`,
`   `,
`                               [  1   1  1  1  1   1 ]`,
`                               [                     ]`,
`                               [ -1  -1  0  0  1  -1 ]`,
`                               [                     ]`,
`                          B := [ -2  -2  0  0  3   1 ]`,
`                               [                     ]`,
`                               [  0   0  1  1  3  -1 ]`,
`                               [                     ]`,
`                               [  1   1  2  2  4   1 ]`,
`   `,
`> ref(B, 5);`,
`   `,
`                             [ 1  1  1  1  1   1 ]`,
`                             [                   ]`,
`                             [ 0  0  1  1  2   0 ]`,
`                             [                   ]`,
`                             [ 0  0  0  0  1   3 ]`,
`                             [                   ]`,
`                             [ 0  0  0  0  0  -4 ]`,
`                             [                   ]`,
`                             [ 0  0  0  0  0  -3 ]`
):