pzmatrix.h

Go to the documentation of this file.

#ifndef _TMATRIXHH_
#define _TMATRIXHH_

#include "pzvec.h"
#include "TPZStream.h"
#include "pzreal.h"
#include "TPZSavable.h"
#include "Hash/TPZHash.h"
#include "pzlog.h"

#include <list>
#include <sstream>

#ifdef _AUTODIFF
#include "fad.h"
#include "tfad.h"
#endif

template<class TVar>
class TPZFMatrix;

#define CLONEDEF(A) virtual TPZMatrix<TVar>*Clone() const override { return new A(*this); }


template<class TVar>
class TPZSolver;

#ifndef USING_MKL
extern "C"{
                double ddot(int *N, double *X, int *INCX, double *Y, int *INCY);
}
#endif

enum DecomposeType {ENoDecompose, ELU, ELUPivot, ECholesky, ELDLt};

enum MatrixOutputFormat {EFormatted, EInputFormat, EMathematicaInput, EMatlabNonZeros, EMatrixMarket};

template<class TVar=REAL>
class TPZMatrix: public TPZSavable

{
public:
    
    //typedef typename TVar;
                TPZMatrix() {
                                fDecomposed = 0;
                                fDefPositive = 0;
                                fRow = 0;
                                fCol = 0;
                }
                
                TPZMatrix<TVar>(const TPZMatrix<TVar>&cp) : TPZRegisterClassId(&TPZMatrix<TVar>::ClassId), fRow(cp.fRow), fCol(cp.fCol), fDecomposed(cp.fDecomposed),fDefPositive(cp.fDefPositive)
                {
                }
                virtual ~TPZMatrix();
                
                virtual TPZMatrix<TVar>*Clone() const = 0;

        virtual int64_t MemoryFootprint() const {
                  std::cout << __PRETTY_FUNCTION__ 
                                    << ": Please, implement me! (class = " << ClassId() 
                            << std::endl;
                  //<< ") (class name = " << ClassName() << ")" << std::endl;
                  //DebugStop();
                  return 0;
                }
    
    template<class TVar2>
    void CopyFrom(TPZMatrix<TVar2> &copy)
    {
        fDecomposed = copy.IsDecomposed();
        fDefPositive = copy.IsDefPositive();
        fRow = copy.Rows();
        fCol = copy.Cols();

    }

                void AutoFill(int64_t nrow, int64_t ncol, int symmetric);
                
                virtual int VerifySymmetry(REAL tol = 1.e-13) const;
                
                virtual int Put(const int64_t row,const int64_t col,const TVar & value );
                virtual const TVar &Get(const int64_t row,const int64_t col ) const;
                
                const TVar &g(const int64_t row, const int64_t col) const {return Get(row,col);}
                
                TVar &operator() (const int64_t row,const int64_t col );
                virtual TVar &s(const int64_t row, const int64_t col);
                TVar &operator()(const int64_t row);
                
                virtual int PutVal(const int64_t /*row*/,const int64_t /*col*/,const TVar & val )
    {
        if(val != ((TVar)(0.))) DebugStop();
        return 0;
    }
    virtual const TVar &GetVal(const int64_t /*row*/, const int64_t /*col*/ ) const;
                
                virtual void Multiply(const TPZFMatrix<TVar>& A,TPZFMatrix<TVar>& res, int opt = 0) const;
                virtual void Add(const TPZMatrix<TVar>& A,TPZMatrix<TVar>& res) const;
                virtual void MultAdd(const TPZFMatrix<TVar> & x,const TPZFMatrix<TVar>& y, TPZFMatrix<TVar>& z,
                                                                                                 const TVar alpha=1., const TVar beta = 0., const int opt = 0) const;
                
                virtual void Residual(const TPZFMatrix<TVar>& x,const TPZFMatrix<TVar>& rhs, TPZFMatrix<TVar>& res ) ;
                virtual void Substract(const TPZMatrix<TVar>& A,TPZMatrix<TVar>& result) const;
                
