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:   PetscFunctionBegin;
 17:   PetscCall(KSPCreate(PetscObjectComm((PetscObject)pc), &jac->ksp));
 18:   PetscCall(KSPSetErrorIfNotConverged(jac->ksp, pc->erroriffailure));
 19:   PetscCall(PetscObjectIncrementTabLevel((PetscObject)jac->ksp, (PetscObject)pc, 1));
 20:   PetscCall(PCGetOptionsPrefix(pc, &prefix));
 21:   PetscCall(KSPSetOptionsPrefix(jac->ksp, prefix));
 22:   PetscCall(KSPAppendOptionsPrefix(jac->ksp, "ksp_"));
 23:   PetscCall(PCGetDM(pc, &dm));
 24:   if (dm) {
 25:     PetscCall(KSPSetDM(jac->ksp, dm));
 26:     PetscCall(KSPSetDMActive(jac->ksp, PETSC_FALSE));
 27:   }
 28:   PetscFunctionReturn(PETSC_SUCCESS);
 29: }

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

 36:   PetscFunctionBegin;
 37:   if (jac->ksp->presolve) {
 38:     PetscCall(VecCopy(x, y));
 39:     PetscCall(KSPSolve(jac->ksp, y, y));
 40:   } else {
 41:     PetscCall(KSPSolve(jac->ksp, x, y));
 42:   }
 43:   PetscCall(KSPCheckSolve(jac->ksp, pc, y));
 44:   PetscCall(KSPGetIterationNumber(jac->ksp, &its));
 45:   jac->its += its;
 46:   PetscFunctionReturn(PETSC_SUCCESS);
 47: }

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

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

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

 72:   PetscFunctionBegin;
 73:   if (jac->ksp->presolve) {
 74:     PetscCall(VecCopy(x, y));
 75:     PetscCall(KSPSolve(jac->ksp, y, y));
 76:   } else {
 77:     PetscCall(KSPSolveTranspose(jac->ksp, x, y));
 78:   }
 79:   PetscCall(KSPCheckSolve(jac->ksp, pc, y));
 80:   PetscCall(KSPGetIterationNumber(jac->ksp, &its));
 81:   jac->its += its;
 82:   PetscFunctionReturn(PETSC_SUCCESS);
 83: }

 85: static PetscErrorCode PCSetUp_KSP(PC pc)
 86: {
 87:   PC_KSP *jac = (PC_KSP *)pc->data;
 88:   Mat     mat;

 90:   PetscFunctionBegin;
 91:   if (!jac->ksp) {
 92:     PetscCall(PCKSPCreateKSP_KSP(pc));
 93:     PetscCall(KSPSetFromOptions(jac->ksp));
 94:   }
 95:   if (pc->useAmat) mat = pc->mat;
 96:   else mat = pc->pmat;
 97:   PetscCall(KSPSetOperators(jac->ksp, mat, pc->pmat));
 98:   PetscCall(KSPSetUp(jac->ksp));
 99:   PetscFunctionReturn(PETSC_SUCCESS);
100: }

102: /* Default destroy, if it has never been setup */
103: static PetscErrorCode PCReset_KSP(PC pc)
104: {
105:   PC_KSP *jac = (PC_KSP *)pc->data;

107:   PetscFunctionBegin;
108:   PetscCall(KSPDestroy(&jac->ksp));
109:   PetscFunctionReturn(PETSC_SUCCESS);
110: }

112: static PetscErrorCode PCDestroy_KSP(PC pc)
113: {
114:   PC_KSP *jac = (PC_KSP *)pc->data;

116:   PetscFunctionBegin;
117:   PetscCall(KSPDestroy(&jac->ksp));
118:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCKSPGetKSP_C", NULL));
119:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCKSPSetKSP_C", NULL));
120:   PetscCall(PetscFree(pc->data));
121:   PetscFunctionReturn(PETSC_SUCCESS);
122: }

124: static PetscErrorCode PCView_KSP(PC pc, PetscViewer viewer)
125: {
126:   PC_KSP   *jac = (PC_KSP *)pc->data;
127:   PetscBool iascii;

129:   PetscFunctionBegin;
130:   if (!jac->ksp) PetscCall(PCKSPCreateKSP_KSP(pc));
131:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
132:   if (iascii) {
133:     if (pc->useAmat) PetscCall(PetscViewerASCIIPrintf(viewer, "  Using Amat (not Pmat) as operator on inner solve\n"));
134:     PetscCall(PetscViewerASCIIPrintf(viewer, "  KSP and PC on KSP preconditioner follow\n"));
135:     PetscCall(PetscViewerASCIIPrintf(viewer, "  ---------------------------------\n"));
136:   }
137:   PetscCall(PetscViewerASCIIPushTab(viewer));
138:   PetscCall(KSPView(jac->ksp, viewer));
139:   PetscCall(PetscViewerASCIIPopTab(viewer));
140:   if (iascii) PetscCall(PetscViewerASCIIPrintf(viewer, "  ---------------------------------\n"));
141:   PetscFunctionReturn(PETSC_SUCCESS);
142: }

