gausselim := proc(A::matrix)
local B,n,m,i,j,c,t,r;
n := rowdim(A); m := coldim(A);
B := map(Normalizer, A); r := 1;
for c to m while r <= n do
    # Search for a provably non-zero pivot element
    i := 0;
    for j from r to n do
        if not Testzero(B[j,c]) then
            if i = 0 then i := j;
            elif length(B[j,c]) < length(B[i,c]) then 
                i := j
            end if; # otherwise B[i,c] is the best
        end if;
    end do;
    if i <> 0 then
        # interchange row i with row r is necessary
        if i <> r then
          for j from c to m do 
            t := B[i,j]; B[i,j] := B[r,j]; B[r,j] := t
          end do
        end if;
        for i from r+1 to n do
          if B[i,c] <> 0 then 
              t := Normalizer(B[i,c]/B[r,c]);
              for j from c+1 to m do 
                B[i,j] := Normalizer(B[i,j]-t*B[r,j]) 
              end do;
              B[i,c] := 0
          end if;
        end do;
        r := r + 1          # go to next row
    end if
end do;              # go to next column
end proc: