LAP_DoubleBandedSymmetricMatrix.h

Go to the documentation of this file.

#ifndef LAP_DoubleBandedSymmetricMatrix_H
#define LAP_DoubleBandedSymmetricMatrix_H

#include "LAP_DoubleSymmetricMatrix.h"
#include "LAP_DoubleBandedStorage.h"

#include "LAP_DoubleFullGeneralMatrix.h"
#include "LAP_DoubleBandedUpperMatrix.h"

#include "LAP_DoubleVector.h"
#include "LAP_IntegerVector.h"

#include "dbandedsymmatrix_functions.h"




DEFINE_SPTR(LAP_DoubleBandedSymmetricMatrix);
class LAP_DoubleBandedSymmetricMatrix: public LAP_DoubleSymmetricMatrix  {
    
    SP_OBJECT(LAP_DoubleBandedSymmetricMatrix);
    
private:
    // ATTRIBUTES

    
    // ASSOCIATIONS
  
    
    
    // METHOD
  
  
public:
    // CONSTRUCTORS
    LAP_DoubleBandedSymmetricMatrix();
    LAP_DoubleBandedSymmetricMatrix(const LAP_DoubleBandedSymmetricMatrix& mat);
  
    // DESTRUCTORS
    virtual ~LAP_DoubleBandedSymmetricMatrix();

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

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

    inline static SP::LAP_DoubleBandedSymmetricMatrix New(const tLVectorIndex& n) {
        SP::LAP_DoubleBandedSymmetricMatrix p=New();
        p->setSize(n);
        return p;
    };
    virtual SP::LAP_DoubleMatrix NewInstance() const {
        return New();
    }
   
    //-----------
    // OPERATORS
    // ----------
   
    virtual double& operator()(const tLVectorIndex& i,const tLVectorIndex& j) {
        return (i<=j)?getBandedStorage()(i,j):getBandedStorage()(j,i);
    };
    
    virtual const double& operator()(const tLVectorIndex& i,const tLVectorIndex& j) const {
        return (i<=j)?getBandedStorage()(i,j):getBandedStorage()(j,i);
    };
    
  
 
    
    //copy operator
    //=============
    
    LAP_DoubleBandedSymmetricMatrix& operator=(const double& s)  {
        getStorage().init(s);
        return (*this);
    };
    LAP_DoubleBandedSymmetricMatrix& operator=(const LAP_DoubleMatrix& s)  {
        copy(s);
        return (*this);
    };
 
    
    //scalar operators
    //================
    
    inline LAP_DoubleBandedSymmetricMatrix& operator^=(const lapack_real&  s) {
        getBandedStorage()^=s;
        return (*this);
    }
    
    inline LAP_DoubleBandedSymmetricMatrix& operator*=(const lapack_real&  s) {
        getBandedStorage()*=s;
        return (*this);
    }
    inline LAP_DoubleBandedSymmetricMatrix& operator/=(const lapack_real&  s) {
        getBandedStorage()/=s;
        return (*this);
    }
 
    inline LAP_DoubleBandedSymmetricMatrix& operator+=(const lapack_real&  s) {
        getBandedStorage()+=s;
        return (*this);
    }
  
  
    inline LAP_DoubleBandedSymmetricMatrix& operator-=(const lapack_real&  s) {
        getBandedStorage()-=s;
        return (*this);
    }
    
    //pack matrix operator
    //=====================
    
    inline LAP_DoubleBandedSymmetricMatrix& operator*=(const LAP_DoubleBandedSymmetricMatrix& s) {
        getBandedStorage()*=s.getBandedStorage();
        return (*this);
    }
    inline LAP_DoubleBandedSymmetricMatrix& operator/=(const LAP_DoubleBandedSymmetricMatrix& s) {
        getBandedStorage()/=s.getBandedStorage();
        return (*this);
    }
    
    inline LAP_DoubleBandedSymmetricMatrix& operator-=(const LAP_DoubleBandedSymmetricMatrix& s) {
        getBandedStorage()-=s.getBandedStorage();
        return (*this);
    }
    inline LAP_DoubleBandedSymmetricMatrix& operator+=(const LAP_DoubleBandedSymmetricMatrix& s) {
        getBandedStorage()+=s.getBandedStorage();
        return (*this);
    }
    
    
    // ----------------
    // MATRIX SETTINGS
    // ----------------

public:
    
   
    
