work/git/mixpp/bdm/math/libDC.h

Go to the documentation of this file.

#ifndef DC_H
#define DC_H

#include <itpp/itbase.h>
#include "../itpp_ext.h"

using namespace itpp;

void dydr( double * r, double *f, double *Dr, double *Df, double *R, int jl, int jh, double *kr, int m, int mx );

//TODO can be done via: dtrtri.f from lapack
mat ltuinv( const mat &L );

class sqmat
{
        public:
                virtual void opupdt ( const vec &v, double w ) =0;

                virtual mat to_mat() const =0;

                virtual void mult_sym ( const mat &C ) =0;
                
                virtual void mult_sym_t ( const mat &C ) =0;


                virtual double logdet() const =0;

                virtual vec sqrt_mult (const vec &v ) const =0;

                virtual double qform (const vec &v ) const =0;

                virtual double invqform (const vec &v ) const =0;

//      //! easy version of the
//      sqmat inv();

                virtual void clear() =0;

                int cols() const {return dim;};

                int rows() const {return dim;};

                virtual ~sqmat(){};
                sqmat(const int dim0): dim(dim0){};
        protected:
                int dim;
};


class fsqmat: public sqmat
{
        protected:
                mat M;
        public:
                void opupdt ( const vec &v, double w );
                mat to_mat() const;
                void mult_sym ( const mat &C);
                void mult_sym_t ( const mat &C);
                void mult_sym ( const mat &C, fsqmat &U) const;
                void mult_sym_t ( const mat &C, fsqmat &U) const;
                void clear();

                fsqmat(); // mat will be initialized OK
                fsqmat(const int dim0); // mat will be initialized OK
                fsqmat ( const mat &M );
                fsqmat ( const vec &d ):sqmat(d.length()){M=diag(d);};

                virtual ~fsqmat(){};


                virtual void inv ( fsqmat &Inv );

                double logdet() const {return log ( det ( M ) );};
                double qform (const  vec &v ) const {return ( v* ( M*v ) );};
                double invqform (const  vec &v ) const {return ( v* ( itpp::inv(M)*v ) );};
                vec sqrt_mult (const vec &v ) const {mat Ch=chol(M); return Ch*v;};

                void add ( const fsqmat &fsq2, double w=1.0 ){M+=fsq2.M;};

                void setD (const vec &nD){M=diag(nD);}
                vec getD (){return diag(M);}
                void setD (const vec &nD, int i){for(int j=i;j<nD.length();j++){M(j,j)=nD(j-i);}} //Fixme can be more general


                fsqmat& operator += ( const fsqmat &A ) {M+=A.M;return *this;};
                fsqmat& operator -= ( const fsqmat &A ) {M-=A.M;return *this;};
                fsqmat& operator *= ( double x ) {M*=x;return *this;};
//              fsqmat& operator = ( const fsqmat &A) {M=A.M; return *this;};
                friend std::ostream &operator<< ( std::ostream &os, const fsqmat &sq );

};

class ldmat: sqmat
{
        public:
                ldmat ( const mat &L, const vec &D );
                ldmat (const mat &V );
                ldmat (const ldmat &V0, const ivec &perm):sqmat(V0.rows()){     ldform(V0.L.get_cols(perm), V0.D);};
                ldmat ( vec D0 );
                ldmat ();
                ldmat(const int dim0);
                
                virtual ~ldmat(){};

                // Reimplementation of compulsory operatios

                void opupdt ( const vec &v, double w );
                mat to_mat() const;
                void mult_sym ( const mat &C);
                void mult_sym_t ( const mat &C);
                void add ( const ldmat &ld2, double w=1.0 );
                double logdet() const;
                double qform (const vec &v ) const;
                double invqform (const vec &v ) const;
                void clear();
                int cols() const;
                int rows() const;
                vec sqrt_mult ( const vec &v ) const;

                
                virtual void inv ( ldmat &Inv ) const;

                void mult_sym ( const mat &C, ldmat &U) const;

                void mult_sym_t ( const mat &C, ldmat &U) const;


                void ldform (const mat &A,const vec &D0 );

                void setD (const vec &nD){D=nD;}
                void setD (const vec &nD, int i){D.replace_mid(i,nD);} //Fixme can be more general
                void setL (const vec &nL){L=nL;}

                const vec& _D() const {return D;}
                const mat& _L() const {return L;}

                ldmat& operator += ( const ldmat &ldA );
                ldmat& operator -= ( const ldmat &ldA );
                ldmat& operator *= ( double x );

                friend std::ostream &operator<< ( std::ostream &os, const ldmat &sq );

        protected:
                vec D;
                mat L;

};

inline ldmat& ldmat::operator += ( const ldmat &ldA )  {this->add ( ldA );return *this;}
inline ldmat& ldmat::operator -= ( const ldmat &ldA )  {this->add ( ldA,-1.0 );return *this;}
inline int ldmat::cols() const {return dim;}
inline int ldmat::rows() const {return dim;}

#endif // DC_H

---