Actual source code: pcmat.c
1: #include <petsc/private/pcimpl.h>
3: typedef enum {
4: PCMATOP_MULT,
5: PCMATOP_MULT_TRANSPOSE,
6: PCMATOP_MULT_HERMITIAN_TRANSPOSE,
7: PCMATOP_SOLVE,
8: PCMATOP_SOLVE_TRANSPOSE,
9: } PCMatOperation;
11: const char *const PCMatOpTypes[] = {"Mult", "MultTranspose", "MultHermitianTranspose", "Solve", "SolveTranspose", NULL};
13: typedef struct _PCMAT {
14: PCMatOperation apply;
15: } PC_Mat;
17: static PetscErrorCode PCApply_Mat(PC pc, Vec x, Vec y)
18: {
19: PC_Mat *pcmat = (PC_Mat *)pc->data;
21: PetscFunctionBegin;
22: switch (pcmat->apply) {
23: case PCMATOP_MULT:
24: PetscCall(MatMult(pc->pmat, x, y));
25: break;
26: case PCMATOP_MULT_TRANSPOSE:
27: PetscCall(MatMultTranspose(pc->pmat, x, y));
28: break;
29: case PCMATOP_SOLVE:
30: PetscCall(MatSolve(pc->pmat, x, y));
31: break;
32: case PCMATOP_SOLVE_TRANSPOSE:
33: PetscCall(MatSolveTranspose(pc->pmat, x, y));
34: break;
35: case PCMATOP_MULT_HERMITIAN_TRANSPOSE:
36: PetscCall(MatMultHermitianTranspose(pc->pmat, x, y));
37: break;
38: }
39: PetscFunctionReturn(PETSC_SUCCESS);
40: }
42: static PetscErrorCode PCMatApply_Mat(PC pc, Mat X, Mat Y)
43: {
44: PC_Mat *pcmat = (PC_Mat *)pc->data;
46: PetscFunctionBegin;
47: switch (pcmat->apply) {
48: case PCMATOP_MULT:
49: PetscCall(MatMatMult(pc->pmat, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &Y));
50: break;
51: case PCMATOP_MULT_TRANSPOSE:
52: PetscCall(MatTransposeMatMult(pc->pmat, X, MAT_REUSE_MATRIX, PETSC_DEFAULT, &Y));
53: break;
54: case PCMATOP_SOLVE:
55: PetscCall(MatMatSolve(pc->pmat, X, Y));
56: break;
57: case PCMATOP_SOLVE_TRANSPOSE:
58: PetscCall(MatMatSolveTranspose(pc->pmat, X, Y));
59: break;
60: case PCMATOP_MULT_HERMITIAN_TRANSPOSE: {
61: Mat W;
63: PetscCall(MatDuplicate(X, MAT_COPY_VALUES, &W));
64: PetscCall(MatConjugate(W));
65: PetscCall(MatTransposeMatMult(pc->pmat, W, MAT_REUSE_MATRIX, PETSC_DEFAULT, &Y));
66: PetscCall(MatConjugate(Y));
67: PetscCall(MatDestroy(&W));
68: break;
69: }
70: }
71: PetscFunctionReturn(PETSC_SUCCESS);
72: }
74: static PetscErrorCode PCApplyTranspose_Mat(PC pc, Vec x, Vec y)
75: {
76: PC_Mat *pcmat = (PC_Mat *)pc->data;
78: PetscFunctionBegin;
79: switch (pcmat->apply) {
80: case PCMATOP_MULT:
81: PetscCall(MatMultTranspose(pc->pmat, x, y));
82: break;
83: case PCMATOP_MULT_TRANSPOSE:
84: PetscCall(MatMult(pc->pmat, x, y));
85: break;
86: case PCMATOP_SOLVE:
87: PetscCall(MatSolveTranspose(pc->pmat, x, y));
88: break;
89: case PCMATOP_SOLVE_TRANSPOSE:
90: PetscCall(MatSolve(pc->pmat, x, y));
91: break;
92: case PCMATOP_MULT_HERMITIAN_TRANSPOSE: {
93: Vec w;
95: PetscCall(VecDuplicate(x, &w));
96: PetscCall(VecCopy(x, w));
97: PetscCall(VecConjugate(w));
98: PetscCall(MatMult(pc->pmat, w, y));
99: PetscCall(VecConjugate(y));
100: PetscCall(VecDestroy(&w));
101: break;
102: }
103: }
104: PetscFunctionReturn(PETSC_SUCCESS);
105: }
107: static PetscErrorCode PCDestroy_Mat(PC pc)
108: {
109: PetscFunctionBegin;
110: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatSetApplyOperation_C", NULL));
111: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatGetApplyOperation_C", NULL));
112: PetscCall(PetscFree(pc->data));
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: /*@
117: PCMatSetApplyOperation - Set which matrix operation of the preconditioning matrix implements `PCApply()` for `PCMAT`.
119: Logically collective
121: Input Parameters:
122: + pc - An instance of `PCMAT`
123: - matop - The selected `MatOperation`
125: Level: intermediate
127: Note:
128: If you have a matrix type that implements an exact inverse that isn't a factorization,
129: you can use `PCMatSetApplyOperation(pc, MATOP_SOLVE)`.
131: .seealso: [](ch_ksp), `PCMAT`, `PCMatGetApplyOperation()`, `PCApply()`, `MatOperation`
132: @*/
133: PetscErrorCode PCMatSetApplyOperation(PC pc, MatOperation matop)
134: {
135: PetscFunctionBegin;
137: PetscTryMethod((PetscObject)pc, "PCMatSetApplyOperation_C", (PC, MatOperation), (pc, matop));
138: PetscFunctionReturn(PETSC_SUCCESS);
139: }
141: /*@
142: PCMatGetApplyOperation - Get which matrix operation of the preconditioning matrix implements `PCApply()` for `PCMAT`.
144: Logically collective
146: Input Parameter:
147: . pc - An instance of `PCMAT`
149: Output Parameter:
150: . matop - The `MatOperation`
152: Level: intermediate
154: .seealso: [](ch_ksp), `PCMAT`, `PCMatSetApplyOperation()`, `PCApply()`, `MatOperation`
155: @*/
156: PetscErrorCode PCMatGetApplyOperation(PC pc, MatOperation *matop)
157: {
158: PetscFunctionBegin;
160: PetscAssertPointer(matop, 2);
161: PetscUseMethod((PetscObject)pc, "PCMatGetApplyOperation_C", (PC, MatOperation *), (pc, matop));
162: PetscFunctionReturn(PETSC_SUCCESS);
163: }
165: static PetscErrorCode PCMatSetApplyOperation_Mat(PC pc, MatOperation matop)
166: {
167: PC_Mat *pcmat = (PC_Mat *)pc->data;
168: PCMatOperation pcmatop;
170: PetscFunctionBegin;
172: // clang-format off
173: #define MATOP_TO_PCMATOP_CASE(var, OP) case MATOP_##OP: (var) = PCMATOP_##OP; break
174: switch (matop) {
175: MATOP_TO_PCMATOP_CASE(pcmatop, MULT);
176: MATOP_TO_PCMATOP_CASE(pcmatop, MULT_TRANSPOSE);
177: MATOP_TO_PCMATOP_CASE(pcmatop, MULT_HERMITIAN_TRANSPOSE);
178: MATOP_TO_PCMATOP_CASE(pcmatop, SOLVE);
179: MATOP_TO_PCMATOP_CASE(pcmatop, SOLVE_TRANSPOSE);
180: default:
181: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Unsupported MatOperation %d for PCMatSetApplyOperation()", (int)matop);
182: }
183: #undef MATOP_TO_PCMATOP_CASE
184: // clang-format on
186: pcmat->apply = pcmatop;
187: PetscFunctionReturn(PETSC_SUCCESS);
188: }
190: static PetscErrorCode PCMatGetApplyOperation_Mat(PC pc, MatOperation *matop_p)
191: {
192: PC_Mat *pcmat = (PC_Mat *)pc->data;
193: MatOperation matop = MATOP_MULT;
195: PetscFunctionBegin;
197: // clang-format off
198: #define PCMATOP_TO_MATOP_CASE(var, OP) case PCMATOP_##OP: (var) = MATOP_##OP; break
199: switch (pcmat->apply) {
200: PCMATOP_TO_MATOP_CASE(matop, MULT);
201: PCMATOP_TO_MATOP_CASE(matop, MULT_TRANSPOSE);
202: PCMATOP_TO_MATOP_CASE(matop, MULT_HERMITIAN_TRANSPOSE);
203: PCMATOP_TO_MATOP_CASE(matop, SOLVE);
204: PCMATOP_TO_MATOP_CASE(matop, SOLVE_TRANSPOSE);
205: }
206: #undef PCMATOP_TO_MATOP_CASE
207: // clang-format on
209: *matop_p = matop;
210: PetscFunctionReturn(PETSC_SUCCESS);
211: }
213: static PetscErrorCode PCView_Mat(PC pc, PetscViewer viewer)
214: {
215: PC_Mat *pcmat = (PC_Mat *)pc->data;
216: PetscBool iascii;
218: PetscFunctionBegin;
219: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
220: if (iascii) { PetscCall(PetscViewerASCIIPrintf(viewer, "PCApply() == Mat%s()\n", PCMatOpTypes[pcmat->apply])); }
221: PetscFunctionReturn(PETSC_SUCCESS);
222: }
224: /*MC
225: PCMAT - A preconditioner obtained by applying an operation of the `pmat` provided in
226: in `PCSetOperators(pc,amat,pmat)` or `KSPSetOperators(ksp,amat,pmat)`. By default, the operation is `MATOP_MULT`,
227: meaning that the `pmat` provides an approximate inverse of `amat`.
228: If some other operation of `pmat` implements the approximate inverse,
229: use `PCMatSetApplyOperation()` to select that operation.
231: Level: intermediate
233: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSHELL`, `MatOperation`, `PCMatSetApplyOperation()`, `PCMatGetApplyOperation()`
234: M*/
236: PETSC_EXTERN PetscErrorCode PCCreate_Mat(PC pc)
237: {
238: PC_Mat *data;
240: PetscFunctionBegin;
241: PetscCall(PetscNew(&data));
242: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatSetApplyOperation_C", PCMatSetApplyOperation_Mat));
243: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMatGetApplyOperation_C", PCMatGetApplyOperation_Mat));
244: data->apply = PCMATOP_MULT;
245: pc->data = data;
246: pc->ops->apply = PCApply_Mat;
247: pc->ops->matapply = PCMatApply_Mat;
248: pc->ops->applytranspose = PCApplyTranspose_Mat;
249: pc->ops->setup = NULL;
250: pc->ops->destroy = PCDestroy_Mat;
251: pc->ops->setfromoptions = NULL;
252: pc->ops->view = PCView_Mat;
253: pc->ops->applyrichardson = NULL;
254: pc->ops->applysymmetricleft = NULL;
255: pc->ops->applysymmetricright = NULL;
256: PetscFunctionReturn(PETSC_SUCCESS);
257: }