LAP_DoubleFullSymmetricMatrix.h

Go to the documentation of this file.

#ifndef LAP_DoubleFullSymmetricMatrix_H
#define LAP_DoubleFullSymmetricMatrix_H

//inheritance classes
#include "LAP_DoubleSymmetricMatrix.h"

//vector classes
#include "LAP_DoubleVector.h"
#include "LAP_IntegerVector.h"

//association classes
#include "LAP_DoubleFullStorage.h"

//C-lapack routines
#include "dfullsymmatrix_functions.h"

//friend classes
#include "LAP_DoubleFullGeneralMatrix.h"
#include "LAP_DoubleFullUpperMatrix.h"
class LAP_DoublePackedSymmetricMatrix;


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

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

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

    inline static SP::LAP_DoubleFullSymmetricMatrix New(const tLVectorIndex& n) {
        SP::LAP_DoubleFullSymmetricMatrix p=New();
        p->setSize(n);
        return p;
    };
    
    inline static SP::LAP_DoubleFullSymmetricMatrix Identity(const tLVectorIndex& nRows) {
        SP::LAP_DoubleFullSymmetricMatrix 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)? getFullStorage()(j,i):getFullStorage()(i,j);
    };
    
    virtual const double& operator()(const tLVectorIndex& i,const tLVectorIndex& j) const {
        return (j<i)? getFullStorage()(j,i):getFullStorage()(i,j);
    };
    
  
 
    
    //copy operator
    //=============
    
    LAP_DoubleFullSymmetricMatrix& operator=(const double& s)  {
        getStorage().init(s);
        return (*this);
    };
    LAP_DoubleFullSymmetricMatrix& operator=(const LAP_DoubleFullSymmetricMatrix& s)  {
        copy(s);
        return (*this);
    };
 
    //scalar operators
    //================
    
    inline LAP_DoubleFullSymmetricMatrix& operator^=(const lapack_real&  s) {
        getFullStorage()^=s;
        return (*this);
    }
    
    inline LAP_DoubleFullSymmetricMatrix& operator*=(const lapack_real&  s) {
        getFullStorage()*=s;
        return (*this);
    }
    inline LAP_DoubleFullSymmetricMatrix& operator/=(const lapack_real&  s) {
        getFullStorage()/=s;
        return (*this);
    }
 
    inline LAP_DoubleFullSymmetricMatrix& operator+=(const lapack_real&  s) {
        getFullStorage()+=s;
        return (*this);
    }
  
  
    inline LAP_DoubleFullSymmetricMatrix& operator-=(const lapack_real&  s) {
        getFullStorage()-=s;
        return (*this);
    }
    
    //pack matrix operator
    //=====================
    
    inline LAP_DoubleFullSymmetricMatrix& operator*=(const LAP_DoubleFullSymmetricMatrix& s) {
        getFullStorage()*=s.getFullStorage();
        return (*this);
    }
    inline LAP_DoubleFullSymmetricMatrix& operator/=(const LAP_DoubleFullSymmetricMatrix& s) {
        getFullStorage()/=s.getFullStorage();
        return (*this);
    }
    
    inline LAP_DoubleFullSymmetricMatrix& operator-=(const LAP_DoubleFullSymmetricMatrix& s) {
        getFullStorage()-=s.getFullStorage();
        return (*this);
    }
    inline LAP_DoubleFullSymmetricMatrix& operator+=(const LAP_DoubleFullSymmetricMatrix& s) {
        getFullStorage()+=s.getFullStorage();
        return (*this);
    }
    LAP_DoubleFullSymmetricMatrix& operator*=(const LAP_DoublePackedSymmetricMatrix& s);
    LAP_DoubleFullSymmetricMatrix& operator/=(const LAP_DoublePackedSymmetricMatrix& s);
    
    LAP_DoubleFullSymmetricMatrix& operator-=(const LAP_DoublePackedSymmetricMatrix& s);
    LAP_DoubleFullSymmetricMatrix& operator+=(const LAP_DoublePackedSymmetricMatrix& s);
    
    
    // ----------------
    // MATRIX SETTINGS
    // ----------------
