Real matrix (Eiffel)
From LiteratePrograms
This class makes use of Matrix (Eiffel), LU Decomposition (Eiffel) and QR Decomposition (Eiffel).
<<real_matrix.e>>= class REAL_MATRIX inherit MATRIX[REAL] redefine determinant end; NUMERIC creation make, make_identity, from_array feature {ANY} -- transforms invert is require rows = columns determinant /= element_type.zero local i : INTEGER lu : LU_DECOMPOSITION[like Current] vector, col : like Current do create lu.make(Current) create vector.make(rows, 1) from i := 1 until i > rows loop vector.set_element(i, 1, 1.0) col := backsub(lu, vector) Current.set_submatrix(1, rows, i, i, col) vector.set_element(i, 1, 0.0) -- for next pass i := i + 1 end end --invert feature {ANY} -- arithmetic infix "/" (other : like Current) : like Current is require else other.rows = columns other.rows = other.columns other.determinant /= element_type.zero do Result := Current * other.inverse end -- divide feature {ANY} -- operators divisible(other : like Current) : BOOLEAN is do if other.rows = other.columns and other.determinant /= element_type.zero then Result := True else Result := False end end -- divisible determinant : REAL is local i : INTEGER decomposition : LU_DECOMPOSITION[like Current] det_l, det_u : REAL do create decomposition.make(Current) from i := 1 det_l := element_type.one det_u := element_type.one until i > rows loop det_l := det_l * (decomposition.lower).element(i,i) det_u := det_u * (decomposition.upper).element(i,i) i := i + 1 end Result := det_l * det_u * decomposition.parity end -- determinant inverse : like Current is require rows = columns determinant /= element_type.zero local i : INTEGER lu : LU_DECOMPOSITION[like Current] vector, col : like Current do create lu.make(Current) create Result.make(rows, columns) create vector.make(rows, 1) from i := 1 until i > rows loop vector.set_element(i, 1, 1.0) col := backsub(lu, vector) Result.set_submatrix(1, rows, i, i, col) vector.set_element(i, 1, 0.0) -- for next pass i := i + 1 end end -- inverse eigenvalues : ARRAY[REAL] is require rows = columns determinant /= element_type.zero local i, iteration_count : INTEGER b, c : REAL a : like Current shifts : ARRAY[REAL] qr : QR_DECOMPOSITION[REAL_MATRIX] eigenvalues_found, failed : BOOLEAN do create Result.with_capacity(0,1) -- handle 2 x 2 case up front if rows = 2 then b := -(element(1,1) + element(2,2)) c := determinant if b^2 - 4*c < 0 then -- complex roots: fail Result := <<1.0/0.0, 1.0/0.0>> end Result.add_last((-b + (b^2 - 4*c).sqrt)/2.0) Result.add_last((-b - (b^2 - 4*c).sqrt)/2.0) else from a := hessenberg_form iteration_count := 0 failed := False eigenvalues_found := False create qr.make_hessenberg(a) until eigenvalues_found or failed loop shifts := a.submatrix(a.rows - 1, a.rows, a.columns - 1, a.columns).eigenvalues if shifts.item(1).is_infinity or shifts.item(2).is_infinity then failed := True else from i := 1 until i > 2 loop a := a - (a.one #* shifts.item(i)) qr.decompose_hessenberg(a) a := qr.r * qr.q + (a.one #* shifts.item(i)) i := i + 1 end if a.element(a.rows, a.columns - 1).abs < 0.0000000000000001 then Result.add_last(a.element(a.rows, a.columns)) if a.rows > 3 then a := a.submatrix(1, a.rows - 1, 1, a.columns - 1) else Result.append_collection(a.submatrix(1,2,1,2).eigenvalues) eigenvalues_found := True end end if a.element(a.rows - 1, a.columns - 2).abs < 0.0000000000000001 then Result.append_collection(a.submatrix(a.rows - 1, a.rows, a.columns - 1, a.columns).eigenvalues) if a.rows > 4 then a := a.submatrix(1, a.rows - 2, 1, a.columns - 2) elseif a.rows = 3 then Result.add_last(a.element(1,1)) eigenvalues_found := True else Result.append_collection(a.submatrix(1,2,1,2).eigenvalues) eigenvalues_found := True end end iteration_count := iteration_count + 1 io.put_string("eigenvalue iteration " + iteration_count.out + "%N") if iteration_count > 30 * rows then failed := True end end -- if minor has real evals end -- loop if failed then create Result.make(1, rows) Result.set_all_with(1.0/0.0) end end -- if 2 x 2 ensure Result.count = rows end -- eigenvalues vector_norm : REAL is require columns = 1 local i : INTEGER do from i := 1 Result := 0.0 until i > rows loop Result := Result + element(i,1)^2 i := i + 1 end Result := Result.sqrt end -- vector_norm feature {NONE, MATRIX_TEST} backsub(lu : LU_DECOMPOSITION[like Current]; vector : like Current) : like Current is local i, j, n : INTEGER sum : REAL x, y, b : like Current do n := rows create x.make(n,1) create y.make(n,1) -- correct for permutation caused by pivot in decomposition b := lu.permutation * vector -- solve lu.lower * y = b by forward substitution y.set_element(1, 1, b.element(1,1)/lu.lower.element(1,1)) from i := 2 until i > n loop from j := 1 sum := 0 until j > i - 1 loop sum := sum + lu.lower.element(i,j) * y.element(j,1) j := j + 1 end y.set_element(i, 1, (b.element(i,1) - sum)/lu.lower.element(i,i)) i := i + 1 end -- solve lu.upper * x = y by back substitution x.set_element(n, 1, y.element(n,1)/lu.upper.element(n,n)) from i := n until i < 1 loop from j := i + 1 sum := 0 until j > n loop sum := sum + lu.upper.element(i,j) * x.element(j,1) j := j + 1 end x.set_element(i, 1, (y.element(i,1) - sum)/lu.upper.element(i,i)) i := i - 1 end Result := x end --backsub hessenberg_form : like Current is require rows = columns local i : INTEGER alpha : REAL v, v_t : like Current u, u_submatrix : like Current do Result := Current.deep_twin from i := 1 until i >= rows - 1 loop v := Result.submatrix(i + 1, 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 u.make_identity(rows) create u_submatrix.make_identity(v.rows) u_submatrix := u_submatrix - (v * v_t) #* 2 u.set_submatrix(i + 1, rows, i + 1, columns, u_submatrix) Result := u * Result u.transpose Result := Result * u end i := i + 1 end end -- hessenberg_form feature {ANY} -- bookkeeping hash_code : INTEGER is local i, j : INTEGER n : INTEGER_8 do from i := 1 n := 0 Result := 0 until i > rows loop from j := 1 until j > columns loop Result := Result + (element(i, j).hash_code * 2^n).hash_code j := j + 1 n := n + 1 end i := i + 1 end end -- hash_code end -- class REAL_MATRIX
| Download code |
