Actual source code: asm.c

  1: /*
  2:   This file defines an additive Schwarz preconditioner for any Mat implementation.

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

  8:        n - total number of true subdomains on all processors
  9:        n_local_true - actual number of subdomains on this processor
 10:        n_local = maximum over all processors of n_local_true
 11: */

 13: #include <petsc/private/pcasmimpl.h>

 15: static PetscErrorCode PCView_ASM(PC pc,PetscViewer viewer)
 16: {
 17:   PC_ASM            *osm = (PC_ASM*)pc->data;
 18:   PetscMPIInt       rank;
 19:   PetscInt          i,bsz;
 20:   PetscBool         iascii,isstring;
 21:   PetscViewer       sviewer;
 22:   PetscViewerFormat format;
 23:   const char        *prefix;

 25:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
 26:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
 27:   if (iascii) {
 28:     char overlaps[256] = "user-defined overlap",blocks[256] = "total subdomain blocks not yet set";
 29:     if (osm->overlap >= 0) PetscSNPrintf(overlaps,sizeof(overlaps),"amount of overlap = %D",osm->overlap);
 30:     if (osm->n > 0) PetscSNPrintf(blocks,sizeof(blocks),"total subdomain blocks = %D",osm->n);
 31:     PetscViewerASCIIPrintf(viewer,"  %s, %s\n",blocks,overlaps);
 32:     PetscViewerASCIIPrintf(viewer,"  restriction/interpolation type - %s\n",PCASMTypes[osm->type]);
 33:     if (osm->dm_subdomains) PetscViewerASCIIPrintf(viewer,"  Additive Schwarz: using DM to define subdomains\n");
 34:     if (osm->loctype != PC_COMPOSITE_ADDITIVE) PetscViewerASCIIPrintf(viewer,"  Additive Schwarz: local solve composition type - %s\n",PCCompositeTypes[osm->loctype]);
 35:     MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
 36:     PetscViewerGetFormat(viewer,&format);
 37:     if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
 38:       if (osm->ksp) {
 39:         PetscViewerASCIIPrintf(viewer,"  Local solver information for first block is in the following KSP and PC objects on rank 0:\n");
 40:         PCGetOptionsPrefix(pc,&prefix);
 41:         PetscViewerASCIIPrintf(viewer,"  Use -%sksp_view ::ascii_info_detail to display information for all blocks\n",prefix?prefix:"");
 42:         PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
 43:         if (rank == 0) {
 44:           PetscViewerASCIIPushTab(viewer);
 45:           KSPView(osm->ksp[0],sviewer);
 46:           PetscViewerASCIIPopTab(viewer);
 47:         }
 48:         PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
 49:       }
 50:     } else {
 51:       PetscViewerASCIIPushSynchronized(viewer);
 52:       PetscViewerASCIISynchronizedPrintf(viewer,"  [%d] number of local blocks = %D\n",(int)rank,osm->n_local_true);
 53:       PetscViewerFlush(viewer);
 54:       PetscViewerASCIIPrintf(viewer,"  Local solver information for each block is in the following KSP and PC objects:\n");
 55:       PetscViewerASCIIPushTab(viewer);
 56:       PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
 57:       PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
 58:       for (i=0; i<osm->n_local_true; i++) {
 59:         ISGetLocalSize(osm->is[i],&bsz);
 60:         PetscViewerASCIISynchronizedPrintf(sviewer,"[%d] local block number %D, size = %D\n",(int)rank,i,bsz);
 61:         KSPView(osm->ksp[i],sviewer);
 62:         PetscViewerASCIISynchronizedPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
 63:       }
 64:       PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
 65:       PetscViewerASCIIPopTab(viewer);
 66:       PetscViewerFlush(viewer);
 67:       PetscViewerASCIIPopSynchronized(viewer);
 68:     }
 69:   } else if (isstring) {
 70:     PetscViewerStringSPrintf(viewer," blocks=%D, overlap=%D, type=%s",osm->n,osm->overlap,PCASMTypes[osm->type]);
 71:     PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
 72:     if (osm->ksp) KSPView(osm->ksp[0],sviewer);
 73:     PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
 74:   }
 75:   return 0;
 76: }

 78: static PetscErrorCode PCASMPrintSubdomains(PC pc)
 79: {
 80:   PC_ASM         *osm = (PC_ASM*)pc->data;
 81:   const char     *prefix;
 82:   char           fname[PETSC_MAX_PATH_LEN+1];
 83:   PetscViewer    viewer, sviewer;
 84:   char           *s;
 85:   PetscInt       i,j,nidx;
 86:   const PetscInt *idx;
 87:   PetscMPIInt    rank, size;

 89:   MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
 90:   MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
 91:   PCGetOptionsPrefix(pc,&prefix);
 92:   PetscOptionsGetString(NULL,prefix,"-pc_asm_print_subdomains",fname,sizeof(fname),NULL);
 93:   if (fname[0] == 0) PetscStrcpy(fname,"stdout");
 94:   PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),fname,&viewer);
 95:   for (i=0; i<osm->n_local; i++) {
 96:     if (i < osm->n_local_true) {
 97:       ISGetLocalSize(osm->is[i],&nidx);
 98:       ISGetIndices(osm->is[i],&idx);
 99:       /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
100: #define len  16*(nidx+1)+512
101:       PetscMalloc1(len,&s);
102:       PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer);
103: #undef len
104:       PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D with overlap:\n", rank, size, i);
105:       for (j=0; j<nidx; j++) {
106:         PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
107:       }
108:       ISRestoreIndices(osm->is[i],&idx);
109:       PetscViewerStringSPrintf(sviewer,"\n");
110:       PetscViewerDestroy(&sviewer);
111:       PetscViewerASCIIPushSynchronized(viewer);
112:       PetscViewerASCIISynchronizedPrintf(viewer, s);
113:       PetscViewerFlush(viewer);
114:       PetscViewerASCIIPopSynchronized(viewer);
115:       PetscFree(s);
116:       if (osm->is_local) {
117:         /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
118: #define len  16*(nidx+1)+512
119:         PetscMalloc1(len, &s);
120:         PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer);
121: #undef len
122:         PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D without overlap:\n", rank, size, i);
123:         ISGetLocalSize(osm->is_local[i],&nidx);
124:         ISGetIndices(osm->is_local[i],&idx);
125:         for (j=0; j<nidx; j++) {
126:           PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
127:         }
128:         ISRestoreIndices(osm->is_local[i],&idx);
129:         PetscViewerStringSPrintf(sviewer,"\n");
130:         PetscViewerDestroy(&sviewer);
131:         PetscViewerASCIIPushSynchronized(viewer);
132:         PetscViewerASCIISynchronizedPrintf(viewer, s);
133:         PetscViewerFlush(viewer);
134:         PetscViewerASCIIPopSynchronized(viewer);
135:         PetscFree(s);
136:       }
137:     } else {
138:       /* Participate in collective viewer calls. */
139:       PetscViewerASCIIPushSynchronized(viewer);
140:       PetscViewerFlush(viewer);
141:       PetscViewerASCIIPopSynchronized(viewer);
142:       /* Assume either all ranks have is_local or none do. */
143:       if (osm->is_local) {
144:         PetscViewerASCIIPushSynchronized(viewer);
145:         PetscViewerFlush(viewer);
146:         PetscViewerASCIIPopSynchronized(viewer);
147:       }
148:     }
149:   }
150:   PetscViewerFlush(viewer);
151:   PetscViewerDestroy(&viewer);
152:   return 0;
153: }

