Actual source code: pcksp.c

  1: #include <petsc/private/pcimpl.h>
  2: #include <petsc/private/kspimpl.h>
  3: #include <petscksp.h>

  5: typedef struct {
  6:   KSP      ksp;
  7:   PetscInt its; /* total number of iterations KSP uses */
  8: } PC_KSP;

 10: static PetscErrorCode PCKSPCreateKSP_KSP(PC pc)
 11: {
 12:   const char *prefix;
 13:   PC_KSP     *jac = (PC_KSP *)pc->data;
 14:   DM          dm;

 16:   KSPCreate(PetscObjectComm((PetscObject)pc), &jac->ksp);
 17:   KSPSetErrorIfNotConverged(jac->ksp, pc->erroriffailure);
 18:   PetscObjectIncrementTabLevel((PetscObject)jac->ksp, (PetscObject)pc, 1);
 19:   PCGetOptionsPrefix(pc, &prefix);
 20:   KSPSetOptionsPrefix(jac->ksp, prefix);
 21:   KSPAppendOptionsPrefix(jac->ksp, "ksp_");
 22:   PCGetDM(pc, &dm);
 23:   if (dm) {
 24:     KSPSetDM(jac->ksp, dm);
 25:     KSPSetDMActive(jac->ksp, PETSC_FALSE);
 26:   }
 27:   return 0;
 28: }

 30: static PetscErrorCode PCApply_KSP(PC pc, Vec x, Vec y)
 31: {
 32:   PetscInt its;
 33:   PC_KSP  *jac = (PC_KSP *)pc->data;

 35:   if (jac->ksp->presolve) {
 36:     VecCopy(x, y);
 37:     KSPSolve(jac->ksp, y, y);
 38:   } else {
 39:     KSPSolve(jac->ksp, x, y);
 40:   }
 41:   KSPCheckSolve(jac->ksp, pc, y);
 42:   KSPGetIterationNumber(jac->ksp, &its);
 43:   jac->its += its;
 44:   return 0;
 45: }

 47: static PetscErrorCode PCMatApply_KSP(PC pc, Mat X, Mat Y)
 48: {
 49:   PetscInt its;
 50:   PC_KSP  *jac = (PC_KSP *)pc->data;

 52:   if (jac->ksp->presolve) {
 53:     MatCopy(X, Y, SAME_NONZERO_PATTERN);
 54:     KSPMatSolve(jac->ksp, Y, Y); /* TODO FIXME: this will fail since KSPMatSolve does not allow inplace solve yet */
 55:   } else {
 56:     KSPMatSolve(jac->ksp, X, Y);
 57:   }
 58:   KSPCheckSolve(jac->ksp, pc, NULL);
 59:   KSPGetIterationNumber(jac->ksp, &its);
 60:   jac->its += its;
 61:   return 0;
 62: }

 64: static PetscErrorCode PCApplyTranspose_KSP(PC pc, Vec x, Vec y)
 65: {
 66:   PetscInt its;
 67:   PC_KSP  *jac = (PC_KSP *)pc->data;

 69:   if (jac->ksp->presolve) {
 70:     VecCopy(x, y);
 71:     KSPSolve(jac->ksp, y, y);
 72:   } else {
 73:     KSPSolveTranspose(jac->ksp, x, y);
 74:   }
 75:   KSPCheckSolve(jac->ksp, pc, y);
 76:   KSPGetIterationNumber(jac->ksp, &its);
 77:   jac->its += its;
 78:   return 0;
 79: }

 81: static PetscErrorCode PCSetUp_KSP(PC pc)
 82: {
 83:   PC_KSP *jac = (PC_KSP *)pc->data;
 84:   Mat     mat;

 86:   if (!jac->ksp) {
 87:     PCKSPCreateKSP_KSP(pc);
 88:     KSPSetFromOptions(jac->ksp);
 89:   }
 90:   if (pc->useAmat) mat = pc->mat;
 91:   else mat = pc->pmat;
 92:   KSPSetOperators(jac->ksp, mat, pc->pmat);
 93:   KSPSetUp(jac->ksp);
 94:   return 0;
 95: }

 97: /* Default destroy, if it has never been setup */
 98: static PetscErrorCode PCReset_KSP(PC pc)
 99: {
100:   PC_KSP *jac = (PC_KSP *)pc->data;

102:   KSPDestroy(&jac->ksp);
103:   return 0;
104: }

