Actual source code: asm.c

petsc-master 2019-05-18
Report Typos and Errors
  1: /*
  2:   This file defines an additive Schwarz preconditioner for any Mat implementation.

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

  8:        n - total number of true subdomains on all processors
  9:        n_local_true - actual number of subdomains on this processor
 10:        n_local = maximum over all processors of n_local_true
 11: */
 12:  #include <petsc/private/pcimpl.h>
 13:  #include <petscdm.h>

 15: typedef struct {
 16:   PetscInt   n, n_local, n_local_true;
 17:   PetscInt   overlap;             /* overlap requested by user */
 18:   KSP        *ksp;                /* linear solvers for each block */
 19:   VecScatter restriction;         /* mapping from global to overlapping (process) subdomain*/
 20:   VecScatter *lrestriction;       /* mapping from subregion to overlapping (process) subdomain */
 21:   VecScatter *lprolongation;      /* mapping from non-overlapping subregion to overlapping (process) subdomain; used for restrict additive version of algorithms */
 22:   Vec        lx, ly;              /* work vectors */
 23:   Vec        *x,*y;               /* work vectors */
 24:   IS         lis;                 /* index set that defines each overlapping multiplicative (process) subdomain */
 25:   IS         *is;                 /* index set that defines each overlapping subdomain */
 26:   IS         *is_local;           /* index set that defines each non-overlapping subdomain, may be NULL */
 27:   Mat        *mat,*pmat;          /* mat is not currently used */
 28:   PCASMType  type;                /* use reduced interpolation, restriction or both */
 29:   PetscBool  type_set;            /* if user set this value (so won't change it for symmetric problems) */
 30:   PetscBool  same_local_solves;   /* flag indicating whether all local solvers are same */
 31:   PetscBool  sort_indices;        /* flag to sort subdomain indices */
 32:   PetscBool  dm_subdomains;       /* whether DM is allowed to define subdomains */
 33:   PCCompositeType loctype;        /* the type of composition for local solves */
 34:   MatType    sub_mat_type;        /* the type of Mat used for subdomain solves (can be MATSAME or NULL) */
 35:   /* For multiplicative solve */
 36:   Mat       *lmats;               /* submatrices for overlapping multiplicative (process) subdomain */
 37: } PC_ASM;

 39: static PetscErrorCode PCView_ASM(PC pc,PetscViewer viewer)
 40: {
 41:   PC_ASM         *osm = (PC_ASM*)pc->data;
 43:   PetscMPIInt    rank;
 44:   PetscInt       i,bsz;
 45:   PetscBool      iascii,isstring;
 46:   PetscViewer    sviewer;

 49:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
 50:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
 51:   if (iascii) {
 52:     char overlaps[256] = "user-defined overlap",blocks[256] = "total subdomain blocks not yet set";
 53:     if (osm->overlap >= 0) {PetscSNPrintf(overlaps,sizeof(overlaps),"amount of overlap = %D",osm->overlap);}
 54:     if (osm->n > 0) {PetscSNPrintf(blocks,sizeof(blocks),"total subdomain blocks = %D",osm->n);}
 55:     PetscViewerASCIIPrintf(viewer,"  %s, %s\n",blocks,overlaps);
 56:     PetscViewerASCIIPrintf(viewer,"  restriction/interpolation type - %s\n",PCASMTypes[osm->type]);
 57:     if (osm->dm_subdomains) {PetscViewerASCIIPrintf(viewer,"  Additive Schwarz: using DM to define subdomains\n");}
 58:     if (osm->loctype != PC_COMPOSITE_ADDITIVE) {PetscViewerASCIIPrintf(viewer,"  Additive Schwarz: local solve composition type - %s\n",PCCompositeTypes[osm->loctype]);}
 59:     MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
 60:     if (osm->same_local_solves) {
 61:       if (osm->ksp) {
 62:         PetscViewerASCIIPrintf(viewer,"  Local solve is same for all blocks, in the following KSP and PC objects:\n");
 63:         PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
 64:         if (!rank) {
 65:           PetscViewerASCIIPushTab(viewer);
 66:           KSPView(osm->ksp[0],sviewer);
 67:           PetscViewerASCIIPopTab(viewer);
 68:         }
 69:         PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
 70:       }
 71:     } else {
 72:       PetscViewerASCIIPushSynchronized(viewer);
 73:       PetscViewerASCIISynchronizedPrintf(viewer,"  [%d] number of local blocks = %D\n",(int)rank,osm->n_local_true);
 74:       PetscViewerFlush(viewer);
 75:       PetscViewerASCIIPrintf(viewer,"  Local solve info for each block is in the following KSP and PC objects:\n");
 76:       PetscViewerASCIIPushTab(viewer);
 77:       PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
 78:       PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
 79:       for (i=0; i<osm->n_local_true; i++) {
 80:         ISGetLocalSize(osm->is[i],&bsz);
 81:         PetscViewerASCIISynchronizedPrintf(sviewer,"[%d] local block number %D, size = %D\n",(int)rank,i,bsz);
 82:         KSPView(osm->ksp[i],sviewer);
 83:         PetscViewerASCIISynchronizedPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
 84:       }
 85:       PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
 86:       PetscViewerASCIIPopTab(viewer);
 87:       PetscViewerFlush(viewer);
 88:       PetscViewerASCIIPopSynchronized(viewer);
 89:     }
 90:   } else if (isstring) {
 91:     PetscViewerStringSPrintf(viewer," blocks=%D, overlap=%D, type=%s",osm->n,osm->overlap,PCASMTypes[osm->type]);
 92:     PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
 93:     if (osm->ksp) {KSPView(osm->ksp[0],sviewer);}
 94:     PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
 95:   }
 96:   return(0);
 97: }

 99: static PetscErrorCode PCASMPrintSubdomains(PC pc)
100: {
101:   PC_ASM         *osm = (PC_ASM*)pc->data;
102:   const char     *prefix;
103:   char           fname[PETSC_MAX_PATH_LEN+1];
104:   PetscViewer    viewer, sviewer;
105:   char           *s;
106:   PetscInt       i,j,nidx;
107:   const PetscInt *idx;
108:   PetscMPIInt    rank, size;

112:   MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
113:   MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
114:   PCGetOptionsPrefix(pc,&prefix);
115:   PetscOptionsGetString(NULL,prefix,"-pc_asm_print_subdomains",fname,PETSC_MAX_PATH_LEN,NULL);
116:   if (fname[0] == 0) { PetscStrcpy(fname,"stdout"); };
117:   PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),fname,&viewer);
118:   for (i=0; i<osm->n_local; i++) {
119:     if (i < osm->n_local_true) {
120:       ISGetLocalSize(osm->is[i],&nidx);
121:       ISGetIndices(osm->is[i],&idx);
122:       /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
123: #define len  16*(nidx+1)+512
124:       PetscMalloc1(len,&s);
125:       PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer);
126: #undef len
127:       PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D with overlap:\n", rank, size, i);
128:       for (j=0; j<nidx; j++) {
129:         PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
130:       }
131:       ISRestoreIndices(osm->is[i],&idx);
132:       PetscViewerStringSPrintf(sviewer,"\n");
133:       PetscViewerDestroy(&sviewer);
134:       PetscViewerASCIIPushSynchronized(viewer);
135:       PetscViewerASCIISynchronizedPrintf(viewer, s);
136:       PetscViewerFlush(viewer);
137:       PetscViewerASCIIPopSynchronized(viewer);
138:       PetscFree(s);
139:       if (osm->is_local) {
140:         /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
141: #define len  16*(nidx+1)+512
142:         PetscMalloc1(len, &s);
143:         PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer);
144: #undef len
145:         PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D without overlap:\n", rank, size, i);
146:         ISGetLocalSize(osm->is_local[i],&nidx);
147:         ISGetIndices(osm->is_local[i],&idx);
148:         for (j=0; j<nidx; j++) {
149:           PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
150:         }
151:         ISRestoreIndices(osm->is_local[i],&idx);
152:         PetscViewerStringSPrintf(sviewer,"\n");
153:         PetscViewerDestroy(&sviewer);
154:         PetscViewerASCIIPushSynchronized(viewer);
155:         PetscViewerASCIISynchronizedPrintf(viewer, s);
156:         PetscViewerFlush(viewer);
157:         PetscViewerASCIIPopSynchronized(viewer);
158:         PetscFree(s);
159:       }
160:     } else {
161:       /* Participate in collective viewer calls. */
162:       PetscViewerASCIIPushSynchronized(viewer);
163:       PetscViewerFlush(viewer);
164:       PetscViewerASCIIPopSynchronized(viewer);
165:       /* Assume either all ranks have is_local or none do. */
166:       if (osm->is_local) {
167:         PetscViewerASCIIPushSynchronized(viewer);
168:         PetscViewerFlush(viewer);
169:         PetscViewerASCIIPopSynchronized(viewer);
170:       }
171:     }
172:   }
173:   PetscViewerFlush(viewer);
174:   PetscViewerDestroy(&viewer);
175:   return(0);
176: }