155: static PetscErrorCode PCSetUp_ASM(PC pc)
156: {
157:   PC_ASM         *osm = (PC_ASM*)pc->data;
158:   PetscBool      flg;
159:   PetscInt       i,m,m_local;
160:   MatReuse       scall = MAT_REUSE_MATRIX;
161:   IS             isl;
162:   KSP            ksp;
163:   PC             subpc;
164:   const char     *prefix,*pprefix;
165:   Vec            vec;
166:   DM             *domain_dm = NULL;

168:   if (!pc->setupcalled) {
169:     PetscInt m;

171:     /* Note: if subdomains have been set either via PCASMSetTotalSubdomains() or via PCASMSetLocalSubdomains(), osm->n_local_true will not be PETSC_DECIDE */
172:     if (osm->n_local_true == PETSC_DECIDE) {
173:       /* no subdomains given */
174:       /* try pc->dm first, if allowed */
175:       if (osm->dm_subdomains && pc->dm) {
176:         PetscInt  num_domains, d;
177:         char      **domain_names;
178:         IS        *inner_domain_is, *outer_domain_is;
179:         DMCreateDomainDecomposition(pc->dm, &num_domains, &domain_names, &inner_domain_is, &outer_domain_is, &domain_dm);
180:         osm->overlap = -1; /* We do not want to increase the overlap of the IS.
181:                               A future improvement of this code might allow one to use
182:                               DM-defined subdomains and also increase the overlap,
183:                               but that is not currently supported */
184:         if (num_domains) {
185:           PCASMSetLocalSubdomains(pc, num_domains, outer_domain_is, inner_domain_is);
186:         }
187:         for (d = 0; d < num_domains; ++d) {
188:           if (domain_names)    PetscFree(domain_names[d]);
189:           if (inner_domain_is) ISDestroy(&inner_domain_is[d]);
190:           if (outer_domain_is) ISDestroy(&outer_domain_is[d]);
191:         }
192:         PetscFree(domain_names);
193:         PetscFree(inner_domain_is);
194:         PetscFree(outer_domain_is);
195:       }
196:       if (osm->n_local_true == PETSC_DECIDE) {
197:         /* still no subdomains; use one subdomain per processor */
198:         osm->n_local_true = 1;
199:       }
200:     }
201:     { /* determine the global and max number of subdomains */
202:       struct {PetscInt max,sum;} inwork,outwork;
203:       PetscMPIInt size;

205:       inwork.max   = osm->n_local_true;
206:       inwork.sum   = osm->n_local_true;
207:       MPIU_Allreduce(&inwork,&outwork,1,MPIU_2INT,MPIU_MAXSUM_OP,PetscObjectComm((PetscObject)pc));
208:       osm->n_local = outwork.max;
209:       osm->n       = outwork.sum;

211:       MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
212:       if (outwork.max == 1 && outwork.sum == size) {
213:         /* osm->n_local_true = 1 on all processes, set this option may enable use of optimized MatCreateSubMatrices() implementation */
214:         MatSetOption(pc->pmat,MAT_SUBMAT_SINGLEIS,PETSC_TRUE);
215:       }
216:     }
217:     if (!osm->is) { /* create the index sets */
218:       PCASMCreateSubdomains(pc->pmat,osm->n_local_true,&osm->is);
219:     }
220:     if (osm->n_local_true > 1 && !osm->is_local) {
221:       PetscMalloc1(osm->n_local_true,&osm->is_local);
222:       for (i=0; i<osm->n_local_true; i++) {
223:         if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
224:           ISDuplicate(osm->is[i],&osm->is_local[i]);
225:           ISCopy(osm->is[i],osm->is_local[i]);
226:         } else {
227:           PetscObjectReference((PetscObject)osm->is[i]);
228:           osm->is_local[i] = osm->is[i];
229:         }
230:       }
231:     }
232:     PCGetOptionsPrefix(pc,&prefix);
233:     flg  = PETSC_FALSE;
234:     PetscOptionsHasName(NULL,prefix,"-pc_asm_print_subdomains",&flg);
235:     if (flg) PCASMPrintSubdomains(pc);

237:     if (osm->overlap > 0) {
238:       /* Extend the "overlapping" regions by a number of steps */
239:       MatIncreaseOverlap(pc->pmat,osm->n_local_true,osm->is,osm->overlap);
240:     }
241:     if (osm->sort_indices) {
242:       for (i=0; i<osm->n_local_true; i++) {
243:         ISSort(osm->is[i]);
244:         if (osm->is_local) {
245:           ISSort(osm->is_local[i]);
246:         }
247:       }
248:     }

250:     if (!osm->ksp) {
251:       /* Create the local solvers */
252:       PetscMalloc1(osm->n_local_true,&osm->ksp);
253:       if (domain_dm) {
254:         PetscInfo(pc,"Setting up ASM subproblems using the embedded DM\n");
255:       }
256:       for (i=0; i<osm->n_local_true; i++) {
257:         KSPCreate(PETSC_COMM_SELF,&ksp);
258:         KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
259:         PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
260:         PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
261:         KSPSetType(ksp,KSPPREONLY);
262:         KSPGetPC(ksp,&subpc);
263:         PCGetOptionsPrefix(pc,&prefix);
264:         KSPSetOptionsPrefix(ksp,prefix);
265:         KSPAppendOptionsPrefix(ksp,"sub_");
266:         if (domain_dm) {
267:           KSPSetDM(ksp, domain_dm[i]);
268:           KSPSetDMActive(ksp, PETSC_FALSE);
269:           DMDestroy(&domain_dm[i]);
270:         }
271:         osm->ksp[i] = ksp;
272:       }
273:       if (domain_dm) {
274:         PetscFree(domain_dm);
275:       }
276:     }

278:     ISConcatenate(PETSC_COMM_SELF, osm->n_local_true, osm->is, &osm->lis);
279:     ISSortRemoveDups(osm->lis);
280:     ISGetLocalSize(osm->lis, &m);

282:     scall = MAT_INITIAL_MATRIX;
283:   } else {
284:     /*
285:        Destroy the blocks from the previous iteration
286:     */
287:     if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
288:       MatDestroyMatrices(osm->n_local_true,&osm->pmat);
289:       scall = MAT_INITIAL_MATRIX;
290:     }
291:   }

293:   /* Destroy previous submatrices of a different type than pc->pmat since MAT_REUSE_MATRIX won't work in that case */
294:   if ((scall == MAT_REUSE_MATRIX) && osm->sub_mat_type) {
295:     if (osm->n_local_true > 0) {
296:       MatDestroySubMatrices(osm->n_local_true,&osm->pmat);
297:     }
298:     scall = MAT_INITIAL_MATRIX;
299:   }

301:   /*
302:      Extract out the submatrices
303:   */
304:   MatCreateSubMatrices(pc->pmat,osm->n_local_true,osm->is,osm->is,scall,&osm->pmat);
305:   if (scall == MAT_INITIAL_MATRIX) {
306:     PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
307:     for (i=0; i<osm->n_local_true; i++) {
308:       PetscLogObjectParent((PetscObject)pc,(PetscObject)osm->pmat[i]);
309:       PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);
310:     }
311:   }

313:   /* Convert the types of the submatrices (if needbe) */
314:   if (osm->sub_mat_type) {
315:     for (i=0; i<osm->n_local_true; i++) {
316:       MatConvert(osm->pmat[i],osm->sub_mat_type,MAT_INPLACE_MATRIX,&(osm->pmat[i]));
317:     }
318:   }

