Actual source code: asm.c

      1: #define PETSCKSP_DLL

      3: /*
      4:   This file defines an additive Schwarz preconditioner for any Mat implementation.

      6:   Note that each processor may have any number of subdomains. But in order to
      7:   deal easily with the VecScatter(), we treat each processor as if it has the
      8:   same number of subdomains.

     10:        n - total number of true subdomains on all processors
     11:        n_local_true - actual number of subdomains on this processor
     12:        n_local = maximum over all processors of n_local_true
     13: */
     14:  #include src/ksp/pc/pcimpl.h

     16: typedef struct {
     17:   PetscInt   n,n_local,n_local_true;
     18:   PetscTruth is_flg;              /* flg set to 1 if the IS created in pcsetup */
     19:   PetscInt   overlap;             /* overlap requested by user */
     20:   KSP        *ksp;               /* linear solvers for each block */
     21:   VecScatter *scat;               /* mapping to subregion */
     22:   Vec        *x,*y;
     23:   IS         *is;                 /* index set that defines each subdomain */
     24:   Mat        *mat,*pmat;          /* mat is not currently used */
     25:   PCASMType  type;                /* use reduced interpolation, restriction or both */
     26:   PetscTruth type_set;            /* if user set this value (so won't change it for symmetric problems) */
     27:   PetscTruth same_local_solves;   /* flag indicating whether all local solvers are same */
     28:   PetscTruth inplace;             /* indicates that the sub-matrices are deleted after
     29:                                      PCSetUpOnBlocks() is done. Similar to inplace
     30:                                      factorization in the case of LU and ILU */
     31: } PC_ASM;

     35: static PetscErrorCode PCView_ASM(PC pc,PetscViewer viewer)
     36: {
     37:   PC_ASM         *jac = (PC_ASM*)pc->data;
     39:   PetscMPIInt    rank;
     40:   PetscInt       i;
     41:   PetscTruth     iascii,isstring;
     42:   PetscViewer    sviewer;

     46:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
     47:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
     48:   if (iascii) {
     49:     if (jac->n > 0) {
     50:       PetscViewerASCIIPrintf(viewer,"  Additive Schwarz: total subdomain blocks = %D, amount of overlap = %D\n",jac->n,jac->overlap);
     51:     } else {
     52:       PetscViewerASCIIPrintf(viewer,"  Additive Schwarz: total subdomain blocks not yet set, amount of overlap = %D\n",jac->overlap);
     53:     }
     54:     PetscViewerASCIIPrintf(viewer,"  Additive Schwarz: restriction/interpolation type - %s\n",PCASMTypes[jac->type]);
     55:     MPI_Comm_rank(pc->comm,&rank);
     56:     if (jac->same_local_solves) {
     57:       PetscViewerASCIIPrintf(viewer,"  Local solve is same for all blocks, in the following KSP and PC objects:\n");
     58:       PetscViewerGetSingleton(viewer,&sviewer);
     59:       if (!rank && jac->ksp) {
     60:         PetscViewerASCIIPushTab(viewer);
     61:         KSPView(jac->ksp[0],sviewer);
     62:         PetscViewerASCIIPopTab(viewer);
     63:       }
     64:       PetscViewerRestoreSingleton(viewer,&sviewer);
     65:     } else {
     66:       PetscViewerASCIIPrintf(viewer,"  Local solve info for each block is in the following KSP and PC objects:\n");
     67:       PetscViewerASCIISynchronizedPrintf(viewer,"[%d] number of local blocks = %D\n",rank,jac->n_local_true);
     68:       PetscViewerASCIIPushTab(viewer);
     69:       for (i=0; i<jac->n_local; i++) {
     70:         PetscViewerGetSingleton(viewer,&sviewer);
     71:         if (i < jac->n_local_true) {
     72:           PetscViewerASCIISynchronizedPrintf(sviewer,"[%d] local block number %D\n",rank,i);
     73:           KSPView(jac->ksp[i],sviewer);
     74:           PetscViewerASCIISynchronizedPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
     75:         }
     76:         PetscViewerRestoreSingleton(viewer,&sviewer);
     77:       }
     78:       PetscViewerASCIIPopTab(viewer);
     79:       PetscViewerFlush(viewer);
     80:     }
     81:   } else if (isstring) {
     82:     PetscViewerStringSPrintf(viewer," blks=%D, overlap=%D, type=%D",jac->n,jac->overlap,jac->type);
     83:     PetscViewerGetSingleton(viewer,&sviewer);
     84:       if (jac->ksp) {KSPView(jac->ksp[0],sviewer);}
     85:     PetscViewerGetSingleton(viewer,&sviewer);
     86:   } else {
     87:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCASM",((PetscObject)viewer)->type_name);
     88:   }
     89:   return(0);
     90: }

     94: static PetscErrorCode PCSetUp_ASM(PC pc)
     95: {
     96:   PC_ASM         *osm  = (PC_ASM*)pc->data;
     98:   PetscInt       i,m,n_local = osm->n_local,n_local_true = osm->n_local_true;
     99:   PetscInt       start,start_val,end_val,sz,bs;
    100:   PetscMPIInt    size;
    101:   MatReuse       scall = MAT_REUSE_MATRIX;
    102:   IS             isl;
    103:   KSP            ksp;
    104:   PC             subpc;
    105:   const char     *prefix,*pprefix;
    106:   Vec            vec;

    109:   MatGetVecs(pc->pmat,&vec,0);
    110:   if (!pc->setupcalled) {
    111:     if (osm->n == PETSC_DECIDE && osm->n_local_true == PETSC_DECIDE) {
    112:       /* no subdomains given, use one per processor */
    113:       osm->n_local_true = osm->n_local = 1;
    114:       MPI_Comm_size(pc->comm,&size);
    115:       osm->n = size;
    116:     } else if (osm->n == PETSC_DECIDE) { /* determine global number of subdomains */
    117:       PetscInt inwork[2],outwork[2];
    118:       inwork[0] = inwork[1] = osm->n_local_true;
    119:       MPI_Allreduce(inwork,outwork,1,MPIU_2INT,PetscMaxSum_Op,pc->comm);
    120:       osm->n_local = outwork[0];
    121:       osm->n       = outwork[1];
    122:     }
    123:     n_local      = osm->n_local;
    124:     n_local_true = osm->n_local_true;
    125:     if (!osm->is){ /* build the index sets */
    126:       PetscMalloc((n_local_true+1)*sizeof(IS **),&osm->is);
    127:       MatGetOwnershipRange(pc->pmat,&start_val,&end_val);
    128:       MatGetBlockSize(pc->pmat,&bs);
    129:       sz    = end_val - start_val;
    130:       start = start_val;
    131:       if (end_val/bs*bs != end_val || start_val/bs*bs != start_val) {
    132:         SETERRQ(PETSC_ERR_ARG_WRONG,"Bad distribution for matrix block size");
    133:       }
    134:       for (i=0; i<n_local_true; i++){
    135:         size       =  ((sz/bs)/n_local_true + (((sz/bs) % n_local_true) > i))*bs;
    136:          ISCreateStride(PETSC_COMM_SELF,size,start,1,&isl);
    137:         start      += size;
    138:         osm->is[i] =  isl;
    139:       }
    140:       osm->is_flg = PETSC_TRUE;
    141:     }

    143:     PetscMalloc((n_local_true+1)*sizeof(KSP **),&osm->ksp);
    144:     PetscMalloc(n_local*sizeof(VecScatter **),&osm->scat);
    145:     PetscMalloc(2*n_local*sizeof(Vec **),&osm->x);
    146:     osm->y = osm->x + n_local;

    148:     /*  Extend the "overlapping" regions by a number of steps  */
    149:     MatIncreaseOverlap(pc->pmat,n_local_true,osm->is,osm->overlap);
    150:     for (i=0; i<n_local_true; i++) {
    151:       ISSort(osm->is[i]);
    152:     }

    154:     /* create the local work vectors and scatter contexts */
    155:     for (i=0; i<n_local_true; i++) {
    156:       ISGetLocalSize(osm->is[i],&m);
    157:       VecCreateSeq(PETSC_COMM_SELF,m,&osm->x[i]);
    158:       VecDuplicate(osm->x[i],&osm->y[i]);
    159:       ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
    160:       VecScatterCreate(vec,osm->is[i],osm->x[i],isl,&osm->scat[i]);
    161:       ISDestroy(isl);
    162:     }
    163:     for (i=n_local_true; i<n_local; i++) {
    164:       VecCreateSeq(PETSC_COMM_SELF,0,&osm->x[i]);
    165:       VecDuplicate(osm->x[i],&osm->y[i]);
    166:       ISCreateStride(PETSC_COMM_SELF,0,0,1,&isl);
    167:       VecScatterCreate(vec,isl,osm->x[i],isl,&osm->scat[i]);
    168:       ISDestroy(isl);
    169:     }

    171:    /*
    172:        Create the local solvers.
    173:     */
    174:     for (i=0; i<n_local_true; i++) {
    175:       KSPCreate(PETSC_COMM_SELF,&ksp);
    176:       PetscLogObjectParent(pc,ksp);
    177:       KSPSetType(ksp,KSPPREONLY);
    178:       KSPGetPC(ksp,&subpc);
    179:       PCGetOptionsPrefix(pc,&prefix);
    180:       KSPSetOptionsPrefix(ksp,prefix);
    181:       KSPAppendOptionsPrefix(ksp,"sub_");
    182:       osm->ksp[i] = ksp;
    183:     }
    184:     scall = MAT_INITIAL_MATRIX;
    185:   } else {
    186:     /*
    187:        Destroy the blocks from the previous iteration
    188:     */
    189:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
    190:       MatDestroyMatrices(osm->n_local_true,&osm->pmat);
    191:       scall = MAT_INITIAL_MATRIX;
    192:     }
    193:   }

    195:   /* extract out the submatrices */
    196:   MatGetSubMatrices(pc->pmat,osm->n_local_true,osm->is,osm->is,scall,&osm->pmat);

    198:   /* Return control to the user so that the submatrices can be modified (e.g., to apply
    199:      different boundary conditions for the submatrices than for the global problem) */
    200:   PCModifySubMatrices(pc,osm->n_local,osm->is,osm->is,osm->pmat,pc->modifysubmatricesP);

    202:   /* loop over subdomains putting them into local ksp */
    203:   PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
    204:   for (i=0; i<n_local_true; i++) {
    205:     PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);
    206:     PetscLogObjectParent(pc,osm->pmat[i]);
    207:     KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i],pc->flag);
    208:     KSPSetFromOptions(osm->ksp[i]);
    209:   }
    210:   VecDestroy(vec);
    211:   return(0);
    212: }

    216: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
    217: {
    218:   PC_ASM         *osm = (PC_ASM*)pc->data;
    220:   PetscInt       i;

    223:   for (i=0; i<osm->n_local_true; i++) {
    224:     KSPSetUp(osm->ksp[i]);
    225:   }
    226:   /*
    227:      If inplace flag is set, then destroy the matrix after the setup
    228:      on blocks is done.
    229:   */
    230:   if (osm->inplace && osm->n_local_true > 0) {
    231:     MatDestroyMatrices(osm->n_local_true,&osm->pmat);
    232:   }
    233:   return(0);
    234: }

    238: static PetscErrorCode PCApply_ASM(PC pc,Vec x,Vec y)
    239: {
    240:   PC_ASM         *osm = (PC_ASM*)pc->data;
    242:   PetscInt       i,n_local = osm->n_local,n_local_true = osm->n_local_true;
    243:   PetscScalar    zero = 0.0;
    244:   ScatterMode    forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;

    247:   /*
    248:        Support for limiting the restriction or interpolation to only local
    249:      subdomain values (leaving the other values 0).
    250:   */
    251:   if (!(osm->type & PC_ASM_RESTRICT)) {
    252:     forward = SCATTER_FORWARD_LOCAL;
    253:     /* have to zero the work RHS since scatter may leave some slots empty */
    254:     for (i=0; i<n_local_true; i++) {
    255:       VecSet(osm->x[i],zero);
    256:     }
    257:   }
    258:   if (!(osm->type & PC_ASM_INTERPOLATE)) {
    259:     reverse = SCATTER_REVERSE_LOCAL;
    260:   }

    262:   for (i=0; i<n_local; i++) {
    263:     VecScatterBegin(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
    264:   }
    265:   VecSet(y,zero);
    266:   /* do the local solves */
    267:   for (i=0; i<n_local_true; i++) {
    268:     VecScatterEnd(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
    269:     KSPSolve(osm->ksp[i],osm->x[i],osm->y[i]);
    270:     VecScatterBegin(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
    271:   }
    272:   /* handle the rest of the scatters that do not have local solves */
    273:   for (i=n_local_true; i<n_local; i++) {
    274:     VecScatterEnd(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
    275:     VecScatterBegin(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
    276:   }
    277:   for (i=0; i<n_local; i++) {
    278:     VecScatterEnd(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
    279:   }
    280:   return(0);
    281: }

    285: static PetscErrorCode PCApplyTranspose_ASM(PC pc,Vec x,Vec y)
    286: {
    287:   PC_ASM         *osm = (PC_ASM*)pc->data;
    289:   PetscInt       i,n_local = osm->n_local,n_local_true = osm->n_local_true;
    290:   PetscScalar    zero = 0.0;
    291:   ScatterMode    forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;

    294:   /*
    295:        Support for limiting the restriction or interpolation to only local
    296:      subdomain values (leaving the other values 0).

    298:        Note: these are reversed from the PCApply_ASM() because we are applying the
    299:      transpose of the three terms
    300:   */
    301:   if (!(osm->type & PC_ASM_INTERPOLATE)) {
    302:     forward = SCATTER_FORWARD_LOCAL;
    303:     /* have to zero the work RHS since scatter may leave some slots empty */
    304:     for (i=0; i<n_local_true; i++) {
    305:       VecSet(osm->x[i],zero);
    306:     }
    307:   }
    308:   if (!(osm->type & PC_ASM_RESTRICT)) {
    309:     reverse = SCATTER_REVERSE_LOCAL;
    310:   }

    312:   for (i=0; i<n_local; i++) {
    313:     VecScatterBegin(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
    314:   }
    315:   VecSet(y,zero);
    316:   /* do the local solves */
    317:   for (i=0; i<n_local_true; i++) {
    318:     VecScatterEnd(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
    319:     KSPSolveTranspose(osm->ksp[i],osm->x[i],osm->y[i]);
    320:     VecScatterBegin(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
    321:   }
    322:   /* handle the rest of the scatters that do not have local solves */
    323:   for (i=n_local_true; i<n_local; i++) {
    324:     VecScatterEnd(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
    325:     VecScatterBegin(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
    326:   }
    327:   for (i=0; i<n_local; i++) {
    328:     VecScatterEnd(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
    329:   }
    330:   return(0);
    331: }

    335: static PetscErrorCode PCDestroy_ASM(PC pc)
    336: {
    337:   PC_ASM         *osm = (PC_ASM*)pc->data;
    339:   PetscInt       i;

    342:   for (i=0; i<osm->n_local; i++) {
    343:     VecScatterDestroy(osm->scat[i]);
    344:     VecDestroy(osm->x[i]);
    345:     VecDestroy(osm->y[i]);
    346:   }
    347:   if (osm->n_local_true > 0 && !osm->inplace && osm->pmat) {
    348:     MatDestroyMatrices(osm->n_local_true,&osm->pmat);
    349:   }
    350:   if (osm->ksp) {
    351:     for (i=0; i<osm->n_local_true; i++) {
    352:       KSPDestroy(osm->ksp[i]);
    353:     }
    354:   }
    355:   if (osm->is_flg) {
    356:     for (i=0; i<osm->n_local_true; i++) {ISDestroy(osm->is[i]);}
    357:     PetscFree(osm->is);
    358:   }
    359:   if (osm->ksp) {PetscFree(osm->ksp);}
    360:   if (osm->scat) {PetscFree(osm->scat);}
    361:   if (osm->x) {PetscFree(osm->x);}
    362:   PetscFree(osm);
    363:   return(0);
    364: }

    368: static PetscErrorCode PCSetFromOptions_ASM(PC pc)
    369: {
    370:   PC_ASM         *osm = (PC_ASM*)pc->data;
    372:   PetscInt       blocks,ovl;
    373:   PetscTruth     flg,set,sym;

    376:   /* set the type to symmetric if matrix is symmetric */
    377:   if (pc->pmat && !osm->type_set) {
    378:     MatIsSymmetricKnown(pc->pmat,&set,&sym);
    379:     if (set && sym) {
    380:       osm->type = PC_ASM_BASIC;
    381:     }
    382:   }
    383:   PetscOptionsHead("Additive Schwarz options");
    384:     PetscOptionsInt("-pc_asm_blocks","Number of subdomains","PCASMSetTotalSubdomains",osm->n,&blocks,&flg);
    385:     if (flg) {PCASMSetTotalSubdomains(pc,blocks,PETSC_NULL); }
    386:     PetscOptionsInt("-pc_asm_overlap","Number of grid points overlap","PCASMSetOverlap",osm->overlap,&ovl,&flg);
    387:     if (flg) {PCASMSetOverlap(pc,ovl); }
    388:     PetscOptionsName("-pc_asm_in_place","Perform matrix factorization inplace","PCASMSetUseInPlace",&flg);
    389:     if (flg) {PCASMSetUseInPlace(pc); }
    390:     PetscOptionsEnum("-pc_asm_type","Type of restriction/extension","PCASMSetType",PCASMTypes,(PetscEnum)osm->type,(PetscEnum*)&osm->type,&flg);
    391:   PetscOptionsTail();
    392:   return(0);
    393: }

    395: /*------------------------------------------------------------------------------------*/

    400: PetscErrorCode PETSCKSP_DLLEXPORT PCASMSetLocalSubdomains_ASM(PC pc,PetscInt n,IS is[])
    401: {
    402:   PC_ASM *osm = (PC_ASM*)pc->data;

    405:   if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Each process must have 0 or more blocks");

    407:   if (pc->setupcalled && n != osm->n_local_true) {
    408:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"PCASMSetLocalSubdomains() should be called before calling PCSetup().");
    409:   }
    410:   if (!pc->setupcalled){
    411:     osm->n_local_true = n;
    412:     osm->is           = is;
    413:   }
    414:   return(0);
    415: }

    421: PetscErrorCode PETSCKSP_DLLEXPORT PCASMSetTotalSubdomains_ASM(PC pc,PetscInt N,IS *is)
    422: {
    423:   PC_ASM         *osm = (PC_ASM*)pc->data;
    425:   PetscMPIInt    rank,size;
    426:   PetscInt       n;

    430:   if (is) SETERRQ(PETSC_ERR_SUP,"Use PCASMSetLocalSubdomains() to set specific index sets\n\
    431: they cannot be set globally yet.");

    433:   /*
    434:      Split the subdomains equally amoung all processors
    435:   */
    436:   MPI_Comm_rank(pc->comm,&rank);
    437:   MPI_Comm_size(pc->comm,&size);
    438:   n = N/size + ((N % size) > rank);
    439:   if (pc->setupcalled && n != osm->n_local_true) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"PCASMSetTotalSubdomains() should be called before PCSetup().");
    440:   if (!pc->setupcalled){
    441:     osm->n_local_true = n;
    442:     osm->is           = 0;
    443:   }
    444:   return(0);
    445: }

    451: PetscErrorCode PETSCKSP_DLLEXPORT PCASMSetOverlap_ASM(PC pc,PetscInt ovl)
    452: {
    453:   PC_ASM *osm;

    456:   if (ovl < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested");

    458:   osm               = (PC_ASM*)pc->data;
    459:   osm->overlap      = ovl;
    460:   return(0);
    461: }

    467: PetscErrorCode PETSCKSP_DLLEXPORT PCASMSetType_ASM(PC pc,PCASMType type)
    468: {
    469:   PC_ASM *osm;

    472:   osm           = (PC_ASM*)pc->data;
    473:   osm->type     = type;
    474:   osm->type_set = PETSC_TRUE;
    475:   return(0);
    476: }

    482: PetscErrorCode PETSCKSP_DLLEXPORT PCASMGetSubKSP_ASM(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
    483: {
    484:   PC_ASM         *jac = (PC_ASM*)pc->data;

    488:   if (jac->n_local_true < 0) {
    489:     SETERRQ(PETSC_ERR_ORDER,"Need to call PCSetUP() on PC (or KSPSetUp() on the outer KSP object) before calling here");
    490:   }

    492:   if (n_local)     *n_local     = jac->n_local_true;
    493:   if (first_local) {
    494:     MPI_Scan(&jac->n_local_true,first_local,1,MPIU_INT,MPI_SUM,pc->comm);
    495:     *first_local -= jac->n_local_true;
    496:   }
    497:   *ksp                         = jac->ksp;
    498:   jac->same_local_solves        = PETSC_FALSE; /* Assume that local solves are now different;
    499:                                       not necessarily true though!  This flag is
    500:                                       used only for PCView_ASM() */
    501:   return(0);
    502: }

    508: PetscErrorCode PETSCKSP_DLLEXPORT PCASMSetUseInPlace_ASM(PC pc)
    509: {
    510:   PC_ASM *dir;

    513:   dir          = (PC_ASM*)pc->data;
    514:   dir->inplace = PETSC_TRUE;
    515:   return(0);
    516: }

    519: /*----------------------------------------------------------------------------*/
    522: /*@
    523:    PCASMSetUseInPlace - Tells the system to destroy the matrix after setup is done.

    525:    Collective on PC

    527:    Input Parameters:
    528: .  pc - the preconditioner context

    530:    Options Database Key:
    531: .  -pc_asm_in_place - Activates in-place factorization

    533:    Note:
    534:    PCASMSetUseInplace() can only be used with the KSP method KSPPREONLY, and
    535:    when the original matrix is not required during the Solve process.
    536:    This destroys the matrix, early thus, saving on memory usage.

    538:    Level: intermediate

    540: .keywords: PC, set, factorization, direct, inplace, in-place, ASM

    542: .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace ()
    543: @*/
    544: PetscErrorCode PETSCKSP_DLLEXPORT PCASMSetUseInPlace(PC pc)
    545: {
    546:   PetscErrorCode ierr,(*f)(PC);

    550:   PetscObjectQueryFunction((PetscObject)pc,"PCASMSetUseInPlace_C",(void (**)(void))&f);
    551:   if (f) {
    552:     (*f)(pc);
    553:   }
    554:   return(0);
    555: }
    556: /*----------------------------------------------------------------------------*/

    560: /*@C
    561:     PCASMSetLocalSubdomains - Sets the local subdomains (for this processor
    562:     only) for the additive Schwarz preconditioner.

    564:     Collective on PC

    566:     Input Parameters:
    567: +   pc - the preconditioner context
    568: .   n - the number of subdomains for this processor (default value = 1)
    569: -   is - the index sets that define the subdomains for this processor
    570:          (or PETSC_NULL for PETSc to determine subdomains)

    572:     Notes:
    573:     The IS numbering is in the parallel, global numbering of the vector.

    575:     By default the ASM preconditioner uses 1 block per processor.

    577:     These index sets cannot be destroyed until after completion of the
    578:     linear solves for which the ASM preconditioner is being used.

    580:     Use PCASMSetTotalSubdomains() to set the subdomains for all processors.

    582:     Level: advanced

    584: .keywords: PC, ASM, set, local, subdomains, additive Schwarz

    586: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
    587:           PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains()
    588: @*/
    589: PetscErrorCode PETSCKSP_DLLEXPORT PCASMSetLocalSubdomains(PC pc,PetscInt n,IS is[])
    590: {
    591:   PetscErrorCode ierr,(*f)(PC,PetscInt,IS[]);

    595:   PetscObjectQueryFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",(void (**)(void))&f);
    596:   if (f) {
    597:     (*f)(pc,n,is);
    598:   }
    599:   return(0);
    600: }

    604: /*@C
    605:     PCASMSetTotalSubdomains - Sets the subdomains for all processor for the
    606:     additive Schwarz preconditioner.  Either all or no processors in the
    607:     PC communicator must call this routine, with the same index sets.

    609:     Collective on PC

    611:     Input Parameters:
    612: +   pc - the preconditioner context
    613: .   n - the number of subdomains for all processors
    614: -   is - the index sets that define the subdomains for all processor
    615:          (or PETSC_NULL for PETSc to determine subdomains)

    617:     Options Database Key:
    618:     To set the total number of subdomain blocks rather than specify the
    619:     index sets, use the option
    620: .    -pc_asm_blocks <blks> - Sets total blocks

    622:     Notes:
    623:     Currently you cannot use this to set the actual subdomains with the argument is.

    625:     By default the ASM preconditioner uses 1 block per processor.

    627:     These index sets cannot be destroyed until after completion of the
    628:     linear solves for which the ASM preconditioner is being used.

    630:     Use PCASMSetLocalSubdomains() to set local subdomains.

    632:     Level: advanced

    634: .keywords: PC, ASM, set, total, global, subdomains, additive Schwarz

    636: .seealso: PCASMSetLocalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
    637:           PCASMCreateSubdomains2D()
    638: @*/
    639: PetscErrorCode PETSCKSP_DLLEXPORT PCASMSetTotalSubdomains(PC pc,PetscInt N,IS *is)
    640: {
    641:   PetscErrorCode ierr,(*f)(PC,PetscInt,IS *);

    645:   PetscObjectQueryFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",(void (**)(void))&f);
    646:   if (f) {
    647:     (*f)(pc,N,is);
    648:   }
    649:   return(0);
    650: }

    654: /*@
    655:     PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
    656:     additive Schwarz preconditioner.  Either all or no processors in the
    657:     PC communicator must call this routine.

    659:     Collective on PC

    661:     Input Parameters:
    662: +   pc  - the preconditioner context
    663: -   ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1)

    665:     Options Database Key:
    666: .   -pc_asm_overlap <ovl> - Sets overlap

    668:     Notes:
    669:     By default the ASM preconditioner uses 1 block per processor.  To use
    670:     multiple blocks per perocessor, see PCASMSetTotalSubdomains() and
    671:     PCASMSetLocalSubdomains() (and the option -pc_asm_blocks <blks>).

    673:     The overlap defaults to 1, so if one desires that no additional
    674:     overlap be computed beyond what may have been set with a call to
    675:     PCASMSetTotalSubdomains() or PCASMSetLocalSubdomains(), then ovl
    676:     must be set to be 0.  In particular, if one does not explicitly set
    677:     the subdomains an application code, then all overlap would be computed
    678:     internally by PETSc, and using an overlap of 0 would result in an ASM
    679:     variant that is equivalent to the block Jacobi preconditioner.

    681:     Note that one can define initial index sets with any overlap via
    682:     PCASMSetTotalSubdomains() or PCASMSetLocalSubdomains(); the routine
    683:     PCASMSetOverlap() merely allows PETSc to extend that overlap further
    684:     if desired.

    686:     Level: intermediate

    688: .keywords: PC, ASM, set, overlap

    690: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
    691:           PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains()
    692: @*/
    693: PetscErrorCode PETSCKSP_DLLEXPORT PCASMSetOverlap(PC pc,PetscInt ovl)
    694: {
    695:   PetscErrorCode ierr,(*f)(PC,PetscInt);

    699:   PetscObjectQueryFunction((PetscObject)pc,"PCASMSetOverlap_C",(void (**)(void))&f);
    700:   if (f) {
    701:     (*f)(pc,ovl);
    702:   }
    703:   return(0);
    704: }

    708: /*@
    709:     PCASMSetType - Sets the type of restriction and interpolation used
    710:     for local problems in the additive Schwarz method.

    712:     Collective on PC

    714:     Input Parameters:
    715: +   pc  - the preconditioner context
    716: -   type - variant of ASM, one of
    717: .vb
    718:       PC_ASM_BASIC       - full interpolation and restriction
    719:       PC_ASM_RESTRICT    - full restriction, local processor interpolation
    720:       PC_ASM_INTERPOLATE - full interpolation, local processor restriction
    721:       PC_ASM_NONE        - local processor restriction and interpolation
    722: .ve

    724:     Options Database Key:
    725: .   -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type

    727:     Level: intermediate

    729: .keywords: PC, ASM, set, type

    731: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
    732:           PCASMCreateSubdomains2D()
    733: @*/
    734: PetscErrorCode PETSCKSP_DLLEXPORT PCASMSetType(PC pc,PCASMType type)
    735: {
    736:   PetscErrorCode ierr,(*f)(PC,PCASMType);

    740:   PetscObjectQueryFunction((PetscObject)pc,"PCASMSetType_C",(void (**)(void))&f);
    741:   if (f) {
    742:     (*f)(pc,type);
    743:   }
    744:   return(0);
    745: }

    749: /*@C
    750:    PCASMGetSubKSP - Gets the local KSP contexts for all blocks on
    751:    this processor.
    752:
    753:    Collective on PC iff first_local is requested

    755:    Input Parameter:
    756: .  pc - the preconditioner context

    758:    Output Parameters:
    759: +  n_local - the number of blocks on this processor or PETSC_NULL
    760: .  first_local - the global number of the first block on this processor or PETSC_NULL,
    761:                  all processors must request or all must pass PETSC_NULL
    762: -  ksp - the array of KSP contexts

    764:    Note:
    765:    After PCASMGetSubKSP() the array of KSPes is not to be freed

    767:    Currently for some matrix implementations only 1 block per processor
    768:    is supported.
    769:
    770:    You must call KSPSetUp() before calling PCASMGetSubKSP().

    772:    Level: advanced

    774: .keywords: PC, ASM, additive Schwarz, get, sub, KSP, context

    776: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap(),
    777:           PCASMCreateSubdomains2D(),
    778: @*/
    779: PetscErrorCode PETSCKSP_DLLEXPORT PCASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
    780: {
    781:   PetscErrorCode ierr,(*f)(PC,PetscInt*,PetscInt*,KSP **);

    785:   PetscObjectQueryFunction((PetscObject)pc,"PCASMGetSubKSP_C",(void (**)(void))&f);
    786:   if (f) {
    787:     (*f)(pc,n_local,first_local,ksp);
    788:   } else {
    789:     SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot get subksp for this type of PC");
    790:   }

    792:  return(0);
    793: }

    795: /* -------------------------------------------------------------------------------------*/
    796: /*MC
    797:    PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
    798:            its own KSP object.

    800:    Options Database Keys:
    801: +  -pc_asm_truelocal - Activates PCASMSetUseTrueLocal()
    802: .  -pc_asm_in_place - Activates in-place factorization
    803: .  -pc_asm_blocks <blks> - Sets total blocks
    804: .  -pc_asm_overlap <ovl> - Sets overlap
    805: -  -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type

    807:      IMPORTANT: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you
    808:       will get a different convergence rate due to the default option of -pc_asm_type restrict. Use
    809:       -pc_asm_type basic to use the standard ASM.

    811:    Notes: Each processor can have one or more blocks, but a block cannot be shared by more
    812:      than one processor. Defaults to one block per processor.

    814:      To set options on the solvers for each block append -sub_ to all the KSP, and PC
    815:         options database keys. For example, -sub_pc_type ilu -sub_pc_ilu_levels 1 -sub_ksp_type preonly
    816:
    817:      To set the options on the solvers seperate for each block call PCASMGetSubKSP()
    818:          and set the options directly on the resulting KSP object (you can access its PC
    819:          with KSPGetPC())

    822:    Level: beginner

    824:    Concepts: additive Schwarz method

    826: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
    827:            PCBJACOBI, PCASMSetUseTrueLocal(), PCASMGetSubKSP(), PCASMSetLocalSubdomains(),
    828:            PCASMSetTotalSubdomains(), PCSetModifySubmatrices(), PCASMSetOverlap(), PCASMSetType(),
    829:            PCASMSetUseInPlace()
    830: M*/

    835: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_ASM(PC pc)
    836: {
    838:   PC_ASM         *osm;

    841:   PetscNew(PC_ASM,&osm);
    842:   PetscLogObjectMemory(pc,sizeof(PC_ASM));
    843:   osm->n                 = PETSC_DECIDE;
    844:   osm->n_local           = 0;
    845:   osm->n_local_true      = PETSC_DECIDE;
    846:   osm->overlap           = 1;
    847:   osm->is_flg            = PETSC_FALSE;
    848:   osm->ksp              = 0;
    849:   osm->scat              = 0;
    850:   osm->is                = 0;
    851:   osm->mat               = 0;
    852:   osm->pmat              = 0;
    853:   osm->type              = PC_ASM_RESTRICT;
    854:   osm->same_local_solves = PETSC_TRUE;
    855:   osm->inplace           = PETSC_FALSE;
    856:   pc->data               = (void*)osm;

    858:   pc->ops->apply             = PCApply_ASM;
    859:   pc->ops->applytranspose    = PCApplyTranspose_ASM;
    860:   pc->ops->setup             = PCSetUp_ASM;
    861:   pc->ops->destroy           = PCDestroy_ASM;
    862:   pc->ops->setfromoptions    = PCSetFromOptions_ASM;
    863:   pc->ops->setuponblocks     = PCSetUpOnBlocks_ASM;
    864:   pc->ops->view              = PCView_ASM;
    865:   pc->ops->applyrichardson   = 0;

    867:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMSetLocalSubdomains_C","PCASMSetLocalSubdomains_ASM",
    868:                     PCASMSetLocalSubdomains_ASM);
    869:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMSetTotalSubdomains_C","PCASMSetTotalSubdomains_ASM",
    870:                     PCASMSetTotalSubdomains_ASM);
    871:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMSetOverlap_C","PCASMSetOverlap_ASM",
    872:                     PCASMSetOverlap_ASM);
    873:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMSetType_C","PCASMSetType_ASM",
    874:                     PCASMSetType_ASM);
    875:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMGetSubKSP_C","PCASMGetSubKSP_ASM",
    876:                     PCASMGetSubKSP_ASM);
    877: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMSetUseInPlace_C","PCASMSetUseInPlace_ASM",
    878:                     PCASMSetUseInPlace_ASM);
    879:   return(0);
    880: }

    886: /*@
    887:    PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
    888:    preconditioner for a two-dimensional problem on a regular grid.

    890:    Not Collective

    892:    Input Parameters:
    893: +  m, n - the number of mesh points in the x and y directions
    894: .  M, N - the number of subdomains in the x and y directions
    895: .  dof - degrees of freedom per node
    896: -  overlap - overlap in mesh lines

    898:    Output Parameters:
    899: +  Nsub - the number of subdomains created
    900: -  is - the array of index sets defining the subdomains

    902:    Note:
    903:    Presently PCAMSCreateSubdomains2d() is valid only for sequential
    904:    preconditioners.  More general related routines are
    905:    PCASMSetTotalSubdomains() and PCASMSetLocalSubdomains().

    907:    Level: advanced

    909: .keywords: PC, ASM, additive Schwarz, create, subdomains, 2D, regular grid

    911: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
    912:           PCASMSetOverlap()
    913: @*/
    914: PetscErrorCode PETSCKSP_DLLEXPORT PCASMCreateSubdomains2D(PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt dof,PetscInt overlap,PetscInt *Nsub,IS **is)
    915: {
    916:   PetscInt       i,j,height,width,ystart,xstart,yleft,yright,xleft,xright,loc_outter;
    918:   PetscInt       nidx,*idx,loc,ii,jj,count;

    921:   if (dof != 1) SETERRQ(PETSC_ERR_SUP," ");

    923:   *Nsub = N*M;
    924:   PetscMalloc((*Nsub)*sizeof(IS **),is);
    925:   ystart = 0;
    926:   loc_outter = 0;
    927:   for (i=0; i<N; i++) {
    928:     height = n/N + ((n % N) > i); /* height of subdomain */
    929:     if (height < 2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Too many N subdomains for mesh dimension n");
    930:     yleft  = ystart - overlap; if (yleft < 0) yleft = 0;
    931:     yright = ystart + height + overlap; if (yright > n) yright = n;
    932:     xstart = 0;
    933:     for (j=0; j<M; j++) {
    934:       width = m/M + ((m % M) > j); /* width of subdomain */
    935:       if (width < 2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Too many M subdomains for mesh dimension m");
    936:       xleft  = xstart - overlap; if (xleft < 0) xleft = 0;
    937:       xright = xstart + width + overlap; if (xright > m) xright = m;
    938:       nidx   = (xright - xleft)*(yright - yleft);
    939:       PetscMalloc(nidx*sizeof(PetscInt),&idx);
    940:       loc    = 0;
    941:       for (ii=yleft; ii<yright; ii++) {
    942:         count = m*ii + xleft;
    943:         for (jj=xleft; jj<xright; jj++) {
    944:           idx[loc++] = count++;
    945:         }
    946:       }
    947:       ISCreateGeneral(PETSC_COMM_SELF,nidx,idx,(*is)+loc_outter++);
    948:       PetscFree(idx);
    949:       xstart += width;
    950:     }
    951:     ystart += height;
    952:   }
    953:   for (i=0; i<*Nsub; i++) { ISSort((*is)[i]); }
    954:   return(0);
    955: }

    959: /*@C
    960:     PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
    961:     only) for the additive Schwarz preconditioner.

    963:     Collective on PC

    965:     Input Parameter:
    966: .   pc - the preconditioner context

    968:     Output Parameters:
    969: +   n - the number of subdomains for this processor (default value = 1)
    970: -   is - the index sets that define the subdomains for this processor
    971:

    973:     Notes:
    974:     The IS numbering is in the parallel, global numbering of the vector.

    976:     Level: advanced

    978: .keywords: PC, ASM, set, local, subdomains, additive Schwarz

    980: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
    981:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubmatrices()
    982: @*/
    983: PetscErrorCode PETSCKSP_DLLEXPORT PCASMGetLocalSubdomains(PC pc,PetscInt *n,IS *is[])
    984: {
    985:   PC_ASM *osm;

    990:   if (!pc->setupcalled) {
    991:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUP() or PCSetUp().");
    992:   }

    994:   osm = (PC_ASM*)pc->data;
    995:   if (n)  *n = osm->n_local_true;
    996:   if (is) *is = osm->is;
    997:   return(0);
    998: }

    1002: /*@C
    1003:     PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
    1004:     only) for the additive Schwarz preconditioner.

    1006:     Collective on PC

    1008:     Input Parameter:
    1009: .   pc - the preconditioner context

    1011:     Output Parameters:
    1012: +   n - the number of matrices for this processor (default value = 1)
    1013: -   mat - the matrices
    1014:

    1016:     Level: advanced

    1018: .keywords: PC, ASM, set, local, subdomains, additive Schwarz, block Jacobi

    1020: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
    1021:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
    1022: @*/
    1023: PetscErrorCode PETSCKSP_DLLEXPORT PCASMGetLocalSubmatrices(PC pc,PetscInt *n,Mat *mat[])
    1024: {
    1025:   PC_ASM *osm;

    1030:   if (!pc->setupcalled) {
    1031:     SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUP() or PCSetUp().");
    1032:   }

    1034:   osm = (PC_ASM*)pc->data;
    1035:   if (n)   *n   = osm->n_local_true;
    1036:   if (mat) *mat = osm->pmat;
    1037:   return(0);
    1038: }