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++)
28 for (int j=0; j < i; j++)
29 sum += y (j) * L(i,j);
30 y (i) = (rhs (i) - sum)/L(i,i);
33 for (int i=0; i < n; i++)
37 Vector &x (out); // using input as return val.
38 for (int i=n-1; i >= 0; i--)
41 for (int j=i+1; j < n; j++)
43 x (i) = (y (i) - sum)/L(i,i);
48 Choleski_decomposition::band_matrix_solve (Vector &out, Vector const &rhs) const
52 assert (n == L.dim());
59 // forward substitution
60 for (int i=0; i < n; i++)
63 for (int j= 0 >? i - b; j < i; j++)
64 sum += y (j) * L(i,j);
65 y (i) = (rhs (i) - sum)/L(i,i);
67 for (int i=0; i < n; i++)
71 Vector &x (out); // using input as return val.
72 for (int i=n-1; i >= 0; i--)
75 for (int j=i+1; j <= i + b&&j < n ; j++)
77 x (i) = (y (i) - sum)/L(i,i);
82 Choleski_decomposition::solve (Vector &x, Vector const &rhs) const
86 band_matrix_solve (x,rhs);
89 full_matrix_solve (x,rhs);
93 Choleski_decomposition::solve (Vector rhs) const
101 Choleski_decomposition::full_matrix_decompose (Matrix const & P)
106 for (int k= 0; k < n; k++)
108 for (int j = 0; j < k; j++)
111 for (int l=0; l < j; l++)
112 sum += L(k,l)*L(j,l)*D(l);
113 L(k,j) = (P(k,j) - sum)/D(j);
117 for (int l=0; l < k; l++)
118 sum += sqr (L(k,l))*D(l);
119 Real d = P(k,k) - sum;
126 Choleski_decomposition::band_matrix_decompose (Matrix const &P)
132 for (int i= 0; i < n; i++)
134 for (int j = 0 >? i - b; j < i; j++)
137 for (int l=0 >? i - b; l < j; l++)
138 sum += L(i,l)*L(j,l)*D(l);
139 L(i,j) = (P(i,j) - sum)/D(j);
143 for (int l=0 >? i - b; l < i; l++)
144 sum += sqr (L(i,l))*D(l);
145 Real d = P(i,i) - sum;
149 assert (L.band_i() == P.band_i ());
156 Standard matrix algorithm.
159 Choleski_decomposition::Choleski_decomposition (Matrix const & P)
160 : L(P.dim()), D(P.dim ())
163 assert ((P-P.transposed()).norm ()/P.norm () < EPS);
166 band_matrix_decompose (P);
168 full_matrix_decompose (P);
172 assert ((original()-P).norm () / P.norm () < EPS);
177 Choleski_decomposition::original() const
181 return L*T*L.transposed();
185 Choleski_decomposition::inverse() const
191 for (int i = 0; i < n; i++)
195 for (int j = 0 ; j<n; j++)
196 invm (i,j) = inv (j);
200 Matrix I1(n), I2(original());
202 assert ((I1-I2*invm).norm()/I2.norm () < EPS);