Actual source code: asm.c

petsc-master 2016-08-28
Report Typos and Errors
  2: /*
  3:   This file defines an additive Schwarz preconditioner for any Mat implementation.

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

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

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

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

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

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

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

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

196:   if (!pc->setupcalled) {

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

203:     /* Note: if subdomains have been set either via PCASMSetTotalSubdomains() or via PCASMSetLocalSubdomains(), osm->n_local_true will not be PETSC_DECIDE */
204:     if (osm->n_local_true == PETSC_DECIDE) {
205:       /* no subdomains given */
206:       /* try pc->dm first, if allowed */
207:       if (osm->dm_subdomains && pc->dm) {
208:         PetscInt  num_domains, d;
209:         char      **domain_names;
210:         IS        *inner_domain_is, *outer_domain_is;
211:         DMCreateDomainDecomposition(pc->dm, &num_domains, &domain_names, &inner_domain_is, &outer_domain_is, &domain_dm);
212:         osm->overlap = -1; /* We do not want to increase the overlap of the IS. 
213:                               A future improvement of this code might allow one to use 
214:                               DM-defined subdomains and also increase the overlap, 
215:                               but that is not currently supported */
216:         if (num_domains) {
217:           PCASMSetLocalSubdomains(pc, num_domains, outer_domain_is, inner_domain_is);
218:         }
219:         for (d = 0; d < num_domains; ++d) {
220:           if (domain_names)    {PetscFree(domain_names[d]);}
221:           if (inner_domain_is) {ISDestroy(&inner_domain_is[d]);}
222:           if (outer_domain_is) {ISDestroy(&outer_domain_is[d]);}
223:         }
224:         PetscFree(domain_names);
225:         PetscFree(inner_domain_is);
226:         PetscFree(outer_domain_is);
227:       }
228:       if (osm->n_local_true == PETSC_DECIDE) {
229:         /* still no subdomains; use one subdomain per processor */
230:         osm->n_local_true = 1;
231:       }
232:     }
233:     { /* determine the global and max number of subdomains */
234:       struct {PetscInt max,sum;} inwork,outwork;
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;
240:     }
241:     if (!osm->is) { /* create the index sets */
242:       PCASMCreateSubdomains(pc->pmat,osm->n_local_true,&osm->is);
243:     }
244:     if (osm->n_local_true > 1 && !osm->is_local) {
245:       PetscMalloc1(osm->n_local_true,&osm->is_local);
246:       for (i=0; i<osm->n_local_true; i++) {
247:         if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
248:           ISDuplicate(osm->is[i],&osm->is_local[i]);
249:           ISCopy(osm->is[i],osm->is_local[i]);
250:         } else {
251:           PetscObjectReference((PetscObject)osm->is[i]);
252:           osm->is_local[i] = osm->is[i];
253:         }
254:       }
255:     }
256:     PCGetOptionsPrefix(pc,&prefix);
257:     flg  = PETSC_FALSE;
258:     PetscOptionsGetBool(NULL,prefix,"-pc_asm_print_subdomains",&flg,NULL);
259:     if (flg) { PCASMPrintSubdomains(pc); }

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

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

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);
309:     }
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:   MatGetSubMatrices(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);
343:     PetscMalloc1(osm->n_local,&osm->restriction);
344:     if (osm->is_local) {PetscMalloc1(osm->n_local,&osm->localization);}
345:     PetscMalloc1(osm->n_local,&osm->prolongation);
346:     PetscMalloc1(osm->n_local,&osm->x);
347:     PetscMalloc1(osm->n_local,&osm->y);
348:     PetscMalloc1(osm->n_local,&osm->y_local);
349:     for (i=0; i<osm->n_local_true; ++i) {
350:       ISGetLocalSize(osm->is[i],&m);
351:       MatCreateVecs(osm->pmat[i],&osm->x[i],NULL);
352:       ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
353:       VecScatterCreate(vec,osm->is[i],osm->x[i],isl,&osm->restriction[i]);
354:       ISDestroy(&isl);
355:       VecDuplicate(osm->x[i],&osm->y[i]);
356:       if (osm->is_local) {
357:         ISLocalToGlobalMapping ltog;
358:         IS                     isll;
359:         const PetscInt         *idx_local;
360:         PetscInt               *idx,nout;

362:         ISLocalToGlobalMappingCreateIS(osm->is[i],&ltog);
363:         ISGetLocalSize(osm->is_local[i],&m_local);
364:         ISGetIndices(osm->is_local[i], &idx_local);
365:         PetscMalloc1(m_local,&idx);
366:         ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx);
367:         ISLocalToGlobalMappingDestroy(&ltog);
368:         if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of is");
369:         ISRestoreIndices(osm->is_local[i], &idx_local);
370:         ISCreateGeneral(PETSC_COMM_SELF,m_local,idx,PETSC_OWN_POINTER,&isll);
371:         ISCreateStride(PETSC_COMM_SELF,m_local,0,1,&isl);
372:         VecCreateSeq(PETSC_COMM_SELF,m_local,&osm->y_local[i]);
373:         VecScatterCreate(osm->y[i],isll,osm->y_local[i],isl,&osm->localization[i]);
374:         ISDestroy(&isll);

376:         VecScatterCreate(vec,osm->is_local[i],osm->y_local[i],isl,&osm->prolongation[i]);
377:         ISDestroy(&isl);
378:       } else {
379:         VecGetLocalSize(vec,&m_local);

381:         osm->y_local[i] = osm->y[i];

383:         PetscObjectReference((PetscObject) osm->y[i]);

385:         osm->prolongation[i] = osm->restriction[i];

387:         PetscObjectReference((PetscObject) osm->restriction[i]);
388:       }
389:     }
390:     for (i=osm->n_local_true; i<osm->n_local; i++) {
391:       VecCreateSeq(PETSC_COMM_SELF,0,&osm->x[i]);
392:       VecDuplicate(osm->x[i],&osm->y[i]);
393:       VecDuplicate(osm->x[i],&osm->y_local[i]);
394:       ISCreateStride(PETSC_COMM_SELF,0,0,1,&isl);
395:       VecScatterCreate(vec,isl,osm->x[i],isl,&osm->restriction[i]);
396:       if (osm->is_local) {
397:         VecScatterCreate(osm->y[i],isl,osm->y_local[i],isl,&osm->localization[i]);
398:         VecScatterCreate(vec,isl,osm->x[i],isl,&osm->prolongation[i]);
399:       } else {
400:         osm->prolongation[i] = osm->restriction[i];
401:         PetscObjectReference((PetscObject) osm->restriction[i]);
402:       }
403:       ISDestroy(&isl);
404:     }
405:     VecDestroy(&vec);
406:   }

