Interface.hpp

#pragma once
 
#include <petscksp.h>
#include <petsc/private/kspimpl.h>
 
#include <amdis/Output.hpp>
 
namespace AMDiS
{
  namespace PETSc
  {
    // KSP routines
 
    inline PetscErrorCode KSPMonitorCyclic(KSP ksp, PetscInt n, PetscReal rnorm, void* ctx)
    {
      if (n % 10 == 0)
        msg("iteration {:>4}: resid {:<.6e}", n, rnorm);
      return 0;
    }
 
    inline PetscErrorCode KSPMonitorNoisy(KSP ksp, PetscInt n, PetscReal rnorm, void* ctx)
    {
      Vec resid;
      KSPBuildResidual(ksp, NULL, NULL, &resid);
 
      PetscReal truenorm = 0.0;
      VecNorm(resid, NORM_2, &truenorm);
      VecDestroy(&resid);
 
      PetscReal bnorm = 0.0;
      VecNorm(ksp->vec_rhs, NORM_2, &bnorm);
 
      msg("iteration {:>4}: resid {:<.6e}, true resid {:<.2e}, rel. resid {:<.2e}", n, rnorm, truenorm, truenorm/bnorm);
      return 0;
    }
 
    inline PetscErrorCode KSPMonitorTrueResidual(KSP ksp, PetscInt n, PetscReal rnorm, void* ctx)
    {
      assert(ctx != nullptr);
#if (PETSC_VERSION_MINOR >= 15)
      return ::KSPMonitorTrueResidual(ksp, n, rnorm, (PetscViewerAndFormat*)(ctx));
#elif (PETSC_VERSION_MINOR >= 7)
      return ::KSPMonitorTrueResidualNorm(ksp, n, rnorm, (PetscViewerAndFormat*)(ctx));
#else
      return ::KSPMonitorTrueResidualNorm(ksp, n, rnorm, ctx);
#endif
    }
 
    inline PetscErrorCode KSPMonitorResidual(KSP ksp, PetscInt n, PetscReal rnorm, void* ctx)
    {
      assert(ctx != nullptr);
#if (PETSC_VERSION_MINOR >= 15)
      return ::KSPMonitorResidual(ksp, n, rnorm, (PetscViewerAndFormat*)(ctx));
#elif (PETSC_VERSION_MINOR >= 7)
      return ::KSPMonitorDefault(ksp, n, rnorm, (PetscViewerAndFormat*)(ctx));
#else
      return ::KSPMonitorDefault(ksp, n, rnorm, ctx);
#endif
    }
 
    template <class Monitor>
    inline PetscErrorCode KSPMonitorSet(KSP ksp, Monitor monitor)
    {
#if (PETSC_VERSION_MINOR >= 7)
      PetscViewerAndFormat *vf;
      PetscErrorCode ierr;
      ierr = ::PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD,PETSC_VIEWER_DEFAULT,&vf);CHKERRQ(ierr);
      ierr = ::KSPMonitorSet(ksp,(PetscErrorCode (*)(KSP,PetscInt,PetscReal,void*))monitor,vf,(PetscErrorCode (*)(void**))PetscViewerAndFormatDestroy);CHKERRQ(ierr);
      return ierr;
#else
      return ::KSPMonitorSet(ksp, monitor, PETSC_NULL, PETSC_NULL);
#endif
    }
 
    inline PetscErrorCode KSPGetOperators(KSP ksp, Mat *Amat, Mat *Pmat)
    {
#if (PETSC_VERSION_MINOR >= 5)
      return ::KSPGetOperators(ksp, Amat, Pmat);
#else
      return ::KSPGetOperators(ksp, Amat, Pmat, PETSC_NULL);
#endif
    }
 
    inline PetscErrorCode KSPSetOperators(KSP ksp, Mat Amat, Mat Pmat)
    {
#if (PETSC_VERSION_MINOR >= 5)
      return ::KSPSetOperators(ksp, Amat, Pmat);
#else
      return ::KSPSetOperators(ksp, Amat, Pmat, SAME_NONZERO_PATTERN);
#endif
    }
 
    inline PetscErrorCode KSPConvergedReasonView(KSP ksp, PetscViewer viewer)
    {
#if (PETSC_VERSION_MINOR >= 14)
      return ::KSPConvergedReasonView(ksp, viewer);
#else
      return ::KSPReasonView(ksp, viewer);
#endif
    }
 
    // Mat routines
 
    inline PetscErrorCode MatCreateVecs(Mat mat, Vec *right, Vec *left)
    {
#if (PETSC_VERSION_MINOR >= 6)
      return ::MatCreateVecs(mat, right, left);
#else
      return ::MatGetVecs(mat, right, left);
#endif
    }
 
    inline PetscErrorCode MatCreateSubMatrix(Mat mat, IS isrow, IS iscol, MatReuse cll, Mat *newmat)
    {
#if (PETSC_VERSION_MINOR >= 8)
      return ::MatCreateSubMatrix(mat, isrow, iscol, cll, newmat);
#else
      return ::MatGetSubMatrix(mat, isrow, iscol, cll, newmat);
#endif
    }
 
    inline PetscErrorCode MatNullSpaceRemove(MatNullSpace sp, Vec vec)
    {
#if (PETSC_VERSION_MINOR >= 5)
      return ::MatNullSpaceRemove(sp, vec);
#else
      return ::MatNullSpaceRemove(sp, vec, PETSC_NULL);
#endif
    }
 
    // PC routines
 
#if (PETSC_VERSION_MINOR >= 9)
    inline PetscErrorCode PCFactorSetMatSolverType(PC pc, MatSolverType stype)
    {
      return ::PCFactorSetMatSolverType(pc, stype);
    }
#else
    template <class MatSolverType>
    inline PetscErrorCode PCFactorSetMatSolverType(PC pc, MatSolverType stype)
    {
      return ::PCFactorSetMatSolverPackage(pc, stype);
    }
#endif
 
    inline PetscErrorCode PetscOptionsView(PetscViewer viewer)
    {
#if (PETSC_VERSION_MINOR >= 7)
      return ::PetscOptionsView(PETSC_NULL, viewer);
#else
      return ::PetscOptionsView(viewer);
#endif
    }
 
    // PETSc routine
 
    inline PetscErrorCode PetscOptionsInsertString(const char in_str[])
    {
#if (PETSC_VERSION_MINOR >= 7)
      return ::PetscOptionsInsertString(PETSC_NULL, in_str);
#else
      return ::PetscOptionsInsertString(in_str);
#endif
    }
 
  } // end namespace PETSc
} // end namespace AMDiS