LAP_DoubleFullGeneralMatrix.h

Go to the documentation of this file.

#ifndef LAP_DoubleFullGeneralMatrix_H
#define LAP_DoubleFullGeneralMatrix_H

#include "LAP_DoubleGeneralMatrix.h"

#include "LAP_DoubleFullStorage.h"

#include "LAP_DoubleVector.h"
#include "LAP_ConstDoubleVector.h"

#include "dfullmatrix_functions.h"



DEFINE_SPTR(LAP_DoubleFullGeneralMatrix);

class LAP_DoubleFullGeneralMatrix: public LAP_DoubleGeneralMatrix  {
    
    SP_OBJECT(LAP_DoubleFullGeneralMatrix);
    
    
    
    // ATTRIBUTES
public:
    static const tFlag QR;
    static const tFlag QRP;

   
    
private:
    //transposed the element
    tBoolean mIsTransposed;
    
    lapack_int *mLUPivots;
    
    // ASSOCIATIONS
 
    
    // METHOD
  
  
public:
    // CONSTRUCTORS
    LAP_DoubleFullGeneralMatrix();
    LAP_DoubleFullGeneralMatrix(const tLVectorIndex& n,const tLVectorIndex& m);
    LAP_DoubleFullGeneralMatrix(const LAP_DoubleFullGeneralMatrix& mat);
  
    // DESTRUCTORS
    virtual ~LAP_DoubleFullGeneralMatrix();

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

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

    inline static SP::LAP_DoubleFullGeneralMatrix New(const int& n,const int& m) {
        SP::LAP_DoubleFullGeneralMatrix p(new LAP_DoubleFullGeneralMatrix(n,m),
                               LAP_DoubleFullGeneralMatrix::Delete());
        p->setThis(p);
        return p;
    };
    inline static SP::LAP_DoubleFullGeneralMatrix New(const tLVectorIndex& n,const tLVectorIndex & m) {
        SP::LAP_DoubleFullGeneralMatrix p(new LAP_DoubleFullGeneralMatrix(n,m),
                               LAP_DoubleFullGeneralMatrix::Delete());
        p->setThis(p);
        return p;
    };

    inline static SP::LAP_DoubleFullGeneralMatrix Identity(const tLVectorIndex& nRows) {
        SP::LAP_DoubleFullGeneralMatrix Id=New(nRows,nRows);

        //get the values
        double *vs=&(*Id.get())[0];
        //set diagonal to 1
        for (tLVectorIndex i=0;i<nRows-1;i++) {
            *vs=1;
            vs+=nRows+1;
            
        }
        *vs=1;
        return Id;
    }
    virtual SP::LAP_DoubleMatrix NewInstance() const {
        return New();
    }
  
    //-----------
    // OPERATORS
    // ----------
    virtual double& operator()(const tLVectorIndex& i,const tLVectorIndex& j) {
        return getFullStorage()(i,j);
    };
    
    virtual const double& operator()(const tLVectorIndex& i,const tLVectorIndex& j) const {
        return getFullStorage()(i,j);
    };

  
   

    //copy operator
    //=============
   
    
    LAP_DoubleFullGeneralMatrix& operator=(const double& s)  {
        getStorage().init(s);
        return (*this);
    };
    LAP_DoubleFullGeneralMatrix& operator=(const LAP_DoubleFullGeneralMatrix& s)  {
        copy(s);
        return (*this);
    };
   
    
    //scalar operators
    //================
    
    inline LAP_DoubleFullGeneralMatrix& operator^=(const lapack_real&  s) {
        getFullStorage()^=s;
        return (*this);
    }
    
    inline LAP_DoubleFullGeneralMatrix& operator*=(const lapack_real&  s) {
        getFullStorage()*=s;
        return (*this);
    }
    inline LAP_DoubleFullGeneralMatrix& operator/=(const lapack_real&  s) {
        getFullStorage()/=s;
        return (*this);
    }
 
    inline LAP_DoubleFullGeneralMatrix& operator+=(const lapack_real&  s) {
        getFullStorage()+=s;
        return (*this);
    }
  
  
    inline LAP_DoubleFullGeneralMatrix& operator-=(const lapack_real&  s) {
        getFullStorage()-=s;
        return (*this);
    }
    
