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 |