178: static PetscErrorCode PCSetUp_ASM(PC pc)
179: {
180:   PC_ASM         *osm = (PC_ASM*)pc->data;
182:   PetscBool      symset,flg;
183:   PetscInt       i,m,m_local;
184:   MatReuse       scall = MAT_REUSE_MATRIX;
185:   IS             isl;
186:   KSP            ksp;
187:   PC             subpc;
188:   const char     *prefix,*pprefix;
189:   Vec            vec;
190:   DM             *domain_dm = NULL;

193:   if (!pc->setupcalled) {
194:     PetscInt m;

196:     if (!osm->type_set) {
197:       MatIsSymmetricKnown(pc->pmat,&symset,&flg);
198:       if (symset && flg) osm->type = PC_ASM_BASIC;
199:     }

201:     /* Note: if subdomains have been set either via PCASMSetTotalSubdomains() or via PCASMSetLocalSubdomains(), osm->n_local_true will not be PETSC_DECIDE */
202:     if (osm->n_local_true == PETSC_DECIDE) {
203:       /* no subdomains given */
204:       /* try pc->dm first, if allowed */
205:       if (osm->dm_subdomains && pc->dm) {
206:         PetscInt  num_domains, d;
207:         char      **domain_names;
208:         IS        *inner_domain_is, *outer_domain_is;
209:         DMCreateDomainDecomposition(pc->dm, &num_domains, &domain_names, &inner_domain_is, &outer_domain_is, &domain_dm);
210:         osm->overlap = -1; /* We do not want to increase the overlap of the IS.
211:                               A future improvement of this code might allow one to use
212:                               DM-defined subdomains and also increase the overlap,
213:                               but that is not currently supported */
214:         if (num_domains) {
215:           PCASMSetLocalSubdomains(pc, num_domains, outer_domain_is, inner_domain_is);
216:         }
217:         for (d = 0; d < num_domains; ++d) {
218:           if (domain_names)    {PetscFree(domain_names[d]);}
219:           if (inner_domain_is) {ISDestroy(&inner_domain_is[d]);}
220:           if (outer_domain_is) {ISDestroy(&outer_domain_is[d]);}
221:         }
222:         PetscFree(domain_names);
223:         PetscFree(inner_domain_is);
224:         PetscFree(outer_domain_is);
225:       }
226:       if (osm->n_local_true == PETSC_DECIDE) {
227:         /* still no subdomains; use one subdomain per processor */
228:         osm->n_local_true = 1;
229:       }
230:     }
231:     { /* determine the global and max number of subdomains */
232:       struct {PetscInt max,sum;} inwork,outwork;
233:       PetscMPIInt size;

235:       inwork.max   = osm->n_local_true;
236:       inwork.sum   = osm->n_local_true;
237:       MPIU_Allreduce(&inwork,&outwork,1,MPIU_2INT,MPIU_MAXSUM_OP,PetscObjectComm((PetscObject)pc));
238:       osm->n_local = outwork.max;
239:       osm->n       = outwork.sum;

241:       MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
242:       if (outwork.max == 1 && outwork.sum == size) {
243:         /* osm->n_local_true = 1 on all processes, set this option may enable use of optimized MatCreateSubMatrices() implementation */
244:         MatSetOption(pc->pmat,MAT_SUBMAT_SINGLEIS,PETSC_TRUE);
245:       }
246:     }
247:     if (!osm->is) { /* create the index sets */
248:       PCASMCreateSubdomains(pc->pmat,osm->n_local_true,&osm->is);
249:     }
250:     if (osm->n_local_true > 1 && !osm->is_local) {
251:       PetscMalloc1(osm->n_local_true,&osm->is_local);
252:       for (i=0; i<osm->n_local_true; i++) {
253:         if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
254:           ISDuplicate(osm->is[i],&osm->is_local[i]);
255:           ISCopy(osm->is[i],osm->is_local[i]);
256:         } else {
257:           PetscObjectReference((PetscObject)osm->is[i]);
258:           osm->is_local[i] = osm->is[i];
259:         }
260:       }
261:     }
262:     PCGetOptionsPrefix(pc,&prefix);
263:     flg  = PETSC_FALSE;
264:     PetscOptionsGetBool(NULL,prefix,"-pc_asm_print_subdomains",&flg,NULL);
265:     if (flg) { PCASMPrintSubdomains(pc); }

267:     if (osm->overlap > 0) {
268:       /* Extend the "overlapping" regions by a number of steps */
269:       MatIncreaseOverlap(pc->pmat,osm->n_local_true,osm->is,osm->overlap);
270:     }
271:     if (osm->sort_indices) {
272:       for (i=0; i<osm->n_local_true; i++) {
273:         ISSort(osm->is[i]);
274:         if (osm->is_local) {
275:           ISSort(osm->is_local[i]);
276:         }
277:       }
278:     }

280:     if (!osm->ksp) {
281:       /* Create the local solvers */
282:       PetscMalloc1(osm->n_local_true,&osm->ksp);
283:       if (domain_dm) {
284:         PetscInfo(pc,"Setting up ASM subproblems using the embedded DM\n");
285:       }
286:       for (i=0; i<osm->n_local_true; i++) {
287:         KSPCreate(PETSC_COMM_SELF,&ksp);
288:         KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
289:         PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
290:         PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
291:         KSPSetType(ksp,KSPPREONLY);
292:         KSPGetPC(ksp,&subpc);
293:         PCGetOptionsPrefix(pc,&prefix);
294:         KSPSetOptionsPrefix(ksp,prefix);
295:         KSPAppendOptionsPrefix(ksp,"sub_");
296:         if (domain_dm) {
297:           KSPSetDM(ksp, domain_dm[i]);
298:           KSPSetDMActive(ksp, PETSC_FALSE);
299:           DMDestroy(&domain_dm[i]);
300:         }
301:         osm->ksp[i] = ksp;
302:       }
303:       if (domain_dm) {
304:         PetscFree(domain_dm);
305:       }
306:     }

308:     ISConcatenate(PETSC_COMM_SELF, osm->n_local_true, osm->is, &osm->lis);
309:     ISSortRemoveDups(osm->lis);
310:     ISGetLocalSize(osm->lis, &m);
311:     VecCreateSeq(PETSC_COMM_SELF, m, &osm->lx);
312:     VecDuplicate(osm->lx, &osm->ly);

314:     scall = MAT_INITIAL_MATRIX;
315:   } else {
316:     /*
317:        Destroy the blocks from the previous iteration
318:     */
319:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
320:       MatDestroyMatrices(osm->n_local_true,&osm->pmat);
321:       scall = MAT_INITIAL_MATRIX;
322:     }
323:   }

325:   /*
326:      Extract out the submatrices
327:   */
328:   MatCreateSubMatrices(pc->pmat,osm->n_local_true,osm->is,osm->is,scall,&osm->pmat);
329:   if (scall == MAT_INITIAL_MATRIX) {
330:     PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
331:     for (i=0; i<osm->n_local_true; i++) {
332:       PetscLogObjectParent((PetscObject)pc,(PetscObject)osm->pmat[i]);
333:       PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);
334:     }
335:   }

337:   /* Convert the types of the submatrices (if needbe) */
338:   if (osm->sub_mat_type) {
339:     for (i=0; i<osm->n_local_true; i++) {
340:       MatConvert(osm->pmat[i],osm->sub_mat_type,MAT_INPLACE_MATRIX,&(osm->pmat[i]));
341:     }
342:   }

