Actual source code: asm.c

petsc-master 2015-02-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>     /*I "petscpc.h" I*/
 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: } PC_ASM;

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

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

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

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

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

188:   if (!pc->setupcalled) {

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

195:     /* Note: if subdomains have been set either via PCASMSetTotalSubdomains() or via PCASMSetLocalSubdomains(), osm->n_local_true will not be PETSC_DECIDE */
196:     if (osm->n_local_true == PETSC_DECIDE) {
197:       /* no subdomains given */
198:       /* try pc->dm first, if allowed */
199:       if (osm->dm_subdomains && pc->dm) {
200:         PetscInt  num_domains, d;
201:         char      **domain_names;
202:         IS        *inner_domain_is, *outer_domain_is;
203:         DMCreateDomainDecomposition(pc->dm, &num_domains, &domain_names, &inner_domain_is, &outer_domain_is, &domain_dm);
204:         if (num_domains) {
205:           PCASMSetLocalSubdomains(pc, num_domains, outer_domain_is, inner_domain_is);
206:         }
207:         for (d = 0; d < num_domains; ++d) {
208:           if (domain_names)    {PetscFree(domain_names[d]);}
209:           if (inner_domain_is) {ISDestroy(&inner_domain_is[d]);}
210:           if (outer_domain_is) {ISDestroy(&outer_domain_is[d]);}
211:         }
212:         PetscFree(domain_names);
213:         PetscFree(inner_domain_is);
214:         PetscFree(outer_domain_is);
215:       }
216:       if (osm->n_local_true == PETSC_DECIDE) {
217:         /* still no subdomains; use one subdomain per processor */
218:         osm->n_local_true = 1;
219:       }
220:     }
221:     { /* determine the global and max number of subdomains */
222:       struct {PetscInt max,sum;} inwork,outwork;
223:       inwork.max   = osm->n_local_true;
224:       inwork.sum   = osm->n_local_true;
225:       MPI_Allreduce(&inwork,&outwork,1,MPIU_2INT,PetscMaxSum_Op,PetscObjectComm((PetscObject)pc));
226:       osm->n_local = outwork.max;
227:       osm->n       = outwork.sum;
228:     }
229:     if (!osm->is) { /* create the index sets */
230:       PCASMCreateSubdomains(pc->pmat,osm->n_local_true,&osm->is);
231:     }
232:     if (osm->n_local_true > 1 && !osm->is_local) {
233:       PetscMalloc1(osm->n_local_true,&osm->is_local);
234:       for (i=0; i<osm->n_local_true; i++) {
235:         if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
236:           ISDuplicate(osm->is[i],&osm->is_local[i]);
237:           ISCopy(osm->is[i],osm->is_local[i]);
238:         } else {
239:           PetscObjectReference((PetscObject)osm->is[i]);
240:           osm->is_local[i] = osm->is[i];
241:         }
242:       }
243:     }
244:     PCGetOptionsPrefix(pc,&prefix);
245:     flg  = PETSC_FALSE;
246:     PetscOptionsGetBool(prefix,"-pc_asm_print_subdomains",&flg,NULL);
247:     if (flg) { PCASMPrintSubdomains(pc); }

249:     if (osm->overlap > 0) {
250:       /* Extend the "overlapping" regions by a number of steps */
251:       MatIncreaseOverlap(pc->pmat,osm->n_local_true,osm->is,osm->overlap);
252:     }
253:     if (osm->sort_indices) {
254:       for (i=0; i<osm->n_local_true; i++) {
255:         ISSort(osm->is[i]);
256:         if (osm->is_local) {
257:           ISSort(osm->is_local[i]);
258:         }
259:       }
260:     }
261:     /* Create the local work vectors and scatter contexts */
262:     MatCreateVecs(pc->pmat,&vec,0);
263:     PetscMalloc1(osm->n_local,&osm->restriction);
264:     if (osm->is_local) {PetscMalloc1(osm->n_local,&osm->localization);}
265:     PetscMalloc1(osm->n_local,&osm->prolongation);
266:     PetscMalloc1(osm->n_local,&osm->x);
267:     PetscMalloc1(osm->n_local,&osm->y);
268:     PetscMalloc1(osm->n_local,&osm->y_local);
269:     for (i=0; i<osm->n_local_true; ++i) {
270:       ISGetLocalSize(osm->is[i],&m);
271:       VecCreateSeq(PETSC_COMM_SELF,m,&osm->x[i]);
272:       ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
273:       VecScatterCreate(vec,osm->is[i],osm->x[i],isl,&osm->restriction[i]);
274:       ISDestroy(&isl);
275:       VecDuplicate(osm->x[i],&osm->y[i]);
276:       if (osm->is_local) {
277:         ISLocalToGlobalMapping ltog;
278:         IS                     isll;
279:         const PetscInt         *idx_local;
280:         PetscInt               *idx,nout;

282:         ISLocalToGlobalMappingCreateIS(osm->is[i],&ltog);
283:         ISGetLocalSize(osm->is_local[i],&m_local);
284:         ISGetIndices(osm->is_local[i], &idx_local);
285:         PetscMalloc1(m_local,&idx);
286:         ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx);
287:         ISLocalToGlobalMappingDestroy(&ltog);
288:         if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of is");
289:         ISRestoreIndices(osm->is_local[i], &idx_local);
290:         ISCreateGeneral(PETSC_COMM_SELF,m_local,idx,PETSC_OWN_POINTER,&isll);
291:         ISCreateStride(PETSC_COMM_SELF,m_local,0,1,&isl);
292:         VecCreateSeq(PETSC_COMM_SELF,m_local,&osm->y_local[i]);
293:         VecScatterCreate(osm->y[i],isll,osm->y_local[i],isl,&osm->localization[i]);
294:         ISDestroy(&isll);

296:         VecScatterCreate(vec,osm->is_local[i],osm->y_local[i],isl,&osm->prolongation[i]);
297:         ISDestroy(&isl);
298:       } else {
299:         VecGetLocalSize(vec,&m_local);

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

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

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

307:         PetscObjectReference((PetscObject) osm->restriction[i]);
308:       }
309:     }
310:     for (i=osm->n_local_true; i<osm->n_local; i++) {
311:       VecCreateSeq(PETSC_COMM_SELF,0,&osm->x[i]);
312:       VecDuplicate(osm->x[i],&osm->y[i]);
313:       VecDuplicate(osm->x[i],&osm->y_local[i]);
314:       ISCreateStride(PETSC_COMM_SELF,0,0,1,&isl);
315:       VecScatterCreate(vec,isl,osm->x[i],isl,&osm->restriction[i]);
316:       if (osm->is_local) {
317:         VecScatterCreate(osm->y[i],isl,osm->y_local[i],isl,&osm->localization[i]);
318:         VecScatterCreate(vec,isl,osm->x[i],isl,&osm->prolongation[i]);
319:       } else {
320:         osm->prolongation[i] = osm->restriction[i];
321:         PetscObjectReference((PetscObject) osm->restriction[i]);
322:       }
323:       ISDestroy(&isl);
324:     }
325:     VecDestroy(&vec);