    //pack matrix operator
    //=====================
    
    inline LAP_DoubleFullGeneralMatrix& operator*=(const LAP_DoubleFullGeneralMatrix& s) {
        getFullStorage()*=s.getFullStorage();
        return (*this);
    }
    inline LAP_DoubleFullGeneralMatrix& operator/=(const LAP_DoubleFullGeneralMatrix& s) {
        getFullStorage()/=s.getFullStorage();
        return (*this);
    }
    
    inline LAP_DoubleFullGeneralMatrix& operator-=(const LAP_DoubleFullGeneralMatrix& s) {
        getFullStorage()-=s.getFullStorage();
        return (*this);
    }
    inline LAP_DoubleFullGeneralMatrix& operator+=(const LAP_DoubleFullGeneralMatrix& s) {
        getFullStorage()+=s.getFullStorage();
        return (*this);
    }
    
    
    //dimension setting
    //=================
    
    virtual void setSize(const tLVectorIndex& n,const tLVectorIndex& p) {
        //set the number of the matrix
        LAP_DoubleGeneralMatrix::setSize(n,p);
        
        //set the storage 
        tLVectorIndex dim=((tLVectorIndex)n)*((tLVectorIndex)p);

        //set the size of the storage
        getFullStorage().setValuesNumber(n,dim);
    
    }
    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_DoubleGeneralMatrix::setSize(n,m/p);
        else
            LAP_DoubleGeneralMatrix::setSize(n,p);

