30 #ifndef _BACKENDPETSC_HPP_
31 #define _BACKENDPETSC_HPP_
33 #include <boost/program_options/variables_map.hpp>
45 namespace po = boost::program_options;
51 #if defined( FEELPP_HAS_PETSC_H )
58 class BackendPetsc :
public Backend<T>
60 typedef Backend<T> super;
64 typedef typename super::value_type value_type;
67 typedef typename super::sparse_matrix_type sparse_matrix_type;
68 typedef typename super::sparse_matrix_ptrtype sparse_matrix_ptrtype;
69 typedef MatrixPetsc<value_type> petsc_sparse_matrix_type;
70 typedef boost::shared_ptr<sparse_matrix_type> petsc_sparse_matrix_ptrtype;
71 typedef MatrixPetscMPI<value_type> petscMPI_sparse_matrix_type;
74 typedef typename sparse_matrix_type::graph_type graph_type;
75 typedef typename sparse_matrix_type::graph_ptrtype graph_ptrtype;
78 typedef typename super::vector_type vector_type;
79 typedef typename super::vector_ptrtype vector_ptrtype;
80 typedef VectorPetsc<value_type> petsc_vector_type;
81 typedef boost::shared_ptr<vector_type> petsc_vector_ptrtype;
82 typedef VectorPetscMPI<value_type> petscMPI_vector_type;
84 typedef typename super::solve_return_type solve_return_type;
85 typedef typename super::nl_solve_return_type nl_solve_return_type;
87 typedef typename super::datamap_type datamap_type;
88 typedef typename super::datamap_ptrtype datamap_ptrtype;
92 BackendPetsc( WorldComm
const& worldComm=Environment::worldComm() )
95 M_solver_petsc( worldComm ),
96 M_nl_solver_petsc( worldComm )
99 BackendPetsc( po::variables_map
const& vm, std::string
const& prefix =
"",
100 WorldComm
const& worldComm=Environment::worldComm() )
102 super( vm, prefix, worldComm ),
103 M_solver_petsc( vm, worldComm ),
104 M_nl_solver_petsc( worldComm )
106 std::string _prefix = prefix;
108 if ( !_prefix.empty() )
114 template<
typename DomainSpace,
typename DualImageSpace>
115 static sparse_matrix_ptrtype newMatrix( DomainSpace
const& Xh,
116 DualImageSpace
const& Yh,
119 auto s = stencil( _test=Yh,_trial=Xh );
121 sparse_matrix_ptrtype mat;
122 if ( Yh->worldComm().globalSize()>1 )
123 mat = sparse_matrix_ptrtype(
new petscMPI_sparse_matrix_type( Yh->dof(),Xh->dof() ) );
125 mat = sparse_matrix_ptrtype(
new petsc_sparse_matrix_type( Yh->dof(),Xh->dof() ) );
127 mat->setMatrixProperties( matrix_properties );
128 mat->init( Yh->nDof(), Xh->nDof(),
129 Yh->nLocalDofWithoutGhost(), Xh->nLocalDofWithoutGhost(),
133 auto nSpace = DomainSpace::nSpaces;
135 std::vector < std::vector<int> > is( nSpace );
140 auto result = boost::make_tuple( cptSpaces,is );
141 boost::fusion::fold( Xh->functionSpaces(), result, computeNDofForEachSpace() );
143 for ( uint i = 0; i<nSpace; i++ )
153 sparse_matrix_ptrtype
162 sparse_matrix_ptrtype mat;
164 if ( this->comm().globalSize()>1 )
165 mat = sparse_matrix_ptrtype(
new petscMPI_sparse_matrix_type );
167 mat = sparse_matrix_ptrtype(
new petsc_sparse_matrix_type );
169 mat->setMatrixProperties( matrix_properties );
170 mat->init( m,n,m_l,n_l,nnz,noz );
174 sparse_matrix_ptrtype
175 newMatrix( datamap_ptrtype
const& domainmap,
176 datamap_ptrtype
const& imagemap,
180 sparse_matrix_ptrtype mat;
182 if ( imagemap->worldComm().globalSize()>1 )
183 mat = sparse_matrix_ptrtype(
new petscMPI_sparse_matrix_type( imagemap,domainmap,imagemap->worldComm() ) );
185 mat = sparse_matrix_ptrtype(
new petsc_sparse_matrix_type( imagemap,domainmap,imagemap->worldComm() ) );
187 mat->setMatrixProperties( matrix_properties );
191 mat->init( imagemap->nDof(), domainmap->nDof(),
192 imagemap->nLocalDofWithoutGhost(), domainmap->nLocalDofWithoutGhost() );
198 sparse_matrix_ptrtype
203 graph_ptrtype
const & graph,
206 sparse_matrix_ptrtype mat;
207 auto const& mapGraphRow = graph->mapRowPtr();
208 auto const& mapGraphCol = graph->mapColPtr();
209 if ( this->comm().globalSize()>1 )
210 mat = sparse_matrix_ptrtype(
new petscMPI_sparse_matrix_type( mapGraphRow,mapGraphCol,mapGraphRow->worldComm() ) );
212 mat = sparse_matrix_ptrtype(
new petsc_sparse_matrix_type( mapGraphRow,mapGraphCol,mapGraphRow->worldComm() ) ) ;
214 mat->setMatrixProperties( matrix_properties );
216 mat->init( mapGraphRow->nDof(), mapGraphCol->nDof(),
217 mapGraphRow->nLocalDofWithoutGhost(), mapGraphCol->nLocalDofWithoutGhost(),
224 sparse_matrix_ptrtype
225 newZeroMatrix( datamap_ptrtype
const& domainmap,
226 datamap_ptrtype
const& imagemap )
228 graph_ptrtype sparsity_graph(
new graph_type( imagemap, domainmap ) );
229 sparsity_graph->zero();
230 sparsity_graph->close();
232 sparse_matrix_ptrtype mat;
233 if ( imagemap->worldComm().globalSize()>1 )
234 mat = sparse_matrix_ptrtype(
new petscMPI_sparse_matrix_type( imagemap,domainmap,imagemap->worldComm() ) );
236 mat = sparse_matrix_ptrtype(
new petsc_sparse_matrix_type( imagemap,domainmap,imagemap->worldComm() ) );
238 mat->init( imagemap->nDof(), domainmap->nDof(),
239 imagemap->nLocalDofWithoutGhost(), domainmap->nLocalDofWithoutGhost(),
245 sparse_matrix_ptrtype
251 graph_ptrtype sparsity_graph(
new graph_type( 0,0,m_l-1,0,n_l-1 ) );
252 sparsity_graph->zero();
253 sparsity_graph->close();
255 sparse_matrix_ptrtype mat;
257 if ( this->comm().globalSize()>1 )
258 mat = sparse_matrix_ptrtype(
new petscMPI_sparse_matrix_type );
260 mat = sparse_matrix_ptrtype(
new petsc_sparse_matrix_type );
263 mat->init( m, n, m_l, n_l, sparsity_graph );
268 template<
typename SpaceT>
269 static vector_ptrtype newVector( SpaceT
const& space )
271 if ( space->worldComm().globalSize()>1 )
272 return vector_ptrtype(
new petscMPI_vector_type( space->dof() ) );
274 return vector_ptrtype(
new petsc_vector_type( space->dof() ) );
277 vector_ptrtype newVector( datamap_ptrtype
const& dm )
279 if ( this->comm().globalSize()>1 )
return vector_ptrtype(
new petscMPI_vector_type( dm ) );
281 else return vector_ptrtype(
new petsc_vector_type( dm ) );
286 return vector_ptrtype(
new petsc_vector_type( n, n_local ) );
290 void set_symmetric(
bool ) {}
293 template <
class Vector>
294 static void applyMatrix( sparse_matrix_type
const& A,
301 void prod( sparse_matrix_type
const& A,
302 vector_type
const& x,
303 vector_type& b )
const
306 petsc_sparse_matrix_type
const& _A =
dynamic_cast<petsc_sparse_matrix_type const&
>( A );
307 petsc_vector_type
const& _x =
dynamic_cast<petsc_vector_type const&
>( x );
308 petsc_vector_type
const& _b =
dynamic_cast<petsc_vector_type const&
>( b );
309 if ( _A.mapCol().worldComm().globalSize() == x.map().worldComm().globalSize() )
312 ierr = MatMult( _A.mat(), _x.vec(), _b.vec() );
313 CHKERRABORT( _A.comm().globalComm(),ierr );
318 auto x_convert = petscMPI_vector_type(_A.mapColPtr());
319 x_convert.duplicateFromOtherPartition(x);
321 ierr = MatMult( _A.mat(), x_convert.vec(), _b.vec() );
322 CHKERRABORT( _A.comm().globalComm(),ierr );
328 solve_return_type solve( sparse_matrix_type
const& A,
330 vector_type
const& b );
332 solve_return_type solve( sparse_matrix_ptrtype
const& A,
333 sparse_matrix_ptrtype
const& B,
335 vector_ptrtype
const& b );
337 template <
class Vector>
338 static value_type dot(
const vector_type& f,
341 value_type result( 0 );
350 SolverLinearPetsc<double> M_solver_petsc;
351 SolverNonLinearPetsc<double> M_nl_solver_petsc;
357 BackendPetsc<T>::~BackendPetsc()
363 BackendPetsc<T>::clear()
365 LOG(INFO) <<
"Deleting linear solver petsc";
366 M_solver_petsc.clear();
367 LOG(INFO) <<
"Deleting non linear solver petsc";
368 M_nl_solver_petsc.clear();
369 LOG(INFO) <<
"Deleting backend petsc";
376 typename BackendPetsc<T>::solve_return_type
377 BackendPetsc<T>::solve( sparse_matrix_ptrtype
const& A,
378 sparse_matrix_ptrtype
const& B,
380 vector_ptrtype
const& b )
382 M_solver_petsc.setPrefix( this->prefix() );
383 M_solver_petsc.setPreconditionerType( this->pcEnumType() );
384 M_solver_petsc.setSolverType( this->kspEnumType() );
385 if (!M_solver_petsc.initialized())
386 M_solver_petsc.attachPreconditioner( this->M_preconditioner );
387 M_solver_petsc.setConstantNullSpace( this->hasConstantNullSpace() );
388 M_solver_petsc.setFieldSplitType( this->fieldSplitEnumType() );
389 M_solver_petsc.setTolerances( _rtolerance=this->rTolerance(),
390 _atolerance=this->aTolerance(),
391 _dtolerance=this->dTolerance(),
392 _maxit = this->maxIterations() );
393 M_solver_petsc.setPrecMatrixStructure( this->precMatrixStructure() );
394 M_solver_petsc.setMatSolverPackageType( this->matSolverPackageEnumType() );
395 M_solver_petsc.setShowKSPMonitor( this->showKSPMonitor() );
396 M_solver_petsc.setShowKSPConvergedReason( this->showKSPConvergedReason() );
398 auto res = M_solver_petsc.solve( *A, *B, *x, *b, this->rTolerance(), this->maxIterations(), this->transpose() );
399 DVLOG(2) <<
"[BackendPetsc::solve] number of iterations : " << res.template get<1>() <<
"\n";
400 DVLOG(2) <<
"[BackendPetsc::solve] residual : " << res.template get<2>() <<
"\n";
402 if ( !res.template get<0>() )
403 LOG(ERROR) <<
"Backend " << this->prefix() <<
" : linear solver failed to converge" << std::endl;
410 typename BackendPetsc<T>::solve_return_type
411 BackendPetsc<T>::solve( sparse_matrix_type
const& A,
413 vector_type
const& b )
415 M_solver_petsc.setPrefix( this->prefix() );
416 M_solver_petsc.setPreconditionerType( this->pcEnumType() );
417 M_solver_petsc.setSolverType( this->kspEnumType() );
418 if (!M_solver_petsc.initialized())
419 M_solver_petsc.attachPreconditioner( this->M_preconditioner );
420 M_solver_petsc.setConstantNullSpace( this->hasConstantNullSpace() );
421 M_solver_petsc.setFieldSplitType( this->fieldSplitEnumType() );
422 M_solver_petsc.setTolerances( _rtolerance=this->rTolerance(),
423 _atolerance=this->aTolerance(),
424 _dtolerance=this->dTolerance(),
425 _maxit = this->maxIterations() );
426 M_solver_petsc.setPrecMatrixStructure( this->precMatrixStructure() );
427 M_solver_petsc.setMatSolverPackageType( this->matSolverPackageEnumType() );
428 M_solver_petsc.setShowKSPMonitor( this->showKSPMonitor() );
429 M_solver_petsc.setShowKSPConvergedReason( this->showKSPConvergedReason() );
431 auto res = M_solver_petsc.solve( A, x, b, this->rTolerance(), this->maxIterations() );
432 DVLOG(2) <<
"[BackendPetsc::solve] number of iterations : " << res.template get<1>() <<
"\n";
433 DVLOG(2) <<
"[BackendPetsc::solve] residual : " << res.template get<2>() <<
"\n";
435 if ( !res.template get<0>() )
436 LOG(ERROR) <<
"Backend " << this->prefix() <<
" : linear solver failed to converge" << std::endl;
443 #endif // FEELPP_HAS_PETSC_H
Definition: feelalg/enums.hpp:67
size_t size_type
Indices (starting from 0)
Definition: feelcore/feel.hpp:319