public:
    void copyUpperToLowerSubMatrix();
public:
    

    virtual void setSize(const tLVectorIndex& n,const tLVectorIndex& p) {

        tLVectorIndex d=(n<p)?n:p;
        //set the number of the matrix
        LAP_DoubleSymmetricMatrix::setSize(d,d);
        
        //set the storage 
        tLVectorIndex dim=((tLVectorIndex)d);
        dim*=dim;
        
        getFullStorage().setValuesNumber(d,dim);
    
    }
    virtual void setSize(const tLVectorIndex & n) {
        setSize(n,n);
    }
    
    virtual void setValues(SP::LAP_DoubleVector v) {
        
        
        tLVectorIndex m=v->getSize();
        tLVectorIndex nRows=(tLVectorIndex) (sqrt(m)+0.5);
        LAP_DoubleSymmetricMatrix::setSize(nRows,nRows);
        
        //set the vector values
        getFullStorage().setValuesPointer(nRows,v);
    };
  
    virtual void setValues(const tLVectorIndex& n,
                           const double* values) {
        tLVectorIndex nRows=(tLVectorIndex) (sqrt(n)+0.5);
        LAP_DoubleSymmetricMatrix::setSize(nRows,nRows);
        //set the storage
        getFullStorage().setValues(nRows,n,values);
    }
   
    //======