        //set the vector values
        getFullStorage().setValuesPointer(n,v);
    };
    virtual void setValues(SP::LAP_DoubleVector v) {
        tLVectorIndex n=v->getSize();
        tLVectorIndex nRows=(tLVectorIndex) (sqrt(n)+0.5);
        LAP_DoubleGeneralMatrix::setSize(nRows,nRows);

        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_DoubleGeneralMatrix::setSize(nRows,nRows);

        getFullStorage().setValues(nRows,n,values);
    }
    

 
    inline void concat(const double& alpha,const LAP_DoubleFullGeneralMatrix& X,
                       const double& beta, const LAP_DoubleFullGeneralMatrix& Y) {
        
        //verify the number of rows must be the same
        if (X.getSize(LAP_Matrix::ROW)!=Y.getSize(LAP_Matrix::ROW))
            throw LAP_Exception("math/linalg/core",
                                "LAP_DoubleFullGeneralMatrix::concat()",
                                "incompatible dimenion of concat matrices");
        
        // the new matrix is of size nRows x (nColsOfX+nColsOfY)
        setSize(X.getSize(ROW),X.getSize(COL)+Y.getSize(COL));
        //merge the values
        getFullStorage().concat(alpha,X.getFullStorage(),
                                beta,Y.getFullStorage());
        
        
    }
    
    // GET
    // ======
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();
    }

    void transpose();
    
    inline tBoolean isTransposed() const {
      return mIsTransposed;
    }
    virtual tLVectorIndex getRowsNumber() const {
        return getFullStorage().getViewedRowsNumber();
    };
    virtual tLVectorIndex getColumnsNumber() const {
        return getFullStorage().getViewedColumnsNumber();
    };
    void add(const double& alpha,
             const LAP_DoubleFullGeneralMatrix& x);
   
    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());
    };

   
   
    inline void getColumnByReference(const tLVectorIndex& j,LAP_DoubleVector& v) {
        getFullStorage().getColumnByReference(j,v);
    }
    inline void getColumnByReference(const tLVectorIndex& j,LAP_ConstDoubleVector& v) const {
        getFullStorage().getColumnByReference(j,v);
    }
    
    virtual void getRow(const tLVectorIndex& i,SP::LAP_DoubleVector v) const {
        if (v.get()==null) throw LAP_Exception("math/linalg,core",
                                               "getRow()",
                                               "can not store the row in a null shared pointer"); 
        getRow(i,*v.get());
    };
    virtual void getRow(const tLVectorIndex& i,LAP_DoubleVector& v) const;
    
    inline void getRowByReference(const tLVectorIndex& i,LAP_DoubleVector& v) {
        getFullStorage().getRowByReference(i,v);
    }
    inline void getRowByReference(const tLVectorIndex& i,LAP_ConstDoubleVector& v) const {
        getFullStorage().getRowByReference(i,v);
    }
    
    
    //column manipulations
    // ====================
    
    void removeColumn(const tLVectorIndex& index);

    void swapColumns(const tLVectorIndex& i,const tLVectorIndex& j) {
        getFullStorage().swapColumns(i,j);
    }
    
    void reverseColumns() {
        getFullStorage().reverseColumns();
    }
    
    
    
   
    // =============
    // Tools methods
    // =============
    
    
    inline tBoolean symmetricTest() const {
        tLVectorIndex k,n=getRowsNumber();
        tLVectorIndex j,p=getColumnsNumber();
        if (n!=p) return false;
        
        tLVectorIncrement inc=getRowIncrement();
        tLVectorIndex ld=getLeadingDimension();
        
        const double *vc_jk,*vc_j=&(*this)(0,0);
        const double *vr_jk,*vr_j=vc_j;
        for (j=0;j<p;j++) {
            vr_jk=vr_j;
            vc_jk=vc_j;
            for (k=0;k<j;k++) {
                if (fabs(*vc_jk-(*vr_jk))>LAP_DoubleMatrixStorage::ZERO_EPSILON) return false;
                vr_jk+=ld;
                vc_jk+=inc;
            }
            vc_j+=ld;
            vr_j+=inc;
        }
        return true;
    };
    inline tBoolean upperTest() const {
        tLVectorIndex i,j,p=getColumnsNumber();
        tLVectorIndex n=getRowsNumber();
        const double *v_ij,*v_jj=&(*this)(0,0);
        const tLVectorIncrement& inc=getRowIncrement();
        tLVectorIndex ld=getLeadingDimension();
        tLVectorIndex ldd=ld+inc;
        for (j=0;j<p;j++) {
            v_ij=v_jj+inc;
            for (i=j+1;i<n;i++) {
                if (fabs(*v_ij)>LAP_DoubleMatrixStorage::ZERO_EPSILON) return false;
                v_ij+=inc;
            }
            v_jj+=ldd;
        }
        return true;
    };
    inline tBoolean lowerTest() const {
        tLVectorIndex i,n=getRowsNumber();
        tLVectorIndex j,p=getColumnsNumber();
        const tLVectorIncrement& inc=getRowIncrement();
        tLVectorIndex ld=getLeadingDimension();
        const double *v_ij,*v_j=&(*this)(0,0);
        for (j=0;j<p;j++) {
            v_ij=v_j+(j+1)*inc;
            for (i=j+1;i<n;i++) {
                if (fabs(*v_ij)>LAP_DoubleMatrixStorage::ZERO_EPSILON) return false;
                v_ij+=inc;
            }
            v_j+=ld;
        }
        return true;
    };
    
    
    // ------------------------
    // 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{
        LAP_DoubleMatrix::vectorProduct(false,X,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!=nRows) {
            throw LAP_Exception("math/linalg/core","LAP_DoubleFullGeneralMatrix::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_DoubleFullGeneralMatrix::product()","incompatible size of Y vector");
        }
        //make the product
        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);

        //get the diemsnion of the vector X
        tLVectorIndex nX=X.getSize();
        
        //verify the dimension of X
        if (nX!=nRows) {
            throw LAP_Exception("math/linalg/core","LAP_DoubleFullGeneralMatrix::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_DoubleFullGeneralMatrix::product()","incompatible size of Y vector");
        }
        //make the product
        vectorProduct(isTrans,
                      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 {
        
        if (getRowIncrement()!=1) {
            cout <<"WARNING:row increment matrix ignored in lapack functions. It is set to 1. \n";
        }
        tLVectorIndex nRows=getRowsNumber();
        tLVectorIndex nCols=getColumnsNumber();
        tLVectorIndex ld= getLeadingDimension();
        
        if (nRows*nCols==0) return;
        if (nCols!=nx) {
            throw LAP_Exception("math/linal/core","LAP_DoubleFullGeneralMatrix::product","incompatible size of X vector");
        }
        if (nRows!=ny) {
            throw LAP_Exception("math/linal/core","LAP_DoubleFullGeneralMatrix::product","incompatible size of Y vector");
        }
        
        if (x==y) {
            throw LAP_Exception("math/linal/core","LAP_DoubleFullGeneralMatrix::product","x & y vector must not have same pointer");
        }
        
        DoubleFullMatrixVectorProduct(nx,incx,x,
                                      isTrans,
                                      nRows,nCols,ld,
                                      &(*this)(0,0),
                                      alpha,beta,
                                      ny,incy,y);
    }
    
    // ----------------
    // matrix product
    // ----------------
public:
    virtual SP::LAP_DoubleMatrix matrixProduct(const LAP_DoubleMatrix& B) const {
        SP::LAP_DoubleFullGeneralMatrix  C=LAP_DoubleFullGeneralMatrix::New();
        matrixProduct(B,*C.get());
        return C;
    }
    
    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(true,false,B,1.,C);
    }
  
  
 
    void matrixProduct(const tBoolean& isLeftSide,
                      const tBoolean& isTrans,
                      const LAP_DoubleFullGeneralMatrix& B,
                      const double& alpha,
                      LAP_DoubleFullGeneralMatrix& C) const {
        tLVectorIndex nBCols=B.getColumnsNumber();
        tLVectorIndex nBRows=B.getRowsNumber();
      
        
        tLVectorIndex nCols=getColumnsNumber();
        tLVectorIndex nRows=getRowsNumber();
        if (isTrans) std::swap(nCols,nRows);
        
        tLVectorIndex nCRows=nRows;
        tLVectorIndex nCCols=nBCols;

        if (!isLeftSide) {
            nCRows=nBRows;
            nCCols=nCols;
            if (nBCols!=nRows) throw LAP_Exception("math/linalg/core","LAPD_DoubleFullGeneralMatrix::product","incompatible size of matrices");
        } else {
            if (nCols!=nBRows) throw LAP_Exception("math/linalg/core","LAPD_DoubleFullGeneralMatrix::product","incompatible size of matrices");
        }
        
        
        C.setSize(nCRows,nCCols);
        if (isLeftSide) {
            matrixProduct(isTrans,(*this),false,B,alpha,0.,C);
        } else {
            matrixProduct(false,B,isTrans,(*this),alpha,0.,C);
        }
    };
  

    
    static inline void matrixProduct(const tBoolean& isTrA,
                                     const LAP_DoubleFullGeneralMatrix& A,
                                     const tBoolean& isTrB,
                                     const LAP_DoubleFullGeneralMatrix& B,
                                     const lapack_real& alpha,
                                     const lapack_real& beta,
                                     LAP_DoubleFullGeneralMatrix& C) {
        if ((A.getRowIncrement()!=1) || (B.getRowIncrement()!=1)) {
            cout << "WARNING : Row increment in matric are ignored in lapack routines...set to 1 \n";
        }
        tLVectorIndex nARows=A.getRowsNumber();
        tLVectorIndex nACols=A.getColumnsNumber();
        tLVectorIndex nBRows=B.getRowsNumber();
        tLVectorIndex nBCols=B.getColumnsNumber();
        if (nARows*nACols==0) return;
        if (nBRows*nBCols==0) return;
        
        matrixProduct(isTrA,nARows,nACols,A.getLeadingDimension(),&A(0,0),
                      isTrB,nBRows,nBCols,B.getLeadingDimension(),&B(0,0),
                      alpha,beta,C);
        
    } 
    static void matrixProduct(const tBoolean& isTrA,
                              const tLVectorIndex& nRowsA,const tLVectorIndex& nColsA,const tLVectorIndex& ldA,const double *A,
                              const tBoolean& isTrB,
                              const tLVectorIndex& nRowsB,const tLVectorIndex& nColsB,const tLVectorIndex& ldB,const double *B,
                              const lapack_real& alpha,
                              const lapack_real& beta,
                              LAP_DoubleFullGeneralMatrix& C);
  
  
   
    // ----------------------
    // eigen values & vectors
    // ------------------------
