MMSD_DoubleSymmetricMatrix.h

Go to the documentation of this file.

#ifndef MMSD_DoubleSymmetricMatrix_H
#define MMSD_DoubleSymmetricMatrix_H

#include "MMSD_Object.h"

#include "MMSD_DoubleVector.h"
#include "LAP_DoubleFullSymmetricMatrix.h"
#include "LAP_DoublePackedSymmetricMatrix.h"

DEFINE_SPTR(MMSD_DoubleFullMatrix);

DEFINE_SPTR(MMSD_DoubleSymmetricMatrix);
DEFINE_SVPTR(MMSD_DoubleSymmetricMatrix);

class MMSD_DoubleSymmetricMatrix: public MMSD_Object  {
    
    SP_OBJECT(MMSD_DoubleSymmetricMatrix);

public:
    static const tFlag PACKED_STORAGE;
    static const tFlag FULL_STORAGE;

    
    
private:
    // ATTRIBUTES
    
    // ASSOCIATIONS

    SP::LAP_DoubleFullSymmetricMatrix mFullStorage;
    SP::LAP_DoublePackedSymmetricMatrix mPackedStorage;
    
   
    
    
    // METHOD
  
  
public:
    // CONSTRUCTORS
    MMSD_DoubleSymmetricMatrix();
  
  
    // DESTRUCTORS
    virtual ~MMSD_DoubleSymmetricMatrix();

    
    
    // -------------
    // NEW METHODS
    // -------------
public:

    inline static SP::MMSD_DoubleSymmetricMatrix New() {
        SP::MMSD_DoubleSymmetricMatrix p(new MMSD_DoubleSymmetricMatrix(),
                                            MMSD_DoubleSymmetricMatrix::Delete());
        p->setThis(p);
        return p;
    };

  
   
   
   
    //-----------
    // OPERATORS
    // ----------
  
    inline double& operator()(const size_t& i,const size_t& j) {
        return (mFullStorage.get()!=null)?(*mFullStorage.get())(i,j):(*mPackedStorage.get())(i,j);  
    };
    
    inline const double& operator()(const size_t& i,const size_t& j) const {
        return (mFullStorage.get()!=null)?(*mFullStorage)(i,j):(*mPackedStorage.get())(i,j);  
    };
    MMSD_DoubleSymmetricMatrix& operator*=(const MMSD_DoubleSymmetricMatrix& s) {
        const LAP_DoubleFullSymmetricMatrix *fsy=dynamic_cast< const LAP_DoubleFullSymmetricMatrix *>(&s);
        const LAP_DoublePackedSymmetricMatrix *psy=dynamic_cast< const LAP_DoublePackedSymmetricMatrix *>(&s);
        if (mFullStorage!=null) {
            if (fsy!=null) (*mFullStorage)*=(*fsy);
            else if (psy!=null) (*mFullStorage)*=(*psy);
        } else if (mPackedStorage!=null) {
            if (fsy!=null) (*mPackedStorage)*=(*fsy);
            else if (psy!=null) (*mPackedStorage)*=(*psy);
        }
        return (*this);
    };
    MMSD_DoubleSymmetricMatrix& operator*=(const lapack_real&  s)  {
        if (mFullStorage.get()!=null) (*mFullStorage.get())*=s;
        else if (mPackedStorage.get()!=null) (*mPackedStorage.get())*=s;
        return (*this);
    };
    inline MMSD_DoubleSymmetricMatrix& operator/=(const lapack_real&  s)  {
        if (mFullStorage.get()!=null) (*mFullStorage.get())/=s;
        else if (mPackedStorage.get()!=null) (*mPackedStorage.get())/=s;
        return (*this);
    };
 
