pzfmatrix.h

Go to the documentation of this file.

#ifndef _TFULLMATRIXH_
#define _TFULLMATRIXH_

#include <math.h>            // for sqrt
#include <string.h>          // for NULL, memset
#include <complex>           // for complex
#include <iostream>          // for operator<<, ostream, cout
#include <list>              // for list
#include <sstream>           // for basic_stringbuf<>::int_type, basic_strin...
#include "pzmatrix.h"                            // for TPZMatrix
#include "pzerror.h"         // for DebugStop
#include "pzmanvector.h"     // for TPZManVector
#include "pzreal.h"          // for TPZFlopCounter, IsZero, REAL, sqrt, fabs
#include "pzvec.h"           // for TPZVec
#include "tpzautopointer.h"  // for TPZAutoPointer
#include <cstdlib>

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

class TPZSavable;
class TPZStream;
template <class TVar> class TPZFMatrix;
template <class TVar> class TPZVerySparseMatrix;


template <class T>
class TPZVec;

#define GETVAL(MAT,rows,row,col) MAT->fElem[((unsigned)col)*rows+row]

#define PUTVAL(MAT,rows,row,col,val) MAT->fElem[((unsigned)col)*rows+row]=val

#define SELECTEL(ptr,rows,row,col) ptr[col*rows+row]


template<class TVar>
TVar Dot(const TPZFMatrix<TVar> &A,const TPZFMatrix<TVar> &B);

template<class TVar>
TVar Norm(const TPZFMatrix<TVar> &A);


template<class TVar=REAL>
class TPZFMatrix: public TPZMatrix<TVar> {
    
public:
    
    TPZFMatrix() : TPZRegisterClassId(&TPZFMatrix<TVar>::ClassId),
    TPZMatrix<TVar>( 0, 0 ), fElem(0),fGiven(0),fSize(0) {}
    TPZFMatrix(const int64_t rows ,const int64_t columns, TVar* buf,const int64_t size);
    TPZFMatrix(const int64_t rows ,const int64_t columns,const TVar & val );
    inline  TPZFMatrix(const int64_t rows ,const int64_t columns = 1) : 
    TPZRegisterClassId(&TPZFMatrix<TVar>::ClassId),
    TPZMatrix<TVar>(rows,columns), fElem(0),fGiven(0),fSize(0) {
        if(rows*columns) fElem = new TVar[rows*columns];
    }
    
    TPZFMatrix(TPZVerySparseMatrix<TVar> const & A);
    
    TPZFMatrix(const TPZFMatrix<TVar> & refmat);
    
    
    
    CLONEDEF(TPZFMatrix<TVar>)
    TPZFMatrix(const TPZMatrix<TVar> & refmat);

                inline TPZFMatrix(const std::initializer_list<TVar>& list);

    inline TPZFMatrix(const std::initializer_list< std::initializer_list<TVar> > &list);

    
    virtual  ~TPZFMatrix();
    
    int64_t MemoryFootprint() const  override
    {
        return (sizeof(TVar)*this->Rows()*this->Cols());
    }
    
    TVar *Adress()
    {
        return fElem;
    }
    
    friend class TPZFMatrix<float>;
    friend class TPZFMatrix<double>;

    template<class TVar2>
    void CopyFrom(TPZFMatrix<TVar2> &orig)
    {
        Resize(orig.Rows(), orig.Cols());
        TPZMatrix<TVar>::CopyFrom(orig);
        int64_t nel = orig.Rows()*orig.Cols();
        for (int64_t el=0; el<nel; el++) {
            fElem[el] = orig.fElem[el];
        }
    }

    virtual void UpdateFrom(TPZAutoPointer<TPZMatrix<TVar> >  mat) override
    {
        TPZMatrix<TVar> *matptr = mat.operator->();
        TPZFMatrix<TVar> *from = dynamic_cast<TPZFMatrix<TVar> *>(matptr);
        if (from) {
            *this = *from;
        }
        else
        {
            std::cout << "TPZMatrix<TVar>::UdateFrom is not implemented\n";
            DebugStop();
        }
    }

    int PutVal(const int64_t row,const int64_t col,const TVar & value ) override;
    const TVar &GetVal(const int64_t row,const int64_t col ) const override;
    
    virtual TVar &s(const int64_t row, const int64_t col) override;
    