344:   if(!pc->setupcalled){
345:     /* Create the local work vectors (from the local matrices) and scatter contexts */
346:     MatCreateVecs(pc->pmat,&vec,0);

348:     if (osm->is_local && (osm->type == PC_ASM_INTERPOLATE || osm->type == PC_ASM_NONE )) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot use interpolate or none PCASMType if is_local was provided to PCASMSetLocalSubdomains()");
349:     if (osm->is_local && osm->type == PC_ASM_RESTRICT && osm->loctype == PC_COMPOSITE_ADDITIVE) {
350:       PetscMalloc1(osm->n_local_true,&osm->lprolongation);
351:     }
352:     PetscMalloc1(osm->n_local_true,&osm->lrestriction);
353:     PetscMalloc1(osm->n_local_true,&osm->x);
354:     PetscMalloc1(osm->n_local_true,&osm->y);

356:     ISGetLocalSize(osm->lis,&m);
357:     ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
358:     VecScatterCreate(vec,osm->lis,osm->lx,isl,&osm->restriction);
359:     ISDestroy(&isl);


362:     for (i=0; i<osm->n_local_true; ++i) {
363:       ISLocalToGlobalMapping ltog;
364:       IS                     isll;
365:       const PetscInt         *idx_is;
366:       PetscInt               *idx_lis,nout;

368:       ISGetLocalSize(osm->is[i],&m);
369:       MatCreateVecs(osm->pmat[i],&osm->x[i],NULL);
370:       VecDuplicate(osm->x[i],&osm->y[i]);

372:       /* generate a scatter from ly to y[i] picking all the overlapping is[i] entries */
373:       ISLocalToGlobalMappingCreateIS(osm->lis,&ltog);
374:       ISGetLocalSize(osm->is[i],&m);
375:       ISGetIndices(osm->is[i], &idx_is);
376:       PetscMalloc1(m,&idx_lis);
377:       ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m,idx_is,&nout,idx_lis);
378:       if (nout != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is not a subset of lis");
379:       ISRestoreIndices(osm->is[i], &idx_is);
380:       ISCreateGeneral(PETSC_COMM_SELF,m,idx_lis,PETSC_OWN_POINTER,&isll);
381:       ISLocalToGlobalMappingDestroy(&ltog);
382:       ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
383:       VecScatterCreate(osm->ly,isll,osm->y[i],isl,&osm->lrestriction[i]);
384:       ISDestroy(&isll);
385:       ISDestroy(&isl);
386:       if (osm->lprolongation) { /* generate a scatter from y[i] to ly picking only the the non-overalapping is_local[i] entries */
387:         ISLocalToGlobalMapping ltog;
388:         IS                     isll,isll_local;
389:         const PetscInt         *idx_local;
390:         PetscInt               *idx1, *idx2, nout;

392:         ISGetLocalSize(osm->is_local[i],&m_local);
393:         ISGetIndices(osm->is_local[i], &idx_local);

395:         ISLocalToGlobalMappingCreateIS(osm->is[i],&ltog);
396:         PetscMalloc1(m_local,&idx1);
397:         ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx1);
398:         ISLocalToGlobalMappingDestroy(&ltog);
399:         if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of is");
400:         ISCreateGeneral(PETSC_COMM_SELF,m_local,idx1,PETSC_OWN_POINTER,&isll);

402:         ISLocalToGlobalMappingCreateIS(osm->lis,&ltog);
403:         PetscMalloc1(m_local,&idx2);
404:         ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx2);
405:         ISLocalToGlobalMappingDestroy(&ltog);
406:         if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of lis");
407:         ISCreateGeneral(PETSC_COMM_SELF,m_local,idx2,PETSC_OWN_POINTER,&isll_local);

409:         ISRestoreIndices(osm->is_local[i], &idx_local);
410:         VecScatterCreate(osm->y[i],isll,osm->ly,isll_local,&osm->lprolongation[i]);

412:         ISDestroy(&isll);
413:         ISDestroy(&isll_local);
414:       }
415:     }
416:     VecDestroy(&vec);
417:   }

419:   if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
420:     IS      *cis;
421:     PetscInt c;

423:     PetscMalloc1(osm->n_local_true, &cis);
424:     for (c = 0; c < osm->n_local_true; ++c) cis[c] = osm->lis;
425:     MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, cis, scall, &osm->lmats);
426:     PetscFree(cis);
427:   }

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

433:   /*
434:      Loop over subdomains putting them into local ksp
435:   */
436:   for (i=0; i<osm->n_local_true; i++) {
437:     KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i]);
438:     if (!pc->setupcalled) {
439:       KSPSetFromOptions(osm->ksp[i]);
440:     }
441:   }
442:   return(0);
443: }

445: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
446: {
447:   PC_ASM             *osm = (PC_ASM*)pc->data;
448:   PetscErrorCode     ierr;
449:   PetscInt           i;
450:   KSPConvergedReason reason;

453:   for (i=0; i<osm->n_local_true; i++) {
454:     KSPSetUp(osm->ksp[i]);
455:     KSPGetConvergedReason(osm->ksp[i],&reason);
456:     if (reason == KSP_DIVERGED_PC_FAILED) {
457:       pc->failedreason = PC_SUBPC_ERROR;
458:     }
459:   }
460:   return(0);
461: }

463: static PetscErrorCode PCApply_ASM(PC pc,Vec x,Vec y)
464: {
465:   PC_ASM         *osm = (PC_ASM*)pc->data;
467:   PetscInt       i,n_local_true = osm->n_local_true;
468:   ScatterMode    forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;

471:   /*
472:      Support for limiting the restriction or interpolation to only local
473:      subdomain values (leaving the other values 0).
474:   */
475:   if (!(osm->type & PC_ASM_RESTRICT)) {
476:     forward = SCATTER_FORWARD_LOCAL;
477:     /* have to zero the work RHS since scatter may leave some slots empty */
478:     VecSet(osm->lx, 0.0);
479:   }
480:   if (!(osm->type & PC_ASM_INTERPOLATE)) {
481:     reverse = SCATTER_REVERSE_LOCAL;
482:   }

484:   if(osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE){
485:     /* zero the global and the local solutions */
486:     VecZeroEntries(y);
487:     VecSet(osm->ly, 0.0);

489:     /* Copy the global RHS to local RHS including the ghost nodes */
490:     VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
491:     VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);

493:     /* Restrict local RHS to the overlapping 0-block RHS */
494:     VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
495:     VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0],  INSERT_VALUES, forward);

497:     /* do the local solves */
498:     for (i = 0; i < n_local_true; ++i) {

500:       /* solve the overlapping i-block */
501:       PetscLogEventBegin(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
502:       KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]);
503:       KSPCheckSolve(osm->ksp[i],pc,osm->y[i]);
504:       PetscLogEventEnd(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);

506:       if (osm->lprolongation) { /* interpolate the non-overalapping i-block solution to the local solution (only for restrictive additive) */
507:         VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
508:         VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
509:       }
510:       else{ /* interpolate the overalapping i-block solution to the local solution */
511:         VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
512:         VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
513:       }

515:       if (i < n_local_true-1) {
516:         /* Restrict local RHS to the overlapping (i+1)-block RHS */
517:         VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
518:         VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);

520:         if ( osm->loctype == PC_COMPOSITE_MULTIPLICATIVE){
521:           /* udpdate the overlapping (i+1)-block RHS using the current local solution */
522:           MatMult(osm->lmats[i+1], osm->ly, osm->y[i+1]);
523:           VecAXPBY(osm->x[i+1],-1.,1., osm->y[i+1]);
524:         }
525:       }
526:     }
527:     /* Add the local solution to the global solution including the ghost nodes */
528:     VecScatterBegin(osm->restriction, osm->ly, y,  ADD_VALUES, reverse);
529:     VecScatterEnd(osm->restriction,  osm->ly, y, ADD_VALUES, reverse);
530:   }else{
531:     SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
532:   }
533:   return(0);
534: }

536: static PetscErrorCode PCApplyTranspose_ASM(PC pc,Vec x,Vec y)
537: {
538:   PC_ASM         *osm = (PC_ASM*)pc->data;
540:   PetscInt       i,n_local_true = osm->n_local_true;
541:   ScatterMode    forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;

544:   /*
545:      Support for limiting the restriction or interpolation to only local
546:      subdomain values (leaving the other values 0).

548:      Note: these are reversed from the PCApply_ASM() because we are applying the
549:      transpose of the three terms
550:   */

552:   if (!(osm->type & PC_ASM_INTERPOLATE)) {
553:     forward = SCATTER_FORWARD_LOCAL;
554:     /* have to zero the work RHS since scatter may leave some slots empty */
555:     VecSet(osm->lx, 0.0);
556:   }
557:   if (!(osm->type & PC_ASM_RESTRICT)) reverse = SCATTER_REVERSE_LOCAL;

559:   /* zero the global and the local solutions */
560:   VecZeroEntries(y);
561:   VecSet(osm->ly, 0.0);

563:   /* Copy the global RHS to local RHS including the ghost nodes */
564:   VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
565:   VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);

567:   /* Restrict local RHS to the overlapping 0-block RHS */
568:   VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
569:   VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0],  INSERT_VALUES, forward);

571:   /* do the local solves */
572:   for (i = 0; i < n_local_true; ++i) {

574:     /* solve the overlapping i-block */
575:     PetscLogEventBegin(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
576:     KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i]);
577:     KSPCheckSolve(osm->ksp[i],pc,osm->y[i]);
578:     PetscLogEventEnd(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);

580:     if (osm->lprolongation) { /* interpolate the non-overalapping i-block solution to the local solution */
581:      VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
582:      VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
583:     }
584:     else{ /* interpolate the overalapping i-block solution to the local solution */
585:       VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
586:       VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
587:     }

