Reduced row echelon form

You are encouraged to solve this task according to the task description, using any language you may know.


Task

Show how to compute the reduced row echelon form
(a.k.a. row canonical form) of a matrix.

The matrix can be stored in any datatype that is convenient
(for most languages, this will probably be a two-dimensional array).

Built-in functions or this pseudocode (from Wikipedia) may be used:

function ToReducedRowEchelonForm(Matrix M) is
    lead := 0
    rowCount := the number of rows in M
    columnCount := the number of columns in M
    for 0 ≤ r < rowCount do
        if columnCount ≤ lead then
            stop
        end if
        i = r
        while M[i, lead] = 0 do
            i = i + 1
            if rowCount = i then
                i = r
                lead = lead + 1
                if columnCount = lead then
                    stop
                end if
            end if
        end while
        Swap rows i and r
        If M[r, lead] is not 0 divide row r by M[r, lead]
        for 0 ≤ i < rowCount do
            if i ≠ r do
                Subtract M[i, lead] multiplied by row r from row i
            end if
        end for
        lead = lead + 1
    end for
end function

For testing purposes, the RREF of this matrix:

 1    2   -1   -4
 2    3   -1   -11
-2    0   -3    22

is:

 1    0    0   -8
 0    1    0    1
 0    0    1   -2



#include once "matmult.bas"
 
sub rowswap( byval M as Matrix, i as uinteger, j as uinteger )
    dim as integer k
    for k = 0 to ubound(M.m, 2)
        swap M.m(j, k), M.m(i, k)
    next k
end sub
 
function rowech(byval M as Matrix) as Matrix
    dim as uinteger lead = 0, rowCount = 1+ubound(M.m, 1), colCount = 1+ubound(M.m, 2)
    dim as uinteger r, i, j
    dim as double K
    for r = 0 to rowCount-1
        if lead >= colCount then exit for
        i = r
        while M.m(i, lead) = 0
            i += 1
            if i = rowCount then
                i = r
                lead += 1
                if lead = colCount then exit for
            endif
        wend
        rowswap M, r, i
        K = M.m(r,lead)
        if K <> 0 then 
            for j = 0 to colCount-1
                M.m(r,j) /= K
            next j
        endif
        for i = 0 to rowCount-1
            if i <> r then
                K = M.m(i, lead)
                for j = 0 to colCount-1
                    M.m(i,j) -= M.m(r,j) * K
                next j
            endif
        next i
        lead += 1
    next r
    return M
end function
 
 
dim as Matrix M = Matrix (3, 4)
dim as Matrix N
 
M.m(0,0) = 1 : M.m(0,1) = 2 : M.m(0,2) = -1 : M.M(0,3) = -4
M.m(1,0) = 2 : M.m(1,1) = 3 : M.m(1,2) = -1 : M.m(1,3) = -11
M.m(2,0) = -2: M.m(2,1) = 0 : M.m(2,2) = -3 : M.m(2,3) = 22
 
dim as integer i, j
 
N = rowech(M)
for i=0 to 2
    for j = 0 to 3
        print N.m(i, j),
    next j
    print
next i