408:   if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
409:     IS      *cis;
410:     PetscInt c;

412:     PetscMalloc1(osm->n_local_true, &cis);
413:     for (c = 0; c < osm->n_local_true; ++c) cis[c] = osm->lis;
414:     MatGetSubMatrices(pc->pmat, osm->n_local_true, osm->is, cis, scall, &osm->lmats);
415:     PetscFree(cis);
416:   }

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

422:   /*
423:      Loop over subdomains putting them into local ksp
424:   */
425:   for (i=0; i<osm->n_local_true; i++) {
426:     KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i]);
427:     if (!pc->setupcalled) {
428:       KSPSetFromOptions(osm->ksp[i]);
429:     }
430:   }
431:   return(0);
432: }

436: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
437: {
438:   PC_ASM             *osm = (PC_ASM*)pc->data;
439:   PetscErrorCode     ierr;
440:   PetscInt           i;
441:   KSPConvergedReason reason;

444:   for (i=0; i<osm->n_local_true; i++) {
445:     KSPSetUp(osm->ksp[i]);
446:     KSPGetConvergedReason(osm->ksp[i],&reason);
447:     if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
448:       pc->failedreason = PC_SUBPC_ERROR;
449:     }
450:   }
451:   return(0);
452: }

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

464:   /*
465:      Support for limiting the restriction or interpolation to only local
466:      subdomain values (leaving the other values 0).
467:   */
468:   if (!(osm->type & PC_ASM_RESTRICT)) {
469:     forward = SCATTER_FORWARD_LOCAL;
470:     /* have to zero the work RHS since scatter may leave some slots empty */
471:     for (i=0; i<n_local_true; i++) {
472:       VecZeroEntries(osm->x[i]);
473:     }
474:   }
475:   if (!(osm->type & PC_ASM_INTERPOLATE)) reverse = SCATTER_REVERSE_LOCAL;