    //dimension setting
    //=================
    
    virtual void setSize(const tLVectorIndex & nRows,const tLVectorIndex& nCols) {
        LAP_DoubleSymmetricMatrix::setSize(nRows,nCols);
        LAP_DoubleBandedStorage &s=getBandedStorage();
        tLVectorIndex nBandsNumber=s.getBandsNumber();
        if (nBandsNumber<=0) {
            if (nRows>=nCols) s.setBandsNumber(nCols,nRows-nCols+1);
            else s.setBandsNumber(nCols,1);
            nBandsNumber=s.getBandsNumber();
        }
        getStorage().setValuesNumber( ((tLVectorIndex) nCols) * nBandsNumber);
    }
    
    virtual void setSize(const tLVectorIndex & n) {
        LAP_DoubleSymmetricMatrix::setSize(n,n);
        LAP_DoubleBandedStorage &s=getBandedStorage();
        tLVectorIndex nBandsNumber=s.getBandsNumber();
        if (nBandsNumber<=0) {
            s.setBandsNumber(n,1);
            nBandsNumber=s.getBandsNumber();
        }
        getStorage().setValuesNumber( ((tLVectorIndex) n) * nBandsNumber);
    }
    
    inline void setBandsNumber(const tLVectorIndex& n) {
        getBandedStorage().setBandsNumber(n,1);
    }
        

    //set the values from vectors
    //===========================


    virtual void setValues(SP::LAP_DoubleVector values) {
        ASSERT_IN(values.get()!=null);

        //compute the numbr of columns
        tLVectorIndex n=values->getSize();
        n/=getBandsNumber();
        LAP_DoubleSymmetricMatrix::setSize(n,n);

        //set the values pointer
        getStorage().setValuesPointer(values);
        
    }
    virtual void setValues(const tLVectorIndex& n,const double* values) {
        tLVectorIndex nRows=n/getBandsNumber();
        LAP_DoubleSymmetricMatrix::setSize(nRows,nRows);
        
        //set the values pointer
        getStorage().setValues(n,values);
    }

    inline tBoolean compile() {
        return getBandedStorage().compile();
    }
    
    //======

public:
    inline LAP_DoubleBandedStorage& getBandedStorage() {
        return *dynamic_cast<LAP_DoubleBandedStorage*>(&getStorage());
    };
    inline const LAP_DoubleBandedStorage& getBandedStorage() const {
        return *dynamic_cast<const LAP_DoubleBandedStorage*>(&getStorage());
    };
public:
    inline  const tLVectorIndex&  getBandsNumber() const {
        return getBandedStorage().getBandsNumber();
    }
    inline  const tLVectorIndex&  getUpperBandsNumber() const {
        return getBandedStorage().getUpperBandsNumber();
    }
    
    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);
        LAP_DoubleBandedSymmetricMatrix::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);
        LAP_DoubleBandedSymmetricMatrix::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;
   

    // Linear methods
    // ==============
    
  
    
   
    // ----------------
    // vector product
    // ----------------
    
public:

    virtual void vectorProduct(const LAP_DoubleVector& X,LAP_DoubleVector& Y) {
        LAP_DoubleMatrix::vectorProduct(X,Y);
    }
    
    virtual void vectorProduct(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 nX=X.getSize();
        
        //verify the dimension of X
        if (nX!=nRows) {
            throw LAP_Exception("math/linalg/core","LAP_DoubleBandedSymmetricMatrix::product()","incompatible size of X vector");
        }
        
        //set the dimension of Y is beta=0
        if (beta==0) {
            Y.setSize(nRows);
        }
        tLVectorIndex nY=Y.getSize();
        if (nY!=nRows) {
            throw LAP_Exception("math/linalg/core","LAP_DoubleBandedSymmetricMatrix::product()","incompatible size of Y vector");
        }
        //make the product
        vectorProduct(nX,X.getIncrement(),&X(0),
                      alpha,
                      beta,
                      nY,Y.getIncrement(),&Y(0));
    };
  
    virtual void vectorProduct(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 nX=X.getSize();
        
        //verify the dimension of X
        if (nX!=nRows) {
            throw LAP_Exception("math/linalg/core","LAP_DoubleBandedSymmetricMatrix::product()","incompatible size of X vector");
        }
        
        //set the dimension of Y is beta=0
        if (beta==0) {
            Y.setSize(nRows);
        }
        tLVectorIndex nY=Y.getSize();
        if (nY!=nRows) {
            throw LAP_Exception("math/linalg/core","LAP_DoubleBandedSymmetricMatrix::product()","incompatible size of Y vector");
        }
        //make the product
        vectorProduct(nX,X.getIncrement(),&X(0),
                      alpha,
                      beta,
                      nY,Y.getIncrement(),&Y(0));
    };
    
    