327:     if (!osm->ksp) {
328:       /* Create the local solvers */
329:       PetscMalloc1(osm->n_local_true,&osm->ksp);
330:       if (domain_dm) {
331:         PetscInfo(pc,"Setting up ASM subproblems using the embedded DM\n");
332:       }
333:       for (i=0; i<osm->n_local_true; i++) {
334:         KSPCreate(PETSC_COMM_SELF,&ksp);
335:         PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
336:         PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
337:         KSPSetType(ksp,KSPPREONLY);
338:         KSPGetPC(ksp,&subpc);
339:         PCGetOptionsPrefix(pc,&prefix);
340:         KSPSetOptionsPrefix(ksp,prefix);
341:         KSPAppendOptionsPrefix(ksp,"sub_");
342:         if (domain_dm) {
343:           KSPSetDM(ksp, domain_dm[i]);
344:           KSPSetDMActive(ksp, PETSC_FALSE);
345:           DMDestroy(&domain_dm[i]);
346:         }
347:         osm->ksp[i] = ksp;
348:       }
349:       if (domain_dm) {
350:         PetscFree(domain_dm);
351:       }
352:     }
353:     scall = MAT_INITIAL_MATRIX;
354:   } else {
355:     /*
356:        Destroy the blocks from the previous iteration
357:     */
358:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
359:       MatDestroyMatrices(osm->n_local_true,&osm->pmat);
360:       scall = MAT_INITIAL_MATRIX;
361:     }
362:   }

364:   /*
365:      Extract out the submatrices
366:   */
367:   MatGetSubMatrices(pc->pmat,osm->n_local_true,osm->is,osm->is,scall,&osm->pmat);
368:   if (scall == MAT_INITIAL_MATRIX) {
369:     PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
370:     for (i=0; i<osm->n_local_true; i++) {
371:       PetscLogObjectParent((PetscObject)pc,(PetscObject)osm->pmat[i]);
372:       PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);
373:     }
374:   }

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

380:   /*
381:      Loop over subdomains putting them into local ksp
382:   */
383:   for (i=0; i<osm->n_local_true; i++) {
384:     KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i]);
385:     if (!pc->setupcalled) {
386:       KSPSetFromOptions(osm->ksp[i]);
387:     }
388:   }
389:   return(0);
390: }

394: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
395: {
396:   PC_ASM         *osm = (PC_ASM*)pc->data;
398:   PetscInt       i;

401:   for (i=0; i<osm->n_local_true; i++) {
402:     KSPSetUp(osm->ksp[i]);
403:   }
404:   return(0);
405: }

