Actual source code: galerkin.c
1: /*
2: Defines a preconditioner defined by R^T S R
3: */
4: #include <petsc/private/pcimpl.h>
5: #include <petscksp.h>
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;
15: static PetscErrorCode PCApply_Galerkin(PC pc, Vec x, Vec y)
16: {
17: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
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: }
35: static PetscErrorCode PCSetUp_Galerkin(PC pc)
36: {
37: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
38: PetscBool a;
39: Vec *xx, *yy;
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: }
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: }
68: static PetscErrorCode PCReset_Galerkin(PC pc)
69: {
70: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
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: }
81: static PetscErrorCode PCDestroy_Galerkin(PC pc)
82: {
83: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
85: PetscFunctionBegin;
86: PetscCall(PCReset_Galerkin(pc));
87: PetscCall(KSPDestroy(&jac->ksp));
88: PetscCall(PetscFree(pc->data));
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: static PetscErrorCode PCView_Galerkin(PC pc, PetscViewer viewer)
93: {
94: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
95: PetscBool iascii;
97: PetscFunctionBegin;
98: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
99: if (iascii) {
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: }
107: static PetscErrorCode PCGalerkinGetKSP_Galerkin(PC pc, KSP *ksp)
108: {
109: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
111: PetscFunctionBegin;
112: *ksp = jac->ksp;
113: PetscFunctionReturn(PETSC_SUCCESS);
114: }
116: static PetscErrorCode PCGalerkinSetRestriction_Galerkin(PC pc, Mat R)
117: {
118: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
120: PetscFunctionBegin;
121: PetscCall(PetscObjectReference((PetscObject)R));
122: PetscCall(MatDestroy(&jac->R));
123: jac->R = R;
124: PetscFunctionReturn(PETSC_SUCCESS);
125: }
127: static PetscErrorCode PCGalerkinSetInterpolation_Galerkin(PC pc, Mat P)
128: {
129: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
131: PetscFunctionBegin;
132: PetscCall(PetscObjectReference((PetscObject)P));
133: PetscCall(MatDestroy(&jac->P));
134: jac->P = P;
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }
138: static PetscErrorCode PCGalerkinSetComputeSubmatrix_Galerkin(PC pc, PetscErrorCode (*computeAsub)(PC, Mat, Mat, Mat *, void *), void *ctx)
139: {
140: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
142: PetscFunctionBegin;
143: jac->computeasub = computeAsub;
144: jac->computeasub_ctx = ctx;
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: /*@
149: PCGalerkinSetRestriction - Sets the restriction operator for the `PCGALERKIN` preconditioner
151: Logically Collective
153: Input Parameters:
154: + pc - the preconditioner context
155: - R - the restriction operator
157: Level: intermediate
159: Note:
160: Either this or `PCGalerkinSetInterpolation()` or both must be called
162: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
163: `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
164: @*/
165: PetscErrorCode PCGalerkinSetRestriction(PC pc, Mat R)
166: {
167: PetscFunctionBegin;
169: PetscTryMethod(pc, "PCGalerkinSetRestriction_C", (PC, Mat), (pc, R));
170: PetscFunctionReturn(PETSC_SUCCESS);
171: }
173: /*@
174: PCGalerkinSetInterpolation - Sets the interpolation operator for the `PCGALERKIN` preconditioner
176: Logically Collective
178: Input Parameters:
179: + pc - the preconditioner context
180: - P - the interpolation operator
182: Level: intermediate
184: Note:
185: Either this or `PCGalerkinSetRestriction()` or both must be called
187: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
188: `PCGalerkinSetRestriction()`, `PCGalerkinGetKSP()`
189: @*/
190: PetscErrorCode PCGalerkinSetInterpolation(PC pc, Mat P)
191: {
192: PetscFunctionBegin;
194: PetscTryMethod(pc, "PCGalerkinSetInterpolation_C", (PC, Mat), (pc, P));
195: PetscFunctionReturn(PETSC_SUCCESS);
196: }
198: /*@C
199: PCGalerkinSetComputeSubmatrix - Provide a routine that will be called to compute the Galerkin submatrix
201: Logically Collective
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`
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
215: Level: intermediate
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.
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.
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()`
228: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
229: `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
230: @*/
231: PetscErrorCode PCGalerkinSetComputeSubmatrix(PC pc, PetscErrorCode (*computeAsub)(PC pc, Mat A, Mat Ap, Mat *cAp, void *ctx), void *ctx)
232: {
233: PetscFunctionBegin;
235: PetscTryMethod(pc, "PCGalerkinSetComputeSubmatrix_C", (PC, PetscErrorCode (*)(PC, Mat, Mat, Mat *, void *), void *), (pc, computeAsub, ctx));
236: PetscFunctionReturn(PETSC_SUCCESS);
237: }
239: /*@
240: PCGalerkinGetKSP - Gets the `KSP` object in the `PCGALERKIN`
242: Not Collective
244: Input Parameter:
245: . pc - the preconditioner context
247: Output Parameter:
248: . ksp - the `KSP` object
250: Level: intermediate
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.
256: .seealso: [](ch_ksp), `PC`, `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
257: `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinSetComputeSubmatrix()`
258: @*/
259: PetscErrorCode PCGalerkinGetKSP(PC pc, KSP *ksp)
260: {
261: PetscFunctionBegin;
263: PetscAssertPointer(ksp, 2);
264: PetscUseMethod(pc, "PCGalerkinGetKSP_C", (PC, KSP *), (pc, ksp));
265: PetscFunctionReturn(PETSC_SUCCESS);
266: }
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;
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: }
283: PetscOptionsHeadBegin(PetscOptionsObject, "Galerkin options");
284: if (jac->ksp) PetscCall(KSPSetFromOptions(jac->ksp));
285: PetscOptionsHeadEnd();
286: PetscFunctionReturn(PETSC_SUCCESS);
287: }
289: /*MC
290: PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T)
292: Level: intermediate
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
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.
307: Should there be a prefix for the inner `KSP`?
309: There is no `KSPSetFromOptions_Galerkin()` that calls `KSPSetFromOptions()` on the inner `KSP`
311: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
312: `PCSHELL`, `PCKSP`, `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
313: M*/
315: PETSC_EXTERN PetscErrorCode PCCreate_Galerkin(PC pc)
316: {
317: PC_Galerkin *jac;
319: PetscFunctionBegin;
320: PetscCall(PetscNew(&jac));
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;
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));
335: pc->data = (void *)jac;
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: }