                virtual void Identity();
                virtual void Transpose(TPZMatrix<TVar>*const T) const;
                
                int Inverse(TPZFMatrix<TVar>&Inv, DecomposeType dec);
                
                TVar MatrixNorm(int p, int64_t numiter = 2000000, REAL tol = 1.e-10) const; // -<
                
                TVar ConditionNumber(int p, int64_t numiter = 2000000, REAL tol = 1.e-10);
                
                virtual void Input(std::istream & in = std::cin );
                
                template <class TT> friend std::istream & operator >> (std::istream& in,TPZMatrix<TT>& A) ;
                
    virtual void Print(std::ostream &out) const 
    {
        std::string name = __PRETTY_FUNCTION__;
        Print(name.c_str(),out);
    }
                virtual void Print(const char *name, std::ostream &out = std::cout ,const MatrixOutputFormat form = EFormatted) const;
                
                //void PrintMath(const char *name, std::ostream &out);
                
                int64_t Rows() const;
                int64_t Cols() const;
                
                inline virtual int64_t Dim() const;
                
                
                virtual int Resize(const int64_t newRows, const int64_t newCols ) {
                                fRow = newRows;
                                fCol = newCols;
                                return 0;
                }
                
                virtual int Redim(const int64_t newRows, const int64_t newCols ) {
                                fRow = newRows;
                                fCol = newCols;
                                return 0;
                }
                
                virtual int Zero(){
                                std::cout << "WARNING! TPZMatrix<TVar>::Zero is called\n";
        DebugStop();
                                return 0; }
                
                virtual int PutSub( const int64_t sRow, const int64_t sCol, const TPZFMatrix<TVar>& Source );
                
                virtual int GetSub( const int64_t sRow, const int64_t sCol, const int64_t rowSize,
                                                                                   const int64_t colSize, TPZFMatrix<TVar>& Target ) const;
                
                virtual int AddSub(const int64_t sRow, const int64_t sCol, const TPZFMatrix<TVar>& Source );
                
                virtual int InsertSub(const int64_t sRow,const int64_t sCol,const int64_t rowSize,const int64_t colSize,
                                                                                                  const int64_t pRow,const int64_t pCol, TPZMatrix<TVar>* Target ) const;
                
                virtual int AddSub(const int64_t sRow,const int64_t sCol,const int64_t rowSize,
                                                                                   const int64_t colSize,const int64_t pRow,const int64_t pCol, TPZMatrix<TVar>* pA ) const;
                
                virtual  void AddKel(TPZFMatrix<TVar>&elmat, TPZVec<int64_t> &destinationindex);

                virtual  void AddKel(TPZFMatrix<TVar>&elmat, TPZVec<int64_t> &sourceindex,  TPZVec<int64_t> &destinationindex);

                virtual void UpdateFrom(TPZAutoPointer<TPZMatrix<TVar> > /* mat*/)    
                {
                                std::cout << "TPZMatrix<TVar>::UdateFrom is not implemented\n";
        DebugStop();
                }
                
                virtual int IsSimetric() const    { return 0; }
                inline int IsSquare() const { return fRow == fCol;}
                
                virtual void Simetrize();
                
                
                virtual int IsDefPositive() const{ return 0; }
                int IsDecomposed() const         { return fDecomposed; }
                
                void SetIsDecomposed(int val) {fDecomposed = (char) val; }
                
                virtual void SolveJacobi(int64_t & numiterations, const TPZFMatrix<TVar>& F, TPZFMatrix<TVar>& result,
                                                                                                                 TPZFMatrix<TVar>* residual, TPZFMatrix<TVar>& scratch, REAL & tol, const int FromCurrent = 0);
                