public:
    virtual void vectorProduct(const tBoolean& isTrans,
                               const tLVectorIndex& nX,const tLVectorIncrement& incX,const double *X,
                               const lapack_real& alpha, const lapack_real& beta,
                               const tLVectorIndex& nY,const tLVectorIncrement& incY,double *Y) const {
        vectorProduct(nX,incX,X,alpha,beta,nY,incY,Y);
    }
    
    virtual void vectorProduct(const tLVectorIndex& nx,const tLVectorIncrement& incx,const double* x,
                               const lapack_real& alpha,
                               const lapack_real& beta,
                               const tLVectorIndex& ny,const tLVectorIncrement&incy,double* y) const {
        
        tLVectorIndex nRows=getRowsNumber();
        tLVectorIndex nBands=getBandsNumber();
        if (nRows==0) return;
        if (nRows!=nx) {
            throw LAP_Exception("math/linal/core","LAP_DoubleBandedSymmetricMatrix::product","incompatible size of X vector");
        }
        if (nRows!=ny) {
            throw LAP_Exception("math/linal/core","LAP_DoubleBandedSymmetricMatrix::product","incompatible size of Y vector");
      }
        if (x==y) {
            throw LAP_Exception("math/linal/core","LAP_DoubleBandedSymmetricMatrix::product","x & y must have different pointer");
        }
        
        DoubleBandedSymmetricMatrixVectorProduct(nx,incx,x,
                                               nRows,nBands,                                          
                                               &(*this)[0],
                                               alpha,beta,
                                               ny,incy,y);
    }
    
    inline void rankSymmetricProduct(const LAP_DoubleFullGeneralMatrix& A) {
        rankSymmetricProduct(1,0,A.getRowsNumber(),A.getColumnsNumber(),
                             A.getRowIncrement(),A.getLeadingDimension(),false,&A(0,0));
    }
    inline void rankSymmetricProduct(const double alpha,const tBoolean& isTransA,const LAP_DoubleFullGeneralMatrix& A) {
        rankSymmetricProduct(alpha,0,A.getRowsNumber(),A.getColumnsNumber(),
                             A.getRowIncrement(),A.getLeadingDimension(),isTransA,&A(0,0));
    }
       
    virtual void rankSymmetricProduct(const double& alpha,const double& beta,
                                      const tLVectorIndex& nARows,const tLVectorIndex& nACols,const tLVectorIncrement& incA,const tLVectorIndex& ldA,
                                      const tBoolean& isTransA,const double* A);

    // ----------------
    // matrix product
    // ----------------
    virtual SP::LAP_DoubleMatrix  matrixProduct(const LAP_DoubleMatrix& B) const {
        SP::LAP_DoubleFullGeneralMatrix C=LAP_DoubleFullGeneralMatrix::New();
       
        const LAP_DoubleFullGeneralMatrix* dgeM=dynamic_cast<const LAP_DoubleFullGeneralMatrix*>(&B);
        if (dgeM!=null) {
            matrixProduct(*dgeM,*C.get());
            return C;
        }
        return C;
    };
   
    void matrixProduct(const LAP_DoubleFullGeneralMatrix& B,
                       LAP_DoubleFullGeneralMatrix& C) const {
        
        C.setSize(getRowsNumber(),B.getColumnsNumber());
        matrixProduct(true,(*this),B,1.0,0,C);
    };
   

    virtual void  matrixProduct(const tBoolean& leftSide,
                                const LAP_DoubleFullGeneralMatrix& B,
                                const lapack_real& alpha,
                                const lapack_real& beta,
                                LAP_DoubleFullGeneralMatrix& C) const {
        
       matrixProduct(leftSide,*this,B,alpha,beta,C);
       
    }
    virtual void  matrixProduct(const tBoolean& leftSide,
                                const tLVectorIndex& nBRows,
                                const tLVectorIndex& nBCols,
                                const tLVectorIncrement& incB,
                                const tLVectorIndex& ldB,
                                const double* B,
                                const lapack_real& alpha,
                                const lapack_real& beta,
                                LAP_DoubleFullGeneralMatrix& C) const {
        ASSERT_IN(B!=&C(0,0));
        
        //verify the dimension
        tLVectorIndex N=getRowsNumber();
        
        tLVectorIndex M=nBRows;
        tLVectorIndex P=nBCols;
        if (!leftSide) std::swap(M,P);
        
        if (N!=M) throw LAP_Exception("math/linalg/core","LAP_DoubleFullSymmetricMatrix::product()","incompatible size of matrix");
        if (beta==0) {
            if (leftSide) C.setSize(N,P);
            else C.setSize(P,N);
        } else {
            if (leftSide && ((C.getRowsNumber()!=N) || (C.getColumnsNumber()!=P))) C.setSize(N,P);
            else if (!leftSide && ((C.getRowsNumber()!=P) || (C.getColumnsNumber()!=N))) C.setSize(P,N);
        }
      
        if (nBRows!=nBCols) throw LAP_Exception("math/linalg/core","LAP_DoubleFullSymmetricMatrix::product()","incompatible size of matrix");
        if (beta==0) {
            if (leftSide) C.setSize(N,P);
            else C.setSize(P,N);
        } else {
            if (leftSide && ((C.getRowsNumber()!=N) || (C.getColumnsNumber()!=P))) C.setSize(N,P);
            else if (!leftSide && ((C.getRowsNumber()!=P) || (C.getColumnsNumber()!=N))) C.setSize(P,N);
        }
        DoubleBandedSymmetricMatrixMatrixProduct(leftSide,
                                                 getRowsNumber(),getBandsNumber(),&(*this)[0],//S
                                                 nBRows,nBCols,incB,ldB,B,//B
                                                 alpha,beta,
                                                 C.getRowsNumber(),C.getColumnsNumber(),
                                                 C.getRowIncrement(),
                                                 C.getLeadingDimension(),&C(0,0));
    }
   
    static  void matrixProduct(const tBoolean& leftSide,
                               const LAP_DoubleBandedSymmetricMatrix& S,
                               const LAP_DoubleFullGeneralMatrix& B,
                               const lapack_real& alpha,
                               const lapack_real& beta,
                               LAP_DoubleFullGeneralMatrix& C) {
        ASSERT_IN(&B!=&C);
        
        //verify the dimension
        tLVectorIndex N=S.getRowsNumber();
        tLVectorIndex M=B.getRowsNumber();
        tLVectorIndex P=B.getColumnsNumber();
        if (!leftSide) std::swap(M,P);
        
        if (N!=M) throw LAP_Exception("math/linalg/core","LAP_DoubleBandedSymmetricMatrix::product()","incompatible size of matrix");
        if (beta==0) {
            if (leftSide) C.setSize(N,P);
            else C.setSize(P,N);
        } else {
            if (leftSide && ((C.getRowsNumber()!=N) || (C.getColumnsNumber()!=P))) C.setSize(N,P);
            else if (!leftSide && ((C.getRowsNumber()!=P) || (C.getColumnsNumber()!=N))) C.setSize(P,N);
        }
        DoubleBandedSymmetricMatrixMatrixProduct(leftSide,
                                                 S.getRowsNumber(),S.getBandsNumber(),&S[0],//S
                                                 B.getRowsNumber(),B.getColumnsNumber(),
                                                 B.getRowIncrement(),
                                                 B.getLeadingDimension(),&B(0,0),//B
                                                 alpha,beta,
                                                 C.getRowsNumber(),C.getColumnsNumber(),
                                                 C.getRowIncrement(),C.getLeadingDimension(),&C(0,0));
        
    }
    
  
    // ----------------------
    // eigen values & vectors
    // ------------------------
