macro( GaussJordan = `Matrix/GaussJordan` );
macro( GE = readlib(`Matrix/GaussianElimination`, ``.gauss.`/MX/GE.m`) );

#
#--> GaussianJordan(M,A,'rank','det')
#
# Complete elimination of 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
#

GaussJordan := proc(M,B,rank,det)
local i,j,c,m,n,r,A;

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

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

    # Forward-elimination: reduce A to unit normal upper triangular form
    A := M[ArrayCoeffs](GE(args));

    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 next fi;

        for i to r-1 do
	    if A[i,c] = R[0] then next fi;
    	    for j from c+1 to m do
		A[i,j] := R[`-`](A[i,j],R[`*`](A[i,c],A[r,j])) od;
	    A[i,c] := R[0]
        od;
        r := r + 1		# go to next row

    od;			# go to next column

    M[Matrix](A,n,m)

end:


save `GJ.m`;
quit
