Actual source code: pcksp.c

petsc-3.10.4 2019-02-26
Report Typos and Errors
  1:  #include <petsc/private/pcimpl.h>
  2:  #include <petscksp.h>

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

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

 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:   return(0);
 23: }

 25: static PetscErrorCode PCApply_KSP(PC pc,Vec x,Vec y)
 26: {
 27:   PetscErrorCode     ierr;
 28:   PetscInt           its;
 29:   PC_KSP             *jac = (PC_KSP*)pc->data;
 30:   KSPConvergedReason reason;

 33:   KSPSolve(jac->ksp,x,y);
 34:   KSPGetConvergedReason(jac->ksp,&reason);
 35:   if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
 36:     pc->failedreason = PC_SUBPC_ERROR;
 37:   }
 38:   KSPGetIterationNumber(jac->ksp,&its);
 39:   jac->its += its;
 40:   return(0);
 41: }

 43: static PetscErrorCode PCApplyTranspose_KSP(PC pc,Vec x,Vec y)
 44: {
 45:   PetscErrorCode     ierr;
 46:   PetscInt           its;
 47:   PC_KSP             *jac = (PC_KSP*)pc->data;
 48:   KSPConvergedReason reason;

 51:   KSPSolveTranspose(jac->ksp,x,y);
 52:   KSPGetConvergedReason(jac->ksp,&reason);
 53:   if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
 54:     pc->failedreason = PC_SUBPC_ERROR;
 55:   }
 56:   KSPGetIterationNumber(jac->ksp,&its);
 57:   jac->its += its;
 58:   return(0);
 59: }

 61: static PetscErrorCode PCSetUp_KSP(PC pc)
 62: {
 64:   PC_KSP         *jac = (PC_KSP*)pc->data;
 65:   Mat            mat;

 68:   if (!jac->ksp) {
 69:     PCKSPCreateKSP_KSP(pc);
 70:     KSPSetFromOptions(jac->ksp);
 71:   }
 72:   if (pc->useAmat) mat = pc->mat;
 73:   else             mat = pc->pmat;
 74:   KSPSetOperators(jac->ksp,mat,pc->pmat);
 75:   KSPSetUp(jac->ksp);
 76:   return(0);
 77: }

 79: /* Default destroy, if it has never been setup */
 80: static PetscErrorCode PCReset_KSP(PC pc)
 81: {
 82:   PC_KSP         *jac = (PC_KSP*)pc->data;

 86:   KSPDestroy(&jac->ksp);
 87:   return(0);
 88: }

 90: static PetscErrorCode PCDestroy_KSP(PC pc)
 91: {
 92:   PC_KSP         *jac = (PC_KSP*)pc->data;

 96:   KSPDestroy(&jac->ksp);
 97:   PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",NULL);
 98:   PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",NULL);
 99:   PetscFree(pc->data);
100:   return(0);
101: }

103: static PetscErrorCode PCView_KSP(PC pc,PetscViewer viewer)
104: {
105:   PC_KSP         *jac = (PC_KSP*)pc->data;
107:   PetscBool      iascii;

110:   if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
111:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
112:   if (iascii) {
113:     if (pc->useAmat) {
114:       PetscViewerASCIIPrintf(viewer,"  Using Amat (not Pmat) as operator on inner solve\n");
115:     }
116:     PetscViewerASCIIPrintf(viewer,"  KSP and PC on KSP preconditioner follow\n");
117:     PetscViewerASCIIPrintf(viewer,"  ---------------------------------\n");
118:   }
119:   PetscViewerASCIIPushTab(viewer);
120:   KSPView(jac->ksp,viewer);
121:   PetscViewerASCIIPopTab(viewer);
122:   if (iascii) {
123:     PetscViewerASCIIPrintf(viewer,"  ---------------------------------\n");
124:   }
125:   return(0);
126: }

128: static PetscErrorCode  PCKSPSetKSP_KSP(PC pc,KSP ksp)
129: {
130:   PC_KSP         *jac = (PC_KSP*)pc->data;

134:   PetscObjectReference((PetscObject)ksp);
135:   KSPDestroy(&jac->ksp);
136:   jac->ksp = ksp;
137:   return(0);
138: }

140: /*@
141:    PCKSPSetKSP - Sets the KSP context for a KSP PC.

143:    Collective on PC

145:    Input Parameter:
146: +  pc - the preconditioner context
147: -  ksp - the KSP solver

149:    Notes:
150:    The PC and the KSP must have the same communicator

152:    Level: advanced

154: .keywords:  PC, KSP, set, context
155: @*/
156: PetscErrorCode  PCKSPSetKSP(PC pc,KSP ksp)
157: {

164:   PetscTryMethod(pc,"PCKSPSetKSP_C",(PC,KSP),(pc,ksp));
165:   return(0);
166: }

168: static PetscErrorCode  PCKSPGetKSP_KSP(PC pc,KSP *ksp)
169: {
170:   PC_KSP         *jac = (PC_KSP*)pc->data;

174:   if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
175:   *ksp = jac->ksp;
176:   return(0);
177: }

179: /*@
180:    PCKSPGetKSP - Gets the KSP context for a KSP PC.

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

184:    Input Parameter:
185: .  pc - the preconditioner context

187:    Output Parameters:
188: .  ksp - the KSP solver

190:    Notes:
191:    You must call KSPSetUp() before calling PCKSPGetKSP().

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

195:    Level: advanced

197: .keywords:  PC, KSP, get, context
198: @*/
199: PetscErrorCode  PCKSPGetKSP(PC pc,KSP *ksp)
200: {

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;

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

224: /* ----------------------------------------------------------------------------------*/

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

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

235:    Level: intermediate

237:    Concepts: inner iteration

239:    Notes:
240:     Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
241:           the incorrect answer) unless you use KSPFGMRES as the other Krylov method

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

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

256: M*/

258: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
259: {
261:   PC_KSP         *jac;

264:   PetscNewLog(pc,&jac);
265:   pc->data = (void*)jac;

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

276:   PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",PCKSPGetKSP_KSP);
277:   PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",PCKSPSetKSP_KSP);
278:   return(0);
279: }