Actual source code: asm.c

petsc-master 2018-07-20
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:       PetscMalloc1(16*(nidx+1)+512, &s);
124:       PetscViewerStringOpen(PETSC_COMM_SELF, s, 16*(nidx+1)+512, &sviewer);
125:       PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D with overlap:\n", rank, size, i);
126:       for (j=0; j<nidx; j++) {
127:         PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
128:       }
129:       ISRestoreIndices(osm->is[i],&idx);
130:       PetscViewerStringSPrintf(sviewer,"\n");
131:       PetscViewerDestroy(&sviewer);
132:       PetscViewerASCIIPushSynchronized(viewer);
133:       PetscViewerASCIISynchronizedPrintf(viewer, s);
134:       PetscViewerFlush(viewer);
135:       PetscViewerASCIIPopSynchronized(viewer);
136:       PetscFree(s);
137:       if (osm->is_local) {
138:         /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
139:         PetscMalloc1(16*(nidx+1)+512, &s);
140:         PetscViewerStringOpen(PETSC_COMM_SELF, s, 16*(nidx+1)+512, &sviewer);
141:         PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D without overlap:\n", rank, size, i);
142:         ISGetLocalSize(osm->is_local[i],&nidx);
143:         ISGetIndices(osm->is_local[i],&idx);
144:         for (j=0; j<nidx; j++) {
145:           PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
146:         }
147:         ISRestoreIndices(osm->is_local[i],&idx);
148:         PetscViewerStringSPrintf(sviewer,"\n");
149:         PetscViewerDestroy(&sviewer);
150:         PetscViewerASCIIPushSynchronized(viewer);
151:         PetscViewerASCIISynchronizedPrintf(viewer, s);
152:         PetscViewerFlush(viewer);
153:         PetscViewerASCIIPopSynchronized(viewer);
154:         PetscFree(s);
155:       }
156:     } else {
157:       /* Participate in collective viewer calls. */
158:       PetscViewerASCIIPushSynchronized(viewer);
159:       PetscViewerFlush(viewer);
160:       PetscViewerASCIIPopSynchronized(viewer);
161:       /* Assume either all ranks have is_local or none do. */
162:       if (osm->is_local) {
163:         PetscViewerASCIIPushSynchronized(viewer);
164:         PetscViewerFlush(viewer);
165:         PetscViewerASCIIPopSynchronized(viewer);
166:       }
167:     }
168:   }
169:   PetscViewerFlush(viewer);
170:   PetscViewerDestroy(&viewer);
171:   return(0);
172: }

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

189:   if (!pc->setupcalled) {
190:     PetscInt m;

192:     if (!osm->type_set) {
193:       MatIsSymmetricKnown(pc->pmat,&symset,&flg);
194:       if (symset && flg) osm->type = PC_ASM_BASIC;
195:     }

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

231:       inwork.max   = osm->n_local_true;
232:       inwork.sum   = osm->n_local_true;
233:       MPIU_Allreduce(&inwork,&outwork,1,MPIU_2INT,MPIU_MAXSUM_OP,PetscObjectComm((PetscObject)pc));
234:       osm->n_local = outwork.max;
235:       osm->n       = outwork.sum;

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

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

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

304:     ISConcatenate(PETSC_COMM_SELF, osm->n_local_true, osm->is, &osm->lis);
305:     ISSortRemoveDups(osm->lis);
306:     ISGetLocalSize(osm->lis, &m);
307:     VecCreateSeq(PETSC_COMM_SELF, m, &osm->lx);
308:     VecDuplicate(osm->lx, &osm->ly);

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

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

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

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

344:     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()");
345:     if (osm->is_local && osm->type == PC_ASM_RESTRICT && osm->loctype == PC_COMPOSITE_ADDITIVE) {
346:       PetscMalloc1(osm->n_local_true,&osm->lprolongation);
347:     }
348:     PetscMalloc1(osm->n_local_true,&osm->lrestriction);
349:     PetscMalloc1(osm->n_local_true,&osm->x);
350:     PetscMalloc1(osm->n_local_true,&osm->y);

352:     ISGetLocalSize(osm->lis,&m);
353:     ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
354:     VecScatterCreate(vec,osm->lis,osm->lx,isl,&osm->restriction);
355:     ISDestroy(&isl);


358:     for (i=0; i<osm->n_local_true; ++i) {
359:       ISLocalToGlobalMapping ltog;
360:       IS                     isll;
361:       const PetscInt         *idx_is;
362:       PetscInt               *idx_lis,nout;

364:       ISGetLocalSize(osm->is[i],&m);
365:       MatCreateVecs(osm->pmat[i],&osm->x[i],NULL);
366:       VecDuplicate(osm->x[i],&osm->y[i]);

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

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

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

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

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

408:         ISDestroy(&isll);
409:         ISDestroy(&isll_local);
410:       }
411:     }
412:     VecDestroy(&vec);
413:   }

