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: if (jac->R) {
21: MatRestrict(jac->R, x, jac->b);
22: } else {
23: MatRestrict(jac->P, x, jac->b);
24: }
25: KSPSolve(jac->ksp, jac->b, jac->x);
26: KSPCheckSolve(jac->ksp, pc, jac->x);
27: if (jac->P) {
28: MatInterpolate(jac->P, jac->x, y);
29: } else {
30: MatInterpolate(jac->R, jac->x, y);
31: }
32: return 0;
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: if (jac->computeasub) {
42: Mat Ap;
43: if (!pc->setupcalled) {
44: (*jac->computeasub)(pc, pc->pmat, NULL, &Ap, jac->computeasub_ctx);
45: KSPSetOperators(jac->ksp, Ap, Ap);
46: MatDestroy(&Ap);
47: } else {
48: KSPGetOperators(jac->ksp, NULL, &Ap);
49: (*jac->computeasub)(pc, pc->pmat, Ap, NULL, jac->computeasub_ctx);
50: }
51: }
53: if (!jac->x) {
54: KSPGetOperatorsSet(jac->ksp, &a, NULL);
56: KSPCreateVecs(jac->ksp, 1, &xx, 1, &yy);
57: jac->x = *xx;
58: jac->b = *yy;
59: PetscFree(xx);
60: PetscFree(yy);
61: }
63: /* should check here that sizes of R/P match size of a */
65: return 0;
66: }
68: static PetscErrorCode PCReset_Galerkin(PC pc)
69: {
70: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
72: MatDestroy(&jac->R);
73: MatDestroy(&jac->P);
74: VecDestroy(&jac->x);
75: VecDestroy(&jac->b);
76: KSPReset(jac->ksp);
77: return 0;
78: }
80: static PetscErrorCode PCDestroy_Galerkin(PC pc)
81: {
82: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
84: PCReset_Galerkin(pc);
85: KSPDestroy(&jac->ksp);
86: PetscFree(pc->data);
87: return 0;
88: }
90: static PetscErrorCode PCView_Galerkin(PC pc, PetscViewer viewer)
91: {
92: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
93: PetscBool iascii;
95: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
96: if (iascii) {
97: PetscViewerASCIIPrintf(viewer, " KSP on Galerkin follow\n");
98: PetscViewerASCIIPrintf(viewer, " ---------------------------------\n");
99: }
100: KSPView(jac->ksp, viewer);
101: return 0;
102: }
104: static PetscErrorCode PCGalerkinGetKSP_Galerkin(PC pc, KSP *ksp)
105: {
106: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
108: *ksp = jac->ksp;
109: return 0;
110: }
112: static PetscErrorCode PCGalerkinSetRestriction_Galerkin(PC pc, Mat R)
113: {
114: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
116: PetscObjectReference((PetscObject)R);
117: MatDestroy(&jac->R);
118: jac->R = R;
119: return 0;
120: }
122: static PetscErrorCode PCGalerkinSetInterpolation_Galerkin(PC pc, Mat P)
123: {
124: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
126: PetscObjectReference((PetscObject)P);
127: MatDestroy(&jac->P);
128: jac->P = P;
129: return 0;
130: }
132: static PetscErrorCode PCGalerkinSetComputeSubmatrix_Galerkin(PC pc, PetscErrorCode (*computeAsub)(PC, Mat, Mat, Mat *, void *), void *ctx)
133: {
134: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
136: jac->computeasub = computeAsub;
137: jac->computeasub_ctx = ctx;
138: return 0;
139: }
141: /*@
142: PCGalerkinSetRestriction - Sets the restriction operator for the `PCGALERKIN` preconditioner
144: Logically Collective
146: Input Parameters:
147: + pc - the preconditioner context
148: - R - the restriction operator
150: Note:
151: Either this or `PCGalerkinSetInterpolation()` or both must be called
153: Level: Intermediate
155: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
156: `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
157: @*/
158: PetscErrorCode PCGalerkinSetRestriction(PC pc, Mat R)
159: {
161: PetscTryMethod(pc, "PCGalerkinSetRestriction_C", (PC, Mat), (pc, R));
162: return 0;
163: }
165: /*@
166: PCGalerkinSetInterpolation - Sets the interpolation operator for the `PCGALERKIN` preconditioner
168: Logically Collective
170: Input Parameters:
171: + pc - the preconditioner context
172: - R - the interpolation operator
174: Note:
175: Either this or `PCGalerkinSetRestriction()` or both must be called
177: Level: Intermediate
179: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
180: `PCGalerkinSetRestriction()`, `PCGalerkinGetKSP()`
181: @*/
182: PetscErrorCode PCGalerkinSetInterpolation(PC pc, Mat P)
183: {
185: PetscTryMethod(pc, "PCGalerkinSetInterpolation_C", (PC, Mat), (pc, P));
186: return 0;
187: }
189: /*@
190: PCGalerkinSetComputeSubmatrix - Provide a routine that will be called to compute the Galerkin submatrix
192: Logically Collective
194: Input Parameters:
195: + pc - the preconditioner context
196: . computeAsub - routine that computes the submatrix from the global matrix
197: - ctx - context used by the routine, or NULL
199: Calling sequence of computeAsub:
200: $ computeAsub(PC pc,Mat A, Mat Ap, Mat *cAP,void *ctx);
202: + PC - the `PCGALERKIN`
203: . A - the matrix in the `PCGALERKIN`
204: . Ap - the computed submatrix from any previous computation, if NULL it has not previously been computed
205: . cAp - the submatrix computed by this routine
206: - ctx - optional user-defined function context
208: Level: Intermediate
210: Notes:
211: Instead of providing this routine you can call `PCGalerkinGetKSP()` and then `KSPSetOperators()` to provide the submatrix,
212: but that will not work for multiple `KSPSolve()`s with different matrices unless you call it for each solve.
214: 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
215: matrix and computes its values in cAp. On each subsequent call the routine should up the Ap matrix.
217: Developer Note:
218: If the user does not call this routine nor call `PCGalerkinGetKSP()` and `KSPSetOperators()` then `PCGALERKIN`
219: could automatically compute the submatrix via calls to `MatGalerkin()` or `MatRARt()`
221: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
222: `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
223: @*/
224: PetscErrorCode PCGalerkinSetComputeSubmatrix(PC pc, PetscErrorCode (*computeAsub)(PC, Mat, Mat, Mat *, void *), void *ctx)
225: {
227: PetscTryMethod(pc, "PCGalerkinSetComputeSubmatrix_C", (PC, PetscErrorCode(*)(PC, Mat, Mat, Mat *, void *), void *), (pc, computeAsub, ctx));
228: return 0;
229: }
231: /*@
232: PCGalerkinGetKSP - Gets the `KSP` object in the `PCGALERKIN`
234: Not Collective
236: Input Parameter:
237: . pc - the preconditioner context
239: Output Parameter:
240: . ksp - the `KSP` object
242: Level: Intermediate
244: Note:
245: Once you have called this routine you can call `KSPSetOperators()` on the resulting ksp to provide the operator for the Galerkin problem,
246: an alternative is to use `PCGalerkinSetComputeSubmatrix()` to provide a routine that computes the submatrix as needed.
248: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PCGALERKIN`,
249: `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinSetComputeSubmatrix()`
250: @*/
251: PetscErrorCode PCGalerkinGetKSP(PC pc, KSP *ksp)
252: {
255: PetscUseMethod(pc, "PCGalerkinGetKSP_C", (PC, KSP *), (pc, ksp));
256: return 0;
257: }
259: static PetscErrorCode PCSetFromOptions_Galerkin(PC pc, PetscOptionItems *PetscOptionsObject)
260: {
261: PC_Galerkin *jac = (PC_Galerkin *)pc->data;
262: const char *prefix;
263: PetscBool flg;
265: KSPGetOptionsPrefix(jac->ksp, &prefix);
266: PetscStrendswith(prefix, "galerkin_", &flg);
267: if (!flg) {
268: PCGetOptionsPrefix(pc, &prefix);
269: KSPSetOptionsPrefix(jac->ksp, prefix);
270: KSPAppendOptionsPrefix(jac->ksp, "galerkin_");
271: }
273: PetscOptionsHeadBegin(PetscOptionsObject, "Galerkin options");
274: if (jac->ksp) KSPSetFromOptions(jac->ksp);
275: PetscOptionsHeadEnd();
276: return 0;
277: }
279: /*MC
280: PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T)
282: Use `PCGalerkinSetRestriction`(pc,R) and/or `PCGalerkinSetInterpolation`(pc,P) followed by `PCGalerkinGetKSP`(pc,&ksp); `KSPSetOperators`(ksp,A,....)
284: Level: intermediate
286: Developer Notes:
287: If `KSPSetOperators()` has not been called on the inner `KSP` then `PCGALERKIN` could use `MatRARt()` or `MatPtAP()` to compute
288: the operators automatically.
290: Should there be a prefix for the inner `KSP`?
292: There is no `KSPSetFromOptions_Galerkin()` that calls `KSPSetFromOptions()` on the inner `KSP`
294: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
295: `PCSHELL`, `PCKSP`, `PCGalerkinSetRestriction()`, `PCGalerkinSetInterpolation()`, `PCGalerkinGetKSP()`
296: M*/
298: PETSC_EXTERN PetscErrorCode PCCreate_Galerkin(PC pc)
299: {
300: PC_Galerkin *jac;
302: PetscNew(&jac);
304: pc->ops->apply = PCApply_Galerkin;
305: pc->ops->setup = PCSetUp_Galerkin;
306: pc->ops->reset = PCReset_Galerkin;
307: pc->ops->destroy = PCDestroy_Galerkin;
308: pc->ops->view = PCView_Galerkin;
309: pc->ops->setfromoptions = PCSetFromOptions_Galerkin;
310: pc->ops->applyrichardson = NULL;
312: KSPCreate(PetscObjectComm((PetscObject)pc), &jac->ksp);
313: KSPSetErrorIfNotConverged(jac->ksp, pc->erroriffailure);
314: PetscObjectIncrementTabLevel((PetscObject)jac->ksp, (PetscObject)pc, 1);
316: pc->data = (void *)jac;
318: PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetRestriction_C", PCGalerkinSetRestriction_Galerkin);
319: PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetInterpolation_C", PCGalerkinSetInterpolation_Galerkin);
320: PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinGetKSP_C", PCGalerkinGetKSP_Galerkin);
321: PetscObjectComposeFunction((PetscObject)pc, "PCGalerkinSetComputeSubmatrix_C", PCGalerkinSetComputeSubmatrix_Galerkin);
322: return 0;
323: }