root/libDC.h @ 12

/*!
 * \file
 * \brief Matrices in decomposed forms (LDL', LU, UDU', etc).
 * \author Vaclav Smidl.
 *
 * -----------------------------------
 * BDM++ - C++ library for Bayesian Decision Making under Uncertainty
 *
 * Using IT++ for numerical operations
 * -----------------------------------
 */

#ifndef DC_H
#define DC_H

#include <itpp/itbase.h>

using namespace itpp;

/*! \brief Virtual class for representation of double symmetric matrices in square-root form.

All operations defined on this class should be optimized for the chosed decomposition.
*/
class sqmat {
public:
        /*!
         * Perfroms a rank-1 update by outer product of vectors: $V = V + w v v'$.
         * @param  v Vector forming the outer product to be added
         * @param  w weight of updating; can be negative

         BLAS-2b operation.
         */
        virtual void opupdt( const vec &v, double w ) =0;

        /*! \brief Conversion to full matrix.
        */

        virtual mat to_mat() =0;

        /*! \brief Inplace symmetric multiplication by a SQUARE matrix $C$, i.e. $V = C*V*C'$
        @param C multiplying matrix,
        @param trans if true, product $V = C'*V*C$ will be computed instead;
        */
        virtual void mult_sym( const mat &C, bool trans=true ) =0;

        
        /*!
        \brief Logarithm of a determinant.
        
        */
        virtual double logdet() =0;

        /*!
        \brief Multiplies square root of $V$ by vector $x$.
        
        Used e.g. in generating normal samples.
        */
        virtual vec sqrt_mult(vec &v) =0;
        
        /*!
        \brief Evaluates quadratic form $x= v'*V*v$;
        
        */
        virtual double qform(vec &v) =0;

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

        //! Clearing matrix so that it corresponds to zeros.
        virtual void clear() =0;
        
        //! Reimplementing common functions of mat: cols().
        virtual int cols() =0;

        //! Reimplementing common functions of mat: cols().
        virtual int rows() =0;

protected:
        int dim;
};


/*! \brief Fake sqmat. This class maps sqmat operations to operations on full matrix.

This class can be used to compare performance of algorithms using decomposed matrices with perormance of the same algorithms using full matrices;
*/
class fsqmat: sqmat {
        void opupdt( const vec &v, double w );
        mat to_mat();
        void mult_sym( const mat &C, bool trans=false );
   void mult_sym( const mat &C, fsqmat &U, bool trans=false );
        void inv(fsqmat &Inv);
        void clear();
        
        
        //! Constructor
        fsqmat(const mat &M);
        
        /*! \brief Matrix inversion preserving the chosen form.

        @param Inv a space where the inverse is stored.
        
        */
        virtual void inv(fsqmat* Inv);
        

};

class ldmat: sqmat {
public:

        //! Construct by copy of L and D.
        ldmat( const mat &L, const vec &D );
        //! Construct by decomposition of full matrix V.
        ldmat( mat V );
        ldmat ();

        // Reimplementation of compulsory operatios

        void opupdt( const vec &v, double w );
        mat to_mat();
   void mult_sym( const mat &C, bool trans=false );
        void add ( const ldmat &ld2, double w=1.0 );
        double logdet();
        double qform(vec &v);
//      sqmat& operator -= ( const sqmat & ld2 );
        void clear();
        int cols();
        int rows();
        vec sqrt_mult(vec &v);

        /*! \brief Matrix inversion preserving the chosen form.

        @param Inv a space where the inverse is stored.
        
        */
        virtual void inv(ldmat &Inv);
        
        /*! \brief Symmetric multiplication of $U$ by a general matrix $C$, result of which is stored in the current class.

        @param Inv a space where the inverse is stored.
        
        */
   void mult_sym( const mat &C, ldmat &U, bool trans=false );

        /*! \brief Transforms general $A'D0A$ into pure $L'DL$
        
        The new decomposition fullfills: $A'*diag(D)*A = self.L'*diag(self.D)*self.L$

        @param A general matrix
        @param D0 general vector
        
        */
   void ldform( mat &A, vec &D0 );

        ldmat& operator += (const ldmat &ldA);
        ldmat& operator -= (const ldmat &ldA);
        ldmat& operator *= (double x);
        
        friend std::ostream &operator<< ( std::ostream &os, ldmat &sq );

protected:
        vec D;
        mat L;

};

//////// Operations:

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(){return L.cols();}
inline int ldmat::rows(){return L.rows();}

#endif // DC_H