    TVar &g(const int64_t row, const int64_t col) const;
    void AddFel(TPZFMatrix<TVar> &rhs,TPZVec<int64_t> &destination);
    void AddFel(TPZFMatrix<TVar> &rhs,TPZVec<int64_t> &source, TPZVec<int64_t> &destination);
    
    
    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  override;
    
    
    static void MultAdd(const TVar *ptr, int64_t rows, int64_t cols, const TPZFMatrix<TVar> &x,const TPZFMatrix<TVar> &y, TPZFMatrix<TVar> &z,
                        const TVar alpha=1.,const TVar beta = 0.,const int opt = 0);
    
    TVar &operator()(const int64_t row,const int64_t col);
    TVar &operator()(const int64_t row);
    virtual TPZFMatrix&operator= (const TPZFMatrix<TVar> &A );
                TPZFMatrix<TVar>& operator= (const std::initializer_list<TVar>& list);
                TPZFMatrix<TVar>& operator= (const std::initializer_list< std::initializer_list<TVar> >& list);
    TPZFMatrix<TVar> operator+  (const TPZFMatrix<TVar> &A ) const;
    TPZFMatrix<TVar> operator-  (const TPZFMatrix<TVar> &A ) const;
    TPZFMatrix<TVar> operator*  ( TPZFMatrix<TVar> A ) const ;
    TPZFMatrix<TVar> &operator+=(const TPZFMatrix<TVar> &A );
    TPZFMatrix<TVar> &operator-=(const TPZFMatrix<TVar> &A );
    
    void ZAXPY(const TVar alpha, const TPZFMatrix<TVar> &p);
    void TimesBetaPlusZ(const TVar beta, const TPZFMatrix<TVar> &z);
    
    TPZFMatrix<TVar> &operator= (const TPZMatrix<TVar> &A );
    
    TPZFMatrix<TVar> &operator= (const TVar val );
    TPZFMatrix<TVar> operator+  (const TVar val ) const;
    TPZFMatrix<TVar> operator-  (const TVar val ) const;
    TPZFMatrix<TVar> operator*  (const TVar val ) const;
    TPZFMatrix<TVar> &operator+=(const TVar val );
    TPZFMatrix<TVar> &operator-=(const TVar val )  { return operator+=( -val ); }
    TPZFMatrix<TVar> &operator*=(const TVar val );
    
    //          TPZFMatrix<TVar> operator-() const;// { return operator*( -1.0 ); }
    
    int Resize(const int64_t newRows,const int64_t wCols ) override;
    
    int SetSize(int64_t newRows, int64_t newCols);
    
    int Remodel(const int64_t newRows,const int64_t wCols );
    
    int Redim(const int64_t newRows,const int64_t newCols ) override;
    
    int Zero() override;
    
#ifdef USING_LAPACK

    void InitializePivot();
#endif

    
    void GramSchmidt(TPZFMatrix<TVar> &Orthog, TPZFMatrix<TVar> &TransfToOrthog);
    
    void DeterminantInverse(TVar &determinant, TPZFMatrix<TVar> &inverse);
    
    void Transpose(TPZMatrix<TVar> *const T) const override;
    
    void Transpose();
    
    /*** @name Solve linear system of equations ***/
    virtual int Decompose_Cholesky() override;
    virtual int Decompose_Cholesky(std::list<int64_t> &singular) override;
    
    virtual int Decompose_LU(std::list<int64_t> &singular) override;
    virtual int Decompose_LU() override;
    
    virtual int Decompose_LDLt() override;
    
    static int Substitution(const TVar *ptr, int64_t rows, TPZFMatrix<TVar> *B);
    
    virtual int Substitution( TPZFMatrix<TVar> *B ) const override;
    
    virtual int Decompose_LU(TPZVec<int> &index);
    
    virtual int Substitution( TPZFMatrix<TVar> *B, const TPZVec<int> &index ) const;
    
    static int Substitution(const TVar *ptr, int64_t rows,  TPZFMatrix<TVar> *B, const TPZVec<int> &index );
    
#ifdef USING_LAPACK

    virtual int Subst_Forward( TPZFMatrix<TVar>* b ) const override;
    
    virtual int Subst_Backward( TPZFMatrix<TVar>* b ) const override;
    
    virtual int Subst_LForward( TPZFMatrix<TVar>* b ) const override;
    
    virtual int Subst_LBackward( TPZFMatrix<TVar>* b ) const override;
    
