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
++) {
24 for (int i
=n
-1; i
>= 0; i
--) {
26 for (int j
=i
+1; j
< n
; j
++)
28 x(i
) = (y(i
) - sum
)/L(i
,i
);
34 Standard matrix algorithm.
37 Choleski_decomposition::Choleski_decomposition(Matrix P
)
38 : L(P
.dim()), D(P
.dim())
41 assert((P
-P
.transposed()).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() < 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() < EPS
);