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
);