QR decomposition (Eiffel)

From LiteratePrograms

This class makes use of Real matrix (Eiffel).

<<qr_decomposition.e>>=
class QR_DECOMPOSITION[X -> REAL_MATRIX] inherit MATRIX_DECOMPOSITION[X]
        rename
            matrix1 as q,
            matrix2 as r
        end
creation
    make,
    make_hessenberg
feature {ANY} -- creation
    make(a : X) is
        require
            a.rows = a.columns
        do
            create q.make(a.rows, a.columns)
            create r.make(a.rows, a.columns)
            decompose(a.deep_twin)
        end
    make_hessenberg(a : X) is
        require
            a.rows = a.columns
            is_hessenberg(a)
        do
            create q.make(a.rows, a.columns)
            create r.make(a.rows, a.columns)
            decompose_hessenberg(a.deep_twin)
        end            
feature {ANY}
    decompose(a : X) is
        local
            i : INTEGER
            alpha : REAL
            v, v_t : X
            q_i, q_i_submatrix : X
        do
            from
                i := 1
                q := a.one
                r := a.deep_twin
            until
                i >= a.rows
            loop
                v := r.submatrix(i, r.rows, i, i)
                alpha := v.vector_norm
                if alpha /= 0.0 then
                    v.set_element(1, 1, v.element(1, 1) - alpha)
                    v := v #* (v.vector_norm^(-1))
                    v_t := v.deep_twin
                    v_t.transpose
                    create q_i.make_identity(a.rows)
                    create q_i_submatrix.make_identity(v.rows)
                    q_i_submatrix := q_i_submatrix - (v * v_t) #* 2
                    q_i.set_submatrix(i, a.rows, i, a.columns, q_i_submatrix)
                    r := q_i * r
                    q := q * q_i
                end
                i := i + 1
            end
        end -- decompose
    decompose_hessenberg(a : X) is
        require
            a.rows = a.columns
            is_hessenberg(a)
        local
            i : INTEGER
            theta : REAL
            omega : X
        do
            from
                i := 1
                q := a.one
                r := a.deep_twin
            until
                i >= a.rows
            loop
                omega := a.one
                theta := (r.element(i + 1, i)/r.element(i,i)).atan
                omega.set_element(i, i, theta.cos)
                omega.set_element(i + 1, i + 1, theta.cos)
                omega.set_element(i, i + 1, theta.sin)
                omega.set_element(i + 1, i, -theta.sin)
                q := omega * q
                r := omega * r
                i := i + 1
            end
            q.transpose
        end -- decompose_hessenberg
feature {NONE}
    is_hessenberg(a : X) : BOOLEAN is
        local
            i, j : INTEGER
        do
            Result := True
            from
                i := 3
            until
                i > a.rows
            loop
                from
                    j := 1
                until
                    j >= a.columns - 1
                loop
                    -- allow a little leeway
                    if a.element(i,j).abs > 0.000000000000001 then
                        Result := False
                    end
                    j := j + 1
                end
                i := i + 1
            end
        end -- is_hessenberg
end -- class QR_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