[262] | 1 | |
---|
[39] | 2 | #include "chmat.h" |
---|
| 3 | |
---|
| 4 | //using std::endl; |
---|
| 5 | |
---|
| 6 | |
---|
| 7 | void chmat::opupdt ( const vec &v, double w ) { |
---|
| 8 | //TODO see cholupdt in lhotse |
---|
| 9 | mat Z; |
---|
| 10 | mat R; |
---|
[76] | 11 | mat V(1,v.length()); |
---|
[434] | 12 | V.set_row(0,v*sqrt(w)); |
---|
[76] | 13 | Z = concat_vertical ( Ch,V ); |
---|
| 14 | qr ( Z,R ); |
---|
[39] | 15 | Ch = R ( 0, Ch.rows()-1, 0, Ch.cols()-1 ); |
---|
| 16 | }; |
---|
[168] | 17 | mat chmat::to_mat() const {mat F=Ch.T() *Ch;return F;}; |
---|
[39] | 18 | void chmat::mult_sym ( const mat &C ) { |
---|
[437] | 19 | it_assert_debug(C.cols()==dim, "Wrong dimension of U"); |
---|
| 20 | if(!qr(Ch*C.T(), Ch)) {it_warning("QR unstable in chmat mult_sym");} |
---|
[39] | 21 | }; |
---|
[76] | 22 | void chmat::mult_sym ( const mat &C , chmat &U ) const { |
---|
[437] | 23 | it_assert_debug(C.cols()==U.dim, "Wrong dimension of U"); |
---|
| 24 | if(!qr(Ch*C.T(), U.Ch)) {it_warning("QR unstable in chmat mult_sym");} |
---|
[39] | 25 | }; |
---|
| 26 | void chmat::mult_sym_t ( const mat &C ) { |
---|
[437] | 27 | it_assert_debug(C.rows()==dim, "Wrong dimension of U"); |
---|
| 28 | if(!qr(Ch*C, Ch)) {it_warning("QR unstable in chmat mult_sym");} |
---|
[39] | 29 | }; |
---|
| 30 | void chmat::mult_sym_t ( const mat &C, chmat &U ) const { |
---|
[437] | 31 | it_assert_debug(C.rows()==U.dim, "Wrong dimension of U"); |
---|
| 32 | if(!qr(Ch*C, U.Ch)) {it_warning("QR unstable in chmat mult_sym");} |
---|
[39] | 33 | }; |
---|
| 34 | double chmat::logdet() const { |
---|
| 35 | double ldet=0.0; int i; |
---|
| 36 | //sum of logs of (possibly negative!) diagonal entries |
---|
[76] | 37 | for ( i=0;i<Ch.rows();i++ ) {ldet+=log ( std::fabs ( Ch ( i,i ) ) );} |
---|
[39] | 38 | return 2*ldet; //compensate for Ch being sqrt() |
---|
| 39 | }; |
---|
| 40 | //TODO can be done more efficiently using BLAS, see triangular matrices |
---|
[108] | 41 | vec chmat::sqrt_mult ( const vec &v ) const {vec pom; pom = Ch*v; return pom;}; |
---|
| 42 | double chmat::qform ( const vec &v ) const {vec pom; pom = Ch*v; return pom*pom;}; |
---|
| 43 | double chmat::invqform ( const vec &v ) const { |
---|
[76] | 44 | vec pom(v.length()); |
---|
| 45 | forward_substitution(Ch.T(),v,pom); |
---|
| 46 | return pom*pom; |
---|
| 47 | }; |
---|
[108] | 48 | void chmat::clear() {Ch.clear();}; |
---|