Actual source code: pcksp.c

      1: #define PETSCKSP_DLL

     3:  #include src/ksp/pc/pcimpl.h
     4:  #include petscksp.h

      6: typedef struct {
      7:   PetscTruth use_true_matrix;       /* use mat rather than pmat in inner linear solve */
      8:   KSP        ksp;
      9:   PetscInt   its;                   /* total number of iterations KSP uses */
     10: } PC_KSP;

     14: static PetscErrorCode PCApply_KSP(PC pc,Vec x,Vec y)
     15: {
     17:   PetscInt       its;
     18:   PC_KSP         *jac = (PC_KSP*)pc->data;

     21:   KSPSolve(jac->ksp,x,y);
     22:   KSPGetIterationNumber(jac->ksp,&its);
     23:   jac->its += its;
     24:   return(0);
     25: }

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

     36:   KSPSolveTranspose(jac->ksp,x,y);
     37:   KSPGetIterationNumber(jac->ksp,&its);
     38:   jac->its += its;
     39:   return(0);
     40: }

     44: static PetscErrorCode PCSetUp_KSP(PC pc)
     45: {
     47:   PC_KSP         *jac = (PC_KSP*)pc->data;
     48:   Mat            mat;

     51:   KSPSetFromOptions(jac->ksp);
     52:   if (jac->use_true_matrix) mat = pc->mat;
     53:   else                      mat = pc->pmat;

     55:   KSPSetOperators(jac->ksp,mat,pc->pmat,pc->flag);
     56:   KSPSetUp(jac->ksp);
     57:   return(0);
     58: }

     60: /* Default destroy, if it has never been setup */
     63: static PetscErrorCode PCDestroy_KSP(PC pc)
     64: {
     65:   PC_KSP         *jac = (PC_KSP*)pc->data;

     69:   KSPDestroy(jac->ksp);
     70:   PetscFree(jac);
     71:   return(0);
     72: }

     76: static PetscErrorCode PCView_KSP(PC pc,PetscViewer viewer)
     77: {
     78:   PC_KSP         *jac = (PC_KSP*)pc->data;
     80:   PetscTruth     iascii;

     83:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
     84:   if (iascii) {
     85:     if (jac->use_true_matrix) {
     86:       PetscViewerASCIIPrintf(viewer,"Using true matrix (not preconditioner matrix) on inner solve\n");
     87:     }
     88:     PetscViewerASCIIPrintf(viewer,"KSP and PC on KSP preconditioner follow\n");
     89:     PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
     90:   } else {
     91:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for this object",((PetscObject)viewer)->type_name);
     92:   }
     93:   PetscViewerASCIIPushTab(viewer);
     94:   KSPView(jac->ksp,viewer);
     95:   PetscViewerASCIIPopTab(viewer);
     96:   if (iascii) {
     97:     PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
     98:   }
     99:   return(0);
    100: }

    104: static PetscErrorCode PCSetFromOptions_KSP(PC pc)
    105: {
    107:   PetscTruth     flg;

    110:   PetscOptionsHead("KSP preconditioner options");
    111:     PetscOptionsName("-pc_ksp_true","Use true matrix to define inner linear system, not preconditioner matrix","PCKSPSetUseTrue",&flg);
    112:     if (flg) {
    113:       PCKSPSetUseTrue(pc);
    114:     }
    115:   PetscOptionsTail();
    116:   return(0);
    117: }

    119: /* ----------------------------------------------------------------------------------*/

    124: PetscErrorCode PETSCKSP_DLLEXPORT PCKSPSetUseTrue_KSP(PC pc)
    125: {
    126:   PC_KSP   *jac;

    129:   jac                  = (PC_KSP*)pc->data;
    130:   jac->use_true_matrix = PETSC_TRUE;
    131:   return(0);
    132: }

    138: PetscErrorCode PETSCKSP_DLLEXPORT PCKSPGetKSP_KSP(PC pc,KSP *ksp)
    139: {
    140:   PC_KSP   *jac;

    143:   jac          = (PC_KSP*)pc->data;
    144:   *ksp        = jac->ksp;
    145:   return(0);
    146: }

    151: /*@
    152:    PCKSPSetUseTrue - Sets a flag to indicate that the true matrix (rather than
    153:    the matrix used to define the preconditioner) is used to compute the
    154:    residual inside the inner solve.

    156:    Collective on PC

    158:    Input Parameters:
    159: .  pc - the preconditioner context

    161:    Options Database Key:
    162: .  -pc_ksp_true - Activates PCKSPSetUseTrue()

    164:    Note:
    165:    For the common case in which the preconditioning and linear
    166:    system matrices are identical, this routine is unnecessary.

    168:    Level: advanced

    170: .keywords:  PC, KSP, set, true, local, flag

    172: .seealso: PCSetOperators(), PCBJacobiSetUseTrueLocal()
    173: @*/
    174: PetscErrorCode PETSCKSP_DLLEXPORT PCKSPSetUseTrue(PC pc)
    175: {
    176:   PetscErrorCode ierr,(*f)(PC);

    180:   PetscObjectQueryFunction((PetscObject)pc,"PCKSPSetUseTrue_C",(void (**)(void))&f);
    181:   if (f) {
    182:     (*f)(pc);
    183:   }
    184:   return(0);
    185: }

    189: /*@
    190:    PCKSPGetKSP - Gets the KSP context for a KSP PC.

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

    194:    Input Parameter:
    195: .  pc - the preconditioner context

    197:    Output Parameters:
    198: .  ksp - the PC solver

    200:    Notes:
    201:    You must call KSPSetUp() before calling PCKSPGetKSP().

    203:    Level: advanced

    205: .keywords:  PC, KSP, get, context
    206: @*/
    207: PetscErrorCode PETSCKSP_DLLEXPORT PCKSPGetKSP(PC pc,KSP *ksp)
    208: {
    209:   PetscErrorCode ierr,(*f)(PC,KSP*);

    214:   if (!pc->setupcalled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must call KSPSetUp first");
    215:   PetscObjectQueryFunction((PetscObject)pc,"PCKSPGetKSP_C",(void (**)(void))&f);
    216:   if (f) {
    217:     (*f)(pc,ksp);
    218:   }
    219:   return(0);
    220: }

    222: /* ----------------------------------------------------------------------------------*/

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

    228:    Options Database Key:
    229: .     -pc_ksp_true - use the matrix that defines the linear system as the matrix for the
    230:                     inner solver, otherwise by default it uses the matrix used to construct
    231:                     the preconditioner (see PCSetOperators())

    233:    Level: intermediate

    235:    Concepts: inner iteration

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

    241: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
    242:            PCSHELL, PCCOMPOSITE, PCKSPUseTrue(), PCKSPGetKSP()

    244: M*/

    249: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_KSP(PC pc)
    250: {
    252:   const char     *prefix;
    253:   PC_KSP         *jac;

    256:   PetscNew(PC_KSP,&jac);
    257:   PetscLogObjectMemory(pc,sizeof(PC_KSP));
    258:   pc->ops->apply              = PCApply_KSP;
    259:   pc->ops->applytranspose     = PCApplyTranspose_KSP;
    260:   pc->ops->setup              = PCSetUp_KSP;
    261:   pc->ops->destroy            = PCDestroy_KSP;
    262:   pc->ops->setfromoptions     = PCSetFromOptions_KSP;
    263:   pc->ops->view               = PCView_KSP;
    264:   pc->ops->applyrichardson    = 0;

    266:   pc->data               = (void*)jac;
    267:   KSPCreate(pc->comm,&jac->ksp);

    269:   PCGetOptionsPrefix(pc,&prefix);
    270:   KSPSetOptionsPrefix(jac->ksp,prefix);
    271:   KSPAppendOptionsPrefix(jac->ksp,"ksp_");
    272:   jac->use_true_matrix = PETSC_FALSE;
    273:   jac->its             = 0;

    275:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCKSPSetUseTrue_C","PCKSPSetUseTrue_KSP",
    276:                     PCKSPSetUseTrue_KSP);
    277:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCKSPGetKSP_C","PCKSPGetKSP_KSP",
    278:                     PCKSPGetKSP_KSP);

    280:   return(0);
    281: }