Actual source code: gasm.c
petsc-3.6.1 2015-08-06
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 distinct global subdomains (set explicitly in PCGASMSetTotalSubdomains() or implicitly PCGASMSetSubdomains() and then 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 PCSetUp_GASM())
10: */
11: #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
12: #include <petscdm.h>
14: typedef struct {
15: PetscInt N,n,nmax;
16: PetscInt overlap; /* overlap requested by user */
17: PCGASMType type; /* use reduced interpolation, restriction or both */
18: PetscBool type_set; /* if user set this value (so won't change it for symmetric problems) */
19: PetscBool same_subdomain_solvers; /* flag indicating whether all local solvers are same */
20: PetscBool sort_indices; /* flag to sort subdomain indices */
21: PetscBool user_subdomains; /* whether the user set explicit subdomain index sets -- keep them on PCReset() */
22: PetscBool dm_subdomains; /* whether DM is allowed to define subdomains */
23: IS *ois; /* index sets that define the outer (conceptually, overlapping) subdomains */
24: IS *iis; /* index sets that define the inner (conceptually, nonoverlapping) subdomains */
25: KSP *ksp; /* linear solvers for each subdomain */
26: Mat *pmat; /* subdomain block matrices */
27: Vec gx,gy; /* Merged work vectors */
28: Vec *x,*y; /* Split work vectors; storage aliases pieces of storage of the above merged vectors. */
29: VecScatter gorestriction; /* merged restriction to disjoint union of outer subdomains */
30: VecScatter girestriction; /* merged restriction to disjoint union of inner subdomains */
31: } PC_GASM;
35: static PetscErrorCode PCGASMComputeGlobalSubdomainNumbering_Private(PC pc,PetscInt **numbering,PetscInt **permutation)
36: {
37: PC_GASM *osm = (PC_GASM*)pc->data;
38: PetscInt i;
42: /* Determine the number of globally-distinct subdomains and compute a global numbering for them. */
43: PetscMalloc2(osm->n,numbering,osm->n,permutation);
44: PetscObjectsListGetGlobalNumbering(PetscObjectComm((PetscObject)pc),osm->n,(PetscObject*)osm->iis,NULL,*numbering);
45: for (i = 0; i < osm->n; ++i) (*permutation)[i] = i;
46: PetscSortIntWithPermutation(osm->n,*numbering,*permutation);
47: return(0);
48: }
52: static PetscErrorCode PCGASMSubdomainView_Private(PC pc, PetscInt i, PetscViewer viewer)
53: {
54: PC_GASM *osm = (PC_GASM*)pc->data;
55: PetscInt j,nidx;
56: const PetscInt *idx;
57: PetscViewer sviewer;
58: char *cidx;
62: if (i < 0 || i > osm->n) SETERRQ2(PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_WRONG, "Invalid subdomain %D: must nonnegative and less than %D", i, osm->n);
63: /* Inner subdomains. */
64: ISGetLocalSize(osm->iis[i], &nidx);
65: /*
66: No more than 15 characters per index plus a space.
67: PetscViewerStringSPrintf requires a string of size at least 2, so use (nidx+1) instead of nidx,
68: in case nidx == 0. That will take care of the space for the trailing '\0' as well.
69: For nidx == 0, the whole string 16 '\0'.
70: */
71: PetscMalloc1(16*(nidx+1)+1, &cidx);
72: PetscViewerStringOpen(PETSC_COMM_SELF, cidx, 16*(nidx+1)+1, &sviewer);
73: ISGetIndices(osm->iis[i], &idx);
74: for (j = 0; j < nidx; ++j) {
75: PetscViewerStringSPrintf(sviewer, "%D ", idx[j]);
76: }
77: ISRestoreIndices(osm->iis[i],&idx);
78: PetscViewerDestroy(&sviewer);
79: PetscViewerASCIIPrintf(viewer, "Inner subdomain:\n");
80: PetscViewerFlush(viewer);
81: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
82: PetscViewerASCIISynchronizedPrintf(viewer, "%s", cidx);
83: PetscViewerFlush(viewer);
84: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
85: PetscViewerASCIIPrintf(viewer, "\n");
86: PetscViewerFlush(viewer);
87: PetscFree(cidx);
88: /* Outer subdomains. */
89: ISGetLocalSize(osm->ois[i], &nidx);
90: /*
91: No more than 15 characters per index plus a space.
92: PetscViewerStringSPrintf requires a string of size at least 2, so use (nidx+1) instead of nidx,
93: in case nidx == 0. That will take care of the space for the trailing '\0' as well.
94: For nidx == 0, the whole string 16 '\0'.
95: */
96: PetscMalloc1(16*(nidx+1)+1, &cidx);
97: PetscViewerStringOpen(PETSC_COMM_SELF, cidx, 16*(nidx+1)+1, &sviewer);
98: ISGetIndices(osm->ois[i], &idx);
99: for (j = 0; j < nidx; ++j) {
100: PetscViewerStringSPrintf(sviewer,"%D ", idx[j]);
101: }
102: PetscViewerDestroy(&sviewer);
103: ISRestoreIndices(osm->ois[i],&idx);
104: PetscViewerASCIIPrintf(viewer, "Outer subdomain:\n");
105: PetscViewerFlush(viewer);
106: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
107: PetscViewerASCIISynchronizedPrintf(viewer, "%s", cidx);
108: PetscViewerFlush(viewer);
109: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
110: PetscViewerASCIIPrintf(viewer, "\n");
111: PetscViewerFlush(viewer);
112: PetscFree(cidx);
113: return(0);
114: }
118: static PetscErrorCode PCGASMPrintSubdomains(PC pc)
119: {
120: PC_GASM *osm = (PC_GASM*)pc->data;
121: const char *prefix;
122: char fname[PETSC_MAX_PATH_LEN+1];
123: PetscInt l, d, count;
124: PetscBool doprint,found;
125: PetscViewer viewer, sviewer = NULL;
126: PetscInt *numbering,*permutation;/* global numbering of locally-supported subdomains and the permutation from the local ordering */
130: PCGetOptionsPrefix(pc,&prefix);
131: doprint = PETSC_FALSE;
132: PetscOptionsGetBool(prefix,"-pc_gasm_print_subdomains",&doprint,NULL);
133: if (!doprint) return(0);
134: PetscOptionsGetString(prefix,"-pc_gasm_print_subdomains",fname,PETSC_MAX_PATH_LEN,&found);
135: if (!found) { PetscStrcpy(fname,"stdout"); };
136: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),fname,&viewer);
137: /*
138: Make sure the viewer has a name. Otherwise this may cause a deadlock or other weird errors when creating a subcomm viewer:
139: the subcomm viewer will attempt to inherit the viewer's name, which, if not set, will be constructed collectively on the comm.
140: */
141: PetscObjectName((PetscObject)viewer);
142: l = 0;
143: PCGASMComputeGlobalSubdomainNumbering_Private(pc,&numbering,&permutation);
144: for (count = 0; count < osm->N; ++count) {
145: /* Now let subdomains go one at a time in the global numbering order and print their subdomain/solver info. */
146: if (l<osm->n) {
147: d = permutation[l]; /* d is the local number of the l-th smallest (in the global ordering) among the locally supported subdomains */
148: if (numbering[d] == count) {
149: PetscViewerGetSubcomm(viewer,((PetscObject)osm->ois[d])->comm, &sviewer);
150: PCGASMSubdomainView_Private(pc,d,sviewer);
151: PetscViewerRestoreSubcomm(viewer,((PetscObject)osm->ois[d])->comm, &sviewer);
152: ++l;
153: }
154: }
155: MPI_Barrier(PetscObjectComm((PetscObject)pc));
156: }
157: PetscFree2(numbering,permutation);
158: PetscViewerDestroy(&viewer);
159: return(0);
160: }
165: static PetscErrorCode PCView_GASM(PC pc,PetscViewer viewer)
166: {
167: PC_GASM *osm = (PC_GASM*)pc->data;
168: const char *prefix;
170: PetscMPIInt rank, size;
171: PetscInt bsz;
172: PetscBool iascii,view_subdomains=PETSC_FALSE;
173: PetscViewer sviewer;
174: PetscInt count, l;
175: char overlap[256] = "user-defined overlap";
176: char gsubdomains[256] = "unknown total number of subdomains";
177: char msubdomains[256] = "unknown max number of local subdomains";
178: PetscInt *numbering,*permutation;/* global numbering of locally-supported subdomains and the permutation from the local ordering */
181: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
182: MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
184: if (osm->overlap >= 0) {
185: PetscSNPrintf(overlap,sizeof(overlap),"requested amount of overlap = %D",osm->overlap);
186: }
187: if (osm->N != PETSC_DETERMINE) {
188: PetscSNPrintf(gsubdomains, sizeof(gsubdomains), "total number of subdomains = %D",osm->N);
189: }
190: if (osm->nmax != PETSC_DETERMINE) {
191: PetscSNPrintf(msubdomains,sizeof(msubdomains),"max number of local subdomains = %D",osm->nmax);
192: }
194: PCGetOptionsPrefix(pc,&prefix);
195: PetscOptionsGetBool(prefix,"-pc_gasm_view_subdomains",&view_subdomains,NULL);
197: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
198: if (iascii) {
199: /*
200: Make sure the viewer has a name. Otherwise this may cause a deadlock when creating a subcomm viewer:
201: the subcomm viewer will attempt to inherit the viewer's name, which, if not set, will be constructed
202: collectively on the comm.
203: */
204: PetscObjectName((PetscObject)viewer);
205: PetscViewerASCIIPrintf(viewer,"Generalized additive Schwarz:\n");
206: PetscViewerASCIIPrintf(viewer,"Restriction/interpolation type: %s\n",PCGASMTypes[osm->type]);
207: PetscViewerASCIIPrintf(viewer,"%s\n",overlap);
208: PetscViewerASCIIPrintf(viewer,"%s\n",gsubdomains);
209: PetscViewerASCIIPrintf(viewer,"%s\n",msubdomains);
210: PetscViewerASCIISynchronizedAllow(viewer,PETSC_TRUE);
211: PetscViewerASCIISynchronizedPrintf(viewer,"[%d|%d] number of locally-supported subdomains = %D\n",rank,size,osm->n);
212: PetscViewerFlush(viewer);
213: PetscViewerASCIISynchronizedAllow(viewer,PETSC_FALSE);
214: /* Cannot take advantage of osm->same_subdomain_solvers without a global numbering of subdomains. */
215: PetscViewerASCIIPrintf(viewer,"Subdomain solver info is as follows:\n");
216: PetscViewerASCIIPushTab(viewer);
217: PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
218: /* Make sure that everybody waits for the banner to be printed. */
219: MPI_Barrier(PetscObjectComm((PetscObject)viewer));
220: /* Now let subdomains go one at a time in the global numbering order and print their subdomain/solver info. */
221: PCGASMComputeGlobalSubdomainNumbering_Private(pc,&numbering,&permutation);
222: l = 0;
223: for (count = 0; count < osm->N; ++count) {
224: PetscMPIInt srank, ssize;
225: if (l<osm->n) {
226: PetscInt d = permutation[l]; /* d is the local number of the l-th smallest (in the global ordering) among the locally supported subdomains */
227: if (numbering[d] == count) {
228: MPI_Comm_size(((PetscObject)osm->ois[d])->comm, &ssize);
229: MPI_Comm_rank(((PetscObject)osm->ois[d])->comm, &srank);
230: PetscViewerGetSubcomm(viewer,((PetscObject)osm->ois[d])->comm, &sviewer);
231: ISGetLocalSize(osm->ois[d],&bsz);
232: PetscViewerASCIISynchronizedAllow(sviewer,PETSC_TRUE);
233: PetscViewerASCIISynchronizedPrintf(sviewer,"[%d|%d] (subcomm [%d|%d]) local subdomain number %D, local size = %D\n",rank,size,srank,ssize,d,bsz);
234: PetscViewerFlush(sviewer);
235: PetscViewerASCIISynchronizedAllow(sviewer,PETSC_FALSE);
236: if (view_subdomains) {
237: PCGASMSubdomainView_Private(pc,d,sviewer);
238: }
239: if (!pc->setupcalled) {
240: PetscViewerASCIIPrintf(sviewer, "Solver not set up yet: PCSetUp() not yet called\n");
241: } else {
242: KSPView(osm->ksp[d],sviewer);
243: }
244: PetscViewerASCIIPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
245: PetscViewerFlush(sviewer);
246: PetscViewerRestoreSubcomm(viewer,((PetscObject)osm->ois[d])->comm, &sviewer);
247: ++l;
248: }
249: }
250: MPI_Barrier(PetscObjectComm((PetscObject)pc));
251: }
252: PetscFree2(numbering,permutation);
253: PetscViewerASCIIPopTab(viewer);
254: }
255: return(0);
256: }
258: PETSC_INTERN PetscErrorCode PCGASMCreateLocalSubdomains(Mat A, PetscInt nloc, IS *iis[]);
262: static PetscErrorCode PCSetUp_GASM(PC pc)
263: {
264: PC_GASM *osm = (PC_GASM*)pc->data;
266: PetscBool symset,flg;
267: PetscInt i;
268: PetscMPIInt rank, size;
269: MatReuse scall = MAT_REUSE_MATRIX;
270: KSP ksp;
271: PC subpc;
272: const char *prefix,*pprefix;
273: Vec x,y;
274: PetscInt oni; /* Number of indices in the i-th local outer subdomain. */
275: const PetscInt *oidxi; /* Indices from the i-th subdomain local outer subdomain. */
276: PetscInt on; /* Number of indices in the disjoint union of local outer subdomains. */
277: PetscInt *oidx; /* Indices in the disjoint union of local outer subdomains. */
278: IS gois; /* Disjoint union the global indices of outer subdomains. */
279: IS goid; /* Identity IS of the size of the disjoint union of outer subdomains. */
280: PetscScalar *gxarray, *gyarray;
281: PetscInt gofirst; /* Start of locally-owned indices in the vectors -- osm->gx,osm->gy -- over the disjoint union of outer subdomains. */
282: PetscInt num_subdomains = 0;
283: DM *subdomain_dm = NULL;
284: char **subdomain_names = NULL;
288: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
289: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
290: if (!pc->setupcalled) {
292: if (!osm->type_set) {
293: MatIsSymmetricKnown(pc->pmat,&symset,&flg);
294: if (symset && flg) osm->type = PC_GASM_BASIC;
295: }
297: if (osm->n == PETSC_DETERMINE) {
298: if (osm->N != PETSC_DETERMINE) {
299: /* No local subdomains given, but the desired number of total subdomains is known, so construct them accordingly. */
300: PCGASMCreateSubdomains(pc->pmat,osm->N,&osm->n,&osm->iis);
301: } else if (osm->dm_subdomains && pc->dm) {
302: /* try pc->dm next, if allowed */
303: PetscInt d;
304: IS *inner_subdomain_is, *outer_subdomain_is;
305: DMCreateDomainDecomposition(pc->dm, &num_subdomains, &subdomain_names, &inner_subdomain_is, &outer_subdomain_is, &subdomain_dm);
306: if (num_subdomains) {
307: PCGASMSetSubdomains(pc, num_subdomains, inner_subdomain_is, outer_subdomain_is);
308: }
309: for (d = 0; d < num_subdomains; ++d) {
310: if (inner_subdomain_is) {ISDestroy(&inner_subdomain_is[d]);}
311: if (outer_subdomain_is) {ISDestroy(&outer_subdomain_is[d]);}
312: }
313: PetscFree(inner_subdomain_is);
314: PetscFree(outer_subdomain_is);
315: } else {
316: /* still no subdomains; use one per processor */
317: osm->nmax = osm->n = 1;
318: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
319: osm->N = size;
320: PCGASMCreateLocalSubdomains(pc->pmat,osm->n,&osm->iis);
321: }
322: }
323: if (!osm->iis) {
324: /*
325: osm->n was set in PCGASMSetSubdomains(), but the actual subdomains have not been supplied.
326: We create the requisite number of local inner subdomains and then expand them into
327: out subdomains, if necessary.
328: */
329: PCGASMCreateLocalSubdomains(pc->pmat,osm->n,&osm->iis);
330: }
331: if (!osm->ois) {
332: /*
333: Initially make outer subdomains the same as inner subdomains. If nonzero additional overlap
334: has been requested, copy the inner subdomains over so they can be modified.
335: */
336: PetscMalloc1(osm->n,&osm->ois);
337: for (i=0; i<osm->n; ++i) {
338: if (osm->overlap > 0) { /* With positive overlap, osm->iis[i] will be modified */
339: ISDuplicate(osm->iis[i],(osm->ois)+i);
340: ISCopy(osm->iis[i],osm->ois[i]);
341: } else {
342: PetscObjectReference((PetscObject)((osm->iis)[i]));
343: osm->ois[i] = osm->iis[i];
344: }
345: }
346: if (osm->overlap > 0) {
347: /* Extend the "overlapping" regions by a number of steps */
348: MatIncreaseOverlap(pc->pmat,osm->n,osm->ois,osm->overlap);
349: }
350: }
352: /* Now the subdomains are defined. Determine their global and max local numbers, if necessary. */
353: if (osm->nmax == PETSC_DETERMINE) {
354: PetscMPIInt inwork,outwork;
355: /* determine global number of subdomains and the max number of local subdomains */
356: inwork = osm->n;
357: MPI_Allreduce(&inwork,&outwork,1,MPI_INT,MPI_MAX,PetscObjectComm((PetscObject)pc));
358: osm->nmax = outwork;
359: }
360: if (osm->N == PETSC_DETERMINE) {
361: /* Determine the number of globally-distinct subdomains and compute a global numbering for them. */
362: PetscObjectsListGetGlobalNumbering(PetscObjectComm((PetscObject)pc),osm->n,(PetscObject*)osm->ois,&osm->N,NULL);
363: }
366: if (osm->sort_indices) {
367: for (i=0; i<osm->n; i++) {
368: ISSort(osm->ois[i]);
369: ISSort(osm->iis[i]);
370: }
371: }
373: PCGetOptionsPrefix(pc,&prefix);
374: PCGASMPrintSubdomains(pc);
376: /*
377: Merge the ISs, create merged vectors and restrictions.
378: */
379: /* Merge outer subdomain ISs and construct a restriction onto the disjoint union of local outer subdomains. */
380: on = 0;
381: for (i=0; i<osm->n; i++) {
382: ISGetLocalSize(osm->ois[i],&oni);
383: on += oni;
384: }
385: PetscMalloc1(on, &oidx);
386: on = 0;
387: for (i=0; i<osm->n; i++) {
388: ISGetLocalSize(osm->ois[i],&oni);
389: ISGetIndices(osm->ois[i],&oidxi);
390: PetscMemcpy(oidx+on, oidxi, sizeof(PetscInt)*oni);
391: ISRestoreIndices(osm->ois[i], &oidxi);
392: on += oni;
393: }
394: ISCreateGeneral(((PetscObject)(pc))->comm, on, oidx, PETSC_OWN_POINTER, &gois);
395: MatCreateVecs(pc->pmat,&x,&y);
396: VecCreateMPI(PetscObjectComm((PetscObject)pc), on, PETSC_DECIDE, &osm->gx);
397: VecDuplicate(osm->gx,&osm->gy);
398: VecGetOwnershipRange(osm->gx, &gofirst, NULL);
399: ISCreateStride(PetscObjectComm((PetscObject)pc),on,gofirst,1, &goid);
400: VecScatterCreate(x,gois,osm->gx,goid, &(osm->gorestriction));
401: VecDestroy(&x);
402: ISDestroy(&gois);
404: /* Merge inner subdomain ISs and construct a restriction onto the disjoint union of local inner subdomains. */
405: {
406: PetscInt ini; /* Number of indices the i-th a local inner subdomain. */
407: PetscInt in; /* Number of indices in the disjoint uniont of local inner subdomains. */
408: PetscInt *iidx; /* Global indices in the merged local inner subdomain. */
409: PetscInt *ioidx; /* Global indices of the disjoint union of inner subdomains within the disjoint union of outer subdomains. */
410: IS giis; /* IS for the disjoint union of inner subdomains. */
411: IS giois; /* IS for the disjoint union of inner subdomains within the disjoint union of outer subdomains. */
412: /**/
413: in = 0;
414: for (i=0; i<osm->n; i++) {
415: ISGetLocalSize(osm->iis[i],&ini);
416: in += ini;
417: }
418: PetscMalloc1(in, &iidx);
419: PetscMalloc1(in, &ioidx);
420: VecGetOwnershipRange(osm->gx,&gofirst, NULL);
421: in = 0;
422: on = 0;
423: for (i=0; i<osm->n; i++) {
424: const PetscInt *iidxi; /* Global indices of the i-th local inner subdomain. */
425: ISLocalToGlobalMapping ltogi; /* Map from global to local indices of the i-th outer local subdomain. */
426: PetscInt *ioidxi; /* Local indices of the i-th local inner subdomain within the local outer subdomain. */
427: PetscInt ioni; /* Number of indices in ioidxi; if ioni != ini the inner subdomain is not a subdomain of the outer subdomain (error). */
428: PetscInt k;
429: ISGetLocalSize(osm->iis[i],&ini);
430: ISGetLocalSize(osm->ois[i],&oni);
431: ISGetIndices(osm->iis[i],&iidxi);
432: PetscMemcpy(iidx+in, iidxi, sizeof(PetscInt)*ini);
433: ISLocalToGlobalMappingCreateIS(osm->ois[i],<ogi);
434: ioidxi = ioidx+in;
435: ISGlobalToLocalMappingApply(ltogi,IS_GTOLM_DROP,ini,iidxi,&ioni,ioidxi);
436: ISLocalToGlobalMappingDestroy(<ogi);
437: ISRestoreIndices(osm->iis[i], &iidxi);
438: if (ioni != ini) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"Inner subdomain %D contains %D indices outside of its outer subdomain", i, ini - ioni);
439: for (k = 0; k < ini; ++k) ioidxi[k] += gofirst+on;
440: in += ini;
441: on += oni;
442: }
443: ISCreateGeneral(PetscObjectComm((PetscObject)pc), in, iidx, PETSC_OWN_POINTER, &giis);
444: ISCreateGeneral(PetscObjectComm((PetscObject)pc), in, ioidx, PETSC_OWN_POINTER, &giois);
445: VecScatterCreate(y,giis,osm->gy,giois,&osm->girestriction);
446: VecDestroy(&y);
447: ISDestroy(&giis);
448: ISDestroy(&giois);
449: }
450: ISDestroy(&goid);
452: /* Create the subdomain work vectors. */
453: PetscMalloc1(osm->n,&osm->x);
454: PetscMalloc1(osm->n,&osm->y);
455: VecGetArray(osm->gx, &gxarray);
456: VecGetArray(osm->gy, &gyarray);
457: for (i=0, on=0; i<osm->n; ++i, on += oni) {
458: PetscInt oNi;
459: ISGetLocalSize(osm->ois[i],&oni);
460: ISGetSize(osm->ois[i],&oNi);
461: VecCreateMPIWithArray(((PetscObject)(osm->ois[i]))->comm,1,oni,oNi,gxarray+on,&osm->x[i]);
462: VecCreateMPIWithArray(((PetscObject)(osm->ois[i]))->comm,1,oni,oNi,gyarray+on,&osm->y[i]);
463: }
464: VecRestoreArray(osm->gx, &gxarray);
465: VecRestoreArray(osm->gy, &gyarray);
466: /* Create the subdomain solvers */
467: PetscMalloc1(osm->n,&osm->ksp);
468: for (i=0; i<osm->n; i++) {
469: char subprefix[PETSC_MAX_PATH_LEN+1];
470: KSPCreate(((PetscObject)(osm->ois[i]))->comm,&ksp);
471: PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
472: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
473: KSPSetType(ksp,KSPPREONLY);
474: KSPGetPC(ksp,&subpc); /* Why do we need this here? */
475: if (subdomain_dm) {
476: KSPSetDM(ksp,subdomain_dm[i]);
477: DMDestroy(subdomain_dm+i);
478: }
479: PCGetOptionsPrefix(pc,&prefix);
480: KSPSetOptionsPrefix(ksp,prefix);
481: if (subdomain_names && subdomain_names[i]) {
482: PetscSNPrintf(subprefix,PETSC_MAX_PATH_LEN,"sub_%s_",subdomain_names[i]);
483: KSPAppendOptionsPrefix(ksp,subprefix);
484: PetscFree(subdomain_names[i]);
485: }
486: KSPAppendOptionsPrefix(ksp,"sub_");
487: osm->ksp[i] = ksp;
488: }
489: PetscFree(subdomain_dm);
490: PetscFree(subdomain_names);
491: scall = MAT_INITIAL_MATRIX;
493: } else { /*if (pc->setupcalled)*/
494: /*
495: Destroy the submatrices from the previous iteration
496: */
497: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
498: MatDestroyMatrices(osm->n,&osm->pmat);
499: scall = MAT_INITIAL_MATRIX;
500: }
501: }
503: /*
504: Extract out the submatrices.
505: */
506: if (size > 1) {
507: MatGetSubMatricesMPI(pc->pmat,osm->n,osm->ois,osm->ois,scall,&osm->pmat);
508: } else {
509: MatGetSubMatrices(pc->pmat,osm->n,osm->ois,osm->ois,scall,&osm->pmat);
510: }
511: if (scall == MAT_INITIAL_MATRIX) {
512: PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
513: for (i=0; i<osm->n; i++) {
514: PetscLogObjectParent((PetscObject)pc,(PetscObject)osm->pmat[i]);
515: PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);
516: }
517: }
519: /* Return control to the user so that the submatrices can be modified (e.g., to apply
520: different boundary conditions for the submatrices than for the global problem) */
521: PCModifySubMatrices(pc,osm->n,osm->ois,osm->ois,osm->pmat,pc->modifysubmatricesP);
523: /*
524: Loop over submatrices putting them into local ksps
525: */
526: for (i=0; i<osm->n; i++) {
527: KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i]);
528: if (!pc->setupcalled) {
529: KSPSetFromOptions(osm->ksp[i]);
530: }
531: }
532: return(0);
533: }
537: static PetscErrorCode PCSetUpOnBlocks_GASM(PC pc)
538: {
539: PC_GASM *osm = (PC_GASM*)pc->data;
541: PetscInt i;
544: for (i=0; i<osm->n; i++) {
545: KSPSetUp(osm->ksp[i]);
546: }
547: return(0);
548: }
552: static PetscErrorCode PCApply_GASM(PC pc,Vec x,Vec y)
553: {
554: PC_GASM *osm = (PC_GASM*)pc->data;
556: PetscInt i;
557: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
560: /*
561: Support for limiting the restriction or interpolation only to the inner
562: subdomain values (leaving the other values 0).
563: */
564: if (!(osm->type & PC_GASM_RESTRICT)) {
565: /* have to zero the work RHS since scatter may leave some slots empty */
566: VecZeroEntries(osm->gx);
567: VecScatterBegin(osm->girestriction,x,osm->gx,INSERT_VALUES,forward);
568: } else {
569: VecScatterBegin(osm->gorestriction,x,osm->gx,INSERT_VALUES,forward);
570: }
571: VecZeroEntries(osm->gy);
572: if (!(osm->type & PC_GASM_RESTRICT)) {
573: VecScatterEnd(osm->girestriction,x,osm->gx,INSERT_VALUES,forward);
574: } else {
575: VecScatterEnd(osm->gorestriction,x,osm->gx,INSERT_VALUES,forward);
576: }
577: /* do the subdomain solves */
578: for (i=0; i<osm->n; ++i) {
579: KSPSolve(osm->ksp[i],osm->x[i],osm->y[i]);
580: }
581: VecZeroEntries(y);
582: if (!(osm->type & PC_GASM_INTERPOLATE)) {
583: VecScatterBegin(osm->girestriction,osm->gy,y,ADD_VALUES,reverse);
584: VecScatterEnd(osm->girestriction,osm->gy,y,ADD_VALUES,reverse); return(0);
585: } else {
586: VecScatterBegin(osm->gorestriction,osm->gy,y,ADD_VALUES,reverse);
587: VecScatterEnd(osm->gorestriction,osm->gy,y,ADD_VALUES,reverse); return(0);
588: }
589: }
593: static PetscErrorCode PCApplyTranspose_GASM(PC pc,Vec x,Vec y)
594: {
595: PC_GASM *osm = (PC_GASM*)pc->data;
597: PetscInt i;
598: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
601: /*
602: Support for limiting the restriction or interpolation to only local
603: subdomain values (leaving the other values 0).
605: Note: these are reversed from the PCApply_GASM() because we are applying the
606: transpose of the three terms
607: */
608: if (!(osm->type & PC_GASM_INTERPOLATE)) {
609: /* have to zero the work RHS since scatter may leave some slots empty */
610: VecZeroEntries(osm->gx);
611: VecScatterBegin(osm->girestriction,x,osm->gx,INSERT_VALUES,forward);
612: } else {
613: VecScatterBegin(osm->gorestriction,x,osm->gx,INSERT_VALUES,forward);
614: }
615: VecZeroEntries(osm->gy);
616: if (!(osm->type & PC_GASM_INTERPOLATE)) {
617: VecScatterEnd(osm->girestriction,x,osm->gx,INSERT_VALUES,forward);
618: } else {
619: VecScatterEnd(osm->gorestriction,x,osm->gx,INSERT_VALUES,forward);
620: }
621: /* do the local solves */
622: for (i=0; i<osm->n; ++i) { /* Note that the solves are local, so we can go to osm->n, rather than osm->nmax. */
623: KSPSolveTranspose(osm->ksp[i],osm->x[i],osm->y[i]);
624: }
625: VecZeroEntries(y);
626: if (!(osm->type & PC_GASM_RESTRICT)) {
627: VecScatterBegin(osm->girestriction,osm->gy,y,ADD_VALUES,reverse);
628: VecScatterEnd(osm->girestriction,osm->gy,y,ADD_VALUES,reverse);
629: } else {
630: VecScatterBegin(osm->gorestriction,osm->gy,y,ADD_VALUES,reverse);
631: VecScatterEnd(osm->gorestriction,osm->gy,y,ADD_VALUES,reverse);
632: }
633: return(0);
634: }
638: static PetscErrorCode PCReset_GASM(PC pc)
639: {
640: PC_GASM *osm = (PC_GASM*)pc->data;
642: PetscInt i;
645: if (osm->ksp) {
646: for (i=0; i<osm->n; i++) {
647: KSPReset(osm->ksp[i]);
648: }
649: }
650: if (osm->pmat) {
651: if (osm->n > 0) {
652: MatDestroyMatrices(osm->n,&osm->pmat);
653: }
654: }
655: if (osm->x) {
656: for (i=0; i<osm->n; i++) {
657: VecDestroy(&osm->x[i]);
658: VecDestroy(&osm->y[i]);
659: }
660: }
661: VecDestroy(&osm->gx);
662: VecDestroy(&osm->gy);
664: VecScatterDestroy(&osm->gorestriction);
665: VecScatterDestroy(&osm->girestriction);
666: if (!osm->user_subdomains) {
667: PCGASMDestroySubdomains(osm->n,&osm->ois,&osm->iis);
668: osm->N = PETSC_DETERMINE;
669: osm->nmax = PETSC_DETERMINE;
670: }
671: return(0);
672: }
676: static PetscErrorCode PCDestroy_GASM(PC pc)
677: {
678: PC_GASM *osm = (PC_GASM*)pc->data;
680: PetscInt i;
683: PCReset_GASM(pc);
685: /* PCReset will not destroy subdomains, if user_subdomains is true. */
686: PCGASMDestroySubdomains(osm->n,&osm->ois,&osm->iis);
688: if (osm->ksp) {
689: for (i=0; i<osm->n; i++) {
690: KSPDestroy(&osm->ksp[i]);
691: }
692: PetscFree(osm->ksp);
693: }
694: PetscFree(osm->x);
695: PetscFree(osm->y);
696: PetscFree(pc->data);
697: return(0);
698: }
702: static PetscErrorCode PCSetFromOptions_GASM(PetscOptions *PetscOptionsObject,PC pc)
703: {
704: PC_GASM *osm = (PC_GASM*)pc->data;
706: PetscInt blocks,ovl;
707: PetscBool symset,flg;
708: PCGASMType gasmtype;
711: /* set the type to symmetric if matrix is symmetric */
712: if (!osm->type_set && pc->pmat) {
713: MatIsSymmetricKnown(pc->pmat,&symset,&flg);
714: if (symset && flg) osm->type = PC_GASM_BASIC;
715: }
716: PetscOptionsHead(PetscOptionsObject,"Generalized additive Schwarz options");
717: PetscOptionsBool("-pc_gasm_use_dm_subdomains","If subdomains aren't set, use DMCreateDomainDecomposition() to define subdomains.","PCGASMSetUseDMSubdomains",osm->dm_subdomains,&osm->dm_subdomains,&flg);
718: PetscOptionsInt("-pc_gasm_total_subdomains","Total number of subdomains across communicator","PCGASMSetTotalSubdomains",osm->N,&blocks,&flg);
719: if (flg) {
720: PCGASMSetTotalSubdomains(pc,blocks);
721: }
722: PetscOptionsInt("-pc_gasm_overlap","Number of overlapping degrees of freedom","PCGASMSetOverlap",osm->overlap,&ovl,&flg);
723: if (flg) {
724: PCGASMSetOverlap(pc,ovl);
725: osm->dm_subdomains = PETSC_FALSE;
726: }
727: flg = PETSC_FALSE;
728: PetscOptionsEnum("-pc_gasm_type","Type of restriction/extension","PCGASMSetType",PCGASMTypes,(PetscEnum)osm->type,(PetscEnum*)&gasmtype,&flg);
729: if (flg) {PCGASMSetType(pc,gasmtype);}
730: PetscOptionsTail();
731: return(0);
732: }
734: /*------------------------------------------------------------------------------------*/
738: /*@
739: PCGASMSetTotalSubdomains - sets the total number of subdomains to use across the
740: communicator.
741: Logically collective on pc
743: Input Parameters:
744: + pc - the preconditioner
745: - N - total number of subdomains
748: Level: beginner
750: .keywords: PC, ASM, DM, set, subdomains, additive Schwarz
752: .seealso: PCGASMSetSubdomains(), PCGASMSetOverlap()
753: PCGASMCreateSubdomains2D()
754: @*/
755: PetscErrorCode PCGASMSetTotalSubdomains(PC pc,PetscInt N)
756: {
757: PC_GASM *osm = (PC_GASM*)pc->data;
758: PetscMPIInt size,rank;
762: if (N < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Total number of subdomains must be 1 or more, got N = %D",N);
763: if (pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCGASMSetTotalSubdomains() should be called before calling PCSetUp().");
765: PCGASMDestroySubdomains(osm->n,&osm->iis,&osm->ois);
766: osm->ois = osm->iis = NULL;
768: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
769: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
770: osm->N = N;
771: osm->n = PETSC_DETERMINE;
772: osm->nmax = PETSC_DETERMINE;
773: osm->dm_subdomains = PETSC_FALSE;
774: return(0);
775: }
779: static PetscErrorCode PCGASMSetSubdomains_GASM(PC pc,PetscInt n,IS iis[],IS ois[])
780: {
781: PC_GASM *osm = (PC_GASM*)pc->data;
783: PetscInt i;
786: if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Each process must have 1 or more subdomains, got n = %D",n);
787: if (pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCGASMSetSubdomains() should be called before calling PCSetUp().");
789: PCGASMDestroySubdomains(osm->n,&osm->iis,&osm->ois);
790: osm->iis = osm->ois = NULL;
791: osm->n = n;
792: osm->N = PETSC_DETERMINE;
793: osm->nmax = PETSC_DETERMINE;
794: if (ois) {
795: PetscMalloc1(n,&osm->ois);
796: for (i=0; i<n; i++) {
797: PetscObjectReference((PetscObject)ois[i]);
798: osm->ois[i] = ois[i];
799: }
800: /*
801: Since the user set the outer subdomains, even if nontrivial overlap was requested via PCGASMSetOverlap(),
802: it will be ignored. To avoid confusion later on (e.g., when viewing the PC), the overlap size is set to -1.
803: */
804: osm->overlap = -1;
805: if (!iis) {
806: PetscMalloc1(n,&osm->iis);
807: for (i=0; i<n; i++) {
808: for (i=0; i<n; i++) {
809: PetscObjectReference((PetscObject)ois[i]);
810: osm->iis[i] = ois[i];
811: }
812: }
813: }
814: }
815: if (iis) {
816: PetscMalloc1(n,&osm->iis);
817: for (i=0; i<n; i++) {
818: PetscObjectReference((PetscObject)iis[i]);
819: osm->iis[i] = iis[i];
820: }
821: if (!ois) {
822: PetscMalloc1(n,&osm->ois);
823: for (i=0; i<n; i++) {
824: for (i=0; i<n; i++) {
825: PetscObjectReference((PetscObject)iis[i]);
826: osm->ois[i] = iis[i];
827: }
828: }
829: /* If necessary, osm->ois[i] will be expanded using the requested overlap size in PCSetUp_GASM(). */
830: }
831: }
832: if (iis) osm->user_subdomains = PETSC_TRUE;
833: return(0);
834: }
839: static PetscErrorCode PCGASMSetOverlap_GASM(PC pc,PetscInt ovl)
840: {
841: PC_GASM *osm = (PC_GASM*)pc->data;
844: if (ovl < 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested");
845: if (pc->setupcalled && ovl != osm->overlap) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCGASMSetOverlap() should be called before PCSetUp().");
846: if (!pc->setupcalled) osm->overlap = ovl;
847: return(0);
848: }
852: static PetscErrorCode PCGASMSetType_GASM(PC pc,PCGASMType type)
853: {
854: PC_GASM *osm = (PC_GASM*)pc->data;
857: osm->type = type;
858: osm->type_set = PETSC_TRUE;
859: return(0);
860: }
864: static PetscErrorCode PCGASMSetSortIndices_GASM(PC pc,PetscBool doSort)
865: {
866: PC_GASM *osm = (PC_GASM*)pc->data;
869: osm->sort_indices = doSort;
870: return(0);
871: }
875: /*
876: FIXME: This routine might need to be modified now that multiple ranks per subdomain are allowed.
877: In particular, it would upset the global subdomain number calculation.
878: */
879: static PetscErrorCode PCGASMGetSubKSP_GASM(PC pc,PetscInt *n,PetscInt *first,KSP **ksp)
880: {
881: PC_GASM *osm = (PC_GASM*)pc->data;
885: if (osm->n < 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Need to call PCSetUP() on PC (or KSPSetUp() on the outer KSP object) before calling here");
887: if (n) *n = osm->n;
888: if (first) {
889: MPI_Scan(&osm->n,first,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
890: *first -= osm->n;
891: }
892: if (ksp) {
893: /* Assume that local solves are now different; not necessarily
894: true, though! This flag is used only for PCView_GASM() */
895: *ksp = osm->ksp;
896: osm->same_subdomain_solvers = PETSC_FALSE;
897: }
898: return(0);
899: } /* PCGASMGetSubKSP_GASM() */
903: /*@C
904: PCGASMSetSubdomains - Sets the subdomains for this processor
905: for the additive Schwarz preconditioner.
907: Collective on pc
909: Input Parameters:
910: + pc - the preconditioner object
911: . n - the number of subdomains for this processor
912: . iis - the index sets that define the inner subdomains (or NULL for PETSc to determine subdomains)
913: - ois - the index sets that define the outer subdomains (or NULL to use the same as iis, or to construct by expanding iis by the requested overlap)
915: Notes:
916: The IS indices use the parallel, global numbering of the vector entries.
917: Inner subdomains are those where the correction is applied.
918: Outer subdomains are those where the residual necessary to obtain the
919: corrections is obtained (see PCGASMType for the use of inner/outer subdomains).
920: Both inner and outer subdomains can extend over several processors.
921: This processor's portion of a subdomain is known as a local subdomain.
923: By default the GASM preconditioner uses 1 (local) subdomain per processor.
926: Level: advanced
928: .keywords: PC, GASM, set, subdomains, additive Schwarz
930: .seealso: PCGASMSetNumSubdomains(), PCGASMSetOverlap(), PCGASMGetSubKSP(),
931: PCGASMCreateSubdomains2D(), PCGASMGetSubdomains()
932: @*/
933: PetscErrorCode PCGASMSetSubdomains(PC pc,PetscInt n,IS iis[],IS ois[])
934: {
935: PC_GASM *osm = (PC_GASM*)pc->data;
940: PetscTryMethod(pc,"PCGASMSetSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,iis,ois));
941: osm->dm_subdomains = PETSC_FALSE;
942: return(0);
943: }
948: /*@
949: PCGASMSetOverlap - Sets the overlap between a pair of subdomains for the
950: additive Schwarz preconditioner. Either all or no processors in the
951: pc communicator must call this routine.
953: Logically Collective on pc
955: Input Parameters:
956: + pc - the preconditioner context
957: - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 0)
959: Options Database Key:
960: . -pc_gasm_overlap <overlap> - Sets overlap
962: Notes:
963: By default the GASM preconditioner uses 1 subdomain per processor. To use
964: multiple subdomain per perocessor or "straddling" subdomains that intersect
965: multiple processors use PCGASMSetSubdomains() (or option -pc_gasm_total_subdomains <n>).
967: The overlap defaults to 0, so if one desires that no additional
968: overlap be computed beyond what may have been set with a call to
969: PCGASMSetSubdomains(), then ovl must be set to be 0. In particular, if one does
970: not explicitly set the subdomains in application code, then all overlap would be computed
971: internally by PETSc, and using an overlap of 0 would result in an GASM
972: variant that is equivalent to the block Jacobi preconditioner.
974: Note that one can define initial index sets with any overlap via
975: PCGASMSetSubdomains(); the routine PCGASMSetOverlap() merely allows
976: PETSc to extend that overlap further, if desired.
978: Level: intermediate
980: .keywords: PC, GASM, set, overlap
982: .seealso: PCGASMSetSubdomains(), PCGASMGetSubKSP(),
983: PCGASMCreateSubdomains2D(), PCGASMGetSubdomains()
984: @*/
985: PetscErrorCode PCGASMSetOverlap(PC pc,PetscInt ovl)
986: {
988: PC_GASM *osm = (PC_GASM*)pc->data;
993: PetscTryMethod(pc,"PCGASMSetOverlap_C",(PC,PetscInt),(pc,ovl));
994: osm->dm_subdomains = PETSC_FALSE;
995: return(0);
996: }
1000: /*@
1001: PCGASMSetType - Sets the type of restriction and interpolation used
1002: for local problems in the additive Schwarz method.
1004: Logically Collective on PC
1006: Input Parameters:
1007: + pc - the preconditioner context
1008: - type - variant of GASM, one of
1009: .vb
1010: PC_GASM_BASIC - full interpolation and restriction
1011: PC_GASM_RESTRICT - full restriction, local processor interpolation
1012: PC_GASM_INTERPOLATE - full interpolation, local processor restriction
1013: PC_GASM_NONE - local processor restriction and interpolation
1014: .ve
1016: Options Database Key:
1017: . -pc_gasm_type [basic,restrict,interpolate,none] - Sets GASM type
1019: Level: intermediate
1021: .keywords: PC, GASM, set, type
1023: .seealso: PCGASMSetSubdomains(), PCGASMGetSubKSP(),
1024: PCGASMCreateSubdomains2D()
1025: @*/
1026: PetscErrorCode PCGASMSetType(PC pc,PCGASMType type)
1027: {
1033: PetscTryMethod(pc,"PCGASMSetType_C",(PC,PCGASMType),(pc,type));
1034: return(0);
1035: }
1039: /*@
1040: PCGASMSetSortIndices - Determines whether subdomain indices are sorted.
1042: Logically Collective on PC
1044: Input Parameters:
1045: + pc - the preconditioner context
1046: - doSort - sort the subdomain indices
1048: Level: intermediate
1050: .keywords: PC, GASM, set, type
1052: .seealso: PCGASMSetSubdomains(), PCGASMGetSubKSP(),
1053: PCGASMCreateSubdomains2D()
1054: @*/
1055: PetscErrorCode PCGASMSetSortIndices(PC pc,PetscBool doSort)
1056: {
1062: PetscTryMethod(pc,"PCGASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));
1063: return(0);
1064: }
1068: /*@C
1069: PCGASMGetSubKSP - Gets the local KSP contexts for all blocks on
1070: this processor.
1072: Collective on PC iff first_local is requested
1074: Input Parameter:
1075: . pc - the preconditioner context
1077: Output Parameters:
1078: + n_local - the number of blocks on this processor or NULL
1079: . first_local - the global number of the first block on this processor or NULL,
1080: all processors must request or all must pass NULL
1081: - ksp - the array of KSP contexts
1083: Note:
1084: After PCGASMGetSubKSP() the array of KSPes is not to be freed
1086: Currently for some matrix implementations only 1 block per processor
1087: is supported.
1089: You must call KSPSetUp() before calling PCGASMGetSubKSP().
1091: Level: advanced
1093: .keywords: PC, GASM, additive Schwarz, get, sub, KSP, context
1095: .seealso: PCGASMSetSubdomains(), PCGASMSetOverlap(),
1096: PCGASMCreateSubdomains2D(),
1097: @*/
1098: PetscErrorCode PCGASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
1099: {
1104: PetscUseMethod(pc,"PCGASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
1105: return(0);
1106: }
1108: /* -------------------------------------------------------------------------------------*/
1109: /*MC
1110: PCGASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1111: its own KSP object.
1113: Options Database Keys:
1114: + -pc_gasm_total_subdomains <n> - Sets total number of local subdomains to be distributed among processors
1115: . -pc_gasm_view_subdomains - activates the printing of subdomain indices in PCView(), -ksp_view or -snes_view
1116: . -pc_gasm_print_subdomains - activates the printing of subdomain indices in PCSetUp()
1117: . -pc_gasm_overlap <ovl> - Sets overlap by which to (automatically) extend local subdomains
1118: - -pc_gasm_type [basic,restrict,interpolate,none] - Sets GASM type
1120: IMPORTANT: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you
1121: will get a different convergence rate due to the default option of -pc_gasm_type restrict. Use
1122: -pc_gasm_type basic to use the standard GASM.
1124: Notes: Each processor can have one or more blocks, but a block cannot be shared by more
1125: than one processor. Defaults to one block per processor.
1127: To set options on the solvers for each block append -sub_ to all the KSP, and PC
1128: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
1130: To set the options on the solvers separate for each block call PCGASMGetSubKSP()
1131: and set the options directly on the resulting KSP object (you can access its PC
1132: with KSPGetPC())
1135: Level: beginner
1137: Concepts: additive Schwarz method
1139: References:
1140: An additive variant of the Schwarz alternating method for the case of many subregions
1141: M Dryja, OB Widlund - Courant Institute, New York University Technical report
1143: Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations,
1144: Barry Smith, Petter Bjorstad, and William Gropp, Cambridge University Press, ISBN 0-521-49589-X.
1146: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
1147: PCBJACOBI, PCGASMGetSubKSP(), PCGASMSetSubdomains(),
1148: PCSetModifySubmatrices(), PCGASMSetOverlap(), PCGASMSetType()
1150: M*/
1154: PETSC_EXTERN PetscErrorCode PCCreate_GASM(PC pc)
1155: {
1157: PC_GASM *osm;
1160: PetscNewLog(pc,&osm);
1162: osm->N = PETSC_DETERMINE;
1163: osm->n = PETSC_DECIDE;
1164: osm->nmax = PETSC_DETERMINE;
1165: osm->overlap = 0;
1166: osm->ksp = 0;
1167: osm->gorestriction = 0;
1168: osm->girestriction = 0;
1169: osm->gx = 0;
1170: osm->gy = 0;
1171: osm->x = 0;
1172: osm->y = 0;
1173: osm->ois = 0;
1174: osm->iis = 0;
1175: osm->pmat = 0;
1176: osm->type = PC_GASM_RESTRICT;
1177: osm->same_subdomain_solvers = PETSC_TRUE;
1178: osm->sort_indices = PETSC_TRUE;
1179: osm->dm_subdomains = PETSC_FALSE;
1181: pc->data = (void*)osm;
1182: pc->ops->apply = PCApply_GASM;
1183: pc->ops->applytranspose = PCApplyTranspose_GASM;
1184: pc->ops->setup = PCSetUp_GASM;
1185: pc->ops->reset = PCReset_GASM;
1186: pc->ops->destroy = PCDestroy_GASM;
1187: pc->ops->setfromoptions = PCSetFromOptions_GASM;
1188: pc->ops->setuponblocks = PCSetUpOnBlocks_GASM;
1189: pc->ops->view = PCView_GASM;
1190: pc->ops->applyrichardson = 0;
1192: PetscObjectComposeFunction((PetscObject)pc,"PCGASMSetSubdomains_C",PCGASMSetSubdomains_GASM);
1193: PetscObjectComposeFunction((PetscObject)pc,"PCGASMSetOverlap_C",PCGASMSetOverlap_GASM);
1194: PetscObjectComposeFunction((PetscObject)pc,"PCGASMSetType_C",PCGASMSetType_GASM);
1195: PetscObjectComposeFunction((PetscObject)pc,"PCGASMSetSortIndices_C",PCGASMSetSortIndices_GASM);
1196: PetscObjectComposeFunction((PetscObject)pc,"PCGASMGetSubKSP_C",PCGASMGetSubKSP_GASM);
1197: return(0);
1198: }
1203: PetscErrorCode PCGASMCreateLocalSubdomains(Mat A, PetscInt nloc, IS *iis[])
1204: {
1205: MatPartitioning mpart;
1206: const char *prefix;
1207: PetscErrorCode (*f)(Mat,MatReuse,Mat*);
1208: PetscMPIInt size;
1209: PetscInt i,j,rstart,rend,bs;
1210: PetscBool isbaij = PETSC_FALSE,foundpart = PETSC_FALSE;
1211: Mat Ad = NULL, adj;
1212: IS ispart,isnumb,*is;
1213: PetscErrorCode ierr;
1216: if (nloc < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of local subdomains must > 0, got nloc = %D",nloc);
1218: /* Get prefix, row distribution, and block size */
1219: MatGetOptionsPrefix(A,&prefix);
1220: MatGetOwnershipRange(A,&rstart,&rend);
1221: MatGetBlockSize(A,&bs);
1222: 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);
1224: /* Get diagonal block from matrix if possible */
1225: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1226: PetscObjectQueryFunction((PetscObject)A,"MatGetDiagonalBlock_C",&f);
1227: if (f) {
1228: MatGetDiagonalBlock(A,&Ad);
1229: } else if (size == 1) {
1230: Ad = A;
1231: }
1232: if (Ad) {
1233: PetscObjectTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij);
1234: if (!isbaij) {PetscObjectTypeCompare((PetscObject)Ad,MATSEQSBAIJ,&isbaij);}
1235: }
1236: if (Ad && nloc > 1) {
1237: PetscBool match,done;
1238: /* Try to setup a good matrix partitioning if available */
1239: MatPartitioningCreate(PETSC_COMM_SELF,&mpart);
1240: PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);
1241: MatPartitioningSetFromOptions(mpart);
1242: PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGCURRENT,&match);
1243: if (!match) {
1244: PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGSQUARE,&match);
1245: }
1246: if (!match) { /* assume a "good" partitioner is available */
1247: PetscInt na;
1248: const PetscInt *ia,*ja;
1249: MatGetRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1250: if (done) {
1251: /* Build adjacency matrix by hand. Unfortunately a call to
1252: MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1253: remove the block-aij structure and we cannot expect
1254: MatPartitioning to split vertices as we need */
1255: PetscInt i,j,len,nnz,cnt,*iia=0,*jja=0;
1256: const PetscInt *row;
1257: nnz = 0;
1258: for (i=0; i<na; i++) { /* count number of nonzeros */
1259: len = ia[i+1] - ia[i];
1260: row = ja + ia[i];
1261: for (j=0; j<len; j++) {
1262: if (row[j] == i) { /* don't count diagonal */
1263: len--; break;
1264: }
1265: }
1266: nnz += len;
1267: }
1268: PetscMalloc1(na+1,&iia);
1269: PetscMalloc1(nnz,&jja);
1270: nnz = 0;
1271: iia[0] = 0;
1272: for (i=0; i<na; i++) { /* fill adjacency */
1273: cnt = 0;
1274: len = ia[i+1] - ia[i];
1275: row = ja + ia[i];
1276: for (j=0; j<len; j++) {
1277: if (row[j] != i) jja[nnz+cnt++] = row[j]; /* if not diagonal */
1278: }
1279: nnz += cnt;
1280: iia[i+1] = nnz;
1281: }
1282: /* Partitioning of the adjacency matrix */
1283: MatCreateMPIAdj(PETSC_COMM_SELF,na,na,iia,jja,NULL,&adj);
1284: MatPartitioningSetAdjacency(mpart,adj);
1285: MatPartitioningSetNParts(mpart,nloc);
1286: MatPartitioningApply(mpart,&ispart);
1287: ISPartitioningToNumbering(ispart,&isnumb);
1288: MatDestroy(&adj);
1289: foundpart = PETSC_TRUE;
1290: }
1291: MatRestoreRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1292: }
1293: MatPartitioningDestroy(&mpart);
1294: }
1295: PetscMalloc1(nloc,&is);
1296: if (!foundpart) {
1298: /* Partitioning by contiguous chunks of rows */
1300: PetscInt mbs = (rend-rstart)/bs;
1301: PetscInt start = rstart;
1302: for (i=0; i<nloc; i++) {
1303: PetscInt count = (mbs/nloc + ((mbs % nloc) > i)) * bs;
1304: ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);
1305: start += count;
1306: }
1308: } else {
1310: /* Partitioning by adjacency of diagonal block */
1312: const PetscInt *numbering;
1313: PetscInt *count,nidx,*indices,*newidx,start=0;
1314: /* Get node count in each partition */
1315: PetscMalloc1(nloc,&count);
1316: ISPartitioningCount(ispart,nloc,count);
1317: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1318: for (i=0; i<nloc; i++) count[i] *= bs;
1319: }
1320: /* Build indices from node numbering */
1321: ISGetLocalSize(isnumb,&nidx);
1322: PetscMalloc1(nidx,&indices);
1323: for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */
1324: ISGetIndices(isnumb,&numbering);
1325: PetscSortIntWithPermutation(nidx,numbering,indices);
1326: ISRestoreIndices(isnumb,&numbering);
1327: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1328: PetscMalloc1(nidx*bs,&newidx);
1329: for (i=0; i<nidx; i++) {
1330: for (j=0; j<bs; j++) newidx[i*bs+j] = indices[i]*bs + j;
1331: }
1332: PetscFree(indices);
1333: nidx *= bs;
1334: indices = newidx;
1335: }
1336: /* Shift to get global indices */
1337: for (i=0; i<nidx; i++) indices[i] += rstart;
1339: /* Build the index sets for each block */
1340: for (i=0; i<nloc; i++) {
1341: ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);
1342: ISSort(is[i]);
1343: start += count[i];
1344: }
1346: PetscFree(count);
1347: PetscFree(indices);
1348: ISDestroy(&isnumb);
1349: ISDestroy(&ispart);
1350: }
1351: *iis = is;
1352: return(0);
1353: }
1357: PETSC_INTERN PetscErrorCode PCGASMCreateStraddlingSubdomains(Mat A,PetscInt N,PetscInt *n,IS *iis[])
1358: {
1359: PetscErrorCode ierr;
1362: MatSubdomainsCreateCoalesce(A,N,n,iis);
1363: return(0);
1364: }
1370: /*@C
1371: PCGASMCreateSubdomains - Creates n index sets defining n nonoverlapping subdomains for the additive
1372: Schwarz preconditioner for a any problem based on its matrix.
1374: Collective
1376: Input Parameters:
1377: + A - The global matrix operator
1378: - N - the number of global subdomains requested
1380: Output Parameters:
1381: + n - the number of subdomains created on this processor
1382: - iis - the array of index sets defining the local inner subdomains (on which the correction is applied)
1384: Level: advanced
1386: Note: When N >= A's communicator size, each subdomain is local -- contained within a single processor.
1387: When N < size, the subdomains are 'straddling' (processor boundaries) and are no longer local.
1388: The resulting subdomains can be use in PCGASMSetSubdomains(pc,n,iss,NULL). The overlapping
1389: outer subdomains will be automatically generated from these according to the requested amount of
1390: overlap; this is currently supported only with local subdomains.
1393: .keywords: PC, GASM, additive Schwarz, create, subdomains, unstructured grid
1395: .seealso: PCGASMSetSubdomains(), PCGASMDestroySubdomains()
1396: @*/
1397: PetscErrorCode PCGASMCreateSubdomains(Mat A,PetscInt N,PetscInt *n,IS *iis[])
1398: {
1399: PetscMPIInt size;
1400: PetscErrorCode ierr;
1406: if (N < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Number of subdomains must be > 0, N = %D",N);
1407: MPI_Comm_size(PetscObjectComm((PetscObject)A),&size);
1408: if (N >= size) {
1409: *n = N/size + (N%size);
1410: PCGASMCreateLocalSubdomains(A,*n,iis);
1411: } else {
1412: PCGASMCreateStraddlingSubdomains(A,N,n,iis);
1413: }
1414: return(0);
1415: }
1419: /*@C
1420: PCGASMDestroySubdomains - Destroys the index sets created with
1421: PCGASMCreateSubdomains() or PCGASMCreateSubdomains2D. Should be
1422: called after setting subdomains with PCGASMSetSubdomains().
1424: Collective
1426: Input Parameters:
1427: + n - the number of index sets
1428: . iis - the array of inner subdomains,
1429: - ois - the array of outer subdomains, can be NULL
1431: Level: intermediate
1433: Notes: this is merely a convenience subroutine that walks each list,
1434: destroys each IS on the list, and then frees the list. At the end the
1435: list pointers are set to NULL.
1437: .keywords: PC, GASM, additive Schwarz, create, subdomains, unstructured grid
1439: .seealso: PCGASMCreateSubdomains(), PCGASMSetSubdomains()
1440: @*/
1441: PetscErrorCode PCGASMDestroySubdomains(PetscInt n,IS **iis,IS **ois)
1442: {
1443: PetscInt i;
1447: if (n <= 0) return(0);
1448: if (iis) {
1450: if (*iis) {
1452: for (i=0; i<n; i++) {
1453: ISDestroy(&(*iis)[i]);
1454: }
1455: PetscFree((*iis));
1456: }
1457: }
1458: if (ois) {
1460: if (*ois) {
1462: for (i=0; i<n; i++) {
1463: ISDestroy(&(*ois)[i]);
1464: }
1465: PetscFree((*ois));
1466: }
1467: }
1468: return(0);
1469: }
1472: #define PCGASMLocalSubdomainBounds2D(M,N,xleft,ylow,xright,yhigh,first,last,xleft_loc,ylow_loc,xright_loc,yhigh_loc,n) \
1473: { \
1474: PetscInt first_row = first/M, last_row = last/M+1; \
1475: /* \
1476: Compute ylow_loc and yhigh_loc so that (ylow_loc,xleft) and (yhigh_loc,xright) are the corners \
1477: of the bounding box of the intersection of the subdomain with the local ownership range (local \
1478: subdomain). \
1479: Also compute xleft_loc and xright_loc as the lower and upper bounds on the first and last rows \
1480: of the intersection. \
1481: */ \
1482: /* ylow_loc is the grid row containing the first element of the local sumbdomain */ \
1483: *ylow_loc = PetscMax(first_row,ylow); \
1484: /* xleft_loc is the offset of first element of the local subdomain within its grid row (might actually be outside the local subdomain) */ \
1485: *xleft_loc = *ylow_loc==first_row ? PetscMax(first%M,xleft) : xleft; \
1486: /* yhigh_loc is the grid row above the last local subdomain element */ \
1487: *yhigh_loc = PetscMin(last_row,yhigh); \
1488: /* xright is the offset of the end of the local subdomain within its grid row (might actually be outside the local subdomain) */ \
1489: *xright_loc = *yhigh_loc==last_row ? PetscMin(xright,last%M) : xright; \
1490: /* Now compute the size of the local subdomain n. */ \
1491: *n = 0; \
1492: if (*ylow_loc < *yhigh_loc) { \
1493: PetscInt width = xright-xleft; \
1494: *n += width*(*yhigh_loc-*ylow_loc-1); \
1495: *n += PetscMin(PetscMax(*xright_loc-xleft,0),width); \
1496: *n -= PetscMin(PetscMax(*xleft_loc-xleft,0), width); \
1497: } \
1498: }
1504: /*@
1505: PCGASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1506: preconditioner for a two-dimensional problem on a regular grid.
1508: Collective
1510: Input Parameters:
1511: + M, N - the global number of grid points in the x and y directions
1512: . Mdomains, Ndomains - the global number of subdomains in the x and y directions
1513: . dof - degrees of freedom per node
1514: - overlap - overlap in mesh lines
1516: Output Parameters:
1517: + Nsub - the number of local subdomains created
1518: . iis - array of index sets defining inner (nonoverlapping) subdomains
1519: - ois - array of index sets defining outer (overlapping, if overlap > 0) subdomains
1522: Level: advanced
1524: .keywords: PC, GASM, additive Schwarz, create, subdomains, 2D, regular grid
1526: .seealso: PCGASMSetSubdomains(), PCGASMGetSubKSP(), PCGASMSetOverlap()
1527: @*/
1528: PetscErrorCode PCGASMCreateSubdomains2D(PC pc,PetscInt M,PetscInt N,PetscInt Mdomains,PetscInt Ndomains,PetscInt dof,PetscInt overlap,PetscInt *nsub,IS **iis,IS **ois)
1529: {
1531: PetscMPIInt size, rank;
1532: PetscInt i, j;
1533: PetscInt maxheight, maxwidth;
1534: PetscInt xstart, xleft, xright, xleft_loc, xright_loc;
1535: PetscInt ystart, ylow, yhigh, ylow_loc, yhigh_loc;
1536: PetscInt x[2][2], y[2][2], n[2];
1537: PetscInt first, last;
1538: PetscInt nidx, *idx;
1539: PetscInt ii,jj,s,q,d;
1540: PetscInt k,kk;
1541: PetscMPIInt color;
1542: MPI_Comm comm, subcomm;
1543: IS **xis = 0, **is = ois, **is_local = iis;
1546: PetscObjectGetComm((PetscObject)pc, &comm);
1547: MPI_Comm_size(comm, &size);
1548: MPI_Comm_rank(comm, &rank);
1549: MatGetOwnershipRange(pc->pmat, &first, &last);
1550: 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) "
1551: "does not respect the number of degrees of freedom per grid point %D", first, last, dof);
1553: /* Determine the number of domains with nonzero intersections with the local ownership range. */
1554: s = 0;
1555: ystart = 0;
1556: for (j=0; j<Ndomains; ++j) {
1557: maxheight = N/Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */
1558: 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);
1559: /* Vertical domain limits with an overlap. */
1560: ylow = PetscMax(ystart - overlap,0);
1561: yhigh = PetscMin(ystart + maxheight + overlap,N);
1562: xstart = 0;
1563: for (i=0; i<Mdomains; ++i) {
1564: maxwidth = M/Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */
1565: 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);
1566: /* Horizontal domain limits with an overlap. */
1567: xleft = PetscMax(xstart - overlap,0);
1568: xright = PetscMin(xstart + maxwidth + overlap,M);
1569: /*
1570: Determine whether this subdomain intersects this processor's ownership range of pc->pmat.
1571: */
1572: PCGASMLocalSubdomainBounds2D(M,N,xleft,ylow,xright,yhigh,first,last,(&xleft_loc),(&ylow_loc),(&xright_loc),(&yhigh_loc),(&nidx));
1573: if (nidx) ++s;
1574: xstart += maxwidth;
1575: } /* for (i = 0; i < Mdomains; ++i) */
1576: ystart += maxheight;
1577: } /* for (j = 0; j < Ndomains; ++j) */
1579: /* Now we can allocate the necessary number of ISs. */
1580: *nsub = s;
1581: PetscMalloc1(*nsub,is);
1582: PetscMalloc1(*nsub,is_local);
1583: s = 0;
1584: ystart = 0;
1585: for (j=0; j<Ndomains; ++j) {
1586: maxheight = N/Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */
1587: 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);
1588: /* Vertical domain limits with an overlap. */
1589: y[0][0] = PetscMax(ystart - overlap,0);
1590: y[0][1] = PetscMin(ystart + maxheight + overlap,N);
1591: /* Vertical domain limits without an overlap. */
1592: y[1][0] = ystart;
1593: y[1][1] = PetscMin(ystart + maxheight,N);
1594: xstart = 0;
1595: for (i=0; i<Mdomains; ++i) {
1596: maxwidth = M/Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */
1597: 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);
1598: /* Horizontal domain limits with an overlap. */
1599: x[0][0] = PetscMax(xstart - overlap,0);
1600: x[0][1] = PetscMin(xstart + maxwidth + overlap,M);
1601: /* Horizontal domain limits without an overlap. */
1602: x[1][0] = xstart;
1603: x[1][1] = PetscMin(xstart+maxwidth,M);
1604: /*
1605: Determine whether this domain intersects this processor's ownership range of pc->pmat.
1606: Do this twice: first for the domains with overlaps, and once without.
1607: During the first pass create the subcommunicators, and use them on the second pass as well.
1608: */
1609: for (q = 0; q < 2; ++q) {
1610: /*
1611: domain limits, (xleft, xright) and (ylow, yheigh) are adjusted
1612: according to whether the domain with an overlap or without is considered.
1613: */
1614: xleft = x[q][0]; xright = x[q][1];
1615: ylow = y[q][0]; yhigh = y[q][1];
1616: PCGASMLocalSubdomainBounds2D(M,N,xleft,ylow,xright,yhigh,first,last,(&xleft_loc),(&ylow_loc),(&xright_loc),(&yhigh_loc),(&nidx));
1617: nidx *= dof;
1618: n[q] = nidx;
1619: /*
1620: Based on the counted number of indices in the local domain *with an overlap*,
1621: construct a subcommunicator of all the processors supporting this domain.
1622: Observe that a domain with an overlap might have nontrivial local support,
1623: while the domain without an overlap might not. Hence, the decision to participate
1624: in the subcommunicator must be based on the domain with an overlap.
1625: */
1626: if (q == 0) {
1627: if (nidx) color = 1;
1628: else color = MPI_UNDEFINED;
1630: MPI_Comm_split(comm, color, rank, &subcomm);
1631: }
1632: /*
1633: Proceed only if the number of local indices *with an overlap* is nonzero.
1634: */
1635: if (n[0]) {
1636: if (q == 0) xis = is;
1637: if (q == 1) {
1638: /*
1639: The IS for the no-overlap subdomain shares a communicator with the overlapping domain.
1640: Moreover, if the overlap is zero, the two ISs are identical.
1641: */
1642: if (overlap == 0) {
1643: (*is_local)[s] = (*is)[s];
1644: PetscObjectReference((PetscObject)(*is)[s]);
1645: continue;
1646: } else {
1647: xis = is_local;
1648: subcomm = ((PetscObject)(*is)[s])->comm;
1649: }
1650: } /* if (q == 1) */
1651: idx = NULL;
1652: PetscMalloc1(nidx,&idx);
1653: if (nidx) {
1654: k = 0;
1655: for (jj=ylow_loc; jj<yhigh_loc; ++jj) {
1656: PetscInt x0 = (jj==ylow_loc) ? xleft_loc : xleft;
1657: PetscInt x1 = (jj==yhigh_loc-1) ? xright_loc : xright;
1658: kk = dof*(M*jj + x0);
1659: for (ii=x0; ii<x1; ++ii) {
1660: for (d = 0; d < dof; ++d) {
1661: idx[k++] = kk++;
1662: }
1663: }
1664: }
1665: }
1666: ISCreateGeneral(subcomm,nidx,idx,PETSC_OWN_POINTER,(*xis)+s);
1667: }/* if (n[0]) */
1668: }/* for (q = 0; q < 2; ++q) */
1669: if (n[0]) ++s;
1670: xstart += maxwidth;
1671: } /* for (i = 0; i < Mdomains; ++i) */
1672: ystart += maxheight;
1673: } /* for (j = 0; j < Ndomains; ++j) */
1674: return(0);
1675: }
1679: /*@C
1680: PCGASMGetSubdomains - Gets the subdomains supported on this processor
1681: for the additive Schwarz preconditioner.
1683: Not Collective
1685: Input Parameter:
1686: . pc - the preconditioner context
1688: Output Parameters:
1689: + n - the number of subdomains for this processor (default value = 1)
1690: . iis - the index sets that define the inner subdomains (without overlap) supported on this processor (can be NULL)
1691: - ois - the index sets that define the outer subdomains (with overlap) supported on this processor (can be NULL)
1694: Notes:
1695: The user is responsible for destroying the ISs and freeing the returned arrays.
1696: The IS numbering is in the parallel, global numbering of the vector.
1698: Level: advanced
1700: .keywords: PC, GASM, get, subdomains, additive Schwarz
1702: .seealso: PCGASMSetOverlap(), PCGASMGetSubKSP(), PCGASMCreateSubdomains2D(),
1703: PCGASMSetSubdomains(), PCGASMGetSubmatrices()
1704: @*/
1705: PetscErrorCode PCGASMGetSubdomains(PC pc,PetscInt *n,IS *iis[],IS *ois[])
1706: {
1707: PC_GASM *osm;
1709: PetscBool match;
1710: PetscInt i;
1714: PetscObjectTypeCompare((PetscObject)pc,PCGASM,&match);
1715: if (!match) SETERRQ2(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Incorrect object type: expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name);
1716: osm = (PC_GASM*)pc->data;
1717: if (n) *n = osm->n;
1718: if (iis) {
1719: PetscMalloc1(osm->n, iis);
1720: }
1721: if (ois) {
1722: PetscMalloc1(osm->n, ois);
1723: }
1724: if (iis || ois) {
1725: for (i = 0; i < osm->n; ++i) {
1726: if (iis) (*iis)[i] = osm->iis[i];
1727: if (ois) (*ois)[i] = osm->ois[i];
1728: }
1729: }
1730: return(0);
1731: }
1735: /*@C
1736: PCGASMGetSubmatrices - Gets the local submatrices (for this processor
1737: only) for the additive Schwarz preconditioner.
1739: Not Collective
1741: Input Parameter:
1742: . pc - the preconditioner context
1744: Output Parameters:
1745: + n - the number of matrices for this processor (default value = 1)
1746: - mat - the matrices
1748: Notes: matrices returned by this routine have the same communicators as the index sets (IS)
1749: used to define subdomains in PCGASMSetSubdomains()
1750: Level: advanced
1752: .keywords: PC, GASM, set, local, subdomains, additive Schwarz, block Jacobi
1754: .seealso: PCGASMSetOverlap(), PCGASMGetSubKSP(),
1755: PCGASMCreateSubdomains2D(), PCGASMSetSubdomains(), PCGASMGetSubdomains()
1756: @*/
1757: PetscErrorCode PCGASMGetSubmatrices(PC pc,PetscInt *n,Mat *mat[])
1758: {
1759: PC_GASM *osm;
1761: PetscBool match;
1767: if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUP() or PCSetUp().");
1768: PetscObjectTypeCompare((PetscObject)pc,PCGASM,&match);
1769: if (!match) SETERRQ2(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name);
1770: osm = (PC_GASM*)pc->data;
1771: if (n) *n = osm->n;
1772: if (mat) *mat = osm->pmat;
1773: return(0);
1774: }
1778: /*@
1779: PCGASMSetUseDMSubdomains - Indicates whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1780: Logically Collective
1782: Input Parameter:
1783: + pc - the preconditioner
1784: - flg - boolean indicating whether to use subdomains defined by the DM
1786: Options Database Key:
1787: . -pc_gasm_dm_subdomains -pc_gasm_overlap -pc_gasm_total_subdomains
1789: Level: intermediate
1791: Notes:
1792: PCGASMSetSubdomains(), PCGASMSetTotalSubdomains() or PCGASMSetOverlap() take precedence over PCGASMSetUseDMSubdomains(),
1793: so setting PCGASMSetSubdomains() with nontrivial subdomain ISs or any of PCGASMSetTotalSubdomains() and PCGASMSetOverlap()
1794: automatically turns the latter off.
1796: .keywords: PC, ASM, DM, set, subdomains, additive Schwarz
1798: .seealso: PCGASMGetUseDMSubdomains(), PCGASMSetSubdomains(), PCGASMSetOverlap()
1799: PCGASMCreateSubdomains2D()
1800: @*/
1801: PetscErrorCode PCGASMSetUseDMSubdomains(PC pc,PetscBool flg)
1802: {
1803: PC_GASM *osm = (PC_GASM*)pc->data;
1805: PetscBool match;
1810: if (pc->setupcalled) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for a setup PC.");
1811: PetscObjectTypeCompare((PetscObject)pc,PCGASM,&match);
1812: if (match) {
1813: if (!osm->user_subdomains && osm->N == PETSC_DETERMINE && osm->overlap < 0) {
1814: osm->dm_subdomains = flg;
1815: }
1816: }
1817: return(0);
1818: }
1822: /*@
1823: PCGASMGetUseDMSubdomains - Returns flag indicating whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1824: Not Collective
1826: Input Parameter:
1827: . pc - the preconditioner
1829: Output Parameter:
1830: . flg - boolean indicating whether to use subdomains defined by the DM
1832: Level: intermediate
1834: .keywords: PC, ASM, DM, set, subdomains, additive Schwarz
1836: .seealso: PCGASMSetUseDMSubdomains(), PCGASMSetOverlap()
1837: PCGASMCreateSubdomains2D()
1838: @*/
1839: PetscErrorCode PCGASMGetUseDMSubdomains(PC pc,PetscBool* flg)
1840: {
1841: PC_GASM *osm = (PC_GASM*)pc->data;
1843: PetscBool match;
1848: PetscObjectTypeCompare((PetscObject)pc,PCGASM,&match);
1849: if (match) {
1850: if (flg) *flg = osm->dm_subdomains;
1851: }
1852: return(0);
1853: }