409: static PetscErrorCode PCApply_ASM(PC pc,Vec x,Vec y)
410: {
411:   PC_ASM         *osm = (PC_ASM*)pc->data;
413:   PetscInt       i,n_local = osm->n_local,n_local_true = osm->n_local_true;
414:   ScatterMode    forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;

417:   /*
418:      Support for limiting the restriction or interpolation to only local
419:      subdomain values (leaving the other values 0).
420:   */
421:   if (!(osm->type & PC_ASM_RESTRICT)) {
422:     forward = SCATTER_FORWARD_LOCAL;
423:     /* have to zero the work RHS since scatter may leave some slots empty */
424:     for (i=0; i<n_local_true; i++) {
425:       VecZeroEntries(osm->x[i]);
426:     }
427:   }
428:   if (!(osm->type & PC_ASM_INTERPOLATE)) reverse = SCATTER_REVERSE_LOCAL;

430:   for (i=0; i<n_local; i++) {
431:     VecScatterBegin(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
432:   }
433:   VecZeroEntries(y);
434:   /* do the local solves */
435:   for (i=0; i<n_local_true; i++) {
436:     VecScatterEnd(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
437:     KSPSolve(osm->ksp[i],osm->x[i],osm->y[i]);
438:     if (osm->localization) {
439:       VecScatterBegin(osm->localization[i],osm->y[i],osm->y_local[i],INSERT_VALUES,forward);
440:       VecScatterEnd(osm->localization[i],osm->y[i],osm->y_local[i],INSERT_VALUES,forward);
441:     }
442:     VecScatterBegin(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
443:   }
444:   /* handle the rest of the scatters that do not have local solves */
445:   for (i=n_local_true; i<n_local; i++) {
446:     VecScatterEnd(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
447:     VecScatterBegin(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
448:   }
449:   for (i=0; i<n_local; i++) {
450:     VecScatterEnd(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
451:   }
452:   return(0);
453: }

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

465:   /*
466:      Support for limiting the restriction or interpolation to only local
467:      subdomain values (leaving the other values 0).

469:      Note: these are reversed from the PCApply_ASM() because we are applying the
470:      transpose of the three terms
471:   */
472:   if (!(osm->type & PC_ASM_INTERPOLATE)) {
473:     forward = SCATTER_FORWARD_LOCAL;
474:     /* have to zero the work RHS since scatter may leave some slots empty */
475:     for (i=0; i<n_local_true; i++) {
476:       VecZeroEntries(osm->x[i]);
477:     }
478:   }
479:   if (!(osm->type & PC_ASM_RESTRICT)) reverse = SCATTER_REVERSE_LOCAL;

481:   for (i=0; i<n_local; i++) {
482:     VecScatterBegin(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
483:   }
484:   VecZeroEntries(y);
485:   /* do the local solves */
486:   for (i=0; i<n_local_true; i++) {
487:     VecScatterEnd(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
488:     KSPSolveTranspose(osm->ksp[i],osm->x[i],osm->y[i]);
489:     if (osm->localization) {
490:       VecScatterBegin(osm->localization[i],osm->y[i],osm->y_local[i],INSERT_VALUES,forward);
491:       VecScatterEnd(osm->localization[i],osm->y[i],osm->y_local[i],INSERT_VALUES,forward);
492:     }
493:     VecScatterBegin(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
494:   }
495:   /* handle the rest of the scatters that do not have local solves */
496:   for (i=n_local_true; i<n_local; i++) {
497:     VecScatterEnd(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
498:     VecScatterBegin(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
499:   }
500:   for (i=0; i<n_local; i++) {
501:     VecScatterEnd(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
502:   }
503:   return(0);
504: }

508: static PetscErrorCode PCReset_ASM(PC pc)
509: {
510:   PC_ASM         *osm = (PC_ASM*)pc->data;
512:   PetscInt       i;

515:   if (osm->ksp) {
516:     for (i=0; i<osm->n_local_true; i++) {
517:       KSPReset(osm->ksp[i]);
518:     }
519:   }
520:   if (osm->pmat) {
521:     if (osm->n_local_true > 0) {
522:       MatDestroyMatrices(osm->n_local_true,&osm->pmat);
523:     }
524:   }
525:   if (osm->restriction) {
526:     for (i=0; i<osm->n_local; i++) {
527:       VecScatterDestroy(&osm->restriction[i]);
528:       if (osm->localization) {VecScatterDestroy(&osm->localization[i]);}
529:       VecScatterDestroy(&osm->prolongation[i]);
530:       VecDestroy(&osm->x[i]);
531:       VecDestroy(&osm->y[i]);
532:       VecDestroy(&osm->y_local[i]);
533:     }
534:     PetscFree(osm->restriction);
535:     if (osm->localization) {PetscFree(osm->localization);}
536:     PetscFree(osm->prolongation);
537:     PetscFree(osm->x);
538:     PetscFree(osm->y);
539:     PetscFree(osm->y_local);
540:   }
541:   PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);

543:   osm->is       = 0;
544:   osm->is_local = 0;
545:   return(0);
546: }

550: static PetscErrorCode PCDestroy_ASM(PC pc)
551: {
552:   PC_ASM         *osm = (PC_ASM*)pc->data;
554:   PetscInt       i;

557:   PCReset_ASM(pc);
558:   if (osm->ksp) {
559:     for (i=0; i<osm->n_local_true; i++) {
560:       KSPDestroy(&osm->ksp[i]);
561:     }
562:     PetscFree(osm->ksp);
563:   }
564:   PetscFree(pc->data);
565:   return(0);
566: }

570: static PetscErrorCode PCSetFromOptions_ASM(PetscOptions *PetscOptionsObject,PC pc)
571: {
572:   PC_ASM         *osm = (PC_ASM*)pc->data;
574:   PetscInt       blocks,ovl;
575:   PetscBool      symset,flg;
576:   PCASMType      asmtype;

579:   /* set the type to symmetric if matrix is symmetric */
580:   if (!osm->type_set && pc->pmat) {
581:     MatIsSymmetricKnown(pc->pmat,&symset,&flg);
582:     if (symset && flg) osm->type = PC_ASM_BASIC;
583:   }
584:   PetscOptionsHead(PetscOptionsObject,"Additive Schwarz options");
585:   PetscOptionsBool("-pc_asm_dm_subdomains","Use DMCreateDomainDecomposition() to define subdomains","PCASMSetDMSubdomains",osm->dm_subdomains,&osm->dm_subdomains,&flg);
586:   PetscOptionsInt("-pc_asm_blocks","Number of subdomains","PCASMSetTotalSubdomains",osm->n,&blocks,&flg);
587:   if (flg) {
588:     PCASMSetTotalSubdomains(pc,blocks,NULL,NULL);
589:     osm->dm_subdomains = PETSC_FALSE;
590:   }
591:   PetscOptionsInt("-pc_asm_overlap","Number of grid points overlap","PCASMSetOverlap",osm->overlap,&ovl,&flg);
592:   if (flg) {
593:     PCASMSetOverlap(pc,ovl);
594:     osm->dm_subdomains = PETSC_FALSE;
595:   }
596:   flg  = PETSC_FALSE;
597:   PetscOptionsEnum("-pc_asm_type","Type of restriction/extension","PCASMSetType",PCASMTypes,(PetscEnum)osm->type,(PetscEnum*)&asmtype,&flg);
598:   if (flg) {PCASMSetType(pc,asmtype); }
599:   PetscOptionsTail();
600:   return(0);
601: }

603: /*------------------------------------------------------------------------------------*/

607: static PetscErrorCode  PCASMSetLocalSubdomains_ASM(PC pc,PetscInt n,IS is[],IS is_local[])
608: {
609:   PC_ASM         *osm = (PC_ASM*)pc->data;
611:   PetscInt       i;

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

617:   if (!pc->setupcalled) {
618:     if (is) {
619:       for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is[i]);}
620:     }
621:     if (is_local) {
622:       for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is_local[i]);}
623:     }
624:     PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);

626:     osm->n_local_true = n;
627:     osm->is           = 0;
628:     osm->is_local     = 0;
629:     if (is) {
630:       PetscMalloc1(n,&osm->is);
631:       for (i=0; i<n; i++) osm->is[i] = is[i];
632:       /* Flag indicating that the user has set overlapping subdomains so PCASM should not increase their size. */
633:       osm->overlap = -1;
634:     }
635:     if (is_local) {
636:       PetscMalloc1(n,&osm->is_local);
637:       for (i=0; i<n; i++) osm->is_local[i] = is_local[i];
638:       if (!is) {
639:         PetscMalloc1(osm->n_local_true,&osm->is);
640:         for (i=0; i<osm->n_local_true; i++) {
641:           if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
642:             ISDuplicate(osm->is_local[i],&osm->is[i]);
643:             ISCopy(osm->is_local[i],osm->is[i]);
644:           } else {
645:             PetscObjectReference((PetscObject)osm->is_local[i]);
646:             osm->is[i] = osm->is_local[i];
647:           }
648:         }
649:       }
650:     }
651:   }
652:   return(0);
653: }

657: static PetscErrorCode  PCASMSetTotalSubdomains_ASM(PC pc,PetscInt N,IS *is,IS *is_local)
658: {
659:   PC_ASM         *osm = (PC_ASM*)pc->data;
661:   PetscMPIInt    rank,size;
662:   PetscInt       n;

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

668:   /*
669:      Split the subdomains equally among all processors
670:   */
671:   MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
672:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
673:   n    = N/size + ((N % size) > rank);
674:   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);
675:   if (pc->setupcalled && n != osm->n_local_true) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCASMSetTotalSubdomains() should be called before PCSetUp().");
676:   if (!pc->setupcalled) {
677:     PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);

679:     osm->n_local_true = n;
680:     osm->is           = 0;
681:     osm->is_local     = 0;
682:   }
683:   return(0);
684: }

688: static PetscErrorCode  PCASMSetOverlap_ASM(PC pc,PetscInt ovl)
689: {
690:   PC_ASM *osm = (PC_ASM*)pc->data;

693:   if (ovl < 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested");
694:   if (pc->setupcalled && ovl != osm->overlap) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetOverlap() should be called before PCSetUp().");
695:   if (!pc->setupcalled) osm->overlap = ovl;
696:   return(0);
697: }

701: static PetscErrorCode  PCASMSetType_ASM(PC pc,PCASMType type)
702: {
703:   PC_ASM *osm = (PC_ASM*)pc->data;

706:   osm->type     = type;
707:   osm->type_set = PETSC_TRUE;
708:   return(0);
709: }

713: static PetscErrorCode  PCASMGetType_ASM(PC pc,PCASMType *type)
714: {
715:   PC_ASM *osm = (PC_ASM*)pc->data;

718:   *type = osm->type;
719:   return(0);
720: }

724: static PetscErrorCode  PCASMSetSortIndices_ASM(PC pc,PetscBool  doSort)
725: {
726:   PC_ASM *osm = (PC_ASM*)pc->data;

729:   osm->sort_indices = doSort;
730:   return(0);
731: }

735: static PetscErrorCode  PCASMGetSubKSP_ASM(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
736: {
737:   PC_ASM         *osm = (PC_ASM*)pc->data;

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

743:   if (n_local) *n_local = osm->n_local_true;
744:   if (first_local) {
745:     MPI_Scan(&osm->n_local_true,first_local,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
746:     *first_local -= osm->n_local_true;
747:   }
748:   if (ksp) {
749:     /* Assume that local solves are now different; not necessarily
750:        true though!  This flag is used only for PCView_ASM() */
751:     *ksp                   = osm->ksp;
752:     osm->same_local_solves = PETSC_FALSE;
753:   }
754:   return(0);
755: }

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

762:     Collective on PC

764:     Input Parameters:
765: +   pc - the preconditioner context
766: .   n - the number of subdomains for this processor (default value = 1)
767: .   is - the index set that defines the subdomains for this processor
768:          (or NULL for PETSc to determine subdomains)
769: -   is_local - the index sets that define the local part of the subdomains for this processor
770:          (or NULL to use the default of 1 subdomain per process)

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

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

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

779:     Level: advanced

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

783: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
784:           PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains()
785: @*/
786: PetscErrorCode  PCASMSetLocalSubdomains(PC pc,PetscInt n,IS is[],IS is_local[])
787: {

792:   PetscTryMethod(pc,"PCASMSetLocalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,is,is_local));
793:   return(0);
794: }

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

803:     Collective on PC

805:     Input Parameters:
806: +   pc - the preconditioner context
807: .   N  - the number of subdomains for all processors
808: .   is - the index sets that define the subdomains for all processors
809:          (or NULL to ask PETSc to compe up with subdomains)
810: -   is_local - the index sets that define the local part of the subdomains for this processor
811:          (or NULL to use the default of 1 subdomain per process)

813:     Options Database Key:
814:     To set the total number of subdomain blocks rather than specify the
815:     index sets, use the option
816: .    -pc_asm_blocks <blks> - Sets total blocks

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

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

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

826:     Use PCASMSetLocalSubdomains() to set local subdomains.

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

830:     Level: advanced

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

834: .seealso: PCASMSetLocalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
835:           PCASMCreateSubdomains2D()
836: @*/
837: PetscErrorCode  PCASMSetTotalSubdomains(PC pc,PetscInt N,IS is[],IS is_local[])
838: {

843:   PetscTryMethod(pc,"PCASMSetTotalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,N,is,is_local));
844:   return(0);
845: }

849: /*@
850:     PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
851:     additive Schwarz preconditioner.  Either all or no processors in the
852:     PC communicator must call this routine.

854:     Logically Collective on PC

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

860:     Options Database Key:
861: .   -pc_asm_overlap <ovl> - Sets overlap

863:     Notes:
864:     By default the ASM preconditioner uses 1 block per processor.  To use
865:     multiple blocks per perocessor, see PCASMSetTotalSubdomains() and
866:     PCASMSetLocalSubdomains() (and the option -pc_asm_blocks <blks>).

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

876:     Note that one can define initial index sets with any overlap via
877:     PCASMSetTotalSubdomains() or PCASMSetLocalSubdomains(); the routine
878:     PCASMSetOverlap() merely allows PETSc to extend that overlap further
879:     if desired.

881:     Level: intermediate

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

885: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
886:           PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains()
887: @*/
888: PetscErrorCode  PCASMSetOverlap(PC pc,PetscInt ovl)
889: {

895:   PetscTryMethod(pc,"PCASMSetOverlap_C",(PC,PetscInt),(pc,ovl));
896:   return(0);
897: }

901: /*@
902:     PCASMSetType - Sets the type of restriction and interpolation used
903:     for local problems in the additive Schwarz method.

905:     Logically Collective on PC

907:     Input Parameters:
908: +   pc  - the preconditioner context
909: -   type - variant of ASM, one of
910: .vb
911:       PC_ASM_BASIC       - full interpolation and restriction
912:       PC_ASM_RESTRICT    - full restriction, local processor interpolation
913:       PC_ASM_INTERPOLATE - full interpolation, local processor restriction
914:       PC_ASM_NONE        - local processor restriction and interpolation
915: .ve

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

920:     Level: intermediate

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

924: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
925:           PCASMCreateSubdomains2D()
926: @*/
927: PetscErrorCode  PCASMSetType(PC pc,PCASMType type)
928: {

934:   PetscTryMethod(pc,"PCASMSetType_C",(PC,PCASMType),(pc,type));
935:   return(0);
936: }

940: /*@
941:     PCASMGetType - Gets the type of restriction and interpolation used
942:     for local problems in the additive Schwarz method.

944:     Logically Collective on PC

946:     Input Parameter:
947: .   pc  - the preconditioner context

949:     Output Parameter:
950: .   type - variant of ASM, one of

952: .vb
953:       PC_ASM_BASIC       - full interpolation and restriction
954:       PC_ASM_RESTRICT    - full restriction, local processor interpolation
955:       PC_ASM_INTERPOLATE - full interpolation, local processor restriction
956:       PC_ASM_NONE        - local processor restriction and interpolation
957: .ve

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

962:     Level: intermediate

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

966: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
967:           PCASMCreateSubdomains2D()
968: @*/
969: PetscErrorCode  PCASMGetType(PC pc,PCASMType *type)
970: {

975:   PetscUseMethod(pc,"PCASMGetType_C",(PC,PCASMType*),(pc,type));
976:   return(0);
977: }

981: /*@
982:     PCASMSetSortIndices - Determines whether subdomain indices are sorted.

984:     Logically Collective on PC

986:     Input Parameters:
987: +   pc  - the preconditioner context
988: -   doSort - sort the subdomain indices

990:     Level: intermediate

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

994: .seealso: PCASMSetLocalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
995:           PCASMCreateSubdomains2D()
996: @*/
997: PetscErrorCode  PCASMSetSortIndices(PC pc,PetscBool doSort)
998: {

1004:   PetscTryMethod(pc,"PCASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));
1005:   return(0);
1006: }

1010: /*@C
1011:    PCASMGetSubKSP - Gets the local KSP contexts for all blocks on
1012:    this processor.

1014:    Collective on PC iff first_local is requested

1016:    Input Parameter:
1017: .  pc - the preconditioner context

1019:    Output Parameters:
1020: +  n_local - the number of blocks on this processor or NULL
1021: .  first_local - the global number of the first block on this processor or NULL,
1022:                  all processors must request or all must pass NULL
1023: -  ksp - the array of KSP contexts

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

1028:    Currently for some matrix implementations only 1 block per processor
1029:    is supported.

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

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

1036:    Level: advanced

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

1040: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap(),
1041:           PCASMCreateSubdomains2D(),
1042: @*/
1043: PetscErrorCode  PCASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
1044: {

1049:   PetscUseMethod(pc,"PCASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
1050:   return(0);
1051: }

1053: /* -------------------------------------------------------------------------------------*/
1054: /*MC
1055:    PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1056:            its own KSP object.

1058:    Options Database Keys:
1059: +  -pc_asm_blocks <blks> - Sets total blocks
1060: .  -pc_asm_overlap <ovl> - Sets overlap
1061: -  -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type

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

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

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

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


1078:    Level: beginner

1080:    Concepts: additive Schwarz method

1082:     References:
1083:     An additive variant of the Schwarz alternating method for the case of many subregions
1084:     M Dryja, OB Widlund - Courant Institute, New York University Technical report

1086:     Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations,
1087:     Barry Smith, Petter Bjorstad, and William Gropp, Cambridge University Press, ISBN 0-521-49589-X.

1089: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
1090:            PCBJACOBI, PCASMGetSubKSP(), PCASMSetLocalSubdomains(),
1091:            PCASMSetTotalSubdomains(), PCSetModifySubmatrices(), PCASMSetOverlap(), PCASMSetType()

1093: M*/

1097: PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1098: {
1100:   PC_ASM         *osm;

1103:   PetscNewLog(pc,&osm);

1105:   osm->n                 = PETSC_DECIDE;
1106:   osm->n_local           = 0;
1107:   osm->n_local_true      = PETSC_DECIDE;
1108:   osm->overlap           = 1;
1109:   osm->ksp               = 0;
1110:   osm->restriction       = 0;
1111:   osm->localization      = 0;
1112:   osm->prolongation      = 0;
1113:   osm->x                 = 0;
1114:   osm->y                 = 0;
1115:   osm->y_local           = 0;
1116:   osm->is                = 0;
1117:   osm->is_local          = 0;
1118:   osm->mat               = 0;
1119:   osm->pmat              = 0;
1120:   osm->type              = PC_ASM_RESTRICT;
1121:   osm->same_local_solves = PETSC_TRUE;
1122:   osm->sort_indices      = PETSC_TRUE;
1123:   osm->dm_subdomains     = PETSC_FALSE;

1125:   pc->data                 = (void*)osm;
1126:   pc->ops->apply           = PCApply_ASM;
1127:   pc->ops->applytranspose  = PCApplyTranspose_ASM;
1128:   pc->ops->setup           = PCSetUp_ASM;
1129:   pc->ops->reset           = PCReset_ASM;
1130:   pc->ops->destroy         = PCDestroy_ASM;
1131:   pc->ops->setfromoptions  = PCSetFromOptions_ASM;
1132:   pc->ops->setuponblocks   = PCSetUpOnBlocks_ASM;
1133:   pc->ops->view            = PCView_ASM;
1134:   pc->ops->applyrichardson = 0;

1136:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",PCASMSetLocalSubdomains_ASM);
1137:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",PCASMSetTotalSubdomains_ASM);
1138:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetOverlap_C",PCASMSetOverlap_ASM);
1139:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetType_C",PCASMSetType_ASM);
1140:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetType_C",PCASMGetType_ASM);
1141:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSortIndices_C",PCASMSetSortIndices_ASM);
1142:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubKSP_C",PCASMGetSubKSP_ASM);
1143:   return(0);
1144: }

1148: /*@C
1149:    PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1150:    preconditioner for a any problem on a general grid.

1152:    Collective

1154:    Input Parameters:
1155: +  A - The global matrix operator
1156: -  n - the number of local blocks

1158:    Output Parameters:
1159: .  outis - the array of index sets defining the subdomains

1161:    Level: advanced

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

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

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

1170: .seealso: PCASMSetLocalSubdomains(), PCASMDestroySubdomains()
1171: @*/
1172: PetscErrorCode  PCASMCreateSubdomains(Mat A, PetscInt n, IS* outis[])
1173: {
1174:   MatPartitioning mpart;
1175:   const char      *prefix;
1176:   PetscErrorCode  (*f)(Mat,Mat*);
1177:   PetscMPIInt     size;
1178:   PetscInt        i,j,rstart,rend,bs;
1179:   PetscBool       isbaij = PETSC_FALSE,foundpart = PETSC_FALSE;
1180:   Mat             Ad     = NULL, adj;
1181:   IS              ispart,isnumb,*is;
1182:   PetscErrorCode  ierr;

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

1189:   /* Get prefix, row distribution, and block size */
1190:   MatGetOptionsPrefix(A,&prefix);
1191:   MatGetOwnershipRange(A,&rstart,&rend);
1192:   MatGetBlockSize(A,&bs);
1193:   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);

1195:   /* Get diagonal block from matrix if possible */
1196:   MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1197:   PetscObjectQueryFunction((PetscObject)A,"MatGetDiagonalBlock_C",&f);
1198:   if (f) {
1199:     MatGetDiagonalBlock(A,&Ad);
1200:   } else if (size == 1) {
1201:     Ad = A;
1202:   }
1203:   if (Ad) {
1204:     PetscObjectTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij);
1205:     if (!isbaij) {PetscObjectTypeCompare((PetscObject)Ad,MATSEQSBAIJ,&isbaij);}
1206:   }
1207:   if (Ad && n > 1) {
1208:     PetscBool match,done;
1209:     /* Try to setup a good matrix partitioning if available */
1210:     MatPartitioningCreate(PETSC_COMM_SELF,&mpart);
1211:     PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);
1212:     MatPartitioningSetFromOptions(mpart);
1213:     PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGCURRENT,&match);
1214:     if (!match) {
1215:       PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGSQUARE,&match);
1216:     }
1217:     if (!match) { /* assume a "good" partitioner is available */
1218:       PetscInt       na;
1219:       const PetscInt *ia,*ja;
1220:       MatGetRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1221:       if (done) {
1222:         /* Build adjacency matrix by hand. Unfortunately a call to
1223:            MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1224:            remove the block-aij structure and we cannot expect
1225:            MatPartitioning to split vertices as we need */
1226:         PetscInt       i,j,len,nnz,cnt,*iia=0,*jja=0;
1227:         const PetscInt *row;
1228:         nnz = 0;
1229:         for (i=0; i<na; i++) { /* count number of nonzeros */
1230:           len = ia[i+1] - ia[i];
1231:           row = ja + ia[i];
1232:           for (j=0; j<len; j++) {
1233:             if (row[j] == i) { /* don't count diagonal */
1234:               len--; break;
1235:             }
1236:           }
1237:           nnz += len;
1238:         }
1239:         PetscMalloc1(na+1,&iia);
1240:         PetscMalloc1(nnz,&jja);
1241:         nnz    = 0;
1242:         iia[0] = 0;
1243:         for (i=0; i<na; i++) { /* fill adjacency */
1244:           cnt = 0;
1245:           len = ia[i+1] - ia[i];
1246:           row = ja + ia[i];
1247:           for (j=0; j<len; j++) {
1248:             if (row[j] != i) { /* if not diagonal */
1249:               jja[nnz+cnt++] = row[j];
1250:             }
1251:           }
1252:           nnz     += cnt;
1253:           iia[i+1] = nnz;
1254:         }
1255:         /* Partitioning of the adjacency matrix */
1256:         MatCreateMPIAdj(PETSC_COMM_SELF,na,na,iia,jja,NULL,&adj);
1257:         MatPartitioningSetAdjacency(mpart,adj);
1258:         MatPartitioningSetNParts(mpart,n);
1259:         MatPartitioningApply(mpart,&ispart);
1260:         ISPartitioningToNumbering(ispart,&isnumb);
1261:         MatDestroy(&adj);
1262:         foundpart = PETSC_TRUE;
1263:       }
1264:       MatRestoreRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1265:     }
1266:     MatPartitioningDestroy(&mpart);
1267:   }

1269:   PetscMalloc1(n,&is);
1270:   *outis = is;

1272:   if (!foundpart) {

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

1276:     PetscInt mbs   = (rend-rstart)/bs;
1277:     PetscInt start = rstart;
1278:     for (i=0; i<n; i++) {
1279:       PetscInt count = (mbs/n + ((mbs % n) > i)) * bs;
1280:       ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);
1281:       start += count;
1282:     }

1284:   } else {

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

1288:     const PetscInt *numbering;
1289:     PetscInt       *count,nidx,*indices,*newidx,start=0;
1290:     /* Get node count in each partition */
1291:     PetscMalloc1(n,&count);
1292:     ISPartitioningCount(ispart,n,count);
1293:     if (isbaij && bs > 1) { /* adjust for the block-aij case */
1294:       for (i=0; i<n; i++) count[i] *= bs;
1295:     }
1296:     /* Build indices from node numbering */
1297:     ISGetLocalSize(isnumb,&nidx);
1298:     PetscMalloc1(nidx,&indices);
1299:     for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */
1300:     ISGetIndices(isnumb,&numbering);
1301:     PetscSortIntWithPermutation(nidx,numbering,indices);
1302:     ISRestoreIndices(isnumb,&numbering);
1303:     if (isbaij && bs > 1) { /* adjust for the block-aij case */
1304:       PetscMalloc1(nidx*bs,&newidx);
1305:       for (i=0; i<nidx; i++) {
1306:         for (j=0; j<bs; j++) newidx[i*bs+j] = indices[i]*bs + j;
1307:       }
1308:       PetscFree(indices);
1309:       nidx   *= bs;
1310:       indices = newidx;
1311:     }
1312:     /* Shift to get global indices */
1313:     for (i=0; i<nidx; i++) indices[i] += rstart;

1315:     /* Build the index sets for each block */
1316:     for (i=0; i<n; i++) {
1317:       ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);
1318:       ISSort(is[i]);
1319:       start += count[i];
1320:     }

1322:     PetscFree(count);
1323:     PetscFree(indices);
1324:     ISDestroy(&isnumb);
1325:     ISDestroy(&ispart);

1327:   }
1328:   return(0);
1329: }

1333: /*@C
1334:    PCASMDestroySubdomains - Destroys the index sets created with
1335:    PCASMCreateSubdomains(). Should be called after setting subdomains
1336:    with PCASMSetLocalSubdomains().

1338:    Collective

1340:    Input Parameters:
1341: +  n - the number of index sets
1342: .  is - the array of index sets
1343: -  is_local - the array of local index sets, can be NULL

1345:    Level: advanced

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

1349: .seealso: PCASMCreateSubdomains(), PCASMSetLocalSubdomains()
1350: @*/
1351: PetscErrorCode  PCASMDestroySubdomains(PetscInt n, IS is[], IS is_local[])
1352: {
1353:   PetscInt       i;

1357:   if (n <= 0) return(0);
1358:   if (is) {
1360:     for (i=0; i<n; i++) { ISDestroy(&is[i]); }
1361:     PetscFree(is);
1362:   }
1363:   if (is_local) {
1365:     for (i=0; i<n; i++) { ISDestroy(&is_local[i]); }
1366:     PetscFree(is_local);
1367:   }
1368:   return(0);
1369: }

1373: /*@
1374:    PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1375:    preconditioner for a two-dimensional problem on a regular grid.

1377:    Not Collective

1379:    Input Parameters:
1380: +  m, n - the number of mesh points in the x and y directions
1381: .  M, N - the number of subdomains in the x and y directions
1382: .  dof - degrees of freedom per node
1383: -  overlap - overlap in mesh lines

1385:    Output Parameters:
1386: +  Nsub - the number of subdomains created
1387: .  is - array of index sets defining overlapping (if overlap > 0) subdomains
1388: -  is_local - array of index sets defining non-overlapping subdomains

1390:    Note:
1391:    Presently PCAMSCreateSubdomains2d() is valid only for sequential
1392:    preconditioners.  More general related routines are
1393:    PCASMSetTotalSubdomains() and PCASMSetLocalSubdomains().

1395:    Level: advanced

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

1399: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1400:           PCASMSetOverlap()
1401: @*/
1402: PetscErrorCode  PCASMCreateSubdomains2D(PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt dof,PetscInt overlap,PetscInt *Nsub,IS **is,IS **is_local)
1403: {
1404:   PetscInt       i,j,height,width,ystart,xstart,yleft,yright,xleft,xright,loc_outer;
1406:   PetscInt       nidx,*idx,loc,ii,jj,count;

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

1411:   *Nsub     = N*M;
1412:   PetscMalloc1(*Nsub,is);
1413:   PetscMalloc1(*Nsub,is_local);
1414:   ystart    = 0;
1415:   loc_outer = 0;
1416:   for (i=0; i<N; i++) {
1417:     height = n/N + ((n % N) > i); /* height of subdomain */
1418:     if (height < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many N subdomains for mesh dimension n");
1419:     yleft  = ystart - overlap; if (yleft < 0) yleft = 0;
1420:     yright = ystart + height + overlap; if (yright > n) yright = n;
1421:     xstart = 0;
1422:     for (j=0; j<M; j++) {
1423:       width = m/M + ((m % M) > j); /* width of subdomain */
1424:       if (width < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many M subdomains for mesh dimension m");
1425:       xleft  = xstart - overlap; if (xleft < 0) xleft = 0;
1426:       xright = xstart + width + overlap; if (xright > m) xright = m;
1427:       nidx   = (xright - xleft)*(yright - yleft);
1428:       PetscMalloc1(nidx,&idx);
1429:       loc    = 0;
1430:       for (ii=yleft; ii<yright; ii++) {
1431:         count = m*ii + xleft;
1432:         for (jj=xleft; jj<xright; jj++) idx[loc++] = count++;
1433:       }
1434:       ISCreateGeneral(PETSC_COMM_SELF,nidx,idx,PETSC_COPY_VALUES,(*is)+loc_outer);
1435:       if (overlap == 0) {
1436:         PetscObjectReference((PetscObject)(*is)[loc_outer]);

1438:         (*is_local)[loc_outer] = (*is)[loc_outer];
1439:       } else {
1440:         for (loc=0,ii=ystart; ii<ystart+height; ii++) {
1441:           for (jj=xstart; jj<xstart+width; jj++) {
1442:             idx[loc++] = m*ii + jj;
1443:           }
1444:         }
1445:         ISCreateGeneral(PETSC_COMM_SELF,loc,idx,PETSC_COPY_VALUES,*is_local+loc_outer);
1446:       }
1447:       PetscFree(idx);
1448:       xstart += width;
1449:       loc_outer++;
1450:     }
1451:     ystart += height;
1452:   }
1453:   for (i=0; i<*Nsub; i++) { ISSort((*is)[i]); }
1454:   return(0);
1455: }

1459: /*@C
1460:     PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
1461:     only) for the additive Schwarz preconditioner.

1463:     Not Collective

1465:     Input Parameter:
1466: .   pc - the preconditioner context

1468:     Output Parameters:
1469: +   n - the number of subdomains for this processor (default value = 1)
1470: .   is - the index sets that define the subdomains for this processor
1471: -   is_local - the index sets that define the local part of the subdomains for this processor (can be NULL)


1474:     Notes:
1475:     The IS numbering is in the parallel, global numbering of the vector.

1477:     Level: advanced

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

1481: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1482:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubmatrices()
1483: @*/
1484: PetscErrorCode  PCASMGetLocalSubdomains(PC pc,PetscInt *n,IS *is[],IS *is_local[])
1485: {
1486:   PC_ASM         *osm;
1488:   PetscBool      match;

1494:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1495:   if (!match) {
1496:     if (n) *n = 0;
1497:     if (is) *is = NULL;
1498:   } else {
1499:     osm = (PC_ASM*)pc->data;
1500:     if (n) *n = osm->n_local_true;
1501:     if (is) *is = osm->is;
1502:     if (is_local) *is_local = osm->is_local;
1503:   }
1504:   return(0);
1505: }

1509: /*@C
1510:     PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1511:     only) for the additive Schwarz preconditioner.

1513:     Not Collective

1515:     Input Parameter:
1516: .   pc - the preconditioner context

1518:     Output Parameters:
1519: +   n - the number of matrices for this processor (default value = 1)
1520: -   mat - the matrices


1523:     Level: advanced

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

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

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

1531: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1532:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains(), PCSetModifySubmatrices()
1533: @*/
1534: PetscErrorCode  PCASMGetLocalSubmatrices(PC pc,PetscInt *n,Mat *mat[])
1535: {
1536:   PC_ASM         *osm;
1538:   PetscBool      match;

1544:   if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUP() or PCSetUp().");
1545:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1546:   if (!match) {
1547:     if (n) *n = 0;
1548:     if (mat) *mat = NULL;
1549:   } else {
1550:     osm = (PC_ASM*)pc->data;
1551:     if (n) *n = osm->n_local_true;
1552:     if (mat) *mat = osm->pmat;
1553:   }
1554:   return(0);
1555: }

1559: /*@
1560:     PCASMSetDMSubdomains - Indicates whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1561:     Logically Collective

1563:     Input Parameter:
1564: +   pc  - the preconditioner
1565: -   flg - boolean indicating whether to use subdomains defined by the DM

1567:     Options Database Key:
1568: .   -pc_asm_dm_subdomains

1570:     Level: intermediate

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

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

1578: .seealso: PCASMGetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1579:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1580: @*/
1581: PetscErrorCode  PCASMSetDMSubdomains(PC pc,PetscBool flg)
1582: {
1583:   PC_ASM         *osm = (PC_ASM*)pc->data;
1585:   PetscBool      match;

1590:   if (pc->setupcalled) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for a setup PC.");
1591:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1592:   if (match) {
1593:     osm->dm_subdomains = flg;
1594:   }
1595:   return(0);
1596: }

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

1604:     Input Parameter:
1605: .   pc  - the preconditioner

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

1610:     Level: intermediate

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

1614: .seealso: PCASMSetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1615:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1616: @*/
1617: PetscErrorCode  PCASMGetDMSubdomains(PC pc,PetscBool* flg)
1618: {
1619:   PC_ASM         *osm = (PC_ASM*)pc->data;
1621:   PetscBool      match;

1626:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1627:   if (match) {
1628:     if (flg) *flg = osm->dm_subdomains;
1629:   }
1630:   return(0);
1631: }