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: }