Actual source code: pcmat.c

  1: #include <petsc/private/pcimpl.h>

  3: typedef enum {
  4:   PCMATOP_UNSPECIFIED,
  5:   PCMATOP_MULT,
  6:   PCMATOP_MULT_TRANSPOSE,
  7:   PCMATOP_MULT_HERMITIAN_TRANSPOSE,
  8:   PCMATOP_SOLVE,
  9:   PCMATOP_SOLVE_TRANSPOSE,
 10: } PCMatOperation;

 12: const char *const PCMatOpTypes[] = {"Unspecified", "Mult", "MultTranspose", "MultHermitianTranspose", "Solve", "SolveTranspose", NULL};

 14: typedef struct _PCMAT {
 15:   PCMatOperation apply;
 16: } PC_Mat;

 18: static PetscErrorCode PCApply_Mat(PC pc, Vec x, Vec y)
 19: {
 20:   PC_Mat *pcmat = (PC_Mat *)pc->data;

 22:   PetscFunctionBegin;
 23:   switch (pcmat->apply) {
 24:   case PCMATOP_MULT:
 25:     PetscCall(MatMult(pc->pmat, x, y));
 26:     break;
 27:   case PCMATOP_MULT_TRANSPOSE:
 28:     PetscCall(MatMultTranspose(pc->pmat, x, y));
 29:     break;
 30:   case PCMATOP_SOLVE:
 31:     PetscCall(MatSolve(pc->pmat, x, y));
 32:     break;
 33:   case PCMATOP_SOLVE_TRANSPOSE:
 34:     PetscCall(MatSolveTranspose(pc->pmat, x, y));
 35:     break;
 36:   case PCMATOP_MULT_HERMITIAN_TRANSPOSE:
 37:     PetscCall(MatMultHermitianTranspose(pc->pmat, x, y));
 38:     break;
 39:   default:
 40:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Unsupported %s case", PCMatOpTypes[pcmat->apply]);
 41:   }
 42:   PetscFunctionReturn(PETSC_SUCCESS);
 43: }

 45: static PetscErrorCode PCSetUp_Mat(PC pc)
 46: {
 47:   PC_Mat *pcmat = (PC_Mat *)pc->data;

 49:   PetscFunctionBegin;
 50:   if (pcmat->apply == PCMATOP_UNSPECIFIED) {
 51:     PetscBool hassolve;
 52:     PetscCall(MatHasOperation(pc->pmat, MATOP_SOLVE, &hassolve));
 53:     if (hassolve) pcmat->apply = PCMATOP_SOLVE;
 54:     else pcmat->apply = PCMATOP_MULT;
 55:   }
 56:   PetscFunctionReturn(PETSC_SUCCESS);
 57: }

 59: static PetscErrorCode PCMatApply_Mat(PC pc, Mat X, Mat Y)
 60: {
 61:   PC_Mat *pcmat = (PC_Mat *)pc->data;
 62:   Mat     W;

 64:   PetscFunctionBegin;
 65:   switch (pcmat->apply) {
 66:   case PCMATOP_MULT:
 67:     PetscCall(MatMatMult(pc->pmat, X, MAT_REUSE_MATRIX, PETSC_CURRENT, &Y));
 68:     break;
 69:   case PCMATOP_MULT_TRANSPOSE:
 70:     PetscCall(MatTransposeMatMult(pc->pmat, X, MAT_REUSE_MATRIX, PETSC_CURRENT, &Y));
 71:     break;
 72:   case PCMATOP_SOLVE:
 73:     PetscCall(MatMatSolve(pc->pmat, X, Y));
 74:     break;
 75:   case PCMATOP_SOLVE_TRANSPOSE:
 76:     PetscCall(MatMatSolveTranspose(pc->pmat, X, Y));
 77:     break;
 78:   case PCMATOP_MULT_HERMITIAN_TRANSPOSE:
 79:     PetscCall(MatDuplicate(X, MAT_COPY_VALUES, &W));
 80:     PetscCall(MatConjugate(W));
 81:     PetscCall(MatTransposeMatMult(pc->pmat, W, MAT_REUSE_MATRIX, PETSC_CURRENT, &Y));
 82:     PetscCall(MatConjugate(Y));
 83:     PetscCall(MatDestroy(&W));
 84:     break;
 85:   default:
 86:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Unsupported %s case", PCMatOpTypes[pcmat->apply]);
 87:   }
 88:   PetscFunctionReturn(PETSC_SUCCESS);
 89: }

 91: static PetscErrorCode PCApplyTranspose_Mat(PC pc, Vec x, Vec y)
 92: {
 93:   PC_Mat *pcmat = (PC_Mat *)pc->data;
 94:   Vec     w;

 96:   PetscFunctionBegin;
 97:   switch (pcmat->apply) {
 98:   case PCMATOP_MULT:
 99:     PetscCall(MatMultTranspose(pc->pmat, x, y));
100:     break;
101:   case PCMATOP_MULT_TRANSPOSE:
102:     PetscCall(MatMult(pc->pmat, x, y));
103:     break;
104:   case PCMATOP_SOLVE:
105:     PetscCall(MatSolveTranspose(pc->pmat, x, y));
106:     break;
107:   case PCMATOP_SOLVE_TRANSPOSE:
108:     PetscCall(MatSolve(pc->pmat, x, y));
109:     break;
110:   case PCMATOP_MULT_HERMITIAN_TRANSPOSE:
111:     PetscCall(VecDuplicate(x, &w));
112:     PetscCall(VecCopy(x, w));
113:     PetscCall(VecConjugate(w));
114:     PetscCall(MatMult(pc->pmat, w, y));
115:     PetscCall(VecConjugate(y));
116:     PetscCall(VecDestroy(&w));
117:     break;
118:   default:
119:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Unsupported %s case", PCMatOpTypes[pcmat->apply]);
120:   }
121:   PetscFunctionReturn(PETSC_SUCCESS);
122: }