    inline MMSD_DoubleSymmetricMatrix& operator+=(const lapack_real&  s)  {
        if (mFullStorage.get()!=null) (*mFullStorage.get())+=s;
        else if (mPackedStorage.get()!=null) (*mPackedStorage.get())+=s;
    
        return (*this);
    };
    inline MMSD_DoubleSymmetricMatrix& operator-=(const lapack_real&  s)  {
        if (mFullStorage.get()!=null) (*mFullStorage.get())-=s;
        else if (mPackedStorage.get()!=null) (*mPackedStorage.get())-=s;
        return (*this);
    };
    inline MMSD_DoubleSymmetricMatrix& operator^=(const lapack_real&  s)  {
        if (mFullStorage.get()!=null) (*mFullStorage.get())^=s;
        else if (mPackedStorage.get()!=null) (*mPackedStorage.get())^=s;
        return (*this);
    };
    
    
    MMSD_DoubleSymmetricMatrix& operator=(const double& s)  {
        if (mFullStorage.get()!=null) (*mFullStorage.get())=s;
        else if (mPackedStorage.get()!=null) (*mPackedStorage.get())=s;
        return (*this);
    };
    MMSD_DoubleSymmetricMatrix& operator=(const MMSD_DoubleSymmetricMatrix& s)  {
        if (mFullStorage.get()!=null) mFullStorage->copy(s.getMatrix());
        else if (mPackedStorage.get()!=null) mPackedStorage->copy(s.getMatrix());
        return (*this);
    };
    MMSD_DoubleSymmetricMatrix& operator=(const LAP_DoubleFullSymmetricMatrix& s)  {
        if (mFullStorage!=null) mFullStorage->copy(s);
        else if (mPackedStorage.get()!=null) mPackedStorage->copy(s);
        return (*this);
    };
    MMSD_DoubleSymmetricMatrix& operator=(const LAP_DoublePackedSymmetricMatrix& s)  {
        if (mFullStorage.get()!=null) mFullStorage->copy(s);
        else if (mPackedStorage.get()!=null) mPackedStorage->copy(s);
        return (*this);
    };
    MMSD_DoubleSymmetricMatrix& operator=(const LAP_DoubleFullGeneralMatrix& s)  {
        if (mFullStorage.get()!=null) mFullStorage->copy(s);
        else if (mPackedStorage.get()!=null) mPackedStorage->copy(s);
        return (*this);
    };
    
 
    
    // ----------------
    // MATRIX SETTINGS
    // ----------------
public:
    virtual void copy(const LAP_DoubleMatrix& x);
    
    void copy(const MMSD_DoubleSymmetricMatrix& x);
    
    void copy(SP::MMSD_DoubleSymmetricMatrix x) {
        if (x.get()!=null) copy(*x.get());
    }

    void copy(const LAP_DoubleFullSymmetricMatrix& x);
    
    void copy(const LAP_DoublePackedSymmetricMatrix& x);

    void copy(const LAP_DoubleFullGeneralMatrix& x);
    
    inline const LAP_DoubleVector&  getValues() const {
        return getMatrix().getStorage().getValues();
    };
    inline void setValues(SP::MMSD_DoubleVector values) {
        getMatrix().setValues(values);
    };
   
    
    inline int getValuesNumber() const {
        return getMatrix().getValuesNumber();
    } 

    inline const LAP_DoubleSymmetricMatrix& getMatrix() const {
        if (mFullStorage.get()!=null) {
            return (*mFullStorage.get());
        } 
        return (*mPackedStorage.get());
        
    };
    
    inline LAP_DoubleSymmetricMatrix& getMatrix()  {
        if (mFullStorage.get()!=null) {
            return (*mFullStorage.get());
        } 
        return (*mPackedStorage.get());
        
        
    };

    inline size_t getSize() const {
        return getMatrix().getSize(LAP_DoubleMatrix::ROW);
    }
    void setStorage(const tFlag& storageType);
    inline tFlag getStorage() const {
        if (mFullStorage.get()!=null) return FULL_STORAGE;
        else return PACKED_STORAGE;
    }
    
    virtual void setSize(const size_t& n) {
        if (mFullStorage.get()!=null) mFullStorage->setSize(n);
        else if (mPackedStorage.get()!=null) mPackedStorage->setSize(n);
        
    }
    virtual void setSize(const size_t& n,const size_t& p) {
        if (mFullStorage.get()!=null) mFullStorage->setSize(n,p);
        else if (mPackedStorage.get()!=null) mPackedStorage->setSize(n,p);
    }

   
    virtual void getColumn(const size_t& j,LAP_DoubleVector& v) const {
        if (mFullStorage.get()!=null) mFullStorage->getColumn(j,v);
        else if (mPackedStorage.get()!=null) mPackedStorage->getColumn(j,v);
    }
    virtual void getColumn(const size_t& j,SP::LAP_DoubleVector v) const {
        if (mFullStorage.get()!=null) mFullStorage->getColumn(j,v);
        else if (mPackedStorage.get()!=null) mPackedStorage->getColumn(j,v);
    };

