- 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 #include "../bdmerror.h"
    00019
    00020 namespace bdm
    00021 {
    00022
    00023 using namespace itpp;
    00024
    00026 void dydr ( double * r, double *f, double *Dr, double *Df, double *R, int jl, int jh, double *kr, int m, int mx );
    00027
    00029 //TODO can be done via: dtrtri.f from lapack
    00030 mat ltuinv ( const mat &L );
    00031
    00036 class sqmat {
    00037 public:
    00045         virtual void opupdt ( const vec &v, double w ) { bdm_error("not implemented"); };
    00046
    00050         virtual mat to_mat() const { bdm_error("not implemented"); return mat(0,0); }
    00051
    00055         virtual void mult_sym ( const mat &C ) { bdm_error("not implemented"); };
    00056
    00060         virtual void mult_sym_t ( const mat &C ) { bdm_error("not implemented"); }
    00061
    00062
    00067         virtual double logdet() const { bdm_error("not implemented"); return 0;};
    00068
    00074         virtual vec sqrt_mult ( const vec &v )  const { bdm_error("not implemented"); return vec(0); };
    00075
    00080         virtual double qform ( const vec &v ) const { bdm_error("not implemented"); return 0; };
    00081
    00086         virtual double invqform ( const vec &v ) const { bdm_error("not implemented"); return 0; };
    00087
    00088 //      //! easy version of the
    00089 //      sqmat inv();
    00090
    00092         virtual void clear() { bdm_error("not implemented"); };
    00093
    00095         int cols() const {
    00096                 return dim;
    00097         };
    00098
    00100         int rows() const {
    00101                 return dim;
    00102         };
    00103
    00105         virtual ~sqmat() {};
    00107         sqmat ( const int dim0 ) : dim ( dim0 ) {};
    00109         sqmat() : dim ( 0 ) {};
    00110 protected:
    00112         int dim;
    00113 };
    00114
    00115
    00120 class fsqmat: public sqmat {
    00121 protected:
    00123         mat M;
    00124 public:
    00125         void opupdt ( const vec &v, double w );
    00126         mat to_mat() const;
    00127         void mult_sym ( const mat &C );
    00128         void mult_sym_t ( const mat &C );
    00130         void mult_sym ( const mat &C, fsqmat &U ) const;
    00132         void mult_sym_t ( const mat &C, fsqmat &U ) const;
    00133         void clear();
    00134
    00136         fsqmat() {}; // mat will be initialized OK
    00138         fsqmat ( const int dim0 ); // mat will be initialized OK
    00140         fsqmat ( const mat &M );
    00141
    00146         fsqmat ( const fsqmat &M, const ivec &perm ) {
    00147                 bdm_error ( "not implemented" );
    00148         }
    00149
    00151         fsqmat ( const vec &d ) : sqmat ( d.length() ) {
    00152                 M = diag ( d );
    00153         };
    00154
    00156         virtual ~fsqmat() {};
    00157
    00158
    00164         void inv ( fsqmat &Inv ) const;
    00165
    00166         double logdet() const {
    00167                 return log ( det ( M ) );
    00168         };
    00169         double qform ( const  vec &v ) const {
    00170                 return ( v* ( M*v ) );
    00171         };
    00172         double invqform ( const  vec &v ) const {
    00173                 return ( v* ( itpp::inv ( M ) *v ) );
    00174         };
    00175         vec sqrt_mult ( const vec &v ) const {
    00176                 mat Ch = chol ( M );
    00177                 return Ch*v;
    00178         };
    00179
    00181         void add ( const fsqmat &fsq2, double w = 1.0 ) {
    00182                 M += fsq2.M;
    00183         };
    00184
    00186         void setD ( const vec &nD ) {
    00187                 M = diag ( nD );
    00188         }
    00190         vec getD () {
    00191                 return diag ( M );
    00192         }
    00194         void setD ( const vec &nD, int i ) {
    00195                 for ( int j = i; j < nD.length(); j++ ) {
    00196                         M ( j, j ) = nD ( j - i );    //Fixme can be more general
    00197                 }
    00198         }
    00199
    00200
    00202         fsqmat& operator += ( const fsqmat &A ) {
    00203                 M += A.M;
    00204                 return *this;
    00205         };
    00207         fsqmat& operator -= ( const fsqmat &A ) {
    00208                 M -= A.M;
    00209                 return *this;
    00210         };
    00212         fsqmat& operator *= ( double x ) {
    00213                 M *= x;
    00214                 return *this;
    00215         };
    00216
    00218         operator mat&() {return M;};
    00219
    00220 //              fsqmat& operator = ( const fsqmat &A) {M=A.M; return *this;};
    00222         friend std::ostream &operator<< ( std::ostream &os, const fsqmat &sq );
    00223
    00224 };
    00225
    00226
    00232 class ldmat: public sqmat {
    00233 public:
    00235         ldmat ( const mat &L, const vec &D );
    00237         ldmat ( const mat &V );
    00239         ldmat ( const ldmat &V0, const ivec &perm ) : sqmat ( V0.rows() ) {
    00240                 ldform ( V0.L.get_cols ( perm ), V0.D );
    00241         };
    00243         ldmat ( vec D0 );
    00245         ldmat ();
    00247         ldmat ( const int dim0 );
    00248
    00250         virtual ~ldmat() {};
    00251
    00252         // Reimplementation of compulsory operatios
    00253
    00254         void opupdt ( const vec &v, double w );
    00255         mat to_mat() const;
    00256         void mult_sym ( const mat &C );
    00257         void mult_sym_t ( const mat &C );
    00259         void add ( const ldmat &ld2, double w = 1.0 );
    00260         double logdet() const;
    00261         double qform ( const vec &v ) const;
    00262         double invqform ( const vec &v ) const;
    00263         void clear();
    00264         vec sqrt_mult ( const vec &v ) const;
    00265
    00266
    00270         void inv ( ldmat &Inv ) const;
    00271
    00276         void mult_sym ( const mat &C, ldmat &U ) const;
    00277
    00282         void mult_sym_t ( const mat &C, ldmat &U ) const;
    00283
    00284
    00291         void ldform ( const mat &A, const vec &D0 );
    00292
    00294         void setD ( const vec &nD ) {
    00295                 D = nD;
    00296         }
    00298         void setD ( const vec &nD, int i ) {
    00299                 D.replace_mid ( i, nD );    //Fixme can be more general
    00300         }
    00302         void setL ( const vec &nL ) {
    00303                 L = nL;
    00304         }
    00305
    00307         const vec& _D() const {
    00308                 return D;
    00309         }
    00311         const mat& _L() const {
    00312                 return L;
    00313         }
    00314
    00316         ldmat& operator += ( const ldmat &ldA );
    00318         ldmat& operator -= ( const ldmat &ldA );
    00320         ldmat& operator *= ( double x );
    00321
    00323         friend std::ostream &operator<< ( std::ostream &os, const ldmat &sq );
    00324
    00325 protected:
    00327         vec D;
    00329         mat L;
    00330
    00331 };
    00332
    00335 inline ldmat& ldmat::operator += ( const ldmat & ldA )  {
    00336         this->add ( ldA );
    00337         return *this;
    00338 }
    00340 inline ldmat& ldmat::operator -= ( const ldmat & ldA )  {
    00341         this->add ( ldA, -1.0 );
    00342         return *this;
    00343 }
    00344
    00345 }
    00346
    00347 #endif // DC_H