589:     if (i < n_local_true-1) {
590:       /* Restrict local RHS to the overlapping (i+1)-block RHS */
591:       VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
592:       VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
593:     }
594:   }
595:   /* Add the local solution to the global solution including the ghost nodes */
596:   VecScatterBegin(osm->restriction, osm->ly, y,  ADD_VALUES, reverse);
597:   VecScatterEnd(osm->restriction,  osm->ly, y, ADD_VALUES, reverse);

599:   return(0);

601: }

603: static PetscErrorCode PCReset_ASM(PC pc)
604: {
605:   PC_ASM         *osm = (PC_ASM*)pc->data;
607:   PetscInt       i;

610:   if (osm->ksp) {
611:     for (i=0; i<osm->n_local_true; i++) {
612:       KSPReset(osm->ksp[i]);
613:     }
614:   }
615:   if (osm->pmat) {
616:     if (osm->n_local_true > 0) {
617:       MatDestroySubMatrices(osm->n_local_true,&osm->pmat);
618:     }
619:   }
620:   if (osm->lrestriction) {
621:     VecScatterDestroy(&osm->restriction);
622:     for (i=0; i<osm->n_local_true; i++) {
623:       VecScatterDestroy(&osm->lrestriction[i]);
624:       if (osm->lprolongation) {VecScatterDestroy(&osm->lprolongation[i]);}
625:       VecDestroy(&osm->x[i]);
626:       VecDestroy(&osm->y[i]);
627:     }
628:     PetscFree(osm->lrestriction);
629:     if (osm->lprolongation) {PetscFree(osm->lprolongation);}
630:     PetscFree(osm->x);
631:     PetscFree(osm->y);

633:   }
634:   PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
635:   ISDestroy(&osm->lis);
636:   VecDestroy(&osm->lx);
637:   VecDestroy(&osm->ly);
638:   if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
639:     MatDestroyMatrices(osm->n_local_true, &osm->lmats);
640:   }

642:   PetscFree(osm->sub_mat_type);

644:   osm->is       = 0;
645:   osm->is_local = 0;
646:   return(0);
647: }

649: static PetscErrorCode PCDestroy_ASM(PC pc)
650: {
651:   PC_ASM         *osm = (PC_ASM*)pc->data;
653:   PetscInt       i;

656:   PCReset_ASM(pc);
657:   if (osm->ksp) {
658:     for (i=0; i<osm->n_local_true; i++) {
659:       KSPDestroy(&osm->ksp[i]);
660:     }
661:     PetscFree(osm->ksp);
662:   }
663:   PetscFree(pc->data);
664:   return(0);
665: }

667: static PetscErrorCode PCSetFromOptions_ASM(PetscOptionItems *PetscOptionsObject,PC pc)
668: {
669:   PC_ASM         *osm = (PC_ASM*)pc->data;
671:   PetscInt       blocks,ovl;
672:   PetscBool      symset,flg;
673:   PCASMType      asmtype;
674:   PCCompositeType loctype;
675:   char           sub_mat_type[256];

678:   /* set the type to symmetric if matrix is symmetric */
679:   if (!osm->type_set && pc->pmat) {
680:     MatIsSymmetricKnown(pc->pmat,&symset,&flg);
681:     if (symset && flg) osm->type = PC_ASM_BASIC;
682:   }
683:   PetscOptionsHead(PetscOptionsObject,"Additive Schwarz options");
684:   PetscOptionsBool("-pc_asm_dm_subdomains","Use DMCreateDomainDecomposition() to define subdomains","PCASMSetDMSubdomains",osm->dm_subdomains,&osm->dm_subdomains,&flg);
685:   PetscOptionsInt("-pc_asm_blocks","Number of subdomains","PCASMSetTotalSubdomains",osm->n,&blocks,&flg);
686:   if (flg) {
687:     PCASMSetTotalSubdomains(pc,blocks,NULL,NULL);
688:     osm->dm_subdomains = PETSC_FALSE;
689:   }
690:   PetscOptionsInt("-pc_asm_local_blocks","Number of local subdomains","PCASMSetLocalSubdomains",osm->n_local_true,&blocks,&flg);
691:   if (flg) {
692:     PCASMSetLocalSubdomains(pc,blocks,NULL,NULL);
693:     osm->dm_subdomains = PETSC_FALSE;
694:   }
695:   PetscOptionsInt("-pc_asm_overlap","Number of grid points overlap","PCASMSetOverlap",osm->overlap,&ovl,&flg);
696:   if (flg) {
697:     PCASMSetOverlap(pc,ovl);
698:     osm->dm_subdomains = PETSC_FALSE;
699:   }
700:   flg  = PETSC_FALSE;
701:   PetscOptionsEnum("-pc_asm_type","Type of restriction/extension","PCASMSetType",PCASMTypes,(PetscEnum)osm->type,(PetscEnum*)&asmtype,&flg);
702:   if (flg) {PCASMSetType(pc,asmtype); }
703:   flg  = PETSC_FALSE;
704:   PetscOptionsEnum("-pc_asm_local_type","Type of local solver composition","PCASMSetLocalType",PCCompositeTypes,(PetscEnum)osm->loctype,(PetscEnum*)&loctype,&flg);
705:   if (flg) {PCASMSetLocalType(pc,loctype); }
706:   PetscOptionsFList("-pc_asm_sub_mat_type","Subsolve Matrix Type","PCASMSetSubMatType",MatList,NULL,sub_mat_type,256,&flg);
707:   if(flg){
708:     PCASMSetSubMatType(pc,sub_mat_type);
709:   }
710:   PetscOptionsTail();
711:   return(0);
712: }

714: /*------------------------------------------------------------------------------------*/

716: static PetscErrorCode  PCASMSetLocalSubdomains_ASM(PC pc,PetscInt n,IS is[],IS is_local[])
717: {
718:   PC_ASM         *osm = (PC_ASM*)pc->data;
720:   PetscInt       i;

723:   if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Each process must have 1 or more blocks, n = %D",n);
724:   if (pc->setupcalled && (n != osm->n_local_true || is)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetLocalSubdomains() should be called before calling PCSetUp().");

726:   if (!pc->setupcalled) {
727:     if (is) {
728:       for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is[i]);}
729:     }
730:     if (is_local) {
731:       for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is_local[i]);}
732:     }
733:     PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);

735:     osm->n_local_true = n;
736:     osm->is           = 0;
737:     osm->is_local     = 0;
738:     if (is) {
739:       PetscMalloc1(n,&osm->is);
740:       for (i=0; i<n; i++) osm->is[i] = is[i];
741:       /* Flag indicating that the user has set overlapping subdomains so PCASM should not increase their size. */
742:       osm->overlap = -1;
743:     }
744:     if (is_local) {
745:       PetscMalloc1(n,&osm->is_local);
746:       for (i=0; i<n; i++) osm->is_local[i] = is_local[i];
747:       if (!is) {
748:         PetscMalloc1(osm->n_local_true,&osm->is);
749:         for (i=0; i<osm->n_local_true; i++) {
750:           if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
751:             ISDuplicate(osm->is_local[i],&osm->is[i]);
752:             ISCopy(osm->is_local[i],osm->is[i]);
753:           } else {
754:             PetscObjectReference((PetscObject)osm->is_local[i]);
755:             osm->is[i] = osm->is_local[i];
756:           }
757:         }
758:       }
759:     }
760:   }
761:   return(0);
762: }

764: static PetscErrorCode  PCASMSetTotalSubdomains_ASM(PC pc,PetscInt N,IS *is,IS *is_local)
765: {
766:   PC_ASM         *osm = (PC_ASM*)pc->data;
768:   PetscMPIInt    rank,size;
769:   PetscInt       n;

772:   if (N < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of total blocks must be > 0, N = %D",N);
773:   if (is || is_local) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Use PCASMSetLocalSubdomains() to set specific index sets\n\they cannot be set globally yet.");

775:   /*
776:      Split the subdomains equally among all processors
777:   */
778:   MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
779:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
780:   n    = N/size + ((N % size) > rank);
781:   if (!n) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Process %d must have at least one block: total processors %d total blocks %D",(int)rank,(int)size,N);
782:   if (pc->setupcalled && n != osm->n_local_true) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCASMSetTotalSubdomains() should be called before PCSetUp().");
783:   if (!pc->setupcalled) {
784:     PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);

786:     osm->n_local_true = n;
787:     osm->is           = 0;
788:     osm->is_local     = 0;
789:   }
790:   return(0);
791: }