    inline size_t getRowsNumber() const {
        if (mFullStorage.get()!=null) return mFullStorage->getRowsNumber();
        else if (mPackedStorage.get()!=null) return mPackedStorage->getRowsNumber();
        return 0;
    };
    inline size_t getSize(const tFlag& d) const {
        if (mFullStorage.get()!=null) return mFullStorage->getSize(d);
        else if (mPackedStorage.get()!=null) return mPackedStorage->getSize(d);
        return 0;
    };
    inline size_t getColumnsNumber() const {
        if (mFullStorage.get() !=null) return mFullStorage->getColumnsNumber();
        else if (mPackedStorage.get()!=null) return mPackedStorage->getColumnsNumber();
        return 0;
    };
    
    virtual void set(const size_t& i,
                     const size_t& j,
                     const double& v) {
        (*this)(i,j)=v;
    };
   
    
    virtual double get(const size_t& i,
                       const size_t& j) const {
        return (*this)(i,j);
    }
    
    virtual void add(const size_t& i,
                     const size_t& j,
                     const double& v) {
       (*this)(i,j)+=v;
    }

    inline void setValues(const size_t& n,const double* values) {
        if (mFullStorage.get()!=null) mFullStorage->setValues(n,values);
        else if (mPackedStorage.get()!=null) mPackedStorage->setValues(n,values);
    }
    
    // ----------------
    // vector product
    // ----------------

public:
  
      
    void product(const LAP_DoubleVector& vect,
                 LAP_DoubleVector& result) const {
        if (mFullStorage.get() !=null) mFullStorage->vectorProduct(vect,result);
        else if (mPackedStorage.get()!=null) mPackedStorage->vectorProduct(vect,result);
    };    
    void product(const LAP_ConstDoubleVector& vect,
                 LAP_DoubleVector& result) const {
        tLVectorIndex n=vect.getSize();
        result.setSize(n);
        if (mFullStorage.get() !=null)
            mFullStorage->vectorProduct(false,
                                        n,vect.getIncrement(),&vect(0),1,0,
                                        n,result.getIncrement(),&result(0));
        else if (mPackedStorage.get()!=null)
            mPackedStorage->vectorProduct(false,
                                          n,vect.getIncrement(),&vect(0),1,0,
                                          n,result.getIncrement(),&result(0));
    }
    
    void product(SPC::LAP_DoubleVector vect,
                 SP::LAP_DoubleVector result) const {
        if (mFullStorage.get() !=null) mFullStorage->vectorProduct(*vect.get(),*result.get());
        else if (mPackedStorage.get()!=null) mPackedStorage->vectorProduct(*vect.get(),*result.get());
    };
   
    virtual void product(const LAP_DoubleVector& vect,
                         const lapack_real& alpha,
                         const lapack_real& beta,
                         LAP_DoubleVector& result) const {
        if (mFullStorage.get() !=null) mFullStorage->vectorProduct(vect,alpha,beta,result);
        else if (mPackedStorage.get()!=null) mPackedStorage->vectorProduct(vect,alpha,beta,result);
    };


    virtual void rankSymmetricProduct(const lapack_real& alpha,const LAP_DoubleVector& vect) {
        if (mFullStorage.get() !=null) {
            mFullStorage->dotProduct(alpha,vect);
        }
        else if (mPackedStorage.get()!=null) mPackedStorage->dotProduct(alpha,vect);
    }
  