public:
    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 tLVectorIndex 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());
    };

    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,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) {
        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_DoubleFullSymmetricMatrix::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_DoubleFullSymmetricMatrix::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_DoubleFullSymmetricMatrix::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_DoubleFullSymmetricMatrix::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 {
        if (getRowIncrement()!=1) {
            cout <<"WARNING:row increment matrix ignored in lapack functions. It is set to 1. \n";
        }
        tLVectorIndex nRows=getRowsNumber();
        tLVectorIndex ld= getLeadingDimension();
        if (nRows==0) return;
        if (nRows!=nx) {
            throw LAP_Exception("math/linal/core","LAP_DoubleFullSymmetricMatrix::product","incompatible size of X vector");
        }
        if (nRows!=ny) {
            throw LAP_Exception("math/linal/core","LAP_DoubleFullSymmetricMatrix::product","incompatible size of Y vector");
      }
        if (x==y) {
            throw LAP_Exception("math/linal/core","LAP_DoubleFullSymmetricMatrix::product","x & y must have different pointer");
        }
        
        DoubleFullSymmetricMatrixVectorProduct(nx,incx,x,
                                               nRows,
                                               ld,
                                               &(*this)(0,0),
                                               alpha,beta,
                                               ny,incy,y);
    }
    
    
    

public:
    virtual void dotProduct(const LAP_DoubleVector& v) {
        dotProduct(1,v);
    };
    virtual void dotProduct(const double& alpha,const LAP_DoubleVector& X) {
        //A is considered to be of size 1 x Nx : X is the first row
        tLVectorIndex nRows=X.getSize();
        setSize(nRows);
        
        DoubleFullSymmetricMatrixRankSymmetricProduct(X.getIncrement(),&X(0),//X
                                                      alpha,//alpha
                                                      nRows,getLeadingDimension(),&(*this)(0,0));
        
    }
    void rankSymmetricProduct(const LAP_DoubleFullGeneralMatrix& A) {
        rankSymmetricProduct(1.,0.,false,A);
    }
    void rankSymmetricProduct(const double& alpha,const tBoolean& isTrans,const LAP_DoubleFullGeneralMatrix& A) {
        rankSymmetricProduct(alpha,0,isTrans,A);
    }
    
    virtual void rankSymmetricProduct(const double& alpha,const double& beta,const tBoolean& isTrans,const LAP_DoubleFullGeneralMatrix& A) {
        //A is considered to be of size 1 x Nx : X is the first row
        tLVectorIndex nARows=A.getRowsNumber();
        tLVectorIndex nACols=A.getColumnsNumber();
        tLVectorIndex ldA=A.getLeadingDimension();
        const tLVectorIncrement &incA=A.getRowIncrement();
        if (nARows*nACols==0) return;
        rankSymmetricProduct(alpha,beta,nARows,nACols,incA,ldA,isTrans,&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) {
        if (alpha==0) {
            (*this)*=beta;
            return;
        }
        tLVectorIndex nRows=nARows;
        tLVectorIndex nCols=nACols;
        if (isTransA) std::swap(nRows,nCols);
        if (beta==0) {
            setSize(nRows);
        } else {
            if (getRowsNumber()!=nRows) {
                throw LAP_Exception("math/linalg/core",
                                    "LAP_DoubleFullSymmetricMatrix::rankSymmetricProduct",
                                    "matrices size incompatible");
            }
        }
        
        DoubleFullSymmetricMatrixKRankSymmetricProduct(nARows,nACols,incA,ldA,isTransA,A,
                                                       alpha,beta,
                                                       0,0,null,
                                                       nRows,getRowIncrement(),getLeadingDimension(),&(*this)(0,0));
    }
    
    void rankSymmetricProduct(const LAP_DoubleVector& D,const LAP_DoubleFullGeneralMatrix& A) {
        rankSymmetricProduct(D,1.,0.,false,A);
    }
    void rankSymmetricProduct(const LAP_DoubleVector& D,const double& alpha,const tBoolean& isTrans,const LAP_DoubleFullGeneralMatrix& A) {
        rankSymmetricProduct(D,alpha,0,isTrans,A);
    }
    
    virtual void rankSymmetricProduct(const LAP_DoubleVector& D,const double& alpha,const double& beta,const tBoolean& isTrans,const LAP_DoubleFullGeneralMatrix& A) {
        if (alpha==0) {
            (*this)*=beta;
            return;
        }
        //A is considered to be of size 1 x Nx : X is the first row
        tLVectorIndex nRows=A.getRowsNumber();
        tLVectorIndex nCols=A.getColumnsNumber();
        if (nRows*nCols==0) return;
        
        if (isTrans) std::swap(nRows,nCols);
        if (beta==0) {
            setSize(nRows);
        } else {
            if (getRowsNumber()!=nRows) {
                throw LAP_Exception("math/linalg/core",
                                    "LAP_DoubleFullSymmetricMatrix::rankSymmetricProduct",
                                    "matrices size incompatible");
            }
        }
        if (nCols!=D.getSize()) { 
                throw LAP_Exception("math/linalg/core",
                                   "LAP_DoubleFullSymmetricMatrix::rankSymmetricProduct",
                                   "matrix & diagonal have incompatible size");
        }
        DoubleFullSymmetricMatrixKRankSymmetricProduct(A.getRowsNumber(),A.getColumnsNumber(),
                                                       A.getRowIncrement(),A.getLeadingDimension(),isTrans,&A(0,0),
                                                       alpha,beta,
                                                       D.getSize(),D.getIncrement(),&D(0),
                                                       nRows,getRowIncrement(),getLeadingDimension(),&(*this)(0,0));
    }
  
   
 

    
    
    // ----------------
    // matrix product
    // ----------------
    virtual SP::LAP_DoubleMatrix  matrixProduct(const LAP_DoubleMatrix& B) const {
        SP::LAP_DoubleFullGeneralMatrix C=LAP_DoubleFullGeneralMatrix::New();
        
        const LAP_DoubleFullSymmetricMatrix* dsyM=dynamic_cast<const LAP_DoubleFullSymmetricMatrix*>(&B);
        if (dsyM!=null) {
            matrixProduct(*dsyM,*C.get());
            return C;
        }

        
        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_DoubleFullSymmetricMatrix& B,
                     LAP_DoubleFullGeneralMatrix& C) const {

        C.setSize(getRowsNumber(),B.getColumnsNumber());
        matrixProduct(true,(*this),B,1.0,0,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 LAP_DoubleFullGeneralMatrix& B,
                               const lapack_real& alpha,
                               const lapack_real& beta,
                               LAP_DoubleFullGeneralMatrix& C) const {
        matrixProduct(true,(*this),B,alpha,beta,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);
        }
        DoubleFullSymmetricMatrixMatrixProduct(leftSide,
                                               getRowsNumber(),getLeadingDimension(),&(*this)(0,0),//S
                                               nBRows,nBCols,ldB,false,B,//B
                                               alpha,beta,
                                               C.getRowsNumber(),C.getColumnsNumber(),C.getLeadingDimension(),&C(0,0));
    }
   
    static  void matrixProduct(const tBoolean& leftSide,
                               const LAP_DoubleFullSymmetricMatrix& 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_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);
        }
        DoubleFullSymmetricMatrixMatrixProduct(leftSide,
                                               S.getRowsNumber(),S.getLeadingDimension(),&S(0,0),//S
                                               B.getRowsNumber(),B.getColumnsNumber(),B.getLeadingDimension(),false,&B(0,0),//B
                                               alpha,beta,
                                               C.getRowsNumber(),C.getColumnsNumber(),C.getLeadingDimension(),&C(0,0));
        
    }
    
    static  void matrixProduct(const tBoolean& leftSide,
                               const LAP_DoubleFullSymmetricMatrix& S1,
                               const LAP_DoubleFullSymmetricMatrix& S2,
                               const lapack_real& alpha,
                               const lapack_real& beta,
                               LAP_DoubleFullGeneralMatrix& C) {  //verify the dimension
        tLVectorIndex N=S1.getRowsNumber();
        tLVectorIndex M=S2.getRowsNumber();
        tLVectorIndex P=S2.getColumnsNumber();
        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);
        }
        DoubleFullSymmetricMatrixMatrixProduct(leftSide,
                                               S1.getRowsNumber(),S1.getLeadingDimension(),&S1(0,0),//S1
                                               S2.getRowsNumber(),S2.getColumnsNumber(),S2.getRowsNumber(),true,&S2(0,0),//S2
                                               alpha,beta,
                                               C.getRowsNumber(),C.getColumnsNumber(),C.getLeadingDimension(),&C(0,0));
    }
    
    

  
    // ----------------------
    // eigen values & vectors
    // ------------------------