320:   if (!pc->setupcalled) {
321:     VecType vtype;

323:     /* Create the local work vectors (from the local matrices) and scatter contexts */
324:     MatCreateVecs(pc->pmat,&vec,NULL);

327:     if (osm->is_local && osm->type == PC_ASM_RESTRICT && osm->loctype == PC_COMPOSITE_ADDITIVE) {
328:       PetscMalloc1(osm->n_local_true,&osm->lprolongation);
329:     }
330:     PetscMalloc1(osm->n_local_true,&osm->lrestriction);
331:     PetscMalloc1(osm->n_local_true,&osm->x);
332:     PetscMalloc1(osm->n_local_true,&osm->y);

334:     ISGetLocalSize(osm->lis,&m);
335:     ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
336:     MatGetVecType(osm->pmat[0],&vtype);
337:     VecCreate(PETSC_COMM_SELF,&osm->lx);
338:     VecSetSizes(osm->lx,m,m);
339:     VecSetType(osm->lx,vtype);
340:     VecDuplicate(osm->lx, &osm->ly);
341:     VecScatterCreate(vec,osm->lis,osm->lx,isl,&osm->restriction);
342:     ISDestroy(&isl);

344:     for (i=0; i<osm->n_local_true; ++i) {
345:       ISLocalToGlobalMapping ltog;
346:       IS                     isll;
347:       const PetscInt         *idx_is;
348:       PetscInt               *idx_lis,nout;

350:       ISGetLocalSize(osm->is[i],&m);
351:       MatCreateVecs(osm->pmat[i],&osm->x[i],NULL);
352:       VecDuplicate(osm->x[i],&osm->y[i]);

354:       /* generate a scatter from ly to y[i] picking all the overlapping is[i] entries */
355:       ISLocalToGlobalMappingCreateIS(osm->lis,&ltog);
356:       ISGetLocalSize(osm->is[i],&m);
357:       ISGetIndices(osm->is[i], &idx_is);
358:       PetscMalloc1(m,&idx_lis);
359:       ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m,idx_is,&nout,idx_lis);
361:       ISRestoreIndices(osm->is[i], &idx_is);
362:       ISCreateGeneral(PETSC_COMM_SELF,m,idx_lis,PETSC_OWN_POINTER,&isll);
363:       ISLocalToGlobalMappingDestroy(&ltog);
364:       ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
365:       VecScatterCreate(osm->ly,isll,osm->y[i],isl,&osm->lrestriction[i]);
366:       ISDestroy(&isll);
367:       ISDestroy(&isl);
368:       if (osm->lprolongation) { /* generate a scatter from y[i] to ly picking only the the non-overlapping is_local[i] entries */
369:         ISLocalToGlobalMapping ltog;
370:         IS                     isll,isll_local;
371:         const PetscInt         *idx_local;
372:         PetscInt               *idx1, *idx2, nout;

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

377:         ISLocalToGlobalMappingCreateIS(osm->is[i],&ltog);
378:         PetscMalloc1(m_local,&idx1);
379:         ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx1);
380:         ISLocalToGlobalMappingDestroy(&ltog);
382:         ISCreateGeneral(PETSC_COMM_SELF,m_local,idx1,PETSC_OWN_POINTER,&isll);

384:         ISLocalToGlobalMappingCreateIS(osm->lis,&ltog);
385:         PetscMalloc1(m_local,&idx2);
386:         ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx2);
387:         ISLocalToGlobalMappingDestroy(&ltog);
389:         ISCreateGeneral(PETSC_COMM_SELF,m_local,idx2,PETSC_OWN_POINTER,&isll_local);

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

394:         ISDestroy(&isll);
395:         ISDestroy(&isll_local);
396:       }
397:     }
398:     VecDestroy(&vec);
399:   }

401:   if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
402:     IS      *cis;
403:     PetscInt c;

405:     PetscMalloc1(osm->n_local_true, &cis);
406:     for (c = 0; c < osm->n_local_true; ++c) cis[c] = osm->lis;
407:     MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, cis, scall, &osm->lmats);
408:     PetscFree(cis);
409:   }

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

415:   /*
416:      Loop over subdomains putting them into local ksp
417:   */
418:   KSPGetOptionsPrefix(osm->ksp[0],&prefix);
419:   for (i=0; i<osm->n_local_true; i++) {
420:     KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i]);
421:     MatSetOptionsPrefix(osm->pmat[i],prefix);
422:     if (!pc->setupcalled) {
423:       KSPSetFromOptions(osm->ksp[i]);
424:     }
425:   }
426:   return 0;
427: }

429: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
430: {
431:   PC_ASM             *osm = (PC_ASM*)pc->data;
432:   PetscInt           i;
433:   KSPConvergedReason reason;

435:   for (i=0; i<osm->n_local_true; i++) {
436:     KSPSetUp(osm->ksp[i]);
437:     KSPGetConvergedReason(osm->ksp[i],&reason);
438:     if (reason == KSP_DIVERGED_PC_FAILED) {
439:       pc->failedreason = PC_SUBPC_ERROR;
440:     }
441:   }
442:   return 0;
443: }

445: static PetscErrorCode PCApply_ASM(PC pc,Vec x,Vec y)
446: {
447:   PC_ASM         *osm = (PC_ASM*)pc->data;
448:   PetscInt       i,n_local_true = osm->n_local_true;
449:   ScatterMode    forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;

451:   /*
452:      support for limiting the restriction or interpolation to only local
453:      subdomain values (leaving the other values 0).
454:   */
455:   if (!(osm->type & PC_ASM_RESTRICT)) {
456:     forward = SCATTER_FORWARD_LOCAL;
457:     /* have to zero the work RHS since scatter may leave some slots empty */
458:     VecSet(osm->lx, 0.0);
459:   }
460:   if (!(osm->type & PC_ASM_INTERPOLATE)) {
461:     reverse = SCATTER_REVERSE_LOCAL;
462:   }

464:   if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE) {
465:     /* zero the global and the local solutions */
466:     VecSet(y, 0.0);
467:     VecSet(osm->ly, 0.0);

469:     /* copy the global RHS to local RHS including the ghost nodes */
470:     VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
471:     VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);

473:     /* restrict local RHS to the overlapping 0-block RHS */
474:     VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
475:     VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);

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

480:       /* solve the overlapping i-block */
481:       PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i],0);
482:       KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]);
483:       KSPCheckSolve(osm->ksp[i], pc, osm->y[i]);
484:       PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0);

486:       if (osm->lprolongation) { /* interpolate the non-overlapping i-block solution to the local solution (only for restrictive additive) */
487:         VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
488:         VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
489:       } else { /* interpolate the overlapping i-block solution to the local solution */
490:         VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
491:         VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
492:       }

494:       if (i < n_local_true-1) {
495:         /* restrict local RHS to the overlapping (i+1)-block RHS */
496:         VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
497:         VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);

499:         if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
500:           /* update the overlapping (i+1)-block RHS using the current local solution */
501:           MatMult(osm->lmats[i+1], osm->ly, osm->y[i+1]);
502:           VecAXPBY(osm->x[i+1],-1.,1., osm->y[i+1]);
503:         }
504:       }
505:     }
506:     /* add the local solution to the global solution including the ghost nodes */
507:     VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
508:     VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
509:   } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
510:   return 0;
511: }