793: static PetscErrorCode  PCASMSetOverlap_ASM(PC pc,PetscInt ovl)
794: {
795:   PC_ASM *osm = (PC_ASM*)pc->data;

798:   if (ovl < 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested");
799:   if (pc->setupcalled && ovl != osm->overlap) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetOverlap() should be called before PCSetUp().");
800:   if (!pc->setupcalled) osm->overlap = ovl;
801:   return(0);
802: }

804: static PetscErrorCode  PCASMSetType_ASM(PC pc,PCASMType type)
805: {
806:   PC_ASM *osm = (PC_ASM*)pc->data;

809:   osm->type     = type;
810:   osm->type_set = PETSC_TRUE;
811:   return(0);
812: }

814: static PetscErrorCode  PCASMGetType_ASM(PC pc,PCASMType *type)
815: {
816:   PC_ASM *osm = (PC_ASM*)pc->data;

819:   *type = osm->type;
820:   return(0);
821: }

823: static PetscErrorCode  PCASMSetLocalType_ASM(PC pc, PCCompositeType type)
824: {
825:   PC_ASM *osm = (PC_ASM *) pc->data;

828:   if (type != PC_COMPOSITE_ADDITIVE && type != PC_COMPOSITE_MULTIPLICATIVE) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Only supports additive or multiplicative as the local type");
829:   osm->loctype = type;
830:   return(0);
831: }

833: static PetscErrorCode  PCASMGetLocalType_ASM(PC pc, PCCompositeType *type)
834: {
835:   PC_ASM *osm = (PC_ASM *) pc->data;

838:   *type = osm->loctype;
839:   return(0);
840: }

842: static PetscErrorCode  PCASMSetSortIndices_ASM(PC pc,PetscBool  doSort)
843: {
844:   PC_ASM *osm = (PC_ASM*)pc->data;

847:   osm->sort_indices = doSort;
848:   return(0);
849: }

851: static PetscErrorCode  PCASMGetSubKSP_ASM(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
852: {
853:   PC_ASM         *osm = (PC_ASM*)pc->data;

857:   if (osm->n_local_true < 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Need to call PCSetUP() on PC (or KSPSetUp() on the outer KSP object) before calling here");

859:   if (n_local) *n_local = osm->n_local_true;
860:   if (first_local) {
861:     MPI_Scan(&osm->n_local_true,first_local,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
862:     *first_local -= osm->n_local_true;
863:   }
864:   if (ksp) {
865:     /* Assume that local solves are now different; not necessarily
866:        true though!  This flag is used only for PCView_ASM() */
867:     *ksp                   = osm->ksp;
868:     osm->same_local_solves = PETSC_FALSE;
869:   }
870:   return(0);
871: }

873: static PetscErrorCode  PCASMGetSubMatType_ASM(PC pc,MatType *sub_mat_type)
874: {
875:   PC_ASM         *osm = (PC_ASM*)pc->data;

880:   *sub_mat_type = osm->sub_mat_type;
881:   return(0);
882: }

884: static PetscErrorCode PCASMSetSubMatType_ASM(PC pc,MatType sub_mat_type)
885: {
886:   PetscErrorCode    ierr;
887:   PC_ASM            *osm = (PC_ASM*)pc->data;

891:   PetscFree(osm->sub_mat_type);
892:   PetscStrallocpy(sub_mat_type,(char**)&osm->sub_mat_type);
893:   return(0);
894: }

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

899:     Collective on PC

901:     Input Parameters:
902: +   pc - the preconditioner context
903: .   n - the number of subdomains for this processor (default value = 1)
904: .   is - the index set that defines the subdomains for this processor
905:          (or NULL for PETSc to determine subdomains)
906: -   is_local - the index sets that define the local part of the subdomains for this processor, not used unless PCASMType is PC_ASM_RESTRICT
907:          (or NULL to not provide these)

909:     Options Database Key:
910:     To set the total number of subdomain blocks rather than specify the
911:     index sets, use the option
912: .    -pc_asm_local_blocks <blks> - Sets local blocks

914:     Notes:
915:     The IS numbering is in the parallel, global numbering of the vector for both is and is_local

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

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

921:     If is_local is provided and PCASMType is PC_ASM_RESTRICT then the solution only over the is_local region is interpolated
922:     back to form the global solution (this is the standard restricted additive Schwarz method)

924:     If the is_local is provided and PCASMType is PC_ASM_INTERPOLATE or PC_ASM_NONE then an error is generated since there is
925:     no code to handle that case.

927:     Level: advanced

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

931: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
932:           PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), PCASMType, PCASMSetType()
933: @*/
934: PetscErrorCode  PCASMSetLocalSubdomains(PC pc,PetscInt n,IS is[],IS is_local[])
935: {

940:   PetscTryMethod(pc,"PCASMSetLocalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,is,is_local));
941:   return(0);
942: }

944: /*@C
945:     PCASMSetTotalSubdomains - Sets the subdomains for all processors for the
946:     additive Schwarz preconditioner.  Either all or no processors in the
947:     PC communicator must call this routine, with the same index sets.

949:     Collective on PC

951:     Input Parameters:
952: +   pc - the preconditioner context
953: .   N  - the number of subdomains for all processors
954: .   is - the index sets that define the subdomains for all processors
955:          (or NULL to ask PETSc to determine the subdomains)
956: -   is_local - the index sets that define the local part of the subdomains for this processor
957:          (or NULL to not provide this information)

959:     Options Database Key:
960:     To set the total number of subdomain blocks rather than specify the
961:     index sets, use the option
962: .    -pc_asm_blocks <blks> - Sets total blocks

964:     Notes:
965:     Currently you cannot use this to set the actual subdomains with the argument is or is_local.

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

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

972:     Use PCASMSetLocalSubdomains() to set local subdomains.

974:     The IS numbering is in the parallel, global numbering of the vector for both is and is_local

976:     Level: advanced

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

980: .seealso: PCASMSetLocalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
981:           PCASMCreateSubdomains2D()
982: @*/
983: PetscErrorCode  PCASMSetTotalSubdomains(PC pc,PetscInt N,IS is[],IS is_local[])
984: {

989:   PetscTryMethod(pc,"PCASMSetTotalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,N,is,is_local));
990:   return(0);
991: }

993: /*@
994:     PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
995:     additive Schwarz preconditioner.  Either all or no processors in the
996:     PC communicator must call this routine.

998:     Logically Collective on PC

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

1004:     Options Database Key:
1005: .   -pc_asm_overlap <ovl> - Sets overlap

1007:     Notes:
1008:     By default the ASM preconditioner uses 1 block per processor.  To use
1009:     multiple blocks per perocessor, see PCASMSetTotalSubdomains() and
1010:     PCASMSetLocalSubdomains() (and the option -pc_asm_blocks <blks>).

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

1020:     The default algorithm used by PETSc to increase overlap is fast, but not scalable,
1021:     use the option -mat_increase_overlap_scalable when the problem and number of processes is large.

1023:     Note that one can define initial index sets with any overlap via
1024:     PCASMSetLocalSubdomains(); the routine
1025:     PCASMSetOverlap() merely allows PETSc to extend that overlap further
1026:     if desired.

1028:     Level: intermediate

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

1032: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1033:           PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), MatIncreaseOverlap()
1034: @*/
1035: PetscErrorCode  PCASMSetOverlap(PC pc,PetscInt ovl)
1036: {

1042:   PetscTryMethod(pc,"PCASMSetOverlap_C",(PC,PetscInt),(pc,ovl));
1043:   return(0);
1044: }

1046: /*@
1047:     PCASMSetType - Sets the type of restriction and interpolation used
1048:     for local problems in the additive Schwarz method.

1050:     Logically Collective on PC

1052:     Input Parameters:
1053: +   pc  - the preconditioner context
1054: -   type - variant of ASM, one of
1055: .vb
1056:       PC_ASM_BASIC       - full interpolation and restriction
1057:       PC_ASM_RESTRICT    - full restriction, local processor interpolation
1058:       PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1059:       PC_ASM_NONE        - local processor restriction and interpolation
1060: .ve

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

1065:     Notes:
1066:     if the is_local arguments are passed to PCASMSetLocalSubdomains() then they are used when PC_ASM_RESTRICT has been selected
1067:     to limit the local processor interpolation

1069:     Level: intermediate

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

1073: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1074:           PCASMCreateSubdomains2D(), PCASMType, PCASMSetLocalType(), PCASMGetLocalType()
1075: @*/
1076: PetscErrorCode  PCASMSetType(PC pc,PCASMType type)
1077: {

1083:   PetscTryMethod(pc,"PCASMSetType_C",(PC,PCASMType),(pc,type));
1084:   return(0);
1085: }

1087: /*@
1088:     PCASMGetType - Gets the type of restriction and interpolation used
1089:     for local problems in the additive Schwarz method.

1091:     Logically Collective on PC

1093:     Input Parameter:
1094: .   pc  - the preconditioner context

1096:     Output Parameter:
1097: .   type - variant of ASM, one of

1099: .vb
1100:       PC_ASM_BASIC       - full interpolation and restriction
1101:       PC_ASM_RESTRICT    - full restriction, local processor interpolation
1102:       PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1103:       PC_ASM_NONE        - local processor restriction and interpolation
1104: .ve

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

1109:     Level: intermediate

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

1113: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1114:           PCASMCreateSubdomains2D(), PCASMType, PCASMSetType(), PCASMSetLocalType(), PCASMGetLocalType()
1115: @*/
1116: PetscErrorCode  PCASMGetType(PC pc,PCASMType *type)
1117: {

1122:   PetscUseMethod(pc,"PCASMGetType_C",(PC,PCASMType*),(pc,type));
1123:   return(0);
1124: }

1126: /*@
1127:   PCASMSetLocalType - Sets the type of composition used for local problems in the additive Schwarz method.

1129:   Logically Collective on PC

1131:   Input Parameters:
1132: + pc  - the preconditioner context
1133: - type - type of composition, one of
1134: .vb
1135:   PC_COMPOSITE_ADDITIVE       - local additive combination
1136:   PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1137: .ve

1139:   Options Database Key:
1140: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type

1142:   Level: intermediate

1144: .seealso: PCASMSetType(), PCASMGetType(), PCASMGetLocalType(), PCASM, PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1145: @*/
1146: PetscErrorCode PCASMSetLocalType(PC pc, PCCompositeType type)
1147: {

1153:   PetscTryMethod(pc, "PCASMSetLocalType_C", (PC, PCCompositeType), (pc, type));
1154:   return(0);
1155: }

1157: /*@
1158:   PCASMGetLocalType - Gets the type of composition used for local problems in the additive Schwarz method.

1160:   Logically Collective on PC

1162:   Input Parameter:
1163: . pc  - the preconditioner context

1165:   Output Parameter:
1166: . type - type of composition, one of
1167: .vb
1168:   PC_COMPOSITE_ADDITIVE       - local additive combination
1169:   PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1170: .ve

1172:   Options Database Key:
1173: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type

1175:   Level: intermediate

1177: .seealso: PCASMSetType(), PCASMGetType(), PCASMSetLocalType(), PCASMCreate(), PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1178: @*/
1179: PetscErrorCode PCASMGetLocalType(PC pc, PCCompositeType *type)
1180: {

1186:   PetscUseMethod(pc, "PCASMGetLocalType_C", (PC, PCCompositeType *), (pc, type));
1187:   return(0);
1188: }

1190: /*@
1191:     PCASMSetSortIndices - Determines whether subdomain indices are sorted.

1193:     Logically Collective on PC

1195:     Input Parameters:
1196: +   pc  - the preconditioner context
1197: -   doSort - sort the subdomain indices

1199:     Level: intermediate

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

1203: .seealso: PCASMSetLocalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1204:           PCASMCreateSubdomains2D()
1205: @*/
1206: PetscErrorCode  PCASMSetSortIndices(PC pc,PetscBool doSort)
1207: {

1213:   PetscTryMethod(pc,"PCASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));
1214:   return(0);
1215: }

1217: /*@C
1218:    PCASMGetSubKSP - Gets the local KSP contexts for all blocks on
1219:    this processor.

1221:    Collective on PC iff first_local is requested

1223:    Input Parameter:
1224: .  pc - the preconditioner context

1226:    Output Parameters:
1227: +  n_local - the number of blocks on this processor or NULL
1228: .  first_local - the global number of the first block on this processor or NULL,
1229:                  all processors must request or all must pass NULL
1230: -  ksp - the array of KSP contexts

1232:    Note:
1233:    After PCASMGetSubKSP() the array of KSPes is not to be freed.

1235:    You must call KSPSetUp() before calling PCASMGetSubKSP().

1237:    Fortran note:
1238:    The output argument 'ksp' must be an array of sufficient length or PETSC_NULL_KSP. The latter can be used to learn the necessary length.

1240:    Level: advanced

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

1244: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap(),
1245:           PCASMCreateSubdomains2D(),
1246: @*/
1247: PetscErrorCode  PCASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
1248: {

1253:   PetscUseMethod(pc,"PCASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
1254:   return(0);
1255: }

1257: /* -------------------------------------------------------------------------------------*/
1258: /*MC
1259:    PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1260:            its own KSP object.

1262:    Options Database Keys:
1263: +  -pc_asm_blocks <blks> - Sets total blocks
1264: .  -pc_asm_overlap <ovl> - Sets overlap
1265: .  -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type, default is restrict
1266: -  -pc_asm_local_type [additive, multiplicative] - Sets ASM type, default is additive

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

1272:    Notes:
1273:     Each processor can have one or more blocks, but a block cannot be shared by more
1274:      than one processor. Use PCGASM for subdomains shared by multiple processes. Defaults to one block per processor.

1276:      To set options on the solvers for each block append -sub_ to all the KSP, and PC
1277:         options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly

1279:      To set the options on the solvers separate for each block call PCASMGetSubKSP()
1280:          and set the options directly on the resulting KSP object (you can access its PC
1281:          with KSPGetPC())

1283:    Level: beginner

1285:    Concepts: additive Schwarz method

1287:     References:
1288: +   1. - M Dryja, OB Widlund, An additive variant of the Schwarz alternating method for the case of many subregions
1289:      Courant Institute, New York University Technical report
1290: -   1. - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations,
1291:     Cambridge University Press.

1293: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
1294:            PCBJACOBI, PCASMGetSubKSP(), PCASMSetLocalSubdomains(), PCASMType, PCASMGetType(), PCASMSetLocalType(), PCASMGetLocalType()
1295:            PCASMSetTotalSubdomains(), PCSetModifySubmatrices(), PCASMSetOverlap(), PCASMSetType(), PCCompositeType

1297: M*/

1299: PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1300: {
1302:   PC_ASM         *osm;

1305:   PetscNewLog(pc,&osm);

1307:   osm->n                 = PETSC_DECIDE;
1308:   osm->n_local           = 0;
1309:   osm->n_local_true      = PETSC_DECIDE;
1310:   osm->overlap           = 1;
1311:   osm->ksp               = 0;
1312:   osm->restriction       = 0;
1313:   osm->lprolongation     = 0;
1314:   osm->lrestriction      = 0;
1315:   osm->x                 = 0;
1316:   osm->y                 = 0;
1317:   osm->is                = 0;
1318:   osm->is_local          = 0;
1319:   osm->mat               = 0;
1320:   osm->pmat              = 0;
1321:   osm->type              = PC_ASM_RESTRICT;
1322:   osm->loctype           = PC_COMPOSITE_ADDITIVE;
1323:   osm->same_local_solves = PETSC_TRUE;
1324:   osm->sort_indices      = PETSC_TRUE;
1325:   osm->dm_subdomains     = PETSC_FALSE;
1326:   osm->sub_mat_type      = NULL;

1328:   pc->data                 = (void*)osm;
1329:   pc->ops->apply           = PCApply_ASM;
1330:   pc->ops->applytranspose  = PCApplyTranspose_ASM;
1331:   pc->ops->setup           = PCSetUp_ASM;
1332:   pc->ops->reset           = PCReset_ASM;
1333:   pc->ops->destroy         = PCDestroy_ASM;
1334:   pc->ops->setfromoptions  = PCSetFromOptions_ASM;
1335:   pc->ops->setuponblocks   = PCSetUpOnBlocks_ASM;
1336:   pc->ops->view            = PCView_ASM;
1337:   pc->ops->applyrichardson = 0;

1339:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",PCASMSetLocalSubdomains_ASM);
1340:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",PCASMSetTotalSubdomains_ASM);
1341:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetOverlap_C",PCASMSetOverlap_ASM);
1342:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetType_C",PCASMSetType_ASM);
1343:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetType_C",PCASMGetType_ASM);
1344:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalType_C",PCASMSetLocalType_ASM);
1345:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetLocalType_C",PCASMGetLocalType_ASM);
1346:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSortIndices_C",PCASMSetSortIndices_ASM);
1347:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubKSP_C",PCASMGetSubKSP_ASM);
1348:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubMatType_C",PCASMGetSubMatType_ASM);
1349:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSubMatType_C",PCASMSetSubMatType_ASM);
1350:   return(0);
1351: }

1353: /*@C
1354:    PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1355:    preconditioner for a any problem on a general grid.

1357:    Collective

1359:    Input Parameters:
1360: +  A - The global matrix operator
1361: -  n - the number of local blocks

1363:    Output Parameters:
1364: .  outis - the array of index sets defining the subdomains

1366:    Level: advanced

1368:    Note: this generates nonoverlapping subdomains; the PCASM will generate the overlap
1369:     from these if you use PCASMSetLocalSubdomains()

1371:     In the Fortran version you must provide the array outis[] already allocated of length n.

1373: .keywords: PC, ASM, additive Schwarz, create, subdomains, unstructured grid

1375: .seealso: PCASMSetLocalSubdomains(), PCASMDestroySubdomains()
1376: @*/
1377: PetscErrorCode  PCASMCreateSubdomains(Mat A, PetscInt n, IS* outis[])
1378: {
1379:   MatPartitioning mpart;
1380:   const char      *prefix;
1381:   PetscInt        i,j,rstart,rend,bs;
1382:   PetscBool       hasop, isbaij = PETSC_FALSE,foundpart = PETSC_FALSE;
1383:   Mat             Ad     = NULL, adj;
1384:   IS              ispart,isnumb,*is;
1385:   PetscErrorCode  ierr;

1390:   if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of local blocks must be > 0, n = %D",n);

1392:   /* Get prefix, row distribution, and block size */
1393:   MatGetOptionsPrefix(A,&prefix);
1394:   MatGetOwnershipRange(A,&rstart,&rend);
1395:   MatGetBlockSize(A,&bs);
1396:   if (rstart/bs*bs != rstart || rend/bs*bs != rend) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"bad row distribution [%D,%D) for matrix block size %D",rstart,rend,bs);

1398:   /* Get diagonal block from matrix if possible */
1399:   MatHasOperation(A,MATOP_GET_DIAGONAL_BLOCK,&hasop);
1400:   if (hasop) {
1401:     MatGetDiagonalBlock(A,&Ad);
1402:   }
1403:   if (Ad) {
1404:     PetscObjectTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij);
1405:     if (!isbaij) {PetscObjectTypeCompare((PetscObject)Ad,MATSEQSBAIJ,&isbaij);}
1406:   }
1407:   if (Ad && n > 1) {
1408:     PetscBool match,done;
1409:     /* Try to setup a good matrix partitioning if available */
1410:     MatPartitioningCreate(PETSC_COMM_SELF,&mpart);
1411:     PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);
1412:     MatPartitioningSetFromOptions(mpart);
1413:     PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGCURRENT,&match);
1414:     if (!match) {
1415:       PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGSQUARE,&match);
1416:     }
1417:     if (!match) { /* assume a "good" partitioner is available */
1418:       PetscInt       na;
1419:       const PetscInt *ia,*ja;
1420:       MatGetRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1421:       if (done) {
1422:         /* Build adjacency matrix by hand. Unfortunately a call to
1423:            MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1424:            remove the block-aij structure and we cannot expect
1425:            MatPartitioning to split vertices as we need */
1426:         PetscInt       i,j,len,nnz,cnt,*iia=0,*jja=0;
1427:         const PetscInt *row;
1428:         nnz = 0;
1429:         for (i=0; i<na; i++) { /* count number of nonzeros */
1430:           len = ia[i+1] - ia[i];
1431:           row = ja + ia[i];
1432:           for (j=0; j<len; j++) {
1433:             if (row[j] == i) { /* don't count diagonal */
1434:               len--; break;
1435:             }
1436:           }
1437:           nnz += len;
1438:         }
1439:         PetscMalloc1(na+1,&iia);
1440:         PetscMalloc1(nnz,&jja);
1441:         nnz    = 0;
1442:         iia[0] = 0;
1443:         for (i=0; i<na; i++) { /* fill adjacency */
1444:           cnt = 0;
1445:           len = ia[i+1] - ia[i];
1446:           row = ja + ia[i];
1447:           for (j=0; j<len; j++) {
1448:             if (row[j] != i) { /* if not diagonal */
1449:               jja[nnz+cnt++] = row[j];
1450:             }
1451:           }
1452:           nnz     += cnt;
1453:           iia[i+1] = nnz;
1454:         }
1455:         /* Partitioning of the adjacency matrix */
1456:         MatCreateMPIAdj(PETSC_COMM_SELF,na,na,iia,jja,NULL,&adj);
1457:         MatPartitioningSetAdjacency(mpart,adj);
1458:         MatPartitioningSetNParts(mpart,n);
1459:         MatPartitioningApply(mpart,&ispart);
1460:         ISPartitioningToNumbering(ispart,&isnumb);
1461:         MatDestroy(&adj);
1462:         foundpart = PETSC_TRUE;
1463:       }
1464:       MatRestoreRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1465:     }
1466:     MatPartitioningDestroy(&mpart);
1467:   }

1469:   PetscMalloc1(n,&is);
1470:   *outis = is;

1472:   if (!foundpart) {

1474:     /* Partitioning by contiguous chunks of rows */

1476:     PetscInt mbs   = (rend-rstart)/bs;
1477:     PetscInt start = rstart;
1478:     for (i=0; i<n; i++) {
1479:       PetscInt count = (mbs/n + ((mbs % n) > i)) * bs;
1480:       ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);
1481:       start += count;
1482:     }

1484:   } else {

1486:     /* Partitioning by adjacency of diagonal block  */

1488:     const PetscInt *numbering;
1489:     PetscInt       *count,nidx,*indices,*newidx,start=0;
1490:     /* Get node count in each partition */
1491:     PetscMalloc1(n,&count);
1492:     ISPartitioningCount(ispart,n,count);
1493:     if (isbaij && bs > 1) { /* adjust for the block-aij case */
1494:       for (i=0; i<n; i++) count[i] *= bs;
1495:     }
1496:     /* Build indices from node numbering */
1497:     ISGetLocalSize(isnumb,&nidx);
1498:     PetscMalloc1(nidx,&indices);
1499:     for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */
1500:     ISGetIndices(isnumb,&numbering);
1501:     PetscSortIntWithPermutation(nidx,numbering,indices);
1502:     ISRestoreIndices(isnumb,&numbering);
1503:     if (isbaij && bs > 1) { /* adjust for the block-aij case */
1504:       PetscMalloc1(nidx*bs,&newidx);
1505:       for (i=0; i<nidx; i++) {
1506:         for (j=0; j<bs; j++) newidx[i*bs+j] = indices[i]*bs + j;
1507:       }
1508:       PetscFree(indices);
1509:       nidx   *= bs;
1510:       indices = newidx;
1511:     }
1512:     /* Shift to get global indices */
1513:     for (i=0; i<nidx; i++) indices[i] += rstart;

1515:     /* Build the index sets for each block */
1516:     for (i=0; i<n; i++) {
1517:       ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);
1518:       ISSort(is[i]);
1519:       start += count[i];
1520:     }

1522:     PetscFree(count);
1523:     PetscFree(indices);
1524:     ISDestroy(&isnumb);
1525:     ISDestroy(&ispart);

1527:   }
1528:   return(0);
1529: }

