LU decomposition (Eiffel)

From LiteratePrograms

This class makes use of Real matrix (Eiffel).

<<lu_decomposition.e>>=
class LU_DECOMPOSITION[X -> REAL_MATRIX] inherit MATRIX_DECOMPOSITION[X]
        rename
            matrix1 as lower,
            matrix2 as upper
        redefine
            recompose
        end
creation
    make
feature {ANY} -- extra properties
    permutation : X
    parity : REAL
feature {ANY} -- creation
    make(a : X) is
        require
            a.rows = a.columns
        do
            create lower.make(a.rows, a.columns)
            create upper.make_identity(a.rows)
            create permutation.make_identity(a.rows)
            parity := 1.0
            decompose(a.deep_twin)
        end
feature {ANY}
    decompose(a : X) is
        local
            i, j, n : INTEGER
            sum : REAL
            element_type : REAL
        do
            pivot(a)
            -- first column of lower triangular, first row of upper triangular
            from
                i := 1
            until
                i > a.rows
            loop
                lower.set_element(i, 1, a.element(i,1))
                upper.set_element(1, i, a.element(1,i)/lower.element(1,1))
                i := i + 1
            end
            -- caluculate the remaining portion
            from
                i := 2
            until
                i > a.rows
            loop
                from
                    j := i
                until
                    j > a.rows
                loop
                    from
                        n := 1
                        sum := element_type.zero
                    until
                        n > i - 1
                    loop
                        sum := sum + lower.element(j,n) * upper.element(n,i)
                        n := n + 1
                    end
                    lower.set_element(j, i, a.element(j,i) - sum)
                    j := j + 1
                end
                from
                    j := i + 1
                until
                    j > a.rows
                loop
                    from
                        n := 1
                        sum := element_type.zero
                    until
                        n > i - 1
                    loop
                        sum := sum + lower.element(i,n) * upper.element(n, j)
                        n := n + 1
                    end
                    upper.set_element(i, j, ((a.element(i,j) - sum)/lower.element(i,i)))
                    j := j + 1
                end
                i := i + 1
            end
        ensure then
            permutation /= Void
        end -- decompose
    recompose : X is
        do
            Result := permutation * lower * upper
        end
feature {NONE}
    pivot(a : X) is
        local
            i, j, max_index : INTEGER
            col_max : REAL
        do
            from
                j := 1
            until
                j > a.columns
            loop
                from
                    i := j
                    col_max := 0
                until
                    i > a.rows
                loop
                    if a.element(i,j).abs >= col_max then
                        col_max := a.element(i,j).abs
                        max_index := i
                    end
                    i := i + 1
                end
                if j /= max_index then
                    a.swap_rows(j, max_index)
                    permutation.swap_rows(j, max_index)
                    parity := parity * -1.0
                end
                j := j + 1
            end
        end -- pivot
end -- class LU_DECOMPOSITION

<<matrix_decomposition.e>>=
deferred class MATRIX_DECOMPOSITION[X -> MATRIX[NUMERIC]]
feature {ANY} -- storage
    matrix1 : X
    matrix2 : X
feature {ANY} -- (de|re)-composition functions
    decompose(a : X) is
        require
            a.rows = a.columns
        deferred
        ensure
            matrix1 /= Void
            matrix2 /= Void
        end -- decompose
    recompose : X is
        require
            matrix1 /= Void
            matrix2 /= Void
        do
            Result := matrix1 * matrix2
        end -- recompose
end -- class MATRIX_DECOMPOSITION

Download code