LAP_DoubleFullUpperMatrix.h

Go to the documentation of this file.

#ifndef LAP_DoubleFullUpperMatrix_H
#define LAP_DoubleFullUpperMatrix_H

//inheritance classes
#include "LAP_DoubleUpperMatrix.h"

//vector classes
#include "LAP_DoubleVector.h"

//association classes
#include "LAP_DoubleFullStorage.h"

//c_lapack routines
#include "dfulluppermatrix_functions.h"
#include "dvector_functions.h"

//friend classes
DEFINE_SPTR(LAP_DoubleFullGeneralMatrix);


DEFINE_SPTR(LAP_DoubleFullUpperMatrix);
class LAP_DoubleFullUpperMatrix: public virtual LAP_DoubleUpperMatrix  {
    
    SP_OBJECT(LAP_DoubleFullUpperMatrix);
    
private:
    // ATTRIBUTES
    
    
    // ASSOCIATIONS
    
    // METHOD
  
  
public:
    // CONSTRUCTORS
    LAP_DoubleFullUpperMatrix();
    LAP_DoubleFullUpperMatrix(const LAP_DoubleFullUpperMatrix& mat);
  
    // DESTRUCTORS
    virtual ~LAP_DoubleFullUpperMatrix();

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

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

    inline static SP::LAP_DoubleFullUpperMatrix New(const tLVectorIndex& n) {
        SP::LAP_DoubleFullUpperMatrix p=New();
        p->setSize(n);
        return p;
    };
    inline static SP::LAP_DoubleFullUpperMatrix Identity(const tLVectorIndex& nRows) {
        SP::LAP_DoubleFullUpperMatrix p=New(nRows);
        
        double *vs=&(*p.get())[0];
        for (tLVectorIndex i=0;i<nRows-1;i++) {
            *vs=1;
            vs+=nRows+1;
            
        }
        *vs=1;
        return p;
    }
    
    virtual SP::LAP_DoubleMatrix NewInstance() const {
        return New();
    }
   
    //-----------
    // OPERATORS
    // ----------
  
    virtual double& operator()(const tLVectorIndex& i,const tLVectorIndex& j) {
        return (j<i)? getNullValue():getFullStorage()(i,j);
    };
    
    virtual const double& operator()(const tLVectorIndex& i,const tLVectorIndex& j) const {
        return (j<i)? getNullValue():getFullStorage()(i,j);
    };
 
   
  
    //copy operator
    //=============
    
    LAP_DoubleFullUpperMatrix& operator=(const double& s)  {
        getStorage().init(s);
        resetLowerSubMatrix();
        return (*this);
    };
    LAP_DoubleFullUpperMatrix& operator=(const LAP_DoubleFullUpperMatrix& s)  {
        copy(s);
        return (*this);
    };
 
    
    //scalar operators
    //================
    
    inline LAP_DoubleFullUpperMatrix& operator^=(const lapack_real&  s) {
        getFullStorage()^=s;
        return (*this);
    }
    
    inline LAP_DoubleFullUpperMatrix& operator*=(const lapack_real&  s) {
        getFullStorage()*=s;
        return (*this);
    }
    inline LAP_DoubleFullUpperMatrix& operator/=(const lapack_real&  s) {
        getFullStorage()/=s;
        return (*this);
    }
 
    inline LAP_DoubleFullUpperMatrix& operator+=(const lapack_real&  s) {
        getFullStorage()+=s;
        return (*this);
    }
  
  
    inline LAP_DoubleFullUpperMatrix& operator-=(const lapack_real&  s) {
        getFullStorage()-=s;
        return (*this);
    }
    
    //pack matrix operator
    //=====================
    
    inline LAP_DoubleFullUpperMatrix& operator*=(const LAP_DoubleFullUpperMatrix& s) {
        getFullStorage()*=s.getFullStorage();
        return (*this);
    }
    inline LAP_DoubleFullUpperMatrix& operator/=(const LAP_DoubleFullUpperMatrix& s) {
        getFullStorage()/=s.getFullStorage();
        return (*this);
    }
    
    inline LAP_DoubleFullUpperMatrix& operator-=(const LAP_DoubleFullUpperMatrix& s) {
        getFullStorage()-=s.getFullStorage();
        return (*this);
    }
    inline LAP_DoubleFullUpperMatrix& operator+=(const LAP_DoubleFullUpperMatrix& s) {
        getFullStorage()+=s.getFullStorage();
        return (*this);
    }
    
