square_mat.h

Go to the documentation of this file.

#ifndef DC_H
#define DC_H


#include "../itpp_ext.h"
#include "../bdmerror.h"

namespace bdm
{

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 ) {};
        sqmat() : dim ( 0 ) {};
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 fsqmat &M, const ivec &perm ) {
                bdm_error ( "not implemented" );
        }

        fsqmat ( const vec &d ) : sqmat ( d.length() ) {
                M = diag ( d );
        };

        virtual ~fsqmat() {};


        void inv ( fsqmat &Inv ) const;

        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;
        };
        
        operator mat&() {return M;};
        
//              fsqmat& operator = ( const fsqmat &A) {M=A.M; return *this;};
        friend std::ostream &operator<< ( std::ostream &os, const fsqmat &sq );

};


class ldmat: public 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();
        vec sqrt_mult ( const vec &v ) const;


        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;
}

}

#endif // DC_H

---