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: }