415:   if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
416:     IS      *cis;
417:     PetscInt c;

419:     PetscMalloc1(osm->n_local_true, &cis);
420:     for (c = 0; c < osm->n_local_true; ++c) cis[c] = osm->lis;
421:     MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, cis, scall, &osm->lmats);
422:     PetscFree(cis);
423:   }

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

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

441: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
442: {
443:   PC_ASM             *osm = (PC_ASM*)pc->data;
444:   PetscErrorCode     ierr;
445:   PetscInt           i;
446:   KSPConvergedReason reason;

449:   for (i=0; i<osm->n_local_true; i++) {
450:     KSPSetUp(osm->ksp[i]);
451:     KSPGetConvergedReason(osm->ksp[i],&reason);
452:     if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
453:       pc->failedreason = PC_SUBPC_ERROR;
454:     }
455:   }
456:   return(0);
457: }

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

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

480:   if(osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE){
481:     /* zero the global and the local solutions */
482:     VecZeroEntries(y);
483:     VecSet(osm->ly, 0.0);

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

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

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

496:       /* solve the overlapping i-block */
497:       PetscLogEventBegin(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
498:       KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]);
499:       PetscLogEventEnd(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);

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

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

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

531: static PetscErrorCode PCApplyTranspose_ASM(PC pc,Vec x,Vec y)
532: {
533:   PC_ASM         *osm = (PC_ASM*)pc->data;
535:   PetscInt       i,n_local_true = osm->n_local_true;
536:   ScatterMode    forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;

539:   /*
540:      Support for limiting the restriction or interpolation to only local
541:      subdomain values (leaving the other values 0).

543:      Note: these are reversed from the PCApply_ASM() because we are applying the
544:      transpose of the three terms
545:   */

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

554:   /* zero the global and the local solutions */
555:   VecZeroEntries(y);
556:   VecSet(osm->ly, 0.0);

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

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

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

569:     /* solve the overlapping i-block */
570:     PetscLogEventBegin(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
571:     KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]);
572:     PetscLogEventEnd(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);

574:     if (osm->lprolongation) { /* interpolate the non-overalapping i-block solution to the local solution */
575:      VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
576:      VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
577:     }
578:     else{ /* interpolate the overalapping i-block solution to the local solution */
579:       VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
580:       VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
581:     }

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

593:   return(0);

595: }

597: static PetscErrorCode PCReset_ASM(PC pc)
598: {
599:   PC_ASM         *osm = (PC_ASM*)pc->data;
601:   PetscInt       i;

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

627:   }
628:   PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
629:   ISDestroy(&osm->lis);
630:   VecDestroy(&osm->lx);
631:   VecDestroy(&osm->ly);
632:   if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
633:     MatDestroyMatrices(osm->n_local_true, &osm->lmats);
634:   }

636:   PetscFree(osm->sub_mat_type);

638:   osm->is       = 0;
639:   osm->is_local = 0;
640:   return(0);
641: }