public:
   
    
    static tBoolean computeEigenValues(LAP_DoubleFullGeneralMatrix&A,
                                       LAP_DoubleVector& U);
    
    virtual tBoolean computeEigenValues(LAP_DoubleVector& U) const {
        LAP_DoubleFullGeneralMatrix A(*this);
        return computeEigenValues(A,U);
    };    
   
    virtual tBoolean computeEigenValueDecomposition(LAP_DoubleVector& U) {
         LAP_DoubleFullGeneralMatrix VR;
         tBoolean succeeds=computeEigenValueDecomposition(*this,U,VR);
         copy(VR);
         return succeeds;
    };
    
     
    virtual tBoolean computeEigenValueDecomposition(LAP_DoubleVector& U,LAP_DoubleMatrix& VR) const {
        LAP_DoubleFullGeneralMatrix C(*this);
        LAP_DoubleFullGeneralMatrix *fVR=dynamic_cast< LAP_DoubleFullGeneralMatrix *>(&VR);
        if (fVR!=null) 
            return computeEigenValueDecomposition(C,U,*fVR);
        else
            return false;
    }


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


    // ------------------------
    // Orthogonalization
    // ------------------------
    inline void orthogonalize() {
        orthogonalize(QRP);
    }
  
    inline void orthogonalize(const tFlag& m) {
        if (m==QR) {
            QRFactorization(*this,true);
        } else if (m==QRP) {
            QRFactorizationWithPivotingColumns(*this,true);
        }
    }