                virtual void SolveSOR(int64_t & numiterations, const TPZFMatrix<TVar>& F, TPZFMatrix<TVar>& result,
                                                                                                  TPZFMatrix<TVar>* residual,TPZFMatrix<TVar>& scratch,const REAL overrelax, REAL & tol,
                                                                                                  const int FromCurrent = 0,const int direction = 1) ;
                virtual void SolveSSOR(int64_t & numiterations,const TPZFMatrix<TVar>& F, TPZFMatrix<TVar>& result,
                                                                                                   TPZFMatrix<TVar>* residual, TPZFMatrix<TVar>& scratch, const REAL overrelax, REAL & tol,
                                                                                                   const int FromCurrent = 0) ;
                
                virtual void SolveCG(int64_t & numiterations, TPZSolver<TVar> & preconditioner,
                                                                                                 const TPZFMatrix<TVar>& F, TPZFMatrix<TVar>& result,
                                                                                                 TPZFMatrix<TVar>* residual, REAL & tol,
                                                                                                 const int FromCurrent = 0) ;
                virtual void SolveBICG(int64_t & numiterations, TPZSolver<TVar> & preconditioner,
                                                                                                   const TPZFMatrix<TVar>& F, TPZFMatrix<TVar>& result,
                                                                                                   REAL & tol) ;
                
                virtual void SolveBICGStab(int64_t & numiterations, TPZSolver<TVar> & preconditioner,
                                                                                                                   const TPZFMatrix<TVar>& F, TPZFMatrix<TVar>& result,
                                                                                                                   TPZFMatrix<TVar>* residual, REAL & tol,
                                                                                                                   const int FromCurrent = 0) ;
                
                virtual void SolveGMRES(int64_t & numiterations, TPZSolver<TVar> & preconditioner,
                                                                                                                TPZFMatrix<TVar>& H, int & numvectors,
                                                                                                                const TPZFMatrix<TVar>& F, TPZFMatrix<TVar>& result,
                                                                                                                TPZFMatrix<TVar>* residual, REAL & tol,const int FromCurrent) ;
                
                virtual void SolveIR(int64_t & numiterations, TPZSolver<TVar> & preconditioner,
                                                                                                 const TPZFMatrix<TVar>& F, TPZFMatrix<TVar>& result,
                                                                                                 TPZFMatrix<TVar>* residual, REAL & tol,
                                                                                                 const int FromCurrent = 0);
                
                virtual bool SolveEigenvaluesJacobi(int64_t &numiterations, REAL &tol, TPZVec<TVar> * Sort = 0);
                
                virtual bool SolveEigensystemJacobi(int64_t &numiterations, REAL & tol, TPZVec<TVar> & Eigenvalues, TPZFMatrix<TVar>& Eigenvectors) const;
                
                virtual int SolveDirect ( TPZFMatrix<TVar>& F , const DecomposeType dt, std::list<int64_t> &singular);
                virtual int SolveDirect ( TPZFMatrix<TVar>& F , const DecomposeType dt);
                
    virtual int Decompose(const DecomposeType dt, std::list<int64_t> &singular)
    {
        switch (dt) {
            case ELU:
                return Decompose_LU(singular);
                break;
            case ELDLt:
                return Decompose_LDLt(singular);
                break;
            case ECholesky:
                return Decompose_Cholesky(singular);
                break;                
            default:
                DebugStop();
                break;
        }
        return -1;
    }
                static TVar ReturnNearestValue(TVar val, TPZVec<TVar> &Vec, TVar tol);
                
                int Solve_LU ( TPZFMatrix<TVar>* B, std::list<int64_t> &singular );
                int Solve_LU ( TPZFMatrix<TVar>* B );
                
                virtual int Solve_Cholesky( TPZFMatrix<TVar>* B);
                int Solve_Cholesky( TPZFMatrix<TVar>* B, std::list<int64_t> &singular );
                int Solve_LDLt    ( TPZFMatrix<TVar>* B, std::list<int64_t> &singular );
                int Solve_LDLt    ( TPZFMatrix<TVar>* B);
                
                virtual int Decompose_LU(std::list<int64_t> &singular);
                virtual int Decompose_LU();
                
