2 const Real EPS = 1e-7; // so sue me. Hard coded
5 Choleski_decomposition::solve(Vector rhs)const
11 // forward substitution
12 for (int i=0; i < n; i++) {
14 for (int j=0; j < i; j++)
16 y(i) = (rhs(i) - sum)/L(i,i);
18 for (int i=0; i < n; i++)
22 Vector &x(rhs); // using input as return val.
23 for (int i=n-1; i >= 0; i--) {
25 for (int j=i+1; j < n; j++)
27 x(i) = (y(i) - sum)/L(i,i);
33 Standard matrix algorithm.
36 Choleski_decomposition::Choleski_decomposition(Matrix P)
37 : L(P.dim()), D(P.dim())
40 assert((P-P.transposed()).norm()/P.norm() < EPS);
43 for (int k= 0; k < n; k++) {
44 for (int j = 0; j < k; j++){
46 for (int l=0; l < j; l++)
47 sum += L(k,l)*L(j,l)*D(l);
48 L(k,j) = (P(k,j) - sum)/D(j);
52 for (int l=0; l < k; l++)
53 sum += sqr(L(k,l))*D(l);
54 Real d = P(k,k) - sum;
59 assert((original()-P).norm() / P.norm() < EPS);
64 Choleski_decomposition::original() const
68 return L*T*L.transposed();
72 Choleski_decomposition::inverse() const
77 for (int i = 0; i < n; i++) {
79 Vector inv(solve(e_i));
80 for (int j = 0 ; j<n; j++)
85 Matrix I1(n), I2(original());
87 assert((I1-original()*invm).norm()/original.norm() < EPS);