Changeset 437

Timestamp:

07/28/09 15:31:57 (15 years ago)

Author:

smidl

Message:

missing functions in chmat - square_matrix test now passes for chmat

Location:

library

Files:

• bdm/math/chmat.cpp (modified) (1 diff)
• bdm/math/chmat.h (modified) (2 diffs)
• tests/square_mat_test.cpp (modified) (1 diff)

library/bdm/math/chmat.cpp

@@ -17 +17 @@
 mat chmat::to_mat() const {mat F=Ch.T() *Ch;return F;};
 void chmat::mult_sym ( const mat &C ) {
-        it_error ( "not implemented" );
+        it_assert_debug(C.cols()==dim, "Wrong dimension of U");
+        if(!qr(Ch*C.T(), Ch)) {it_warning("QR unstable in chmat mult_sym");}
 };
 void chmat::mult_sym ( const mat &C , chmat &U ) const {
-        mat Z=C*Ch;
-        U.Ch= chol(Z*Z.T());
+        it_assert_debug(C.cols()==U.dim, "Wrong dimension of U");
+        if(!qr(Ch*C.T(), U.Ch)) {it_warning("QR unstable in chmat mult_sym");}
 };
 void chmat::mult_sym_t ( const mat &C ) {
-        it_error ( "not implemented" );
+        it_assert_debug(C.rows()==dim, "Wrong dimension of U");
+        if(!qr(Ch*C, Ch)) {it_warning("QR unstable in chmat mult_sym");}
 };
 void chmat::mult_sym_t ( const mat &C, chmat &U ) const {
-        it_error ( "not implemented" );
+        it_assert_debug(C.rows()==U.dim, "Wrong dimension of U");
+        if(!qr(Ch*C, U.Ch)) {it_warning("QR unstable in chmat mult_sym");}
 };
 double chmat::logdet() const {

library/bdm/math/chmat.h

@@ -40 +40 @@
         void clear();
         //! add another chmat \c A2 with weight \c w.
-        void add ( const chmat &A2, double w=1.0 ) {};
+        void add ( const chmat &A2, double w=1.0 ) {
+                it_assert_debug(dim==A2.dim,"Matrices of unequal dimension");
+                mat pre=concat_vertical(Ch,sqrt(w)*A2.Ch); 
+                mat post=zeros(pre.rows(),pre.cols()); 
+                if(!qr(pre,post)) {it_warning("Unstable QR in chmat add");}
+                Ch=post(0,dim-1,0,dim-1);
+        };
         //!Inversion in the same form, i.e. cholesky
         void inv ( chmat &Inv ) const   { ( Inv.Ch ) = itpp::inv ( Ch ).T();}; //Fixme: can be more efficient
@@ -80 +86 @@
         //! Operators
         chmat& operator += ( const chmat &A2 );
-        chmat& operator -= ( const chmat &A2 ) ;
-        chmat& operator * ( const chmat &A2 );
-        chmat& operator * ( const double &d ){Ch*d; return *this;};
+        chmat& operator -= ( const chmat &A2 );
+        chmat& operator * ( const double &d ){Ch*sqrt(d); return *this;};
         chmat& operator = ( const chmat &A2 ){Ch=A2.Ch;dim=A2.dim;return *this;}
+        chmat& operator *= ( double x ){Ch*=sqrt(x); return *this;};
 };
 

library/tests/square_mat_test.cpp

@@ -100 +100 @@
 
 TEST(test_chmat) {
-    test_square_matrix_minimum<chmat>(epsilon);
+    test_square_matrix<chmat>(epsilon);
 }