Actual source code: gasm.c
petsc-3.3-p7 2013-05-11
1: /*
2: This file defines an "generalized" additive Schwarz preconditioner for any Mat implementation.
3: In this version each processor may intersect multiple subdomains and any subdomain may
4: intersect multiple processors. Intersections of subdomains with processors are called *local
5: subdomains*.
7: N - total number of local subdomains on all processors (set in PCGASMSetTotalSubdomains() or calculated in PCSetUp_GASM())
8: n - actual number of local subdomains on this processor (set in PCGASMSetSubdomains() or calculated in PCGASMSetTotalSubdomains())
9: nmax - maximum number of local subdomains per processor (calculated in PCGASMSetTotalSubdomains() or in PCSetUp_GASM())
10: */
11: #include <petsc-private/pcimpl.h> /*I "petscpc.h" I*/
13: typedef struct {
14: PetscInt N,n,nmax;
15: PetscInt overlap; /* overlap requested by user */
16: KSP *ksp; /* linear solvers for each block */
17: Vec gx,gy; /* Merged work vectors */
18: Vec *x,*y; /* Split work vectors; storage aliases pieces of storage of the above merged vectors. */
19: VecScatter gorestriction; /* merged restriction to disjoint union of outer subdomains */
20: VecScatter girestriction; /* merged restriction to disjoint union of inner subdomains */
21: IS *ois; /* index sets that define the outer (conceptually, overlapping) subdomains */
22: IS *iis; /* index sets that define the inner (conceptually, nonoverlapping) subdomains */
23: Mat *pmat; /* subdomain block matrices */
24: PCGASMType type; /* use reduced interpolation, restriction or both */
25: PetscBool create_local; /* whether the autocreated subdomains are local or not. */
26: PetscBool type_set; /* if user set this value (so won't change it for symmetric problems) */
27: PetscBool same_subdomain_solvers; /* flag indicating whether all local solvers are same */
28: PetscBool sort_indices; /* flag to sort subdomain indices */
29: } PC_GASM;
33: static PetscErrorCode PCGASMSubdomainView_Private(PC pc, PetscInt i, PetscViewer viewer)
34: {
35: PC_GASM *osm = (PC_GASM*)pc->data;
36: PetscInt j,nidx;
37: const PetscInt *idx;
38: PetscViewer sviewer;
39: char *cidx;
43: if(i < 0 || i > osm->n) SETERRQ2(((PetscObject)viewer)->comm, PETSC_ERR_ARG_WRONG, "Invalid subdomain %D: must nonnegative and less than %D", i, osm->n);
44: /* Inner subdomains. */
45: ISGetLocalSize(osm->iis[i], &nidx);
46: /*
47: No more than 15 characters per index plus a space.
48: PetscViewerStringSPrintf requires a string of size at least 2, so use (nidx+1) instead of nidx,
49: in case nidx == 0. That will take care of the space for the trailing '\0' as well.
50: For nidx == 0, the whole string 16 '\0'.
51: */
52: PetscMalloc(sizeof(char)*(16*(nidx+1)+1), &cidx);
53: PetscViewerStringOpen(PETSC_COMM_SELF, cidx, 16*(nidx+1)+1, &sviewer);
54: ISGetIndices(osm->iis[i], &idx);
55: for(j = 0; j < nidx; ++j) {
56: PetscViewerStringSPrintf(sviewer, "%D ", idx[j]);
57: }
58: ISRestoreIndices(osm->iis[i],&idx);
59: PetscViewerDestroy(&sviewer);
60: PetscViewerASCIIPrintf(viewer, "Inner subdomain:\n");
61: PetscViewerFlush(viewer);
62: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
63: PetscViewerASCIISynchronizedPrintf(viewer, "%s", cidx);
64: PetscViewerFlush(viewer);
65: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
66: PetscViewerASCIIPrintf(viewer, "\n");
67: PetscViewerFlush(viewer);
68: PetscFree(cidx);
69: /* Outer subdomains. */
70: ISGetLocalSize(osm->ois[i], &nidx);
71: /*
72: No more than 15 characters per index plus a space.
73: PetscViewerStringSPrintf requires a string of size at least 2, so use (nidx+1) instead of nidx,
74: in case nidx == 0. That will take care of the space for the trailing '\0' as well.
75: For nidx == 0, the whole string 16 '\0'.
76: */
77: PetscMalloc(sizeof(char)*(16*(nidx+1)+1), &cidx);
78: PetscViewerStringOpen(PETSC_COMM_SELF, cidx, 16*(nidx+1)+1, &sviewer);
79: ISGetIndices(osm->ois[i], &idx);
80: for(j = 0; j < nidx; ++j) {
81: PetscViewerStringSPrintf(sviewer,"%D ", idx[j]);
82: }
83: PetscViewerDestroy(&sviewer);
84: ISRestoreIndices(osm->ois[i],&idx);
85: PetscViewerASCIIPrintf(viewer, "Outer subdomain:\n");
86: PetscViewerFlush(viewer);
87: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
88: PetscViewerASCIISynchronizedPrintf(viewer, "%s", cidx);
89: PetscViewerFlush(viewer);
90: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
91: PetscViewerASCIIPrintf(viewer, "\n");
92: PetscViewerFlush(viewer);
93: PetscFree(cidx);
95: return(0);
96: }
100: static PetscErrorCode PCGASMPrintSubdomains(PC pc)
101: {
102: PC_GASM *osm = (PC_GASM*)pc->data;
103: const char *prefix;
104: char fname[PETSC_MAX_PATH_LEN+1];
105: PetscInt i, l, d, count, gcount, *permutation, *numbering;
106: PetscBool found;
107: PetscViewer viewer, sviewer = PETSC_NULL;
111: PetscMalloc2(osm->n, PetscInt, &permutation, osm->n, PetscInt, &numbering);
112: for(i = 0; i < osm->n; ++i) permutation[i] = i;
113: PetscObjectsGetGlobalNumbering(((PetscObject)pc)->comm, osm->n, (PetscObject*)osm->ois, &gcount, numbering);
114: PetscSortIntWithPermutation(osm->n, numbering, permutation);
115: PCGetOptionsPrefix(pc,&prefix);
116: PetscOptionsGetString(prefix,"-pc_gasm_print_subdomains",fname,PETSC_MAX_PATH_LEN,&found);
117: if (!found) { PetscStrcpy(fname,"stdout"); };
118: PetscViewerASCIIOpen(((PetscObject)pc)->comm,fname,&viewer);
119: /*
120: Make sure the viewer has a name. Otherwise this may cause a deadlock or other weird errors when creating a subcomm viewer:
121: the subcomm viewer will attempt to inherit the viewer's name, which, if not set, will be constructed collectively on the comm.
122: */
123: PetscObjectName((PetscObject)viewer);
124: l = 0;
125: for(count = 0; count < gcount; ++count) {
126: /* Now let subdomains go one at a time in the global numbering order and print their subdomain/solver info. */
127: if(l<osm->n){
128: d = permutation[l]; /* d is the local number of the l-th smallest (in the global ordering) among the locally supported subdomains */
129: if(numbering[d] == count) {
130: PetscViewerGetSubcomm(viewer,((PetscObject)osm->ois[d])->comm, &sviewer);
131: PCGASMSubdomainView_Private(pc,d,sviewer);
132: PetscViewerRestoreSubcomm(viewer,((PetscObject)osm->ois[d])->comm, &sviewer);
133: ++l;
134: }
135: }
136: MPI_Barrier(((PetscObject)pc)->comm);
137: }
138: PetscViewerDestroy(&viewer);
139: PetscFree2(permutation,numbering);
140: return(0);
141: }
146: static PetscErrorCode PCView_GASM(PC pc,PetscViewer viewer)
147: {
148: PC_GASM *osm = (PC_GASM*)pc->data;
149: const char *prefix;
151: PetscMPIInt rank, size;
152: PetscInt i,bsz;
153: PetscBool iascii,view_subdomains=PETSC_FALSE;
154: PetscViewer sviewer;
155: PetscInt count, l, gcount, *numbering, *permutation;
156: char overlap[256] = "user-defined overlap";
157: char gsubdomains[256] = "unknown total number of subdomains";
158: char lsubdomains[256] = "unknown number of local subdomains";
159: char msubdomains[256] = "unknown max number of local subdomains";
161: MPI_Comm_size(((PetscObject)pc)->comm, &size);
162: MPI_Comm_rank(((PetscObject)pc)->comm, &rank);
165: PetscMalloc2(osm->n, PetscInt, &permutation, osm->n, PetscInt, &numbering);
166: for(i = 0; i < osm->n; ++i) permutation[i] = i;
167: PetscObjectsGetGlobalNumbering(((PetscObject)pc)->comm, osm->n, (PetscObject*)osm->ois, &gcount, numbering);
168: PetscSortIntWithPermutation(osm->n, numbering, permutation);
170: if(osm->overlap >= 0) {
171: PetscSNPrintf(overlap,sizeof(overlap),"requested amount of overlap = %D",osm->overlap);
172: }
173: PetscSNPrintf(gsubdomains, sizeof(gsubdomains), "total number of subdomains = %D",gcount);
174: if(osm->N > 0) {
175: PetscSNPrintf(lsubdomains, sizeof(gsubdomains), "number of local subdomains = %D",osm->N);
176: }
177: if(osm->nmax > 0){
178: PetscSNPrintf(msubdomains,sizeof(msubdomains),"max number of local subdomains = %D",osm->nmax);
179: }
181: PCGetOptionsPrefix(pc,&prefix);
182: PetscOptionsGetBool(prefix,"-pc_gasm_view_subdomains",&view_subdomains,PETSC_NULL);
184: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
185: if (iascii) {
186: /*
187: Make sure the viewer has a name. Otherwise this may cause a deadlock when creating a subcomm viewer:
188: the subcomm viewer will attempt to inherit the viewer's name, which, if not set, will be constructed
189: collectively on the comm.
190: */
191: PetscObjectName((PetscObject)viewer);
192: PetscViewerASCIIPrintf(viewer,"Generalized additive Schwarz:\n");
193: PetscViewerASCIIPrintf(viewer,"Restriction/interpolation type: %s\n",PCGASMTypes[osm->type]);
194: PetscViewerASCIIPrintf(viewer,"%s\n",overlap);
195: PetscViewerASCIIPrintf(viewer,"%s\n",gsubdomains);
196: PetscViewerASCIIPrintf(viewer,"%s\n",lsubdomains);
197: PetscViewerASCIIPrintf(viewer,"%s\n",msubdomains);
198: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
199: PetscViewerASCIISynchronizedPrintf(viewer,"[%d:%d] number of locally-supported subdomains = %D\n",(int)rank,(int)size,osm->n);
200: PetscViewerFlush(viewer);
201: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
202: /* Cannot take advantage of osm->same_subdomain_solvers without a global numbering of subdomains. */
203: PetscViewerASCIIPrintf(viewer,"Subdomain solver info is as follows:\n");
204: PetscViewerASCIIPushTab(viewer);
205: PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
206: /* Make sure that everybody waits for the banner to be printed. */
207: MPI_Barrier(((PetscObject)viewer)->comm);
208: /* Now let subdomains go one at a time in the global numbering order and print their subdomain/solver info. */
209: l = 0;
210: for(count = 0; count < gcount; ++count) {
211: PetscMPIInt srank, ssize;
212: if(l<osm->n){
213: PetscInt d = permutation[l]; /* d is the local number of the l-th smallest (in the global ordering) among the locally supported subdomains */
214: if(numbering[d] == count) {
215: MPI_Comm_size(((PetscObject)osm->ois[d])->comm, &ssize);
216: MPI_Comm_rank(((PetscObject)osm->ois[d])->comm, &srank);
217: PetscViewerGetSubcomm(viewer,((PetscObject)osm->ois[d])->comm, &sviewer);
218: ISGetLocalSize(osm->ois[d],&bsz);
219: PetscViewerASCIISynchronizedAllow(sviewer,PETSC_TRUE);
220: PetscViewerASCIISynchronizedPrintf(sviewer,"[%D:%D] (subcomm [%D:%D]) local subdomain number %D, local size = %D\n",(int)rank,(int)size,(int)srank,(int)ssize,d,bsz);
221: PetscViewerFlush(sviewer);
222: PetscViewerASCIISynchronizedAllow(sviewer,PETSC_FALSE);
223: if(view_subdomains) {
224: PCGASMSubdomainView_Private(pc,d,sviewer);
225: }
226: if(!pc->setupcalled) {
227: PetscViewerASCIIPrintf(sviewer, "Solver not set up yet: PCSetUp() not yet called\n");
228: }
229: else {
230: KSPView(osm->ksp[d],sviewer);
231: }
232: PetscViewerASCIIPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
233: PetscViewerFlush(sviewer);
234: PetscViewerRestoreSubcomm(viewer,((PetscObject)osm->ois[d])->comm, &sviewer);
235: ++l;
236: }
237: }
238: MPI_Barrier(((PetscObject)pc)->comm);
239: }
240: PetscViewerASCIIPopTab(viewer);
241: } else {
242: SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCGASM",((PetscObject)viewer)->type_name);
243: }
244: PetscFree2(permutation,numbering);
245: return(0);
246: }
254: static PetscErrorCode PCSetUp_GASM(PC pc)
255: {
256: PC_GASM *osm = (PC_GASM*)pc->data;
258: PetscBool symset,flg;
259: PetscInt i;
260: PetscMPIInt rank, size;
261: MatReuse scall = MAT_REUSE_MATRIX;
262: KSP ksp;
263: PC subpc;
264: const char *prefix,*pprefix;
265: Vec x,y;
266: PetscInt oni; /* Number of indices in the i-th local outer subdomain. */
267: const PetscInt *oidxi; /* Indices from the i-th subdomain local outer subdomain. */
268: PetscInt on; /* Number of indices in the disjoint union of local outer subdomains. */
269: PetscInt *oidx; /* Indices in the disjoint union of local outer subdomains. */
270: IS gois; /* Disjoint union the global indices of outer subdomains. */
271: IS goid; /* Identity IS of the size of the disjoint union of outer subdomains. */
272: PetscScalar *gxarray, *gyarray;
273: PetscInt gofirst; /* Start of locally-owned indices in the vectors -- osm->gx,osm->gy --
274: over the disjoint union of outer subdomains. */
275: DM *domain_dm = PETSC_NULL;
278: MPI_Comm_size(((PetscObject)pc)->comm,&size);
279: MPI_Comm_rank(((PetscObject)pc)->comm,&rank);
280: if (!pc->setupcalled) {
282: if (!osm->type_set) {
283: MatIsSymmetricKnown(pc->pmat,&symset,&flg);
284: if (symset && flg) { osm->type = PC_GASM_BASIC; }
285: }
287: /*
288: If subdomains have been set, then the local number of subdomains, osm->n, is NOT PETSC_DECIDE and is at least 1.
289: The total number of subdomains, osm->N is not necessarily set, might be PETSC_DECIDE, and then will have to be calculated from osm->n.
290: */
291: if (osm->n == PETSC_DECIDE) {
292: /* no subdomains given */
293: /* try pc->dm first */
294: if(pc->dm) {
295: char ddm_name[1024];
296: DM ddm;
297: PetscBool flg;
298: PetscInt num_domains, d;
299: char **domain_names;
300: IS *inner_domain_is, *outer_domain_is;
301: /* Allow the user to request a decomposition DM by name */
302: PetscStrncpy(ddm_name, "", 1024);
303: PetscOptionsString("-pc_gasm_decomposition","Name of the DM defining the composition", "PCSetDM", ddm_name, ddm_name,1024,&flg);
304: if(flg) {
305: DMCreateDomainDecompositionDM(pc->dm, ddm_name, &ddm);
306: if(!ddm) {
307: SETERRQ1(((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONGSTATE, "Uknown DM decomposition name %s", ddm_name);
308: }
309: PetscInfo(pc,"Using decomposition DM defined using options database\n");
310: PCSetDM(pc,ddm);
311: }
312: DMCreateDomainDecomposition(pc->dm, &num_domains, &domain_names, &inner_domain_is, &outer_domain_is, &domain_dm);
313: if(num_domains) {
314: PCGASMSetSubdomains(pc, num_domains, inner_domain_is, outer_domain_is);
315: }
316: for(d = 0; d < num_domains; ++d) {
317: PetscFree(domain_names[d]);
318: ISDestroy(&inner_domain_is[d]);
319: ISDestroy(&outer_domain_is[d]);
320: }
321: PetscFree(domain_names);
322: PetscFree(inner_domain_is);
323: PetscFree(outer_domain_is);
324: }
325: if (osm->n == PETSC_DECIDE) { /* still no subdomains; use one per processor */
326: osm->nmax = osm->n = 1;
327: MPI_Comm_size(((PetscObject)pc)->comm,&size);
328: osm->N = size;
329: }
330: }
331: if (osm->N == PETSC_DECIDE) {
332: PetscInt inwork[2], outwork[2];
333: /* determine global number of subdomains and the max number of local subdomains */
334: inwork[0] = inwork[1] = osm->n;
335: MPI_Allreduce(inwork,outwork,1,MPIU_2INT,PetscMaxSum_Op,((PetscObject)pc)->comm);
336: osm->nmax = outwork[0];
337: osm->N = outwork[1];
338: }
339: if (!osm->iis){
340: /*
341: The local number of subdomains was set in PCGASMSetTotalSubdomains() or PCGASMSetSubdomains(),
342: but the actual subdomains have not been supplied (in PCGASMSetSubdomains()).
343: We create the requisite number of inner subdomains on PETSC_COMM_SELF (for now).
344: */
345: PCGASMCreateLocalSubdomains(pc->pmat,osm->overlap,osm->n,&osm->iis,&osm->ois);
346: }
348: PCGetOptionsPrefix(pc,&prefix);
349: flg = PETSC_FALSE;
350: PetscOptionsGetBool(prefix,"-pc_gasm_print_subdomains",&flg,PETSC_NULL);
351: if (flg) { PCGASMPrintSubdomains(pc); }
353: if (osm->sort_indices) {
354: for (i=0; i<osm->n; i++) {
355: ISSort(osm->ois[i]);
356: ISSort(osm->iis[i]);
357: }
358: }
359: /*
360: Merge the ISs, create merged vectors and restrictions.
361: */
362: /* Merge outer subdomain ISs and construct a restriction onto the disjoint union of local outer subdomains. */
363: on = 0;
364: for (i=0; i<osm->n; i++) {
365: ISGetLocalSize(osm->ois[i],&oni);
366: on += oni;
367: }
368: PetscMalloc(on*sizeof(PetscInt), &oidx);
369: on = 0;
370: for (i=0; i<osm->n; i++) {
371: ISGetLocalSize(osm->ois[i],&oni);
372: ISGetIndices(osm->ois[i],&oidxi);
373: PetscMemcpy(oidx+on, oidxi, sizeof(PetscInt)*oni);
374: ISRestoreIndices(osm->ois[i], &oidxi);
375: on += oni;
376: }
377: ISCreateGeneral(((PetscObject)(pc))->comm, on, oidx, PETSC_OWN_POINTER, &gois);
378: MatGetVecs(pc->pmat,&x,&y);
379: VecCreateMPI(((PetscObject)pc)->comm, on, PETSC_DECIDE, &osm->gx);
380: VecDuplicate(osm->gx,&osm->gy);
381: VecGetOwnershipRange(osm->gx, &gofirst, PETSC_NULL);
382: ISCreateStride(((PetscObject)pc)->comm,on,gofirst,1, &goid);
383: VecScatterCreate(x,gois,osm->gx,goid, &(osm->gorestriction));
384: VecDestroy(&x);
385: ISDestroy(&gois);
386: /* Merge inner subdomain ISs and construct a restriction onto the disjoint union of local inner subdomains. */
387: { PetscInt ini; /* Number of indices the i-th a local inner subdomain. */
388: PetscInt in; /* Number of indices in the disjoint uniont of local inner subdomains. */
389: PetscInt *iidx; /* Global indices in the merged local inner subdomain. */
390: PetscInt *ioidx; /* Global indices of the disjoint union of inner subdomains within the disjoint union of outer subdomains. */
391: IS giis; /* IS for the disjoint union of inner subdomains. */
392: IS giois; /* IS for the disjoint union of inner subdomains within the disjoint union of outer subdomains. */
393: /**/
394: in = 0;
395: for (i=0; i<osm->n; i++) {
396: ISGetLocalSize(osm->iis[i],&ini);
397: in += ini;
398: }
399: PetscMalloc(in*sizeof(PetscInt), &iidx);
400: PetscMalloc(in*sizeof(PetscInt), &ioidx);
401: VecGetOwnershipRange(osm->gx,&gofirst, PETSC_NULL);
402: in = 0;
403: on = 0;
404: for (i=0; i<osm->n; i++) {
405: const PetscInt *iidxi; /* Global indices of the i-th local inner subdomain. */
406: ISLocalToGlobalMapping ltogi; /* Map from global to local indices of the i-th outer local subdomain. */
407: PetscInt *ioidxi; /* Local indices of the i-th local inner subdomain within the local outer subdomain. */
408: PetscInt ioni; /* Number of indices in ioidxi; if ioni != ini the inner subdomain is not a subdomain of the outer subdomain (error). */
409: PetscInt k;
410: ISGetLocalSize(osm->iis[i],&ini);
411: ISGetLocalSize(osm->ois[i],&oni);
412: ISGetIndices(osm->iis[i],&iidxi);
413: PetscMemcpy(iidx+in, iidxi, sizeof(PetscInt)*ini);
414: ISLocalToGlobalMappingCreateIS(osm->ois[i],<ogi);
415: ioidxi = ioidx+in;
416: ISGlobalToLocalMappingApply(ltogi,IS_GTOLM_DROP,ini,iidxi,&ioni,ioidxi);
417: ISLocalToGlobalMappingDestroy(<ogi);
418: ISRestoreIndices(osm->iis[i], &iidxi);
419: if (ioni != ini) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inner subdomain %D contains %D indices outside of its outer subdomain", i, ini - ioni);
420: for(k = 0; k < ini; ++k) {
421: ioidxi[k] += gofirst+on;
422: }
423: in += ini;
424: on += oni;
425: }
426: ISCreateGeneral(((PetscObject)pc)->comm, in, iidx, PETSC_OWN_POINTER, &giis);
427: ISCreateGeneral(((PetscObject)pc)->comm, in, ioidx, PETSC_OWN_POINTER, &giois);
428: VecScatterCreate(y,giis,osm->gy,giois,&osm->girestriction);
429: VecDestroy(&y);
430: ISDestroy(&giis);
431: ISDestroy(&giois);
432: }
433: ISDestroy(&goid);
434: /* Create the subdomain work vectors. */
435: PetscMalloc(osm->n*sizeof(Vec),&osm->x);
436: PetscMalloc(osm->n*sizeof(Vec),&osm->y);
437: VecGetArray(osm->gx, &gxarray);
438: VecGetArray(osm->gy, &gyarray);
439: for (i=0, on=0; i<osm->n; ++i, on += oni) {
440: PetscInt oNi;
441: ISGetLocalSize(osm->ois[i],&oni);
442: ISGetSize(osm->ois[i],&oNi);
443: VecCreateMPIWithArray(((PetscObject)(osm->ois[i]))->comm,1,oni,oNi,gxarray+on,&osm->x[i]);
444: VecCreateMPIWithArray(((PetscObject)(osm->ois[i]))->comm,1,oni,oNi,gyarray+on,&osm->y[i]);
445: }
446: VecRestoreArray(osm->gx, &gxarray);
447: VecRestoreArray(osm->gy, &gyarray);
448: /* Create the local solvers */
449: PetscMalloc(osm->n*sizeof(KSP *),&osm->ksp);
450: for (i=0; i<osm->n; i++) { /* KSPs are local */
451: KSPCreate(((PetscObject)(osm->ois[i]))->comm,&ksp);
452: PetscLogObjectParent(pc,ksp);
453: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
454: KSPSetType(ksp,KSPPREONLY);
455: KSPGetPC(ksp,&subpc);
456: PCGetOptionsPrefix(pc,&prefix);
457: KSPSetOptionsPrefix(ksp,prefix);
458: KSPAppendOptionsPrefix(ksp,"sub_");
459: osm->ksp[i] = ksp;
460: }
461: scall = MAT_INITIAL_MATRIX;
463: }/*if(!pc->setupcalled)*/
464: else {
465: /*
466: Destroy the blocks from the previous iteration
467: */
468: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
469: MatDestroyMatrices(osm->n,&osm->pmat);
470: scall = MAT_INITIAL_MATRIX;
471: }
472: }
474: /*
475: Extract out the submatrices.
476: */
477: if(size > 1) {
478: MatGetSubMatricesParallel(pc->pmat,osm->n,osm->ois, osm->ois,scall,&osm->pmat);
479: }
480: else {
481: MatGetSubMatrices(pc->pmat,osm->n,osm->ois, osm->ois,scall,&osm->pmat);
482: }
483: if (scall == MAT_INITIAL_MATRIX) {
484: PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
485: for (i=0; i<osm->n; i++) {
486: PetscLogObjectParent(pc,osm->pmat[i]);
487: PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);
488: }
489: }
490:
491: /* Return control to the user so that the submatrices can be modified (e.g., to apply
492: different boundary conditions for the submatrices than for the global problem) */
493: PCModifySubMatrices(pc,osm->n,osm->ois,osm->ois,osm->pmat,pc->modifysubmatricesP);
495: /*
496: Loop over submatrices putting them into local ksps
497: */
498: for (i=0; i<osm->n; i++) {
499: KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i],pc->flag);
500: if (!pc->setupcalled) {
501: KSPSetFromOptions(osm->ksp[i]);
502: }
503: }
505: return(0);
506: }
510: static PetscErrorCode PCSetUpOnBlocks_GASM(PC pc)
511: {
512: PC_GASM *osm = (PC_GASM*)pc->data;
514: PetscInt i;
517: for (i=0; i<osm->n; i++) {
518: KSPSetUp(osm->ksp[i]);
519: }
520: return(0);
521: }
525: static PetscErrorCode PCApply_GASM(PC pc,Vec x,Vec y)
526: {
527: PC_GASM *osm = (PC_GASM*)pc->data;
529: PetscInt i;
530: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
533: /*
534: Support for limiting the restriction or interpolation only to the inner
535: subdomain values (leaving the other values 0).
536: */
537: if(!(osm->type & PC_GASM_RESTRICT)) {
538: /* have to zero the work RHS since scatter may leave some slots empty */
539: VecZeroEntries(osm->gx);
540: VecScatterBegin(osm->girestriction,x,osm->gx,INSERT_VALUES,forward);
541: }
542: else {
543: VecScatterBegin(osm->gorestriction,x,osm->gx,INSERT_VALUES,forward);
544: }
545: VecZeroEntries(osm->gy);
546: if(!(osm->type & PC_GASM_RESTRICT)) {
547: VecScatterEnd(osm->girestriction,x,osm->gx,INSERT_VALUES,forward);
548: }
549: else {
550: VecScatterEnd(osm->gorestriction,x,osm->gx,INSERT_VALUES,forward);
551: }
552: /* do the subdomain solves */
553: for (i=0; i<osm->n; ++i) {
554: KSPSolve(osm->ksp[i],osm->x[i],osm->y[i]);
555: }
556: VecZeroEntries(y);
557: if(!(osm->type & PC_GASM_INTERPOLATE)) {
558: VecScatterBegin(osm->girestriction,osm->gy,y,ADD_VALUES,reverse);
559: VecScatterEnd(osm->girestriction,osm->gy,y,ADD_VALUES,reverse); return(0);
560: }
561: else {
562: VecScatterBegin(osm->gorestriction,osm->gy,y,ADD_VALUES,reverse);
563: VecScatterEnd(osm->gorestriction,osm->gy,y,ADD_VALUES,reverse); return(0);
564: }
565: }
569: static PetscErrorCode PCApplyTranspose_GASM(PC pc,Vec x,Vec y)
570: {
571: PC_GASM *osm = (PC_GASM*)pc->data;
573: PetscInt i;
574: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
577: /*
578: Support for limiting the restriction or interpolation to only local
579: subdomain values (leaving the other values 0).
581: Note: these are reversed from the PCApply_GASM() because we are applying the
582: transpose of the three terms
583: */
584: if (!(osm->type & PC_GASM_INTERPOLATE)) {
585: /* have to zero the work RHS since scatter may leave some slots empty */
586: VecZeroEntries(osm->gx);
587: VecScatterBegin(osm->girestriction,x,osm->gx,INSERT_VALUES,forward);
588: }
589: else {
590: VecScatterBegin(osm->gorestriction,x,osm->gx,INSERT_VALUES,forward);
591: }
592: VecZeroEntries(osm->gy);
593: if (!(osm->type & PC_GASM_INTERPOLATE)) {
594: VecScatterEnd(osm->girestriction,x,osm->gx,INSERT_VALUES,forward);
595: }
596: else {
597: VecScatterEnd(osm->gorestriction,x,osm->gx,INSERT_VALUES,forward);
598: }
599: /* do the local solves */
600: for (i=0; i<osm->n; ++i) { /* Note that the solves are local, so we can go to osm->n, rather than osm->nmax. */
601: KSPSolveTranspose(osm->ksp[i],osm->x[i],osm->y[i]);
602: }
603: VecZeroEntries(y);
604: if (!(osm->type & PC_GASM_RESTRICT)) {
605: VecScatterBegin(osm->girestriction,osm->gy,y,ADD_VALUES,reverse);
606: VecScatterEnd(osm->girestriction,osm->gy,y,ADD_VALUES,reverse);
607: }
608: else {
609: VecScatterBegin(osm->gorestriction,osm->gy,y,ADD_VALUES,reverse);
610: VecScatterEnd(osm->gorestriction,osm->gy,y,ADD_VALUES,reverse);
611: }
613: return(0);
614: }
618: static PetscErrorCode PCReset_GASM(PC pc)
619: {
620: PC_GASM *osm = (PC_GASM*)pc->data;
622: PetscInt i;
625: if (osm->ksp) {
626: for (i=0; i<osm->n; i++) {
627: KSPReset(osm->ksp[i]);
628: }
629: }
630: if (osm->pmat) {
631: if (osm->n > 0) {
632: MatDestroyMatrices(osm->n,&osm->pmat);
633: }
634: }
635: if (osm->x) {
636: for (i=0; i<osm->n; i++) {
637: VecDestroy(&osm->x[i]);
638: VecDestroy(&osm->y[i]);
639: }
640: }
641: VecDestroy(&osm->gx);
642: VecDestroy(&osm->gy);
643:
644: VecScatterDestroy(&osm->gorestriction);
645: VecScatterDestroy(&osm->girestriction);
646: if (osm->ois) {
647: PCGASMDestroySubdomains(osm->n,osm->ois,osm->iis);
648: osm->ois = 0;
649: osm->iis = 0;
650: }
651: return(0);
652: }
656: static PetscErrorCode PCDestroy_GASM(PC pc)
657: {
658: PC_GASM *osm = (PC_GASM*)pc->data;
660: PetscInt i;
663: PCReset_GASM(pc);
664: if (osm->ksp) {
665: for (i=0; i<osm->n; i++) {
666: KSPDestroy(&osm->ksp[i]);
667: }
668: PetscFree(osm->ksp);
669: }
670: PetscFree(osm->x);
671: PetscFree(osm->y);
672: PetscFree(pc->data);
673: return(0);
674: }
678: static PetscErrorCode PCSetFromOptions_GASM(PC pc) {
679: PC_GASM *osm = (PC_GASM*)pc->data;
681: PetscInt blocks,ovl;
682: PetscBool symset,flg;
683: PCGASMType gasmtype;
686: /* set the type to symmetric if matrix is symmetric */
687: if (!osm->type_set && pc->pmat) {
688: MatIsSymmetricKnown(pc->pmat,&symset,&flg);
689: if (symset && flg) { osm->type = PC_GASM_BASIC; }
690: }
691: PetscOptionsHead("Generalized additive Schwarz options");
692: PetscOptionsInt("-pc_gasm_total_subdomains","Total number of subdomains across communicator","PCGASMSetTotalSubdomains",osm->n,&blocks,&flg);
693: osm->create_local = PETSC_TRUE;
694: PetscOptionsBool("-pc_gasm_subdomains_create_local","Whether to make autocreated subdomains local (true by default)","PCGASMSetTotalSubdomains",osm->create_local,&osm->create_local,&flg);
695: if(!osm->create_local) SETERRQ(((PetscObject)pc)->comm, PETSC_ERR_SUP, "No support for autocreation of nonlocal subdomains yet.");
696:
697: if (flg) {PCGASMSetTotalSubdomains(pc,blocks,osm->create_local); }
698: PetscOptionsInt("-pc_gasm_overlap","Number of overlapping degrees of freedom","PCGASMSetOverlap",osm->overlap,&ovl,&flg);
699: if (flg) {PCGASMSetOverlap(pc,ovl); }
700: flg = PETSC_FALSE;
701: PetscOptionsEnum("-pc_gasm_type","Type of restriction/extension","PCGASMSetType",PCGASMTypes,(PetscEnum)osm->type,(PetscEnum*)&gasmtype,&flg);
702: if (flg) {PCGASMSetType(pc,gasmtype); }
703: PetscOptionsTail();
704: return(0);
705: }
707: /*------------------------------------------------------------------------------------*/
709: EXTERN_C_BEGIN
712: PetscErrorCode PCGASMSetSubdomains_GASM(PC pc,PetscInt n,IS is_local[],IS is[])
713: {
714: PC_GASM *osm = (PC_GASM*)pc->data;
716: PetscInt i;
719: if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Each process must have 1 or more subdomains, n = %D",n);
720: if (pc->setupcalled && (n != osm->n || is_local)) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"PCGASMSetSubdomains() should be called before calling PCSetUp().");
722: if (!pc->setupcalled) {
723: osm->n = n;
724: osm->ois = 0;
725: osm->iis = 0;
726: if (is) {
727: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is[i]);}
728: }
729: if (is_local) {
730: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is_local[i]);}
731: }
732: if (osm->ois) {
733: PCGASMDestroySubdomains(osm->n,osm->iis,osm->ois);
734: }
735: if (is) {
736: PetscMalloc(n*sizeof(IS),&osm->ois);
737: for (i=0; i<n; i++) { osm->ois[i] = is[i]; }
738: /* Flag indicating that the user has set outer subdomains, so PCGASM should not increase their size. */
739: osm->overlap = -1;
740: }
741: if (is_local) {
742: PetscMalloc(n*sizeof(IS),&osm->iis);
743: for (i=0; i<n; i++) { osm->iis[i] = is_local[i]; }
744: }
745: }
746: return(0);
747: }
748: EXTERN_C_END
750: EXTERN_C_BEGIN
753: PetscErrorCode PCGASMSetTotalSubdomains_GASM(PC pc,PetscInt N, PetscBool create_local) {
754: PC_GASM *osm = (PC_GASM*)pc->data;
756: PetscMPIInt rank,size;
757: PetscInt n;
758: PetscInt Nmin, Nmax;
760: if(!create_local) SETERRQ(((PetscObject)pc)->comm, PETSC_ERR_SUP, "No suppor for autocreation of nonlocal subdomains.");
761: if (N < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Total number of subdomains must be > 0, N = %D",N);
762: MPI_Allreduce(&N,&Nmin,1,MPI_INT,MPI_MIN,((PetscObject)pc)->comm);
763: MPI_Allreduce(&N,&Nmax,1,MPI_INT,MPI_MAX,((PetscObject)pc)->comm);
764: if(Nmin != Nmax)
765: SETERRQ2(((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONG, "All processors must use the same number of subdomains. min(N) = %D != %D = max(N)", Nmin, Nmax);
767: osm->create_local = create_local;
768: /*
769: Split the subdomains equally among all processors
770: */
771: MPI_Comm_rank(((PetscObject)pc)->comm,&rank);
772: MPI_Comm_size(((PetscObject)pc)->comm,&size);
773: n = N/size + ((N % size) > rank);
774: if (!n) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Process %d must have at least one subdomain: total processors %d total blocks %D",(int)rank,(int)size,N);
775: if (pc->setupcalled && n != osm->n) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCGASMSetTotalSubdomains() should be called before PCSetUp().");
776: if (!pc->setupcalled) {
777: if (osm->ois) {
778: PCGASMDestroySubdomains(osm->n,osm->iis,osm->ois);
779: }
780: osm->N = N;
781: osm->n = n;
782: osm->nmax = N/size + ((N%size)?1:0);
783: osm->ois = 0;
784: osm->iis = 0;
785: }
786: return(0);
787: }
788: EXTERN_C_END
790: EXTERN_C_BEGIN
793: PetscErrorCode PCGASMSetOverlap_GASM(PC pc,PetscInt ovl)
794: {
795: PC_GASM *osm = (PC_GASM*)pc->data;
798: if (ovl < 0) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested");
799: if (pc->setupcalled && ovl != osm->overlap) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"PCGASMSetOverlap() should be called before PCSetUp().");
800: if (!pc->setupcalled) {
801: osm->overlap = ovl;
802: }
803: return(0);
804: }
805: EXTERN_C_END
807: EXTERN_C_BEGIN
810: PetscErrorCode PCGASMSetType_GASM(PC pc,PCGASMType type)
811: {
812: PC_GASM *osm = (PC_GASM*)pc->data;
815: osm->type = type;
816: osm->type_set = PETSC_TRUE;
817: return(0);
818: }
819: EXTERN_C_END
821: EXTERN_C_BEGIN
824: PetscErrorCode PCGASMSetSortIndices_GASM(PC pc,PetscBool doSort)
825: {
826: PC_GASM *osm = (PC_GASM*)pc->data;
829: osm->sort_indices = doSort;
830: return(0);
831: }
832: EXTERN_C_END
834: EXTERN_C_BEGIN
837: /*
838: FIX: This routine might need to be modified once multiple ranks per subdomain are allowed.
839: In particular, it would upset the global subdomain number calculation.
840: */
841: PetscErrorCode PCGASMGetSubKSP_GASM(PC pc,PetscInt *n,PetscInt *first,KSP **ksp)
842: {
843: PC_GASM *osm = (PC_GASM*)pc->data;
847: if (osm->n < 1) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ORDER,"Need to call PCSetUP() on PC (or KSPSetUp() on the outer KSP object) before calling here");
849: if (n) {
850: *n = osm->n;
851: }
852: if (first) {
853: MPI_Scan(&osm->n,first,1,MPIU_INT,MPI_SUM,((PetscObject)pc)->comm);
854: *first -= osm->n;
855: }
856: if (ksp) {
857: /* Assume that local solves are now different; not necessarily
858: true, though! This flag is used only for PCView_GASM() */
859: *ksp = osm->ksp;
860: osm->same_subdomain_solvers = PETSC_FALSE;
861: }
862: return(0);
863: }/* PCGASMGetSubKSP_GASM() */
864: EXTERN_C_END
869: /*@C
870: PCGASMSetSubdomains - Sets the subdomains for this processor
871: for the additive Schwarz preconditioner.
873: Collective on PC
875: Input Parameters:
876: + pc - the preconditioner context
877: . n - the number of subdomains for this processor
878: . iis - the index sets that define this processor's local inner subdomains
879: (or PETSC_NULL for PETSc to determine subdomains)
880: - ois- the index sets that define this processor's local outer subdomains
881: (or PETSC_NULL to use the same as iis)
883: Notes:
884: The IS indices use the parallel, global numbering of the vector entries.
885: Inner subdomains are those where the correction is applied.
886: Outer subdomains are those where the residual necessary to obtain the
887: corrections is obtained (see PCGASMType for the use of inner/outer subdomains).
888: Both inner and outer subdomains can extend over several processors.
889: This processor's portion of a subdomain is known as a local subdomain.
891: By default the GASM preconditioner uses 1 (local) subdomain per processor.
892: Use PCGASMSetTotalSubdomains() to set the total number of subdomains across
893: all processors that PCGASM will create automatically, and to specify whether
894: they should be local or not.
897: Level: advanced
899: .keywords: PC, GASM, set, subdomains, additive Schwarz
901: .seealso: PCGASMSetTotalSubdomains(), PCGASMSetOverlap(), PCGASMGetSubKSP(),
902: PCGASMCreateSubdomains2D(), PCGASMGetSubdomains()
903: @*/
904: PetscErrorCode PCGASMSetSubdomains(PC pc,PetscInt n,IS iis[],IS ois[])
905: {
910: PetscTryMethod(pc,"PCGASMSetSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,iis,ois));
911: return(0);
912: }
916: /*@C
917: PCGASMSetTotalSubdomains - Sets the total number of subdomains to use in the generalized additive
918: Schwarz preconditioner. The number of subdomains is cumulative across all processors in pc's
919: communicator. Either all or no processors in the PC communicator must call this routine with
920: the same N. The subdomains will be created automatically during PCSetUp().
922: Collective on PC
924: Input Parameters:
925: + pc - the preconditioner context
926: . N - the total number of subdomains cumulative across all processors
927: - create_local - whether the subdomains to be created are to be local
929: Options Database Key:
930: To set the total number of subdomains and let PCGASM autocreate them, rather than specify the index sets, use the following options:
931: + -pc_gasm_total_subdomains <n> - sets the total number of subdomains to be autocreated by PCGASM
932: - -pc_gasm_subdomains_create_local <true|false> - whether autocreated subdomains should be local or not (default is true)
934: By default the GASM preconditioner uses 1 subdomain per processor.
937: Use PCGASMSetSubdomains() to set subdomains explicitly or to set different numbers
938: of subdomains per processor.
940: Level: advanced
942: .keywords: PC, GASM, set, total, global, subdomains, additive Schwarz
944: .seealso: PCGASMSetSubdomains(), PCGASMSetOverlap(), PCGASMGetSubKSP(),
945: PCGASMCreateSubdomains2D()
946: @*/
947: PetscErrorCode PCGASMSetTotalSubdomains(PC pc,PetscInt N, PetscBool create_local)
948: {
953: PetscTryMethod(pc,"PCGASMSetTotalSubdomains_C",(PC,PetscInt,PetscBool),(pc,N,create_local));
954: return(0);
955: }
959: /*@
960: PCGASMSetOverlap - Sets the overlap between a pair of subdomains for the
961: additive Schwarz preconditioner. Either all or no processors in the
962: PC communicator must call this routine.
964: Logically Collective on PC
966: Input Parameters:
967: + pc - the preconditioner context
968: - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1)
970: Options Database Key:
971: . -pc_gasm_overlap <overlap> - Sets overlap
973: Notes:
974: By default the GASM preconditioner uses 1 subdomain per processor. To use
975: multiple subdomain per perocessor, see PCGASMSetTotalSubdomains() or
976: PCGASMSetSubdomains() (and the option -pc_gasm_total_subdomains <n>).
978: The overlap defaults to 1, so if one desires that no additional
979: overlap be computed beyond what may have been set with a call to
980: PCGASMSetTotalSubdomains() or PCGASMSetSubdomains(), then ovl
981: must be set to be 0. In particular, if one does not explicitly set
982: the subdomains in application code, then all overlap would be computed
983: internally by PETSc, and using an overlap of 0 would result in an GASM
984: variant that is equivalent to the block Jacobi preconditioner.
986: Note that one can define initial index sets with any overlap via
987: PCGASMSetSubdomains(); the routine PCGASMSetOverlap() merely allows
988: PETSc to extend that overlap further, if desired.
990: Level: intermediate
992: .keywords: PC, GASM, set, overlap
994: .seealso: PCGASMSetTotalSubdomains(), PCGASMSetSubdomains(), PCGASMGetSubKSP(),
995: PCGASMCreateSubdomains2D(), PCGASMGetSubdomains()
996: @*/
997: PetscErrorCode PCGASMSetOverlap(PC pc,PetscInt ovl)
998: {
1004: PetscTryMethod(pc,"PCGASMSetOverlap_C",(PC,PetscInt),(pc,ovl));
1005: return(0);
1006: }
1010: /*@
1011: PCGASMSetType - Sets the type of restriction and interpolation used
1012: for local problems in the additive Schwarz method.
1014: Logically Collective on PC
1016: Input Parameters:
1017: + pc - the preconditioner context
1018: - type - variant of GASM, one of
1019: .vb
1020: PC_GASM_BASIC - full interpolation and restriction
1021: PC_GASM_RESTRICT - full restriction, local processor interpolation
1022: PC_GASM_INTERPOLATE - full interpolation, local processor restriction
1023: PC_GASM_NONE - local processor restriction and interpolation
1024: .ve
1026: Options Database Key:
1027: . -pc_gasm_type [basic,restrict,interpolate,none] - Sets GASM type
1029: Level: intermediate
1031: .keywords: PC, GASM, set, type
1033: .seealso: PCGASMSetTotalSubdomains(), PCGASMSetSubdomains(), PCGASMGetSubKSP(),
1034: PCGASMCreateSubdomains2D()
1035: @*/
1036: PetscErrorCode PCGASMSetType(PC pc,PCGASMType type)
1037: {
1043: PetscTryMethod(pc,"PCGASMSetType_C",(PC,PCGASMType),(pc,type));
1044: return(0);
1045: }
1049: /*@
1050: PCGASMSetSortIndices - Determines whether subdomain indices are sorted.
1052: Logically Collective on PC
1054: Input Parameters:
1055: + pc - the preconditioner context
1056: - doSort - sort the subdomain indices
1058: Level: intermediate
1060: .keywords: PC, GASM, set, type
1062: .seealso: PCGASMSetSubdomains(), PCGASMSetTotalSubdomains(), PCGASMGetSubKSP(),
1063: PCGASMCreateSubdomains2D()
1064: @*/
1065: PetscErrorCode PCGASMSetSortIndices(PC pc,PetscBool doSort)
1066: {
1072: PetscTryMethod(pc,"PCGASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));
1073: return(0);
1074: }
1078: /*@C
1079: PCGASMGetSubKSP - Gets the local KSP contexts for all blocks on
1080: this processor.
1081:
1082: Collective on PC iff first_local is requested
1084: Input Parameter:
1085: . pc - the preconditioner context
1087: Output Parameters:
1088: + n_local - the number of blocks on this processor or PETSC_NULL
1089: . first_local - the global number of the first block on this processor or PETSC_NULL,
1090: all processors must request or all must pass PETSC_NULL
1091: - ksp - the array of KSP contexts
1093: Note:
1094: After PCGASMGetSubKSP() the array of KSPes is not to be freed
1096: Currently for some matrix implementations only 1 block per processor
1097: is supported.
1098:
1099: You must call KSPSetUp() before calling PCGASMGetSubKSP().
1101: Level: advanced
1103: .keywords: PC, GASM, additive Schwarz, get, sub, KSP, context
1105: .seealso: PCGASMSetTotalSubdomains(), PCGASMSetSubdomains(), PCGASMSetOverlap(),
1106: PCGASMCreateSubdomains2D(),
1107: @*/
1108: PetscErrorCode PCGASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
1109: {
1114: PetscUseMethod(pc,"PCGASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
1115: return(0);
1116: }
1118: /* -------------------------------------------------------------------------------------*/
1119: /*MC
1120: PCGASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1121: its own KSP object.
1123: Options Database Keys:
1124: + -pc_gasm_total_block_count <n> - Sets total number of local subdomains (known as blocks) to be distributed among processors
1125: . -pc_gasm_view_subdomains - activates the printing of subdomain indices in PCView(), -ksp_view or -snes_view
1126: . -pc_gasm_print_subdomains - activates the printing of subdomain indices in PCSetUp()
1127: . -pc_gasm_overlap <ovl> - Sets overlap by which to (automatically) extend local subdomains
1128: - -pc_gasm_type [basic,restrict,interpolate,none] - Sets GASM type
1130: IMPORTANT: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you
1131: will get a different convergence rate due to the default option of -pc_gasm_type restrict. Use
1132: -pc_gasm_type basic to use the standard GASM.
1134: Notes: Each processor can have one or more blocks, but a block cannot be shared by more
1135: than one processor. Defaults to one block per processor.
1137: To set options on the solvers for each block append -sub_ to all the KSP, and PC
1138: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
1139:
1140: To set the options on the solvers separate for each block call PCGASMGetSubKSP()
1141: and set the options directly on the resulting KSP object (you can access its PC
1142: with KSPGetPC())
1145: Level: beginner
1147: Concepts: additive Schwarz method
1149: References:
1150: An additive variant of the Schwarz alternating method for the case of many subregions
1151: M Dryja, OB Widlund - Courant Institute, New York University Technical report
1153: Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations,
1154: Barry Smith, Petter Bjorstad, and William Gropp, Cambridge University Press, ISBN 0-521-49589-X.
1156: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
1157: PCBJACOBI, PCGASMSetUseTrueLocal(), PCGASMGetSubKSP(), PCGASMSetSubdomains(),
1158: PCGASMSetTotalSubdomains(), PCSetModifySubmatrices(), PCGASMSetOverlap(), PCGASMSetType()
1160: M*/
1162: EXTERN_C_BEGIN
1165: PetscErrorCode PCCreate_GASM(PC pc)
1166: {
1168: PC_GASM *osm;
1171: PetscNewLog(pc,PC_GASM,&osm);
1172: osm->N = PETSC_DECIDE;
1173: osm->n = PETSC_DECIDE;
1174: osm->nmax = 0;
1175: osm->overlap = 1;
1176: osm->ksp = 0;
1177: osm->gorestriction = 0;
1178: osm->girestriction = 0;
1179: osm->gx = 0;
1180: osm->gy = 0;
1181: osm->x = 0;
1182: osm->y = 0;
1183: osm->ois = 0;
1184: osm->iis = 0;
1185: osm->pmat = 0;
1186: osm->type = PC_GASM_RESTRICT;
1187: osm->same_subdomain_solvers = PETSC_TRUE;
1188: osm->sort_indices = PETSC_TRUE;
1190: pc->data = (void*)osm;
1191: pc->ops->apply = PCApply_GASM;
1192: pc->ops->applytranspose = PCApplyTranspose_GASM;
1193: pc->ops->setup = PCSetUp_GASM;
1194: pc->ops->reset = PCReset_GASM;
1195: pc->ops->destroy = PCDestroy_GASM;
1196: pc->ops->setfromoptions = PCSetFromOptions_GASM;
1197: pc->ops->setuponblocks = PCSetUpOnBlocks_GASM;
1198: pc->ops->view = PCView_GASM;
1199: pc->ops->applyrichardson = 0;
1201: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMSetSubdomains_C","PCGASMSetSubdomains_GASM",
1202: PCGASMSetSubdomains_GASM);
1203: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMSetTotalSubdomains_C","PCGASMSetTotalSubdomains_GASM",
1204: PCGASMSetTotalSubdomains_GASM);
1205: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMSetOverlap_C","PCGASMSetOverlap_GASM",
1206: PCGASMSetOverlap_GASM);
1207: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMSetType_C","PCGASMSetType_GASM",
1208: PCGASMSetType_GASM);
1209: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMSetSortIndices_C","PCGASMSetSortIndices_GASM",
1210: PCGASMSetSortIndices_GASM);
1211: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCGASMGetSubKSP_C","PCGASMGetSubKSP_GASM",
1212: PCGASMGetSubKSP_GASM);
1213: return(0);
1214: }
1215: EXTERN_C_END
1220: /*@C
1221: PCGASMCreateLocalSubdomains - Creates n local index sets for the overlapping
1222: Schwarz preconditioner for a any problem based on its matrix.
1224: Collective
1226: Input Parameters:
1227: + A - The global matrix operator
1228: . overlap - amount of overlap in outer subdomains
1229: - n - the number of local subdomains
1231: Output Parameters:
1232: + iis - the array of index sets defining the local inner subdomains (on which the correction is applied)
1233: - ois - the array of index sets defining the local outer subdomains (on which the residual is computed)
1235: Level: advanced
1237: Note: this generates n nonoverlapping local inner subdomains on PETSC_COMM_SELF;
1238: PCGASM will generate the overlap from these if you use them in PCGASMSetSubdomains() and set a
1239: nonzero overlap with PCGASMSetOverlap()
1241: In the Fortran version you must provide the array outis[] already allocated of length n.
1243: .keywords: PC, GASM, additive Schwarz, create, subdomains, unstructured grid
1245: .seealso: PCGASMSetSubdomains(), PCGASMDestroySubdomains()
1246: @*/
1247: PetscErrorCode PCGASMCreateLocalSubdomains(Mat A, PetscInt overlap, PetscInt n, IS* iis[], IS* ois[])
1248: {
1249: MatPartitioning mpart;
1250: const char *prefix;
1251: PetscErrorCode (*f)(Mat,MatReuse,Mat*);
1252: PetscMPIInt size;
1253: PetscInt i,j,rstart,rend,bs;
1254: PetscBool isbaij = PETSC_FALSE,foundpart = PETSC_FALSE;
1255: Mat Ad = PETSC_NULL, adj;
1256: IS ispart,isnumb,*is;
1257: PetscErrorCode ierr;
1262: if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of local blocks must be > 0, n = %D",n);
1264: /* Get prefix, row distribution, and block size */
1265: MatGetOptionsPrefix(A,&prefix);
1266: MatGetOwnershipRange(A,&rstart,&rend);
1267: MatGetBlockSize(A,&bs);
1268: 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);
1270: /* Get diagonal block from matrix if possible */
1271: MPI_Comm_size(((PetscObject)A)->comm,&size);
1272: PetscObjectQueryFunction((PetscObject)A,"MatGetDiagonalBlock_C",(void (**)(void))&f);
1273: if (f) {
1274: MatGetDiagonalBlock(A,&Ad);
1275: } else if (size == 1) {
1276: Ad = A;
1277: }
1278: if (Ad) {
1279: PetscObjectTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij);
1280: if (!isbaij) {PetscObjectTypeCompare((PetscObject)Ad,MATSEQSBAIJ,&isbaij);}
1281: }
1282: if (Ad && n > 1) {
1283: PetscBool match,done;
1284: /* Try to setup a good matrix partitioning if available */
1285: MatPartitioningCreate(PETSC_COMM_SELF,&mpart);
1286: PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);
1287: MatPartitioningSetFromOptions(mpart);
1288: PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGCURRENT,&match);
1289: if (!match) {
1290: PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGSQUARE,&match);
1291: }
1292: if (!match) { /* assume a "good" partitioner is available */
1293: PetscInt na,*ia,*ja;
1294: MatGetRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1295: if (done) {
1296: /* Build adjacency matrix by hand. Unfortunately a call to
1297: MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1298: remove the block-aij structure and we cannot expect
1299: MatPartitioning to split vertices as we need */
1300: PetscInt i,j,*row,len,nnz,cnt,*iia=0,*jja=0;
1301: nnz = 0;
1302: for (i=0; i<na; i++) { /* count number of nonzeros */
1303: len = ia[i+1] - ia[i];
1304: row = ja + ia[i];
1305: for (j=0; j<len; j++) {
1306: if (row[j] == i) { /* don't count diagonal */
1307: len--; break;
1308: }
1309: }
1310: nnz += len;
1311: }
1312: PetscMalloc((na+1)*sizeof(PetscInt),&iia);
1313: PetscMalloc((nnz)*sizeof(PetscInt),&jja);
1314: nnz = 0;
1315: iia[0] = 0;
1316: for (i=0; i<na; i++) { /* fill adjacency */
1317: cnt = 0;
1318: len = ia[i+1] - ia[i];
1319: row = ja + ia[i];
1320: for (j=0; j<len; j++) {
1321: if (row[j] != i) { /* if not diagonal */
1322: jja[nnz+cnt++] = row[j];
1323: }
1324: }
1325: nnz += cnt;
1326: iia[i+1] = nnz;
1327: }
1328: /* Partitioning of the adjacency matrix */
1329: MatCreateMPIAdj(PETSC_COMM_SELF,na,na,iia,jja,PETSC_NULL,&adj);
1330: MatPartitioningSetAdjacency(mpart,adj);
1331: MatPartitioningSetNParts(mpart,n);
1332: MatPartitioningApply(mpart,&ispart);
1333: ISPartitioningToNumbering(ispart,&isnumb);
1334: MatDestroy(&adj);
1335: foundpart = PETSC_TRUE;
1336: }
1337: MatRestoreRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1338: }
1339: MatPartitioningDestroy(&mpart);
1340: }
1341: PetscMalloc(n*sizeof(IS),&is);
1342: if (!foundpart) {
1344: /* Partitioning by contiguous chunks of rows */
1346: PetscInt mbs = (rend-rstart)/bs;
1347: PetscInt start = rstart;
1348: for (i=0; i<n; i++) {
1349: PetscInt count = (mbs/n + ((mbs % n) > i)) * bs;
1350: ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);
1351: start += count;
1352: }
1354: } else {
1356: /* Partitioning by adjacency of diagonal block */
1358: const PetscInt *numbering;
1359: PetscInt *count,nidx,*indices,*newidx,start=0;
1360: /* Get node count in each partition */
1361: PetscMalloc(n*sizeof(PetscInt),&count);
1362: ISPartitioningCount(ispart,n,count);
1363: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1364: for (i=0; i<n; i++) count[i] *= bs;
1365: }
1366: /* Build indices from node numbering */
1367: ISGetLocalSize(isnumb,&nidx);
1368: PetscMalloc(nidx*sizeof(PetscInt),&indices);
1369: for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */
1370: ISGetIndices(isnumb,&numbering);
1371: PetscSortIntWithPermutation(nidx,numbering,indices);
1372: ISRestoreIndices(isnumb,&numbering);
1373: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1374: PetscMalloc(nidx*bs*sizeof(PetscInt),&newidx);
1375: for (i=0; i<nidx; i++)
1376: for (j=0; j<bs; j++)
1377: newidx[i*bs+j] = indices[i]*bs + j;
1378: PetscFree(indices);
1379: nidx *= bs;
1380: indices = newidx;
1381: }
1382: /* Shift to get global indices */
1383: for (i=0; i<nidx; i++) indices[i] += rstart;
1385: /* Build the index sets for each block */
1386: for (i=0; i<n; i++) {
1387: ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);
1388: ISSort(is[i]);
1389: start += count[i];
1390: }
1392: PetscFree(count);
1393: PetscFree(indices);
1394: ISDestroy(&isnumb);
1395: ISDestroy(&ispart);
1396: }
1397: *iis = is;
1398: if(!ois) return(0);
1399: /*
1400: Initially make outer subdomains the same as inner subdomains. If nonzero additional overlap
1401: has been requested, copy the inner subdomains over so they can be modified.
1402: */
1403: PetscMalloc(n*sizeof(IS),ois);
1404: for (i=0; i<n; ++i) {
1405: if (overlap > 0) { /* With positive overlap, (*iis)[i] will be modified */
1406: ISDuplicate((*iis)[i],(*ois)+i);
1407: ISCopy((*iis)[i],(*ois)[i]);
1408: } else {
1409: PetscObjectReference((PetscObject)(*iis)[i]);
1410: (*ois)[i] = (*iis)[i];
1411: }
1412: }
1413: if (overlap > 0) {
1414: /* Extend the "overlapping" regions by a number of steps */
1415: MatIncreaseOverlap(A,n,*ois,overlap);
1416: }
1417: return(0);
1418: }
1422: /*@C
1423: PCGASMDestroySubdomains - Destroys the index sets created with
1424: PCGASMCreateLocalSubdomains() or PCGASMCreateSubdomains2D. Should be
1425: called after setting subdomains with PCGASMSetSubdomains().
1427: Collective
1429: Input Parameters:
1430: + n - the number of index sets
1431: . iis - the array of inner subdomains,
1432: - ois - the array of outer subdomains, can be PETSC_NULL
1434: Level: intermediate
1436: Notes: this is merely a convenience subroutine that walks each list,
1437: destroys each IS on the list, and then frees the list.
1439: .keywords: PC, GASM, additive Schwarz, create, subdomains, unstructured grid
1441: .seealso: PCGASMCreateLocalSubdomains(), PCGASMSetSubdomains()
1442: @*/
1443: PetscErrorCode PCGASMDestroySubdomains(PetscInt n, IS iis[], IS ois[])
1444: {
1445: PetscInt i;
1448: if (n <= 0) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"n must be > 0: n = %D",n);
1450: for (i=0; i<n; i++) { ISDestroy(&iis[i]); }
1451: PetscFree(iis);
1452: if (ois) {
1453: for (i=0; i<n; i++) { ISDestroy(&ois[i]); }
1454: PetscFree(ois);
1455: }
1456: return(0);
1457: }
1460: #define PCGASMLocalSubdomainBounds2D(M,N,xleft,ylow,xright,yhigh,first,last,xleft_loc,ylow_loc,xright_loc,yhigh_loc,n) \
1461: { \
1462: PetscInt first_row = first/M, last_row = last/M+1; \
1463: /* \
1464: Compute ylow_loc and yhigh_loc so that (ylow_loc,xleft) and (yhigh_loc,xright) are the corners \
1465: of the bounding box of the intersection of the subdomain with the local ownership range (local \
1466: subdomain). \
1467: Also compute xleft_loc and xright_loc as the lower and upper bounds on the first and last rows \
1468: of the intersection. \
1469: */ \
1470: /* ylow_loc is the grid row containing the first element of the local sumbdomain */ \
1471: *ylow_loc = PetscMax(first_row,ylow); \
1472: /* xleft_loc is the offset of first element of the local subdomain within its grid row (might actually be outside the local subdomain) */ \
1473: *xleft_loc = *ylow_loc==first_row?PetscMax(first%M,xleft):xleft; \
1474: /* yhigh_loc is the grid row above the last local subdomain element */ \
1475: *yhigh_loc = PetscMin(last_row,yhigh); \
1476: /* xright is the offset of the end of the local subdomain within its grid row (might actually be outside the local subdomain) */ \
1477: *xright_loc = *yhigh_loc==last_row?PetscMin(xright,last%M):xright; \
1478: /* Now compute the size of the local subdomain n. */ \
1479: *n = 0; \
1480: if(*ylow_loc < *yhigh_loc) { \
1481: PetscInt width = xright-xleft; \
1482: *n += width*(*yhigh_loc-*ylow_loc-1); \
1483: *n += PetscMin(PetscMax(*xright_loc-xleft,0),width); \
1484: *n -= PetscMin(PetscMax(*xleft_loc-xleft,0), width); \
1485: }\
1486: }
1492: /*@
1493: PCGASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1494: preconditioner for a two-dimensional problem on a regular grid.
1496: Collective
1498: Input Parameters:
1499: + M, N - the global number of grid points in the x and y directions
1500: . Mdomains, Ndomains - the global number of subdomains in the x and y directions
1501: . dof - degrees of freedom per node
1502: - overlap - overlap in mesh lines
1504: Output Parameters:
1505: + Nsub - the number of local subdomains created
1506: . iis - array of index sets defining inner (nonoverlapping) subdomains
1507: - ois - array of index sets defining outer (overlapping, if overlap > 0) subdomains
1510: Level: advanced
1512: .keywords: PC, GASM, additive Schwarz, create, subdomains, 2D, regular grid
1514: .seealso: PCGASMSetTotalSubdomains(), PCGASMSetSubdomains(), PCGASMGetSubKSP(),
1515: PCGASMSetOverlap()
1516: @*/
1517: PetscErrorCode PCGASMCreateSubdomains2D(PC pc, PetscInt M,PetscInt N,PetscInt Mdomains,PetscInt Ndomains,PetscInt dof,PetscInt overlap, PetscInt *nsub,IS **iis,IS **ois)
1518: {
1520: PetscMPIInt size, rank;
1521: PetscInt i, j;
1522: PetscInt maxheight, maxwidth;
1523: PetscInt xstart, xleft, xright, xleft_loc, xright_loc;
1524: PetscInt ystart, ylow, yhigh, ylow_loc, yhigh_loc;
1525: PetscInt x[2][2], y[2][2], n[2];
1526: PetscInt first, last;
1527: PetscInt nidx, *idx;
1528: PetscInt ii,jj,s,q,d;
1529: PetscInt k,kk;
1530: PetscMPIInt color;
1531: MPI_Comm comm, subcomm;
1532: IS **xis = 0, **is = ois, **is_local = iis;
1535: PetscObjectGetComm((PetscObject)pc, &comm);
1536: MPI_Comm_size(comm, &size);
1537: MPI_Comm_rank(comm, &rank);
1538: MatGetOwnershipRange(pc->pmat, &first, &last);
1539: if (first%dof || last%dof) SETERRQ3(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Matrix row partitioning unsuitable for domain decomposition: local row range (%D,%D) "
1540: "does not respect the number of degrees of freedom per grid point %D", first, last, dof);
1542: /* Determine the number of domains with nonzero intersections with the local ownership range. */
1543: s = 0;
1544: ystart = 0;
1545: for (j=0; j<Ndomains; ++j) {
1546: maxheight = N/Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */
1547: if (maxheight < 2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many %D subdomains in the vertical directon for mesh height %D", Ndomains, N);
1548: /* Vertical domain limits with an overlap. */
1549: ylow = PetscMax(ystart - overlap,0);
1550: yhigh = PetscMin(ystart + maxheight + overlap,N);
1551: xstart = 0;
1552: for (i=0; i<Mdomains; ++i) {
1553: maxwidth = M/Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */
1554: if (maxwidth < 2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many %D subdomains in the horizontal direction for mesh width %D", Mdomains, M);
1555: /* Horizontal domain limits with an overlap. */
1556: xleft = PetscMax(xstart - overlap,0);
1557: xright = PetscMin(xstart + maxwidth + overlap,M);
1558: /*
1559: Determine whether this subdomain intersects this processor's ownership range of pc->pmat.
1560: */
1561: PCGASMLocalSubdomainBounds2D(M,N,xleft,ylow,xright,yhigh,first,last,(&xleft_loc),(&ylow_loc),(&xright_loc),(&yhigh_loc),(&nidx));
1562: if(nidx) {
1563: ++s;
1564: }
1565: xstart += maxwidth;
1566: }/* for(i = 0; i < Mdomains; ++i) */
1567: ystart += maxheight;
1568: }/* for(j = 0; j < Ndomains; ++j) */
1569: /* Now we can allocate the necessary number of ISs. */
1570: *nsub = s;
1571: PetscMalloc((*nsub)*sizeof(IS*),is);
1572: PetscMalloc((*nsub)*sizeof(IS*),is_local);
1573: s = 0;
1574: ystart = 0;
1575: for (j=0; j<Ndomains; ++j) {
1576: maxheight = N/Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */
1577: if (maxheight < 2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many %D subdomains in the vertical directon for mesh height %D", Ndomains, N);
1578: /* Vertical domain limits with an overlap. */
1579: y[0][0] = PetscMax(ystart - overlap,0);
1580: y[0][1] = PetscMin(ystart + maxheight + overlap,N);
1581: /* Vertical domain limits without an overlap. */
1582: y[1][0] = ystart;
1583: y[1][1] = PetscMin(ystart + maxheight,N);
1584: xstart = 0;
1585: for (i=0; i<Mdomains; ++i) {
1586: maxwidth = M/Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */
1587: if (maxwidth < 2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many %D subdomains in the horizontal direction for mesh width %D", Mdomains, M);
1588: /* Horizontal domain limits with an overlap. */
1589: x[0][0] = PetscMax(xstart - overlap,0);
1590: x[0][1] = PetscMin(xstart + maxwidth + overlap,M);
1591: /* Horizontal domain limits without an overlap. */
1592: x[1][0] = xstart;
1593: x[1][1] = PetscMin(xstart+maxwidth,M);
1594: /*
1595: Determine whether this domain intersects this processor's ownership range of pc->pmat.
1596: Do this twice: first for the domains with overlaps, and once without.
1597: During the first pass create the subcommunicators, and use them on the second pass as well.
1598: */
1599: for(q = 0; q < 2; ++q) {
1600: /*
1601: domain limits, (xleft, xright) and (ylow, yheigh) are adjusted
1602: according to whether the domain with an overlap or without is considered.
1603: */
1604: xleft = x[q][0]; xright = x[q][1];
1605: ylow = y[q][0]; yhigh = y[q][1];
1606: PCGASMLocalSubdomainBounds2D(M,N,xleft,ylow,xright,yhigh,first,last,(&xleft_loc),(&ylow_loc),(&xright_loc),(&yhigh_loc),(&nidx));
1607: nidx *= dof;
1608: n[q] = nidx;
1609: /*
1610: Based on the counted number of indices in the local domain *with an overlap*,
1611: construct a subcommunicator of all the processors supporting this domain.
1612: Observe that a domain with an overlap might have nontrivial local support,
1613: while the domain without an overlap might not. Hence, the decision to participate
1614: in the subcommunicator must be based on the domain with an overlap.
1615: */
1616: if (q == 0) {
1617: if(nidx) {
1618: color = 1;
1619: } else {
1620: color = MPI_UNDEFINED;
1621: }
1622: MPI_Comm_split(comm, color, rank, &subcomm);
1623: }
1624: /*
1625: Proceed only if the number of local indices *with an overlap* is nonzero.
1626: */
1627: if (n[0]) {
1628: if(q == 0) {
1629: xis = is;
1630: }
1631: if (q == 1) {
1632: /*
1633: The IS for the no-overlap subdomain shares a communicator with the overlapping domain.
1634: Moreover, if the overlap is zero, the two ISs are identical.
1635: */
1636: if (overlap == 0) {
1637: (*is_local)[s] = (*is)[s];
1638: PetscObjectReference((PetscObject)(*is)[s]);
1639: continue;
1640: } else {
1641: xis = is_local;
1642: subcomm = ((PetscObject)(*is)[s])->comm;
1643: }
1644: }/* if(q == 1) */
1645: idx = PETSC_NULL;
1646: PetscMalloc(nidx*sizeof(PetscInt),&idx);
1647: if(nidx) {
1648: k = 0;
1649: for (jj=ylow_loc; jj<yhigh_loc; ++jj) {
1650: PetscInt x0 = (jj==ylow_loc)?xleft_loc:xleft;
1651: PetscInt x1 = (jj==yhigh_loc-1)?xright_loc:xright;
1652: kk = dof*(M*jj + x0);
1653: for (ii=x0; ii<x1; ++ii) {
1654: for(d = 0; d < dof; ++d) {
1655: idx[k++] = kk++;
1656: }
1657: }
1658: }
1659: }
1660: ISCreateGeneral(subcomm,nidx,idx,PETSC_OWN_POINTER,(*xis)+s);
1661: }/* if(n[0]) */
1662: }/* for(q = 0; q < 2; ++q) */
1663: if(n[0]) {
1664: ++s;
1665: }
1666: xstart += maxwidth;
1667: }/* for(i = 0; i < Mdomains; ++i) */
1668: ystart += maxheight;
1669: }/* for(j = 0; j < Ndomains; ++j) */
1670: return(0);
1671: }
1675: /*@C
1676: PCGASMGetSubdomains - Gets the subdomains supported on this processor
1677: for the additive Schwarz preconditioner.
1679: Not Collective
1681: Input Parameter:
1682: . pc - the preconditioner context
1684: Output Parameters:
1685: + n - the number of subdomains for this processor (default value = 1)
1686: . iis - the index sets that define the inner subdomains (without overlap) supported on this processor (can be PETSC_NULL)
1687: - ois - the index sets that define the outer subdomains (with overlap) supported on this processor (can be PETSC_NULL)
1688:
1690: Notes:
1691: The user is responsible for destroying the ISs and freeing the returned arrays.
1692: The IS numbering is in the parallel, global numbering of the vector.
1694: Level: advanced
1696: .keywords: PC, GASM, get, subdomains, additive Schwarz
1698: .seealso: PCGASMSetTotalSubdomains(), PCGASMSetOverlap(), PCGASMGetSubKSP(),
1699: PCGASMCreateSubdomains2D(), PCGASMSetSubdomains(), PCGASMGetSubmatrices()
1700: @*/
1701: PetscErrorCode PCGASMGetSubdomains(PC pc,PetscInt *n,IS *iis[],IS *ois[])
1702: {
1703: PC_GASM *osm;
1705: PetscBool match;
1706: PetscInt i;
1709: PetscObjectTypeCompare((PetscObject)pc,PCGASM,&match);
1710: if (!match)
1711: SETERRQ2(((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONG, "Incorrect object type: expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name);
1712: osm = (PC_GASM*)pc->data;
1713: if (n) *n = osm->n;
1714: if(iis) {
1715: PetscMalloc(osm->n*sizeof(IS), iis);
1716: }
1717: if(ois) {
1718: PetscMalloc(osm->n*sizeof(IS), ois);
1719: }
1720: if(iis || ois) {
1721: for(i = 0; i < osm->n; ++i) {
1722: if(iis) (*iis)[i] = osm->iis[i];
1723: if(ois) (*ois)[i] = osm->ois[i];
1724: }
1725: }
1726: return(0);
1727: }
1731: /*@C
1732: PCGASMGetSubmatrices - Gets the local submatrices (for this processor
1733: only) for the additive Schwarz preconditioner.
1735: Not Collective
1737: Input Parameter:
1738: . pc - the preconditioner context
1740: Output Parameters:
1741: + n - the number of matrices for this processor (default value = 1)
1742: - mat - the matrices
1743:
1744: Notes: matrices returned by this routine have the same communicators as the index sets (IS)
1745: used to define subdomains in PCGASMSetSubdomains(), or PETSC_COMM_SELF, if the
1746: subdomains were defined using PCGASMSetTotalSubdomains().
1747: Level: advanced
1749: .keywords: PC, GASM, set, local, subdomains, additive Schwarz, block Jacobi
1751: .seealso: PCGASMSetTotalSubdomain(), PCGASMSetOverlap(), PCGASMGetSubKSP(),
1752: PCGASMCreateSubdomains2D(), PCGASMSetSubdomains(), PCGASMGetSubdomains()
1753: @*/
1754: PetscErrorCode PCGASMGetSubmatrices(PC pc,PetscInt *n,Mat *mat[])
1755: {
1756: PC_GASM *osm;
1758: PetscBool match;
1764: if (!pc->setupcalled) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUP() or PCSetUp().");
1765: PetscObjectTypeCompare((PetscObject)pc,PCGASM,&match);
1766: if (!match) SETERRQ2(((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONG, "Expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name);
1767: osm = (PC_GASM*)pc->data;
1768: if (n) *n = osm->n;
1769: if (mat) *mat = osm->pmat;
1771: return(0);
1772: }