Actual source code: fieldsplit.c

      1: #define PETSCKSP_DLL

      3: /*

      5: */
     6:  #include src/ksp/pc/pcimpl.h

      8: typedef struct _PC_FieldSplitLink *PC_FieldSplitLink;
      9: struct _PC_FieldSplitLink {
     10:   KSP               ksp;
     11:   Vec               x,y;
     12:   PetscInt          nfields;
     13:   PetscInt          *fields;
     14:   VecScatter        sctx;
     15:   PC_FieldSplitLink next;
     16: };

     18: typedef struct {
     19:   PCCompositeType   type;              /* additive or multiplicative */
     20:   PetscTruth        defaultsplit;
     21:   PetscInt          bs;
     22:   PetscInt          nsplits;
     23:   Vec               *x,*y,w1,w2;
     24:   Mat               *pmat;
     25:   IS                *is;
     26:   PC_FieldSplitLink head;
     27: } PC_FieldSplit;

     31: static PetscErrorCode PCView_FieldSplit(PC pc,PetscViewer viewer)
     32: {
     33:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
     34:   PetscErrorCode    ierr;
     35:   PetscTruth        iascii;
     36:   PetscInt          i,j;
     37:   PC_FieldSplitLink ilink = jac->head;

     40:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
     41:   if (iascii) {
     42:     PetscViewerASCIIPrintf(viewer,"  FieldSplit with %s composition: total splits = %D",PCCompositeTypes[jac->type],jac->nsplits);
     43:     PetscViewerASCIIPrintf(viewer,"  Solver info for each split is in the following KSP objects:\n");
     44:     PetscViewerASCIIPushTab(viewer);
     45:     for (i=0; i<jac->nsplits; i++) {
     46:       PetscViewerASCIIPrintf(viewer,"Split number %D Fields ",i);
     47:       PetscViewerASCIIUseTabs(viewer,PETSC_FALSE);
     48:       for (j=0; j<ilink->nfields; j++) {
     49:         if (j > 0) {
     50:           PetscViewerASCIIPrintf(viewer,",");
     51:         }
     52:         PetscViewerASCIIPrintf(viewer," %D",ilink->fields[j]);
     53:       }
     54:       PetscViewerASCIIPrintf(viewer,"\n");
     55:       PetscViewerASCIIUseTabs(viewer,PETSC_TRUE);
     56:       KSPView(ilink->ksp,viewer);
     57:       ilink = ilink->next;
     58:     }
     59:     PetscViewerASCIIPopTab(viewer);
     60:   } else {
     61:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCFieldSplit",((PetscObject)viewer)->type_name);
     62:   }
     63:   return(0);
     64: }

     68: static PetscErrorCode PCFieldSplitSetDefaults(PC pc)
     69: {
     70:   PC_FieldSplit     *jac  = (PC_FieldSplit*)pc->data;
     71:   PetscErrorCode    ierr;
     72:   PC_FieldSplitLink ilink = jac->head;
     73:   PetscInt          i;
     74:   PetscTruth        flg = PETSC_FALSE,*fields;

     77:   PetscOptionsGetTruth(pc->prefix,"-pc_fieldsplit_default",&flg,PETSC_NULL);
     78:   if (!ilink || flg) {
     79:     PetscLogInfo((pc,"PCFieldSplitSetDefaults: Using default splitting of fields\n"));
     80:     if (jac->bs <= 0) {
     81:       MatGetBlockSize(pc->pmat,&jac->bs);
     82:     }
     83:     PetscMalloc(jac->bs*sizeof(PetscTruth),&fields);
     84:     PetscMemzero(fields,jac->bs*sizeof(PetscTruth));
     85:     while (ilink) {
     86:       for (i=0; i<ilink->nfields; i++) {
     87:         fields[ilink->fields[i]] = PETSC_TRUE;
     88:       }
     89:       ilink = ilink->next;
     90:     }
     91:     jac->defaultsplit = PETSC_TRUE;
     92:     for (i=0; i<jac->bs; i++) {
     93:       if (!fields[i]) {
     94:         PCFieldSplitSetFields(pc,1,&i);
     95:       } else {
     96:         jac->defaultsplit = PETSC_FALSE;
     97:       }
     98:     }
     99:   }
    100:   return(0);
    101: }

    106: static PetscErrorCode PCSetUp_FieldSplit(PC pc)
    107: {
    108:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
    109:   PetscErrorCode    ierr;
    110:   PC_FieldSplitLink ilink;
    111:   PetscInt          i,nsplit;
    112:   MatStructure      flag = pc->flag;

    115:   PCFieldSplitSetDefaults(pc);
    116:   nsplit = jac->nsplits;
    117:   ilink   = jac->head;

    119:   /* get the matrices for each split */
    120:   if (!jac->is) {
    121:     PetscInt rstart,rend,nslots,bs;

    123:     MatGetBlockSize(pc->pmat,&bs);
    124:     MatGetOwnershipRange(pc->pmat,&rstart,&rend);
    125:     nslots = (rend - rstart)/bs;
    126:     PetscMalloc(nsplit*sizeof(IS),&jac->is);
    127:     for (i=0; i<nsplit; i++) {
    128:       if (jac->defaultsplit) {
    129:         ISCreateStride(pc->comm,nslots,rstart+i,nsplit,&jac->is[i]);
    130:       } else {
    131:         PetscInt   *ii,j,k,nfields = ilink->nfields,*fields = ilink->fields;
    132:         PetscTruth sorted;
    133:         PetscMalloc(ilink->nfields*nslots*sizeof(PetscInt),&ii);
    134:         for (j=0; j<nslots; j++) {
    135:           for (k=0; k<nfields; k++) {
    136:             ii[nfields*j + k] = rstart + bs*j + fields[k];
    137:           }
    138:         }
    139:         ISCreateGeneral(pc->comm,nslots*nfields,ii,&jac->is[i]);
    140:         ISSorted(jac->is[i],&sorted);
    141:         if (!sorted) SETERRQ(PETSC_ERR_USER,"Fields must be sorted when creating split");
    142:         PetscFree(ii);
    143:         ilink = ilink->next;
    144:       }
    145:     }
    146:   }
    147:
    148:   if (!jac->pmat) {
    149:     MatGetSubMatrices(pc->pmat,nsplit,jac->is,jac->is,MAT_INITIAL_MATRIX,&jac->pmat);
    150:   } else {
    151:     MatGetSubMatrices(pc->pmat,nsplit,jac->is,jac->is,MAT_REUSE_MATRIX,&jac->pmat);
    152:   }

    154:   /* set up the individual PCs */
    155:   i    = 0;
    156:   ilink = jac->head;
    157:   while (ilink) {
    158:     KSPSetOperators(ilink->ksp,jac->pmat[i],jac->pmat[i],flag);
    159:     KSPSetFromOptions(ilink->ksp);
    160:     KSPSetUp(ilink->ksp);
    161:     i++;
    162:     ilink = ilink->next;
    163:   }
    164:
    165:   /* create work vectors for each split */
    166:   if (!jac->x) {
    167:     Vec xtmp;
    168:     PetscMalloc2(nsplit,Vec,&jac->x,nsplit,Vec,&jac->y);
    169:     ilink = jac->head;
    170:     for (i=0; i<nsplit; i++) {
    171:       Mat A;
    172:       KSPGetOperators(ilink->ksp,PETSC_NULL,&A,PETSC_NULL);
    173:       MatGetVecs(A,&ilink->x,&ilink->y);
    174:       jac->x[i] = ilink->x;
    175:       jac->y[i] = ilink->y;
    176:       ilink      = ilink->next;
    177:     }
    178:     /* compute scatter contexts needed by multiplicative versions and non-default splits */
    179:
    180:     ilink = jac->head;
    181:     MatGetVecs(pc->pmat,&xtmp,PETSC_NULL);
    182:     for (i=0; i<nsplit; i++) {
    183:       VecScatterCreate(xtmp,jac->is[i],jac->x[i],PETSC_NULL,&ilink->sctx);
    184:       ilink = ilink->next;
    185:     }
    186:     VecDestroy(xtmp);
    187:   }
    188:   return(0);
    189: }

    191: #define FieldSplitSplitSolveAdd(ilink,xx,yy) \
    192:     (VecScatterBegin(xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD,ilink->sctx) || \
    193:      VecScatterEnd(xx,ilink->x,INSERT_VALUES,SCATTER_FORWARD,ilink->sctx) || \
    194:      KSPSolve(ilink->ksp,ilink->x,ilink->y) || \
    195:      VecScatterBegin(ilink->y,yy,ADD_VALUES,SCATTER_REVERSE,ilink->sctx) || \
    196:      VecScatterEnd(ilink->y,yy,ADD_VALUES,SCATTER_REVERSE,ilink->sctx))

    200: static PetscErrorCode PCApply_FieldSplit(PC pc,Vec x,Vec y)
    201: {
    202:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
    203:   PetscErrorCode    ierr;
    204:   PC_FieldSplitLink ilink = jac->head;
    205:   PetscScalar       zero = 0.0,mone = -1.0;

    208:   if (jac->type == PC_COMPOSITE_ADDITIVE) {
    209:     if (jac->defaultsplit) {
    210:       VecStrideGatherAll(x,jac->x,INSERT_VALUES);
    211:       while (ilink) {
    212:         KSPSolve(ilink->ksp,ilink->x,ilink->y);
    213:         ilink = ilink->next;
    214:       }
    215:       VecStrideScatterAll(jac->y,y,INSERT_VALUES);
    216:     } else {
    217:       PetscInt    i = 0;

    219:       VecSet(y,zero);
    220:       while (ilink) {
    221:         FieldSplitSplitSolveAdd(ilink,x,y);
    222:         ilink = ilink->next;
    223:         i++;
    224:       }
    225:     }
    226:   } else {
    227:     if (!jac->w1) {
    228:       VecDuplicate(x,&jac->w1);
    229:       VecDuplicate(x,&jac->w2);
    230:     }
    231:     VecSet(y,zero);
    232:     FieldSplitSplitSolveAdd(ilink,x,y);
    233:     while (ilink->next) {
    234:       ilink = ilink->next;
    235:       MatMult(pc->pmat,y,jac->w1);
    236:       VecWAXPY(jac->w2,mone,jac->w1,x);
    237:       FieldSplitSplitSolveAdd(ilink,jac->w2,y);
    238:     }
    239:   }
    240:   return(0);
    241: }

    245: static PetscErrorCode PCDestroy_FieldSplit(PC pc)
    246: {
    247:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
    248:   PetscErrorCode    ierr;
    249:   PC_FieldSplitLink ilink = jac->head,next;

    252:   while (ilink) {
    253:     KSPDestroy(ilink->ksp);
    254:     if (ilink->x) {VecDestroy(ilink->x);}
    255:     if (ilink->y) {VecDestroy(ilink->y);}
    256:     if (ilink->sctx) {VecScatterDestroy(ilink->sctx);}
    257:     next = ilink->next;
    258:     PetscFree2(ilink,ilink->fields);
    259:     ilink = next;
    260:   }
    261:   if (jac->x) {PetscFree2(jac->x,jac->y);}
    262:   if (jac->pmat) {MatDestroyMatrices(jac->nsplits,&jac->pmat);}
    263:   if (jac->is) {
    264:     PetscInt i;
    265:     for (i=0; i<jac->nsplits; i++) {ISDestroy(jac->is[i]);}
    266:     PetscFree(jac->is);
    267:   }
    268:   if (jac->w1) {VecDestroy(jac->w1);}
    269:   if (jac->w2) {VecDestroy(jac->w2);}
    270:   PetscFree(jac);
    271:   return(0);
    272: }

    276: static PetscErrorCode PCSetFromOptions_FieldSplit(PC pc)
    277: /*   This does not call KSPSetFromOptions() on the subksp's, see PCSetFromOptionsBJacobi/ASM() */
    278: {
    280:   PetscInt       i = 0,nfields,fields[12];
    281:   PetscTruth     flg;
    282:   char           optionname[128];
    283:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

    286:   PetscOptionsHead("FieldSplit options");
    287:   PetscOptionsEnum("-pc_fieldsplit_type","Type of composition","PCFieldSplitSetType",PCCompositeTypes,(PetscEnum)jac->type,(PetscEnum*)&jac->type,&flg);
    288:   while (PETSC_TRUE) {
    289:     sprintf(optionname,"-pc_fieldsplit_%d_fields",(int)i++);
    290:     nfields = 12;
    291:     PetscOptionsIntArray(optionname,"Fields in this split","PCFieldSplitSetFields",fields,&nfields,&flg);
    292:     if (!flg) break;
    293:     if (!nfields) SETERRQ(PETSC_ERR_USER,"Cannot list zero fields");
    294:     PCFieldSplitSetFields(pc,nfields,fields);
    295:   }
    296:   PetscOptionsTail();
    297:   return(0);
    298: }

    300: /*------------------------------------------------------------------------------------*/

    305: PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetFields_FieldSplit(PC pc,PetscInt n,PetscInt *fields)
    306: {
    307:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
    308:   PetscErrorCode    ierr;
    309:   PC_FieldSplitLink ilink,next = jac->head;
    310:   char              prefix[128];

    313:   if (n <= 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative number of fields requested");
    314:   PetscMalloc2(1,struct _PC_FieldSplitLink,&ilink,n,PetscInt,&ilink->fields);
    315:   PetscMemcpy(ilink->fields,fields,n*sizeof(PetscInt));
    316:   ilink->nfields = n;
    317:   ilink->next    = PETSC_NULL;
    318:   KSPCreate(pc->comm,&ilink->ksp);
    319:   KSPSetType(ilink->ksp,KSPPREONLY);

    321:   if (pc->prefix) {
    322:     sprintf(prefix,"%sfieldsplit_%d_",pc->prefix,(int)jac->nsplits);
    323:   } else {
    324:     sprintf(prefix,"fieldsplit_%d_",(int)jac->nsplits);
    325:   }
    326:   KSPSetOptionsPrefix(ilink->ksp,prefix);

    328:   if (!next) {
    329:     jac->head = ilink;
    330:   } else {
    331:     while (next->next) {
    332:       next = next->next;
    333:     }
    334:     next->next = ilink;
    335:   }
    336:   jac->nsplits++;
    337:   return(0);
    338: }

    345: PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP_FieldSplit(PC pc,PetscInt *n,KSP **subksp)
    346: {
    347:   PC_FieldSplit     *jac = (PC_FieldSplit*)pc->data;
    348:   PetscErrorCode    ierr;
    349:   PetscInt          cnt = 0;
    350:   PC_FieldSplitLink ilink = jac->head;

    353:   PetscMalloc(jac->nsplits*sizeof(KSP*),subksp);
    354:   while (ilink) {
    355:     (*subksp)[cnt++] = ilink->ksp;
    356:     ilink = ilink->next;
    357:   }
    358:   if (cnt != jac->nsplits) SETERRQ2(PETSC_ERR_PLIB,"Corrupt PCFIELDSPLIT object: number splits in linked list %D in object %D",cnt,jac->nsplits);
    359:   *n = jac->nsplits;
    360:   return(0);
    361: }

    366: /*@
    367:     PCFieldSplitSetFields - Sets the fields for one particular split in the field split preconditioner

    369:     Collective on PC

    371:     Input Parameters:
    372: +   pc  - the preconditioner context
    373: .   n - the number of fields in this split
    374: .   fields - the fields in this split

    376:     Level: intermediate

    378: .seealso: PCFieldSplitGetSubKSP(), PCFIELDSPLIT

    380: @*/
    381: PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetFields(PC pc,PetscInt n, PetscInt *fields)
    382: {
    383:   PetscErrorCode ierr,(*f)(PC,PetscInt,PetscInt *);

    387:   PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitSetFields_C",(void (**)(void))&f);
    388:   if (f) {
    389:     (*f)(pc,n,fields);
    390:   }
    391:   return(0);
    392: }

    396: /*@C
    397:    PCFieldSplitGetSubKSP - Gets the KSP contexts for all splits
    398:
    399:    Collective on KSP

    401:    Input Parameter:
    402: .  pc - the preconditioner context

    404:    Output Parameters:
    405: +  n - the number of split
    406: -  pc - the array of KSP contexts

    408:    Note:
    409:    After PCFieldSplitGetSubKSP() the array of KSPs IS to be freed

    411:    You must call KSPSetUp() before calling PCFieldSplitGetSubKSP().

    413:    Level: advanced

    415: .seealso: PCFIELDSPLIT
    416: @*/
    417: PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitGetSubKSP(PC pc,PetscInt *n,KSP *subksp[])
    418: {
    419:   PetscErrorCode ierr,(*f)(PC,PetscInt*,KSP **);

    424:   PetscObjectQueryFunction((PetscObject)pc,"PCFieldSplitGetSubKSP_C",(void (**)(void))&f);
    425:   if (f) {
    426:     (*f)(pc,n,subksp);
    427:   } else {
    428:     SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot get subksp for this type of PC");
    429:   }
    430:   return(0);
    431: }

    436: PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetType_FieldSplit(PC pc,PCCompositeType type)
    437: {
    438:   PC_FieldSplit  *jac = (PC_FieldSplit*)pc->data;

    441:   jac->type = type;
    442:   return(0);
    443: }

    448: /*@C
    449:    PCFieldSplitSetType - Sets the type of fieldsplit preconditioner.
    450:
    451:    Collective on PC

    453:    Input Parameter:
    454: .  pc - the preconditioner context
    455: .  type - PC_COMPOSITE_ADDITIVE (default), PC_COMPOSITE_MULTIPLICATIVE

    457:    Options Database Key:
    458: .  -pc_fieldsplit_type <type: one of multiplicative, additive, special> - Sets fieldsplit preconditioner type

    460:    Level: Developer

    462: .keywords: PC, set, type, composite preconditioner, additive, multiplicative

    464: .seealso: PCCompositeSetType()

    466: @*/
    467: PetscErrorCode PETSCKSP_DLLEXPORT PCFieldSplitSetType(PC pc,PCCompositeType type)
    468: {
    469:   PetscErrorCode ierr,(*f)(PC,PCCompositeType);

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

    480: /* -------------------------------------------------------------------------------------*/
    481: /*MC
    482:    PCFIELDSPLIT - Preconditioner created by combining seperate preconditioners for individual
    483:                   fields or groups of fields

    486:      To set options on the solvers for each block append -sub_ to all the PC
    487:         options database keys. For example, -sub_pc_type ilu -sub_pc_ilu_levels 1
    488:
    489:      To set the options on the solvers seperate for each block call PCFieldSplitGetSubKSP()
    490:          and set the options directly on the resulting KSP object

    492:    Level: intermediate

    494:    Options Database Keys:
    495: +   -pc_splitfield_%d_fields <a,b,..> - indicates the fields to be used in the %d'th split
    496: .   -pc_splitfield_default - automatically add any fields to additional splits that have not
    497:                               been supplied explicitly by -pc_splitfield_%d_fields
    498: -   -pc_splitfield_type <additive,multiplicative>

    500:    Concepts: physics based preconditioners

    502: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
    503:            PCFieldSplitGetSubKSP(), PCFieldSplitSetFields(),PCFieldSplitSetType()
    504: M*/

    509: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_FieldSplit(PC pc)
    510: {
    512:   PC_FieldSplit  *jac;

    515:   PetscNew(PC_FieldSplit,&jac);
    516:   PetscLogObjectMemory(pc,sizeof(PC_FieldSplit));
    517:   jac->bs      = -1;
    518:   jac->nsplits = 0;
    519:   jac->type    = PC_COMPOSITE_ADDITIVE;
    520:   pc->data     = (void*)jac;

    522:   pc->ops->apply             = PCApply_FieldSplit;
    523:   pc->ops->setup             = PCSetUp_FieldSplit;
    524:   pc->ops->destroy           = PCDestroy_FieldSplit;
    525:   pc->ops->setfromoptions    = PCSetFromOptions_FieldSplit;
    526:   pc->ops->view              = PCView_FieldSplit;
    527:   pc->ops->applyrichardson   = 0;

    529:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitGetSubKSP_C","PCFieldSplitGetSubKSP_FieldSplit",
    530:                     PCFieldSplitGetSubKSP_FieldSplit);
    531:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetFields_C","PCFieldSplitSetFields_FieldSplit",
    532:                     PCFieldSplitSetFields_FieldSplit);
    533:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCFieldSplitSetType_C","PCFieldSplitSetType_FieldSplit",
    534:                     PCFieldSplitSetType_FieldSplit);
    535:   return(0);
    536: }