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) {
126:       PetscViewerASCIIPrintf(viewer,"  Using Amat (not Pmat) as operator on inner solve\n");
127:     }
128:     PetscViewerASCIIPrintf(viewer,"  KSP and PC on KSP preconditioner follow\n");
129:     PetscViewerASCIIPrintf(viewer,"  ---------------------------------\n");
130:   }
131:   PetscViewerASCIIPushTab(viewer);
132:   KSPView(jac->ksp,viewer);
133:   PetscViewerASCIIPopTab(viewer);
134:   if (iascii) {
135:     PetscViewerASCIIPrintf(viewer,"  ---------------------------------\n");
136:   }
137:   return 0;
138: }

140: static PetscErrorCode  PCKSPSetKSP_KSP(PC pc,KSP ksp)
141: {
142:   PC_KSP         *jac = (PC_KSP*)pc->data;

144:   PetscObjectReference((PetscObject)ksp);
145:   KSPDestroy(&jac->ksp);
146:   jac->ksp = ksp;
147:   return 0;
148: }

150: /*@
151:    PCKSPSetKSP - Sets the KSP context for a KSP PC.

153:    Collective on PC

155:    Input Parameters:
156: +  pc - the preconditioner context
157: -  ksp - the KSP solver

159:    Notes:
160:    The PC and the KSP must have the same communicator

162:    Level: advanced

164: @*/
165: PetscErrorCode  PCKSPSetKSP(PC pc,KSP ksp)
166: {
170:   PetscTryMethod(pc,"PCKSPSetKSP_C",(PC,KSP),(pc,ksp));
171:   return 0;
172: }

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

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

183: /*@
184:    PCKSPGetKSP - Gets the KSP context for a KSP PC.

186:    Not Collective but KSP returned is parallel if PC was parallel

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

191:    Output Parameters:
192: .  ksp - the KSP solver

194:    Notes:
195:    You must call KSPSetUp() before calling PCKSPGetKSP().

197:    If the PC is not a PCKSP object it raises an error

199:    Level: advanced

201: @*/
202: PetscErrorCode  PCKSPGetKSP(PC pc,KSP *ksp)
203: {
206:   PetscUseMethod(pc,"PCKSPGetKSP_C",(PC,KSP*),(pc,ksp));
207:   return 0;
208: }

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

214:   PetscOptionsHead(PetscOptionsObject,"PC KSP options");
215:   if (jac->ksp) {
216:     KSPSetFromOptions(jac->ksp);
217:    }
218:   PetscOptionsTail();
219:   return 0;
220: }

222: /* ----------------------------------------------------------------------------------*/

224: /*MC
225:      PCKSP -    Defines a preconditioner that can consist of any KSP solver.
226:                  This allows, for example, embedding a Krylov method inside a preconditioner.

228:    Options Database Key:
229: .     -pc_use_amat - use the matrix that defines the linear system, Amat as the matrix for the
230:                     inner solver, otherwise by default it uses the matrix used to construct
231:                     the preconditioner, Pmat (see PCSetOperators())

233:    Level: intermediate

235:    Notes:
236:     The application of an inexact Krylov solve is a nonlinear operation. Thus, performing a solve with KSP is,
237:     in general, a nonlinear operation, so PCKSP is in general a nonlinear preconditioner.
238:     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.

240:    Developer Notes:
241:     If the outer Krylov method has a nonzero initial guess it will compute a new residual based on that initial guess
242:     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
243:     except Richardson this is neccessary since Krylov methods, even the flexible ones, need to "see" the result of the action of the preconditioner on the
244:     input (current residual) vector, the action of the preconditioner cannot depend also on some other vector (the "initial guess"). For
245:     KSPRICHARDSON it is possible to provide a PCApplyRichardson_PCKSP() that short circuits returning to the KSP object at each iteration to compute the
246:     residual, see for example PCApplyRichardson_SOR(). We do not implement a PCApplyRichardson_PCKSP()  because (1) using a KSP directly inside a Richardson
247:     is not an efficient algorithm anyways and (2) implementing it for its > 1 would essentially require that we implement Richardson (reimplementing the
248:     Richardson code) inside the PCApplyRichardson_PCKSP() leading to duplicate code.

250: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
251:            PCSHELL, PCCOMPOSITE, PCSetUseAmat(), PCKSPGetKSP()

253: M*/

255: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
256: {
257:   PC_KSP         *jac;

259:   PetscNewLog(pc,&jac);
260:   pc->data = (void*)jac;

262:   PetscMemzero(pc->ops,sizeof(struct _PCOps));
263:   pc->ops->apply           = PCApply_KSP;
264:   pc->ops->matapply        = PCMatApply_KSP;
265:   pc->ops->applytranspose  = PCApplyTranspose_KSP;
266:   pc->ops->setup           = PCSetUp_KSP;
267:   pc->ops->reset           = PCReset_KSP;
268:   pc->ops->destroy         = PCDestroy_KSP;
269:   pc->ops->setfromoptions  = PCSetFromOptions_KSP;
270:   pc->ops->view            = PCView_KSP;

272:   PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",PCKSPGetKSP_KSP);
273:   PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",PCKSPSetKSP_KSP);
274:   return 0;
275: }