    virtual int Subst_Diag( TPZFMatrix<TVar>* b ) const override;
#endif
    
#ifdef USING_LAPACK
    /*** @name Solve eigenvalues ***/
    virtual int SolveEigenProblem(TPZVec < std::complex<double> > &w, TPZFMatrix < std::complex<double> > &eigenVectors);
    virtual int SolveEigenProblem(TPZVec < std::complex<double> > &w);

    virtual int SolveGeneralisedEigenProblem(TPZFMatrix< TVar > &B , TPZVec < std::complex<double> > &w, TPZFMatrix < std::complex<double> > &eigenVectors);
    virtual int SolveGeneralisedEigenProblem(TPZFMatrix< TVar > &B , TPZVec < std::complex<double> > &w);    
#endif
    
    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;
    
    operator const TVar*() const { return fElem; }
    
    static void PrintStatic(const TVar *ptr, int64_t rows, int64_t cols, const char *name, std::ostream& out,const MatrixOutputFormat form);
    
private:
    
    static int Error(const char *msg1,const char *msg2=0 );
    int Clear() override;
    
    TVar *fElem;
    TVar *fGiven;
    int64_t fSize;
    
#ifdef USING_LAPACK
    TPZManVector<int,5> fPivot;
    
    TPZVec<TVar> fWork;
#endif
};

template<class TVar>
inline TPZFMatrix<TVar>::TPZFMatrix(const int64_t rows,const int64_t cols,TVar * buf,const int64_t sz)
: TPZRegisterClassId(&TPZFMatrix<TVar>::ClassId),TPZMatrix<TVar>( rows, cols ), 
fElem(buf),fGiven(buf),fSize(sz) {
    int64_t size = rows * cols;
    if(size == 0)
    {
        fElem = NULL;
    }
    else if(size > sz)
    {
        fElem=new TVar[size];
#ifndef NODEBUG
        if ( fElem == NULL && size) Error( "Constructor <memory allocation error>." );
#endif
    }
}
template<class TVar>
inline TPZFMatrix<TVar>::TPZFMatrix(const int64_t rows,const int64_t cols,const TVar & val )
: TPZRegisterClassId(&TPZFMatrix<TVar>::ClassId), TPZMatrix<TVar>( rows, cols ), 
 fElem(0), fGiven(0), fSize(0) {
    int64_t size = rows * cols;
    if(!size) return;
    fElem=new TVar[size];
#ifdef PZDEBUG
    if ( fElem == NULL && size) Error( "Constructor <memory allocation error>." );
#endif
    for(int64_t i=0;i<size;i++) fElem[i] = val;
}

template<class TVar>
inline TPZFMatrix<TVar>::TPZFMatrix(const std::initializer_list<TVar>& list) 
: TPZRegisterClassId(&TPZFMatrix<TVar>::ClassId), TPZMatrix<TVar>(list.size(), 1)
{
                if (list.size() > 0) {
                                this->fElem = new TVar[list.size()];
                }

                auto it = list.begin();
                auto it_end = list.end();
                TVar* aux = fElem;
                for (; it != it_end; it++, aux++)
                                *aux = *it;
}

template< class TVar >
inline TPZFMatrix<TVar>::TPZFMatrix(const std::initializer_list<std::initializer_list<TVar>> &list)
: TPZRegisterClassId(&TPZFMatrix<TVar>::ClassId)
{
    this->fRow = list.size();

    auto row_it = list.begin();
    auto row_it_end = list.end();

#ifdef PZDEBUG    
    bool col_n_found = false;
    for (auto it = row_it; it != row_it_end; it++) {
        if (!col_n_found) {
            this->fCol = it->size();
            col_n_found = true;
        } else {
            if (this->fCol != it->size())
                Error("TPZFMatrix constructor: inconsistent number of columns in initializer list");
        }
    }
#else
    this->fCol = row_it->size();
#endif

                this->fElem = new TVar[this->fRow * this->fCol];
    
    for (uint32_t row_n = 0; row_it != row_it_end; row_it++, row_n++) {
        auto col_it = row_it->begin();
        auto col_it_end = row_it->end();
        for (uint32_t col_n = 0; col_it != col_it_end; col_it++, col_n++) {
            this->fElem[col_n * this->fRow + row_n] = *col_it;
        }
    }
}

template<class TVar>
inline TPZFMatrix<TVar> operator*(TVar val, const TPZFMatrix<TVar> &A)
{
    return A*val;
}


