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: /*@
143: PCGalerkinSetRestriction - Sets the restriction operator for the "Galerkin-type" preconditioner
145: Logically Collective on PC
147: Input Parameters:
148: + pc - the preconditioner context
149: - R - the restriction operator
151: Notes:
152: Either this or PCGalerkinSetInterpolation() or both must be called
154: Level: Intermediate
156: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
157: PCGalerkinSetInterpolation(), PCGalerkinGetKSP()
159: @*/
160: PetscErrorCode PCGalerkinSetRestriction(PC pc,Mat R)
161: {
163: PetscTryMethod(pc,"PCGalerkinSetRestriction_C",(PC,Mat),(pc,R));
164: return 0;
165: }
167: /*@
168: PCGalerkinSetInterpolation - Sets the interpolation operator for the "Galerkin-type" preconditioner
170: Logically Collective on PC
172: Input Parameters:
173: + pc - the preconditioner context
174: - R - the interpolation operator
176: Notes:
177: Either this or PCGalerkinSetRestriction() or both must be called
179: Level: Intermediate
181: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
182: PCGalerkinSetRestriction(), PCGalerkinGetKSP()
184: @*/
185: PetscErrorCode PCGalerkinSetInterpolation(PC pc,Mat P)
186: {
188: PetscTryMethod(pc,"PCGalerkinSetInterpolation_C",(PC,Mat),(pc,P));
189: return 0;
190: }
192: /*@
193: PCGalerkinSetComputeSubmatrix - Provide a routine that will be called to compute the Galerkin submatrix
195: Logically Collective
197: Input Parameters:
198: + pc - the preconditioner context
199: . computeAsub - routine that computes the submatrix from the global matrix
200: - ctx - context used by the routine, or NULL
202: Calling sequence of computeAsub:
203: $ computeAsub(PC pc,Mat A, Mat Ap, Mat *cAP,void *ctx);
205: + PC - the Galerkin PC
206: . A - the matrix in the Galerkin PC
207: . Ap - the computed submatrix from any previous computation, if NULL it has not previously been computed
208: . cAp - the submatrix computed by this routine
209: - ctx - optional user-defined function context
211: Level: Intermediate
213: Notes:
214: Instead of providing this routine you can call PCGalerkinGetKSP() and then KSPSetOperators() to provide the submatrix,
215: but that will not work for multiple KSPSolves with different matrices unless you call it for each solve.
217: 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
218: matrix and computes its values in cAp. On each subsequent call the routine should up the Ap matrix.
220: Developer Notes:
221: If the user does not call this routine nor call PCGalerkinGetKSP() and KSPSetOperators() then PCGalerkin could
222: could automatically compute the submatrix via calls to MatGalerkin() or MatRARt()
224: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
225: PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinGetKSP()
227: @*/
228: PetscErrorCode PCGalerkinSetComputeSubmatrix(PC pc,PetscErrorCode (*computeAsub)(PC,Mat,Mat,Mat*,void*),void *ctx)
229: {
231: PetscTryMethod(pc,"PCGalerkinSetComputeSubmatrix_C",(PC,PetscErrorCode (*)(PC,Mat,Mat,Mat*,void*),void*),(pc,computeAsub,ctx));
232: return 0;
233: }
235: /*@
236: PCGalerkinGetKSP - Gets the KSP object in the Galerkin PC.
238: Not Collective
240: Input Parameter:
241: . pc - the preconditioner context
243: Output Parameters:
244: . ksp - the KSP object
246: Level: Intermediate
248: Notes:
249: Once you have called this routine you can call KSPSetOperators() on the resulting ksp to provide the operator for the Galerkin problem,
250: an alternative is to use PCGalerkinSetComputeSubmatrix() to provide a routine that computes the submatrix as needed.
252: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
253: PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinSetComputeSubmatrix()
255: @*/
256: PetscErrorCode PCGalerkinGetKSP(PC pc,KSP *ksp)
257: {
260: PetscUseMethod(pc,"PCGalerkinGetKSP_C",(PC,KSP*),(pc,ksp));
261: return 0;
262: }
264: static PetscErrorCode PCSetFromOptions_Galerkin(PetscOptionItems *PetscOptionsObject,PC pc)
265: {
266: PC_Galerkin *jac = (PC_Galerkin*)pc->data;
267: const char *prefix;
268: PetscBool flg;
270: KSPGetOptionsPrefix(jac->ksp,&prefix);
271: PetscStrendswith(prefix,"galerkin_",&flg);
272: if (!flg) {
273: PCGetOptionsPrefix(pc,&prefix);
274: KSPSetOptionsPrefix(jac->ksp,prefix);
275: KSPAppendOptionsPrefix(jac->ksp,"galerkin_");
276: }
278: PetscOptionsHead(PetscOptionsObject,"Galerkin options");
279: if (jac->ksp) {
280: KSPSetFromOptions(jac->ksp);
281: }
282: PetscOptionsTail();
283: return 0;
284: }
286: /* -------------------------------------------------------------------------------------------*/
288: /*MC
289: PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T)
291: $ Use PCGalerkinSetRestriction(pc,R) and/or PCGalerkinSetInterpolation(pc,P) followed by
292: $ PCGalerkinGetKSP(pc,&ksp); KSPSetOperators(ksp,A,....)
294: Level: intermediate
296: Developer Note: If KSPSetOperators() has not been called on the inner KSP then PCGALERKIN could use MatRARt() or MatPtAP() to compute
297: the operators automatically.
298: Should there be a prefix for the inner KSP.
299: There is no KSPSetFromOptions_Galerkin() that calls KSPSetFromOptions() on the inner KSP
301: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
302: PCSHELL, PCKSP, PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinGetKSP()
304: M*/
306: PETSC_EXTERN PetscErrorCode PCCreate_Galerkin(PC pc)
307: {
308: PC_Galerkin *jac;
310: PetscNewLog(pc,&jac);
312: pc->ops->apply = PCApply_Galerkin;
313: pc->ops->setup = PCSetUp_Galerkin;
314: pc->ops->reset = PCReset_Galerkin;
315: pc->ops->destroy = PCDestroy_Galerkin;
316: pc->ops->view = PCView_Galerkin;
317: pc->ops->setfromoptions = PCSetFromOptions_Galerkin;
318: pc->ops->applyrichardson = NULL;
320: KSPCreate(PetscObjectComm((PetscObject)pc),&jac->ksp);
321: KSPSetErrorIfNotConverged(jac->ksp,pc->erroriffailure);
322: PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);
324: pc->data = (void*)jac;
326: PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetRestriction_C",PCGalerkinSetRestriction_Galerkin);
327: PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetInterpolation_C",PCGalerkinSetInterpolation_Galerkin);
328: PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinGetKSP_C",PCGalerkinGetKSP_Galerkin);
329: PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetComputeSubmatrix_C",PCGalerkinSetComputeSubmatrix_Galerkin);
330: return 0;
331: }