124: static PetscErrorCode PCDestroy_Mat(PC pc)
125: {
126:   PetscFunctionBegin;
127:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatSetApplyOperation_C", NULL));
128:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatGetApplyOperation_C", NULL));
129:   PetscCall(PetscFree(pc->data));
130:   PetscFunctionReturn(PETSC_SUCCESS);
131: }

133: /*@
134:   PCMatSetApplyOperation - Set which matrix operation of the preconditioning matrix implements `PCApply()` for `PCMAT`.

136:   Logically collective

138:   Input Parameters:
139: + pc    - An instance of `PCMAT`
140: - matop - The selected `MatOperation`

142:   Level: intermediate

144:   Note:
145:   If you have a matrix type that implements an exact inverse that isn't a factorization,
146:   you can use `PCMatSetApplyOperation(pc, MATOP_SOLVE)`.

148: .seealso: [](ch_ksp), `PCMAT`, `PCMatGetApplyOperation()`, `PCApply()`, `MatOperation`
149: @*/
150: PetscErrorCode PCMatSetApplyOperation(PC pc, MatOperation matop)
151: {
152:   PetscFunctionBegin;
154:   PetscTryMethod(pc, "PCMatSetApplyOperation_C", (PC, MatOperation), (pc, matop));
155:   PetscFunctionReturn(PETSC_SUCCESS);
156: }

158: /*@
159:   PCMatGetApplyOperation - Get which matrix operation of the preconditioning matrix implements `PCApply()` for `PCMAT`.

161:   Logically collective

163:   Input Parameter:
164: . pc - An instance of `PCMAT`

166:   Output Parameter:
167: . matop - The `MatOperation`

169:   Level: intermediate

171: .seealso: [](ch_ksp), `PCMAT`, `PCMatSetApplyOperation()`, `PCApply()`, `MatOperation`
172: @*/
173: PetscErrorCode PCMatGetApplyOperation(PC pc, MatOperation *matop)
174: {
175:   PetscFunctionBegin;
177:   PetscAssertPointer(matop, 2);
178:   PetscUseMethod(pc, "PCMatGetApplyOperation_C", (PC, MatOperation *), (pc, matop));
179:   PetscFunctionReturn(PETSC_SUCCESS);
180: }

