4 Choleski_decomposition::solve(Vector rhs
)const
10 // forward substitution
11 for (int i
=0; i
< n
; i
++) {
13 for (int j
=0; j
< i
; j
++)
15 y(i
) = (rhs(i
) - sum
)/L(i
,i
);
17 for (int i
=0; i
< n
; i
++) {
22 for (int i
=n
-1; i
>= 0; i
--) {
24 for (int j
=i
+1; j
< n
; j
++)
26 x(i
) = (y(i
) - sum
)/L(i
,i
);
31 Choleski_decomposition::Choleski_decomposition(Matrix P
)
32 : L(P
.dim()), D(P
.dim())
35 assert((P
-P
.transposed()).norm() < EPS
);
37 for (int k
= 0; k
< n
; k
++) {
38 for (int j
= 0; j
< k
; j
++){
40 for (int l
=0; l
< j
; l
++)
41 sum
+= L(k
,l
)*L(j
,l
)*D(l
);
42 L(k
,j
) = (P(k
,j
) - sum
)/D(j
);
46 for (int l
=0; l
< k
; l
++)
47 sum
+= sqr(L(k
,l
))*D(l
);
48 Real d
= P(k
,k
) - sum
;
53 assert((original()-P
).norm() < EPS
);
58 Choleski_decomposition::original() const
62 return L
*T
*L
.transposed();
66 Choleski_decomposition::inverse() const
71 for (int i
= 0; i
< n
; i
++) {
73 Vector
inv(solve(e_i
));
74 for (int j
= 0 ; j
<n
; j
++)
79 Matrix
I1(n
), I2(original());
81 assert((I1
-original()*invm
).norm() < EPS
);