Actual source code: pcksp.c

petsc-3.11.2 2019-05-18
Report Typos and Errors
  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: {
 13:   const char     *prefix;
 14:   PC_KSP         *jac = (PC_KSP*)pc->data;

 17:   KSPCreate(PetscObjectComm((PetscObject)pc),&jac->ksp);
 18:   KSPSetErrorIfNotConverged(jac->ksp,pc->erroriffailure);
 19:   PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);
 20:   PCGetOptionsPrefix(pc,&prefix);
 21:   KSPSetOptionsPrefix(jac->ksp,prefix);
 22:   KSPAppendOptionsPrefix(jac->ksp,"ksp_");
 23:   return(0);
 24: }

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

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

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

 52:   if (jac->ksp->presolve) {
 53:     VecCopy(x,y);
 54:     KSPSolve(jac->ksp,y,y);
 55:   } else {
 56:     KSPSolveTranspose(jac->ksp,x,y);
 57:   }
 58:   KSPCheckSolve(jac->ksp,pc,y);
 59:   KSPGetIterationNumber(jac->ksp,&its);
 60:   jac->its += its;
 61:   return(0);
 62: }

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

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

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

 89:   KSPDestroy(&jac->ksp);
 90:   return(0);
 91: }

 93: static PetscErrorCode PCDestroy_KSP(PC pc)
 94: {
 95:   PC_KSP         *jac = (PC_KSP*)pc->data;

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

106: static PetscErrorCode PCView_KSP(PC pc,PetscViewer viewer)
107: {
108:   PC_KSP         *jac = (PC_KSP*)pc->data;
110:   PetscBool      iascii;

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

131: static PetscErrorCode  PCKSPSetKSP_KSP(PC pc,KSP ksp)
132: {
133:   PC_KSP         *jac = (PC_KSP*)pc->data;

137:   PetscObjectReference((PetscObject)ksp);
138:   KSPDestroy(&jac->ksp);
139:   jac->ksp = ksp;
140:   return(0);
141: }

143: /*@
144:    PCKSPSetKSP - Sets the KSP context for a KSP PC.

146:    Collective on PC

148:    Input Parameter:
149: +  pc - the preconditioner context
150: -  ksp - the KSP solver

152:    Notes:
153:    The PC and the KSP must have the same communicator

155:    Level: advanced

157: .keywords:  PC, KSP, set, context
158: @*/
159: PetscErrorCode  PCKSPSetKSP(PC pc,KSP ksp)
160: {

167:   PetscTryMethod(pc,"PCKSPSetKSP_C",(PC,KSP),(pc,ksp));
168:   return(0);
169: }

171: static PetscErrorCode  PCKSPGetKSP_KSP(PC pc,KSP *ksp)
172: {
173:   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 KSP PC.

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

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

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

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

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

198:    Level: advanced

200: .keywords:  PC, KSP, get, context
201: @*/
202: PetscErrorCode  PCKSPGetKSP(PC pc,KSP *ksp)
203: {

209:   PetscUseMethod(pc,"PCKSPGetKSP_C",(PC,KSP*),(pc,ksp));
210:   return(0);
211: }

213: static PetscErrorCode PCSetFromOptions_KSP(PetscOptionItems *PetscOptionsObject,PC pc)
214: {
215:   PC_KSP         *jac = (PC_KSP*)pc->data;

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

227: /* ----------------------------------------------------------------------------------*/

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

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

238:    Level: intermediate

240:    Concepts: inner iteration

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

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

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

259: M*/

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

267:   PetscNewLog(pc,&jac);
268:   pc->data = (void*)jac;

270:   PetscMemzero(pc->ops,sizeof(struct _PCOps));
271:   pc->ops->apply           = PCApply_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:   PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",PCKSPGetKSP_KSP);
280:   PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",PCKSPSetKSP_KSP);
281:   return(0);
282: }