144: static PetscErrorCode PCKSPSetKSP_KSP(PC pc, KSP ksp)
145: {
146:   PC_KSP *jac = (PC_KSP *)pc->data;

148:   PetscFunctionBegin;
149:   PetscCall(PetscObjectReference((PetscObject)ksp));
150:   PetscCall(KSPDestroy(&jac->ksp));
151:   jac->ksp = ksp;
152:   PetscFunctionReturn(PETSC_SUCCESS);
153: }

155: /*@
156:    PCKSPSetKSP - Sets the `KSP` context for a `PCKSP`.

158:    Collective

160:    Input Parameters:
161: +  pc - the preconditioner context
162: -  ksp - the `KSP` solver

164:    Level: advanced

166:    Notes:
167:    The `PC` and the `KSP` must have the same communicator

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

171: .seealso: `PCKSP`, `PCKSPGetKSP()`
172: @*/
173: PetscErrorCode PCKSPSetKSP(PC pc, KSP ksp)
174: {
175:   PetscFunctionBegin;
178:   PetscCheckSameComm(pc, 1, ksp, 2);
179:   PetscTryMethod(pc, "PCKSPSetKSP_C", (PC, KSP), (pc, ksp));
180:   PetscFunctionReturn(PETSC_SUCCESS);
181: }

183: static PetscErrorCode PCKSPGetKSP_KSP(PC pc, KSP *ksp)
184: {
185:   PC_KSP *jac = (PC_KSP *)pc->data;

187:   PetscFunctionBegin;
188:   if (!jac->ksp) PetscCall(PCKSPCreateKSP_KSP(pc));
189:   *ksp = jac->ksp;
190:   PetscFunctionReturn(PETSC_SUCCESS);
191: }

193: /*@
194:    PCKSPGetKSP - Gets the `KSP` context for a `PCKSP`.

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

198:    Input Parameter:
199: .  pc - the preconditioner context

201:    Output Parameter:
202: .  ksp - the `KSP` solver

204:    Note:
205:    If the `PC` is not a `PCKSP` object it raises an error

207:    Level: advanced

209: .seealso: `PCKSP`, `PCKSPSetKSP()`
210: @*/
211: PetscErrorCode PCKSPGetKSP(PC pc, KSP *ksp)
212: {
213:   PetscFunctionBegin;
216:   PetscUseMethod(pc, "PCKSPGetKSP_C", (PC, KSP *), (pc, ksp));
217:   PetscFunctionReturn(PETSC_SUCCESS);
218: }

220: static PetscErrorCode PCSetFromOptions_KSP(PC pc, PetscOptionItems *PetscOptionsObject)
221: {
222:   PC_KSP *jac = (PC_KSP *)pc->data;

224:   PetscFunctionBegin;
225:   PetscOptionsHeadBegin(PetscOptionsObject, "PC KSP options");
226:   if (jac->ksp) PetscCall(KSPSetFromOptions(jac->ksp));
227:   PetscOptionsHeadEnd();
228:   PetscFunctionReturn(PETSC_SUCCESS);
229: }

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

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

240:    Level: intermediate

242:    Note:
243:     The application of an inexact Krylov solve is a nonlinear operation. Thus, performing a solve with `KSP` is,
244:     in general, a nonlinear operation, so `PCKSP` is in general a nonlinear preconditioner.
245:     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.

247:    Developer Note:
248:     If the outer Krylov method has a nonzero initial guess it will compute a new residual based on that initial guess
249:     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
250:     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
251:     input (current residual) vector, the action of the preconditioner cannot depend also on some other vector (the "initial guess"). For
252:     `KSPRICHARDSON` it is possible to provide a `PCApplyRichardson_PCKSP()` that short circuits returning to the `KSP` object at each iteration to compute the
253:     residual, see for example `PCApplyRichardson_SOR()`. We do not implement a `PCApplyRichardson_PCKSP()`  because (1) using a `KSP` directly inside a Richardson
254:     is not an efficient algorithm anyways and (2) implementing it for its > 1 would essentially require that we implement Richardson (reimplementing the
255:     Richardson code) inside the `PCApplyRichardson_PCKSP()` leading to duplicate code.

257: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
258:           `PCSHELL`, `PCCOMPOSITE`, `PCSetUseAmat()`, `PCKSPGetKSP()`, `KSPFGMRES`, `KSPGCR`, `KSPFCG`
259: M*/

261: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
262: {
263:   PC_KSP *jac;

265:   PetscFunctionBegin;
266:   PetscCall(PetscNew(&jac));
267:   pc->data = (void *)jac;

269:   PetscCall(PetscMemzero(pc->ops, sizeof(struct _PCOps)));
270:   pc->ops->apply          = PCApply_KSP;
271:   pc->ops->matapply       = PCMatApply_KSP;
272:   pc->ops->applytranspose = PCApplyTranspose_KSP;
273:   pc->ops->setup          = PCSetUp_KSP;
274:   pc->ops->reset          = PCReset_KSP;
275:   pc->ops->destroy        = PCDestroy_KSP;
276:   pc->ops->setfromoptions = PCSetFromOptions_KSP;
277:   pc->ops->view           = PCView_KSP;

279:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCKSPGetKSP_C", PCKSPGetKSP_KSP));
280:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCKSPSetKSP_C", PCKSPSetKSP_KSP));
281:   PetscFunctionReturn(PETSC_SUCCESS);
282: }