CDMATH/linearsolver/inc/SparseMatrixPetsc.hxx

   1 /*
   2  * SparseMatrixPetsc.hxx
   3  *
   4  *  Created on: 03/11/2017. 2017
   5  *      Author: mndjinga
   6  */
   7
   8 #ifndef SOURCE_DIRECTORY__BASE_INC_SparseMatrixPetsc_HXX_
   9 #define SOURCE_DIRECTORY__BASE_INC_SparseMatrixPetsc_HXX_
  10
  11 #include <iostream>
  12 #include <vector>
  13
  14 #include "MEDCouplingUMesh.hxx"
  15 #include "Vector.hxx"
  16 #include <petsc.h>
  17
  18 #include <slepceps.h>
  19 #include <slepcsvd.h>
  20
  21 class SparseMatrixPetsc: public GenericMatrix {
  22
  23 public:
  24         SparseMatrixPetsc();
  25         virtual ~SparseMatrixPetsc();
  26
  27         /**
  28          * constructor for a block sparse matrix
  29          * @param numberOfRows : The number of rows
  30          * @param numberOfColumns : The number of columns
  31          */
  32         SparseMatrixPetsc( int numberOfRows, int numberOfColumns ) ;
  33
  34         /**
  35          * constructor for a block sparse matrix with number of non zero coefficients given
  36          * @param numberOfRows : The number of rows
  37          * @param numberOfColumns : The number of columns
  38          * @param nnz : The maximum number of nonzeros coefficients per line (or an upper bound)
  39          */
  40         SparseMatrixPetsc( int numberOfRows, int numberOfColumns, int nnz ) ;
  41
  42         /**
  43          * constructor for a sparse matrix with block structure
  44          * @param blockSize : The block size
  45          * @param numberOfRows : The number of rows
  46          * @param numberOfColumns : The number of columns
  47          * @param nnz : The maximum number of nonzeros coefficients per line (or an upper bound)
  48      * @comment blockSize should always divide numberOfRows and numberOfColumns
  49          */
  50     SparseMatrixPetsc( int blockSize, int numberOfRows, int numberOfColumns, int nnz );
  51
  52         /**
  53          * constructor by copy
  54          * @param SparseMatrixPetsc : The SparseMatrixPetsc object to be copied
  55          */
  56         SparseMatrixPetsc ( const SparseMatrixPetsc& sparseMatrixPetsc ) ;
  57
  58         /**
  59          * constructor with data
  60          * @param Petsc matris : mat
  61          */
  62         SparseMatrixPetsc(Mat mat);
  63
  64         SparseMatrixPetsc transpose() const ;
  65
  66         double operator()( int i, int j ) const ;
  67
  68         void setValue( int i, int j, double value ) ;
  69         void addValue( int i, int j, double value ) ;
  70
  71         void setValue( int i, int j, Matrix M ) ;
  72         void addValue( int i, int j, Matrix M ) ;
  73
  74         void setValuesBlocked( int i, int j, Matrix M ) ;
  75         void addValuesBlocked( int i, int j, Matrix M ) ;
  76
  77         SparseMatrixPetsc& operator+= (const SparseMatrixPetsc& SparseMatrixPetsc) ;
  78
  79         SparseMatrixPetsc& operator-= (const SparseMatrixPetsc& SparseMatrixPetsc) ;
  80
  81         SparseMatrixPetsc& operator*= (double value) ;
  82
  83         SparseMatrixPetsc& operator/= (double value) ;
  84
  85         SparseMatrixPetsc& operator*= (const SparseMatrixPetsc& matrix) ;
  86
  87         Vector operator* (const Vector& vector) const ;
  88
  89         const SparseMatrixPetsc& operator= ( const SparseMatrixPetsc& SparseMatrixPetsc ) ;
  90
  91         friend SparseMatrixPetsc operator+ (const SparseMatrixPetsc& SparseMatrixPetsc1, const SparseMatrixPetsc& SparseMatrixPetsc2);
  92
  93         friend SparseMatrixPetsc operator- (const SparseMatrixPetsc& SparseMatrixPetsc1, const SparseMatrixPetsc& SparseMatrixPetsc2);
  94
  95         friend SparseMatrixPetsc operator* (double value , const SparseMatrixPetsc& SparseMatrixPetsc ) ;
  96
  97         friend SparseMatrixPetsc operator* (const SparseMatrixPetsc& SparseMatrixPetsc, double value ) ;
  98
  99         friend SparseMatrixPetsc operator/ (const SparseMatrixPetsc& SparseMatrixPetsc, double value) ;
 100
 101         friend SparseMatrixPetsc operator*(const SparseMatrixPetsc& M, const SparseMatrixPetsc& N) ;
 102
 103         void viewMatrix() const ;
 104         double getMatrixCoeff(int i, int j) const;
 105
 106         bool containsPetscMatrix() const;
 107         Mat getPetscMatrix() const;
 108
 109     void diagonalShift(double lambda);
 110     void zeroEntries();//sets the matrix coefficients to zero
 111
 112     std::vector< double > getEigenvalues(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
 113     std::vector< Vector > getEigenvectors(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
 114     MEDCoupling::DataArrayDouble * getEigenvectorsDataArrayDouble(int nev, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
 115     std::vector< double > getSingularValues(int nsv, SVDWhich which=SVD_SMALLEST, double tol=1e-6) const;
 116     double getConditionNumber(bool isSingular=false, double tol=1e-6) const;
 117
 118     bool isSymmetric(double tol=1.e-6) const ;
 119
 120     void  leftDiagonalScale(Vector v);
 121     void rightDiagonalScale(Vector v);
 122
 123 private:
 124         Mat _mat;
 125
 126         int _numberOfNonZeros ;//The maximum number of nonzeros coefficients per line (or an upper bound)
 127
 128         int computeSpectrum(int nev, double ** valP, double ***vecP, EPSWhich which=EPS_SMALLEST_MAGNITUDE, double tol=1e-6) const;
 129         int computeSVD     (int nsv, double ** valS,                 SVDWhich which=SVD_SMALLEST          , double tol=1e-6) const;
 130
 131         Vector vecToVector(const Vec& vec) const ;
 132         Vec vectorToVec( const Vector& myVector ) const ;
 133 };
 134
 135
 136 #endif /* SOURCE_DIRECTORY__BASE_INC_SparseMatrixPetsc_HXX_ */