    // ----------------
    // matrix product
    // ----------------
    tBoolean product(const MMSD_DoubleSymmetricMatrix& B,
                     LAP_DoubleFullGeneralMatrix& C) const {
        C.setSize(getMatrix().getRowsNumber(),B.getMatrix().getColumnsNumber());
        if (mFullStorage.get() !=null) {
            const LAP_DoubleFullSymmetricMatrix* fB=dynamic_cast<const LAP_DoubleFullSymmetricMatrix*>(&B);
            if (fB!=null) {
                mFullStorage->matrixProduct(true,*mFullStorage,*fB,1.,0.,C);
                return true;
            }
            const LAP_DoublePackedSymmetricMatrix* pB=dynamic_cast<const LAP_DoublePackedSymmetricMatrix*>(&B);
            if (pB!=null) {
                pB->matrixProduct(false,*mFullStorage,1.,0.,C);
                return true;
            }
            return false;
        }
        else if (mPackedStorage.get()!=null) {
            const LAP_DoubleFullSymmetricMatrix* fB=dynamic_cast<const LAP_DoubleFullSymmetricMatrix*>(&B);
            if (fB!=null) {
                mPackedStorage->matrixProduct(true,*fB,1.,0.,C);
                return true;
            }
            const LAP_DoublePackedSymmetricMatrix* pB=dynamic_cast<const LAP_DoublePackedSymmetricMatrix*>(&B);
            if (pB!=null) {
                mPackedStorage->matrixProduct(*pB,C);
                return true;
            }
            return false;
        }
        return false;
        
    };
    void product(const LAP_DoubleFullGeneralMatrix& B,
                 LAP_DoubleFullGeneralMatrix& C) const {
        
        if (mFullStorage.get()!=null) mFullStorage->matrixProduct(B,C);
        else if (mPackedStorage.get()!=null) mPackedStorage->matrixProduct(B,C);
    };
    virtual void  product(const LAP_DoubleFullGeneralMatrix& B,
                          const lapack_real& alpha,
                          const lapack_real& beta,
                          LAP_DoubleFullGeneralMatrix& C) const {
        if (mFullStorage.get() !=null) mFullStorage->matrixProduct(B,alpha,beta,C);
        else  if (mPackedStorage.get()!=null) mPackedStorage->matrixProduct(true,B,alpha,beta,C);
    };

    virtual  tBoolean product(const tBoolean& leftSide,
                              const LAP_DoubleFullGeneralMatrix& B,
                              const lapack_real& alpha,
                              const lapack_real& beta,
                              LAP_DoubleFullGeneralMatrix& C) const {
        if (mFullStorage.get() !=null) {
            mFullStorage->matrixProduct(leftSide,B,alpha,beta,C);
            return true;
        }
        else if (mPackedStorage.get()!=null) {
            mPackedStorage->matrixProduct(leftSide,B,alpha,beta,C);
            return true;
        }
        return false;
        
    }
   
    static  void product(const tBoolean& leftSide,
                         const MMSD_DoubleSymmetricMatrix& A,
                         const LAP_DoubleFullGeneralMatrix& B,
                         const lapack_real& alpha,
                         const lapack_real& beta,
                         LAP_DoubleFullGeneralMatrix& C) {
        const LAP_DoubleFullSymmetricMatrix *fsy=dynamic_cast<const LAP_DoubleFullSymmetricMatrix *>(&A.getMatrix());
        const LAP_DoublePackedSymmetricMatrix *psy=dynamic_cast<const LAP_DoublePackedSymmetricMatrix *>(&A.getMatrix());
        
        if (fsy!=null) LAP_DoubleFullSymmetricMatrix::matrixProduct(leftSide,*fsy,B,alpha,beta,C);
        if (psy!=null) LAP_DoublePackedSymmetricMatrix::matrixProduct(leftSide,*psy,B,alpha,beta,C);
    };
   
    

  
    // ----------------------
    // eigen values & vectors
    // ------------------------
public:
    virtual tBoolean computeEigenValues(LAP_DoubleVector& U) const {
        if (mFullStorage.get()!=null) {
            LAP_DoubleFullSymmetricMatrix C(*mFullStorage);
            return LAP_DoubleFullSymmetricMatrix::computeEigenValues(C,U);
        } else if (mPackedStorage.get()!=null) {
            LAP_DoublePackedSymmetricMatrix C(*mPackedStorage);
            return LAP_DoublePackedSymmetricMatrix::computeEigenValues(C,U);
        }
        return false;
    };


  
public:
  

    
    virtual tBoolean computeEigenValueDecomposition(LAP_DoubleVector& U,
                                                    LAP_DoubleFullGeneralMatrix& VR) const {
        if (mFullStorage.get()!=null) {
            return mFullStorage->computeEigenValueDecomposition(U,VR);
        } else if (mPackedStorage.get()!=null) {
            return mPackedStorage->computeEigenValueDecomposition(U,VR);
        }
        return false;
    };
  
    
    static tBoolean computeEigenValueDecomposition(MMSD_DoubleSymmetricMatrix& A,
                                                   LAP_DoubleVector& U,
                                                   LAP_DoubleFullGeneralMatrix& VR) {
        LAP_DoublePackedSymmetricMatrix *psy=dynamic_cast<LAP_DoublePackedSymmetricMatrix *>(&A.getMatrix());
        if (psy!=null) {
            return LAP_DoublePackedSymmetricMatrix::computeEigenValueDecomposition(*psy,U,VR);
        }
        LAP_DoubleFullSymmetricMatrix *fsy=dynamic_cast<LAP_DoubleFullSymmetricMatrix *>(&A.getMatrix());
        if (fsy!=null) {
            return LAP_DoubleFullSymmetricMatrix::computeEigenValueDecomposition(*fsy,U,VR);
        }
        return false;
    };
    

    
    tBoolean computeEigenValueDecomposition(LAP_DoubleVector& U,
                                            LAP_DoubleFullGeneralMatrix& VR,
                                            const int& il,const int& iu) const {
        if (mFullStorage.get()!=null) {
            LAP_DoubleFullSymmetricMatrix C(*mFullStorage);
            return mFullStorage->computeEigenValueDecomposition(C,U,VR,il,iu);
        } else if (mPackedStorage.get()!=null) {
            LAP_DoublePackedSymmetricMatrix C(*mPackedStorage);
            return mPackedStorage->computeEigenValueDecomposition(C,U,VR,il,iu);
        }
        return false;
    };
    
