30 #ifndef __MatrixPetsc_H
31 #define __MatrixPetsc_H 1
33 #include <feel/feelconfig.h>
35 #if defined(FEELPP_HAS_PETSC_H)
49 #ifndef USE_COMPLEX_NUMBERS
51 # include <petscversion.h>
52 # if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR <= 1)
53 # include <petscsles.h>
55 # include <petscksp.h>
59 # include <petscversion.h>
60 # if (PETSC_VERSION_MAJOR == 2) && (PETSC_VERSION_MINOR <= 1)
61 # include <petscsles.h>
63 # include <petscksp.h>
71 template<
typename T>
class VectorPetsc;
86 class MatrixPetsc :
public MatrixSparse<T>
88 typedef MatrixSparse<T> super;
94 static const bool is_row_major =
true;
102 typedef typename super::value_type value_type;
103 typedef typename super::real_type real_type;
105 typedef std::vector<std::set<size_type> > pattern_type;
107 typedef typename super::graph_type graph_type;
108 typedef typename super::graph_ptrtype graph_ptrtype;
110 typedef typename super::datamap_type datamap_type;
111 typedef typename super::datamap_ptrtype datamap_ptrtype;
136 MatrixPetsc( datamap_ptrtype
const& dmRow, datamap_ptrtype
const& dmCol, WorldComm
const& worldComm=Environment::worldComm() );
146 MatrixPetsc ( Mat m );
147 MatrixPetsc ( Mat m, datamap_ptrtype
const& dmRow, datamap_ptrtype
const& dmCol );
148 MatrixPetsc ( MatrixSparse<value_type>
const& M, IS& isrow, IS& iscol );
149 MatrixPetsc ( MatrixSparse<value_type>
const& M, std::vector<int>
const& rowIndex, std::vector<int>
const& colIndex );
174 value_type operator () (
const size_type i,
181 MatrixPetsc&
operator=( MatrixSparse<value_type>
const& M );
248 graph_ptrtype
const& graph );
253 void setIndexSplit( std::vector< std::vector<size_type> >
const &indexSplit );
260 this->setInitialized(
false );
264 void fill( pattern_type
const& patt )
307 real_type l1Norm ()
const;
319 real_type linftyNorm ()
const;
329 const value_type& value );
341 const value_type& value );
349 void addMatrix (
const ublas::matrix<value_type> &dm,
350 const std::vector<size_type> &rows,
351 const std::vector<size_type> &cols );
359 void addMatrix (
int* rows,
int nrows,
360 int* cols,
int ncols,
367 void addMatrix (
const ublas::matrix<value_type> &dm,
368 const std::vector<size_type> &dof_indices )
370 this->addMatrix ( dm, dof_indices, dof_indices );
384 void addMatrix (
const T a, MatrixSparse<T> &X );
391 void matMatMult ( MatrixSparse<T>
const& In, MatrixSparse<T> &Res );
397 void scale( T
const a );
402 void diagonal ( Vector<value_type>& dest )
const;
410 void transpose( MatrixSparse<value_type>& Mt )
const;
415 virtual void symmetricPart( MatrixSparse<value_type>& Ms )
const;
424 FEELPP_ASSERT ( M_mat != NULL ).error(
"null petsc matrix" );
429 FEELPP_ASSERT ( M_mat != NULL ).warn(
"null petsc matrix" );
439 void printMatlab(
const std::string name=
"NULL" )
const;
445 energy( Vector<value_type>
const& __v,
446 Vector<value_type>
const& __u,
447 bool transpose =
false )
const;
456 void zeroRows( std::vector<int>
const& rows, Vector<value_type>
const& values, Vector<value_type>& rhs, Context
const& on_context );
461 void updateBlockMat( boost::shared_ptr<MatrixSparse<T> > m, std::vector<size_type> start_i, std::vector<size_type> start_j );
464 void updatePCFieldSplit( PC & pc );
466 std::vector<IS>
const& petscSplitIS()
const {
return M_petscIS; }
467 std::map<PC*,bool > & mapSplitPC() {
return M_mapPC; }
469 std::vector<PetscInt> ia() {
return M_ia; }
470 std::vector<PetscInt> ja() {
return M_ja; }
478 void zeroEntriesDiagonal();
483 MatrixPetsc( MatrixPetsc
const & );
495 std::vector<IS> M_petscIS;
497 std::map<PC*,bool > M_mapPC;
503 const bool M_destroy_mat_on_exit;
504 std::vector<PetscInt> M_ia,M_ja;
510 class MatrixPetscMPI :
public MatrixPetsc<T>
512 typedef MatrixPetsc<T> super;
516 typedef typename super::graph_type graph_type;
517 typedef typename super::graph_ptrtype graph_ptrtype;
518 typedef typename super::value_type value_type;
520 typedef typename super::datamap_type datamap_type;
521 typedef typename super::datamap_ptrtype datamap_ptrtype;
525 MatrixPetscMPI( datamap_ptrtype
const& dmRow, datamap_ptrtype
const& dmCol, WorldComm
const& worldComm=Environment::worldComm() );
527 MatrixPetscMPI( Mat m, datamap_ptrtype
const& dmRow, datamap_ptrtype
const& dmCol );
548 graph_ptrtype
const& graph );
560 const value_type& value );
564 const value_type& value );
566 void addMatrix(
const ublas::matrix<value_type>& dm,
567 const std::vector<size_type>& rows,
568 const std::vector<size_type>& cols );
570 void addMatrix(
int* rows,
int nrows,
571 int* cols,
int ncols,
574 void addMatrix(
const T a, MatrixSparse<T> &X );
580 void zeroRows( std::vector<int>
const& rows,
581 Vector<value_type>
const& values,
582 Vector<value_type>& rhs,
583 Context
const& on_context );
585 value_type energy( Vector<value_type>
const& __v,
586 Vector<value_type>
const& __u,
587 bool transpose =
false )
const;
591 void addMatrixSameNonZeroPattern(
const T a, MatrixSparse<T> &X );