public:
    virtual tBoolean computeEigenValues(LAP_DoubleVector& U) const {
        LAP_DoubleBandedSymmetricMatrix C(*this);
        return computeEigenValues(C,U);
    };
    static tBoolean computeEigenValues(LAP_DoubleBandedSymmetricMatrix&A,LAP_DoubleVector& U) {
        tLVectorIndex nRows=A.getRowsNumber();
        tLVectorIndex nCols=A.getColumnsNumber();
        
        if (nRows>nCols)
            throw LAP_Exception("math/linal/core","LAP_DoubleBandedSymmetricMatrix::computeEigenValues","matrix columns number must be greater than its rows number");
       
        U.resetView();
        U.setSize(nRows);
        
        return DoubleBandedSymmetricMatrixEigenValues(nRows,A.getBandsNumber(),&A[0],nRows,&U(0),1,null);
        
    }

  public:
    tBoolean computeEigenValueDecomposition(LAP_DoubleVector& U,LAP_DoubleFullGeneralMatrix&W) const {
        LAP_DoubleBandedSymmetricMatrix A;
        A.copy(*this);
        return computeEigenValueDecomposition(A,U,W);  
    };

private:
    static tBoolean computeEigenValueDecomposition(LAP_DoubleBandedSymmetricMatrix&A,
                                                   LAP_DoubleVector& U,
                                                   LAP_DoubleFullGeneralMatrix&W) {
        tLVectorIndex nRows=A.getRowsNumber();;
      
        U.resetView();
        U.setSize(nRows);
        W.setSize(nRows,nRows);
        return DoubleBandedSymmetricMatrixEigenValues(nRows,A.getBandsNumber(),&A[0],nRows,&U(0),
                                                      W.getLeadingDimension(),&W(0,0));
    }

    
   
   