106: static PetscErrorCode PCDestroy_KSP(PC pc)
107: {
108:   PC_KSP *jac = (PC_KSP *)pc->data;

110:   KSPDestroy(&jac->ksp);
111:   PetscObjectComposeFunction((PetscObject)pc, "PCKSPGetKSP_C", NULL);
112:   PetscObjectComposeFunction((PetscObject)pc, "PCKSPSetKSP_C", NULL);
113:   PetscFree(pc->data);
114:   return 0;
115: }

117: static PetscErrorCode PCView_KSP(PC pc, PetscViewer viewer)
118: {
119:   PC_KSP   *jac = (PC_KSP *)pc->data;
120:   PetscBool iascii;

122:   if (!jac->ksp) PCKSPCreateKSP_KSP(pc);
123:   PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
124:   if (iascii) {
125:     if (pc->useAmat) PetscViewerASCIIPrintf(viewer, "  Using Amat (not Pmat) as operator on inner solve\n");
126:     PetscViewerASCIIPrintf(viewer, "  KSP and PC on KSP preconditioner follow\n");
127:     PetscViewerASCIIPrintf(viewer, "  ---------------------------------\n");
128:   }
129:   PetscViewerASCIIPushTab(viewer);
130:   KSPView(jac->ksp, viewer);
131:   PetscViewerASCIIPopTab(viewer);
132:   if (iascii) PetscViewerASCIIPrintf(viewer, "  ---------------------------------\n");
133:   return 0;
134: }

136: static PetscErrorCode PCKSPSetKSP_KSP(PC pc, KSP ksp)
137: {
138:   PC_KSP *jac = (PC_KSP *)pc->data;

140:   PetscObjectReference((PetscObject)ksp);
141:   KSPDestroy(&jac->ksp);
142:   jac->ksp = ksp;
143:   return 0;
144: }

146: /*@
147:    PCKSPSetKSP - Sets the `KSP` context for a `PCKSP`.

149:    Collective

151:    Input Parameters:
152: +  pc - the preconditioner context
153: -  ksp - the `KSP` solver

155:    Level: advanced

157:    Notes:
158:    The `PC` and the `KSP` must have the same communicator

160:    This would rarely be used, the standard usage is to call `PCKSPGetKSP()` and then change options on that `KSP`

162: .seealso: `PCKSP`, `PCKSPGetKSP()`
163: @*/
164: PetscErrorCode PCKSPSetKSP(PC pc, KSP ksp)
165: {
169:   PetscTryMethod(pc, "PCKSPSetKSP_C", (PC, KSP), (pc, ksp));
170:   return 0;
171: }

173: static PetscErrorCode PCKSPGetKSP_KSP(PC pc, KSP *ksp)
174: {
175:   PC_KSP *jac = (PC_KSP *)pc->data;

177:   if (!jac->ksp) PCKSPCreateKSP_KSP(pc);
178:   *ksp = jac->ksp;
179:   return 0;
180: }

182: /*@
183:    PCKSPGetKSP - Gets the `KSP` context for a `PCKSP`.

185:    Not Collective but ksp returned is parallel if pc was parallel

187:    Input Parameter:
188: .  pc - the preconditioner context

190:    Output Parameter:
191: .  ksp - the `KSP` solver

193:    Note:
194:    If the `PC` is not a `PCKSP` object it raises an error

196:    Level: advanced

198: .seealso: `PCKSP`, `PCKSPSetKSP()`
199: @*/
200: PetscErrorCode PCKSPGetKSP(PC pc, KSP *ksp)
201: {
204:   PetscUseMethod(pc, "PCKSPGetKSP_C", (PC, KSP *), (pc, ksp));
205:   return 0;
206: }