643: static PetscErrorCode PCDestroy_ASM(PC pc)
644: {
645:   PC_ASM         *osm = (PC_ASM*)pc->data;
647:   PetscInt       i;

650:   PCReset_ASM(pc);
651:   if (osm->ksp) {
652:     for (i=0; i<osm->n_local_true; i++) {
653:       KSPDestroy(&osm->ksp[i]);
654:     }
655:     PetscFree(osm->ksp);
656:   }
657:   PetscFree(pc->data);
658:   return(0);
659: }

661: static PetscErrorCode PCSetFromOptions_ASM(PetscOptionItems *PetscOptionsObject,PC pc)
662: {
663:   PC_ASM         *osm = (PC_ASM*)pc->data;
665:   PetscInt       blocks,ovl;
666:   PetscBool      symset,flg;
667:   PCASMType      asmtype;
668:   PCCompositeType loctype;
669:   char           sub_mat_type[256];

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

708: /*------------------------------------------------------------------------------------*/

710: static PetscErrorCode  PCASMSetLocalSubdomains_ASM(PC pc,PetscInt n,IS is[],IS is_local[])
711: {
712:   PC_ASM         *osm = (PC_ASM*)pc->data;
714:   PetscInt       i;

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

720:   if (!pc->setupcalled) {
721:     if (is) {
722:       for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is[i]);}
723:     }
724:     if (is_local) {
725:       for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is_local[i]);}
726:     }
727:     PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);

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

758: static PetscErrorCode  PCASMSetTotalSubdomains_ASM(PC pc,PetscInt N,IS *is,IS *is_local)
759: {
760:   PC_ASM         *osm = (PC_ASM*)pc->data;
762:   PetscMPIInt    rank,size;
763:   PetscInt       n;

766:   if (N < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of total blocks must be > 0, N = %D",N);
767:   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.");

769:   /*
770:      Split the subdomains equally among all processors
771:   */
772:   MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
773:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
774:   n    = N/size + ((N % size) > rank);
775:   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);
776:   if (pc->setupcalled && n != osm->n_local_true) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCASMSetTotalSubdomains() should be called before PCSetUp().");
777:   if (!pc->setupcalled) {
778:     PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);

780:     osm->n_local_true = n;
781:     osm->is           = 0;
782:     osm->is_local     = 0;
783:   }
784:   return(0);
785: }

