1 | |
---|
2 | #include "chmat.h" |
---|
3 | |
---|
4 | //using std::endl; |
---|
5 | |
---|
6 | namespace bdm |
---|
7 | { |
---|
8 | |
---|
9 | void chmat::opupdt ( const vec &v, double w ) { |
---|
10 | //TODO see cholupdt in lhotse |
---|
11 | mat Z; |
---|
12 | mat R; |
---|
13 | mat V ( 1, v.length() ); |
---|
14 | V.set_row ( 0, v*sqrt ( w ) ); |
---|
15 | Z = concat_vertical ( Ch, V ); |
---|
16 | qr ( Z, R ); |
---|
17 | Ch = R ( 0, Ch.rows() - 1, 0, Ch.cols() - 1 ); |
---|
18 | } |
---|
19 | |
---|
20 | mat chmat::to_mat() const { |
---|
21 | mat F = Ch.T() * Ch; |
---|
22 | return F; |
---|
23 | } |
---|
24 | |
---|
25 | void chmat::mult_sym ( const mat &C ) { |
---|
26 | bdm_assert_debug ( C.rows() == dim, "Wrong dimension of U" ); |
---|
27 | if ( !qr ( Ch*C.T(), Ch ) ) { |
---|
28 | bdm_warning ( "QR unstable in chmat mult_sym" ); |
---|
29 | } |
---|
30 | } |
---|
31 | |
---|
32 | void chmat::mult_sym ( const mat &C , chmat &U ) const { |
---|
33 | bdm_assert_debug ( C.rows() == U.dim, "Wrong dimension of U" ); |
---|
34 | if ( !qr ( Ch*C.T(), U.Ch ) ) { |
---|
35 | bdm_warning ( "QR unstable in chmat mult_sym" ); |
---|
36 | } |
---|
37 | } |
---|
38 | |
---|
39 | void chmat::mult_sym_t ( const mat &C ) { |
---|
40 | bdm_assert_debug ( C.cols() == dim, "Wrong dimension of U" ); |
---|
41 | if ( !qr ( Ch*C, Ch ) ) { |
---|
42 | bdm_warning ( "QR unstable in chmat mult_sym" ); |
---|
43 | } |
---|
44 | } |
---|
45 | |
---|
46 | void chmat::mult_sym_t ( const mat &C, chmat &U ) const { |
---|
47 | bdm_assert_debug ( C.cols() == U.dim, "Wrong dimension of U" ); |
---|
48 | if ( !qr ( Ch*C, U.Ch ) ) { |
---|
49 | bdm_warning ( "QR unstable in chmat mult_sym" ); |
---|
50 | } |
---|
51 | } |
---|
52 | |
---|
53 | double chmat::logdet() const { |
---|
54 | double ldet = 0.0; |
---|
55 | int i; |
---|
56 | //sum of logs of (possibly negative!) diagonal entries |
---|
57 | for ( i = 0; i < Ch.rows(); i++ ) { |
---|
58 | ldet += log ( std::fabs ( Ch ( i, i ) ) ); |
---|
59 | } |
---|
60 | return 2*ldet; //compensate for Ch being sqrt() |
---|
61 | } |
---|
62 | |
---|
63 | //TODO can be done more efficiently using BLAS, see triangular matrices |
---|
64 | vec chmat::sqrt_mult ( const vec &v ) const { |
---|
65 | vec pom; |
---|
66 | pom = Ch * v; |
---|
67 | return pom; |
---|
68 | } |
---|
69 | |
---|
70 | double chmat::qform ( const vec &v ) const { |
---|
71 | vec pom; |
---|
72 | pom = Ch * v; |
---|
73 | return pom*pom; |
---|
74 | } |
---|
75 | |
---|
76 | double chmat::invqform ( const vec &v ) const { |
---|
77 | vec pom ( v.length() ); |
---|
78 | forward_substitution ( Ch.T(), v, pom ); |
---|
79 | return pom*pom; |
---|
80 | } |
---|
81 | |
---|
82 | void chmat::clear() { |
---|
83 | Ch.clear(); |
---|
84 | } |
---|
85 | |
---|
86 | } |
---|