- Main Page
- Related Pages
- Classes
- Files

- File List
- File Members

# 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
    00019 namespace bdm
    00020 {
    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
    00049         virtual mat to_mat() const = 0;
    00050
    00054         virtual void mult_sym ( const mat &C ) = 0;
    00055
    00059         virtual void mult_sym_t ( const mat &C ) = 0;
    00060
    00061
    00066         virtual double logdet() const = 0;
    00067
    00073         virtual vec sqrt_mult ( const vec &v ) const = 0;
    00074
    00079         virtual double qform ( const vec &v ) const = 0;
    00080
    00085         virtual double invqform ( const vec &v ) const = 0;
    00086
    00087 //      //! easy version of the
    00088 //      sqmat inv();
    00089
    00091         virtual void clear() = 0;
    00092
    00094         int cols() const {
    00095                 return dim;
    00096         };
    00097
    00099         int rows() const {
    00100                 return dim;
    00101         };
    00102
    00104         virtual ~sqmat() {};
    00106         sqmat ( const int dim0 ) : dim ( dim0 ) {};
    00108         sqmat() : dim ( 0 ) {};
    00109 protected:
    00111         int dim;
    00112 };
    00113
    00114
    00119 class fsqmat: public sqmat {
    00120 protected:
    00122         mat M;
    00123 public:
    00124         void opupdt ( const vec &v, double w );
    00125         mat to_mat() const;
    00126         void mult_sym ( const mat &C );
    00127         void mult_sym_t ( const mat &C );
    00129         void mult_sym ( const mat &C, fsqmat &U ) const;
    00131         void mult_sym_t ( const mat &C, fsqmat &U ) const;
    00132         void clear();
    00133
    00135         fsqmat() {}; // mat will be initialized OK
    00137         fsqmat ( const int dim0 ); // mat will be initialized OK
    00139         fsqmat ( const mat &M );
    00141         fsqmat ( const fsqmat &M, const ivec &perm ) : sqmat ( M.rows() ) {
    00142                 it_error ( "not implemneted" );
    00143         };
    00145         fsqmat ( const vec &d ) : sqmat ( d.length() ) {
    00146                 M = diag ( d );
    00147         };
    00148
    00150         virtual ~fsqmat() {};
    00151
    00152
    00158         void inv ( fsqmat &Inv ) const;
    00159
    00160         double logdet() const {
    00161                 return log ( det ( M ) );
    00162         };
    00163         double qform ( const  vec &v ) const {
    00164                 return ( v* ( M*v ) );
    00165         };
    00166         double invqform ( const  vec &v ) const {
    00167                 return ( v* ( itpp::inv ( M ) *v ) );
    00168         };
    00169         vec sqrt_mult ( const vec &v ) const {
    00170                 mat Ch = chol ( M );
    00171                 return Ch*v;
    00172         };
    00173
    00175         void add ( const fsqmat &fsq2, double w = 1.0 ) {
    00176                 M += fsq2.M;
    00177         };
    00178
    00180         void setD ( const vec &nD ) {
    00181                 M = diag ( nD );
    00182         }
    00184         vec getD () {
    00185                 return diag ( M );
    00186         }
    00188         void setD ( const vec &nD, int i ) {
    00189                 for ( int j = i; j < nD.length(); j++ ) {
    00190                         M ( j, j ) = nD ( j - i );    //Fixme can be more general
    00191                 }
    00192         }
    00193
    00194
    00196         fsqmat& operator += ( const fsqmat &A ) {
    00197                 M += A.M;
    00198                 return *this;
    00199         };
    00201         fsqmat& operator -= ( const fsqmat &A ) {
    00202                 M -= A.M;
    00203                 return *this;
    00204         };
    00206         fsqmat& operator *= ( double x ) {
    00207                 M *= x;
    00208                 return *this;
    00209         };
    00210 //              fsqmat& operator = ( const fsqmat &A) {M=A.M; return *this;};
    00212         friend std::ostream &operator<< ( std::ostream &os, const fsqmat &sq );
    00213
    00214 };
    00215
    00221 class ldmat: public sqmat {
    00222 public:
    00224         ldmat ( const mat &L, const vec &D );
    00226         ldmat ( const mat &V );
    00228         ldmat ( const ldmat &V0, const ivec &perm ) : sqmat ( V0.rows() ) {
    00229                 ldform ( V0.L.get_cols ( perm ), V0.D );
    00230         };
    00232         ldmat ( vec D0 );
    00234         ldmat ();
    00236         ldmat ( const int dim0 );
    00237
    00239         virtual ~ldmat() {};
    00240
    00241         // Reimplementation of compulsory operatios
    00242
    00243         void opupdt ( const vec &v, double w );
    00244         mat to_mat() const;
    00245         void mult_sym ( const mat &C );
    00246         void mult_sym_t ( const mat &C );
    00248         void add ( const ldmat &ld2, double w = 1.0 );
    00249         double logdet() const;
    00250         double qform ( const vec &v ) const;
    00251         double invqform ( const vec &v ) const;
    00252         void clear();
    00253         vec sqrt_mult ( const vec &v ) const;
    00254
    00255
    00259         void inv ( ldmat &Inv ) const;
    00260
    00265         void mult_sym ( const mat &C, ldmat &U ) const;
    00266
    00271         void mult_sym_t ( const mat &C, ldmat &U ) const;
    00272
    00273
    00280         void ldform ( const mat &A, const vec &D0 );
    00281
    00283         void setD ( const vec &nD ) {
    00284                 D = nD;
    00285         }
    00287         void setD ( const vec &nD, int i ) {
    00288                 D.replace_mid ( i, nD );    //Fixme can be more general
    00289         }
    00291         void setL ( const vec &nL ) {
    00292                 L = nL;
    00293         }
    00294
    00296         const vec& _D() const {
    00297                 return D;
    00298         }
    00300         const mat& _L() const {
    00301                 return L;
    00302         }
    00303
    00305         ldmat& operator += ( const ldmat &ldA );
    00307         ldmat& operator -= ( const ldmat &ldA );
    00309         ldmat& operator *= ( double x );
    00310
    00312         friend std::ostream &operator<< ( std::ostream &os, const ldmat &sq );
    00313
    00314 protected:
    00316         vec D;
    00318         mat L;
    00319
    00320 };
    00321
    00324 inline ldmat& ldmat::operator += ( const ldmat & ldA )  {
    00325         this->add ( ldA );
    00326         return *this;
    00327 }
    00329 inline ldmat& ldmat::operator -= ( const ldmat & ldA )  {
    00330         this->add ( ldA, -1.0 );
    00331         return *this;
    00332 }
    00333
    00334 }
    00335
    00336 #endif // DC_H