787: static PetscErrorCode  PCASMSetOverlap_ASM(PC pc,PetscInt ovl)
788: {
789:   PC_ASM *osm = (PC_ASM*)pc->data;

792:   if (ovl < 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested");
793:   if (pc->setupcalled && ovl != osm->overlap) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetOverlap() should be called before PCSetUp().");
794:   if (!pc->setupcalled) osm->overlap = ovl;
795:   return(0);
796: }

798: static PetscErrorCode  PCASMSetType_ASM(PC pc,PCASMType type)
799: {
800:   PC_ASM *osm = (PC_ASM*)pc->data;

803:   osm->type     = type;
804:   osm->type_set = PETSC_TRUE;
805:   return(0);
806: }

808: static PetscErrorCode  PCASMGetType_ASM(PC pc,PCASMType *type)
809: {
810:   PC_ASM *osm = (PC_ASM*)pc->data;

813:   *type = osm->type;
814:   return(0);
815: }

817: static PetscErrorCode  PCASMSetLocalType_ASM(PC pc, PCCompositeType type)
818: {
819:   PC_ASM *osm = (PC_ASM *) pc->data;

822:   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");
823:   osm->loctype = type;
824:   return(0);
825: }

827: static PetscErrorCode  PCASMGetLocalType_ASM(PC pc, PCCompositeType *type)
828: {
829:   PC_ASM *osm = (PC_ASM *) pc->data;

832:   *type = osm->loctype;
833:   return(0);
834: }

836: static PetscErrorCode  PCASMSetSortIndices_ASM(PC pc,PetscBool  doSort)
837: {
838:   PC_ASM *osm = (PC_ASM*)pc->data;

841:   osm->sort_indices = doSort;
842:   return(0);
843: }

845: static PetscErrorCode  PCASMGetSubKSP_ASM(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
846: {
847:   PC_ASM         *osm = (PC_ASM*)pc->data;

851:   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");

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

867: static PetscErrorCode  PCASMGetSubMatType_ASM(PC pc,MatType *sub_mat_type)
868: {
869:   PC_ASM         *osm = (PC_ASM*)pc->data;

874:   *sub_mat_type = osm->sub_mat_type;
875:   return(0);
876: }

878: static PetscErrorCode PCASMSetSubMatType_ASM(PC pc,MatType sub_mat_type)
879: {
880:   PetscErrorCode    ierr;
881:   PC_ASM            *osm = (PC_ASM*)pc->data;

885:   PetscFree(osm->sub_mat_type);
886:   PetscStrallocpy(sub_mat_type,(char**)&osm->sub_mat_type);
887:   return(0);
888: }

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

893:     Collective on PC

895:     Input Parameters:
896: +   pc - the preconditioner context
897: .   n - the number of subdomains for this processor (default value = 1)
898: .   is - the index set that defines the subdomains for this processor
899:          (or NULL for PETSc to determine subdomains)
900: -   is_local - the index sets that define the local part of the subdomains for this processor, not used unless PCASMType is PC_ASM_RESTRICT
901:          (or NULL to not provide these)

903:     Options Database Key:
904:     To set the total number of subdomain blocks rather than specify the
905:     index sets, use the option
906: .    -pc_asm_local_blocks <blks> - Sets local blocks

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

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

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

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

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

921:     Level: advanced

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

925: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
926:           PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), PCASMType, PCASMSetType()
927: @*/
928: PetscErrorCode  PCASMSetLocalSubdomains(PC pc,PetscInt n,IS is[],IS is_local[])
929: {

934:   PetscTryMethod(pc,"PCASMSetLocalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,is,is_local));
935:   return(0);
936: }

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

943:     Collective on PC

945:     Input Parameters:
946: +   pc - the preconditioner context
947: .   N  - the number of subdomains for all processors
948: .   is - the index sets that define the subdomains for all processors
949:          (or NULL to ask PETSc to determine the subdomains)
950: -   is_local - the index sets that define the local part of the subdomains for this processor
951:          (or NULL to not provide this information)

953:     Options Database Key:
954:     To set the total number of subdomain blocks rather than specify the
955:     index sets, use the option
956: .    -pc_asm_blocks <blks> - Sets total blocks

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

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

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

966:     Use PCASMSetLocalSubdomains() to set local subdomains.

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

970:     Level: advanced

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

974: .seealso: PCASMSetLocalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
975:           PCASMCreateSubdomains2D()
976: @*/
977: PetscErrorCode  PCASMSetTotalSubdomains(PC pc,PetscInt N,IS is[],IS is_local[])
978: {

983:   PetscTryMethod(pc,"PCASMSetTotalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,N,is,is_local));
984:   return(0);
985: }

987: /*@
988:     PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
989:     additive Schwarz preconditioner.  Either all or no processors in the
990:     PC communicator must call this routine.

992:     Logically Collective on PC

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

998:     Options Database Key:
999: .   -pc_asm_overlap <ovl> - Sets overlap

1001:     Notes:
1002:     By default the ASM preconditioner uses 1 block per processor.  To use
1003:     multiple blocks per perocessor, see PCASMSetTotalSubdomains() and
1004:     PCASMSetLocalSubdomains() (and the option -pc_asm_blocks <blks>).

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

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

1017:     Note that one can define initial index sets with any overlap via
1018:     PCASMSetLocalSubdomains(); the routine
1019:     PCASMSetOverlap() merely allows PETSc to extend that overlap further
1020:     if desired.

1022:     Level: intermediate

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

1026: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1027:           PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), MatIncreaseOverlap()
1028: @*/
1029: PetscErrorCode  PCASMSetOverlap(PC pc,PetscInt ovl)
1030: {

1036:   PetscTryMethod(pc,"PCASMSetOverlap_C",(PC,PetscInt),(pc,ovl));
1037:   return(0);
1038: }

1040: /*@
1041:     PCASMSetType - Sets the type of restriction and interpolation used
1042:     for local problems in the additive Schwarz method.

1044:     Logically Collective on PC

1046:     Input Parameters:
1047: +   pc  - the preconditioner context
1048: -   type - variant of ASM, one of
1049: .vb
1050:       PC_ASM_BASIC       - full interpolation and restriction
1051:       PC_ASM_RESTRICT    - full restriction, local processor interpolation
1052:       PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1053:       PC_ASM_NONE        - local processor restriction and interpolation
1054: .ve

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

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

1063:     Level: intermediate

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

1067: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1068:           PCASMCreateSubdomains2D(), PCASMType, PCASMSetLocalType(), PCASMGetLocalType()
1069: @*/
1070: PetscErrorCode  PCASMSetType(PC pc,PCASMType type)
1071: {

1077:   PetscTryMethod(pc,"PCASMSetType_C",(PC,PCASMType),(pc,type));
1078:   return(0);
1079: }

1081: /*@
1082:     PCASMGetType - Gets the type of restriction and interpolation used
1083:     for local problems in the additive Schwarz method.

1085:     Logically Collective on PC

1087:     Input Parameter:
1088: .   pc  - the preconditioner context

1090:     Output Parameter:
1091: .   type - variant of ASM, one of

1093: .vb
1094:       PC_ASM_BASIC       - full interpolation and restriction
1095:       PC_ASM_RESTRICT    - full restriction, local processor interpolation
1096:       PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1097:       PC_ASM_NONE        - local processor restriction and interpolation
1098: .ve

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

1103:     Level: intermediate

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

1107: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1108:           PCASMCreateSubdomains2D(), PCASMType, PCASMSetType(), PCASMSetLocalType(), PCASMGetLocalType()
1109: @*/
1110: PetscErrorCode  PCASMGetType(PC pc,PCASMType *type)
1111: {

1116:   PetscUseMethod(pc,"PCASMGetType_C",(PC,PCASMType*),(pc,type));
1117:   return(0);
1118: }

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

1123:   Logically Collective on PC

1125:   Input Parameters:
1126: + pc  - the preconditioner context
1127: - type - type of composition, one of
1128: .vb
1129:   PC_COMPOSITE_ADDITIVE       - local additive combination
1130:   PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1131: .ve

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

1136:   Level: intermediate

1138: .seealso: PCASMSetType(), PCASMGetType(), PCASMGetLocalType(), PCASM, PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1139: @*/
1140: PetscErrorCode PCASMSetLocalType(PC pc, PCCompositeType type)
1141: {

1147:   PetscTryMethod(pc, "PCASMSetLocalType_C", (PC, PCCompositeType), (pc, type));
1148:   return(0);
1149: }

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

1154:   Logically Collective on PC

1156:   Input Parameter:
1157: . pc  - the preconditioner context

1159:   Output Parameter:
1160: . type - type of composition, one of
1161: .vb
1162:   PC_COMPOSITE_ADDITIVE       - local additive combination
1163:   PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1164: .ve

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

1169:   Level: intermediate

1171: .seealso: PCASMSetType(), PCASMGetType(), PCASMSetLocalType(), PCASMCreate(), PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1172: @*/
1173: PetscErrorCode PCASMGetLocalType(PC pc, PCCompositeType *type)
1174: {

1180:   PetscUseMethod(pc, "PCASMGetLocalType_C", (PC, PCCompositeType *), (pc, type));
1181:   return(0);
1182: }

1184: /*@
1185:     PCASMSetSortIndices - Determines whether subdomain indices are sorted.

1187:     Logically Collective on PC

1189:     Input Parameters:
1190: +   pc  - the preconditioner context
1191: -   doSort - sort the subdomain indices

1193:     Level: intermediate

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

1197: .seealso: PCASMSetLocalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1198:           PCASMCreateSubdomains2D()
1199: @*/
1200: PetscErrorCode  PCASMSetSortIndices(PC pc,PetscBool doSort)
1201: {

1207:   PetscTryMethod(pc,"PCASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));
1208:   return(0);
1209: }

1211: /*@C
1212:    PCASMGetSubKSP - Gets the local KSP contexts for all blocks on
1213:    this processor.

1215:    Collective on PC iff first_local is requested

1217:    Input Parameter:
1218: .  pc - the preconditioner context

1220:    Output Parameters:
1221: +  n_local - the number of blocks on this processor or NULL
1222: .  first_local - the global number of the first block on this processor or NULL,
1223:                  all processors must request or all must pass NULL
1224: -  ksp - the array of KSP contexts

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

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

1231:    Fortran note:
1232:    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.

1234:    Level: advanced

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

1238: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap(),
1239:           PCASMCreateSubdomains2D(),
1240: @*/
1241: PetscErrorCode  PCASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
1242: {

1247:   PetscUseMethod(pc,"PCASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
1248:   return(0);
1249: }

1251: /* -------------------------------------------------------------------------------------*/
1252: /*MC
1253:    PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1254:            its own KSP object.

1256:    Options Database Keys:
1257: +  -pc_asm_blocks <blks> - Sets total blocks
1258: .  -pc_asm_overlap <ovl> - Sets overlap
1259: .  -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type, default is restrict
1260: -  -pc_asm_local_type [additive, multiplicative] - Sets ASM type, default is additive

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

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

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

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

1277:    Level: beginner

1279:    Concepts: additive Schwarz method

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

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

1291: M*/

1293: PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1294: {
1296:   PC_ASM         *osm;

1299:   PetscNewLog(pc,&osm);

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

1322:   pc->data                 = (void*)osm;
1323:   pc->ops->apply           = PCApply_ASM;
1324:   pc->ops->applytranspose  = PCApplyTranspose_ASM;
1325:   pc->ops->setup           = PCSetUp_ASM;
1326:   pc->ops->reset           = PCReset_ASM;
1327:   pc->ops->destroy         = PCDestroy_ASM;
1328:   pc->ops->setfromoptions  = PCSetFromOptions_ASM;
1329:   pc->ops->setuponblocks   = PCSetUpOnBlocks_ASM;
1330:   pc->ops->view            = PCView_ASM;
1331:   pc->ops->applyrichardson = 0;

1333:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",PCASMSetLocalSubdomains_ASM);
1334:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",PCASMSetTotalSubdomains_ASM);
1335:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetOverlap_C",PCASMSetOverlap_ASM);
1336:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetType_C",PCASMSetType_ASM);
1337:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetType_C",PCASMGetType_ASM);
1338:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalType_C",PCASMSetLocalType_ASM);
1339:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetLocalType_C",PCASMGetLocalType_ASM);
1340:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSortIndices_C",PCASMSetSortIndices_ASM);
1341:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubKSP_C",PCASMGetSubKSP_ASM);
1342:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubMatType_C",PCASMGetSubMatType_ASM);
1343:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSubMatType_C",PCASMSetSubMatType_ASM);
1344:   return(0);
1345: }

1347: /*@C
1348:    PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1349:    preconditioner for a any problem on a general grid.

1351:    Collective

1353:    Input Parameters:
1354: +  A - The global matrix operator
1355: -  n - the number of local blocks

1357:    Output Parameters:
1358: .  outis - the array of index sets defining the subdomains

1360:    Level: advanced

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

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

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

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

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

1386:   /* Get prefix, row distribution, and block size */
1387:   MatGetOptionsPrefix(A,&prefix);
1388:   MatGetOwnershipRange(A,&rstart,&rend);
1389:   MatGetBlockSize(A,&bs);
1390:   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);

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

1463:   PetscMalloc1(n,&is);
1464:   *outis = is;

1466:   if (!foundpart) {

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

1470:     PetscInt mbs   = (rend-rstart)/bs;
1471:     PetscInt start = rstart;
1472:     for (i=0; i<n; i++) {
1473:       PetscInt count = (mbs/n + ((mbs % n) > i)) * bs;
1474:       ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);
1475:       start += count;
1476:     }

1478:   } else {

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

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

1509:     /* Build the index sets for each block */
1510:     for (i=0; i<n; i++) {
1511:       ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);
1512:       ISSort(is[i]);
1513:       start += count[i];
1514:     }

1516:     PetscFree(count);
1517:     PetscFree(indices);
1518:     ISDestroy(&isnumb);
1519:     ISDestroy(&ispart);

1521:   }
1522:   return(0);
1523: }

1525: /*@C
1526:    PCASMDestroySubdomains - Destroys the index sets created with
1527:    PCASMCreateSubdomains(). Should be called after setting subdomains
1528:    with PCASMSetLocalSubdomains().

1530:    Collective

1532:    Input Parameters:
1533: +  n - the number of index sets
1534: .  is - the array of index sets
1535: -  is_local - the array of local index sets, can be NULL

1537:    Level: advanced

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

1541: .seealso: PCASMCreateSubdomains(), PCASMSetLocalSubdomains()
1542: @*/
1543: PetscErrorCode  PCASMDestroySubdomains(PetscInt n, IS is[], IS is_local[])
1544: {
1545:   PetscInt       i;

1549:   if (n <= 0) return(0);
1550:   if (is) {
1552:     for (i=0; i<n; i++) { ISDestroy(&is[i]); }
1553:     PetscFree(is);
1554:   }
1555:   if (is_local) {
1557:     for (i=0; i<n; i++) { ISDestroy(&is_local[i]); }
1558:     PetscFree(is_local);
1559:   }
1560:   return(0);
1561: }

1563: /*@
1564:    PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1565:    preconditioner for a two-dimensional problem on a regular grid.

1567:    Not Collective

1569:    Input Parameters:
1570: +  m, n - the number of mesh points in the x and y directions
1571: .  M, N - the number of subdomains in the x and y directions
1572: .  dof - degrees of freedom per node
1573: -  overlap - overlap in mesh lines

1575:    Output Parameters:
1576: +  Nsub - the number of subdomains created
1577: .  is - array of index sets defining overlapping (if overlap > 0) subdomains
1578: -  is_local - array of index sets defining non-overlapping subdomains

1580:    Note:
1581:    Presently PCAMSCreateSubdomains2d() is valid only for sequential
1582:    preconditioners.  More general related routines are
1583:    PCASMSetTotalSubdomains() and PCASMSetLocalSubdomains().

1585:    Level: advanced

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

1589: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1590:           PCASMSetOverlap()
1591: @*/
1592: PetscErrorCode  PCASMCreateSubdomains2D(PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt dof,PetscInt overlap,PetscInt *Nsub,IS **is,IS **is_local)
1593: {
1594:   PetscInt       i,j,height,width,ystart,xstart,yleft,yright,xleft,xright,loc_outer;
1596:   PetscInt       nidx,*idx,loc,ii,jj,count;

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

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

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

1647: /*@C
1648:     PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
1649:     only) for the additive Schwarz preconditioner.

1651:     Not Collective

1653:     Input Parameter:
1654: .   pc - the preconditioner context

1656:     Output Parameters:
1657: +   n - the number of subdomains for this processor (default value = 1)
1658: .   is - the index sets that define the subdomains for this processor
1659: -   is_local - the index sets that define the local part of the subdomains for this processor (can be NULL)


1662:     Notes:
1663:     The IS numbering is in the parallel, global numbering of the vector.

1665:     Level: advanced

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

1669: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1670:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubmatrices()
1671: @*/
1672: PetscErrorCode  PCASMGetLocalSubdomains(PC pc,PetscInt *n,IS *is[],IS *is_local[])
1673: {
1674:   PC_ASM         *osm = (PC_ASM*)pc->data;
1676:   PetscBool      match;

1682:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1683:   if (!match) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PC is not a PCASM");
1684:   if (n) *n = osm->n_local_true;
1685:   if (is) *is = osm->is;
1686:   if (is_local) *is_local = osm->is_local;
1687:   return(0);
1688: }

1690: /*@C
1691:     PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1692:     only) for the additive Schwarz preconditioner.

1694:     Not Collective

1696:     Input Parameter:
1697: .   pc - the preconditioner context

1699:     Output Parameters:
1700: +   n - the number of matrices for this processor (default value = 1)
1701: -   mat - the matrices

1703:     Level: advanced

1705:     Notes:
1706:     Call after PCSetUp() (or KSPSetUp()) but before PCApply() (or KSPApply()) and before PCSetUpOnBlocks())

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

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

1712: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1713:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains(), PCSetModifySubmatrices()
1714: @*/
1715: PetscErrorCode  PCASMGetLocalSubmatrices(PC pc,PetscInt *n,Mat *mat[])
1716: {
1717:   PC_ASM         *osm;
1719:   PetscBool      match;

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

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

1741:     Logically Collective

1743:     Input Parameter:
1744: +   pc  - the preconditioner
1745: -   flg - boolean indicating whether to use subdomains defined by the DM

1747:     Options Database Key:
1748: .   -pc_asm_dm_subdomains

1750:     Level: intermediate

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

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

1758: .seealso: PCASMGetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1759:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1760: @*/
1761: PetscErrorCode  PCASMSetDMSubdomains(PC pc,PetscBool flg)
1762: {
1763:   PC_ASM         *osm = (PC_ASM*)pc->data;
1765:   PetscBool      match;

1770:   if (pc->setupcalled) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for a setup PC.");
1771:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1772:   if (match) {
1773:     osm->dm_subdomains = flg;
1774:   }
1775:   return(0);
1776: }

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

1782:     Input Parameter:
1783: .   pc  - the preconditioner

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

1788:     Level: intermediate

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

1792: .seealso: PCASMSetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1793:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1794: @*/
1795: PetscErrorCode  PCASMGetDMSubdomains(PC pc,PetscBool* flg)
1796: {
1797:   PC_ASM         *osm = (PC_ASM*)pc->data;
1799:   PetscBool      match;

1804:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1805:   if (match) {
1806:     if (flg) *flg = osm->dm_subdomains;
1807:   }
1808:   return(0);
1809: }

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

1814:    Not Collective

1816:    Input Parameter:
1817: .  pc - the PC

1819:    Output Parameter:
1820: .  -pc_asm_sub_mat_type - name of matrix type

1822:    Level: advanced

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

1826: .seealso: PCASMSetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1827: @*/
1828: PetscErrorCode  PCASMGetSubMatType(PC pc,MatType *sub_mat_type){

1831:   PetscTryMethod(pc,"PCASMGetSubMatType_C",(PC,MatType*),(pc,sub_mat_type));
1832:   return(0);
1833: }

1835: /*@
1836:      PCASMSetSubMatType - Set the type of matrix used for ASM subsolves

1838:    Collective on Mat

1840:    Input Parameters:
1841: +  pc             - the PC object
1842: -  sub_mat_type   - matrix type

1844:    Options Database Key:
1845: .  -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.

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

1850:   Level: advanced

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

1854: .seealso: PCASMGetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1855: @*/
1856: PetscErrorCode PCASMSetSubMatType(PC pc,MatType sub_mat_type)
1857: {

1860:   PetscTryMethod(pc,"PCASMSetSubMatType_C",(PC,MatType),(pc,sub_mat_type));
1861:   return(0);
1862: }