square_mat.h
Go to the documentation of this file.00001 00013 #ifndef DC_H 00014 #define DC_H 00015 00016 00017 #include "../itpp_ext.h" 00018 #include "../bdmerror.h" 00019 00020 namespace bdm { 00021 00022 using namespace itpp; 00023 00025 void dydr ( double * r, double *f, double *Dr, double *Df, double *R, int jl, int jh, double *kr, int m, int mx ); 00026 00028 //TODO can be done via: dtrtri.f from lapack 00029 mat ltuinv ( const mat &L ); 00030 00035 class sqmat { 00036 public: 00044 virtual void opupdt ( const vec &v, double w ) = 0; 00045 00048 virtual mat to_mat() const = 0; 00049 00053 virtual void mult_sym ( const mat &C ) = 0; 00054 00058 virtual void mult_sym_t ( const mat &C ) = 0; 00059 00064 virtual double logdet() const = 0; 00065 00071 virtual vec sqrt_mult ( const vec &v ) const = 0; 00072 00077 virtual double qform ( const vec &v ) const = 0; 00078 00083 virtual double invqform ( const vec &v ) const = 0; 00084 00085 // //! easy version of the 00086 // sqmat inv(); 00087 00089 virtual void clear() = 0; 00090 00092 int cols() const { 00093 return dim; 00094 }; 00095 00097 int rows() const { 00098 return dim; 00099 }; 00100 00102 virtual ~sqmat() {}; 00104 sqmat ( const int dim0 ) : dim ( dim0 ) {}; 00106 sqmat() : dim ( 0 ) {}; 00107 protected: 00109 int dim; 00110 }; 00111 00112 00117 class fsqmat: public sqmat { 00118 protected: 00120 mat M; 00121 public: 00122 void opupdt ( const vec &v, double w ); 00123 mat to_mat() const; 00124 void mult_sym ( const mat &C ); 00125 void mult_sym_t ( const mat &C ); 00127 void mult_sym ( const mat &C, fsqmat &U ) const; 00129 void mult_sym_t ( const mat &C, fsqmat &U ) const; 00130 void clear(); 00131 00133 fsqmat() {}; // mat will be initialized OK 00135 fsqmat ( const int dim0 ); // mat will be initialized OK 00137 fsqmat ( const mat &M ); 00138 00143 fsqmat ( const fsqmat &M, const ivec &perm ) { 00144 bdm_error ( "not implemented" ); 00145 } 00146 00148 fsqmat ( const vec &d ) : sqmat ( d.length() ) { 00149 M = diag ( d ); 00150 }; 00151 00153 virtual ~fsqmat() {}; 00154 00155 00161 void inv ( fsqmat &Inv ) const; 00162 00163 double logdet() const { 00164 return log ( det ( M ) ); 00165 }; 00166 double qform ( const vec &v ) const { 00167 return ( v* ( M*v ) ); 00168 }; 00169 double invqform ( const vec &v ) const { 00170 return ( v* ( itpp::inv ( M ) *v ) ); 00171 }; 00172 vec sqrt_mult ( const vec &v ) const { 00173 mat Ch = chol ( M ); 00174 return Ch*v; 00175 }; 00176 00178 void add ( const fsqmat &fsq2, double w = 1.0 ) { 00179 M += fsq2.M; 00180 }; 00181 00183 void setD ( const vec &nD ) { 00184 M = diag ( nD ); 00185 } 00187 vec getD () { 00188 return diag ( M ); 00189 } 00191 void setD ( const vec &nD, int i ) { 00192 for ( int j = i; j < nD.length(); j++ ) { 00193 M ( j, j ) = nD ( j - i ); //Fixme can be more general 00194 } 00195 } 00196 00197 00199 fsqmat& operator += ( const fsqmat &A ) { 00200 M += A.M; 00201 return *this; 00202 }; 00204 fsqmat& operator -= ( const fsqmat &A ) { 00205 M -= A.M; 00206 return *this; 00207 }; 00209 fsqmat& operator *= ( double x ) { 00210 M *= x; 00211 return *this; 00212 }; 00213 00215 operator mat&() { 00216 return M; 00217 }; 00218 00219 // fsqmat& operator = ( const fsqmat &A) {M=A.M; return *this;}; 00221 friend std::ostream &operator<< ( std::ostream &os, const fsqmat &sq ); 00223 mat & _M ( ) { 00224 return M; 00225 }; 00226 00227 }; 00228 00229 00235 class ldmat: public sqmat { 00236 public: 00238 ldmat ( const mat &L, const vec &D ); 00240 ldmat ( const mat &V ); 00242 ldmat ( const ldmat &V0, const ivec &perm ) : sqmat ( V0.rows() ) { 00243 ldform ( V0.L.get_cols ( perm ), V0.D ); 00244 }; 00246 ldmat ( vec D0 ); 00248 ldmat (); 00250 ldmat ( const int dim0 ); 00251 00253 virtual ~ldmat() {}; 00254 00255 // Reimplementation of compulsory operatios 00256 00257 void opupdt ( const vec &v, double w ); 00258 mat to_mat() const; 00259 void mult_sym ( const mat &C ); 00260 void mult_sym_t ( const mat &C ); 00262 void add ( const ldmat &ld2, double w = 1.0 ); 00263 double logdet() const; 00264 double qform ( const vec &v ) const; 00265 double invqform ( const vec &v ) const; 00266 void clear(); 00267 vec sqrt_mult ( const vec &v ) const; 00268 00269 00273 void inv ( ldmat &Inv ) const; 00274 00279 void mult_sym ( const mat &C, ldmat &U ) const; 00280 00285 void mult_sym_t ( const mat &C, ldmat &U ) const; 00286 00287 00294 void ldform ( const mat &A, const vec &D0 ); 00295 00297 void setD ( const vec &nD ) { 00298 D = nD; 00299 } 00301 void setD ( const vec &nD, int i ) { 00302 D.replace_mid ( i, nD ); //Fixme can be more general 00303 } 00305 void setL ( const vec &nL ) { 00306 L = nL; 00307 } 00308 00310 const vec& _D() const { 00311 return D; 00312 } 00314 const mat& _L() const { 00315 return L; 00316 } 00318 vec& __D() { 00319 return D; 00320 } 00322 mat& __L() { 00323 return L; 00324 } 00325 void validate() { 00326 bdm_assert(L.rows()==D.length(),"Incompatible L and D in ldmat"); 00327 dim= L.rows(); 00328 } 00329 00331 ldmat& operator += ( const ldmat &ldA ); 00333 ldmat& operator -= ( const ldmat &ldA ); 00335 ldmat& operator *= ( double x ); 00336 00338 friend std::ostream &operator<< ( std::ostream &os, const ldmat &sq ); 00339 00340 protected: 00342 vec D; 00344 mat L; 00345 00346 }; 00347 00350 inline ldmat& ldmat::operator += ( const ldmat & ldA ) { 00351 this->add ( ldA ); 00352 return *this; 00353 } 00355 inline ldmat& ldmat::operator -= ( const ldmat & ldA ) { 00356 this->add ( ldA, -1.0 ); 00357 return *this; 00358 } 00359 00360 } 00361 00362 #endif // DC_H
Generated on 2 Dec 2013 for mixpp by
1.4.7