    static tBoolean computeEigenValueDecomposition(MMSD_DoubleSymmetricMatrix&A,
                                                   LAP_DoubleVector& U,
                                                   LAP_DoubleFullGeneralMatrix& Z,
                                                   const lapack_int &il,const lapack_int& ul) {
        LAP_DoublePackedSymmetricMatrix *psy=dynamic_cast<LAP_DoublePackedSymmetricMatrix *>(&A.getMatrix());
        if (psy!=null) {
            return LAP_DoublePackedSymmetricMatrix::computeEigenValueDecomposition(*psy,U,Z,il,ul);
        }
        LAP_DoubleFullSymmetricMatrix *fsy=dynamic_cast<LAP_DoubleFullSymmetricMatrix *>(&A.getMatrix());
        if (fsy!=null) {
            return LAP_DoubleFullSymmetricMatrix::computeEigenValueDecomposition(*fsy,U,Z,il,ul);
        }
        return false;
    };
    
  
    
    
    
    // ------------------------
    // operation on matrix
    // ------------------------
public:
    virtual tReal norm2() const {
        if (mFullStorage.get()!=null) return mFullStorage->norm2();
        else if (mPackedStorage.get()!=null) return mPackedStorage->norm2();
        return 0;
    }

    virtual tReal norm2(LAP_DoubleVector& v) const {
        if (mFullStorage.get()!=null) return mFullStorage->norm2(v);
        else if (mPackedStorage.get()!=null) return mPackedStorage->norm2(v);
        return 0;
    };

  
    
    virtual tReal sum(const tFlag& d,LAP_DoubleVector& s) const {
        if (mFullStorage.get()!=null) return mFullStorage->sum(d,s);
        else if (mPackedStorage.get()!=null) return mPackedStorage->sum(d,s);
        return 0;
        
    };
      
    virtual tReal sum(const tFlag& d,const int& index) const {
        if (mFullStorage.get()!=null) return mFullStorage->sum(d,index);
        else if (mPackedStorage.get()!=null) return mPackedStorage->sum(d,index);
        return 0;
        
    };
    
    virtual tReal trace() const {
        if (mFullStorage.get()!=null) return mFullStorage->trace();
        else if (mPackedStorage.get()!=null) return mPackedStorage->trace();
        return 0;
    };

    // -----------------
    // Orthogonalization
    // ------------------
public:
    void orthogonalTransformation(const int& n,
                                  const double& errorMax,
                                  MMSD_DoubleSymmetricMatrix& H);
private:
    void FStepForFluryGautshiAlgorithm(const MMSD_DoubleSymmetricMatrix& A,
                                       const int& n,
                                       const double& errorMax,
                                       LAP_DoubleFullSymmetricMatrix& H);
    void GStepForFluryGautshiAlgorithm(const MMSD_DoubleFullMatrix& T,
                                       const int& n,
                                       const double& errorMax,
                                       MMSD_DoubleFullMatrix& Q);
    
    // -----------------------
    // STRING Representation
    // ----------------------
public:
    virtual tString toString() const {
        if (mFullStorage.get()!=null) return mFullStorage->toString();
        else if (mPackedStorage.get()!=null) return mPackedStorage->toString();
        return "";
    };
  
};
#endif