Interface.hpp 5.34 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
#pragma once

#include <petscksp.h>
#include <petsc/private/kspimpl.h>

#include <amdis/Output.hpp>

namespace AMDiS
{
  namespace PETSc
  {
    // KSP routines

    /// Print the residual norm at every 10th iteration.
    inline PetscErrorCode KSPMonitorCyclic(KSP ksp, PetscInt n, PetscReal rnorm, void* ctx)
    {
      if (n % 10 == 0)
        msg("iteration {:>4}: resid {:<.6e}", n, rnorm);
      return 0;
    }

    /// Prints the residual norm, the true residual norm, and the relative residual norm at each iteration.
    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;
    }

    /// Prints the true residual norm as well as the preconditioned residual norm at each iteration of an iterative solver.
    inline PetscErrorCode KSPMonitorTrueResidualNorm(KSP ksp, PetscInt n, PetscReal rnorm, void* ctx)
    {
      assert(ctx != nullptr);
#if (PETSC_VERSION_MINOR >= 7)
      return ::KSPMonitorTrueResidualNorm(ksp, n, rnorm, (PetscViewerAndFormat*)(ctx));
#else
      return ::KSPMonitorTrueResidualNorm(ksp, n, rnorm, ctx);
#endif
    }

    /// Print the residual norm at each iteration of an iterative solver.
    inline PetscErrorCode KSPMonitorDefault(KSP ksp, PetscInt n, PetscReal rnorm, void* ctx)
    {
      assert(ctx != nullptr);
#if (PETSC_VERSION_MINOR >= 7)
      return ::KSPMonitorDefault(ksp, n, rnorm, (PetscViewerAndFormat*)(ctx));
#else
      return ::KSPMonitorDefault(ksp, n, rnorm, ctx);
#endif
    }

    /// Sets an additional function to be called at every iteration to monitor the residual/error etc.
    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
    }

    /// Gets the matrix associated with the linear system and a (possibly) different one associated with the preconditioner.
    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
    }

    /// Sets the matrix associated with the linear system and a (possibly) different one associated with the preconditioner.
    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
    }

    // Mat routines

    /// Get vector(s) compatible with the matrix, i.e. with the same parallel layout
    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
    }

    /// Gets a single submatrix on the same number of processors as the original matrix.
    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
    }

    /// Removes all the components of a null space from a vector.
    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)
    /// Sets the software that is used to perform the factorization
    inline PetscErrorCode PCFactorSetMatSolverType(PC pc, MatSolverType stype)
    {
      return ::PCFactorSetMatSolverType(pc, stype);
    }
#else
    /// Sets the software that is used to perform the factorization
    template <class MatSolverType>
    inline PetscErrorCode PCFactorSetMatSolverType(PC pc, MatSolverType stype)
    {
      return ::PCFactorSetMatSolverPackage(pc, stype);
    }
#endif

    /// Prints the options that have been loaded. This is useful for debugging purposes.
    inline PetscErrorCode PetscOptionsView(PetscViewer viewer)
    {
#if (PETSC_VERSION_MINOR >= 7)
      return ::PetscOptionsView(PETSC_NULL, viewer);
#else
      return ::PetscOptionsView(viewer);
#endif
    }

    // PETSc routine

    /// Inserts options into the database from a string
    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