1 /*
2 Defines a preconditioner defined by R^T S R
3 */
4 #include <petsc/private/pcimpl.h>
5 #include <petscksp.h> /*I "petscksp.h" I*/
6
7 typedef struct {
8 KSP ksp;
9 Mat R, P;
10 Vec b, x;
11 PetscErrorCode (*computeasub)(PC, Mat, Mat, Mat *, void *);
12 void *computeasub_ctx;
13 } PC_Galerkin;
14
PCApply_Galerkin(PC pc,Vec x,Vec y)15 static PetscErrorCode PCApply_Galerkin(PC pc, Vec x, Vec y)
16 {
17 PC_Galerkin *jac = (PC_Galerkin *)pc->data;
18
19 PetscFunctionBegin;
20 if (jac->R) {
21 PetscCall(MatRestrict(jac->R, x, jac->b));
22 } else {
23 PetscCall(MatRestrict(jac->P, x, jac->b));
24 }
25 PetscCall(KSPSolve(jac->ksp, jac->b, jac->x));
26 PetscCall(KSPCheckSolve(jac->ksp, pc, jac->x));
27 if (jac->P) {
28 PetscCall(MatInterpolate(jac->P, jac->x, y));
29 } else {
30 PetscCall(MatInterpolate(jac->R, jac->x, y));
31 }
32 PetscFunctionReturn(PETSC_SUCCESS);
33 }
34
PCSetUp_Galerkin(PC pc)35 static PetscErrorCode PCSetUp_Galerkin(PC pc)
36 {
37 PC_Galerkin *jac = (PC_Galerkin *)pc->data;
38 PetscBool a;
39 Vec *xx, *yy;
40
41 PetscFunctionBegin;
42 if (jac->computeasub) {
43 Mat Ap;
44 if (!pc->setupcalled) {
45 PetscCall((*jac->computeasub)(pc, pc->pmat, NULL, &Ap, jac->computeasub_ctx));
46 PetscCall(KSPSetOperators(jac->ksp, Ap, Ap));
47 PetscCall(MatDestroy(&Ap));
48 } else {
49 PetscCall(KSPGetOperators(jac->ksp, NULL, &Ap));
50 PetscCall((*jac->computeasub)(pc, pc->pmat, Ap, NULL, jac->computeasub_ctx));
51 }
52 }
53
54 if (!jac->x) {
55 PetscCall(KSPGetOperatorsSet(jac->ksp, &a, NULL));
56 PetscCheck(a, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set operator of PCGALERKIN KSP with PCGalerkinGetKSP()/KSPSetOperators()");
57 PetscCall(KSPCreateVecs(jac->ksp, 1, &xx, 1, &yy));
58 jac->x = *xx;
59 jac->b = *yy;
60 PetscCall(PetscFree(xx));
61 PetscCall(PetscFree(yy));
62 }
63 PetscCheck(jac->R || jac->P, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set restriction or interpolation of PCGALERKIN with PCGalerkinSetRestriction()/Interpolation()");
64 /* should check here that sizes of R/P match size of a */
65 PetscFunctionReturn(PETSC_SUCCESS);
66 }
67
PCReset_Galerkin(PC pc)68 static PetscErrorCode PCReset_Galerkin(PC pc)
69 {
70 PC_Galerkin *jac = (PC_Galerkin *)pc->data;
71
72 PetscFunctionBegin;
73 PetscCall(MatDestroy(&jac->R));
74 PetscCall(MatDestroy(&jac->P));
75 PetscCall(VecDestroy(&jac->x));
76 PetscCall(VecDestroy(&jac->b));
77 PetscCall(KSPReset(jac->ksp));
78 PetscFunctionReturn(PETSC_SUCCESS);
79 }
80
PCDestroy_Galerkin(PC pc)81 static PetscErrorCode PCDestroy_Galerkin(PC pc)
82 {
83 PC_Galerkin *jac = (PC_Galerkin *)pc->data;
84
85 PetscFunctionBegin;
86 PetscCall(PCReset_Galerkin(pc));
87 PetscCall(KSPDestroy(&jac->ksp));
88 PetscCall(PetscFree(pc->data));
89 PetscFunctionReturn(PETSC_SUCCESS);
90 }
91
PCView_Galerkin(PC pc,PetscViewer viewer)92 static PetscErrorCode PCView_Galerkin(PC pc, PetscViewer viewer)
93 {
94 PC_Galerkin *jac = (PC_Galerkin *)pc->data;
95 PetscBool isascii;
96
97 PetscFunctionBegin;
98 PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
99 if (isascii) {
100 PetscCall(PetscViewerASCIIPrintf(viewer, " KSP on Galerkin follow\n"));
101 PetscCall(PetscViewerASCIIPrintf(viewer, " ---------------------------------\n"));
102 }
103 PetscCall(KSPView(jac->ksp, viewer));
104 PetscFunctionReturn(PETSC_SUCCESS);
105 }
106
PCGalerkinGetKSP_Galerkin(PC pc,KSP * ksp)107 static PetscErrorCode PCGalerkinGetKSP_Galerkin(PC pc, KSP *ksp)
108 {
109 PC_Galerkin *jac = (PC_Galerkin *)pc->data;
110
111 PetscFunctionBegin;
112 *ksp = jac->ksp;
113 PetscFunctionReturn(PETSC_SUCCESS);
114 }
115
PCGalerkinSetRestriction_Galerkin(PC pc,Mat R)116 static PetscErrorCode PCGalerkinSetRestriction_Galerkin(PC pc, Mat R)
117 {
118 PC_Galerkin *jac = (PC_Galerkin *)pc->data;
119
120 PetscFunctionBegin;
121 PetscCall(PetscObjectReference((PetscObject)R));
122 PetscCall(MatDestroy(&jac->R));
123 jac->R = R;
124 PetscFunctionReturn(PETSC_SUCCESS);
125 }
126
PCGalerkinSetInterpolation_Galerkin(PC pc,Mat P)127 static PetscErrorCode PCGalerkinSetInterpolation_Galerkin(PC pc, Mat P)
128 {
129 PC_Galerkin *jac = (PC_Galerkin *)pc->data;
130
131 PetscFunctionBegin;
132 PetscCall(PetscObjectReference((PetscObject)P));
133 PetscCall(MatDestroy(&jac->P));
134 jac->P = P;
135 PetscFunctionReturn(PETSC_SUCCESS);
136 }
137
PCGalerkinSetComputeSubmatrix_Galerkin(PC pc,PetscErrorCode (* computeAsub)(PC,Mat,Mat,Mat *,void *),PetscCtx ctx)138 static PetscErrorCode PCGalerkinSetComputeSubmatrix_Galerkin(PC pc, PetscErrorCode (*computeAsub)(PC, Mat, Mat, Mat *, void *), PetscCtx ctx)
139 {
140 PC_Galerkin *jac = (PC_Galerkin *)pc->data;
141
142 PetscFunctionBegin;
143 jac->computeasub = computeAsub;
144 jac->computeasub_ctx = ctx;
145 PetscFunctionReturn(PETSC_SUCCESS);
146 }
147
148 /*@
149 PCGalerkinSetRestriction - Sets the restriction operator for the `PCGALERKIN` preconditioner
150
151 Logically Collective
152
153 Input Parameters:
154 + pc - the preconditioner context
155 - R - the restriction operator
156
157 Level: intermediate
158
159 Note:
160 Either this or `PCGalerkinSetInterpolation()` or both must be called
161
162 .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
163 `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
164 @*/
PCGalerkinSetRestriction(PC pc,Mat R)165 PetscErrorCode PCGalerkinSetRestriction(PC pc, Mat R)
166 {
167 PetscFunctionBegin;
168 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
169 PetscTryMethod(pc, "PCGalerkinSetRestriction_C", (PC, Mat), (pc, R));
170 PetscFunctionReturn(PETSC_SUCCESS);
171 }
172
173 /*@
174 PCGalerkinSetInterpolation - Sets the interpolation operator for the `PCGALERKIN` preconditioner
175
176 Logically Collective
177
178 Input Parameters:
179 + pc - the preconditioner context
180 - P - the interpolation operator
181
182 Level: intermediate
183
184 Note:
185 Either this or `PCGalerkinSetRestriction()` or both must be called
186
187 .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
188 `PCGalerkinSetRestriction()`, `PCGalerkinGetKSP()`
189 @*/
PCGalerkinSetInterpolation(PC pc,Mat P)190 PetscErrorCode PCGalerkinSetInterpolation(PC pc, Mat P)
191 {
192 PetscFunctionBegin;
193 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
194 PetscTryMethod(pc, "PCGalerkinSetInterpolation_C", (PC, Mat), (pc, P));
195 PetscFunctionReturn(PETSC_SUCCESS);
196 }
197
198 /*@C
199 PCGalerkinSetComputeSubmatrix - Provide a routine that will be called to compute the Galerkin submatrix
200
201 Logically Collective
202
203 Input Parameters:
204 + pc - the preconditioner context
205 . computeAsub - routine that computes the submatrix from the global matrix
206 - ctx - context used by the routine, or `NULL`
207
208 Calling sequence of `computeAsub`:
209 + pc - the `PCGALERKIN` preconditioner
210 . A - the matrix in the `PCGALERKIN`
211 . Ap - the computed submatrix from any previous computation, if `NULL` it has not previously been computed
212 . cAp - the submatrix computed by this routine
213 - ctx - optional user-defined function context
214
215 Level: intermediate
216
217 Notes:
218 Instead of providing this routine you can call `PCGalerkinGetKSP()` and then `KSPSetOperators()` to provide the submatrix,
219 but that will not work for multiple `KSPSolve()`s with different matrices unless you call it for each solve.
220
221 This routine is called each time the outer matrix is changed. In the first call the Ap argument is `NULL` and the routine should create the
222 matrix and computes its values in cAp. On each subsequent call the routine should up the Ap matrix.
223
224 Developer Notes:
225 If the user does not call this routine nor call `PCGalerkinGetKSP()` and `KSPSetOperators()` then `PCGALERKIN`
226 could automatically compute the submatrix via calls to `MatGalerkin()` or `MatRARt()`
227
228 .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
229 `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
230 @*/
PCGalerkinSetComputeSubmatrix(PC pc,PetscErrorCode (* computeAsub)(PC pc,Mat A,Mat Ap,Mat * cAp,PetscCtx ctx),PetscCtx ctx)231 PetscErrorCode PCGalerkinSetComputeSubmatrix(PC pc, PetscErrorCode (*computeAsub)(PC pc, Mat A, Mat Ap, Mat *cAp, PetscCtx ctx), PetscCtx ctx)
232 {
233 PetscFunctionBegin;
234 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
235 PetscTryMethod(pc, "PCGalerkinSetComputeSubmatrix_C", (PC, PetscErrorCode (*)(PC, Mat, Mat, Mat *, void *), void *), (pc, computeAsub, ctx));
236 PetscFunctionReturn(PETSC_SUCCESS);
237 }
238
239 /*@
240 PCGalerkinGetKSP - Gets the `KSP` object in the `PCGALERKIN`
241
242 Not Collective
243
244 Input Parameter:
245 . pc - the preconditioner context
246
247 Output Parameter:
248 . ksp - the `KSP` object
249
250 Level: intermediate
251
252 Note:
253 Once you have called this routine you can call `KSPSetOperators()` on the resulting `KSP` to provide the operator for the Galerkin problem,
254 an alternative is to use `PCGalerkinSetComputeSubmatrix()` to provide a routine that computes the submatrix as needed.
255
256 .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
257 `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinSetComputeSubmatrix()`
258 @*/
PCGalerkinGetKSP(PC pc,KSP * ksp)259 PetscErrorCode PCGalerkinGetKSP(PC pc, KSP *ksp)
260 {
261 PetscFunctionBegin;
262 PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
263 PetscAssertPointer(ksp, 2);
264 PetscUseMethod(pc, "PCGalerkinGetKSP_C", (PC, KSP *), (pc, ksp));
265 PetscFunctionReturn(PETSC_SUCCESS);
266 }
267
PCSetFromOptions_Galerkin(PC pc,PetscOptionItems PetscOptionsObject)268 static PetscErrorCode PCSetFromOptions_Galerkin(PC pc, PetscOptionItems PetscOptionsObject)
269 {
270 PC_Galerkin *jac = (PC_Galerkin *)pc->data;
271 const char *prefix;
272 PetscBool flg;
273
274 PetscFunctionBegin;
275 PetscCall(KSPGetOptionsPrefix(jac->ksp, &prefix));
276 PetscCall(PetscStrendswith(prefix, "galerkin_", &flg));
277 if (!flg) {
278 PetscCall(PCGetOptionsPrefix(pc, &prefix));
279 PetscCall(KSPSetOptionsPrefix(jac->ksp, prefix));
280 PetscCall(KSPAppendOptionsPrefix(jac->ksp, "galerkin_"));
281 }
282
283 PetscOptionsHeadBegin(PetscOptionsObject, "Galerkin options");
284 if (jac->ksp) PetscCall(KSPSetFromOptions(jac->ksp));
285 PetscOptionsHeadEnd();
286 PetscFunctionReturn(PETSC_SUCCESS);
287 }
288
289 /*MC
290 PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T)
291
292 Level: intermediate
293
294 Note:
295 Use
296 .vb
297 `PCGalerkinSetRestriction`(pc,R) and/or `PCGalerkinSetInterpolation`(pc,P)
298 `PCGalerkinGetKSP`(pc,&ksp);
299 `KSPSetOperators`(ksp,A,....)
300 ...
301 .ve
302
303 Developer Notes:
304 If `KSPSetOperators()` has not been called on the inner `KSP` then `PCGALERKIN` could use `MatRARt()` or `MatPtAP()` to compute
305 the operators automatically.
306
307 Should there be a prefix for the inner `KSP`?
308
309 There is no `KSPSetFromOptions_Galerkin()` that calls `KSPSetFromOptions()` on the inner `KSP`
310
311 .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
312 `PCSHELL`, `PCKSP`, `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
313 M*/
314
PCCreate_Galerkin(PC pc)315 PETSC_EXTERN PetscErrorCode PCCreate_Galerkin(PC pc)
316 {
317 PC_Galerkin *jac;
318
319 PetscFunctionBegin;
320 PetscCall(PetscNew(&jac));
321
322 pc->ops->apply = PCApply_Galerkin;
323 pc->ops->setup = PCSetUp_Galerkin;
324 pc->ops->reset = PCReset_Galerkin;
325 pc->ops->destroy = PCDestroy_Galerkin;
326 pc->ops->view = PCView_Galerkin;
327 pc->ops->setfromoptions = PCSetFromOptions_Galerkin;
328 pc->ops->applyrichardson = NULL;
329
330 PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->ksp));
331 PetscCall(KSPSetNestLevel(jac->ksp, pc->kspnestlevel));
332 PetscCall(KSPSetErrorIfNotConverged(jac->ksp, pc->erroriffailure));
333 PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->ksp, (PetscObject)pc, 1));
334
335 pc->data = (void *)jac;
336
337 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetRestriction_C", PCGalerkinSetRestriction_Galerkin));
338 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetInterpolation_C", PCGalerkinSetInterpolation_Galerkin));
339 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinGetKSP_C", PCGalerkinGetKSP_Galerkin));
340 PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetComputeSubmatrix_C", PCGalerkinSetComputeSubmatrix_Galerkin));
341 PetscFunctionReturn(PETSC_SUCCESS);
342 }
343