Actual source code: galerkin.c

      2: /*
      3:       Defines a preconditioner defined by R^T S R
      4: */
      5: #include <petsc/private/pcimpl.h>
      6: #include <petscksp.h>

      8: typedef struct {
      9:   KSP            ksp;
     10:   Mat            R,P;
     11:   Vec            b,x;
     12:   PetscErrorCode (*computeasub)(PC,Mat,Mat,Mat*,void*);
     13:   void           *computeasub_ctx;
     14: } PC_Galerkin;

     16: static PetscErrorCode PCApply_Galerkin(PC pc,Vec x,Vec y)
     17: {
     19:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;

     22:   if (jac->R) {
     23:     MatRestrict(jac->R,x,jac->b);
     24:   } else {
     25:     MatRestrict(jac->P,x,jac->b);
     26:   }
     27:   KSPSolve(jac->ksp,jac->b,jac->x);
     28:   KSPCheckSolve(jac->ksp,pc,jac->x);
     29:   if (jac->P) {
     30:     MatInterpolate(jac->P,jac->x,y);
     31:   } else {
     32:     MatInterpolate(jac->R,jac->x,y);
     33:   }
     34:   return(0);
     35: }

     37: static PetscErrorCode PCSetUp_Galerkin(PC pc)
     38: {
     40:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;
     41:   PetscBool      a;
     42:   Vec            *xx,*yy;

     45:   if (jac->computeasub) {
     46:     Mat Ap;
     47:     if (!pc->setupcalled) {
     48:       (*jac->computeasub)(pc,pc->pmat,NULL,&Ap,jac->computeasub_ctx);
     49:       KSPSetOperators(jac->ksp,Ap,Ap);
     50:       MatDestroy(&Ap);
     51:     } else {
     52:       KSPGetOperators(jac->ksp,NULL,&Ap);
     53:       (*jac->computeasub)(pc,pc->pmat,Ap,NULL,jac->computeasub_ctx);
     54:     }
     55:   }

     57:   if (!jac->x) {
     58:     KSPGetOperatorsSet(jac->ksp,&a,NULL);
     59:     if (!a) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set operator of PCGALERKIN KSP with PCGalerkinGetKSP()/KSPSetOperators()");
     60:     KSPCreateVecs(jac->ksp,1,&xx,1,&yy);
     61:     jac->x = *xx;
     62:     jac->b = *yy;
     63:     PetscFree(xx);
     64:     PetscFree(yy);
     65:   }
     66:   if (!jac->R && !jac->P) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set restriction or interpolation of PCGALERKIN with PCGalerkinSetRestriction()/Interpolation()");
     67:   /* should check here that sizes of R/P match size of a */

     69:   return(0);
     70: }

     72: static PetscErrorCode PCReset_Galerkin(PC pc)
     73: {
     74:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;

     78:   MatDestroy(&jac->R);
     79:   MatDestroy(&jac->P);
     80:   VecDestroy(&jac->x);
     81:   VecDestroy(&jac->b);
     82:   KSPReset(jac->ksp);
     83:   return(0);
     84: }

     86: static PetscErrorCode PCDestroy_Galerkin(PC pc)
     87: {
     88:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;

     92:   PCReset_Galerkin(pc);
     93:   KSPDestroy(&jac->ksp);
     94:   PetscFree(pc->data);
     95:   return(0);
     96: }

     98: static PetscErrorCode PCView_Galerkin(PC pc,PetscViewer viewer)
     99: {
    100:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;
    102:   PetscBool      iascii;

    105:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
    106:   if (iascii) {
    107:     PetscViewerASCIIPrintf(viewer,"  KSP on Galerkin follow\n");
    108:     PetscViewerASCIIPrintf(viewer,"  ---------------------------------\n");
    109:   }
    110:   KSPView(jac->ksp,viewer);
    111:   return(0);
    112: }

    114: static PetscErrorCode  PCGalerkinGetKSP_Galerkin(PC pc,KSP *ksp)
    115: {
    116:   PC_Galerkin *jac = (PC_Galerkin*)pc->data;

    119:   *ksp = jac->ksp;
    120:   return(0);
    121: }

    123: static PetscErrorCode  PCGalerkinSetRestriction_Galerkin(PC pc,Mat R)
    124: {
    125:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;

    129:   PetscObjectReference((PetscObject)R);
    130:   MatDestroy(&jac->R);
    131:   jac->R = R;
    132:   return(0);
    133: }

    135: static PetscErrorCode  PCGalerkinSetInterpolation_Galerkin(PC pc,Mat P)
    136: {
    137:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;

    141:   PetscObjectReference((PetscObject)P);
    142:   MatDestroy(&jac->P);
    143:   jac->P = P;
    144:   return(0);
    145: }

    147: static PetscErrorCode  PCGalerkinSetComputeSubmatrix_Galerkin(PC pc,PetscErrorCode (*computeAsub)(PC,Mat,Mat,Mat*,void*),void *ctx)
    148: {
    149:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;

    152:   jac->computeasub     = computeAsub;
    153:   jac->computeasub_ctx = ctx;
    154:   return(0);
    155: }

    157: /* -------------------------------------------------------------------------------- */
    158: /*@
    159:    PCGalerkinSetRestriction - Sets the restriction operator for the "Galerkin-type" preconditioner

    161:    Logically Collective on PC

    163:    Input Parameter:
    164: +  pc - the preconditioner context
    165: -  R - the restriction operator

    167:    Notes:
    168:     Either this or PCGalerkinSetInterpolation() or both must be called

    170:    Level: Intermediate

    172: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
    173:            PCGalerkinSetInterpolation(), PCGalerkinGetKSP()

    175: @*/
    176: PetscErrorCode  PCGalerkinSetRestriction(PC pc,Mat R)
    177: {

    182:   PetscTryMethod(pc,"PCGalerkinSetRestriction_C",(PC,Mat),(pc,R));
    183:   return(0);
    184: }

    186: /*@
    187:    PCGalerkinSetInterpolation - Sets the interpolation operator for the "Galerkin-type" preconditioner

    189:    Logically Collective on PC

    191:    Input Parameter:
    192: +  pc - the preconditioner context
    193: -  R - the interpolation operator

    195:    Notes:
    196:     Either this or PCGalerkinSetRestriction() or both must be called

    198:    Level: Intermediate

    200: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
    201:            PCGalerkinSetRestriction(), PCGalerkinGetKSP()

    203: @*/
    204: PetscErrorCode  PCGalerkinSetInterpolation(PC pc,Mat P)
    205: {

    210:   PetscTryMethod(pc,"PCGalerkinSetInterpolation_C",(PC,Mat),(pc,P));
    211:   return(0);
    212: }

    214: /*@
    215:    PCGalerkinSetComputeSubmatrix - Provide a routine that will be called to compute the Galerkin submatrix

    217:    Logically Collective

    219:    Input Parameter:
    220: +  pc - the preconditioner context
    221: .  computeAsub - routine that computes the submatrix from the global matrix
    222: -  ctx - context used by the routine, or NULL

    224:    Calling sequence of computeAsub:
    225: $    computeAsub(PC pc,Mat A, Mat Ap, Mat *cAP,void *ctx);

    227: +  PC - the Galerkin PC
    228: .  A - the matrix in the Galerkin PC
    229: .  Ap - the computed submatrix from any previous computation, if NULL it has not previously been computed
    230: .  cAp - the submatrix computed by this routine
    231: -  ctx - optional user-defined function context

    233:    Level: Intermediate

    235:    Notes:
    236:     Instead of providing this routine you can call PCGalerkinGetKSP() and then KSPSetOperators() to provide the submatrix,
    237:           but that will not work for multiple KSPSolves with different matrices unless you call it for each solve.

    239:           This routine is called each time the outer matrix is changed. In the first call the Ap argument is NULL and the routine should create the
    240:           matrix and computes its values in cAp. On each subsequent call the routine should up the Ap matrix.

    242:    Developer Notes:
    243:     If the user does not call this routine nor call PCGalerkinGetKSP() and KSPSetOperators() then PCGalerkin could
    244:                     could automatically compute the submatrix via calls to MatGalerkin() or MatRARt()

    246: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
    247:            PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinGetKSP()

    249: @*/
    250: PetscErrorCode  PCGalerkinSetComputeSubmatrix(PC pc,PetscErrorCode (*computeAsub)(PC,Mat,Mat,Mat*,void*),void *ctx)
    251: {

    256:   PetscTryMethod(pc,"PCGalerkinSetComputeSubmatrix_C",(PC,PetscErrorCode (*)(PC,Mat,Mat,Mat*,void*),void*),(pc,computeAsub,ctx));
    257:   return(0);
    258: }

    260: /*@
    261:    PCGalerkinGetKSP - Gets the KSP object in the Galerkin PC.

    263:    Not Collective

    265:    Input Parameter:
    266: .  pc - the preconditioner context

    268:    Output Parameters:
    269: .  ksp - the KSP object

    271:    Level: Intermediate

    273:    Notes:
    274:     Once you have called this routine you can call KSPSetOperators() on the resulting ksp to provide the operator for the Galerkin problem,
    275:           an alternative is to use PCGalerkinSetComputeSubmatrix() to provide a routine that computes the submatrix as needed.

    277: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PCGALERKIN,
    278:            PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinSetComputeSubmatrix()

    280: @*/
    281: PetscErrorCode  PCGalerkinGetKSP(PC pc,KSP *ksp)
    282: {

    288:   PetscUseMethod(pc,"PCGalerkinGetKSP_C",(PC,KSP*),(pc,ksp));
    289:   return(0);
    290: }

    292: static PetscErrorCode PCSetFromOptions_Galerkin(PetscOptionItems *PetscOptionsObject,PC pc)
    293: {
    294:   PC_Galerkin    *jac = (PC_Galerkin*)pc->data;
    296:   const char     *prefix;
    297:   PetscBool      flg;

    300:   KSPGetOptionsPrefix(jac->ksp,&prefix);
    301:   PetscStrendswith(prefix,"galerkin_",&flg);
    302:   if (!flg) {
    303:     PCGetOptionsPrefix(pc,&prefix);
    304:     KSPSetOptionsPrefix(jac->ksp,prefix);
    305:     KSPAppendOptionsPrefix(jac->ksp,"galerkin_");
    306:   }

    308:   PetscOptionsHead(PetscOptionsObject,"Galerkin options");
    309:   if (jac->ksp) {
    310:     KSPSetFromOptions(jac->ksp);
    311:   }
    312:   PetscOptionsTail();
    313:   return(0);
    314: }

    316: /* -------------------------------------------------------------------------------------------*/

    318: /*MC
    319:      PCGALERKIN - Build (part of) a preconditioner by P S R (where P is often R^T)

    321: $   Use PCGalerkinSetRestriction(pc,R) and/or PCGalerkinSetInterpolation(pc,P) followed by
    322: $   PCGalerkinGetKSP(pc,&ksp); KSPSetOperators(ksp,A,....)

    324:    Level: intermediate

    326:    Developer Note: If KSPSetOperators() has not been called on the inner KSP then PCGALERKIN could use MatRARt() or MatPtAP() to compute
    327:                    the operators automatically.
    328:                    Should there be a prefix for the inner KSP.
    329:                    There is no KSPSetFromOptions_Galerkin() that calls KSPSetFromOptions() on the inner KSP

    331: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
    332:            PCSHELL, PCKSP, PCGalerkinSetRestriction(), PCGalerkinSetInterpolation(), PCGalerkinGetKSP()

    334: M*/

    336: PETSC_EXTERN PetscErrorCode PCCreate_Galerkin(PC pc)
    337: {
    339:   PC_Galerkin    *jac;

    342:   PetscNewLog(pc,&jac);

    344:   pc->ops->apply           = PCApply_Galerkin;
    345:   pc->ops->setup           = PCSetUp_Galerkin;
    346:   pc->ops->reset           = PCReset_Galerkin;
    347:   pc->ops->destroy         = PCDestroy_Galerkin;
    348:   pc->ops->view            = PCView_Galerkin;
    349:   pc->ops->setfromoptions  = PCSetFromOptions_Galerkin;
    350:   pc->ops->applyrichardson = NULL;

    352:   KSPCreate(PetscObjectComm((PetscObject)pc),&jac->ksp);
    353:   KSPSetErrorIfNotConverged(jac->ksp,pc->erroriffailure);
    354:   PetscObjectIncrementTabLevel((PetscObject)jac->ksp,(PetscObject)pc,1);

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

    358:   PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetRestriction_C",PCGalerkinSetRestriction_Galerkin);
    359:   PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetInterpolation_C",PCGalerkinSetInterpolation_Galerkin);
    360:   PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinGetKSP_C",PCGalerkinGetKSP_Galerkin);
    361:   PetscObjectComposeFunction((PetscObject)pc,"PCGalerkinSetComputeSubmatrix_C",PCGalerkinSetComputeSubmatrix_Galerkin);
    362:   return(0);
    363: }