    // ----------------
    // MATRIX SETTINGS
    // ----------------
private:
    void resetLowerSubMatrix();


public:
    virtual tBoolean copyLower(const LAP_DoubleFullGeneralMatrix& f);
    

    //set the values from vectors
    //===========================
    
    virtual void setSize(const tLVectorIndex& n,const tLVectorIndex& p) {

       
        //set the number of the matrix
        LAP_DoubleUpperMatrix::setSize(n,p);
        
        //set the storage 
        getFullStorage().setValuesNumber(n,n*p);
    
    }
    virtual void setSize(const tLVectorIndex & n) {
        setSize(n,n);
    }
    virtual void setValues(const tLVectorIndex& n,
                           const tLVectorIndex& p,
                           SP::LAP_DoubleVector v) {
        
        
        tLVectorIndex m=v->getSize();
        if (n*p!=m)
            LAP_DoubleUpperMatrix::setSize(n,m/p);
        else
            LAP_DoubleUpperMatrix::setSize(n,p);

        //set the vector values
        getFullStorage().setValuesPointer(n,v);
    };
    virtual void setValues(SP::LAP_DoubleVector v) {
       
        tLVectorIndex m=v->getSize();
        tLVectorIndex nRows=(sqrt(m)+0.5);
        LAP_DoubleUpperMatrix::setSize(nRows,nRows);
        
        //set the vector values
        getFullStorage().setValuesPointer(nRows,v);
    };
    virtual void setValues(const tLVectorIndex& n,
                           const tLVectorIndex& p,
                           const double* values) {
        //set the size
        setSize(n,p);
        
        //set the storage
        getFullStorage().setValues(n,n*p,values);
    }
    virtual void setValues(const tLVectorIndex& n,
                           const double* values) {
    
        tLVectorIndex nRows=(tLVectorIndex) (sqrt(n)+0.5);
        LAP_DoubleUpperMatrix::setSize(nRows,nRows);
        //set the storage
        getFullStorage().setValues(nRows,n,values);
    }
   

    //======
private:
    inline LAP_DoubleFullStorage& getFullStorage() {
        return *dynamic_cast<LAP_DoubleFullStorage*>(&getStorage());
    };
    inline const LAP_DoubleFullStorage& getFullStorage() const {
        return *dynamic_cast<const LAP_DoubleFullStorage*>(&getStorage());
    };
    
public:
    inline tLVectorIncrement getRowIncrement() const {
        return getFullStorage().getRowIncrement();
    }
    inline tLVectorIndex getLeadingDimension() const {
        return getFullStorage().getLeadingDimension();
    }
    virtual tLVectorIndex getRowsNumber() const {
        return getFullStorage().getViewedRowsNumber();
    };
    virtual tLVectorIndex getColumnsNumber() const {
        return getFullStorage().getViewedColumnsNumber();
    };
    virtual void getColumn(const tLVectorIndex& j,LAP_DoubleVector& v) const;
    virtual void getColumn(const tLVectorIndex& j,SP::LAP_DoubleVector v) const {
        ASSERT_IN(v.get()!=null);
        getColumn(j,*v.get());
    };

   
  
    
    virtual void getRow(const tLVectorIndex& i,LAP_DoubleVector& v) const;
    
    virtual void getRow(const tLVectorIndex& i,SP::LAP_DoubleVector v) const {
        ASSERT_IN(v.get()!=null);
        getRow(i,*v.get());
    };
    
   
    

  
    // ------------------------
    // operation on matrix
    // ------------------------
public:
    
    
    virtual tReal norm2() const;
    
    virtual tReal norm2(LAP_DoubleVector& v) const;
    
 
   
    
    virtual tReal sum(const tFlag& d,LAP_DoubleVector& s) const;
    
    virtual tReal sum(const tFlag& d,const tLVectorIndex& index) const;
    
    virtual tReal max(const tFlag& d,const tLVectorIndex& index) const;
  
    void indexMin(tLVectorIndex& i,tLVectorIndex& j) const;

    
   