public:
   
    static tBoolean computeEigenValues(LAP_DoubleFullSymmetricMatrix&A,LAP_DoubleVector& U) {
        tLVectorIndex nRows=A.getRowsNumber();
        tLVectorIndex nCols=A.getColumnsNumber();
        // The leading dimension of the array A. 
        lapack_int lda = A.getLeadingDimension();
        
        if (nRows>nCols)
            throw LAP_Exception("math/linal/core","LAP_DoubleFullSymmetricMatrix::computeEigenValues","matrix columns number must be greater than its rows number");
        if ( A.getRowIncrement() != 1 )
            throw(LAP_Exception("math/linal/core","LAP_DoubleFullSymmetricMatrix::computeEigenValues","matrix is  non-contiguous."));
        
        U.resetView();
        U.setSize(nRows);
        
        return DoubleFullSymmetricMatrixEigenValues(nRows,nCols,lda,&A(0,0),nRows,&U(0),false);
        
    }

    virtual tBoolean computeEigenValues(LAP_DoubleVector& U) const {
        LAP_DoubleFullSymmetricMatrix C(*this);
        return computeEigenValues(C,U);
    };

private:
    static tBoolean computeEigenValueDecomposition(LAP_DoubleFullSymmetricMatrix&A,
                                                   LAP_DoubleVector& U) {
        tLVectorIndex nRows=A.getRowsNumber();
        tLVectorIndex nCols=A.getColumnsNumber();
        // The leading dimension of the array A. 
        lapack_int lda = A.getLeadingDimension();
        
        if (nRows>nCols)
            throw LAP_Exception("math/linal/core","LAP_DoubleFullSymmetricMatrix::computeEigenValuesDecomposition","matrix columns number must be greater than its rows number");
        if ( A.getRowIncrement() != 1 )
            throw(LAP_Exception("math/linal/core","LAP_DoubleFullSymmetricMatrix::computeEigenValuesDecomposition","matrix is  non-contiguous."));
        
        U.resetView();
        U.setSize(nRows);
        
        return DoubleFullSymmetricMatrixEigenValues(nRows,nCols,lda,&A(0,0),nRows,&U(0),true);
    }
