2 choleski.cc -- implement Choleski_decomposition
4 source file of the Flower Library
6 (c) 1997 Han-Wen Nienhuys <hanwen@stack.nl>
10 const Real EPS = 1e-7; // so sue me. Hard coded
12 // for testing new Matrix_storage.
16 Choleski_decomposition::full_matrix_solve (Vector &out, Vector const &rhs)const
19 assert (n == L.dim());
24 // forward substitution
25 for (int i=0; i < n; i++) {
27 for (int j=0; j < i; j++)
28 sum += y (j) * L(i,j);
29 y (i) = (rhs (i) - sum)/L(i,i);
32 for (int i=0; i < n; i++)
36 Vector &x (out); // using input as return val.
37 for (int i=n-1; i >= 0; i--) {
39 for (int j=i+1; j < n; j++)
41 x (i) = (y (i) - sum)/L(i,i);
46 Choleski_decomposition::band_matrix_solve (Vector &out, Vector const &rhs)const
50 assert (n == L.dim());
57 // forward substitution
58 for (int i=0; i < n; i++) {
60 for (int j= 0 >? i - b; j < i; j++)
61 sum += y (j) * L(i,j);
62 y (i) = (rhs (i) - sum)/L(i,i);
64 for (int i=0; i < n; i++)
68 Vector &x (out); // using input as return val.
69 for (int i=n-1; i >= 0; i--) {
71 for (int j=i+1; j <= i + b&&j < n ; j++)
73 x (i) = (y (i) - sum)/L(i,i);
78 Choleski_decomposition::solve (Vector &x, Vector const &rhs)const
81 band_matrix_solve (x,rhs);
83 full_matrix_solve (x,rhs);
87 Choleski_decomposition::solve (Vector rhs)const
95 Choleski_decomposition::full_matrix_decompose (Matrix const & P)
100 for (int k= 0; k < n; k++) {
101 for (int j = 0; j < k; j++){
103 for (int l=0; l < j; l++)
104 sum += L(k,l)*L(j,l)*D(l);
105 L(k,j) = (P(k,j) - sum)/D(j);
109 for (int l=0; l < k; l++)
110 sum += sqr (L(k,l))*D(l);
111 Real d = P(k,k) - sum;
118 Choleski_decomposition::band_matrix_decompose (Matrix const &P)
124 for (int i= 0; i < n; i++) {
125 for (int j = 0 >? i - b; j < i; j++){
127 for (int l=0 >? i - b; l < j; l++)
128 sum += L(i,l)*L(j,l)*D(l);
129 L(i,j) = (P(i,j) - sum)/D(j);
133 for (int l=0 >? i - b; l < i; l++)
134 sum += sqr (L(i,l))*D(l);
135 Real d = P(i,i) - sum;
139 assert ( L.band_i() == P.band_i ());
146 Standard matrix algorithm.
149 Choleski_decomposition::Choleski_decomposition (Matrix const & P)
150 : L(P.dim()), D(P.dim ())
153 assert ((P-P.transposed()).norm ()/P.norm () < EPS);
156 band_matrix_decompose (P);
158 full_matrix_decompose (P);
162 assert ((original()-P).norm () / P.norm () < EPS);
167 Choleski_decomposition::original() const
171 return L*T*L.transposed();
175 Choleski_decomposition::inverse() const
181 for (int i = 0; i < n; i++) {
184 for (int j = 0 ; j<n; j++)
185 invm (i,j) = inv (j);
189 Matrix I1(n), I2(original());
191 assert ((I1-I2*invm).norm()/I2.norm () < EPS);