#############################################################
#
#       This program assumes that the matrix Mat is kxn, 
#       of full rank. It returns the product of elementary 
#       matrices needed to put the matrix in row-reduced echelon 
#       form. It also returns a column permutation needed 
#       to make the left-most kxk block the identity.
#
#       code is based on an old PutStandardForm program
#       which was used in GUAVA 1.7
#
# 12-12-2004, wdj
###############################################################

PutStandardFormTransform:= function(Mat, F) 
    local Elem,k, n, row, i, j, h, hp, s, zero, P, T, M;
    M:=List(Mat,ShallowCopy);
    k := Length(M);   # not the word length!
    n := Length(M[1]);
    T:=List(IdentityMat(k,F),ShallowCopy);
    zero := Zero(F);
    P := ();
    for j in [1..k] do
        if M[j][j] =zero then
            i := j+1;
            while (i <= k) and (M[i][j] = zero) do
##  this moves down the rows until it gets
##  to a first non-zero entry M[i][j]
                i := i + 1;
            od;
            if i <= k then
     Elem:=List(IdentityMat(k,F),ShallowCopy);
     Elem[i][i]:=zero;
     Elem[j][j]:=zero;
     Elem[i][j]:=One(F);
     Elem[j][i]:=One(F);
     T:=Elem*T;
                row := M[j];
                M[j] := M[i];
                M[i] := row;
##  this swaps row i and row j
            else
                h := j;
                while M[j][h] = zero do
## stops when M[j][h]<>0
                    h := h + 1;
                    if h > n then h := 1; fi;
                od;
                for i in [1..k] do
                    M[i] := Permuted(M[i],(j,h)); 
##  swaps i-th col of M by h-th col 
                od;   
                P := P*(j,h);
            fi;
        fi; #### M[j][j] =zero 
     Elem:=List(IdentityMat(k,F),ShallowCopy);
     Elem[j][j]:=1/M[j][j];
     T:=Elem*T;
        M[j] := M[j]/M[j][j];
##  normalizes j-th row of M so that its
##  leading non-zero entry is a 1
        for i in [1..k] do
            if i <> j then
                if M[i][j] <> zero then
     Elem:=List(IdentityMat(k,F),ShallowCopy);
     Elem[i][j]:=-M[i][j];
     T:=Elem*T;
                    M[i] := M[i]-M[i][j]*M[j];
##  kills off all entries of column j below M[i][j]
##  when i<>j
                fi;
            fi;
        od;
    od;
    return [P,T,M];
end;


#####################################################################