/*******************************/
/*** Operator*( TPZMatrix<TVar> & ) ***/
template<class TVar>
inline TPZFMatrix<TVar> TPZFMatrix<TVar>::operator*( TPZFMatrix<TVar> A ) const {
    if ( this->Cols() != A.Rows() )
        Error( "Operator* <matrixs with incompatible dimensions>" );
    
    TPZFMatrix<TVar> res;
    res.Redim( this->Rows(), A.Cols() );
    MultAdd(A,A,res,1.,0.,0);
    return( res );
}

/**************/
/*** PutVal ***/
template<class TVar>
inline int TPZFMatrix<TVar>::PutVal(const int64_t row, const int64_t col,const TVar & value ) {
    fElem[ ((unsigned)col) * this->Rows() + row ] = value;
    return( 1 );
}

/******************/
/*** Destructor ***/
template<class TVar>
inline TPZFMatrix<TVar>::~TPZFMatrix() {
    if(fElem && fElem != fGiven) delete[]( fElem );
    fElem = 0;
    fSize = 0;
}


/**************/
/*** GetVal ***/
template<class TVar>
inline const TVar &TPZFMatrix<TVar>::GetVal( const int64_t row, const int64_t col ) const {
#ifdef PZDEBUG
    if(row >=  this->Rows() || row<0 || col >=  this->Cols() || col<0) {
        Error("TPZFMatrix::operator() "," Index out of bounds");
        DebugStop();
    }
#endif
    return( fElem[ col*this->fRow + row ] );
}

template<class TVar>
inline TVar &TPZFMatrix<TVar>::operator()( const int64_t row, const int64_t col) {
#ifndef NODEBUG
    if(row >=  this->Rows() || row<0 || col >=  this->Cols() || col<0) {
        Error("TPZFMatrix<TVar>::operator() "," Index out of bounds");
        DebugStop();
    }
#endif
    return *(this->fElem+col*this->fRow+row);
}

template<class TVar>
inline TVar &TPZFMatrix<TVar>::s(const int64_t row, const int64_t col) {
    // verificando se o elemento a inserir esta dentro da matriz
    return operator()(row,col);
}

template<class TVar>
inline TVar &TPZFMatrix<TVar>::g( const int64_t row, const int64_t col) const {
#ifdef PZDEBUG
    if(row >=  this->Rows() || row<0 || col >=  this->Cols() || col<0) {
        Error("TPZFMatrix<TVar>::operator() "," Index out of bounds");
        DebugStop();
    }
#endif
    return *(this->fElem+col*this->fRow+row);
}

template<class TVar>
inline TVar &TPZFMatrix<TVar>::operator()(const int64_t row) {
#ifdef PZDEBUG
    if(row >=  this->Rows() || row<0) {
        Error("TPZFMatrix<TVar>::operator() "," Index out of bounds");
        DebugStop();
    }
#endif
    return *(this->fElem+row);
}

template<class TVar>
inline int TPZFMatrix<TVar>::Redim(const int64_t newRows,const int64_t newCols) {
    int64_t newsize = newRows*newCols;
    int64_t size = this->fRow*this->fCol;
    if ( newsize == size) {
        this->fRow = newRows;
        this->fCol = newCols;
        Zero();
        return( 1 );
    }
    if(this->fElem && this->fElem != this->fGiven) delete []this->fElem;
    
    if(this->fGiven && newsize <= this->fSize) {
        this->fElem = this->fGiven;
    } else if(newsize == 0) {
        this->fElem = NULL;
    } else {
        this->fElem = new TVar[ newsize ] ;
    }
#ifndef NODEBUG
    if (newsize && this->fElem == NULL )
        Error( "Resize <memory allocation error>." );
#endif
    
    this->fRow  = newRows;
    this->fCol  = newCols;
    
    Zero();
    
    return( 1 );
}

/***************/
/****Zero*******/

template<class TVar>
inline int TPZFMatrix<TVar>::Zero() {
    int64_t size = this->fRow * this->fCol * sizeof(TVar);
    memset(((void*)this->fElem),'\0',size);
    this->fDecomposed = 0;
    return( 1 );
}

/**************************/
/*** Operations Global ****/

inline int64_t Norm(const TPZFMatrix<int64_t> &A) {
    return (int64_t)sqrt((REAL)Dot(A,A));
}

inline int Norm(const TPZFMatrix<int> &A) {
    return (int)sqrt((REAL)Dot(A,A));
}

inline float Norm(const TPZFMatrix<float> &A) {
    return sqrt(Dot(A,A));
}

inline double Norm(const TPZFMatrix<double> &A) {
    return sqrt(Dot(A,A));
}