513: static PetscErrorCode PCMatApply_ASM(PC pc,Mat X,Mat Y)
514: {
515:   PC_ASM         *osm = (PC_ASM*)pc->data;
516:   Mat            Z,W;
517:   Vec            x;
518:   PetscInt       i,m,N;
519:   ScatterMode    forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;

522:   /*
523:      support for limiting the restriction or interpolation to only local
524:      subdomain values (leaving the other values 0).
525:   */
526:   if (!(osm->type & PC_ASM_RESTRICT)) {
527:     forward = SCATTER_FORWARD_LOCAL;
528:     /* have to zero the work RHS since scatter may leave some slots empty */
529:     VecSet(osm->lx, 0.0);
530:   }
531:   if (!(osm->type & PC_ASM_INTERPOLATE)) {
532:     reverse = SCATTER_REVERSE_LOCAL;
533:   }
534:   VecGetLocalSize(osm->x[0], &m);
535:   MatGetSize(X, NULL, &N);
536:   MatCreateSeqDense(PETSC_COMM_SELF, m, N, NULL, &Z);
537:   if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE) {
538:     /* zero the global and the local solutions */
539:     MatZeroEntries(Y);
540:     VecSet(osm->ly, 0.0);

542:     for (i = 0; i < N; ++i) {
543:       MatDenseGetColumnVecRead(X, i, &x);
544:       /* copy the global RHS to local RHS including the ghost nodes */
545:       VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
546:       VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
547:       MatDenseRestoreColumnVecRead(X, i, &x);

549:       MatDenseGetColumnVecWrite(Z, i, &x);
550:       /* restrict local RHS to the overlapping 0-block RHS */
551:       VecScatterBegin(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward);
552:       VecScatterEnd(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward);
553:       MatDenseRestoreColumnVecWrite(Z, i, &x);
554:     }
555:     MatCreateSeqDense(PETSC_COMM_SELF, m, N, NULL, &W);
556:     /* solve the overlapping 0-block */
557:     PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[0], Z, W, 0);
558:     KSPMatSolve(osm->ksp[0], Z, W);
559:     KSPCheckSolve(osm->ksp[0], pc, NULL);
560:     PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[0], Z, W,0);
561:     MatDestroy(&Z);

563:     for (i = 0; i < N; ++i) {
564:       VecSet(osm->ly, 0.0);
565:       MatDenseGetColumnVecRead(W, i, &x);
566:       if (osm->lprolongation) { /* interpolate the non-overlapping 0-block solution to the local solution (only for restrictive additive) */
567:         VecScatterBegin(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward);
568:         VecScatterEnd(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward);
569:       } else { /* interpolate the overlapping 0-block solution to the local solution */
570:         VecScatterBegin(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse);
571:         VecScatterEnd(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse);
572:       }
573:       MatDenseRestoreColumnVecRead(W, i, &x);

575:       MatDenseGetColumnVecWrite(Y, i, &x);
576:       /* add the local solution to the global solution including the ghost nodes */
577:       VecScatterBegin(osm->restriction, osm->ly, x, ADD_VALUES, reverse);
578:       VecScatterEnd(osm->restriction, osm->ly, x, ADD_VALUES, reverse);
579:       MatDenseRestoreColumnVecWrite(Y, i, &x);
580:     }
581:     MatDestroy(&W);
582:   } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
583:   return 0;
584: }

586: static PetscErrorCode PCApplyTranspose_ASM(PC pc,Vec x,Vec y)
587: {
588:   PC_ASM         *osm = (PC_ASM*)pc->data;
589:   PetscInt       i,n_local_true = osm->n_local_true;
590:   ScatterMode    forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;

592:   /*
593:      Support for limiting the restriction or interpolation to only local
594:      subdomain values (leaving the other values 0).

596:      Note: these are reversed from the PCApply_ASM() because we are applying the
597:      transpose of the three terms
598:   */

600:   if (!(osm->type & PC_ASM_INTERPOLATE)) {
601:     forward = SCATTER_FORWARD_LOCAL;
602:     /* have to zero the work RHS since scatter may leave some slots empty */
603:     VecSet(osm->lx, 0.0);
604:   }
605:   if (!(osm->type & PC_ASM_RESTRICT)) reverse = SCATTER_REVERSE_LOCAL;

607:   /* zero the global and the local solutions */
608:   VecSet(y, 0.0);
609:   VecSet(osm->ly, 0.0);

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

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

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

622:     /* solve the overlapping i-block */
623:     PetscLogEventBegin(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
624:     KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i]);
625:     KSPCheckSolve(osm->ksp[i],pc,osm->y[i]);
626:     PetscLogEventEnd(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);

628:     if (osm->lprolongation) { /* interpolate the non-overlapping i-block solution to the local solution */
629:       VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
630:       VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
631:     } else { /* interpolate the overlapping i-block solution to the local solution */
632:       VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
633:       VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
634:     }

636:     if (i < n_local_true-1) {
637:       /* Restrict local RHS to the overlapping (i+1)-block RHS */
638:       VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
639:       VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
640:     }
641:   }
642:   /* Add the local solution to the global solution including the ghost nodes */
643:   VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
644:   VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);

646:   return 0;

648: }

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

655:   if (osm->ksp) {
656:     for (i=0; i<osm->n_local_true; i++) {
657:       KSPReset(osm->ksp[i]);
658:     }
659:   }
660:   if (osm->pmat) {
661:     if (osm->n_local_true > 0) {
662:       MatDestroySubMatrices(osm->n_local_true,&osm->pmat);
663:     }
664:   }
665:   if (osm->lrestriction) {
666:     VecScatterDestroy(&osm->restriction);
667:     for (i=0; i<osm->n_local_true; i++) {
668:       VecScatterDestroy(&osm->lrestriction[i]);
669:       if (osm->lprolongation) VecScatterDestroy(&osm->lprolongation[i]);
670:       VecDestroy(&osm->x[i]);
671:       VecDestroy(&osm->y[i]);
672:     }
673:     PetscFree(osm->lrestriction);
674:     if (osm->lprolongation) PetscFree(osm->lprolongation);
675:     PetscFree(osm->x);
676:     PetscFree(osm->y);

678:   }
679:   PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
680:   ISDestroy(&osm->lis);
681:   VecDestroy(&osm->lx);
682:   VecDestroy(&osm->ly);
683:   if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
684:     MatDestroyMatrices(osm->n_local_true, &osm->lmats);
685:   }

687:   PetscFree(osm->sub_mat_type);

689:   osm->is       = NULL;
690:   osm->is_local = NULL;
691:   return 0;
692: }

694: static PetscErrorCode PCDestroy_ASM(PC pc)
695: {
696:   PC_ASM         *osm = (PC_ASM*)pc->data;
697:   PetscInt       i;

699:   PCReset_ASM(pc);
700:   if (osm->ksp) {
701:     for (i=0; i<osm->n_local_true; i++) {
702:       KSPDestroy(&osm->ksp[i]);
703:     }
704:     PetscFree(osm->ksp);
705:   }
706:   PetscFree(pc->data);

708:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",NULL);
709:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",NULL);
710:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetOverlap_C",NULL);
711:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetType_C",NULL);
712:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetType_C",NULL);
713:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalType_C",NULL);
714:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetLocalType_C",NULL);
715:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSortIndices_C",NULL);
716:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubKSP_C",NULL);
717:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubMatType_C",NULL);
718:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSubMatType_C",NULL);
719:   return 0;
720: }