1531: /*@C
1532:    PCASMDestroySubdomains - Destroys the index sets created with
1533:    PCASMCreateSubdomains(). Should be called after setting subdomains
1534:    with PCASMSetLocalSubdomains().

1536:    Collective

1538:    Input Parameters:
1539: +  n - the number of index sets
1540: .  is - the array of index sets
1541: -  is_local - the array of local index sets, can be NULL

1543:    Level: advanced

1545: .keywords: PC, ASM, additive Schwarz, create, subdomains, unstructured grid

1547: .seealso: PCASMCreateSubdomains(), PCASMSetLocalSubdomains()
1548: @*/
1549: PetscErrorCode  PCASMDestroySubdomains(PetscInt n, IS is[], IS is_local[])
1550: {
1551:   PetscInt       i;

1555:   if (n <= 0) return(0);
1556:   if (is) {
1558:     for (i=0; i<n; i++) { ISDestroy(&is[i]); }
1559:     PetscFree(is);
1560:   }
1561:   if (is_local) {
1563:     for (i=0; i<n; i++) { ISDestroy(&is_local[i]); }
1564:     PetscFree(is_local);
1565:   }
1566:   return(0);
1567: }

1569: /*@
1570:    PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1571:    preconditioner for a two-dimensional problem on a regular grid.

1573:    Not Collective

1575:    Input Parameters:
1576: +  m, n - the number of mesh points in the x and y directions
1577: .  M, N - the number of subdomains in the x and y directions
1578: .  dof - degrees of freedom per node
1579: -  overlap - overlap in mesh lines

1581:    Output Parameters:
1582: +  Nsub - the number of subdomains created
1583: .  is - array of index sets defining overlapping (if overlap > 0) subdomains
1584: -  is_local - array of index sets defining non-overlapping subdomains

1586:    Note:
1587:    Presently PCAMSCreateSubdomains2d() is valid only for sequential
1588:    preconditioners.  More general related routines are
1589:    PCASMSetTotalSubdomains() and PCASMSetLocalSubdomains().

1591:    Level: advanced

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

1595: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1596:           PCASMSetOverlap()
1597: @*/
1598: PetscErrorCode  PCASMCreateSubdomains2D(PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt dof,PetscInt overlap,PetscInt *Nsub,IS **is,IS **is_local)
1599: {
1600:   PetscInt       i,j,height,width,ystart,xstart,yleft,yright,xleft,xright,loc_outer;
1602:   PetscInt       nidx,*idx,loc,ii,jj,count;

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

1607:   *Nsub     = N*M;
1608:   PetscMalloc1(*Nsub,is);
1609:   PetscMalloc1(*Nsub,is_local);
1610:   ystart    = 0;
1611:   loc_outer = 0;
1612:   for (i=0; i<N; i++) {
1613:     height = n/N + ((n % N) > i); /* height of subdomain */
1614:     if (height < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many N subdomains for mesh dimension n");
1615:     yleft  = ystart - overlap; if (yleft < 0) yleft = 0;
1616:     yright = ystart + height + overlap; if (yright > n) yright = n;
1617:     xstart = 0;
1618:     for (j=0; j<M; j++) {
1619:       width = m/M + ((m % M) > j); /* width of subdomain */
1620:       if (width < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many M subdomains for mesh dimension m");
1621:       xleft  = xstart - overlap; if (xleft < 0) xleft = 0;
1622:       xright = xstart + width + overlap; if (xright > m) xright = m;
1623:       nidx   = (xright - xleft)*(yright - yleft);
1624:       PetscMalloc1(nidx,&idx);
1625:       loc    = 0;
1626:       for (ii=yleft; ii<yright; ii++) {
1627:         count = m*ii + xleft;
1628:         for (jj=xleft; jj<xright; jj++) idx[loc++] = count++;
1629:       }
1630:       ISCreateGeneral(PETSC_COMM_SELF,nidx,idx,PETSC_COPY_VALUES,(*is)+loc_outer);
1631:       if (overlap == 0) {
1632:         PetscObjectReference((PetscObject)(*is)[loc_outer]);

1634:         (*is_local)[loc_outer] = (*is)[loc_outer];
1635:       } else {
1636:         for (loc=0,ii=ystart; ii<ystart+height; ii++) {
1637:           for (jj=xstart; jj<xstart+width; jj++) {
1638:             idx[loc++] = m*ii + jj;
1639:           }
1640:         }
1641:         ISCreateGeneral(PETSC_COMM_SELF,loc,idx,PETSC_COPY_VALUES,*is_local+loc_outer);
1642:       }
1643:       PetscFree(idx);
1644:       xstart += width;
1645:       loc_outer++;
1646:     }
1647:     ystart += height;
1648:   }
1649:   for (i=0; i<*Nsub; i++) { ISSort((*is)[i]); }
1650:   return(0);
1651: }

1653: /*@C
1654:     PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
1655:     only) for the additive Schwarz preconditioner.

1657:     Not Collective

1659:     Input Parameter:
1660: .   pc - the preconditioner context

1662:     Output Parameters:
1663: +   n - the number of subdomains for this processor (default value = 1)
1664: .   is - the index sets that define the subdomains for this processor
1665: -   is_local - the index sets that define the local part of the subdomains for this processor (can be NULL)


1668:     Notes:
1669:     The IS numbering is in the parallel, global numbering of the vector.

1671:     Level: advanced

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

1675: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1676:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubmatrices()
1677: @*/
1678: PetscErrorCode  PCASMGetLocalSubdomains(PC pc,PetscInt *n,IS *is[],IS *is_local[])
1679: {
1680:   PC_ASM         *osm = (PC_ASM*)pc->data;
1682:   PetscBool      match;

1688:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1689:   if (!match) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PC is not a PCASM");
1690:   if (n) *n = osm->n_local_true;
1691:   if (is) *is = osm->is;
1692:   if (is_local) *is_local = osm->is_local;
1693:   return(0);
1694: }

1696: /*@C
1697:     PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1698:     only) for the additive Schwarz preconditioner.

1700:     Not Collective

1702:     Input Parameter:
1703: .   pc - the preconditioner context

1705:     Output Parameters:
1706: +   n - the number of matrices for this processor (default value = 1)
1707: -   mat - the matrices

1709:     Level: advanced

1711:     Notes:
1712:     Call after PCSetUp() (or KSPSetUp()) but before PCApply() (or KSPApply()) and before PCSetUpOnBlocks())

1714:            Usually one would use PCSetModifySubmatrices() to change the submatrices in building the preconditioner.

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

1718: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1719:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains(), PCSetModifySubmatrices()
1720: @*/
1721: PetscErrorCode  PCASMGetLocalSubmatrices(PC pc,PetscInt *n,Mat *mat[])
1722: {
1723:   PC_ASM         *osm;
1725:   PetscBool      match;

1731:   if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUP() or PCSetUp().");
1732:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1733:   if (!match) {
1734:     if (n) *n = 0;
1735:     if (mat) *mat = NULL;
1736:   } else {
1737:     osm = (PC_ASM*)pc->data;
1738:     if (n) *n = osm->n_local_true;
1739:     if (mat) *mat = osm->pmat;
1740:   }
1741:   return(0);
1742: }

1744: /*@
1745:     PCASMSetDMSubdomains - Indicates whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.

1747:     Logically Collective

1749:     Input Parameter:
1750: +   pc  - the preconditioner
1751: -   flg - boolean indicating whether to use subdomains defined by the DM

1753:     Options Database Key:
1754: .   -pc_asm_dm_subdomains

1756:     Level: intermediate

1758:     Notes:
1759:     PCASMSetTotalSubdomains() and PCASMSetOverlap() take precedence over PCASMSetDMSubdomains(),
1760:     so setting either of the first two effectively turns the latter off.

1762: .keywords: PC, ASM, DM, set, subdomains, additive Schwarz

1764: .seealso: PCASMGetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1765:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1766: @*/
1767: PetscErrorCode  PCASMSetDMSubdomains(PC pc,PetscBool flg)
1768: {
1769:   PC_ASM         *osm = (PC_ASM*)pc->data;
1771:   PetscBool      match;

1776:   if (pc->setupcalled) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for a setup PC.");
1777:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1778:   if (match) {
1779:     osm->dm_subdomains = flg;
1780:   }
1781:   return(0);
1782: }

1784: /*@
1785:     PCASMGetDMSubdomains - Returns flag indicating whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1786:     Not Collective

1788:     Input Parameter:
1789: .   pc  - the preconditioner

1791:     Output Parameter:
1792: .   flg - boolean indicating whether to use subdomains defined by the DM

1794:     Level: intermediate

1796: .keywords: PC, ASM, DM, set, subdomains, additive Schwarz

1798: .seealso: PCASMSetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1799:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1800: @*/
1801: PetscErrorCode  PCASMGetDMSubdomains(PC pc,PetscBool* flg)
1802: {
1803:   PC_ASM         *osm = (PC_ASM*)pc->data;
1805:   PetscBool      match;

1810:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1811:   if (match) {
1812:     if (flg) *flg = osm->dm_subdomains;
1813:   }
1814:   return(0);
1815: }

1817: /*@
1818:      PCASMGetSubMatType - Gets the matrix type used for ASM subsolves, as a string.

1820:    Not Collective

1822:    Input Parameter:
1823: .  pc - the PC

1825:    Output Parameter:
1826: .  -pc_asm_sub_mat_type - name of matrix type

1828:    Level: advanced

1830: .keywords: PC, PCASM, MatType, set

1832: .seealso: PCASMSetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1833: @*/
1834: PetscErrorCode  PCASMGetSubMatType(PC pc,MatType *sub_mat_type){

1837:   PetscTryMethod(pc,"PCASMGetSubMatType_C",(PC,MatType*),(pc,sub_mat_type));
1838:   return(0);
1839: }

1841: /*@
1842:      PCASMSetSubMatType - Set the type of matrix used for ASM subsolves

1844:    Collective on Mat

1846:    Input Parameters:
1847: +  pc             - the PC object
1848: -  sub_mat_type   - matrix type

1850:    Options Database Key:
1851: .  -pc_asm_sub_mat_type  <sub_mat_type> - Sets the matrix type used for subsolves, for example, seqaijviennacl. If you specify a base name like aijviennacl, the corresponding sequential type is assumed.

1853:    Notes:
1854:    See "${PETSC_DIR}/include/petscmat.h" for available types

1856:   Level: advanced

1858: .keywords: PC, PCASM, MatType, set

1860: .seealso: PCASMGetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1861: @*/
1862: PetscErrorCode PCASMSetSubMatType(PC pc,MatType sub_mat_type)
1863: {

1866:   PetscTryMethod(pc,"PCASMSetSubMatType_C",(PC,MatType),(pc,sub_mat_type));
1867:   return(0);
1868: }