public:
    tBoolean computeEigenValueDecomposition(LAP_DoubleVector& U) {
        return computeEigenValueDecomposition(*this,U);  
    };
    
    inline tBoolean computeEigenValueDecomposition(LAP_DoubleVector& U,LAP_DoubleFullGeneralMatrix& VR) const {
        LAP_DoubleFullSymmetricMatrix SA;
        SA.copy(*this);
        tBoolean succeeds=computeEigenValueDecomposition(SA,U);
        // SA is changed into eigen vector matrix
        // get its values
        SP::LAP_DoubleVector v;
        SA.getStorage().getValuesPointer(v);
        // transert it in VP matrix
        VR.setValues((tLVectorIndex)getRowsNumber(),(tLVectorIndex)getRowsNumber(),v);
        // SA is deleted
        return succeeds;
    }
   
   
    static tBoolean computeEigenValueDecomposition(LAP_DoubleFullSymmetricMatrix&A,
                                                   LAP_DoubleVector& U,
                                                   LAP_DoubleFullGeneralMatrix& VR) {
        // compute the eigen values U & vectors in VR
        tBoolean succeeds=A.computeEigenValueDecomposition(U);
        SP::LAP_DoubleVector v;
        A.getStorage().getValuesPointer(v);
        // transfert it in VP matrix
        VR.setValues(A.getRowsNumber(),A.getRowsNumber(),v);
        return succeeds;
    };
public:
   

    virtual tBoolean computeEigenValueDecomposition(LAP_DoubleVector& U,
                                                    LAP_DoubleMatrix& VR) const {
        LAP_DoubleFullGeneralMatrix *fVR=dynamic_cast<LAP_DoubleFullGeneralMatrix*>(&VR);
        if (fVR!=null) {
            LAP_DoubleFullSymmetricMatrix SA;
            SA.copy(*this);
            return computeEigenValueDecomposition(SA,U,*fVR);
        } else {
            return false;
        }
    };

    
    tBoolean computeEigenValueDecomposition(LAP_DoubleVector& U,
                                            LAP_DoubleFullGeneralMatrix& VR,
                                            const int& il,const int& iu) const {
        LAP_DoubleFullSymmetricMatrix C(*this);
        return computeEigenValueDecomposition(C,U,VR,il,iu);
    };
    static tBoolean computeEigenValueDecomposition(LAP_DoubleFullSymmetricMatrix&A,
                                                   LAP_DoubleVector& U,
                                                   LAP_DoubleFullGeneralMatrix& Z,
                                                   const lapack_int &il,const lapack_int& ul) {
        tLVectorIndex nRows=A.getRowsNumber();
        tLVectorIndex nCols=A.getColumnsNumber();
        // The leading dimension of the array A. 
        lapack_int lda = A.getLeadingDimension();
        
        if (nRows>nCols)
            throw LAP_Exception("math/linal/core","LAP_DoubleFullSymmetricMatrix::computeEigenDecomposition","matrix columns number must be greater than its rows number");
        if ( A.getRowIncrement() != 1 )
            throw(LAP_Exception("math/linal/core","LAP_DoubleFullSymmetricMatrix::computeEigenDecomposition","matrix is  non-contiguous."));
        
        //eigen values
        int nEigens=ul-il+1;
        U.resetView();
        U.setSize(nEigens);
        
        //eigen vectors
        Z.setSize(nRows,nEigens);
        lapack_int ldz = Z.getLeadingDimension();
        
        return DoubleFullSymmetricMatrixBoundedEigenValues(nRows,nCols,lda,&A(0,0),//A
                                                           nEigens,&U(0),
                                                           true,il,ul,0,0,ldz,&Z(0,0));
    }
    
    // -----------------------
    // linear system solver
    // -----------------------
    tBoolean inverse() {
        return inverse(*this);
    };
    