722: static PetscErrorCode PCSetFromOptions_ASM(PetscOptionItems *PetscOptionsObject,PC pc)
723: {
724:   PC_ASM         *osm = (PC_ASM*)pc->data;
725:   PetscInt       blocks,ovl;
726:   PetscBool      flg;
727:   PCASMType      asmtype;
728:   PCCompositeType loctype;
729:   char           sub_mat_type[256];

731:   PetscOptionsHead(PetscOptionsObject,"Additive Schwarz options");
732:   PetscOptionsBool("-pc_asm_dm_subdomains","Use DMCreateDomainDecomposition() to define subdomains","PCASMSetDMSubdomains",osm->dm_subdomains,&osm->dm_subdomains,&flg);
733:   PetscOptionsInt("-pc_asm_blocks","Number of subdomains","PCASMSetTotalSubdomains",osm->n,&blocks,&flg);
734:   if (flg) {
735:     PCASMSetTotalSubdomains(pc,blocks,NULL,NULL);
736:     osm->dm_subdomains = PETSC_FALSE;
737:   }
738:   PetscOptionsInt("-pc_asm_local_blocks","Number of local subdomains","PCASMSetLocalSubdomains",osm->n_local_true,&blocks,&flg);
739:   if (flg) {
740:     PCASMSetLocalSubdomains(pc,blocks,NULL,NULL);
741:     osm->dm_subdomains = PETSC_FALSE;
742:   }
743:   PetscOptionsInt("-pc_asm_overlap","Number of grid points overlap","PCASMSetOverlap",osm->overlap,&ovl,&flg);
744:   if (flg) {
745:     PCASMSetOverlap(pc,ovl);
746:     osm->dm_subdomains = PETSC_FALSE;
747:   }
748:   flg  = PETSC_FALSE;
749:   PetscOptionsEnum("-pc_asm_type","Type of restriction/extension","PCASMSetType",PCASMTypes,(PetscEnum)osm->type,(PetscEnum*)&asmtype,&flg);
750:   if (flg) PCASMSetType(pc,asmtype);
751:   flg  = PETSC_FALSE;
752:   PetscOptionsEnum("-pc_asm_local_type","Type of local solver composition","PCASMSetLocalType",PCCompositeTypes,(PetscEnum)osm->loctype,(PetscEnum*)&loctype,&flg);
753:   if (flg) PCASMSetLocalType(pc,loctype);
754:   PetscOptionsFList("-pc_asm_sub_mat_type","Subsolve Matrix Type","PCASMSetSubMatType",MatList,NULL,sub_mat_type,256,&flg);
755:   if (flg) {
756:     PCASMSetSubMatType(pc,sub_mat_type);
757:   }
758:   PetscOptionsTail();
759:   return 0;
760: }

762: /*------------------------------------------------------------------------------------*/

764: static PetscErrorCode  PCASMSetLocalSubdomains_ASM(PC pc,PetscInt n,IS is[],IS is_local[])
765: {
766:   PC_ASM         *osm = (PC_ASM*)pc->data;
767:   PetscInt       i;


772:   if (!pc->setupcalled) {
773:     if (is) {
774:       for (i=0; i<n; i++) PetscObjectReference((PetscObject)is[i]);
775:     }
776:     if (is_local) {
777:       for (i=0; i<n; i++) PetscObjectReference((PetscObject)is_local[i]);
778:     }
779:     PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);

781:     osm->n_local_true = n;
782:     osm->is           = NULL;
783:     osm->is_local     = NULL;
784:     if (is) {
785:       PetscMalloc1(n,&osm->is);
786:       for (i=0; i<n; i++) osm->is[i] = is[i];
787:       /* Flag indicating that the user has set overlapping subdomains so PCASM should not increase their size. */
788:       osm->overlap = -1;
789:     }
790:     if (is_local) {
791:       PetscMalloc1(n,&osm->is_local);
792:       for (i=0; i<n; i++) osm->is_local[i] = is_local[i];
793:       if (!is) {
794:         PetscMalloc1(osm->n_local_true,&osm->is);
795:         for (i=0; i<osm->n_local_true; i++) {
796:           if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
797:             ISDuplicate(osm->is_local[i],&osm->is[i]);
798:             ISCopy(osm->is_local[i],osm->is[i]);
799:           } else {
800:             PetscObjectReference((PetscObject)osm->is_local[i]);
801:             osm->is[i] = osm->is_local[i];
802:           }
803:         }
804:       }
805:     }
806:   }
807:   return 0;
808: }

810: static PetscErrorCode  PCASMSetTotalSubdomains_ASM(PC pc,PetscInt N,IS *is,IS *is_local)
811: {
812:   PC_ASM         *osm = (PC_ASM*)pc->data;
813:   PetscMPIInt    rank,size;
814:   PetscInt       n;


819:   /*
820:      Split the subdomains equally among all processors
821:   */
822:   MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
823:   MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
824:   n    = N/size + ((N % size) > rank);
827:   if (!pc->setupcalled) {
828:     PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);

830:     osm->n_local_true = n;
831:     osm->is           = NULL;
832:     osm->is_local     = NULL;
833:   }
834:   return 0;
835: }

837: static PetscErrorCode  PCASMSetOverlap_ASM(PC pc,PetscInt ovl)
838: {
839:   PC_ASM *osm = (PC_ASM*)pc->data;

843:   if (!pc->setupcalled) osm->overlap = ovl;
844:   return 0;
845: }

847: static PetscErrorCode  PCASMSetType_ASM(PC pc,PCASMType type)
848: {
849:   PC_ASM *osm = (PC_ASM*)pc->data;

851:   osm->type     = type;
852:   osm->type_set = PETSC_TRUE;
853:   return 0;
854: }

856: static PetscErrorCode  PCASMGetType_ASM(PC pc,PCASMType *type)
857: {
858:   PC_ASM *osm = (PC_ASM*)pc->data;

860:   *type = osm->type;
861:   return 0;
862: }

864: static PetscErrorCode  PCASMSetLocalType_ASM(PC pc, PCCompositeType type)
865: {
866:   PC_ASM *osm = (PC_ASM *) pc->data;

869:   osm->loctype = type;
870:   return 0;
871: }

873: static PetscErrorCode  PCASMGetLocalType_ASM(PC pc, PCCompositeType *type)
874: {
875:   PC_ASM *osm = (PC_ASM *) pc->data;

877:   *type = osm->loctype;
878:   return 0;
879: }

881: static PetscErrorCode  PCASMSetSortIndices_ASM(PC pc,PetscBool  doSort)
882: {
883:   PC_ASM *osm = (PC_ASM*)pc->data;

885:   osm->sort_indices = doSort;
886:   return 0;
887: }

