Actual source code: pcksp.c

petsc-3.7.4 2016-10-02
Report Typos and Errors
  2: #include <petsc/private/pcimpl.h>
  3: #include <petscksp.h>            /*I "petscksp.h" I*/

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


 13: static PetscErrorCode  PCKSPCreateKSP_KSP(PC pc)
 14: {
 16:   const char     *prefix;
 17:   PC_KSP         *jac = (PC_KSP*)pc->data;

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

 29: #include <petsc/private/kspimpl.h>
 32: static PetscErrorCode PCApply_KSP(PC pc,Vec x,Vec y)
 33: {
 35:   PetscInt       its;
 36:   PC_KSP         *jac = (PC_KSP*)pc->data;

 39:   KSPSolve(jac->ksp,x,y);
 40:   KSPGetIterationNumber(jac->ksp,&its);
 41:   jac->its += its;
 42:   if (jac->ksp->reason == KSP_DIVERGED_PCSETUP_FAILED) {
 43:     pc->failedreason = PC_SUBPC_ERROR;
 44:   }
 45:   return(0);
 46: }

 50: static PetscErrorCode PCApplyTranspose_KSP(PC pc,Vec x,Vec y)
 51: {
 53:   PetscInt       its;
 54:   PC_KSP         *jac = (PC_KSP*)pc->data;

 57:   KSPSolveTranspose(jac->ksp,x,y);
 58:   KSPGetIterationNumber(jac->ksp,&its);
 59:   jac->its += its;
 60:   return(0);
 61: }

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

 72:   if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
 73:   KSPSetFromOptions(jac->ksp);
 74:   if (pc->useAmat) mat = pc->mat;
 75:   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 */
 85: static PetscErrorCode PCReset_KSP(PC pc)
 86: {
 87:   PC_KSP         *jac = (PC_KSP*)pc->data;

 91:   if (jac->ksp) {KSPReset(jac->ksp);}
 92:   return(0);
 93: }

 97: static PetscErrorCode PCDestroy_KSP(PC pc)
 98: {
 99:   PC_KSP         *jac = (PC_KSP*)pc->data;

103:   PCReset_KSP(pc);
104:   KSPDestroy(&jac->ksp);
105:   PetscFree(pc->data);
106:   return(0);
107: }

111: static PetscErrorCode PCView_KSP(PC pc,PetscViewer viewer)
112: {
113:   PC_KSP         *jac = (PC_KSP*)pc->data;
115:   PetscBool      iascii;

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

138: static PetscErrorCode  PCKSPGetKSP_KSP(PC pc,KSP *ksp)
139: {
140:   PC_KSP         *jac = (PC_KSP*)pc->data;

144:   if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
145:   *ksp = jac->ksp;
146:   return(0);
147: }

151: /*@
152:    PCKSPGetKSP - Gets the KSP context for a KSP PC.

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

156:    Input Parameter:
157: .  pc - the preconditioner context

159:    Output Parameters:
160: .  ksp - the PC solver

162:    Notes:
163:    You must call KSPSetUp() before calling PCKSPGetKSP().

165:    If the PC is not a PCKSP object then a NULL is returned

167:    Level: advanced

169: .keywords:  PC, KSP, get, context
170: @*/
171: PetscErrorCode  PCKSPGetKSP(PC pc,KSP *ksp)
172: {

178:   *ksp = NULL;
179:   PetscTryMethod(pc,"PCKSPGetKSP_C",(PC,KSP*),(pc,ksp));
180:   return(0);
181: }

183: /* ----------------------------------------------------------------------------------*/

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

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

194:    Level: intermediate

196:    Concepts: inner iteration

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

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

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

213: M*/

217: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
218: {
220:   PC_KSP         *jac;

223:   PetscNewLog(pc,&jac);

225:   pc->ops->apply           = PCApply_KSP;
226:   pc->ops->applytranspose  = PCApplyTranspose_KSP;
227:   pc->ops->setup           = PCSetUp_KSP;
228:   pc->ops->reset           = PCReset_KSP;
229:   pc->ops->destroy         = PCDestroy_KSP;
230:   pc->ops->setfromoptions  = 0;
231:   pc->ops->view            = PCView_KSP;
232:   pc->ops->applyrichardson = 0;

234:   pc->data = (void*)jac;


237:   jac->its             = 0;
238:   PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",PCKSPGetKSP_KSP);
239:   return(0);
240: }