Actual source code: composite.c

      1: #define PETSCKSP_DLL

      3: /*
      4:       Defines a preconditioner that can consist of a collection of PCs
      5: */
     6:  #include src/ksp/pc/pcimpl.h
     7:  #include petscksp.h

      9: typedef struct _PC_CompositeLink *PC_CompositeLink;
     10: struct _PC_CompositeLink {
     11:   PC               pc;
     12:   PC_CompositeLink next;
     13: };
     14:
     15: typedef struct {
     16:   PC_CompositeLink head;
     17:   PCCompositeType  type;
     18:   Vec              work1;
     19:   Vec              work2;
     20:   PetscScalar      alpha;
     21:   PetscTruth       use_true_matrix;
     22: } PC_Composite;

     26: static PetscErrorCode PCApply_Composite_Multiplicative(PC pc,Vec x,Vec y)
     27: {
     28:   PetscErrorCode   ierr;
     29:   PC_Composite     *jac = (PC_Composite*)pc->data;
     30:   PC_CompositeLink next = jac->head;
     31:   PetscScalar      one = 1.0,mone = -1.0;
     32:   Mat              mat = pc->pmat;

     35:   if (!next) {
     36:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC()");
     37:   }
     38:   if (next->next && !jac->work2) { /* allocate second work vector */
     39:     VecDuplicate(jac->work1,&jac->work2);
     40:   }
     41:   PCApply(next->pc,x,y);
     42:   if (jac->use_true_matrix) mat = pc->mat;
     43:   while (next->next) {
     44:     next = next->next;
     45:     MatMult(mat,y,jac->work1);
     46:     VecWAXPY(jac->work2,mone,jac->work1,x);
     47:     PCApply(next->pc,jac->work2,jac->work1);
     48:     VecAXPY(y,one,jac->work1);
     49:   }

     51:   return(0);
     52: }

     54: /*
     55:     This is very special for a matrix of the form alpha I + R + S
     56: where first preconditioner is built from alpha I + S and second from
     57: alpha I + R
     58: */
     61: static PetscErrorCode PCApply_Composite_Special(PC pc,Vec x,Vec y)
     62: {
     63:   PetscErrorCode   ierr;
     64:   PC_Composite     *jac = (PC_Composite*)pc->data;
     65:   PC_CompositeLink next = jac->head;

     68:   if (!next) {
     69:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC()");
     70:   }
     71:   if (!next->next || next->next->next) {
     72:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Special composite preconditioners requires exactly two PCs");
     73:   }

     75:   PCApply(next->pc,x,jac->work1);
     76:   PCApply(next->next->pc,jac->work1,y);
     77:   return(0);
     78: }

     82: static PetscErrorCode PCApply_Composite_Additive(PC pc,Vec x,Vec y)
     83: {
     84:   PetscErrorCode   ierr;
     85:   PC_Composite     *jac = (PC_Composite*)pc->data;
     86:   PC_CompositeLink next = jac->head;
     87:   PetscScalar      one = 1.0;

     90:   if (!next) {
     91:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"No composite preconditioners supplied via PCCompositeAddPC()");
     92:   }
     93:   PCApply(next->pc,x,y);
     94:   while (next->next) {
     95:     next = next->next;
     96:     PCApply(next->pc,x,jac->work1);
     97:     VecAXPY(y,one,jac->work1);
     98:   }
     99:   return(0);
    100: }

    104: static PetscErrorCode PCSetUp_Composite(PC pc)
    105: {
    106:   PetscErrorCode   ierr;
    107:   PC_Composite     *jac = (PC_Composite*)pc->data;
    108:   PC_CompositeLink next = jac->head;

    111:   if (!jac->work1) {
    112:    MatGetVecs(pc->pmat,&jac->work1,0);
    113:   }
    114:   while (next) {
    115:     PCSetOperators(next->pc,pc->mat,pc->pmat,pc->flag);
    116:     next = next->next;
    117:   }

    119:   return(0);
    120: }

    124: static PetscErrorCode PCDestroy_Composite(PC pc)
    125: {
    126:   PC_Composite     *jac = (PC_Composite*)pc->data;
    127:   PetscErrorCode   ierr;
    128:   PC_CompositeLink next = jac->head;

    131:   while (next) {
    132:     PCDestroy(next->pc);
    133:     next = next->next;
    134:   }

    136:   if (jac->work1) {VecDestroy(jac->work1);}
    137:   if (jac->work2) {VecDestroy(jac->work2);}
    138:   PetscFree(jac);
    139:   return(0);
    140: }

    144: static PetscErrorCode PCSetFromOptions_Composite(PC pc)
    145: {
    146:   PC_Composite     *jac = (PC_Composite*)pc->data;
    147:   PetscErrorCode   ierr;
    148:   PetscInt         nmax = 8,i;
    149:   PC_CompositeLink next;
    150:   char             *pcs[8];
    151:   PetscTruth       flg;

    154:   PetscOptionsHead("Composite preconditioner options");
    155:   PetscOptionsEnum("-pc_composite_type","Type of composition","PCCompositeSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
    156:     PetscOptionsName("-pc_composite_true","Use true matrix for inner solves","PCCompositeSetUseTrue",&flg);
    157:     if (flg) {
    158:       PCCompositeSetUseTrue(pc);
    159:     }
    160:     PetscOptionsStringArray("-pc_composite_pcs","List of composite solvers","PCCompositeAddPC",pcs,&nmax,&flg);
    161:     if (flg) {
    162:       for (i=0; i<nmax; i++) {
    163:         PCCompositeAddPC(pc,pcs[i]);
    164:       }
    165:     }
    166:   PetscOptionsTail();

    168:   next = jac->head;
    169:   while (next) {
    170:     PCSetFromOptions(next->pc);
    171:     next = next->next;
    172:   }
    173:   return(0);
    174: }

    178: static PetscErrorCode PCView_Composite(PC pc,PetscViewer viewer)
    179: {
    180:   PC_Composite     *jac = (PC_Composite*)pc->data;
    181:   PetscErrorCode   ierr;
    182:   PC_CompositeLink next = jac->head;
    183:   PetscTruth       iascii;

    186:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    187:   if (iascii) {
    188:     PetscViewerASCIIPrintf(viewer,"Composite PC type - %s\n",PCCompositeTypes[jac->type]);
    189:     PetscViewerASCIIPrintf(viewer,"PCs on composite preconditioner follow\n");
    190:     PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
    191:   } else {
    192:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCComposite",((PetscObject)viewer)->type_name);
    193:   }
    194:   if (iascii) {
    195:     PetscViewerASCIIPushTab(viewer);
    196:   }
    197:   while (next) {
    198:     PCView(next->pc,viewer);
    199:     next = next->next;
    200:   }
    201:   if (iascii) {
    202:     PetscViewerASCIIPopTab(viewer);
    203:     PetscViewerASCIIPrintf(viewer,"---------------------------------\n");
    204:   }
    205:   return(0);
    206: }

    208: /* ------------------------------------------------------------------------------*/

    213: PetscErrorCode PETSCKSP_DLLEXPORT PCCompositeSpecialSetAlpha_Composite(PC pc,PetscScalar alpha)
    214: {
    215:   PC_Composite *jac = (PC_Composite*)pc->data;
    217:   jac->alpha = alpha;
    218:   return(0);
    219: }

    225: PetscErrorCode PETSCKSP_DLLEXPORT PCCompositeSetType_Composite(PC pc,PCCompositeType type)
    226: {
    228:   if (type == PC_COMPOSITE_ADDITIVE) {
    229:     pc->ops->apply = PCApply_Composite_Additive;
    230:   } else if (type ==  PC_COMPOSITE_MULTIPLICATIVE) {
    231:     pc->ops->apply = PCApply_Composite_Multiplicative;
    232:   } else if (type ==  PC_COMPOSITE_SPECIAL) {
    233:     pc->ops->apply = PCApply_Composite_Special;
    234:   } else {
    235:     SETERRQ(PETSC_ERR_ARG_WRONG,"Unkown composite preconditioner type");
    236:   }
    237:   return(0);
    238: }

    244: PetscErrorCode PETSCKSP_DLLEXPORT PCCompositeAddPC_Composite(PC pc,PCType type)
    245: {
    246:   PC_Composite     *jac;
    247:   PC_CompositeLink next,ilink;
    248:   PetscErrorCode   ierr;
    249:   PetscInt         cnt = 0;
    250:   const char       *prefix;
    251:   char             newprefix[8];

    254:   PetscNew(struct _PC_CompositeLink,&ilink);
    255:   ilink->next = 0;
    256:   PCCreate(pc->comm,&ilink->pc);

    258:   jac  = (PC_Composite*)pc->data;
    259:   next = jac->head;
    260:   if (!next) {
    261:     jac->head = ilink;
    262:   } else {
    263:     cnt++;
    264:     while (next->next) {
    265:       next = next->next;
    266:       cnt++;
    267:     }
    268:     next->next = ilink;
    269:   }
    270:   PCGetOptionsPrefix(pc,&prefix);
    271:   PCSetOptionsPrefix(ilink->pc,prefix);
    272:   sprintf(newprefix,"sub_%d_",(int)cnt);
    273:   PCAppendOptionsPrefix(ilink->pc,newprefix);
    274:   /* type is set after prefix, because some methods may modify prefix, e.g. pcksp */
    275:   PCSetType(ilink->pc,type);

    277:   return(0);
    278: }

    284: PetscErrorCode PETSCKSP_DLLEXPORT PCCompositeGetPC_Composite(PC pc,PetscInt n,PC *subpc)
    285: {
    286:   PC_Composite     *jac;
    287:   PC_CompositeLink next;
    288:   PetscInt         i;

    291:   jac  = (PC_Composite*)pc->data;
    292:   next = jac->head;
    293:   for (i=0; i<n; i++) {
    294:     if (!next->next) {
    295:       SETERRQ(PETSC_ERR_ARG_INCOMP,"Not enough PCs in composite preconditioner");
    296:     }
    297:     next = next->next;
    298:   }
    299:   *subpc = next->pc;
    300:   return(0);
    301: }

    307: PetscErrorCode PETSCKSP_DLLEXPORT PCCompositeSetUseTrue_Composite(PC pc)
    308: {
    309:   PC_Composite   *jac;

    312:   jac                  = (PC_Composite*)pc->data;
    313:   jac->use_true_matrix = PETSC_TRUE;
    314:   return(0);
    315: }

    318: /* -------------------------------------------------------------------------------- */
    321: /*@C
    322:    PCCompositeSetType - Sets the type of composite preconditioner.
    323:
    324:    Collective on PC

    326:    Input Parameter:
    327: .  pc - the preconditioner context
    328: .  type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE, PC_COMPOSITE_SPECIAL

    330:    Options Database Key:
    331: .  -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type

    333:    Level: Developer

    335: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
    336: @*/
    337: PetscErrorCode PETSCKSP_DLLEXPORT PCCompositeSetType(PC pc,PCCompositeType type)
    338: {
    339:   PetscErrorCode ierr,(*f)(PC,PCCompositeType);

    343:   PetscObjectQueryFunction((PetscObject)pc,"PCCompositeSetType_C",(void (**)(void))&f);
    344:   if (f) {
    345:     (*f)(pc,type);
    346:   }
    347:   return(0);
    348: }

    352: /*@C
    353:    PCCompositeSpecialSetAlpha - Sets alpha for the special composite preconditioner
    354:      for alphaI + R + S
    355:
    356:    Collective on PC

    358:    Input Parameter:
    359: +  pc - the preconditioner context
    360: -  alpha - scale on identity

    362:    Level: Developer

    364: .keywords: PC, set, type, composite preconditioner, additive, multiplicative
    365: @*/
    366: PetscErrorCode PETSCKSP_DLLEXPORT PCCompositeSpecialSetAlpha(PC pc,PetscScalar alpha)
    367: {
    368:   PetscErrorCode ierr,(*f)(PC,PetscScalar);

    372:   PetscObjectQueryFunction((PetscObject)pc,"PCCompositeSpecialSetAlpha_C",(void (**)(void))&f);
    373:   if (f) {
    374:     (*f)(pc,alpha);
    375:   }
    376:   return(0);
    377: }

    381: /*@C
    382:    PCCompositeAddPC - Adds another PC to the composite PC.
    383:
    384:    Collective on PC

    386:    Input Parameters:
    387: .  pc - the preconditioner context
    388: .  type - the type of the new preconditioner

    390:    Level: Developer

    392: .keywords: PC, composite preconditioner, add
    393: @*/
    394: PetscErrorCode PETSCKSP_DLLEXPORT PCCompositeAddPC(PC pc,PCType type)
    395: {
    396:   PetscErrorCode ierr,(*f)(PC,PCType);

    400:   PetscObjectQueryFunction((PetscObject)pc,"PCCompositeAddPC_C",(void (**)(void))&f);
    401:   if (f) {
    402:     (*f)(pc,type);
    403:   }
    404:   return(0);
    405: }

    409: /*@C
    410:    PCCompositeGetPC - Gets one of the PC objects in the composite PC.
    411:
    412:    Not Collective

    414:    Input Parameter:
    415: .  pc - the preconditioner context
    416: .  n - the number of the pc requested

    418:    Output Parameters:
    419: .  subpc - the PC requested

    421:    Level: Developer

    423: .keywords: PC, get, composite preconditioner, sub preconditioner

    425: .seealso: PCCompositeAddPC()
    426: @*/
    427: PetscErrorCode PETSCKSP_DLLEXPORT PCCompositeGetPC(PC pc,PetscInt n,PC *subpc)
    428: {
    429:   PetscErrorCode ierr,(*f)(PC,PetscInt,PC *);

    434:   PetscObjectQueryFunction((PetscObject)pc,"PCCompositeGetPC_C",(void (**)(void))&f);
    435:   if (f) {
    436:     (*f)(pc,n,subpc);
    437:   } else {
    438:     SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot get pc, not composite type");
    439:   }
    440:   return(0);
    441: }

    445: /*@
    446:    PCCompositeSetUseTrue - Sets a flag to indicate that the true matrix (rather than
    447:                       the matrix used to define the preconditioner) is used to compute
    448:                       the residual when the multiplicative scheme is used.

    450:    Collective on PC

    452:    Input Parameters:
    453: .  pc - the preconditioner context

    455:    Options Database Key:
    456: .  -pc_composite_true - Activates PCCompositeSetUseTrue()

    458:    Note:
    459:    For the common case in which the preconditioning and linear
    460:    system matrices are identical, this routine is unnecessary.

    462:    Level: Developer

    464: .keywords: PC, composite preconditioner, set, true, flag

    466: .seealso: PCSetOperators(), PCBJacobiSetUseTrueLocal(), PCKSPSetUseTrue()
    467: @*/
    468: PetscErrorCode PETSCKSP_DLLEXPORT PCCompositeSetUseTrue(PC pc)
    469: {
    470:   PetscErrorCode ierr,(*f)(PC);

    474:   PetscObjectQueryFunction((PetscObject)pc,"PCCompositeSetUseTrue_C",(void (**)(void))&f);
    475:   if (f) {
    476:     (*f)(pc);
    477:   }
    478:   return(0);
    479: }

    481: /* -------------------------------------------------------------------------------------------*/

    483: /*MC
    484:      PCCOMPOSITE - Build a preconditioner by composing together several preconditioners

    486:    Options Database Keys:
    487: .  -pc_composite_type <type: one of multiplicative, additive, special> - Sets composite preconditioner type
    488: .  -pc_composite_true - Activates PCCompositeSetUseTrue()

    490:    Level: intermediate

    492:    Concepts: composing solvers

    494:    Notes: To use a Krylov method inside the composite preconditioner, set the PCType of one or more
    495:           inner PCs to be PCKSP.
    496:           Using a Krylov method inside another Krylov method can be dangerous (you get divergence or
    497:           the incorrect answer) unless you use KSPFGMRES as the outter Krylov method

    500: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
    501:            PCSHELL, PCKSP, PCCompositeSetType(), PCCompositeSpecialSetAlpha(), PCCompositeAddPC(),
    502:            PCCompositeGetPC(), PCCompositeSetUseTrue()

    504: M*/

    509: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_Composite(PC pc)
    510: {
    512:   PC_Composite   *jac;

    515:   PetscNew(PC_Composite,&jac);
    516:   PetscLogObjectMemory(pc,sizeof(PC_Composite));
    517:   pc->ops->apply              = PCApply_Composite_Additive;
    518:   pc->ops->setup              = PCSetUp_Composite;
    519:   pc->ops->destroy            = PCDestroy_Composite;
    520:   pc->ops->setfromoptions     = PCSetFromOptions_Composite;
    521:   pc->ops->view               = PCView_Composite;
    522:   pc->ops->applyrichardson    = 0;

    524:   pc->data               = (void*)jac;
    525:   jac->type              = PC_COMPOSITE_ADDITIVE;
    526:   jac->work1             = 0;
    527:   jac->work2             = 0;
    528:   jac->head              = 0;

    530:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeSetType_C","PCCompositeSetType_Composite",
    531:                     PCCompositeSetType_Composite);
    532:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeAddPC_C","PCCompositeAddPC_Composite",
    533:                     PCCompositeAddPC_Composite);
    534:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeGetPC_C","PCCompositeGetPC_Composite",
    535:                     PCCompositeGetPC_Composite);
    536:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeSetUseTrue_C","PCCompositeSetUseTrue_Composite",
    537:                     PCCompositeSetUseTrue_Composite);
    538:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCCompositeSpecialSetAlpha_C","PCCompositeSpecialSetAlpha_Composite",
    539:                     PCCompositeSpecialSetAlpha_Composite);

    541:   return(0);
    542: }