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