Task

Invert matrix   A   using Gauss-Jordan method.

A   being an   n × n   matrix.




with javascript_semantics
function ToReducedRowEchelonForm(sequence M)
integer lead = 1,
        rowCount = length(M),
        columnCount = length(M[1]),
        i
    for r=1 to rowCount do
        if lead>=columnCount then exit end if
        i = r
        while M[i][lead]=0 do
            i += 1
            if i=rowCount then
                i = r
                lead += 1
                if lead=columnCount then exit end if
            end if
        end while
        object mr = sq_div(M[i],M[i][lead])
        M[i] = M[r]
        M[r] = mr
        for j=1 to rowCount do
            if j!=r then
                M[j] = sq_sub(M[j],sq_mul(M[j][lead],M[r]))
            end if
        end for
        lead += 1
    end for
    return M
end function
--</end of copy>
 
function inverse(sequence mat)
    integer len = length(mat)
    sequence aug = repeat(repeat(0,2*len),len)
    for i=1 to len do
        if length(mat[i])!=len then ?9/0 end if -- "Not a square matrix"
        for j=1 to len do
            aug[i][j] = mat[i][j]
        end for
        -- augment by identity matrix to right
        aug[i][i + len] = 1
    end for
    aug = ToReducedRowEchelonForm(aug)
    sequence inv = repeat(repeat(0,len),len)
    -- remove identity matrix to left
    for i=1 to len do
        for j=len+1 to 2*len do
            inv[i][j-len] = aug[i][j]
        end for
    end for
    return inv
end function
 
constant test = {{ 2, -1,  0},
                 {-1,  2, -1},
                 { 0, -1,  2}}
pp(inverse(test),{pp_Nest,1})
