pzysmp.h

Go to the documentation of this file.

#ifndef YSMPMATH
#define YSMPMATH

#include "pz_config.h"

#ifdef USING_BLAS
#ifdef USING_MKL
#include <mkl.h>
#elif MACOSX
#include <Accelerate/Accelerate.h>
#else
#ifdef MACOSX
#include <Accelerate/Accelerate.h>
#else
extern "C"{
     #include "cblas.h"
     };
#endif
#endif
#endif

template<class TVar>
class TPZVerySparseMatrix;

#include "pzmatrix.h"
#include "pzfmatrix.h"

#ifdef USING_MKL
#include "TPZPardisoControl.h"
#endif

template<class TVar>
class TPZFYsmpMatrix : public TPZMatrix<TVar> {
                
#ifdef USING_MKL
    friend class TPZPardisoControl<TVar>;
#endif
    
                public :
                
                struct TPZMThread {
                                const TPZFYsmpMatrix<TVar> *target;
                                int64_t fFirsteq;
                                int64_t fLasteq;
                                const TPZFMatrix<TVar> *fX;
                                TPZFMatrix<TVar> *fZ;
                                TVar fAlpha;
                                int fOpt;
                };
                
private:
                
                static void * ExecuteMT(void *entrydata);
                
public:
    
    TPZFYsmpMatrix();
    
    TPZFYsmpMatrix(const int64_t rows,const int64_t cols );
                
                TPZFYsmpMatrix(const TPZVerySparseMatrix<TVar> &cp);
                
                TPZFYsmpMatrix &operator=(const TPZFYsmpMatrix<TVar> &copy);
                
                TPZFYsmpMatrix &operator=(const TPZVerySparseMatrix<TVar> &cp);
    
                CLONEDEF(TPZFYsmpMatrix)
                
                virtual ~TPZFYsmpMatrix();      
                
    void AutoFill(int64_t nrow, int64_t ncol, int symmetric);
    

    
                virtual const TVar &GetVal(const int64_t row,const int64_t col ) const override;
                
                int64_t NumTerms()
                {
                                return fIA[this->Rows()];
                }
                
                int PutVal(const int64_t row, const int64_t col, const TVar &Value) override;
                
                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;
                
                virtual void MultAddMT(const TPZFMatrix<TVar> &x,const TPZFMatrix<TVar> &y, TPZFMatrix<TVar> &z,
                                                                                                   const TVar alpha=1.,const TVar beta = 0., const int opt = 0);
                
                virtual int GetSub(const int64_t sRow,const int64_t sCol,const int64_t rowSize,
                                                                                   const int64_t colSize, TPZFMatrix<TVar> & A ) const override;
                
                void GetSub(const TPZVec<int64_t> &indices,TPZFMatrix<TVar> &block) const override;
                
                virtual void SetData( int64_t *IA, int64_t *JA, TVar *A );
    
    virtual void SetData( TPZVec<int64_t> &IA, TPZVec<int64_t> &JA, TPZVec<TVar> &A );
                
                virtual void Print(const char *title, std::ostream &out = std::cout , const MatrixOutputFormat form = EFormatted) const override;
                
                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)  override;
                
                void SolveSOR(int64_t &numiterations, const TPZFMatrix<TVar> &rhs, TPZFMatrix<TVar> &x,
                                                                  TPZFMatrix<TVar> *residual, TPZFMatrix<TVar> &scratch,
                                                                  const REAL overrelax, REAL &tol,
                                                                  const int FromCurrent = 0,const int direction = 1 )  override;
                // @}
                
                virtual void AddKelOld(
                                                                                                   TPZFMatrix<TVar> & elmat 
                                                                                                   , TPZVec < int > & destinationindex 
                                                                                                   );    
                
                virtual void AddKel(TPZFMatrix<TVar> & elmat, TPZVec<int64_t> & destinationindex) override;
                
                virtual void AddKel(TPZFMatrix<TVar> & elmat, TPZVec<int64_t> & sourceindex, TPZVec<int64_t> & destinationindex) override;
                
                void MultiplyDummy(TPZFYsmpMatrix<TVar> & B, TPZFYsmpMatrix<TVar> & Res);
                
                virtual int Zero() override;
                
                virtual int Decompose_LU(std::list<int64_t> &singular) override;
                virtual int Decompose_LU() override;
                
                
                virtual int Substitution( TPZFMatrix<TVar> * B ) const override;
                
                
    public:
int ClassId() const override;

private:
                
                void ComputeDiagonal();
                
                /*
                 * @brief Perform row update of the sparse matrix
                 */
                void RowLUUpdate(int64_t sourcerow, int64_t destrow);
                
protected:
                TPZVec<int64_t>  fIA;
                TPZVec<int64_t>  fJA;
                TPZVec<TVar> fA;
                
                TPZVec<TVar> fDiag;
                
                int   fSymmetric;
                
#ifdef USING_MKL    
    TPZPardisoControl<TVar> fPardisoControl;
#endif
protected:
                
                void InitializeData();
};


template<class TVar>
inline void TPZFYsmpMatrix<TVar>::SetData( int64_t *IA, int64_t *JA, TVar *A ) {
                // Pass the data to the class.
    int nel = this->Rows()+1;
    fIA.resize(nel);
    memccpy(&fIA[0], IA, nel, sizeof(int64_t));
    int64_t nval = fIA[nel-1];
    fJA.resize(nval);
    memccpy(&fJA[0], JA, nval, sizeof(int64_t));
    fA.resize(nval);
    memccpy(&fA[0], A, nval, sizeof(TVar));
                ComputeDiagonal();
}

template<class TVar>
inline void TPZFYsmpMatrix<TVar>::SetData( TPZVec<int64_t> &IA, TPZVec<int64_t> &JA, TPZVec<TVar> &A ){
    
    if (IA.size() != this->Rows() + 1 ) {
        DebugStop();
    }
    
    if (JA.size() != IA[this->Rows()]) {
        DebugStop();
    }
    
    fIA = IA;
    fJA = JA;
    fA = A;
}

#endif