Actual source code: pcksp.c

petsc-3.5.2 2014-09-08
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:   PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);
 22:   PCGetOptionsPrefix(pc,&prefix);
 23:   KSPSetOptionsPrefix(jac->ksp,prefix);
 24:   KSPAppendOptionsPrefix(jac->ksp,"ksp_");
 25:   return(0);
 26: }

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

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

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

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

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

 68:   if (!jac->ksp) {PCKSPCreateKSP_KSP(pc);}
 69:   KSPSetFromOptions(jac->ksp);
 70:   if (pc->useAmat) mat = pc->mat;
 71:   else             mat = pc->pmat;

 73:   KSPSetOperators(jac->ksp,mat,pc->pmat);
 74:   KSPSetUp(jac->ksp);
 75:   return(0);
 76: }

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

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

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

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

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

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

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

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

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

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

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

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

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

161:    Level: advanced

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

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

177: /* ----------------------------------------------------------------------------------*/

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

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

188:    Level: intermediate

190:    Concepts: inner iteration

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

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

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

202: M*/

206: PETSC_EXTERN PetscErrorCode PCCreate_KSP(PC pc)
207: {
209:   PC_KSP         *jac;

212:   PetscNewLog(pc,&jac);

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

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


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