# procedures myinverse, mydet;  Maple V.4, packable

refpkg['myinverse'] := proc(A::matrix)
  # compute inverse matrix by Gauss-Jordan elimination (p. 59)
  local n, i,k, AA, Ainv, alpha;
  option  `Copyright (c) 1997 by Erich Kaltofen`;

  AA := copy(A); # so that argument remains unchanged
  n := linalg[rowdim](AA);
  if n <> linalg[coldim](AA)
     then ERROR(`arg not a square matrix`, A);
  fi;
  Ainv := array(identity, 1..n, 1..n);
  for i from 1 to n do
      # find a pivot element in column i
      for k from i to n do
         if AA[k,i] <> 0 then break; # out of for loop
         fi;
      od; # end pivot search
      if k <= n then
         # found a pivot element
         if i<> k then
            # place pivot element in i-th row
            if printlevel > 2 then print(`Elem step: exchange rows `, i, k); fi;
            AA := linalg[swaprow](AA, i, k);
            Ainv := linalg[swaprow](Ainv, i, k);
            if printlevel > 2 and i<>k then print(`AA=`,AA, `T=`,Ainv); fi;
         fi;
         # now AA[i,i] <> 0. First normalize AA[i,i] to 1
         alpha := 1/AA[i,i];
         if printlevel > 2 then print(`Elem step: mult row `, i, ` by `, alpha); fi;
         AA := linalg[mulrow](AA, i, alpha);
         Ainv := linalg[mulrow](Ainv, i, alpha);
         if printlevel > 2 then print(`AA=`,AA, `T=`,Ainv); fi;
         # Second, eliminate all other entries in col i
         for k from 1 to n do
             if k<>i then
                # set AA[k,i] to 0
                alpha := -AA[k,i];
                if printlevel > 2 then print(`Elem step: add `, alpha, ` times row `,
                                              i, ` to row `, k); fi;
                AA := linalg[addrow](AA, i, k, alpha);
                Ainv := linalg[addrow](Ainv, i, k, alpha);
                if printlevel > 2 then print(`AA=`,AA, `T=`,Ainv); fi;
             fi;
         od;
         if printlevel > 2 then print(`Finished col `, i); fi;
      else
         # all elements in i-th col below i-1-st row are 0
         ERROR(`arg not a non-singular matrix`, A);
      fi;
  od;
  RETURN(evalm(Ainv));
end:

refpkg['mydet'] := proc(A::matrix)
  # compute the determinant of A by Gaussian elimination (p. 94)
  local n, i,k, AA, det, alpha;
  option  `Copyright (c) 1997 by Erich Kaltofen`;

  AA := copy(A); # so that argument remains unchanged
  n := linalg[rowdim](AA);
  if n <> linalg[coldim](AA)
     then ERROR(`arg not a square matrix`, A);
  fi;
  det := 1;
  for i from 1 to n do
      # find a pivot element in column i
      for k from i to n do
         if AA[k,i] <> 0 then break; # out of for loop
         fi;
      od; # end pivot search
      if k <= n then
         # found a pivot element
         # place pivot element in i-th row
         if i<>k then
            if printlevel > 2 then print(`Elem step: exchange rows `, i, k); fi;
            AA := linalg[swaprow](AA, i, k);
            det := -det;
         fi;
         # now AA[i,i] <> 0.
         # Eliminate all entries in col i and rows i+1,...,n
         for k from i+1 to n do
             # set AA[k,i] to 0
             alpha := simplify(-AA[k,i]/AA[i,i]);
             if printlevel > 2 then print(`Elem step: add `, alpha, ` times row `,
                                           i, ` to row `, k); fi;
             AA := linalg[addrow](AA, i, k, alpha);
         od;
         det := simplify(det*AA[i,i]);
         if printlevel > 2 then print(AA, det); fi;
      else
         # all elements in i-th col below i-1-st row are 0
         RETURN(0);
      fi;
  od;
  RETURN(det);
end: