Actual source code: pcksp.c

petsc-3.14.0 2020-09-29
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;
 15:   DM             dm;

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

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

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

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

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

 70: static PetscErrorCode PCApplyTranspose_KSP(PC pc,Vec x,Vec y)
 71: {
 72:   PetscErrorCode     ierr;
 73:   PetscInt           its;
 74:   PC_KSP             *jac = (PC_KSP*)pc->data;

 77:   if (jac->ksp->presolve) {
 78:     VecCopy(x,y);
 79:     KSPSolve(jac->ksp,y,y);
 80:   } else {
 81:     KSPSolveTranspose(jac->ksp,x,y);
 82:   }
 83:   KSPCheckSolve(jac->ksp,pc,y);
 84:   KSPGetIterationNumber(jac->ksp,&its);
 85:   jac->its += its;
 86:   return(0);
 87: }

 89: static PetscErrorCode PCSetUp_KSP(PC pc)
 90: {
 92:   PC_KSP         *jac = (PC_KSP*)pc->data;
 93:   Mat            mat;

 96:   if (!jac->ksp) {
 97:     PCKSPCreateKSP_KSP(pc);
 98:     KSPSetFromOptions(jac->ksp);
 99:   }
100:   if (pc->useAmat) mat = pc->mat;
101:   else             mat = pc->pmat;
102:   KSPSetOperators(jac->ksp,mat,pc->pmat);
103:   KSPSetUp(jac->ksp);
104:   return(0);
105: }

107: /* Default destroy, if it has never been setup */
108: static PetscErrorCode PCReset_KSP(PC pc)
109: {
110:   PC_KSP         *jac = (PC_KSP*)pc->data;

114:   KSPDestroy(&jac->ksp);
115:   return(0);
116: }

118: static PetscErrorCode PCDestroy_KSP(PC pc)
119: {
120:   PC_KSP         *jac = (PC_KSP*)pc->data;

124:   KSPDestroy(&jac->ksp);
125:   PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",NULL);
126:   PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",NULL);
127:   PetscFree(pc->data);
128:   return(0);
129: }

131: static PetscErrorCode PCView_KSP(PC pc,PetscViewer viewer)
132: {
133:   PC_KSP         *jac = (PC_KSP*)pc->data;
135:   PetscBool      iascii;

138:   if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
139:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
140:   if (iascii) {
141:     if (pc->useAmat) {
142:       PetscViewerASCIIPrintf(viewer,"  Using Amat (not Pmat) as operator on inner solve\n");
143:     }
144:     PetscViewerASCIIPrintf(viewer,"  KSP and PC on KSP preconditioner follow\n");
145:     PetscViewerASCIIPrintf(viewer,"  ---------------------------------\n");
146:   }
147:   PetscViewerASCIIPushTab(viewer);
148:   KSPView(jac->ksp,viewer);
149:   PetscViewerASCIIPopTab(viewer);
150:   if (iascii) {
151:     PetscViewerASCIIPrintf(viewer,"  ---------------------------------\n");
152:   }
153:   return(0);
154: }

156: static PetscErrorCode  PCKSPSetKSP_KSP(PC pc,KSP ksp)
157: {
158:   PC_KSP         *jac = (PC_KSP*)pc->data;

162:   PetscObjectReference((PetscObject)ksp);
163:   KSPDestroy(&jac->ksp);
164:   jac->ksp = ksp;
165:   return(0);
166: }

168: /*@
169:    PCKSPSetKSP - Sets the KSP context for a KSP PC.

171:    Collective on PC

173:    Input Parameter:
174: +  pc - the preconditioner context
175: -  ksp - the KSP solver

177:    Notes:
178:    The PC and the KSP must have the same communicator

180:    Level: advanced

182: @*/
183: PetscErrorCode  PCKSPSetKSP(PC pc,KSP ksp)
184: {

191:   PetscTryMethod(pc,"PCKSPSetKSP_C",(PC,KSP),(pc,ksp));
192:   return(0);
193: }

195: static PetscErrorCode  PCKSPGetKSP_KSP(PC pc,KSP *ksp)
196: {
197:   PC_KSP         *jac = (PC_KSP*)pc->data;

201:   if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
202:   *ksp = jac->ksp;
203:   return(0);
204: }

206: /*@
207:    PCKSPGetKSP - Gets the KSP context for a KSP PC.

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

211:    Input Parameter:
212: .  pc - the preconditioner context

214:    Output Parameters:
215: .  ksp - the KSP solver

217:    Notes:
218:    You must call KSPSetUp() before calling PCKSPGetKSP().

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

222:    Level: advanced

224: @*/
225: PetscErrorCode  PCKSPGetKSP(PC pc,KSP *ksp)
226: {

232:   PetscUseMethod(pc,"PCKSPGetKSP_C",(PC,KSP*),(pc,ksp));
233:   return(0);
234: }

236: static PetscErrorCode PCSetFromOptions_KSP(PetscOptionItems *PetscOptionsObject,PC pc)
237: {
238:   PC_KSP         *jac = (PC_KSP*)pc->data;

242:   PetscOptionsHead(PetscOptionsObject,"PC KSP options");
243:   if (jac->ksp) {
244:     KSPSetFromOptions(jac->ksp);
245:    }
246:   PetscOptionsTail();
247:   return(0);
248: }

250: /* ----------------------------------------------------------------------------------*/

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

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

261:    Level: intermediate

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

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

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

280: M*/

282: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
283: {
285:   PC_KSP         *jac;

288:   PetscNewLog(pc,&jac);
289:   pc->data = (void*)jac;

291:   PetscMemzero(pc->ops,sizeof(struct _PCOps));
292:   pc->ops->apply           = PCApply_KSP;
293:   pc->ops->matapply        = PCMatApply_KSP;
294:   pc->ops->applytranspose  = PCApplyTranspose_KSP;
295:   pc->ops->setup           = PCSetUp_KSP;
296:   pc->ops->reset           = PCReset_KSP;
297:   pc->ops->destroy         = PCDestroy_KSP;
298:   pc->ops->setfromoptions  = PCSetFromOptions_KSP;
299:   pc->ops->view            = PCView_KSP;

301:   PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",PCKSPGetKSP_KSP);
302:   PetscObjectComposeFunction((PetscObject)pc,"PCKSPSetKSP_C",PCKSPSetKSP_KSP);
303:   return(0);
304: }