Changeset 26
Legend:
- Unmodified
- Added
- Removed
-
bdm/math/libDC.cpp
r24 r26 15 15 void fsqmat::opupdt ( const vec &v, double w ) {M+=outer_product ( v,v*w );}; 16 16 mat fsqmat::to_mat() {return M;}; 17 void fsqmat::mult_sym ( const mat &C, bool trans ) {M=C *M*C.T();}; 18 void fsqmat::mult_sym ( const mat &C, fsqmat &U, bool trans ) 19 { 20 mat UM= ( C *(M*C.T()) ); 21 U=UM; 22 }; 17 void fsqmat::mult_sym ( const mat &C) {M=C *M*C.T();}; 18 void fsqmat::mult_sym_t ( const mat &C) {M=C.T() *M*C;}; 19 void fsqmat::mult_sym ( const mat &C, fsqmat &U) { U.M = ( C *(M*C.T()) );}; 20 void fsqmat::mult_sym_t ( const mat &C, fsqmat &U) { U.M = ( C.T() *(M*C) );}; 23 21 void fsqmat::inv ( fsqmat &Inv ) {mat IM = itpp::inv ( M ); Inv=IM;}; 24 22 void fsqmat::clear() {M.clear();}; … … 28 26 M=M0;dim=M0.cols(); 29 27 }; 28 30 29 fsqmat::fsqmat() {}; 30 31 fsqmat::fsqmat(const int dim0): M(dim0,dim0) {}; 32 33 std::ostream &operator<< ( std::ostream &os, fsqmat &ld ) { 34 os << ld.M << endl; 35 return os; 36 } 31 37 32 38 … … 41 47 } 42 48 49 ldmat::ldmat(const int dim0): D(dim0),L(dim0,dim0) { 50 dim = dim0; 51 } 52 43 53 ldmat::ldmat(const vec D0) { 44 54 D = D0; … … 47 57 } 48 58 49 ldmat::ldmat( const mat V ) {59 ldmat::ldmat( const mat &V ) { 50 60 //TODO check if correct!! Based on heuristic observation of lu() 51 61 52 62 dim = V.cols(); 53 mat F;54 vec D0;55 63 it_assert_debug( dim == V.rows(),"ldmat::ldmat matrix V is not square!" ); 64 65 // L and D will be allocated by ldform() 56 66 57 //decompose V in cholesky 58 D0 = ones(dim); 59 60 chol(V,F); 61 // L and D will be allocated by ldform() 62 this->ldform(F,D0); 67 //Chol is unstable 68 this->ldform(chol(V),ones(dim)); 69 // this->ldform(ul(V),ones(dim)); 63 70 } 64 71 … … 135 142 } 136 143 137 void ldmat::mult_sym( const mat &C, bool trans ) { 138 139 //TODO better 140 141 it_assert_debug( C.cols()==L.cols(), "ldmat::mult_sym wrong input argument" ); 142 mat Ct=C; 143 144 if ( trans==false ) { // return C*this*C' 145 Ct *= this->to_mat(); 146 Ct *= C.transpose(); 147 } else { // return C'*this*C 148 Ct = C.transpose(); 149 Ct *= this->to_mat(); 150 Ct *= C; 151 } 152 153 ldmat Lnew=ldmat( Ct ); 154 L = Lnew.L; 155 D = Lnew.D; 156 } 157 158 void ldmat::mult_sym( const mat &C, ldmat &U, bool trans ) { 159 160 //TODO better 161 162 //TODO input test 163 164 mat Ct=C; 165 166 if ( trans==false ) { // return C*this*C' 167 Ct *= U.to_mat(); 168 Ct *= C.transpose(); 169 } else { // return C'*this*C 170 Ct = C.transpose(); 171 Ct *= U.to_mat(); 172 Ct *= C; 173 } 174 175 ldmat Lnew=ldmat( Ct ); 176 L = Lnew.L; 177 D = Lnew.D; 178 } 179 180 double ldmat::logdet() { 144 void ldmat::mult_sym( const mat &C) { 145 mat A = L*C.T(); 146 this->ldform(A,D); 147 } 148 149 void ldmat::mult_sym_t( const mat &C) { 150 mat A = L*C; 151 this->ldform(A,D); 152 } 153 154 void ldmat::mult_sym( const mat &C, ldmat &U) { 155 mat A=L*C.T(); //could be done more efficiently using BLAS 156 U.ldform(A,D); 157 } 158 159 void ldmat::mult_sym_t( const mat &C, ldmat &U) { 160 mat A=L*C; 161 /* vec nD=zeros(U.rows()); 162 nD.replace_mid(0, D); //I case that D < nD*/ 163 U.ldform(A,D); 164 } 165 166 167 double ldmat::logdet() const { 181 168 double ldet = 0.0; 182 169 int i; … … 218 205 } 219 206 220 void ldmat::ldform( mat &A,vec &D0 )207 void ldmat::ldform(const mat &A,const vec &D0 ) 221 208 { 222 209 int m = A.rows(); … … 261 248 cc++; 262 249 L.set_row( n-cc, L.get_row( m+i ) ); 263 L.set_row( m+i, 0);250 L.set_row( m+i,zeros(L.cols()) ); 264 251 D( m+i )=0; L( m+i,i )=1; 265 252 L.set_submatrix( n-cc,m+i-1,i,i,0 ); … … 304 291 jj = D.length()-1-n+ii; 305 292 D(jj) = 0; 306 L.set_row(jj, 0);293 L.set_row(jj,zeros(L.cols())); //TODO: set_row accepts Num_T 307 294 L(jj,jj)=1; 308 295 } -
bdm/math/libDC.h
r24 r26 41 41 /*! \brief Inplace symmetric multiplication by a SQUARE matrix $C$, i.e. $V = C*V*C'$ 42 42 @param C multiplying matrix, 43 @param trans if true, product $V = C'*V*C$ will be computed instead; 44 */ 45 virtual void mult_sym ( const mat &C, bool trans=true ) =0; 43 */ 44 virtual void mult_sym ( const mat &C ) =0; 45 46 /*! \brief Inplace symmetric multiplication by a SQUARE transpose of matrix $C$, i.e. $V = C'*V*C$ 47 @param C multiplying matrix, 48 */ 49 virtual void mult_sym_t ( const mat &C ) =0; 46 50 47 51 … … 50 54 51 55 */ 52 virtual double logdet() =0;56 virtual double logdet() const =0; 53 57 54 58 /*! … … 93 97 void opupdt ( const vec &v, double w ); 94 98 mat to_mat() ; 95 void mult_sym ( const mat &C, bool trans=false ); 96 void mult_sym ( const mat &C, fsqmat &U, bool trans=false ); 99 void mult_sym ( const mat &C); 100 void mult_sym_t ( const mat &C); 101 void mult_sym ( const mat &C, fsqmat &U); 102 void mult_sym_t ( const mat &C, fsqmat &U); 97 103 void clear(); 98 104 99 105 //! Default initialization 100 106 fsqmat(); // mat will be initialized OK 107 //! Default initialization with proper size 108 fsqmat(const int dim0); // mat will be initialized OK 101 109 //! Constructor 102 110 fsqmat ( const mat &M ); … … 109 117 virtual void inv ( fsqmat &Inv ); 110 118 111 double logdet() {return log ( det ( M ) );};119 double logdet() const {return log ( det ( M ) );}; 112 120 double qform (const vec &v ) {return ( v* ( M*v ) );}; 113 121 vec sqrt_mult (const vec &v ) {it_error ( "not implemented" );return v;}; … … 117 125 fsqmat& operator *= ( double x ) {M*=x;return *this;}; 118 126 // fsqmat& operator = ( const fsqmat &A) {M=A.M; return *this;}; 127 128 friend std::ostream &operator<< ( std::ostream &os, fsqmat &sq ); 129 119 130 }; 120 131 … … 126 137 ldmat ( const mat &L, const vec &D ); 127 138 //! Construct by decomposition of full matrix V. 128 ldmat ( matV );139 ldmat (const mat &V ); 129 140 //! Construct diagonal matrix with diagonal D0 130 141 ldmat ( vec D0 ); 142 //!Default constructor 131 143 ldmat (); 144 //! Default initialization with proper size 145 ldmat(const int dim0); 132 146 133 147 // Reimplementation of compulsory operatios … … 135 149 void opupdt ( const vec &v, double w ); 136 150 mat to_mat(); 137 void mult_sym ( const mat &C, bool trans=false ); 151 void mult_sym ( const mat &C); 152 void mult_sym_t ( const mat &C); 138 153 void add ( const ldmat &ld2, double w=1.0 ); 139 double logdet() ;154 double logdet() const; 140 155 double qform (const vec &v ); 141 156 // sqmat& operator -= ( const sqmat & ld2 ); 142 157 void clear(); 143 int cols() ;144 int rows() ;158 int cols() const; 159 int rows() const; 145 160 vec sqrt_mult ( const vec &v ); 146 161 … … 154 169 /*! \brief Symmetric multiplication of $U$ by a general matrix $C$, result of which is stored in the current class. 155 170 156 @param Inv a space where the inverse is stored. 157 158 */ 159 void mult_sym ( const mat &C, ldmat &U, bool trans=false ); 160 161 /*! \brief Transforms general $A'D0A$ into pure $L'DL$ 171 @param U a space where the inverse is stored. 172 173 */ 174 void mult_sym ( const mat &C, ldmat &U); 175 176 /*! \brief Symmetric multiplication of $U$ by a transpose of a general matrix $C$, result of which is stored in the current class. 177 178 @param U a space where the inverse is stored. 179 180 */ 181 void mult_sym_t ( const mat &C, ldmat &U); 182 183 184 /*! \brief Transforms general $A'D0 A$ into pure $L'DL$ 162 185 163 186 The new decomposition fullfills: $A'*diag(D)*A = self.L'*diag(self.D)*self.L$ … … 167 190 168 191 */ 169 void ldform ( mat &A,vec &D0 );192 void ldform (const mat &A,const vec &D0 ); 170 193 171 194 ldmat& operator += ( const ldmat &ldA ); … … 185 208 inline ldmat& ldmat::operator += ( const ldmat &ldA ) {this->add ( ldA );return *this;} 186 209 inline ldmat& ldmat::operator -= ( const ldmat &ldA ) {this->add ( ldA,-1.0 );return *this;} 187 inline int ldmat::cols() {return dim;}188 inline int ldmat::rows() {return dim;}210 inline int ldmat::cols() const {return dim;} 211 inline int ldmat::rows() const {return dim;} 189 212 190 213 #endif // DC_H
