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