208: static PetscErrorCode PCSetFromOptions_KSP(PC pc, PetscOptionItems *PetscOptionsObject)
209: {
210:   PC_KSP *jac = (PC_KSP *)pc->data;

212:   PetscOptionsHeadBegin(PetscOptionsObject, "PC KSP options");
213:   if (jac->ksp) KSPSetFromOptions(jac->ksp);
214:   PetscOptionsHeadEnd();
215:   return 0;
216: }

218: /*MC
219:      PCKSP -    Defines a preconditioner as any `KSP` solver.
220:                  This allows, for example, embedding a Krylov method inside a preconditioner.

222:    Options Database Key:
223: .   -pc_use_amat - use the matrix that defines the linear system, Amat as the matrix for the
224:                     inner solver, otherwise by default it uses the matrix used to construct
225:                     the preconditioner, Pmat (see `PCSetOperators()`)

227:    Level: intermediate

229:    Note:
230:     The application of an inexact Krylov solve is a nonlinear operation. Thus, performing a solve with `KSP` is,
231:     in general, a nonlinear operation, so `PCKSP` is in general a nonlinear preconditioner.
232:     Thus, one can see divergence or an incorrect answer unless using a flexible Krylov method (e.g. `KSPFGMRES`, `KSPGCR`, or `KSPFCG`) for the outer Krylov solve.

234:    Developer Note:
235:     If the outer Krylov method has a nonzero initial guess it will compute a new residual based on that initial guess
236:     and pass that as the right hand side into this `KSP` (and hence this `KSP` will always have a zero initial guess). For all outer Krylov methods
237:     except Richardson this is necessary since Krylov methods, even the flexible ones, need to "see" the result of the action of the preconditioner on the
238:     input (current residual) vector, the action of the preconditioner cannot depend also on some other vector (the "initial guess"). For
239:     `KSPRICHARDSON` it is possible to provide a `PCApplyRichardson_PCKSP()` that short circuits returning to the `KSP` object at each iteration to compute the
240:     residual, see for example `PCApplyRichardson_SOR()`. We do not implement a `PCApplyRichardson_PCKSP()`  because (1) using a `KSP` directly inside a Richardson
241:     is not an efficient algorithm anyways and (2) implementing it for its > 1 would essentially require that we implement Richardson (reimplementing the
242:     Richardson code) inside the `PCApplyRichardson_PCKSP()` leading to duplicate code.

244: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
245:           `PCSHELL`, `PCCOMPOSITE`, `PCSetUseAmat()`, `PCKSPGetKSP()`, `KSPFGMRES`, `KSPGCR`, `KSPFCG`
246: M*/

248: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
249: {
250:   PC_KSP *jac;

252:   PetscNew(&jac);
253:   pc->data = (void *)jac;

255:   PetscMemzero(pc->ops, sizeof(struct _PCOps));
256:   pc->ops->apply          = PCApply_KSP;
257:   pc->ops->matapply       = PCMatApply_KSP;
258:   pc->ops->applytranspose = PCApplyTranspose_KSP;
259:   pc->ops->setup          = PCSetUp_KSP;
260:   pc->ops->reset          = PCReset_KSP;
261:   pc->ops->destroy        = PCDestroy_KSP;
262:   pc->ops->setfromoptions = PCSetFromOptions_KSP;
263:   pc->ops->view           = PCView_KSP;

265:   PetscObjectComposeFunction((PetscObject)pc, "PCKSPGetKSP_C", PCKSPGetKSP_KSP);
266:   PetscObjectComposeFunction((PetscObject)pc, "PCKSPSetKSP_C", PCKSPSetKSP_KSP);
267:   return 0;
268: }