                virtual int Decompose_Cholesky() ;
                virtual int Decompose_Cholesky(std::list<int64_t> &singular) ;
                
                virtual int Decompose_LDLt(std::list<int64_t> &singular);
                virtual int Decompose_LDLt();
                
                virtual int Substitution( TPZFMatrix<TVar>* B ) const;
                
                virtual int Subst_Forward( TPZFMatrix<TVar>* b ) const;
                
                virtual int Subst_Backward( TPZFMatrix<TVar>* b ) const;
                
                virtual int Subst_LForward( TPZFMatrix<TVar>* b ) const;
                
                virtual int Subst_LBackward( TPZFMatrix<TVar>* b ) const;
                
                virtual int Subst_Diag( TPZFMatrix<TVar>* b ) const;
                
        public:
int ClassId() const override;

        
                void Read(TPZStream &buf, void *context) override;
                
                void Write(TPZStream &buf, int withclassid) const override;
                

                virtual bool Compare(TPZSavable *copy, bool override = false) override;

                virtual bool Compare(TPZSavable *copy, bool override = false) const override;
                
                virtual void GetSub(const TPZVec<int64_t> &indices,TPZFMatrix<TVar>&block) const;
                
    bool CompareValues(TPZMatrix<TVar>&M, TVar tol);
                
protected:
                
                void PrepareZ(const TPZFMatrix<TVar>& y, TPZFMatrix<TVar>& z,const TVar beta,const int opt) const;
                
                inline  TPZMatrix<TVar>(const int64_t row,const int64_t col ) : TPZRegisterClassId(&TPZMatrix<TVar>::ClassId)
                { fRow = row; fCol = col;fDefPositive=0; fDecomposed = 0;}
                
public:
                static int Error(const char *msg ,const char *msg2 = 0);
                
protected:
                virtual int Clear() { return 0; }
                
                static void Swap(int64_t *a, int64_t *b);
                int64_t fRow;
                int64_t fCol;
                char  fDecomposed;
                char fDefPositive;
                static TVar gZero;
};

template <class TVar>
TVar TPZMatrix<TVar>::gZero = TVar(0);

template<class TVar>
std::ostream & operator<<(std::ostream& out, const TPZMatrix<TVar> & A);

/******** Inline ********/


template<class TVar>
inline int64_t TPZMatrix<TVar>::Rows() const {
                return fRow;
}


template<class TVar>
inline int64_t TPZMatrix<TVar>::Cols() const {
                return fCol;
}

#ifdef _AUTODIFF
template<class TVar>
inline void TPZMatrix<TVar>::Residual(const TPZFMatrix<TVar>& x,const TPZFMatrix<TVar>& rhs, TPZFMatrix<TVar>& res )  {
                DebugStop();
}


template<>
inline void TPZMatrix<double>::Residual(const TPZFMatrix<double>& x,const TPZFMatrix<double>& rhs, TPZFMatrix<double>& res )  {
                MultAdd( x, rhs, res, ((double)-1.0), ((double)1.0) );
}

#else
template<class TVar>
inline void TPZMatrix<TVar>::Residual(const TPZFMatrix<TVar>& x,const TPZFMatrix<TVar>& rhs, TPZFMatrix<TVar>& res )  {
                MultAdd( x, rhs, res, ((TVar)-1.0), ((TVar)1.0) );
}
#endif

/***********/
/*** Put ***/

template<class TVar>
inline int TPZMatrix<TVar>::Put(const int64_t row,const int64_t col,const TVar & value ) {
                // verificando se o elemento a inserir esta dentro da matriz
#ifdef PZDEBUG
                if ( row >= Rows() || col >= Cols() || row <0 || col < 0 ) {
                                std::stringstream sout;
                                sout << "TPZMatrix<TVar>::Put" << " Index out of range row = " << row << " col = " << col << " Rows() " << Rows() << " Cols() " << Cols() << std::endl;
                                Error(sout.str().c_str());
                }
#endif
                
                return( PutVal( row, col, value ) );
}



