Changeset 76

Timestamp:

04/18/08 14:01:20 (17 years ago)

Author:

smidl

Message:

oprava chol + invqform

Location:

bdm/math

Files:

• chmat.cpp (modified) (2 diffs)
• chmat.h (modified) (3 diffs)

bdm/math/chmat.cpp

@@ -11 +11 @@
         mat Z;
         mat R;
-        mat Q;
-        Z = concat_vertical ( Ch,mat(w*v) );
-         qr ( Z,Q,R );
+        mat V(1,v.length());
+        V.set_row(0,v*w);
+        Z = concat_vertical ( Ch,V );
+        qr ( Z,R );
         Ch = R ( 0, Ch.rows()-1, 0, Ch.cols()-1 );
 };
-inline mat chmat::to_mat(){mat F=Ch.T()*Ch;return F;};
+inline mat chmat::to_mat() {mat F=Ch.T() *Ch;return F;};
 void chmat::mult_sym ( const mat &C ) {
-it_error("not implemented");
+        it_error ( "not implemented" );
 };
-void chmat::mult_sym ( const mat &C , chmat &U) const {
-it_error("not implemented");
+void chmat::mult_sym ( const mat &C , chmat &U ) const {
+        it_error ( "not implemented" );
 };
 void chmat::mult_sym_t ( const mat &C ) {
-it_error("not implemented");
+        it_error ( "not implemented" );
 };
 void chmat::mult_sym_t ( const mat &C, chmat &U ) const {
-it_error("not implemented");
+        it_error ( "not implemented" );
 };
 double chmat::logdet() const {
         double ldet=0.0; int i;
         //sum of logs of (possibly negative!) diagonal entries
-        for (i=0;i<Ch.rows();i++ ) {ldet+=log ( std::fabs(Ch ( i,i )) );}
+        for ( i=0;i<Ch.rows();i++ ) {ldet+=log ( std::fabs ( Ch ( i,i ) ) );}
         return 2*ldet; //compensate for Ch being sqrt()
 };
@@ -38 +39 @@
 inline vec chmat::sqrt_mult ( const vec &v ) const {vec pom; pom = Ch*v; return pom;};
 inline double chmat::qform ( const vec &v ) const {vec pom; pom = Ch*v; return pom*pom;};
+inline double chmat::invqform ( const vec &v ) const {
+        vec pom(v.length());
+        forward_substitution(Ch.T(),v,pom);
+        return pom*pom;
+};
 inline void chmat::clear() {Ch.clear();};

bdm/math/chmat.h

@@ -21 +21 @@
 /*! \brief Symmetric matrix stored in square root decomposition using upper cholesky
 
-This matrix represent $A=Ch Ch'$ where only the upper triangle is stored;
+This matrix represent \f$A=Ch' Ch\f$ where only the upper triangle \f$Ch\f$ is stored;
 
 */
@@ -39 +39 @@
         vec sqrt_mult ( const vec &v ) const;
         double qform ( const vec &v ) const;
+        double invqform ( const vec &v ) const;
         void clear();
         void add ( const chmat &A2, double w=1.0 ) {};
         //!Inversion in the same form, i.e. cholesky
-        void inv ( chmat &Inv ) const   { ( Inv.Ch ) = itpp::inv ( Ch );}; //ltuinv is wrong?
+        void inv ( chmat &Inv ) const   { ( Inv.Ch ) = itpp::inv ( Ch ).T();}; //Fixme: can be more efficient
         ;
 //      void inv ( mat &Inv );
@@ -52 +53 @@
         //! Default constructor
         chmat ( const vec &v) : sqmat ( v.length() ),Ch ( diag(sqrt(v)) ) {};
+        //! Copy constructor
+        chmat ( const chmat &Ch0) : sqmat ( Ch0.dim),Ch ( Ch0.dim,Ch0.dim ) {Ch=Ch0.Ch;};
         //! Default constructor (m3k:cholform)
         chmat ( const mat &M ) : sqmat ( M.rows() ),Ch ( M.rows(),M.cols() ) {