889: static PetscErrorCode  PCASMGetSubKSP_ASM(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
890: {
891:   PC_ASM         *osm = (PC_ASM*)pc->data;


895:   if (n_local) *n_local = osm->n_local_true;
896:   if (first_local) {
897:     MPI_Scan(&osm->n_local_true,first_local,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
898:     *first_local -= osm->n_local_true;
899:   }
900:   if (ksp) *ksp   = osm->ksp;
901:   return 0;
902: }

904: static PetscErrorCode  PCASMGetSubMatType_ASM(PC pc,MatType *sub_mat_type)
905: {
906:   PC_ASM         *osm = (PC_ASM*)pc->data;

910:   *sub_mat_type = osm->sub_mat_type;
911:   return 0;
912: }

914: static PetscErrorCode PCASMSetSubMatType_ASM(PC pc,MatType sub_mat_type)
915: {
916:   PC_ASM            *osm = (PC_ASM*)pc->data;

919:   PetscFree(osm->sub_mat_type);
920:   PetscStrallocpy(sub_mat_type,(char**)&osm->sub_mat_type);
921:   return 0;
922: }

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

927:     Collective on pc

929:     Input Parameters:
930: +   pc - the preconditioner context
931: .   n - the number of subdomains for this processor (default value = 1)
932: .   is - the index set that defines the subdomains for this processor
933:          (or NULL for PETSc to determine subdomains)
934: -   is_local - the index sets that define the local part of the subdomains for this processor, not used unless PCASMType is PC_ASM_RESTRICT
935:          (or NULL to not provide these)

937:     Options Database Key:
938:     To set the total number of subdomain blocks rather than specify the
939:     index sets, use the option
940: .    -pc_asm_local_blocks <blks> - Sets local blocks

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

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

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

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

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

955:     Level: advanced

957: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
958:           PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), PCASMType, PCASMSetType()
959: @*/
960: PetscErrorCode  PCASMSetLocalSubdomains(PC pc,PetscInt n,IS is[],IS is_local[])
961: {
963:   PetscTryMethod(pc,"PCASMSetLocalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,is,is_local));
964:   return 0;
965: }

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

972:     Collective on pc

974:     Input Parameters:
975: +   pc - the preconditioner context
976: .   N  - the number of subdomains for all processors
977: .   is - the index sets that define the subdomains for all processors
978:          (or NULL to ask PETSc to determine the subdomains)
979: -   is_local - the index sets that define the local part of the subdomains for this processor
980:          (or NULL to not provide this information)

982:     Options Database Key:
983:     To set the total number of subdomain blocks rather than specify the
984:     index sets, use the option
985: .    -pc_asm_blocks <blks> - Sets total blocks

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

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

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

995:     Use PCASMSetLocalSubdomains() to set local subdomains.

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

999:     Level: advanced

1001: .seealso: PCASMSetLocalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1002:           PCASMCreateSubdomains2D()
1003: @*/
1004: PetscErrorCode  PCASMSetTotalSubdomains(PC pc,PetscInt N,IS is[],IS is_local[])
1005: {
1007:   PetscTryMethod(pc,"PCASMSetTotalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,N,is,is_local));
1008:   return 0;
1009: }

1011: /*@
1012:     PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
1013:     additive Schwarz preconditioner.  Either all or no processors in the
1014:     PC communicator must call this routine.

1016:     Logically Collective on pc

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

1022:     Options Database Key:
1023: .   -pc_asm_overlap <ovl> - Sets overlap

1025:     Notes:
1026:     By default the ASM preconditioner uses 1 block per processor.  To use
1027:     multiple blocks per perocessor, see PCASMSetTotalSubdomains() and
1028:     PCASMSetLocalSubdomains() (and the option -pc_asm_blocks <blks>).

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

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

1041:     Note that one can define initial index sets with any overlap via
1042:     PCASMSetLocalSubdomains(); the routine
1043:     PCASMSetOverlap() merely allows PETSc to extend that overlap further
1044:     if desired.

1046:     Level: intermediate

1048: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1049:           PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), MatIncreaseOverlap()
1050: @*/
1051: PetscErrorCode  PCASMSetOverlap(PC pc,PetscInt ovl)
1052: {
1055:   PetscTryMethod(pc,"PCASMSetOverlap_C",(PC,PetscInt),(pc,ovl));
1056:   return 0;
1057: }

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:     Notes:
1079:     if the is_local arguments are passed to PCASMSetLocalSubdomains() then they are used when PC_ASM_RESTRICT has been selected
1080:     to limit the local processor interpolation

1082:     Level: intermediate

1084: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1085:           PCASMCreateSubdomains2D(), PCASMType, PCASMSetLocalType(), PCASMGetLocalType()
1086: @*/
1087: PetscErrorCode  PCASMSetType(PC pc,PCASMType type)
1088: {
1091:   PetscTryMethod(pc,"PCASMSetType_C",(PC,PCASMType),(pc,type));
1092:   return 0;
1093: }

1095: /*@
1096:     PCASMGetType - Gets the type of restriction and interpolation used
1097:     for local problems in the additive Schwarz method.

1099:     Logically Collective on pc

1101:     Input Parameter:
1102: .   pc  - the preconditioner context

1104:     Output Parameter:
1105: .   type - variant of ASM, one of

1107: .vb
1108:       PC_ASM_BASIC       - full interpolation and restriction
1109:       PC_ASM_RESTRICT    - full restriction, local processor interpolation
1110:       PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1111:       PC_ASM_NONE        - local processor restriction and interpolation
1112: .ve

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

1117:     Level: intermediate

1119: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1120:           PCASMCreateSubdomains2D(), PCASMType, PCASMSetType(), PCASMSetLocalType(), PCASMGetLocalType()
1121: @*/
1122: PetscErrorCode  PCASMGetType(PC pc,PCASMType *type)
1123: {
1125:   PetscUseMethod(pc,"PCASMGetType_C",(PC,PCASMType*),(pc,type));
1126:   return 0;
1127: }

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

1132:   Logically Collective on pc

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

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

1145:   Level: intermediate

1147: .seealso: PCASMSetType(), PCASMGetType(), PCASMGetLocalType(), PCASM, PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1148: @*/
1149: PetscErrorCode PCASMSetLocalType(PC pc, PCCompositeType type)
1150: {
1153:   PetscTryMethod(pc, "PCASMSetLocalType_C", (PC, PCCompositeType), (pc, type));
1154:   return 0;
1155: }

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

1160:   Logically Collective on pc

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

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

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

1175:   Level: intermediate

1177: .seealso: PCASMSetType(), PCASMGetType(), PCASMSetLocalType(), PCASMCreate(), PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1178: @*/
1179: PetscErrorCode PCASMGetLocalType(PC pc, PCCompositeType *type)
1180: {
1183:   PetscUseMethod(pc, "PCASMGetLocalType_C", (PC, PCCompositeType *), (pc, type));
1184:   return 0;
1185: }

1187: /*@
1188:     PCASMSetSortIndices - Determines whether subdomain indices are sorted.

1190:     Logically Collective on pc

1192:     Input Parameters:
1193: +   pc  - the preconditioner context
1194: -   doSort - sort the subdomain indices

1196:     Level: intermediate

1198: .seealso: PCASMSetLocalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1199:           PCASMCreateSubdomains2D()
1200: @*/
1201: PetscErrorCode  PCASMSetSortIndices(PC pc,PetscBool doSort)
1202: {
1205:   PetscTryMethod(pc,"PCASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));
1206:   return 0;
1207: }

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

1213:    Collective on pc iff first_local is requested

1215:    Input Parameter:
1216: .  pc - the preconditioner context

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

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

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

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

1232:    Level: advanced

1234: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap(),
1235:           PCASMCreateSubdomains2D(),
1236: @*/
1237: PetscErrorCode  PCASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
1238: {
1240:   PetscUseMethod(pc,"PCASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
1241:   return 0;
1242: }

1244: /* -------------------------------------------------------------------------------------*/
1245: /*MC
1246:    PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1247:            its own KSP object.

1249:    Options Database Keys:
1250: +  -pc_asm_blocks <blks> - Sets total blocks
1251: .  -pc_asm_overlap <ovl> - Sets overlap
1252: .  -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type, default is restrict
1253: -  -pc_asm_local_type [additive, multiplicative] - Sets ASM type, default is additive

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

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

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

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

1270:    Level: beginner

1272:     References:
1273: +   * - M Dryja, OB Widlund, An additive variant of the Schwarz alternating method for the case of many subregions
1274:      Courant Institute, New York University Technical report
1275: -   * - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations,
1276:     Cambridge University Press.

1278: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
1279:            PCBJACOBI, PCASMGetSubKSP(), PCASMSetLocalSubdomains(), PCASMType, PCASMGetType(), PCASMSetLocalType(), PCASMGetLocalType()
1280:            PCASMSetTotalSubdomains(), PCSetModifySubMatrices(), PCASMSetOverlap(), PCASMSetType(), PCCompositeType

1282: M*/

1284: PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1285: {
1286:   PC_ASM         *osm;

1288:   PetscNewLog(pc,&osm);

1290:   osm->n                 = PETSC_DECIDE;
1291:   osm->n_local           = 0;
1292:   osm->n_local_true      = PETSC_DECIDE;
1293:   osm->overlap           = 1;
1294:   osm->ksp               = NULL;
1295:   osm->restriction       = NULL;
1296:   osm->lprolongation     = NULL;
1297:   osm->lrestriction      = NULL;
1298:   osm->x                 = NULL;
1299:   osm->y                 = NULL;
1300:   osm->is                = NULL;
1301:   osm->is_local          = NULL;
1302:   osm->mat               = NULL;
1303:   osm->pmat              = NULL;
1304:   osm->type              = PC_ASM_RESTRICT;
1305:   osm->loctype           = PC_COMPOSITE_ADDITIVE;
1306:   osm->sort_indices      = PETSC_TRUE;
1307:   osm->dm_subdomains     = PETSC_FALSE;
1308:   osm->sub_mat_type      = NULL;

1310:   pc->data                 = (void*)osm;
1311:   pc->ops->apply           = PCApply_ASM;
1312:   pc->ops->matapply        = PCMatApply_ASM;
1313:   pc->ops->applytranspose  = PCApplyTranspose_ASM;
1314:   pc->ops->setup           = PCSetUp_ASM;
1315:   pc->ops->reset           = PCReset_ASM;
1316:   pc->ops->destroy         = PCDestroy_ASM;
1317:   pc->ops->setfromoptions  = PCSetFromOptions_ASM;
1318:   pc->ops->setuponblocks   = PCSetUpOnBlocks_ASM;
1319:   pc->ops->view            = PCView_ASM;
1320:   pc->ops->applyrichardson = NULL;

1322:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",PCASMSetLocalSubdomains_ASM);
1323:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",PCASMSetTotalSubdomains_ASM);
1324:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetOverlap_C",PCASMSetOverlap_ASM);
1325:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetType_C",PCASMSetType_ASM);
1326:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetType_C",PCASMGetType_ASM);
1327:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalType_C",PCASMSetLocalType_ASM);
1328:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetLocalType_C",PCASMGetLocalType_ASM);
1329:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSortIndices_C",PCASMSetSortIndices_ASM);
1330:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubKSP_C",PCASMGetSubKSP_ASM);
1331:   PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubMatType_C",PCASMGetSubMatType_ASM);
1332:   PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSubMatType_C",PCASMSetSubMatType_ASM);
1333:   return 0;
1334: }

1336: /*@C
1337:    PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1338:    preconditioner for any problem on a general grid.

1340:    Collective

1342:    Input Parameters:
1343: +  A - The global matrix operator
1344: -  n - the number of local blocks

1346:    Output Parameters:
1347: .  outis - the array of index sets defining the subdomains

1349:    Level: advanced

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

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

1356: .seealso: PCASMSetLocalSubdomains(), PCASMDestroySubdomains()
1357: @*/
1358: PetscErrorCode  PCASMCreateSubdomains(Mat A, PetscInt n, IS* outis[])
1359: {
1360:   MatPartitioning mpart;
1361:   const char      *prefix;
1362:   PetscInt        i,j,rstart,rend,bs;
1363:   PetscBool       hasop, isbaij = PETSC_FALSE,foundpart = PETSC_FALSE;
1364:   Mat             Ad     = NULL, adj;
1365:   IS              ispart,isnumb,*is;


1371:   /* Get prefix, row distribution, and block size */
1372:   MatGetOptionsPrefix(A,&prefix);
1373:   MatGetOwnershipRange(A,&rstart,&rend);
1374:   MatGetBlockSize(A,&bs);

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

1448:   PetscMalloc1(n,&is);
1449:   *outis = is;

1451:   if (!foundpart) {

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

1455:     PetscInt mbs   = (rend-rstart)/bs;
1456:     PetscInt start = rstart;
1457:     for (i=0; i<n; i++) {
1458:       PetscInt count = (mbs/n + ((mbs % n) > i)) * bs;
1459:       ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);
1460:       start += count;
1461:     }

1463:   } else {

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

1467:     const PetscInt *numbering;
1468:     PetscInt       *count,nidx,*indices,*newidx,start=0;
1469:     /* Get node count in each partition */
1470:     PetscMalloc1(n,&count);
1471:     ISPartitioningCount(ispart,n,count);
1472:     if (isbaij && bs > 1) { /* adjust for the block-aij case */
1473:       for (i=0; i<n; i++) count[i] *= bs;
1474:     }
1475:     /* Build indices from node numbering */
1476:     ISGetLocalSize(isnumb,&nidx);
1477:     PetscMalloc1(nidx,&indices);
1478:     for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */
1479:     ISGetIndices(isnumb,&numbering);
1480:     PetscSortIntWithPermutation(nidx,numbering,indices);
1481:     ISRestoreIndices(isnumb,&numbering);
1482:     if (isbaij && bs > 1) { /* adjust for the block-aij case */
1483:       PetscMalloc1(nidx*bs,&newidx);
1484:       for (i=0; i<nidx; i++) {
1485:         for (j=0; j<bs; j++) newidx[i*bs+j] = indices[i]*bs + j;
1486:       }
1487:       PetscFree(indices);
1488:       nidx   *= bs;
1489:       indices = newidx;
1490:     }
1491:     /* Shift to get global indices */
1492:     for (i=0; i<nidx; i++) indices[i] += rstart;

1494:     /* Build the index sets for each block */
1495:     for (i=0; i<n; i++) {
1496:       ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);
1497:       ISSort(is[i]);
1498:       start += count[i];
1499:     }

1501:     PetscFree(count);
1502:     PetscFree(indices);
1503:     ISDestroy(&isnumb);
1504:     ISDestroy(&ispart);

1506:   }
1507:   return 0;
1508: }

1510: /*@C
1511:    PCASMDestroySubdomains - Destroys the index sets created with
1512:    PCASMCreateSubdomains(). Should be called after setting subdomains
1513:    with PCASMSetLocalSubdomains().

1515:    Collective

1517:    Input Parameters:
1518: +  n - the number of index sets
1519: .  is - the array of index sets
1520: -  is_local - the array of local index sets, can be NULL

1522:    Level: advanced

1524: .seealso: PCASMCreateSubdomains(), PCASMSetLocalSubdomains()
1525: @*/
1526: PetscErrorCode  PCASMDestroySubdomains(PetscInt n, IS is[], IS is_local[])
1527: {
1528:   PetscInt       i;

1530:   if (n <= 0) return 0;
1531:   if (is) {
1533:     for (i=0; i<n; i++) ISDestroy(&is[i]);
1534:     PetscFree(is);
1535:   }
1536:   if (is_local) {
1538:     for (i=0; i<n; i++) ISDestroy(&is_local[i]);
1539:     PetscFree(is_local);
1540:   }
1541:   return 0;
1542: }

1544: /*@
1545:    PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1546:    preconditioner for a two-dimensional problem on a regular grid.

1548:    Not Collective

1550:    Input Parameters:
1551: +  m   - the number of mesh points in the x direction
1552: .  n   - the number of mesh points in the y direction
1553: .  M   - the number of subdomains in the x direction
1554: .  N   - the number of subdomains in the y direction
1555: .  dof - degrees of freedom per node
1556: -  overlap - overlap in mesh lines

1558:    Output Parameters:
1559: +  Nsub - the number of subdomains created
1560: .  is - array of index sets defining overlapping (if overlap > 0) subdomains
1561: -  is_local - array of index sets defining non-overlapping subdomains

1563:    Note:
1564:    Presently PCAMSCreateSubdomains2d() is valid only for sequential
1565:    preconditioners.  More general related routines are
1566:    PCASMSetTotalSubdomains() and PCASMSetLocalSubdomains().

1568:    Level: advanced

1570: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1571:           PCASMSetOverlap()
1572: @*/
1573: PetscErrorCode  PCASMCreateSubdomains2D(PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt dof,PetscInt overlap,PetscInt *Nsub,IS **is,IS **is_local)
1574: {
1575:   PetscInt       i,j,height,width,ystart,xstart,yleft,yright,xleft,xright,loc_outer;
1576:   PetscInt       nidx,*idx,loc,ii,jj,count;


1580:   *Nsub     = N*M;
1581:   PetscMalloc1(*Nsub,is);
1582:   PetscMalloc1(*Nsub,is_local);
1583:   ystart    = 0;
1584:   loc_outer = 0;
1585:   for (i=0; i<N; i++) {
1586:     height = n/N + ((n % N) > i); /* height of subdomain */
1588:     yleft  = ystart - overlap; if (yleft < 0) yleft = 0;
1589:     yright = ystart + height + overlap; if (yright > n) yright = n;
1590:     xstart = 0;
1591:     for (j=0; j<M; j++) {
1592:       width = m/M + ((m % M) > j); /* width of subdomain */
1594:       xleft  = xstart - overlap; if (xleft < 0) xleft = 0;
1595:       xright = xstart + width + overlap; if (xright > m) xright = m;
1596:       nidx   = (xright - xleft)*(yright - yleft);
1597:       PetscMalloc1(nidx,&idx);
1598:       loc    = 0;
1599:       for (ii=yleft; ii<yright; ii++) {
1600:         count = m*ii + xleft;
1601:         for (jj=xleft; jj<xright; jj++) idx[loc++] = count++;
1602:       }
1603:       ISCreateGeneral(PETSC_COMM_SELF,nidx,idx,PETSC_COPY_VALUES,(*is)+loc_outer);
1604:       if (overlap == 0) {
1605:         PetscObjectReference((PetscObject)(*is)[loc_outer]);

1607:         (*is_local)[loc_outer] = (*is)[loc_outer];
1608:       } else {
1609:         for (loc=0,ii=ystart; ii<ystart+height; ii++) {
1610:           for (jj=xstart; jj<xstart+width; jj++) {
1611:             idx[loc++] = m*ii + jj;
1612:           }
1613:         }
1614:         ISCreateGeneral(PETSC_COMM_SELF,loc,idx,PETSC_COPY_VALUES,*is_local+loc_outer);
1615:       }
1616:       PetscFree(idx);
1617:       xstart += width;
1618:       loc_outer++;
1619:     }
1620:     ystart += height;
1621:   }
1622:   for (i=0; i<*Nsub; i++) ISSort((*is)[i]);
1623:   return 0;
1624: }

1626: /*@C
1627:     PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
1628:     only) for the additive Schwarz preconditioner.

1630:     Not Collective

1632:     Input Parameter:
1633: .   pc - the preconditioner context

1635:     Output Parameters:
1636: +   n - if requested, the number of subdomains for this processor (default value = 1)
1637: .   is - if requested, the index sets that define the subdomains for this processor
1638: -   is_local - if requested, the index sets that define the local part of the subdomains for this processor (can be NULL)

1640:     Notes:
1641:     The IS numbering is in the parallel, global numbering of the vector.

1643:     Level: advanced

1645: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1646:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubmatrices()
1647: @*/
1648: PetscErrorCode  PCASMGetLocalSubdomains(PC pc,PetscInt *n,IS *is[],IS *is_local[])
1649: {
1650:   PC_ASM         *osm = (PC_ASM*)pc->data;
1651:   PetscBool      match;

1657:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1659:   if (n) *n = osm->n_local_true;
1660:   if (is) *is = osm->is;
1661:   if (is_local) *is_local = osm->is_local;
1662:   return 0;
1663: }

1665: /*@C
1666:     PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1667:     only) for the additive Schwarz preconditioner.

1669:     Not Collective

1671:     Input Parameter:
1672: .   pc - the preconditioner context

1674:     Output Parameters:
1675: +   n - if requested, the number of matrices for this processor (default value = 1)
1676: -   mat - if requested, the matrices

1678:     Level: advanced

1680:     Notes:
1681:     Call after PCSetUp() (or KSPSetUp()) but before PCApply() and before PCSetUpOnBlocks())

1683:            Usually one would use PCSetModifySubMatrices() to change the submatrices in building the preconditioner.

1685: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1686:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains(), PCSetModifySubMatrices()
1687: @*/
1688: PetscErrorCode  PCASMGetLocalSubmatrices(PC pc,PetscInt *n,Mat *mat[])
1689: {
1690:   PC_ASM         *osm;
1691:   PetscBool      match;

1697:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1698:   if (!match) {
1699:     if (n) *n = 0;
1700:     if (mat) *mat = NULL;
1701:   } else {
1702:     osm = (PC_ASM*)pc->data;
1703:     if (n) *n = osm->n_local_true;
1704:     if (mat) *mat = osm->pmat;
1705:   }
1706:   return 0;
1707: }

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

1712:     Logically Collective

1714:     Input Parameters:
1715: +   pc  - the preconditioner
1716: -   flg - boolean indicating whether to use subdomains defined by the DM

1718:     Options Database Key:
1719: .   -pc_asm_dm_subdomains <bool> - use subdomains defined by the DM

1721:     Level: intermediate

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

1727: .seealso: PCASMGetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1728:           PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1729: @*/
1730: PetscErrorCode  PCASMSetDMSubdomains(PC pc,PetscBool flg)
1731: {
1732:   PC_ASM         *osm = (PC_ASM*)pc->data;
1733:   PetscBool      match;

1738:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1739:   if (match) {
1740:     osm->dm_subdomains = flg;
1741:   }
1742:   return 0;
1743: }

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

1749:     Input Parameter:
1750: .   pc  - the preconditioner

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

1755:     Level: intermediate

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

1767:   PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1768:   if (match) *flg = osm->dm_subdomains;
1769:   else *flg = PETSC_FALSE;
1770:   return 0;
1771: }

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

1776:    Not Collective

1778:    Input Parameter:
1779: .  pc - the PC

1781:    Output Parameter:
1782: .  -pc_asm_sub_mat_type - name of matrix type

1784:    Level: advanced

1786: .seealso: PCASMSetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1787: @*/
1788: PetscErrorCode  PCASMGetSubMatType(PC pc,MatType *sub_mat_type)
1789: {
1791:   PetscTryMethod(pc,"PCASMGetSubMatType_C",(PC,MatType*),(pc,sub_mat_type));
1792:   return 0;
1793: }

1795: /*@
1796:      PCASMSetSubMatType - Set the type of matrix used for ASM subsolves

1798:    Collective on Mat

1800:    Input Parameters:
1801: +  pc             - the PC object
1802: -  sub_mat_type   - matrix type

1804:    Options Database Key:
1805: .  -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.

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

1810:   Level: advanced

1812: .seealso: PCASMGetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1813: @*/
1814: PetscErrorCode PCASMSetSubMatType(PC pc,MatType sub_mat_type)
1815: {
1817:   PetscTryMethod(pc,"PCASMSetSubMatType_C",(PC,MatType),(pc,sub_mat_type));
1818:   return 0;
1819: }