public:
    
    inline static tBoolean inverse(LAP_DoubleFullSymmetricMatrix& A) {
        if ( A.getRowIncrement() != 1 )
        throw(LAP_Exception("math/linal/core","LAP_DoubleFullSymmetricMatrix::inverseMatrix","matrix is  non-contiguous."));
        tLVectorIndex nRows=A.getRowsNumber();
        LAP_IntegerVector ipiv;
        ipiv.setSize(nRows);
        if (BunchKaufmanFactorization(A,ipiv))
            return DoubleFullSymmetricMatrixInverse(nRows,A.getLeadingDimension(),&A(0,0),&ipiv(0));
        return false;
    };
    inline static tBoolean inverse(LAP_DoubleFullSymmetricMatrix& A,const LAP_IntegerVector& ipiv) {
        if ( A.getRowIncrement() != 1 )
        throw(LAP_Exception("math/linal/core","LAP_DoubleFullSymmetricMatrix::inverseMatrix","matrix is  non-contiguous."));
        return DoubleFullSymmetricMatrixInverse(A.getRowsNumber(),A.getLeadingDimension(),&A(0,0),&ipiv(0));
    };

    // Transformation
    // ==============
private:
    static SP::LAP_DoubleFullUpperMatrix choleskyFactorization(LAP_DoubleFullSymmetricMatrix& A) {
        if ( A.getRowIncrement() != 1 )
            throw(LAP_Exception("math/linal/core","LAP_DoubleFullSymmetricMatrix::choleskyFactorization","matrix is  non-contiguous."));
        tLVectorIndex nRows=A.getRowsNumber();
        tBoolean succeeds=DoubleFullSymmetricMatrixCholeskyFactorization(nRows,A.getLeadingDimension(),&A(0,0));
        SP::LAP_DoubleFullUpperMatrix U=LAP_DoubleFullUpperMatrix::New();
        SP::LAP_DoubleVector vs;
        A.getStorage().getValuesPointer(vs);
        if (succeeds) U->setValues(vs);
        return U;
    }
    static tBoolean BunchKaufmanFactorization(LAP_DoubleFullSymmetricMatrix& A,LAP_IntegerVector& ipiv) {
        if ( A.getRowIncrement() != 1 )
            throw(LAP_Exception("math/linal/core","LAP_DoubleFullSymmetricMatrix::BunchKaufmanFactorization","matrix is  non-contiguous."));
        tLVectorIndex nRows=A.getRowsNumber();
        ipiv.setSize(nRows);
        tBoolean succeeds=DoubleFullSymmetricMatrixBunchKaufmanFactorization(nRows,A.getLeadingDimension(),&A(0,0),&ipiv(0));
        return succeeds;
       
    }
public:

    virtual  SP::LAP_DoubleUpperMatrix choleskyFactorization() {
        return LAP_DoubleFullSymmetricMatrix::choleskyFactorization(*this);
    }
    virtual  tBoolean BunchKaufmanFactorization(LAP_IntegerVector& ipiv) {
        
        return LAP_DoubleFullSymmetricMatrix::BunchKaufmanFactorization(*this,ipiv);
    }
   
  
};
#endif