macro( GaussianElimination = `Matrix/GaussianElimination` );

#
#--> GaussianElimination(M,A,'rank','det')
#
# Ordinary Gaussian elimination on the Matrix A in M where M is a Matrix
# domain over a field.
# Output the reduced matrix and optionally the rank and determinant
#
# Author: MBM 1990
#

GaussianElimination := proc(M,B,rank,det)
local d,i,j,k,c,m,n,r,s,t,A;

    R := M[CoefficientRing];
    if not hasCategory(R,Field) then
        ERROR(`Coefficient ring must be a field`) fi;

    if not M[Type](B) then
	ERROR(`1st argument must be a matrix`) fi;

    n := M[Rows](B);
    m := M[Cols](B);

    if nargs > 2 and not type(rank,name) then
	ERROR(`1st optional argument must be a name`) fi;

    if nargs > 3 and not type(det,name) then
	ERROR(`2nd optional argument must be a name`) fi;

    if nargs > 3 and n <> m then
	ERROR(`cannot compute determinant of a non-square matrix`) fi;

    A := M[ArrayCoeffs](B);

    d := R[1];
    for i to n do
	for j to m while A[i,j] = R[0] do od;
	if j > m or A[i,j] = R[1] then next fi;
	d := R[`*`](d,A[i,j]);
	s := R[Inv](A[i,j]);
	for k from j+1 to m do A[i,k] := R[`*`](s,A[i,k]) od;
	A[i,j] := R[1];
    od;
	
    r := 1;

    for c to m while r <= n do

        for i from r to n while A[i,c] = R[0] do od;
        if i > n then d := R[0]; next fi;

        if i <> r then
	    d := R[`-`](d);
	    # interchange row i with row r
	    for j from c to m do t := A[i,j]; A[i,j] := A[r,j]; A[r,j] := t od
        fi;

        for i from r+1 to n do
	    if A[i,c] = 0 then next fi;
	    A[i,c] := 0;
    	    for j from c+1 to m do A[i,j] := R[`-`](A[i,j],A[r,j]) od;
	    for j from c+1 to m while A[i,j] = R[0] do od;
	    if j > m then next fi;
	    d := R[`*`](A[i,j],d);
	    s := R[Inv](A[i,j]);
	    for k from j+1 to m do A[i,k] := R[`*`](s,A[i,k]) od;
	    A[i,j] := 1;
        od;

        r := r + 1		# go to next row

    od;			# go to next column

    if nargs>2 then rank := r-1 fi;
    if nargs>3 then det := d fi;
    M[Matrix](A,n,m)

end:

save `GE.m`;
quit