    // ----------------
    // vector product
    // ----------------
    
public:
    virtual void vectorProduct(const LAP_DoubleVector& X,LAP_DoubleVector& Y) const {
        vectorProduct(false,X,1.,0.,Y);
    }
    virtual void vectorProduct(const tBoolean& isTrans,const LAP_DoubleVector& X,LAP_DoubleVector& Y) const {
        vectorProduct(isTrans,X,1.,0.,Y);
    }
    void vectorProduct(const tBoolean& isTrans,
                       const LAP_DoubleVector& X,
                       const lapack_real& alpha,
                       const lapack_real& beta,
                       LAP_DoubleVector& Y) const {
        //get the dimenson of Matrix
        tLVectorIndex nRows=getRowsNumber();
        tLVectorIndex nCols=getColumnsNumber();
        if (isTrans) std::swap(nRows,nCols);
        
        tLVectorIndex nX=X.getSize();
        
        //verify the dimension of X
        if (nX!=nCols) {
            throw LAP_Exception("math/linalg/core","LAP_DoubleFullUpperMatrix::product()","incompatible size of X vector");
        }
       
        if (beta==0) {
            Y.setSize(nRows);
        }
        tLVectorIndex nY=Y.getSize();
        if (nY!=nRows) {
            throw LAP_Exception("math/linalg/core","LAP_DoubleFullUpperMatrix::product()","incompatible size of Y vector");
        }
      
        vectorProduct(isTrans,
                      nX,X.getIncrement(),&X(0),
                      alpha,
                      beta,
                      nY,Y.getIncrement(),&Y(0));
    };

    void vectorProduct(const tBoolean& isTrans,
                       const LAP_ConstDoubleVector& X,
                       const lapack_real& alpha,
                       const lapack_real& beta,
                       LAP_DoubleVector& Y) const {
        //get the dimenson of Matrix
        tLVectorIndex nRows=getRowsNumber();
        tLVectorIndex nCols=getColumnsNumber();
        if (isTrans) std::swap(nRows,nCols);
        
        tLVectorIndex nX=X.getSize();
        
        //verify the dimension of X
        if (nX!=nCols) {
            throw LAP_Exception("math/linalg/core","LAP_DoubleFullUpperMatrix::product()","incompatible size of X vector");
        }
       
        if (beta==0) {
            Y.setSize(nRows);
        }
        tLVectorIndex nY=Y.getSize();
        if (nY!=nRows) {
            throw LAP_Exception("math/linalg/core","LAP_DoubleFullUpperMatrix::product()","incompatible size of Y vector");
        }
      
        vectorProduct(isTrans,
                      nX,X.getIncrement(),&X(0),
                      alpha,
                      beta,
                      nY,Y.getIncrement(),&Y(0));
    };
    virtual void vectorProduct(const tBoolean& isTrans,LAP_DoubleVector& X) const {
        vectorProduct(isTrans,1,X);
    }
    /* \brief compute X:=alpha.op(This).X
     * @param isTrans : if (isTrans) op(A)=A else op(A)=tA
     * @param alpha: alpha value
     * @param X: X 
     */
    inline void vectorProduct(const tBoolean& isTrans,const double& alpha,LAP_DoubleVector& X) const {
        if (alpha==0) X=0;
        else {
            tLVectorIndex nRows=getRowsNumber();
            tLVectorIndex nCols=getColumnsNumber();
            if (nRows==nCols) {
                //X:=op(This).X
                vectorProduct(isTrans,nRows,X.getIncrement(),&X(0));
                //X*=alpha
                if (alpha!=1) X*=alpha;
            }
            else {
                LAP_DoubleVector Y;
                Y.setSize(nRows);
                vectorProduct(isTrans,
                              X.getSize(),X.getIncrement(),&X(0),
                              alpha,0.,
                              Y.getSize(),Y.getIncrement(),&Y(0));
                X.copy(Y);
            }
        }
        
    };
   
   
    
    
public:
   
    inline void vectorProduct(const tBoolean& isTrans,const tLVectorIndex& nX,const tLVectorIncrement& incX,double* X) const {
        if (getRowIncrement()!=1) {
            cout <<"WARNING:row increment matrix ignored in lapack functions. It is set to 1. \n";
        }
        tLVectorIndex nRows=getRowsNumber();
        tLVectorIndex nCols=getColumnsNumber();
        if (nRows==0) return;
        
        if (nRows!=nCols) {
            throw LAP_Exception("math/linalg/core","LAP_DoubleFullUpperMatrix::product()",
                                "X=A.X needs to have a squared matrix");
        }
        if (nCols!=nX) {
            throw LAP_Exception("math/linal/core","LAP_DoubleFullUpperMatrix::product","incompatible size of X vector");
        }
        
        DoubleFullUpperMatrixVectorProduct(isTrans,
                                           nRows,
                                           getLeadingDimension(),
                                           &(*this)(0,0),
                                           nX,incX,X);
    
    }
    void vectorProduct(const tBoolean& isTrans,
                       const tLVectorIndex& nX,const tLVectorIncrement& incX,const double *x,
                       const double& alpha,const double& beta,
                       const tLVectorIndex& nY,const tLVectorIncrement& incY,double *y) const;
    
    
public:
         
