1#ifndef DUNE_FEM_PETSCCOMMON_HH
2#define DUNE_FEM_PETSCCOMMON_HH
13#include <dune/common/parallel/communication.hh>
14#include <dune/common/stdstreams.hh>
15#include <dune/common/exceptions.hh>
23#define PETSC_HAVE_BROKEN_RECURSIVE_MACRO
42 class Exception :
public Dune::Exception {};
47 typedef ::PetscErrorCode ErrorCode;
55 inline static ErrorCode ErrorCheckHelper ( ErrorCode errorCode ) { CHKERRQ( errorCode );
return 0; }
57 inline ErrorCode ErrorHandler ( MPI_Comm comm,
int line,
const char *function,
const char *file,
58#
if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 5
61 ErrorCode errorCode, PetscErrorType p,
const char *message,
void *context )
63 std::ostringstream msgout;
64 msgout <<
"PETSc Error in the PETSc function '" << function <<
"' at " << file <<
":" << line <<
":";
67 PetscErrorMessage( errorCode, &text, PETSC_NULL );
69 msgout <<
" '" << text <<
"'. Error message: '" << message <<
"'";
71 msgout <<
" Unable to retrieve error text";
73 DUNE_THROW( Exception, msgout.str() );
78 inline static void ErrorCheck ( ErrorCode errorCode )
82 DUNE_THROW( Exception,
"Generic PETSC exception" );
90 inline bool initialize(
const bool verbose,
int &argc,
char **&args,
const char file[] = 0 ,
const char help[] = 0,
bool ownHandler =
true )
92 bool wasInitializedHere = false ;
93 PetscBool petscInitialized = PETSC_FALSE;
96 ::PetscInitialized( &petscInitialized );
98 if( ! petscInitialized )
100 ::PetscInitialize( &argc, &args, file, help );
101 wasInitializedHere =
true;
109 dverb <<
"INFORMATION: Setting up an own error handler for PETSc errors. If you want the default PETSc handler,\n"
110 <<
"INFORMATION: set the last argument of Dune::Petsc::initialize(...) to 'false'.\n";
112 ::PetscPushErrorHandler( &ErrorHandler, 0 );
114 return wasInitializedHere;
120 inline void finalize ()
123 ::PetscPopErrorHandler();
124 PetscBool finalized = PETSC_FALSE;
125 ErrorCheck( ::PetscFinalized( &finalized ) );
133 template <
class Comm>
134 inline auto castToPetscComm(
const Comm& comm )
138 if constexpr ( std::is_same< Dune::No_Comm, Comm > :: value ||
139 std::is_same< Dune::Communication< No_Comm >, Comm >::value )
141 return PETSC_COMM_SELF;
155 template <
class Comm>
156 inline static void KSPCreate (
const Comm& comm, KSP *inksp )
158 ErrorCheck( ::KSPCreate( castToPetscComm( comm ), inksp ) );
161 inline static void KSPDestroy ( KSP *ksp ) {
162#if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 2
163 ErrorCheck( ::KSPDestroy( *ksp ) );
165 ErrorCheck( ::KSPDestroy( ksp ) );
168 inline static void KSPGetPC ( KSP ksp, PC *pc ) { ErrorCheck( ::KSPGetPC( ksp, pc ) ); }
169 inline static void KSPSetFromOptions ( KSP ksp ) { ErrorCheck( ::KSPSetFromOptions( ksp ) ); }
170 inline static void KSPSetUp ( KSP ksp ) { ErrorCheck( ::KSPSetUp( ksp ) ); }
171 inline static void KSPSetType ( KSP ksp,
const KSPType type ) { ErrorCheck( ::KSPSetType( ksp, type ) ); }
172 inline static void KSPGMRESSetRestart ( KSP ksp, PetscInt restart ) { ErrorCheck( ::KSPGMRESSetRestart( ksp, restart ) ); }
174 template <
class Comm>
175 inline static void KSPView (
const Comm& comm,
179 ErrorCheck( ::KSPView( ksp, viewer ) );
182 template <
class Comm>
183 inline static void KSPView (
const Comm& comm,
186 KSPView( comm, ksp, PETSC_VIEWER_STDOUT_( castToPetscComm(comm) ) );
189 inline static void KSPMonitorSet (KSP ksp, PetscErrorCode (*monitor)(KSP,PetscInt,PetscReal,
void*),
190#
if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 2
191 void *mctx,PetscErrorCode (*monitordestroy)(
void*)
193 void *mctx,PetscErrorCode (*monitordestroy)(
void**)
197 ErrorCheck( ::KSPMonitorSet( ksp, monitor, mctx, monitordestroy ) );
199 inline static void KSPGetIterationNumber( KSP ksp, PetscInt* its )
200 { ErrorCheck( ::KSPGetIterationNumber( ksp, its ) ); }
201 inline static void KSPGetConvergedReason(KSP ksp, KSPConvergedReason *reason)
202 { ErrorCheck( ::KSPGetConvergedReason( ksp, reason ) ); }
205#if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 5
206 inline static void KSPSetOperators (KSP ksp, Mat Amat, Mat Pmat, MatStructure flag ) { ErrorCheck( ::KSPSetOperators( ksp, Amat, Pmat, flag ) ); }
208 inline static void KSPSetOperators (KSP ksp, Mat Amat, Mat Pmat ) { ErrorCheck( ::KSPSetOperators( ksp, Amat, Pmat ) ); }
210 inline static void KSPSetTolerances ( KSP ksp, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt maxits )
211 { ErrorCheck( ::KSPSetTolerances( ksp, rtol, abstol, dtol, maxits ) ); }
212 inline static void KSPSetInitialGuessNonzero( KSP ksp, PetscBool flg ) { ErrorCheck( ::KSPSetInitialGuessNonzero( ksp, flg ) ); };
213 inline static void KSPSolve ( KSP ksp, Vec b, Vec x ) { ErrorCheck( ::KSPSolve( ksp, b, x ) ); }
214 inline static void KSPSetPC ( KSP ksp, PC pc ) { ErrorCheck( ::KSPSetPC( ksp, pc ) ); }
217 template <
class Comm>
218 inline static void PCCreate (
const Comm& comm, PC* pc)
220 ErrorCheck( ::PCCreate( castToPetscComm( comm ), pc ) );
223 inline static void PCDestroy ( PC* pc) {
224#if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 2
225 ErrorCheck( ::PCDestroy( *pc ) );
227 ErrorCheck( ::PCDestroy( pc ) );
230 inline static void PCSetType ( PC pc,
const PCType type ) { ErrorCheck( ::PCSetType( pc, type ) ); }
231 inline static void PCSetFromOptions ( PC pc ) { ErrorCheck( ::PCSetFromOptions( pc ) ); }
232 inline static void PCSetUp ( PC pc ) { ErrorCheck( ::PCSetUp( pc ) ); }
233 inline static void PCFactorSetLevels( PC pc, PetscInt level ) { ErrorCheck( ::PCFactorSetLevels( pc, level ) ); }
234 inline static void PCSORSetOmega( PC pc, PetscReal omega ) { ErrorCheck( ::PCSORSetOmega( pc, omega ) ); }
235#if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 9
236 inline static void PCFactorSetMatSolverPackage( PC pc,
const MatSolverPackage type )
238 ErrorCheck( ::PCFactorSetMatSolverPackage( pc, type ) );
241 inline static void PCFactorSetMatSolverPackage( PC pc,
const MatSolverType type )
243 ErrorCheck( ::PCFactorSetMatSolverType( pc, type ) );
246 inline static void PCHYPRESetType( PC pc,
const char* type )
248 ErrorCheck( ::PCHYPRESetType( pc, type ) );
252 inline static void MatAssemblyBegin ( Mat mat, MatAssemblyType type ) { ErrorCheck( ::MatAssemblyBegin( mat, type ) ); }
253 inline static void MatAssemblyEnd ( Mat mat, MatAssemblyType type ) { ErrorCheck( ::MatAssemblyEnd( mat, type ) ); }
254 inline static void MatAssembled( Mat mat, PetscBool* assembled ) { ErrorCheck( ::MatAssembled( mat, assembled ) ); }
256 template <
class Comm>
257 inline static void MatCreate (
const Comm& comm, Mat *A )
259 ErrorCheck( ::MatCreate( castToPetscComm( comm ), A) );
262 template <
class Comm>
263 inline static void MatCreateBlockMat (
const Comm& comm,
265 PetscInt m, PetscInt n, PetscInt bs, PetscInt nz, PetscInt* nnz )
267 ErrorCheck( ::MatCreateBlockMat( castToPetscComm( comm ), n, m, bs, nz, nnz, A) );
269 inline static void MatDestroy ( Mat *A ) {
270 #if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 2
271 ErrorCheck( ::MatDestroy( *A ) );
273 ErrorCheck( ::MatDestroy( A ) );
276 inline static void MatSetUp( Mat mat )
278 ErrorCheck( ::MatSetUp(mat));
280 inline static void MatSetUp( Mat mat, PetscInt bs, PetscInt nz )
284 ErrorCheck( ::MatSeqAIJSetPreallocation(mat,nz,PETSC_NULL) );
285 ErrorCheck( ::MatMPIAIJSetPreallocation(mat,nz,PETSC_NULL,nz/2,PETSC_NULL) );
289 ErrorCheck( ::MatSeqBAIJSetPreallocation(mat,bs,nz,PETSC_NULL) );
290 ErrorCheck( ::MatMPIBAIJSetPreallocation(mat,bs,nz,PETSC_NULL,nz/2,PETSC_NULL) );
293 ErrorCheck( ::MatSetOption(mat, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE) );
297 ErrorCheck( ::MatSetUp(mat));
299 inline static void MatSetUp( Mat mat, PetscInt bs,
const PetscInt *d_nnz,
const PetscInt *o_nnz )
303 ErrorCheck( ::MatSeqAIJSetPreallocation(mat,0,d_nnz ) );
304 ErrorCheck( ::MatMPIAIJSetPreallocation(mat,0,d_nnz,5,o_nnz) );
308 ErrorCheck( ::MatSeqBAIJSetPreallocation(mat,bs,0,d_nnz ) );
309 ErrorCheck( ::MatMPIBAIJSetPreallocation(mat,bs,0,d_nnz,5,PETSC_NULL) );
312 ErrorCheck( ::MatSetOption(mat, MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE) );
313 ErrorCheck( ::MatSetUp(mat));
316 inline static void MatGetOwnershipRange ( Mat mat, PetscInt *m, PetscInt* n ) { ErrorCheck( ::MatGetOwnershipRange( mat, m, n ) ); }
317 inline static void MatGetSize ( Mat mat, PetscInt *m, PetscInt* n ) { ErrorCheck( ::MatGetSize( mat, m, n ) ); }
318 inline static void MatMult ( Mat mat, Vec x, Vec y ) { ErrorCheck( ::MatMult( mat, x, y ) ); }
319 inline static void MatSetBlockSize ( Mat A, PetscInt bs ) { ErrorCheck( ::MatSetBlockSize( A, bs ) ); }
320 inline static void MatSetSizes ( Mat A, PetscInt m, PetscInt n, PetscInt M, PetscInt N ) { ErrorCheck( ::MatSetSizes( A, m, n, M, N ) ); }
321 inline static void MatSetFromOptions ( Mat B ) { ErrorCheck( ::MatSetFromOptions( B ) ); }
322 inline static void MatSetType ( Mat mat,
const MatType matype ) { ErrorCheck( ::MatSetType( mat, matype ) ); }
324 inline static void MatSetValue ( Mat v, PetscInt i, PetscInt j, PetscScalar va, InsertMode mode )
330 ErrorCheck( ::MatSetValue( v, i, j, va, mode ) );
334 inline static void MatSetValues ( Mat mat, PetscInt m,
const PetscInt idxm[], PetscInt n,
const PetscInt idxn[],
const PetscScalar v[], InsertMode addv )
340 ErrorCheck( ::MatSetValues( mat, m, idxm, n, idxn, v, addv ) );
344 inline static void MatSetValuesBlocked ( Mat mat, PetscInt m,
const PetscInt idxm[], PetscInt n,
const PetscInt idxn[],
const PetscScalar v[], InsertMode addv )
350 ErrorCheck( ::MatSetValuesBlocked( mat, m, idxm, n, idxn, v, addv ) );
354 inline static void MatGetValues ( Mat mat, PetscInt m,
const PetscInt idxm[], PetscInt n,
const PetscInt idxn[], PetscScalar v[] )
355 { ErrorCheck( ::MatGetValues( mat, m, idxm, n, idxn, v ) ); }
357 inline static void MatZeroRows ( Mat mat, PetscInt m,
const PetscInt idxm[],
const PetscScalar v )
363 ErrorCheck( ::MatZeroRows( mat, m, idxm, v, 0, 0 ) );
367 inline static void MatView ( Mat mat, PetscViewer viewer ) { ErrorCheck( ::MatView( mat, viewer ) ); }
368 inline static void MatZeroEntries ( Mat mat )
374 ErrorCheck( ::MatZeroEntries( mat ) );
378 inline static void PetscBarrier ( PetscObject obj ) { ErrorCheck( ::PetscBarrier( obj ) ); }
379 inline static void PetscFinalize () { ErrorCheck( ::PetscFinalize() ); }
380 inline static void PetscInitialize(
int *argc,
char ***args,
const char file[],
const char help[] ) { ErrorCheck( ::PetscInitialize( argc, args, file, help ) ); }
381 inline static void PetscViewerASCIIOpen ( MPI_Comm comm,
const char name[], PetscViewer *lab ) { ErrorCheck( ::PetscViewerASCIIOpen( comm, name, lab ) ); }
382 inline static void PetscViewerDestroy ( PetscViewer *viewer ) {
383 #if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 2
384 ErrorCheck( ::PetscViewerDestroy( *viewer ) );
386 ErrorCheck( ::PetscViewerDestroy( viewer ) );
389 inline static void PetscViewerSetFormat ( PetscViewer viewer, PetscViewerFormat format )
391 #if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 5
392 ErrorCheck( ::PetscViewerSetFormat( viewer, format ) );
394 ErrorCheck( ::PetscViewerPushFormat( viewer, format ) );
397 inline static void VecAssemblyBegin ( Vec vec ) { ErrorCheck( ::VecAssemblyBegin( vec ) ); }
398 inline static void VecAssemblyEnd ( Vec vec ) { ErrorCheck( ::VecAssemblyEnd( vec ) ); }
399 inline static void VecAXPY ( Vec y, PetscScalar alpha, Vec x) { ErrorCheck( ::VecAXPY( y, alpha, x ) ); }
400 inline static void VecCopy ( Vec x, Vec y ) { ErrorCheck( ::VecCopy( x, y ) ); }
402 template <
class Comm>
403 inline static void VecCreate (
const Comm& comm, Vec *vec )
405 ErrorCheck( ::VecCreate( castToPetscComm( comm ), vec ) );
408 template <
class Comm>
409 inline static void VecCreateGhost (
const Comm& comm, PetscInt n, PetscInt N, PetscInt nghost,
const PetscInt ghosts[], Vec *vv )
410 { ErrorCheck( ::VecCreateGhost( castToPetscComm( comm ), n, N, nghost, ghosts, vv ) ); }
412 template <
class Comm>
413 inline static void VecCreateGhostBlock (
const Comm& comm, PetscInt bs, PetscInt n, PetscInt N, PetscInt nghost,
const PetscInt ghosts[], Vec *vv )
415#if PETSC_VERSION_MAJOR < 3 || (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR <= 2)
416 std::unique_ptr< PetscInt[] > ghostsCopy(
new PetscInt[ nghost * bs ] );
417 for( PetscInt i = 0; i < nghost; ++i )
419 for( PetscInt j = 0; j < bs; ++j )
420 ghostsCopy[ i*bs + j ] = ghosts[ i ]*bs + j;
422 VecCreateGhost( n, N, nghost * bs, ghostsCopy.get(), vv );
424 ErrorCheck( ::VecCreateGhostBlock( castToPetscComm( comm ), bs, n, N, nghost, ghosts, vv ) );
428 inline static void VecDestroy ( Vec *v ) {
429 #if PETSC_VERSION_MAJOR <= 3 && PETSC_VERSION_MINOR < 2
430 ErrorCheck( ::VecDestroy( *v ) );
432 ErrorCheck( ::VecDestroy( v ) );
435 inline static void VecDot ( Vec x, Vec y, PetscScalar *val ) { ErrorCheck( ::VecDot( x, y, val ) ); }
436 inline static void VecDuplicate ( Vec v, Vec *newv ) { ErrorCheck( ::VecDuplicate( v, newv ) ); }
437 inline static void VecGetBlockSize ( Vec x, PetscInt* bs ) { ErrorCheck( ::VecGetBlockSize( x, bs ) ); }
438 inline static void VecGetLocalSize ( Vec x, PetscInt *size ) { ErrorCheck( ::VecGetLocalSize( x, size ) ); }
439 inline static void VecGetOwnershipRange ( Vec x, PetscInt *low, PetscInt *high ) { ErrorCheck( ::VecGetOwnershipRange( x, low, high ) ); }
440 inline static void VecGetSize ( Vec x, PetscInt *size ) { ErrorCheck( ::VecGetSize( x, size ) ); }
441 inline static void VecGetValues ( Vec x, PetscInt ni,
const PetscInt ix[], PetscScalar y[] ) { ErrorCheck( ::VecGetValues( x, ni, ix, y ) ); }
442 inline static void VecGhostGetLocalForm ( Vec g, Vec *l ) { ErrorCheck( ::VecGhostGetLocalForm( g, l ) ); }
443 inline static void VecGhostRestoreLocalForm ( Vec g, Vec *l ) { ErrorCheck( ::VecGhostRestoreLocalForm( g, l ) ); }
444 inline static void VecGhostUpdateBegin ( Vec g, InsertMode insertmode, ScatterMode scattermode ) { ErrorCheck( ::VecGhostUpdateBegin( g, insertmode, scattermode ) ); }
445 inline static void VecGhostUpdateEnd ( Vec g, InsertMode insertmode, ScatterMode scattermode ) { ErrorCheck( ::VecGhostUpdateEnd( g, insertmode, scattermode ) ); }
446 inline static void VecNorm ( Vec x, NormType type, PetscReal *val ) { ErrorCheck( ::VecNorm( x, type, val ) ); }
447 inline static void VecScale ( Vec x, PetscScalar alpha ) { ErrorCheck( ::VecScale( x, alpha ) ); }
448 inline static void VecSet ( Vec x, PetscScalar alpha ) { ErrorCheck( ::VecSet( x, alpha ) ); }
449 inline static void VecSetBlockSize ( Vec x, PetscInt bs ) { ErrorCheck( ::VecSetBlockSize( x, bs ) ); }
450 inline static void VecSetFromOptions ( Vec vec ) { ErrorCheck( ::VecSetFromOptions( vec ) ); }
451 inline static void VecSetType ( Vec vec,
const VecType method ) { ErrorCheck( ::VecSetType( vec, method ) ); }
452 inline static void VecSetSizes ( Vec v, PetscInt n, PetscInt N ) { ErrorCheck( ::VecSetSizes( v, n, N ) ); }
453 inline static void VecSetValue ( Vec v,
int row, PetscScalar value, InsertMode mode ) { ErrorCheck( ::VecSetValue( v, row, value, mode ) ); }
454 inline static void VecSetValuesBlocked ( Vec v, PetscInt ni,
const PetscInt xi[],
const PetscScalar values[], InsertMode mode ) { ErrorCheck( ::VecSetValuesBlocked( v, ni, xi, values, mode ) ); }
455 inline static void VecView ( Vec vec, PetscViewer viewer ) { ErrorCheck( ::VecView( vec, viewer ) ); }
Definition: bindguard.hh:11