/***********/
/*** Get ***/

template<class TVar>
inline const TVar &TPZMatrix<TVar>::Get(const int64_t row, const int64_t col ) const {
                // verificando se o elemento pedido esta dentro da matriz
#ifdef PZDEBUG
                if ( (row >= Rows()) || (col >= Cols()) || row <0 || col <0 ) {
                                Error("TPZMatrix::Get", "Index out of range");
        DebugStop();
                }
#endif
                return( GetVal( row, col ) );
}

template<class TVar>
inline TVar &TPZMatrix<TVar>::operator()(const int64_t row, const int64_t col) {
                // verificando se o elemento a inserir esta dentro da matriz
#ifndef NODEBUG
                if ( (row >= Rows()) || (col >= Cols()) || row <0 || col<0 ) {
                                Error("TPZMatrix<TVar>::Operator()","Index out of range");
        DebugStop();
                }
#endif
                return s(row,col);
}

template<class TVar>
inline TVar &TPZMatrix<TVar>::s(const int64_t row, const int64_t col) {
                // verificando se o elemento a inserir esta dentro da matriz
    DebugStop();
    throw "TPZMatrix<TVar>::s not implemented\n";
    
}

template<class TVar>
inline TVar &TPZMatrix<TVar>::operator()(const int64_t row) {
                return operator()(row,0);
}

template<class TVar>
inline int64_t TPZMatrix<TVar>::Dim() const{
                if ( IsSquare() ) return Rows();
                Error( "matrix is not square" );
                return ( 0 );
}
//***Solve LU ***/

template<class TVar>
inline int TPZMatrix<TVar>::Solve_LU( TPZFMatrix<TVar>*B, std::list<int64_t> &singular) {
                if ( IsSimetric() )
        Error( "LU decomposition is a not symetric decomposition" );
                return ( ( !Decompose_LU(singular) )?  0 : Substitution( B )  );
}

template<class TVar>
inline int TPZMatrix<TVar>::Solve_LU( TPZFMatrix<TVar>*B ) {
                if ( IsSimetric() )
        Error( "LU decomposition is a not symetric decomposition" );
                return ( ( !Decompose_LU() )?  0 : Substitution( B )  );
}
/**********************/
/*** Solve Cholesky ***/
//
//  Se nao conseguir resolver por Cholesky retorna 0 e a matriz
//   sera' modificada (seu valor perdera' o sentido).
//
template<class TVar>
inline int TPZMatrix<TVar>::Solve_Cholesky( TPZFMatrix<TVar>* B )
{
                return(
                                   ( !Decompose_Cholesky() )?  0 :( Subst_Forward( B ) && Subst_Backward( B ) )
                                   );
}

template<class TVar>
inline int TPZMatrix<TVar>::Solve_Cholesky( TPZFMatrix<TVar>* B, std::list<int64_t> &singular ) {
                return(
                                   ( !Decompose_Cholesky(singular) )?  0 :( Subst_Forward( B ) && Subst_Backward( B ) )
                                   );
}

/******************/
/*** Solve LDLt ***/

template<class TVar>
inline int TPZMatrix<TVar>::Solve_LDLt( TPZFMatrix<TVar>* B ) {
                
                return(
                                   ( !Decompose_LDLt() )? 0 :
                                   ( Subst_LForward( B ) && Subst_Diag( B ) && Subst_LBackward( B ) )
                                   );
}



template<class TVar>
inline void
TPZMatrix<TVar>::Swap( int64_t *a, int64_t *b )
{
                int64_t aux = *a;
                *a = *b;
                *b = aux;
}

template<class TVar>
int TPZMatrix<TVar>::ClassId() const{
    return Hash("TPZMatrix") ^ ClassIdOrHash<TVar>()<<1;
}

#include "pzfmatrix.h"
#include "pzsolve.h"

#endif