182: static PetscErrorCode PCMatSetApplyOperation_Mat(PC pc, MatOperation matop)
183: {
184:   PC_Mat        *pcmat = (PC_Mat *)pc->data;
185:   PCMatOperation pcmatop;

187:   PetscFunctionBegin;
188:   // clang-format off
189: #define MATOP_TO_PCMATOP_CASE(var, OP) case MATOP_##OP: (var) = PCMATOP_##OP; break
190:   switch (matop) {
191:   MATOP_TO_PCMATOP_CASE(pcmatop, MULT);
192:   MATOP_TO_PCMATOP_CASE(pcmatop, MULT_TRANSPOSE);
193:   MATOP_TO_PCMATOP_CASE(pcmatop, MULT_HERMITIAN_TRANSPOSE);
194:   MATOP_TO_PCMATOP_CASE(pcmatop, SOLVE);
195:   MATOP_TO_PCMATOP_CASE(pcmatop, SOLVE_TRANSPOSE);
196:   default:
197:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Unsupported MatOperation %d for PCMatSetApplyOperation()", (int)matop);
198:   }
199: #undef MATOP_TO_PCMATOP_CASE
200:   // clang-format on

202:   pcmat->apply = pcmatop;
203:   PetscFunctionReturn(PETSC_SUCCESS);
204: }

206: static PetscErrorCode PCMatGetApplyOperation_Mat(PC pc, MatOperation *matop_p)
207: {
208:   PC_Mat      *pcmat = (PC_Mat *)pc->data;
209:   MatOperation matop = MATOP_MULT;

211:   PetscFunctionBegin;
212:   if (!pc->setupcalled) PetscCall(PCSetUp(pc));

214:   // clang-format off
215: #define PCMATOP_TO_MATOP_CASE(var, OP) case PCMATOP_##OP: (var) = MATOP_##OP; break
216:   switch (pcmat->apply) {
217:   PCMATOP_TO_MATOP_CASE(matop, MULT);
218:   PCMATOP_TO_MATOP_CASE(matop, MULT_TRANSPOSE);
219:   PCMATOP_TO_MATOP_CASE(matop, MULT_HERMITIAN_TRANSPOSE);
220:   PCMATOP_TO_MATOP_CASE(matop, SOLVE);
221:   PCMATOP_TO_MATOP_CASE(matop, SOLVE_TRANSPOSE);
222:   default:
223:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Unsupported %s case", PCMatOpTypes[pcmat->apply]);
224:   }
225: #undef PCMATOP_TO_MATOP_CASE
226:   // clang-format on

228:   *matop_p = matop;
229:   PetscFunctionReturn(PETSC_SUCCESS);
230: }

232: static PetscErrorCode PCView_Mat(PC pc, PetscViewer viewer)
233: {
234:   PC_Mat   *pcmat = (PC_Mat *)pc->data;
235:   PetscBool iascii;

237:   PetscFunctionBegin;
238:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
239:   if (iascii) { PetscCall(PetscViewerASCIIPrintf(viewer, "PCApply() == Mat%s()\n", PCMatOpTypes[pcmat->apply])); }
240:   PetscFunctionReturn(PETSC_SUCCESS);
241: }

243: /*MC
244:      PCMAT - A preconditioner obtained by applying an operation of the `pmat` provided in
245:              in `PCSetOperators(pc,amat,pmat)` or `KSPSetOperators(ksp,amat,pmat)`.  By default, the operation is `MATOP_MULT`,
246:              meaning that the `pmat` provides an approximate inverse of `amat`.
247:              If some other operation of `pmat` implements the approximate inverse,
248:              use `PCMatSetApplyOperation()` to select that operation.

250:    Level: intermediate

252: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSHELL`, `MatOperation`, `PCMatSetApplyOperation()`, `PCMatGetApplyOperation()`
253: M*/

255: PETSC_EXTERN PetscErrorCode PCCreate_Mat(PC pc)
256: {
257:   PC_Mat *data;

259:   PetscFunctionBegin;
260:   PetscCall(PetscNew(&data));
261:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatSetApplyOperation_C", PCMatSetApplyOperation_Mat));
262:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatGetApplyOperation_C", PCMatGetApplyOperation_Mat));
263:   pc->data                     = data;
264:   pc->ops->apply               = PCApply_Mat;
265:   pc->ops->matapply            = PCMatApply_Mat;
266:   pc->ops->applytranspose      = PCApplyTranspose_Mat;
267:   pc->ops->setup               = PCSetUp_Mat;
268:   pc->ops->destroy             = PCDestroy_Mat;
269:   pc->ops->setfromoptions      = NULL;
270:   pc->ops->view                = PCView_Mat;
271:   pc->ops->applyrichardson     = NULL;
272:   pc->ops->applysymmetricleft  = NULL;
273:   pc->ops->applysymmetricright = NULL;
274:   PetscFunctionReturn(PETSC_SUCCESS);
275: }