Actual source code: pcksp.c

petsc-3.6.0 2015-06-09
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: }

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

 38:   KSPSetInitialGuessNonzero(jac->ksp,pc->nonzero_guess);
 39:   KSPSolve(jac->ksp,x,y);
 40:   KSPGetIterationNumber(jac->ksp,&its);
 41:   jac->its += its;
 42:   return(0);
 43: }

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

 54:   KSPSolveTranspose(jac->ksp,x,y);
 55:   KSPGetIterationNumber(jac->ksp,&its);
 56:   jac->its += its;
 57:   return(0);
 58: }

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

 69:   if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
 70:   KSPSetFromOptions(jac->ksp);
 71:   if (pc->useAmat) mat = pc->mat;
 72:   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 */
 82: static PetscErrorCode PCReset_KSP(PC pc)
 83: {
 84:   PC_KSP         *jac = (PC_KSP*)pc->data;

 88:   if (jac->ksp) {KSPReset(jac->ksp);}
 89:   return(0);
 90: }

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

100:   PCReset_KSP(pc);
101:   KSPDestroy(&jac->ksp);
102:   PetscFree(pc->data);
103:   return(0);
104: }

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

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

135: static PetscErrorCode  PCKSPGetKSP_KSP(PC pc,KSP *ksp)
136: {
137:   PC_KSP         *jac = (PC_KSP*)pc->data;

141:   if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
142:   *ksp = jac->ksp;
143:   return(0);
144: }

148: /*@
149:    PCKSPGetKSP - Gets the KSP context for a KSP PC.

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

153:    Input Parameter:
154: .  pc - the preconditioner context

156:    Output Parameters:
157: .  ksp - the PC solver

159:    Notes:
160:    You must call KSPSetUp() before calling PCKSPGetKSP().

162:    Level: advanced

164: .keywords:  PC, KSP, get, context
165: @*/
166: PetscErrorCode  PCKSPGetKSP(PC pc,KSP *ksp)
167: {

173:   *ksp = NULL;
174:   PetscTryMethod(pc,"PCKSPGetKSP_C",(PC,KSP*),(pc,ksp));
175:   return(0);
176: }

178: /* ----------------------------------------------------------------------------------*/

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

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

189:    Level: intermediate

191:    Concepts: inner iteration

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

196:    Developer Notes: PCApply_KSP() uses the flag set by PCSetInitialGuessNonzero(), I think this is totally wrong, because it is then not
197:      using this inner KSP as a preconditioner (that is a linear operator applied to some vector), it is actually just using
198:      the inner KSP just like the outer KSP.

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

203: M*/

207: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
208: {
210:   PC_KSP         *jac;

213:   PetscNewLog(pc,&jac);

215:   pc->ops->apply           = PCApply_KSP;
216:   pc->ops->applytranspose  = PCApplyTranspose_KSP;
217:   pc->ops->setup           = PCSetUp_KSP;
218:   pc->ops->reset           = PCReset_KSP;
219:   pc->ops->destroy         = PCDestroy_KSP;
220:   pc->ops->setfromoptions  = 0;
221:   pc->ops->view            = PCView_KSP;
222:   pc->ops->applyrichardson = 0;

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


227:   jac->its             = 0;
228:   PetscObjectComposeFunction((PetscObject)pc,"PCKSPGetKSP_C",PCKSPGetKSP_KSP);
229:   return(0);
230: }