private:
    static tBoolean QRFactorization(LAP_DoubleFullGeneralMatrix& A,const tBoolean& onlyQ);
    static tBoolean QRFactorizationWithPivotingColumns(LAP_DoubleFullGeneralMatrix& A,const tBoolean& onlyQ);
    
    // ------------------------
    // factorization
    // ------------------------

  
    inline int LUFactorization() {
        if (mLUPivots!=null) delete[] mLUPivots;
        tLVectorIndex nRows=getRowsNumber();
        tLVectorIndex nCols=getColumnsNumber();
        int N=(nRows>nCols)?nRows:nCols;
        mLUPivots=new lapack_int[N];
        return LUFactorization(*this,mLUPivots);
    }
    
    static int LUFactorization(LAP_DoubleFullGeneralMatrix&A,lapack_int *ipiv ) {
        
        tLVectorIndex nRows=A.getRowsNumber();
        tLVectorIndex nCols=A.getColumnsNumber();
        tLVectorIndex ld=A.getLeadingDimension();
        return DoubleFullMatrixLUFactorization(nRows,nCols,ld,&A(0,0),ipiv);
    }
   
    
  
    
    
    // -----------------------
    // linear system solvers
    // -----------------------
public:
    virtual tBoolean inverse();
    
    
public:
    tBoolean solve(const tBoolean& isTr,LAP_DoubleFullGeneralMatrix& X) const {
        LAP_DoubleFullGeneralMatrix Ac;
        Ac.copy(*this);
        tLVectorIndex nRows=getRowsNumber();
        tLVectorIndex nCols=getColumnsNumber();
        int N=(nRows>nCols)?nRows:nCols;
        lapack_int *ipiv=new lapack_int[N];
        int info=LUFactorization(Ac,ipiv);
        if (info!=0) {
            delete[] ipiv;
            return false;
        }
        tBoolean succeeds=solve(Ac,ipiv,isTr,X);
        delete[] ipiv;
        return succeeds;
    }
    tBoolean LUSolve(const tBoolean& isTr,LAP_DoubleFullGeneralMatrix& X)  {
        int info=LUFactorization();
        if (info!=0) return false;
        
        return solve(*this,mLUPivots,isTr,X);
    }
    

    static tBoolean solve(const LAP_DoubleFullGeneralMatrix& A,const lapack_int* ipiv,const tBoolean& isTr,LAP_DoubleFullGeneralMatrix& X);
   
    
   
};
#endif