inline long double Norm(const TPZFMatrix<long double> &A) {
    return sqrt(Dot(A,A));
}

inline float Norm(const TPZFMatrix< std::complex <float> > &A) {
    return sqrt(Dot(A,A).real());
}

inline double Norm(const TPZFMatrix< std::complex <double> > &A) {
    return sqrt(Dot(A,A).real());
}

inline long double Norm(const TPZFMatrix< std::complex <long double> > &A) {
    return sqrt(Dot(A,A).real());
}

#ifdef _AUTODIFF
inline float Norm(const TPZFMatrix< Fad <float> > &A) {
    return (float)TPZExtractVal::val(sqrt(Dot(A,A)));
}

inline double Norm(const TPZFMatrix< Fad <double> > &A) {
    return TPZExtractVal::val(sqrt(Dot(A,A)));
}

inline long double Norm(const TPZFMatrix< Fad <long double> > &A) {
    return TPZExtractVal::val(sqrt(Dot(A,A)));
}
#endif

inline TPZFlopCounter Norm(const TPZFMatrix<TPZFlopCounter> &A)
{
    return sqrt(Dot(A, A));
}
template<int N, class TVar=REAL>
class TPZFNMatrix : public TPZFMatrix<TVar> {
    
    TVar fBuf[N+1];
    
public:
    /*
     * @brief Constructor which does not initialize the data. \n
     * WARNING : this class will dynamically allocate memory if the template parameter N is smaller than row*col
     * @param row Number of rows
     * @param col Number of cols
     */
    inline TPZFNMatrix(int64_t row, int64_t col) : TPZRegisterClassId(&TPZFNMatrix<N,TVar>::ClassId), 
    TPZFMatrix<TVar>(row,col,fBuf,N) {}
    
    inline TPZFNMatrix() : TPZRegisterClassId(&TPZFNMatrix<N,TVar>::ClassId), TPZFMatrix<TVar>(0,0,fBuf,N)
    {
    }
    
    inline TPZFNMatrix(const TPZFMatrix<TVar> &copy) : TPZRegisterClassId(&TPZFNMatrix<N,TVar>::ClassId), 
    TPZFMatrix<TVar>(0,0,fBuf,N)
    {
        *this = copy;
    }
    
    inline TPZFNMatrix(const TPZFNMatrix<N,TVar> &copy) : TPZRegisterClassId(&TPZFNMatrix<N,TVar>::ClassId), 
    TPZFMatrix<TVar>(0,0,fBuf,N)
    {
        *this = copy;
    }
    
    virtual ~TPZFNMatrix()
    {
    }
    
    CLONEDEF(TPZFNMatrix)
    /*
     * @brief Constructor which initializes the data. \n
     * WARNING : this class will dynamically allocate memory if the template parameter N is smaller than row*col
     * @param row Number of rows
     * @param col Number of cols
     * @param val initial value of the matrix elements
     */
    inline  TPZFNMatrix(int64_t row, int64_t col, const TVar &val) : TPZRegisterClassId(&TPZFNMatrix<N,TVar>::ClassId), 
    TPZFMatrix<TVar>(row,col,fBuf,N) {
        TPZFMatrix<TVar>::operator=(val);
    }
    
    inline  TPZFMatrix<TVar> &operator=(const TPZFMatrix<TVar> &copy)  override {
        return TPZFMatrix<TVar>::operator=(copy);
    }
    inline  TPZFNMatrix<N, TVar> &operator=(const TPZFNMatrix<N, TVar> &copy) {
        TPZFMatrix<TVar>::operator=(copy);
        return *this;
    }
    
    int ClassId() const override;
    void Read(TPZStream &buf, void *context) override;
    void Write(TPZStream &buf, int withclassid) const override;
    
};

template<int N, class TVar>
int TPZFNMatrix<N, TVar>::ClassId() const{
    return Hash("TPZFNMatrix") ^ TPZFMatrix<TVar>::ClassId() << 1 ^ (N << 2);
}

template<int N, class TVar>
void TPZFNMatrix<N, TVar>::Read(TPZStream& buf, void* context) {
    TPZFMatrix<TVar>::Read(buf, context);
    buf.Read(fBuf, N+1);
}

template<int N, class TVar>
void TPZFNMatrix<N, TVar>::Write(TPZStream& buf, int withclassid) const {
    TPZFMatrix<TVar>::Write(buf, withclassid);
    buf.Write(fBuf, N+1);
}

#endif