477:   switch (osm->loctype)
478:   {
479:   case PC_COMPOSITE_ADDITIVE:
480:     for (i=0; i<n_local; i++) {
481:       VecScatterBegin(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
482:     }
483:     VecZeroEntries(y);
484:     /* do the local solves */
485:     for (i=0; i<n_local_true; i++) {
486:       VecScatterEnd(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
487:       KSPSolve(osm->ksp[i],osm->x[i],osm->y[i]);
488:       if (osm->localization) {
489:         VecScatterBegin(osm->localization[i],osm->y[i],osm->y_local[i],INSERT_VALUES,forward);
490:         VecScatterEnd(osm->localization[i],osm->y[i],osm->y_local[i],INSERT_VALUES,forward);
491:       }
492:       VecScatterBegin(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
493:     }
494:     /* handle the rest of the scatters that do not have local solves */
495:     for (i=n_local_true; i<n_local; i++) {
496:       VecScatterEnd(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
497:       VecScatterBegin(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
498:     }
499:     for (i=0; i<n_local; i++) {
500:       VecScatterEnd(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
501:     }
502:     break;
503:   case PC_COMPOSITE_MULTIPLICATIVE:
504:     VecZeroEntries(y);
505:     /* do the local solves */
506:     for (i = 0; i < n_local_true; ++i) {
507:       if (i > 0) {
508:         /* Update rhs */
509:         VecScatterBegin(osm->restriction[i], osm->lx, osm->x[i], INSERT_VALUES, forward);
510:         VecScatterEnd(osm->restriction[i], osm->lx, osm->x[i], INSERT_VALUES, forward);
511:       } else {
512:         VecZeroEntries(osm->x[i]);
513:       }
514:       VecScatterBegin(osm->restriction[i], x, osm->x[i], ADD_VALUES, forward);
515:       VecScatterEnd(osm->restriction[i], x, osm->x[i], ADD_VALUES, forward);
516:       KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]);
517:       if (osm->localization) {
518:         VecScatterBegin(osm->localization[i], osm->y[i], osm->y_local[i], INSERT_VALUES, forward);
519:         VecScatterEnd(osm->localization[i], osm->y[i], osm->y_local[i], INSERT_VALUES, forward);
520:       }
521:       VecScatterBegin(osm->prolongation[i], osm->y_local[i], y, ADD_VALUES, reverse);
522:       VecScatterEnd(osm->prolongation[i], osm->y_local[i], y, ADD_VALUES, reverse);
523:       if (i < n_local_true-1) {
524:         VecSet(osm->ly, 0.0);
525:         VecScatterBegin(osm->prolongation[i], osm->y_local[i], osm->ly, INSERT_VALUES, reverse);
526:         VecScatterEnd(osm->prolongation[i], osm->y_local[i], osm->ly, INSERT_VALUES, reverse);
527:         VecScale(osm->ly, -1.0);
528:         MatMult(osm->lmats[i+1], osm->ly, osm->y[i+1]);
529:         VecScatterBegin(osm->restriction[i+1], osm->y[i+1], osm->lx, INSERT_VALUES, reverse);
530:         VecScatterEnd(osm->restriction[i+1], osm->y[i+1], osm->lx, INSERT_VALUES, reverse);
531:       }
532:     }
533:     /* handle the rest of the scatters that do not have local solves */
534:     for (i = n_local_true; i < n_local; ++i) {
535:       VecScatterBegin(osm->restriction[i], x, osm->x[i], INSERT_VALUES, forward);
536:       VecScatterEnd(osm->restriction[i], x, osm->x[i], INSERT_VALUES, forward);
537:       VecScatterBegin(osm->prolongation[i], osm->y_local[i], y, ADD_VALUES, reverse);
538:       VecScatterEnd(osm->prolongation[i], osm->y_local[i], y, ADD_VALUES, reverse);
539:     }
540:     break;
541:   default: SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
542:   }
543:   return(0);
544: }

548: static PetscErrorCode PCApplyTranspose_ASM(PC pc,Vec x,Vec y)
549: {
550:   PC_ASM         *osm = (PC_ASM*)pc->data;
552:   PetscInt       i,n_local = osm->n_local,n_local_true = osm->n_local_true;
553:   ScatterMode    forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;

556:   /*
557:      Support for limiting the restriction or interpolation to only local
558:      subdomain values (leaving the other values 0).

560:      Note: these are reversed from the PCApply_ASM() because we are applying the
561:      transpose of the three terms
562:   */
563:   if (!(osm->type & PC_ASM_INTERPOLATE)) {
564:     forward = SCATTER_FORWARD_LOCAL;
565:     /* have to zero the work RHS since scatter may leave some slots empty */
566:     for (i=0; i<n_local_true; i++) {
567:       VecZeroEntries(osm->x[i]);
568:     }
569:   }
570:   if (!(osm->type & PC_ASM_RESTRICT)) reverse = SCATTER_REVERSE_LOCAL;

572:   for (i=0; i<n_local; i++) {
573:     VecScatterBegin(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
574:   }
575:   VecZeroEntries(y);
576:   /* do the local solves */
577:   for (i=0; i<n_local_true; i++) {
578:     VecScatterEnd(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
579:     KSPSolveTranspose(osm->ksp[i],osm->x[i],osm->y[i]);
580:     if (osm->localization) {
581:       VecScatterBegin(osm->localization[i],osm->y[i],osm->y_local[i],INSERT_VALUES,forward);
582:       VecScatterEnd(osm->localization[i],osm->y[i],osm->y_local[i],INSERT_VALUES,forward);
583:     }
584:     VecScatterBegin(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
585:   }
586:   /* handle the rest of the scatters that do not have local solves */
587:   for (i=n_local_true; i<n_local; i++) {
588:     VecScatterEnd(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
589:     VecScatterBegin(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
590:   }
591:   for (i=0; i<n_local; i++) {
592:     VecScatterEnd(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
593:   }
594:   return(0);
595: }

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

606:   if (osm->ksp) {
607:     for (i=0; i<osm->n_local_true; i++) {
608:       KSPReset(osm->ksp[i]);
609:     }
610:   }
611:   if (osm->pmat) {
612:     if (osm->n_local_true > 0) {
613:       MatDestroyMatrices(osm->n_local_true,&osm->pmat);
614:     }
615:   }
616:   if (osm->restriction) {
617:     for (i=0; i<osm->n_local; i++) {
618:       VecScatterDestroy(&osm->restriction[i]);
619:       if (osm->localization) {VecScatterDestroy(&osm->localization[i]);}
620:       VecScatterDestroy(&osm->prolongation[i]);
621:       VecDestroy(&osm->x[i]);
622:       VecDestroy(&osm->y[i]);
623:       VecDestroy(&osm->y_local[i]);
624:     }
625:     PetscFree(osm->restriction);
626:     if (osm->localization) {PetscFree(osm->localization);}
627:     PetscFree(osm->prolongation);
628:     PetscFree(osm->x);
629:     PetscFree(osm->y);
630:     PetscFree(osm->y_local);
631:   }
632:   PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
633:   if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
634:     ISDestroy(&osm->lis);
635:     MatDestroyMatrices(osm->n_local_true, &osm->lmats);
636:     VecDestroy(&osm->lx);
637:     VecDestroy(&osm->ly);
638:   }

640:   PetscFree(osm->sub_mat_type);

642:   osm->is       = 0;
643:   osm->is_local = 0;
644:   return(0);
645: }

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: }

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

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

711: /*------------------------------------------------------------------------------------*/

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

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

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

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

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

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

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

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

796: static PetscErrorCode  PCASMSetOverlap_ASM(PC pc,PetscInt ovl)
797: {
798:   PC_ASM *osm = (PC_ASM*)pc->data;

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

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

814:   osm->type     = type;
815:   osm->type_set = PETSC_TRUE;
816:   return(0);
817: }

821: static PetscErrorCode  PCASMGetType_ASM(PC pc,PCASMType *type)
822: {
823:   PC_ASM *osm = (PC_ASM*)pc->data;

826:   *type = osm->type;
827:   return(0);
828: }

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

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

843: static PetscErrorCode  PCASMGetLocalType_ASM(PC pc, PCCompositeType *type)
844: {
845:   PC_ASM *osm = (PC_ASM *) pc->data;

848:   *type = osm->loctype;
849:   return(0);
850: }

854: static PetscErrorCode  PCASMSetSortIndices_ASM(PC pc,PetscBool  doSort)
855: {
856:   PC_ASM *osm = (PC_ASM*)pc->data;

859:   osm->sort_indices = doSort;
860:   return(0);
861: }

865: static PetscErrorCode  PCASMGetSubKSP_ASM(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
866: {
867:   PC_ASM         *osm = (PC_ASM*)pc->data;

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

873:   if (n_local) *n_local = osm->n_local_true;
874:   if (first_local) {
875:     MPI_Scan(&osm->n_local_true,first_local,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
876:     *first_local -= osm->n_local_true;
877:   }
878:   if (ksp) {
879:     /* Assume that local solves are now different; not necessarily
880:        true though!  This flag is used only for PCView_ASM() */
881:     *ksp                   = osm->ksp;
882:     osm->same_local_solves = PETSC_FALSE;
883:   }
884:   return(0);
885: }

889: static PetscErrorCode  PCASMGetSubMatType_ASM(PC pc,MatType *sub_mat_type)
890: {
891:   PC_ASM         *osm = (PC_ASM*)pc->data;

896:   *sub_mat_type = osm->sub_mat_type;
897:   return(0);
898: }

902: static PetscErrorCode PCASMSetSubMatType_ASM(PC pc,MatType sub_mat_type)
903: {
904:   PetscErrorCode    ierr;
905:   PC_ASM            *osm = (PC_ASM*)pc->data;

909:   PetscFree(osm->sub_mat_type);
910:   PetscStrallocpy(sub_mat_type,(char**)&osm->sub_mat_type);
911:   return(0);
912: }

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

919:     Collective on PC

921:     Input Parameters:
922: +   pc - the preconditioner context
923: .   n - the number of subdomains for this processor (default value = 1)
924: .   is - the index set that defines the subdomains for this processor
925:          (or NULL for PETSc to determine subdomains)
926: -   is_local - the index sets that define the local part of the subdomains for this processor
927:          (or NULL to use the default of 1 subdomain per process)

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

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

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

936:     Level: advanced

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

940: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
941:           PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains()
942: @*/
943: PetscErrorCode  PCASMSetLocalSubdomains(PC pc,PetscInt n,IS is[],IS is_local[])
944: {

949:   PetscTryMethod(pc,"PCASMSetLocalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,is,is_local));
950:   return(0);
951: }

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

960:     Collective on PC

962:     Input Parameters:
963: +   pc - the preconditioner context
964: .   N  - the number of subdomains for all processors
965: .   is - the index sets that define the subdomains for all processors
966:          (or NULL to ask PETSc to compe up with subdomains)
967: -   is_local - the index sets that define the local part of the subdomains for this processor
968:          (or NULL to use the default of 1 subdomain per process)

970:     Options Database Key:
971:     To set the total number of subdomain blocks rather than specify the
972:     index sets, use the option
973: .    -pc_asm_blocks <blks> - Sets total blocks

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

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

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

983:     Use PCASMSetLocalSubdomains() to set local subdomains.

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

987:     Level: advanced

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

991: .seealso: PCASMSetLocalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
992:           PCASMCreateSubdomains2D()
993: @*/
994: PetscErrorCode  PCASMSetTotalSubdomains(PC pc,PetscInt N,IS is[],IS is_local[])
995: {

1000:   PetscTryMethod(pc,"PCASMSetTotalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,N,is,is_local));
1001:   return(0);
1002: }

1006: /*@
1007:     PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
1008:     additive Schwarz preconditioner.  Either all or no processors in the
1009:     PC communicator must call this routine. If MatIncreaseOverlap is used,
1010:     use option -mat_increase_overlap when the problem size large.

1012:     Logically Collective on PC

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

1018:     Options Database Key:
1019: .   -pc_asm_overlap <ovl> - Sets overlap

1021:     Notes:
1022:     By default the ASM preconditioner uses 1 block per processor.  To use
1023:     multiple blocks per perocessor, see PCASMSetTotalSubdomains() and
1024:     PCASMSetLocalSubdomains() (and the option -pc_asm_blocks <blks>).

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

1034:     Note that one can define initial index sets with any overlap via
1035:     PCASMSetTotalSubdomains() or PCASMSetLocalSubdomains(); the routine
1036:     PCASMSetOverlap() merely allows PETSc to extend that overlap further
1037:     if desired.

1039:     Level: intermediate

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

1043: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1044:           PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains()
1045: @*/
1046: PetscErrorCode  PCASMSetOverlap(PC pc,PetscInt ovl)
1047: {

1053:   PetscTryMethod(pc,"PCASMSetOverlap_C",(PC,PetscInt),(pc,ovl));
1054:   return(0);
1055: }

1059: /*@
1060:     PCASMSetType - Sets the type of restriction and interpolation used
1061:     for local problems in the additive Schwarz method.

1063:     Logically Collective on PC

1065:     Input Parameters:
1066: +   pc  - the preconditioner context
1067: -   type - variant of ASM, one of
1068: .vb
1069:       PC_ASM_BASIC       - full interpolation and restriction
1070:       PC_ASM_RESTRICT    - full restriction, local processor interpolation
1071:       PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1072:       PC_ASM_NONE        - local processor restriction and interpolation
1073: .ve

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

1078:     Level: intermediate

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

1082: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1083:           PCASMCreateSubdomains2D()
1084: @*/
1085: PetscErrorCode  PCASMSetType(PC pc,PCASMType type)
1086: {

1092:   PetscTryMethod(pc,"PCASMSetType_C",(PC,PCASMType),(pc,type));
1093:   return(0);
1094: }

1098: /*@
1099:     PCASMGetType - Gets the type of restriction and interpolation used
1100:     for local problems in the additive Schwarz method.

1102:     Logically Collective on PC

1104:     Input Parameter:
1105: .   pc  - the preconditioner context

1107:     Output Parameter:
1108: .   type - variant of ASM, one of

1110: .vb
1111:       PC_ASM_BASIC       - full interpolation and restriction
1112:       PC_ASM_RESTRICT    - full restriction, local processor interpolation
1113:       PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1114:       PC_ASM_NONE        - local processor restriction and interpolation
1115: .ve

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

1120:     Level: intermediate

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

1124: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1125:           PCASMCreateSubdomains2D()
1126: @*/
1127: PetscErrorCode  PCASMGetType(PC pc,PCASMType *type)
1128: {

1133:   PetscUseMethod(pc,"PCASMGetType_C",(PC,PCASMType*),(pc,type));
1134:   return(0);
1135: }

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

1142:   Logically Collective on PC

1144:   Input Parameters:
1145: + pc  - the preconditioner context
1146: - type - type of composition, one of
1147: .vb
1148:   PC_COMPOSITE_ADDITIVE       - local additive combination
1149:   PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1150: .ve

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

1155:   Level: intermediate

1157: .seealso: PCASMSetType(), PCASMGetType(), PCASMGetLocalType(), PCASMCreate()
1158: @*/
1159: PetscErrorCode PCASMSetLocalType(PC pc, PCCompositeType type)
1160: {

1166:   PetscTryMethod(pc, "PCASMSetLocalType_C", (PC, PCCompositeType), (pc, type));
1167:   return(0);
1168: }

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

1175:   Logically Collective on PC

1177:   Input Parameter:
1178: . pc  - the preconditioner context

1180:   Output Parameter:
1181: . type - type of composition, one of
1182: .vb
1183:   PC_COMPOSITE_ADDITIVE       - local additive combination
1184:   PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1185: .ve

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

1190:   Level: intermediate

1192: .seealso: PCASMSetType(), PCASMGetType(), PCASMSetLocalType(), PCASMCreate()
1193: @*/
1194: PetscErrorCode PCASMGetLocalType(PC pc, PCCompositeType *type)
1195: {

1201:   PetscUseMethod(pc, "PCASMGetLocalType_C", (PC, PCCompositeType *), (pc, type));
1202:   return(0);
1203: }

1207: /*@
1208:     PCASMSetSortIndices - Determines whether subdomain indices are sorted.

1210:     Logically Collective on PC

1212:     Input Parameters:
1213: +   pc  - the preconditioner context
1214: -   doSort - sort the subdomain indices

1216:     Level: intermediate

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

1220: .seealso: PCASMSetLocalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1221:           PCASMCreateSubdomains2D()
1222: @*/
1223: PetscErrorCode  PCASMSetSortIndices(PC pc,PetscBool doSort)
1224: {

1230:   PetscTryMethod(pc,"PCASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));
1231:   return(0);
1232: }

1236: /*@C
1237:    PCASMGetSubKSP - Gets the local KSP contexts for all blocks on
1238:    this processor.

1240:    Collective on PC iff first_local is requested

1242:    Input Parameter:
1243: .  pc - the preconditioner context

1245:    Output Parameters:
1246: +  n_local - the number of blocks on this processor or NULL
1247: .  first_local - the global number of the first block on this processor or NULL,
1248:                  all processors must request or all must pass NULL
1249: -  ksp - the array of KSP contexts

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

1254:    Currently for some matrix implementations only 1 block per processor
1255:    is supported.

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

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

1262:    Level: advanced

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

1266: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap(),
1267:           PCASMCreateSubdomains2D(),
1268: @*/
1269: PetscErrorCode  PCASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
1270: {

1275:   PetscUseMethod(pc,"PCASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
1276:   return(0);
1277: }

1279: /* -------------------------------------------------------------------------------------*/
1280: /*MC
1281:    PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1282:            its own KSP object.

1284:    Options Database Keys:
1285: +  -pc_asm_blocks <blks> - Sets total blocks
1286: .  -pc_asm_overlap <ovl> - Sets overlap
1287: -  -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type

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

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

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

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


1304:    Level: beginner

1306:    Concepts: additive Schwarz method

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

1314: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
1315:            PCBJACOBI, PCASMGetSubKSP(), PCASMSetLocalSubdomains(),
1316:            PCASMSetTotalSubdomains(), PCSetModifySubmatrices(), PCASMSetOverlap(), PCASMSetType()

1318: M*/

1322: PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1323: {
1325:   PC_ASM         *osm;

1328:   PetscNewLog(pc,&osm);

1330:   osm->n                 = PETSC_DECIDE;
1331:   osm->n_local           = 0;
1332:   osm->n_local_true      = PETSC_DECIDE;
1333:   osm->overlap           = 1;
1334:   osm->ksp               = 0;
1335:   osm->restriction       = 0;
1336:   osm->localization      = 0;
1337:   osm->prolongation      = 0;
1338:   osm->x                 = 0;
1339:   osm->y                 = 0;
1340:   osm->y_local           = 0;
1341:   osm->is                = 0;
1342:   osm->is_local          = 0;
1343:   osm->mat               = 0;
1344:   osm->pmat              = 0;
1345:   osm->type              = PC_ASM_RESTRICT;
1346:   osm->loctype           = PC_COMPOSITE_ADDITIVE;
1347:   osm->same_local_solves = PETSC_TRUE;
1348:   osm->sort_indices      = PETSC_TRUE;
1349:   osm->dm_subdomains     = PETSC_FALSE;
1350:   osm->sub_mat_type      = NULL;

1352:   pc->data                 = (void*)osm;
1353:   pc->ops->apply           = PCApply_ASM;
1354:   pc->ops->applytranspose  = PCApplyTranspose_ASM;
1355:   pc->ops->setup           = PCSetUp_ASM;
1356:   pc->ops->reset           = PCReset_ASM;
1357:   pc->ops->destroy         = PCDestroy_ASM;
1358:   pc->ops->setfromoptions  = PCSetFromOptions_ASM;
1359:   pc->ops->setuponblocks   = PCSetUpOnBlocks_ASM;
1360:   pc->ops->view            = PCView_ASM;
1361:   pc->ops->applyrichardson = 0;

1363:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",PCASMSetLocalSubdomains_ASM);
1364:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",PCASMSetTotalSubdomains_ASM);
1365:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetOverlap_C",PCASMSetOverlap_ASM);
1366:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetType_C",PCASMSetType_ASM);
1367:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetType_C",PCASMGetType_ASM);
1368:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalType_C",PCASMSetLocalType_ASM);
1369:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetLocalType_C",PCASMGetLocalType_ASM);
1370:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSortIndices_C",PCASMSetSortIndices_ASM);
1371:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubKSP_C",PCASMGetSubKSP_ASM);
1372:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubMatType_C",PCASMGetSubMatType_ASM);
1373:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSubMatType_C",PCASMSetSubMatType_ASM);
1374:   return(0);
1375: }

1379: /*@C
1380:    PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1381:    preconditioner for a any problem on a general grid.

1383:    Collective

1385:    Input Parameters:
1386: +  A - The global matrix operator
1387: -  n - the number of local blocks

1389:    Output Parameters:
1390: .  outis - the array of index sets defining the subdomains

1392:    Level: advanced

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

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

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

1401: .seealso: PCASMSetLocalSubdomains(), PCASMDestroySubdomains()
1402: @*/
1403: PetscErrorCode  PCASMCreateSubdomains(Mat A, PetscInt n, IS* outis[])
1404: {
1405:   MatPartitioning mpart;
1406:   const char      *prefix;
1407:   PetscErrorCode  (*f)(Mat,Mat*);
1408:   PetscMPIInt     size;
1409:   PetscInt        i,j,rstart,rend,bs;
1410:   PetscBool       isbaij = PETSC_FALSE,foundpart = PETSC_FALSE;
1411:   Mat             Ad     = NULL, adj;
1412:   IS              ispart,isnumb,*is;
1413:   PetscErrorCode  ierr;

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

1420:   /* Get prefix, row distribution, and block size */
1421:   MatGetOptionsPrefix(A,&prefix);
1422:   MatGetOwnershipRange(A,&rstart,&rend);
1423:   MatGetBlockSize(A,&bs);
1424:   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);

1426:   /* Get diagonal block from matrix if possible */
1427:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1428:   PetscObjectQueryFunction((PetscObject)A,"MatGetDiagonalBlock_C",&f);
1429:   if (f) {
1430:     MatGetDiagonalBlock(A,&Ad);
1431:   } else if (size == 1) {
1432:     Ad = A;
1433:   }
1434:   if (Ad) {
1435:     PetscObjectTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij);
1436:     if (!isbaij) {PetscObjectTypeCompare((PetscObject)Ad,MATSEQSBAIJ,&isbaij);}
1437:   }
1438:   if (Ad && n > 1) {
1439:     PetscBool match,done;
1440:     /* Try to setup a good matrix partitioning if available */
1441:     MatPartitioningCreate(PETSC_COMM_SELF,&mpart);
1442:     PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);
1443:     MatPartitioningSetFromOptions(mpart);
1444:     PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGCURRENT,&match);
1445:     if (!match) {
1446:       PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGSQUARE,&match);
1447:     }
1448:     if (!match) { /* assume a "good" partitioner is available */
1449:       PetscInt       na;
1450:       const PetscInt *ia,*ja;
1451:       MatGetRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1452:       if (done) {
1453:         /* Build adjacency matrix by hand. Unfortunately a call to
1454:            MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1455:            remove the block-aij structure and we cannot expect
1456:            MatPartitioning to split vertices as we need */
1457:         PetscInt       i,j,len,nnz,cnt,*iia=0,*jja=0;
1458:         const PetscInt *row;
1459:         nnz = 0;
1460:         for (i=0; i<na; i++) { /* count number of nonzeros */
1461:           len = ia[i+1] - ia[i];
1462:           row = ja + ia[i];
1463:           for (j=0; j<len; j++) {
1464:             if (row[j] == i) { /* don't count diagonal */
1465:               len--; break;
1466:             }
1467:           }
1468:           nnz += len;
1469:         }
1470:         PetscMalloc1(na+1,&iia);
1471:         PetscMalloc1(nnz,&jja);
1472:         nnz    = 0;
1473:         iia[0] = 0;
1474:         for (i=0; i<na; i++) { /* fill adjacency */
1475:           cnt = 0;
1476:           len = ia[i+1] - ia[i];
1477:           row = ja + ia[i];
1478:           for (j=0; j<len; j++) {
1479:             if (row[j] != i) { /* if not diagonal */
1480:               jja[nnz+cnt++] = row[j];
1481:             }
1482:           }
1483:           nnz     += cnt;
1484:           iia[i+1] = nnz;
1485:         }
1486:         /* Partitioning of the adjacency matrix */
1487:         MatCreateMPIAdj(PETSC_COMM_SELF,na,na,iia,jja,NULL,&adj);
1488:         MatPartitioningSetAdjacency(mpart,adj);
1489:         MatPartitioningSetNParts(mpart,n);
1490:         MatPartitioningApply(mpart,&ispart);
1491:         ISPartitioningToNumbering(ispart,&isnumb);
1492:         MatDestroy(&adj);
1493:         foundpart = PETSC_TRUE;
1494:       }
1495:       MatRestoreRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1496:     }
1497:     MatPartitioningDestroy(&mpart);
1498:   }

1500:   PetscMalloc1(n,&is);
1501:   *outis = is;

1503:   if (!foundpart) {

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

1507:     PetscInt mbs   = (rend-rstart)/bs;
1508:     PetscInt start = rstart;
1509:     for (i=0; i<n; i++) {
1510:       PetscInt count = (mbs/n + ((mbs % n) > i)) * bs;
1511:       ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);
1512:       start += count;
1513:     }

1515:   } else {

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

1519:     const PetscInt *numbering;
1520:     PetscInt       *count,nidx,*indices,*newidx,start=0;
1521:     /* Get node count in each partition */
1522:     PetscMalloc1(n,&count);
1523:     ISPartitioningCount(ispart,n,count);
1524:     if (isbaij && bs > 1) { /* adjust for the block-aij case */
1525:       for (i=0; i<n; i++) count[i] *= bs;
1526:     }
1527:     /* Build indices from node numbering */
1528:     ISGetLocalSize(isnumb,&nidx);
1529:     PetscMalloc1(nidx,&indices);
1530:     for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */
1531:     ISGetIndices(isnumb,&numbering);
1532:     PetscSortIntWithPermutation(nidx,numbering,indices);
1533:     ISRestoreIndices(isnumb,&numbering);
1534:     if (isbaij && bs > 1) { /* adjust for the block-aij case */
1535:       PetscMalloc1(nidx*bs,&newidx);
1536:       for (i=0; i<nidx; i++) {
1537:         for (j=0; j<bs; j++) newidx[i*bs+j] = indices[i]*bs + j;
1538:       }
1539:       PetscFree(indices);
1540:       nidx   *= bs;
1541:       indices = newidx;
1542:     }
1543:     /* Shift to get global indices */
1544:     for (i=0; i<nidx; i++) indices[i] += rstart;

1546:     /* Build the index sets for each block */
1547:     for (i=0; i<n; i++) {
1548:       ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);
1549:       ISSort(is[i]);
1550:       start += count[i];
1551:     }

1553:     PetscFree(count);
1554:     PetscFree(indices);
1555:     ISDestroy(&isnumb);
1556:     ISDestroy(&ispart);

1558:   }
1559:   return(0);
1560: }

1564: /*@C
1565:    PCASMDestroySubdomains - Destroys the index sets created with
1566:    PCASMCreateSubdomains(). Should be called after setting subdomains
1567:    with PCASMSetLocalSubdomains().

1569:    Collective

1571:    Input Parameters:
1572: +  n - the number of index sets
1573: .  is - the array of index sets
1574: -  is_local - the array of local index sets, can be NULL

1576:    Level: advanced

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

1580: .seealso: PCASMCreateSubdomains(), PCASMSetLocalSubdomains()
1581: @*/
1582: PetscErrorCode  PCASMDestroySubdomains(PetscInt n, IS is[], IS is_local[])
1583: {
1584:   PetscInt       i;

1588:   if (n <= 0) return(0);
1589:   if (is) {
1591:     for (i=0; i<n; i++) { ISDestroy(&is[i]); }
1592:     PetscFree(is);
1593:   }
1594:   if (is_local) {
1596:     for (i=0; i<n; i++) { ISDestroy(&is_local[i]); }
1597:     PetscFree(is_local);
1598:   }
1599:   return(0);
1600: }

1604: /*@
1605:    PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1606:    preconditioner for a two-dimensional problem on a regular grid.

1608:    Not Collective

1610:    Input Parameters:
1611: +  m, n - the number of mesh points in the x and y directions
1612: .  M, N - the number of subdomains in the x and y directions
1613: .  dof - degrees of freedom per node
1614: -  overlap - overlap in mesh lines

1616:    Output Parameters:
1617: +  Nsub - the number of subdomains created
1618: .  is - array of index sets defining overlapping (if overlap > 0) subdomains
1619: -  is_local - array of index sets defining non-overlapping subdomains

1621:    Note:
1622:    Presently PCAMSCreateSubdomains2d() is valid only for sequential
1623:    preconditioners.  More general related routines are
1624:    PCASMSetTotalSubdomains() and PCASMSetLocalSubdomains().

1626:    Level: advanced

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

1630: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1631:           PCASMSetOverlap()
1632: @*/
1633: PetscErrorCode  PCASMCreateSubdomains2D(PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt dof,PetscInt overlap,PetscInt *Nsub,IS **is,IS **is_local)
1634: {
1635:   PetscInt       i,j,height,width,ystart,xstart,yleft,yright,xleft,xright,loc_outer;
1637:   PetscInt       nidx,*idx,loc,ii,jj,count;

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

1642:   *Nsub     = N*M;
1643:   PetscMalloc1(*Nsub,is);
1644:   PetscMalloc1(*Nsub,is_local);
1645:   ystart    = 0;
1646:   loc_outer = 0;
1647:   for (i=0; i<N; i++) {
1648:     height = n/N + ((n % N) > i); /* height of subdomain */
1649:     if (height < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many N subdomains for mesh dimension n");
1650:     yleft  = ystart - overlap; if (yleft < 0) yleft = 0;
1651:     yright = ystart + height + overlap; if (yright > n) yright = n;
1652:     xstart = 0;
1653:     for (j=0; j<M; j++) {
1654:       width = m/M + ((m % M) > j); /* width of subdomain */
1655:       if (width < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many M subdomains for mesh dimension m");
1656:       xleft  = xstart - overlap; if (xleft < 0) xleft = 0;
1657:       xright = xstart + width + overlap; if (xright > m) xright = m;
1658:       nidx   = (xright - xleft)*(yright - yleft);
1659:       PetscMalloc1(nidx,&idx);
1660:       loc    = 0;
1661:       for (ii=yleft; ii<yright; ii++) {
1662:         count = m*ii + xleft;
1663:         for (jj=xleft; jj<xright; jj++) idx[loc++] = count++;
1664:       }
1665:       ISCreateGeneral(PETSC_COMM_SELF,nidx,idx,PETSC_COPY_VALUES,(*is)+loc_outer);
1666:       if (overlap == 0) {
1667:         PetscObjectReference((PetscObject)(*is)[loc_outer]);

1669:         (*is_local)[loc_outer] = (*is)[loc_outer];
1670:       } else {
1671:         for (loc=0,ii=ystart; ii<ystart+height; ii++) {
1672:           for (jj=xstart; jj<xstart+width; jj++) {
1673:             idx[loc++] = m*ii + jj;
1674:           }
1675:         }
1676:         ISCreateGeneral(PETSC_COMM_SELF,loc,idx,PETSC_COPY_VALUES,*is_local+loc_outer);
1677:       }
1678:       PetscFree(idx);
1679:       xstart += width;
1680:       loc_outer++;
1681:     }
1682:     ystart += height;
1683:   }
1684:   for (i=0; i<*Nsub; i++) { ISSort((*is)[i]); }
1685:   return(0);
1686: }

1690: /*@C
1691:     PCASMGetLocalSubdomains - Gets the local subdomains (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 subdomains for this processor (default value = 1)
1701: .   is - the index sets that define the subdomains for this processor
1702: -   is_local - the index sets that define the local part of the subdomains for this processor (can be NULL)


1705:     Notes:
1706:     The IS numbering is in the parallel, global numbering of the vector.

1708:     Level: advanced

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

1712: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1713:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubmatrices()
1714: @*/
1715: PetscErrorCode  PCASMGetLocalSubdomains(PC pc,PetscInt *n,IS *is[],IS *is_local[])
1716: {
1717:   PC_ASM         *osm = (PC_ASM*)pc->data;
1719:   PetscBool      match;

1725:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1726:   if (!match) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PC is not a PCASM");
1727:   if (n) *n = osm->n_local_true;
1728:   if (is) *is = osm->is;
1729:   if (is_local) *is_local = osm->is_local;
1730:   return(0);
1731: }

1735: /*@C
1736:     PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1737:     only) for the additive Schwarz preconditioner.

1739:     Not Collective

1741:     Input Parameter:
1742: .   pc - the preconditioner context

1744:     Output Parameters:
1745: +   n - the number of matrices for this processor (default value = 1)
1746: -   mat - the matrices


1749:     Level: advanced

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

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

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

1757: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1758:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains(), PCSetModifySubmatrices()
1759: @*/
1760: PetscErrorCode  PCASMGetLocalSubmatrices(PC pc,PetscInt *n,Mat *mat[])
1761: {
1762:   PC_ASM         *osm;
1764:   PetscBool      match;

1770:   if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUP() or PCSetUp().");
1771:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1772:   if (!match) {
1773:     if (n) *n = 0;
1774:     if (mat) *mat = NULL;
1775:   } else {
1776:     osm = (PC_ASM*)pc->data;
1777:     if (n) *n = osm->n_local_true;
1778:     if (mat) *mat = osm->pmat;
1779:   }
1780:   return(0);
1781: }

1785: /*@
1786:     PCASMSetDMSubdomains - Indicates whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1787:     Logically Collective

1789:     Input Parameter:
1790: +   pc  - the preconditioner
1791: -   flg - boolean indicating whether to use subdomains defined by the DM

1793:     Options Database Key:
1794: .   -pc_asm_dm_subdomains

1796:     Level: intermediate

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

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

1804: .seealso: PCASMGetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1805:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1806: @*/
1807: PetscErrorCode  PCASMSetDMSubdomains(PC pc,PetscBool flg)
1808: {
1809:   PC_ASM         *osm = (PC_ASM*)pc->data;
1811:   PetscBool      match;

1816:   if (pc->setupcalled) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for a setup PC.");
1817:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1818:   if (match) {
1819:     osm->dm_subdomains = flg;
1820:   }
1821:   return(0);
1822: }

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

1830:     Input Parameter:
1831: .   pc  - the preconditioner

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

1836:     Level: intermediate

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

1840: .seealso: PCASMSetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1841:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1842: @*/
1843: PetscErrorCode  PCASMGetDMSubdomains(PC pc,PetscBool* flg)
1844: {
1845:   PC_ASM         *osm = (PC_ASM*)pc->data;
1847:   PetscBool      match;

1852:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1853:   if (match) {
1854:     if (flg) *flg = osm->dm_subdomains;
1855:   }
1856:   return(0);
1857: }


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

1865:    Not Collective

1867:    Input Parameter:
1868: .  pc - the PC

1870:    Output Parameter:
1871: .  sub_mat_type - name of matrix type

1873:    Level: advanced

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

1877: .seealso: PCASMSetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1878: @*/
1879: PetscErrorCode  PCASMGetSubMatType(PC pc,MatType *sub_mat_type){

1882:   PetscTryMethod(pc,"PCASMGetSubMatType_C",(PC,MatType*),(pc,sub_mat_type));
1883:   return(0);
1884: }

1888: /*@
1889:  PCASMSetSubMatType - Set the type of matrix used for ASM subsolves

1891:    Collective on Mat

1893:    Input Parameters:
1894: +  pc             - the PC object
1895: -  sub_mat_type   - matrix type

1897:    Options Database Key:
1898: .  -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.

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

1903:   Level: advanced

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

1907: .seealso: PCASMGetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1908: @*/
1909: PetscErrorCode PCASMSetSubMatType(PC pc,MatType sub_mat_type)
1910: {

1913:   PetscTryMethod(pc,"PCASMSetSubMatType_C",(PC,MatType),(pc,sub_mat_type));
1914:   return(0);
1915: }