public:
    tBoolean computeEigenValueDecomposition(LAP_DoubleVector& U,
                                            LAP_DoubleFullGeneralMatrix& VR,
                                            const int& il,const int& iu) const {
        LAP_DoubleBandedSymmetricMatrix C(*this);
        LAP_DoubleFullGeneralMatrix Q;
        return computeEigenValueDecomposition(C,U,VR,Q,il,iu);
    };
    
   
    static tBoolean computeEigenValueDecomposition(LAP_DoubleBandedSymmetricMatrix&A,
                                                   LAP_DoubleVector& U,
                                                   LAP_DoubleFullGeneralMatrix& Z,
                                                   LAP_DoubleFullGeneralMatrix& Q,
                                                   const lapack_int &il,const lapack_int& ul) {
        tLVectorIndex nRows=A.getRowsNumber();
        
        //eigen values
        int nEigens=ul-il+1;
        U.resetView();
        U.setSize(nEigens);
        
        //eigen vectors
        Z.setSize(nRows,nEigens);

        //orthogonlaisation matrix
        Q.setSize(nRows,nRows);
        
        return DoubleBandedSymmetricMatrixBoundedEigenValues(nRows,A.getBandsNumber(),&A[0],//A
                                                             nEigens,&U(0),
                                                             il,ul,0,0,
                                                             Q.getLeadingDimension(),&Q(0,0),
                                                             Z.getLeadingDimension(),&Z(0,0));
    }
    
   
    
private:
    static SP::LAP_DoubleBandedUpperMatrix choleskyFactorization(LAP_DoubleBandedSymmetricMatrix& A) {
       
        tLVectorIndex nRows=A.getRowsNumber();
        tBoolean succeeds=DoubleBandedSymmetricMatrixCholeskyFactorization(nRows,A.getBandsNumber(),&A[0]);
        SP::LAP_DoubleBandedUpperMatrix U=LAP_DoubleBandedUpperMatrix::New();
        U->setBandsNumber(A.getBandsNumber());
        SP::LAP_DoubleVector vs;
        A.getStorage().getValuesPointer(vs);
        if (succeeds) U->setValues(vs);
        return U;
    }
public:

    virtual  SP::LAP_DoubleUpperMatrix choleskyFactorization() {
        return LAP_DoubleBandedSymmetricMatrix::choleskyFactorization(*this);
    }
   
  
};
#endif