    virtual void rankProduct(const tBoolean& isTransU,
                             const LAP_DoubleVector& Diag,
                             LAP_DoubleVector& X) const {
        
       
        tLVectorIndex nX=X.getSize();
        if (getRowIncrement()!=1) {
            cout <<"WARNING:row increment matrix ignored in lapack functions. It is set to 1. \n";
        }
        if (nX!=getRowsNumber()) {
            throw LAP_Exception("math/linal/core"," LAP_DoubleFullUpperMatrix:rankDProduct","incompatible dimension of vector");
        }
        if (nX!=getColumnsNumber()) {
            throw LAP_Exception("math/linal/core"," LAP_DoubleFullUpperMatrix:rankDProduct","incompatible dimension of vector");
        }
        if (nX!=Diag.getSize()) {
            throw LAP_Exception("math/linal/core"," LAP_DoubleFullUpperMatrix:rankDProduct","incompatible dimension of vectors");
        }
        DoubleFullUpperMatrixVectorProduct(isTransU,
                                           getRowsNumber(),
                                           getLeadingDimension(),
                                           &(*this)(0,0),
                                           X.getSize(),X.getIncrement(),&X(0));
        
        X*=Diag;
        
        DoubleFullUpperMatrixVectorProduct(!isTransU,
                                           getRowsNumber(),
                                           getLeadingDimension(),
                                           &(*this)(0,0),
                                           X.getSize(),X.getIncrement(),&X(0));
    }
    
    virtual void rankProduct(const tBoolean& isTransU,LAP_DoubleVector& X) const{
        
        tLVectorIndex nX=X.getSize();
        
        if (getRowIncrement()!=1) {
            cout <<"WARNING:row increment matrix ignored in lapack functions. It is set to 1. \n";
        }
        if (nX!=getRowsNumber()) {
            throw LAP_Exception("math/linal/core"," LAP_DoubleFullUpperMatrix:rankProduct","incompatible dimension of vector");
        }
        if (nX!=getColumnsNumber()) {
            throw LAP_Exception("math/linal/core"," LAP_DoubleFullUpperMatrix:rankProduct","incompatible dimension of vector");
        }
        DoubleFullUpperMatrixVectorProduct(isTransU,
                                           getRowsNumber(),
                                           getLeadingDimension(),
                                           &(*this)(0,0),
                                           X.getSize(),X.getIncrement(),&X(0));
        DoubleFullUpperMatrixVectorProduct(!isTransU,
                                           getRowsNumber(),
                                           getLeadingDimension(),
                                           &(*this)(0,0),
                                           X.getSize(),X.getIncrement(),&X(0));
        
    }
    
    
  
    
    
public:

  
    // ----------------
    // matrix product
    // ----------------
    virtual SP::LAP_DoubleMatrix matrixProduct(const LAP_DoubleMatrix& B) const;
    

    // ----------------
    tBoolean matrixProduct(const LAP_DoubleMatrix& B,
                           LAP_DoubleFullGeneralMatrix& C) const;
    
    
    
    //product with full matrix
    //-------------------------
    
    void matrixProduct(const LAP_DoubleFullGeneralMatrix& B,
                       LAP_DoubleFullGeneralMatrix& C) const {
        
        matrixProduct(false,false,B,1.,C);
    };
     
    void matrixProduct(const LAP_DoubleFullUpperMatrix& B,
                       LAP_DoubleFullGeneralMatrix& C) const;
    
        
    void matrixProduct(const tBoolean& isLeft,
                       const tBoolean& isTrans,
                       const LAP_DoubleFullGeneralMatrix& B,
                       const double& alpha,
                       LAP_DoubleFullGeneralMatrix& C) const;
    

