1 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
2
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;
11
12 const char *const PCMatOpTypes[] = {"Unspecified", "Mult", "MultTranspose", "MultHermitianTranspose", "Solve", "SolveTranspose", NULL};
13
14 typedef struct _PCMAT {
15 PCMatOperation apply;
16 } PC_Mat;
17
PCApply_Mat(PC pc,Vec x,Vec y)18 static PetscErrorCode PCApply_Mat(PC pc, Vec x, Vec y)
19 {
20 PC_Mat *pcmat = (PC_Mat *)pc->data;
21
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 }
44
PCSetUp_Mat(PC pc)45 static PetscErrorCode PCSetUp_Mat(PC pc)
46 {
47 PC_Mat *pcmat = (PC_Mat *)pc->data;
48
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 }
58
PCMatApply_Mat(PC pc,Mat X,Mat Y)59 static PetscErrorCode PCMatApply_Mat(PC pc, Mat X, Mat Y)
60 {
61 PC_Mat *pcmat = (PC_Mat *)pc->data;
62 Mat W;
63
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 }
90
PCApplyTranspose_Mat(PC pc,Vec x,Vec y)91 static PetscErrorCode PCApplyTranspose_Mat(PC pc, Vec x, Vec y)
92 {
93 PC_Mat *pcmat = (PC_Mat *)pc->data;
94 Vec w;
95
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 }
123
PCDestroy_Mat(PC pc)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 }
132
133 /*@
134 PCMatSetApplyOperation - Set which matrix operation of the matrix implements `PCApply()` for `PCMAT`.
135
136 Logically collective
137
138 Input Parameters:
139 + pc - An instance of `PCMAT`
140 - matop - The selected `MatOperation`
141
142 Level: intermediate
143
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)`.
147
148 .seealso: [](ch_ksp), `PCMAT`, `PCMatGetApplyOperation()`, `PCApply()`, `MatOperation`
149 @*/
PCMatSetApplyOperation(PC pc,MatOperation matop)150 PetscErrorCode PCMatSetApplyOperation(PC pc, MatOperation matop)
151 {
152 PetscFunctionBegin;
153 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
154 PetscTryMethod(pc, "PCMatSetApplyOperation_C", (PC, MatOperation), (pc, matop));
155 PetscFunctionReturn(PETSC_SUCCESS);
156 }
157
158 /*@
159 PCMatGetApplyOperation - Get which matrix operation of the matrix implements `PCApply()` for `PCMAT`.
160
161 Logically collective
162
163 Input Parameter:
164 . pc - An instance of `PCMAT`
165
166 Output Parameter:
167 . matop - The `MatOperation`
168
169 Level: intermediate
170
171 .seealso: [](ch_ksp), `PCMAT`, `PCMatSetApplyOperation()`, `PCApply()`, `MatOperation`
172 @*/
PCMatGetApplyOperation(PC pc,MatOperation * matop)173 PetscErrorCode PCMatGetApplyOperation(PC pc, MatOperation *matop)
174 {
175 PetscFunctionBegin;
176 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
177 PetscAssertPointer(matop, 2);
178 PetscUseMethod(pc, "PCMatGetApplyOperation_C", (PC, MatOperation *), (pc, matop));
179 PetscFunctionReturn(PETSC_SUCCESS);
180 }
181
PCMatSetApplyOperation_Mat(PC pc,MatOperation matop)182 static PetscErrorCode PCMatSetApplyOperation_Mat(PC pc, MatOperation matop)
183 {
184 PC_Mat *pcmat = (PC_Mat *)pc->data;
185 PCMatOperation pcmatop;
186
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
201
202 pcmat->apply = pcmatop;
203 PetscFunctionReturn(PETSC_SUCCESS);
204 }
205
PCMatGetApplyOperation_Mat(PC pc,MatOperation * matop_p)206 static PetscErrorCode PCMatGetApplyOperation_Mat(PC pc, MatOperation *matop_p)
207 {
208 PC_Mat *pcmat = (PC_Mat *)pc->data;
209 MatOperation matop = MATOP_MULT;
210
211 PetscFunctionBegin;
212 if (!pc->setupcalled) PetscCall(PCSetUp(pc));
213
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
227
228 *matop_p = matop;
229 PetscFunctionReturn(PETSC_SUCCESS);
230 }
231
PCView_Mat(PC pc,PetscViewer viewer)232 static PetscErrorCode PCView_Mat(PC pc, PetscViewer viewer)
233 {
234 PC_Mat *pcmat = (PC_Mat *)pc->data;
235 PetscBool isascii;
236
237 PetscFunctionBegin;
238 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
239 if (isascii) PetscCall(PetscViewerASCIIPrintf(viewer, "PCApply() == Mat%s()\n", PCMatOpTypes[pcmat->apply]));
240 PetscFunctionReturn(PETSC_SUCCESS);
241 }
242
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.
249
250 Level: intermediate
251
252 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCSHELL`, `MatOperation`, `PCMatSetApplyOperation()`, `PCMatGetApplyOperation()`
253 M*/
254
PCCreate_Mat(PC pc)255 PETSC_EXTERN PetscErrorCode PCCreate_Mat(PC pc)
256 {
257 PC_Mat *data;
258
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 }
276