    virtual void matrixProduct(const tBoolean& isLeft,
                               const tBoolean& isTrans,
                               const double& alpha,
                               LAP_DoubleFullGeneralMatrix& A) const {
        matrixProduct(isLeft,isTrans,*this,alpha,A);
    }
    virtual void matrixProduct(const tBoolean& isTrans,
                               const double& alpha,
                               LAP_DoubleFullGeneralMatrix& A) const {
        matrixProduct(true,isTrans,*this,alpha,A);
    }
  
    
    
    
    static void matrixProduct(const tBoolean& isLeft,
                              const tBoolean& isTransA,
                              const LAP_DoubleFullUpperMatrix& A,
                              const lapack_real& alpha,
                              LAP_DoubleFullGeneralMatrix& B);
    
    
    
  
    // ----------------------
    // eigen values & vectors
    // ------------------------
public:
    
    
    static tBoolean computeEigenValues(const LAP_DoubleFullUpperMatrix&A,
                                       LAP_DoubleVector& U) {
        tLVectorIndex i,n=A.getRowsNumber();
        const tLVectorIncrement &incA=A.getRowIncrement();
        tLVectorIndex ldd=incA+A.getLeadingDimension();
        
        U.setSize(n);
        double *U_i=&U(0);
        const tLVectorIncrement &incU=U.getIncrement();
        const double *A_ii=&A(0,0);
        for (i=0;i<n;i++) {
            *U_i=*A_ii;
            U_i+=incU;
            A_ii+=ldd;
        }
        return true;
    }
    
    virtual tBoolean computeEigenValues(LAP_DoubleVector& U) const {
        return computeEigenValues(*this,U);
    };    
   
    virtual tBoolean computeEigenValueDecomposition(LAP_DoubleVector& U);
    
     
    virtual tBoolean computeEigenValueDecomposition(LAP_DoubleVector& U,LAP_DoubleMatrix& VR) const;
    virtual tBoolean computeEigenValueDecomposition(LAP_DoubleVector& U,LAP_DoubleMatrix& VR,LAP_DoubleMatrix& VL) const;
    

    static tBoolean computeEigenValueDecomposition(LAP_DoubleFullUpperMatrix&A,
                                                   LAP_DoubleVector& U,
                                                   LAP_DoubleFullGeneralMatrix& VR);
    
    static tBoolean computeEigenValueDecomposition(LAP_DoubleFullUpperMatrix&A,
                                                   LAP_DoubleVector& U,
                                                   LAP_DoubleFullGeneralMatrix& VR,
                                                   LAP_DoubleFullGeneralMatrix& VL);



  
    // -----------------------
    // linear system solver
    // -----------------------
    void solve(const tBoolean& isTrans,
               LAP_DoubleVector& B) const {
        tLVectorIndex nB=B.getSize();
        if (nB!=0)
            solve(isTrans,nB,B.getIncrement(),&B(0));
    }
    
    virtual void solve(const tBoolean& isTrans,
                       const tLVectorIndex& nB,const tLVectorIncrement& incB,double *B) const {
        tLVectorIndex nRows=getRowsNumber();
        tLVectorIndex nCols=getColumnsNumber();
        if (nRows==0) return;
        if (nRows!=nB) throw LAP_Exception("math/linal/core"," LAP_DoubleFullUpperMatrix:solve","incompatible dimension of vector");
        if (nRows!=nCols)
            throw LAP_Exception("math/linal/core"," LAP_DoubleFullUpperMatrix:solve","only be used with squared matrix");


        DoubleFullUpperMatrixVectorSolve(isTrans,nRows,getLeadingDimension(),&(*this)(0,0),nB,incB,B);
    }
    
    virtual tBoolean solve(const tBoolean& isLeft,
                           const tBoolean& isTrans,
                           const double &alpha,
                           LAP_DoubleFullGeneralMatrix& X) const {
        return solve(isLeft,isTrans,*this,alpha,X);
    }

    virtual tBoolean solve(const tBoolean& isTrans,
                           LAP_DoubleFullGeneralMatrix& X) const {
        return solve(true,isTrans,*this,1.,X);
    }
    
    static tBoolean solve(const tBoolean& isLeft,
                          const tBoolean& isTransA,
                          const LAP_DoubleFullUpperMatrix& A,
                          const double& alpha,
                          LAP_DoubleFullGeneralMatrix& B);
    
    tBoolean inverse() {
        return inverse(*this);
    };
    
public:
    
    inline static tBoolean inverse(LAP_DoubleFullUpperMatrix& A) {
        
        if ( A.getRowIncrement() != 1 )  
            cout<< "math/linal/core "<<"LAP_DoubleFullUpperMatrix::inverseMatrix() "<<"matrix is  non-contiguous."<<endl;
        
        
        return DoubleFullUpperMatrixInverse(A.getRowsNumber(),A.getLeadingDimension(),&A(0,0));
        
    }
    

  
};
#endif