Actual source code: gasm.c
1: /*
2: This file defines an "generalized" additive Schwarz preconditioner for any Mat implementation.
3: In this version each MPI rank may intersect multiple subdomains and any subdomain may
4: intersect multiple MPI ranks. Intersections of subdomains with MPI ranks 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 rank (set in PCGASMSetSubdomains() or calculated in PCGASMSetTotalSubdomains())
9: nmax - maximum number of local subdomains per rank (calculated in PCSetUp_GASM())
10: */
11: #include <petsc/private/pcimpl.h>
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: PetscBool hierarchicalpartitioning;
24: IS *ois; /* index sets that define the outer (conceptually, overlapping) subdomains */
25: IS *iis; /* index sets that define the inner (conceptually, nonoverlapping) subdomains */
26: KSP *ksp; /* linear solvers for each subdomain */
27: Mat *pmat; /* subdomain block matrices */
28: Vec gx, gy; /* Merged work vectors */
29: Vec *x, *y; /* Split work vectors; storage aliases pieces of storage of the above merged vectors. */
30: VecScatter gorestriction; /* merged restriction to disjoint union of outer subdomains */
31: VecScatter girestriction; /* merged restriction to disjoint union of inner subdomains */
32: VecScatter pctoouter;
33: IS permutationIS;
34: Mat permutationP;
35: Mat pcmat;
36: Vec pcx, pcy;
37: } PC_GASM;
39: static PetscErrorCode PCGASMComputeGlobalSubdomainNumbering_Private(PC pc, PetscInt **numbering, PetscInt **permutation)
40: {
41: PC_GASM *osm = (PC_GASM *)pc->data;
42: PetscInt i;
44: /* Determine the number of globally-distinct subdomains and compute a global numbering for them. */
45: PetscMalloc2(osm->n, numbering, osm->n, permutation);
46: PetscObjectsListGetGlobalNumbering(PetscObjectComm((PetscObject)pc), osm->n, (PetscObject *)osm->iis, NULL, *numbering);
47: for (i = 0; i < osm->n; ++i) (*permutation)[i] = i;
48: PetscSortIntWithPermutation(osm->n, *numbering, *permutation);
49: return 0;
50: }
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;
61: /* Inner subdomains. */
62: ISGetLocalSize(osm->iis[i], &nidx);
63: /*
64: No more than 15 characters per index plus a space.
65: PetscViewerStringSPrintf requires a string of size at least 2, so use (nidx+1) instead of nidx,
66: in case nidx == 0. That will take care of the space for the trailing '\0' as well.
67: For nidx == 0, the whole string 16 '\0'.
68: */
69: #define len 16 * (nidx + 1) + 1
70: PetscMalloc1(len, &cidx);
71: PetscViewerStringOpen(PETSC_COMM_SELF, cidx, len, &sviewer);
72: #undef len
73: ISGetIndices(osm->iis[i], &idx);
74: for (j = 0; j < nidx; ++j) PetscViewerStringSPrintf(sviewer, "%" PetscInt_FMT " ", idx[j]);
75: ISRestoreIndices(osm->iis[i], &idx);
76: PetscViewerDestroy(&sviewer);
77: PetscViewerASCIIPrintf(viewer, "Inner subdomain:\n");
78: PetscViewerFlush(viewer);
79: PetscViewerASCIIPushSynchronized(viewer);
80: PetscViewerASCIISynchronizedPrintf(viewer, "%s", cidx);
81: PetscViewerFlush(viewer);
82: PetscViewerASCIIPopSynchronized(viewer);
83: PetscViewerASCIIPrintf(viewer, "\n");
84: PetscViewerFlush(viewer);
85: PetscFree(cidx);
86: /* Outer subdomains. */
87: ISGetLocalSize(osm->ois[i], &nidx);
88: /*
89: No more than 15 characters per index plus a space.
90: PetscViewerStringSPrintf requires a string of size at least 2, so use (nidx+1) instead of nidx,
91: in case nidx == 0. That will take care of the space for the trailing '\0' as well.
92: For nidx == 0, the whole string 16 '\0'.
93: */
94: #define len 16 * (nidx + 1) + 1
95: PetscMalloc1(len, &cidx);
96: PetscViewerStringOpen(PETSC_COMM_SELF, cidx, len, &sviewer);
97: #undef len
98: ISGetIndices(osm->ois[i], &idx);
99: for (j = 0; j < nidx; ++j) PetscViewerStringSPrintf(sviewer, "%" PetscInt_FMT " ", idx[j]);
100: PetscViewerDestroy(&sviewer);
101: ISRestoreIndices(osm->ois[i], &idx);
102: PetscViewerASCIIPrintf(viewer, "Outer subdomain:\n");
103: PetscViewerFlush(viewer);
104: PetscViewerASCIIPushSynchronized(viewer);
105: PetscViewerASCIISynchronizedPrintf(viewer, "%s", cidx);
106: PetscViewerFlush(viewer);
107: PetscViewerASCIIPopSynchronized(viewer);
108: PetscViewerASCIIPrintf(viewer, "\n");
109: PetscViewerFlush(viewer);
110: PetscFree(cidx);
111: return 0;
112: }
114: static PetscErrorCode PCGASMPrintSubdomains(PC pc)
115: {
116: PC_GASM *osm = (PC_GASM *)pc->data;
117: const char *prefix;
118: char fname[PETSC_MAX_PATH_LEN + 1];
119: PetscInt l, d, count;
120: PetscBool found;
121: PetscViewer viewer, sviewer = NULL;
122: PetscInt *numbering, *permutation; /* global numbering of locally-supported subdomains and the permutation from the local ordering */
124: PCGetOptionsPrefix(pc, &prefix);
125: PetscOptionsHasName(NULL, prefix, "-pc_gasm_print_subdomains", &found);
126: if (!found) return 0;
127: PetscOptionsGetString(NULL, prefix, "-pc_gasm_print_subdomains", fname, sizeof(fname), &found);
128: if (!found) PetscStrcpy(fname, "stdout");
129: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc), fname, &viewer);
130: /*
131: Make sure the viewer has a name. Otherwise this may cause a deadlock or other weird errors when creating a subcomm viewer:
132: the subcomm viewer will attempt to inherit the viewer's name, which, if not set, will be constructed collectively on the comm.
133: */
134: PetscObjectName((PetscObject)viewer);
135: l = 0;
136: PCGASMComputeGlobalSubdomainNumbering_Private(pc, &numbering, &permutation);
137: for (count = 0; count < osm->N; ++count) {
138: /* Now let subdomains go one at a time in the global numbering order and print their subdomain/solver info. */
139: if (l < osm->n) {
140: d = permutation[l]; /* d is the local number of the l-th smallest (in the global ordering) among the locally supported subdomains */
141: if (numbering[d] == count) {
142: PetscViewerGetSubViewer(viewer, ((PetscObject)osm->ois[d])->comm, &sviewer);
143: PCGASMSubdomainView_Private(pc, d, sviewer);
144: PetscViewerRestoreSubViewer(viewer, ((PetscObject)osm->ois[d])->comm, &sviewer);
145: ++l;
146: }
147: }
148: MPI_Barrier(PetscObjectComm((PetscObject)pc));
149: }
150: PetscFree2(numbering, permutation);
151: PetscViewerDestroy(&viewer);
152: return 0;
153: }
155: static PetscErrorCode PCView_GASM(PC pc, PetscViewer viewer)
156: {
157: PC_GASM *osm = (PC_GASM *)pc->data;
158: const char *prefix;
159: PetscMPIInt rank, size;
160: PetscInt bsz;
161: PetscBool iascii, view_subdomains = PETSC_FALSE;
162: PetscViewer sviewer;
163: PetscInt count, l;
164: char overlap[256] = "user-defined overlap";
165: char gsubdomains[256] = "unknown total number of subdomains";
166: char msubdomains[256] = "unknown max number of local subdomains";
167: PetscInt *numbering, *permutation; /* global numbering of locally-supported subdomains and the permutation from the local ordering */
169: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
170: MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
172: if (osm->overlap >= 0) PetscSNPrintf(overlap, sizeof(overlap), "requested amount of overlap = %" PetscInt_FMT, osm->overlap);
173: if (osm->N != PETSC_DETERMINE) PetscSNPrintf(gsubdomains, sizeof(gsubdomains), "total number of subdomains = %" PetscInt_FMT, osm->N);
174: if (osm->nmax != PETSC_DETERMINE) PetscSNPrintf(msubdomains, sizeof(msubdomains), "max number of local subdomains = %" PetscInt_FMT, osm->nmax);
176: PCGetOptionsPrefix(pc, &prefix);
177: PetscOptionsGetBool(NULL, prefix, "-pc_gasm_view_subdomains", &view_subdomains, NULL);
179: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
180: if (iascii) {
181: /*
182: Make sure the viewer has a name. Otherwise this may cause a deadlock when creating a subcomm viewer:
183: the subcomm viewer will attempt to inherit the viewer's name, which, if not set, will be constructed
184: collectively on the comm.
185: */
186: PetscObjectName((PetscObject)viewer);
187: PetscViewerASCIIPrintf(viewer, " Restriction/interpolation type: %s\n", PCGASMTypes[osm->type]);
188: PetscViewerASCIIPrintf(viewer, " %s\n", overlap);
189: PetscViewerASCIIPrintf(viewer, " %s\n", gsubdomains);
190: PetscViewerASCIIPrintf(viewer, " %s\n", msubdomains);
191: PetscViewerASCIIPushSynchronized(viewer);
192: PetscViewerASCIISynchronizedPrintf(viewer, " [%d|%d] number of locally-supported subdomains = %" PetscInt_FMT "\n", rank, size, osm->n);
193: PetscViewerFlush(viewer);
194: PetscViewerASCIIPopSynchronized(viewer);
195: /* Cannot take advantage of osm->same_subdomain_solvers without a global numbering of subdomains. */
196: PetscViewerASCIIPrintf(viewer, " Subdomain solver info is as follows:\n");
197: PetscViewerASCIIPushTab(viewer);
198: PetscViewerASCIIPrintf(viewer, " - - - - - - - - - - - - - - - - - -\n");
199: /* Make sure that everybody waits for the banner to be printed. */
200: MPI_Barrier(PetscObjectComm((PetscObject)viewer));
201: /* Now let subdomains go one at a time in the global numbering order and print their subdomain/solver info. */
202: PCGASMComputeGlobalSubdomainNumbering_Private(pc, &numbering, &permutation);
203: l = 0;
204: for (count = 0; count < osm->N; ++count) {
205: PetscMPIInt srank, ssize;
206: if (l < osm->n) {
207: PetscInt d = permutation[l]; /* d is the local number of the l-th smallest (in the global ordering) among the locally supported subdomains */
208: if (numbering[d] == count) {
209: MPI_Comm_size(((PetscObject)osm->ois[d])->comm, &ssize);
210: MPI_Comm_rank(((PetscObject)osm->ois[d])->comm, &srank);
211: PetscViewerGetSubViewer(viewer, ((PetscObject)osm->ois[d])->comm, &sviewer);
212: ISGetLocalSize(osm->ois[d], &bsz);
213: PetscViewerASCIISynchronizedPrintf(sviewer, " [%d|%d] (subcomm [%d|%d]) local subdomain number %" PetscInt_FMT ", local size = %" PetscInt_FMT "\n", rank, size, srank, ssize, d, bsz);
214: PetscViewerFlush(sviewer);
215: if (view_subdomains) PCGASMSubdomainView_Private(pc, d, sviewer);
216: if (!pc->setupcalled) {
217: PetscViewerASCIIPrintf(sviewer, " Solver not set up yet: PCSetUp() not yet called\n");
218: } else {
219: KSPView(osm->ksp[d], sviewer);
220: }
221: PetscViewerASCIIPrintf(sviewer, " - - - - - - - - - - - - - - - - - -\n");
222: PetscViewerFlush(sviewer);
223: PetscViewerRestoreSubViewer(viewer, ((PetscObject)osm->ois[d])->comm, &sviewer);
224: ++l;
225: }
226: }
227: MPI_Barrier(PetscObjectComm((PetscObject)pc));
228: }
229: PetscFree2(numbering, permutation);
230: PetscViewerASCIIPopTab(viewer);
231: PetscViewerFlush(viewer);
232: /* this line is needed to match the extra PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
233: PetscViewerASCIIPopSynchronized(viewer);
234: }
235: return 0;
236: }
238: PETSC_INTERN PetscErrorCode PCGASMCreateLocalSubdomains(Mat A, PetscInt nloc, IS *iis[]);
240: PetscErrorCode PCGASMSetHierarchicalPartitioning(PC pc)
241: {
242: PC_GASM *osm = (PC_GASM *)pc->data;
243: MatPartitioning part;
244: MPI_Comm comm;
245: PetscMPIInt size;
246: PetscInt nlocalsubdomains, fromrows_localsize;
247: IS partitioning, fromrows, isn;
248: Vec outervec;
250: PetscObjectGetComm((PetscObject)pc, &comm);
251: MPI_Comm_size(comm, &size);
252: /* we do not need a hierarchical partitioning when
253: * the total number of subdomains is consistent with
254: * the number of MPI tasks.
255: * For the following cases, we do not need to use HP
256: * */
257: if (osm->N == PETSC_DETERMINE || osm->N >= size || osm->N == 1) return 0;
259: nlocalsubdomains = size / osm->N;
260: osm->n = 1;
261: MatPartitioningCreate(comm, &part);
262: MatPartitioningSetAdjacency(part, pc->pmat);
263: MatPartitioningSetType(part, MATPARTITIONINGHIERARCH);
264: MatPartitioningHierarchicalSetNcoarseparts(part, osm->N);
265: MatPartitioningHierarchicalSetNfineparts(part, nlocalsubdomains);
266: MatPartitioningSetFromOptions(part);
267: /* get new rank owner number of each vertex */
268: MatPartitioningApply(part, &partitioning);
269: ISBuildTwoSided(partitioning, NULL, &fromrows);
270: ISPartitioningToNumbering(partitioning, &isn);
271: ISDestroy(&isn);
272: ISGetLocalSize(fromrows, &fromrows_localsize);
273: MatPartitioningDestroy(&part);
274: MatCreateVecs(pc->pmat, &outervec, NULL);
275: VecCreateMPI(comm, fromrows_localsize, PETSC_DETERMINE, &(osm->pcx));
276: VecDuplicate(osm->pcx, &(osm->pcy));
277: VecScatterCreate(osm->pcx, NULL, outervec, fromrows, &(osm->pctoouter));
278: MatCreateSubMatrix(pc->pmat, fromrows, fromrows, MAT_INITIAL_MATRIX, &(osm->permutationP));
279: PetscObjectReference((PetscObject)fromrows);
280: osm->permutationIS = fromrows;
281: osm->pcmat = pc->pmat;
282: PetscObjectReference((PetscObject)osm->permutationP);
283: pc->pmat = osm->permutationP;
284: VecDestroy(&outervec);
285: ISDestroy(&fromrows);
286: ISDestroy(&partitioning);
287: osm->n = PETSC_DETERMINE;
288: return 0;
289: }
291: static PetscErrorCode PCSetUp_GASM(PC pc)
292: {
293: PC_GASM *osm = (PC_GASM *)pc->data;
294: PetscInt i, nInnerIndices, nTotalInnerIndices;
295: PetscMPIInt rank, size;
296: MatReuse scall = MAT_REUSE_MATRIX;
297: KSP ksp;
298: PC subpc;
299: const char *prefix, *pprefix;
300: Vec x, y;
301: PetscInt oni; /* Number of indices in the i-th local outer subdomain. */
302: const PetscInt *oidxi; /* Indices from the i-th subdomain local outer subdomain. */
303: PetscInt on; /* Number of indices in the disjoint union of local outer subdomains. */
304: PetscInt *oidx; /* Indices in the disjoint union of local outer subdomains. */
305: IS gois; /* Disjoint union the global indices of outer subdomains. */
306: IS goid; /* Identity IS of the size of the disjoint union of outer subdomains. */
307: PetscScalar *gxarray, *gyarray;
308: PetscInt gostart; /* Start of locally-owned indices in the vectors -- osm->gx,osm->gy -- over the disjoint union of outer subdomains. */
309: PetscInt num_subdomains = 0;
310: DM *subdomain_dm = NULL;
311: char **subdomain_names = NULL;
312: PetscInt *numbering;
314: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
315: MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
316: if (!pc->setupcalled) {
317: /* use a hierarchical partitioning */
318: if (osm->hierarchicalpartitioning) PCGASMSetHierarchicalPartitioning(pc);
319: if (osm->n == PETSC_DETERMINE) {
320: if (osm->N != PETSC_DETERMINE) {
321: /* No local subdomains given, but the desired number of total subdomains is known, so construct them accordingly. */
322: PCGASMCreateSubdomains(pc->pmat, osm->N, &osm->n, &osm->iis);
323: } else if (osm->dm_subdomains && pc->dm) {
324: /* try pc->dm next, if allowed */
325: PetscInt d;
326: IS *inner_subdomain_is, *outer_subdomain_is;
327: DMCreateDomainDecomposition(pc->dm, &num_subdomains, &subdomain_names, &inner_subdomain_is, &outer_subdomain_is, &subdomain_dm);
328: if (num_subdomains) PCGASMSetSubdomains(pc, num_subdomains, inner_subdomain_is, outer_subdomain_is);
329: for (d = 0; d < num_subdomains; ++d) {
330: if (inner_subdomain_is) ISDestroy(&inner_subdomain_is[d]);
331: if (outer_subdomain_is) ISDestroy(&outer_subdomain_is[d]);
332: }
333: PetscFree(inner_subdomain_is);
334: PetscFree(outer_subdomain_is);
335: } else {
336: /* still no subdomains; use one per rank */
337: osm->nmax = osm->n = 1;
338: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
339: osm->N = size;
340: PCGASMCreateLocalSubdomains(pc->pmat, osm->n, &osm->iis);
341: }
342: }
343: if (!osm->iis) {
344: /*
345: osm->n was set in PCGASMSetSubdomains(), but the actual subdomains have not been supplied.
346: We create the requisite number of local inner subdomains and then expand them into
347: out subdomains, if necessary.
348: */
349: PCGASMCreateLocalSubdomains(pc->pmat, osm->n, &osm->iis);
350: }
351: if (!osm->ois) {
352: /*
353: Initially make outer subdomains the same as inner subdomains. If nonzero additional overlap
354: has been requested, copy the inner subdomains over so they can be modified.
355: */
356: PetscMalloc1(osm->n, &osm->ois);
357: for (i = 0; i < osm->n; ++i) {
358: if (osm->overlap > 0 && osm->N > 1) { /* With positive overlap, osm->iis[i] will be modified */
359: ISDuplicate(osm->iis[i], (osm->ois) + i);
360: ISCopy(osm->iis[i], osm->ois[i]);
361: } else {
362: PetscObjectReference((PetscObject)((osm->iis)[i]));
363: osm->ois[i] = osm->iis[i];
364: }
365: }
366: if (osm->overlap > 0 && osm->N > 1) {
367: /* Extend the "overlapping" regions by a number of steps */
368: MatIncreaseOverlapSplit(pc->pmat, osm->n, osm->ois, osm->overlap);
369: }
370: }
372: /* Now the subdomains are defined. Determine their global and max local numbers, if necessary. */
373: if (osm->nmax == PETSC_DETERMINE) {
374: PetscMPIInt inwork, outwork;
375: /* determine global number of subdomains and the max number of local subdomains */
376: inwork = osm->n;
377: MPIU_Allreduce(&inwork, &outwork, 1, MPI_INT, MPI_MAX, PetscObjectComm((PetscObject)pc));
378: osm->nmax = outwork;
379: }
380: if (osm->N == PETSC_DETERMINE) {
381: /* Determine the number of globally-distinct subdomains and compute a global numbering for them. */
382: PetscObjectsListGetGlobalNumbering(PetscObjectComm((PetscObject)pc), osm->n, (PetscObject *)osm->ois, &osm->N, NULL);
383: }
385: if (osm->sort_indices) {
386: for (i = 0; i < osm->n; i++) {
387: ISSort(osm->ois[i]);
388: ISSort(osm->iis[i]);
389: }
390: }
391: PCGetOptionsPrefix(pc, &prefix);
392: PCGASMPrintSubdomains(pc);
394: /*
395: Merge the ISs, create merged vectors and restrictions.
396: */
397: /* Merge outer subdomain ISs and construct a restriction onto the disjoint union of local outer subdomains. */
398: on = 0;
399: for (i = 0; i < osm->n; i++) {
400: ISGetLocalSize(osm->ois[i], &oni);
401: on += oni;
402: }
403: PetscMalloc1(on, &oidx);
404: on = 0;
405: /* Merge local indices together */
406: for (i = 0; i < osm->n; i++) {
407: ISGetLocalSize(osm->ois[i], &oni);
408: ISGetIndices(osm->ois[i], &oidxi);
409: PetscArraycpy(oidx + on, oidxi, oni);
410: ISRestoreIndices(osm->ois[i], &oidxi);
411: on += oni;
412: }
413: ISCreateGeneral(((PetscObject)(pc))->comm, on, oidx, PETSC_OWN_POINTER, &gois);
414: nTotalInnerIndices = 0;
415: for (i = 0; i < osm->n; i++) {
416: ISGetLocalSize(osm->iis[i], &nInnerIndices);
417: nTotalInnerIndices += nInnerIndices;
418: }
419: VecCreateMPI(((PetscObject)(pc))->comm, nTotalInnerIndices, PETSC_DETERMINE, &x);
420: VecDuplicate(x, &y);
422: VecCreateMPI(PetscObjectComm((PetscObject)pc), on, PETSC_DECIDE, &osm->gx);
423: VecDuplicate(osm->gx, &osm->gy);
424: VecGetOwnershipRange(osm->gx, &gostart, NULL);
425: ISCreateStride(PetscObjectComm((PetscObject)pc), on, gostart, 1, &goid);
426: /* gois might indices not on local */
427: VecScatterCreate(x, gois, osm->gx, goid, &(osm->gorestriction));
428: PetscMalloc1(osm->n, &numbering);
429: PetscObjectsListGetGlobalNumbering(PetscObjectComm((PetscObject)pc), osm->n, (PetscObject *)osm->ois, NULL, numbering);
430: VecDestroy(&x);
431: ISDestroy(&gois);
433: /* Merge inner subdomain ISs and construct a restriction onto the disjoint union of local inner subdomains. */
434: {
435: PetscInt ini; /* Number of indices the i-th a local inner subdomain. */
436: PetscInt in; /* Number of indices in the disjoint union of local inner subdomains. */
437: PetscInt *iidx; /* Global indices in the merged local inner subdomain. */
438: PetscInt *ioidx; /* Global indices of the disjoint union of inner subdomains within the disjoint union of outer subdomains. */
439: IS giis; /* IS for the disjoint union of inner subdomains. */
440: IS giois; /* IS for the disjoint union of inner subdomains within the disjoint union of outer subdomains. */
441: PetscScalar *array;
442: const PetscInt *indices;
443: PetscInt k;
444: on = 0;
445: for (i = 0; i < osm->n; i++) {
446: ISGetLocalSize(osm->ois[i], &oni);
447: on += oni;
448: }
449: PetscMalloc1(on, &iidx);
450: PetscMalloc1(on, &ioidx);
451: VecGetArray(y, &array);
452: /* set communicator id to determine where overlap is */
453: in = 0;
454: for (i = 0; i < osm->n; i++) {
455: ISGetLocalSize(osm->iis[i], &ini);
456: for (k = 0; k < ini; ++k) array[in + k] = numbering[i];
457: in += ini;
458: }
459: VecRestoreArray(y, &array);
460: VecScatterBegin(osm->gorestriction, y, osm->gy, INSERT_VALUES, SCATTER_FORWARD);
461: VecScatterEnd(osm->gorestriction, y, osm->gy, INSERT_VALUES, SCATTER_FORWARD);
462: VecGetOwnershipRange(osm->gy, &gostart, NULL);
463: VecGetArray(osm->gy, &array);
464: on = 0;
465: in = 0;
466: for (i = 0; i < osm->n; i++) {
467: ISGetLocalSize(osm->ois[i], &oni);
468: ISGetIndices(osm->ois[i], &indices);
469: for (k = 0; k < oni; k++) {
470: /* skip overlapping indices to get inner domain */
471: if (PetscRealPart(array[on + k]) != numbering[i]) continue;
472: iidx[in] = indices[k];
473: ioidx[in++] = gostart + on + k;
474: }
475: ISRestoreIndices(osm->ois[i], &indices);
476: on += oni;
477: }
478: VecRestoreArray(osm->gy, &array);
479: ISCreateGeneral(PetscObjectComm((PetscObject)pc), in, iidx, PETSC_OWN_POINTER, &giis);
480: ISCreateGeneral(PetscObjectComm((PetscObject)pc), in, ioidx, PETSC_OWN_POINTER, &giois);
481: VecScatterCreate(y, giis, osm->gy, giois, &osm->girestriction);
482: VecDestroy(&y);
483: ISDestroy(&giis);
484: ISDestroy(&giois);
485: }
486: ISDestroy(&goid);
487: PetscFree(numbering);
489: /* Create the subdomain work vectors. */
490: PetscMalloc1(osm->n, &osm->x);
491: PetscMalloc1(osm->n, &osm->y);
492: VecGetArray(osm->gx, &gxarray);
493: VecGetArray(osm->gy, &gyarray);
494: for (i = 0, on = 0; i < osm->n; ++i, on += oni) {
495: PetscInt oNi;
496: ISGetLocalSize(osm->ois[i], &oni);
497: /* on a sub communicator */
498: ISGetSize(osm->ois[i], &oNi);
499: VecCreateMPIWithArray(((PetscObject)(osm->ois[i]))->comm, 1, oni, oNi, gxarray + on, &osm->x[i]);
500: VecCreateMPIWithArray(((PetscObject)(osm->ois[i]))->comm, 1, oni, oNi, gyarray + on, &osm->y[i]);
501: }
502: VecRestoreArray(osm->gx, &gxarray);
503: VecRestoreArray(osm->gy, &gyarray);
504: /* Create the subdomain solvers */
505: PetscMalloc1(osm->n, &osm->ksp);
506: for (i = 0; i < osm->n; i++) {
507: char subprefix[PETSC_MAX_PATH_LEN + 1];
508: KSPCreate(((PetscObject)(osm->ois[i]))->comm, &ksp);
509: KSPSetErrorIfNotConverged(ksp, pc->erroriffailure);
510: PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1);
511: KSPSetType(ksp, KSPPREONLY);
512: KSPGetPC(ksp, &subpc); /* Why do we need this here? */
513: if (subdomain_dm) {
514: KSPSetDM(ksp, subdomain_dm[i]);
515: DMDestroy(subdomain_dm + i);
516: }
517: PCGetOptionsPrefix(pc, &prefix);
518: KSPSetOptionsPrefix(ksp, prefix);
519: if (subdomain_names && subdomain_names[i]) {
520: PetscSNPrintf(subprefix, PETSC_MAX_PATH_LEN, "sub_%s_", subdomain_names[i]);
521: KSPAppendOptionsPrefix(ksp, subprefix);
522: PetscFree(subdomain_names[i]);
523: }
524: KSPAppendOptionsPrefix(ksp, "sub_");
525: osm->ksp[i] = ksp;
526: }
527: PetscFree(subdomain_dm);
528: PetscFree(subdomain_names);
529: scall = MAT_INITIAL_MATRIX;
530: } else { /* if (pc->setupcalled) */
531: /*
532: Destroy the submatrices from the previous iteration
533: */
534: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
535: MatDestroyMatrices(osm->n, &osm->pmat);
536: scall = MAT_INITIAL_MATRIX;
537: }
538: if (osm->permutationIS) {
539: MatCreateSubMatrix(pc->pmat, osm->permutationIS, osm->permutationIS, scall, &osm->permutationP);
540: PetscObjectReference((PetscObject)osm->permutationP);
541: osm->pcmat = pc->pmat;
542: pc->pmat = osm->permutationP;
543: }
544: }
546: /*
547: Extract the submatrices.
548: */
549: if (size > 1) {
550: MatCreateSubMatricesMPI(pc->pmat, osm->n, osm->ois, osm->ois, scall, &osm->pmat);
551: } else {
552: MatCreateSubMatrices(pc->pmat, osm->n, osm->ois, osm->ois, scall, &osm->pmat);
553: }
554: if (scall == MAT_INITIAL_MATRIX) {
555: PetscObjectGetOptionsPrefix((PetscObject)pc->pmat, &pprefix);
556: for (i = 0; i < osm->n; i++) { PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i], pprefix); }
557: }
559: /* Return control to the user so that the submatrices can be modified (e.g., to apply
560: different boundary conditions for the submatrices than for the global problem) */
561: PCModifySubMatrices(pc, osm->n, osm->ois, osm->ois, osm->pmat, pc->modifysubmatricesP);
563: /*
564: Loop over submatrices putting them into local ksps
565: */
566: for (i = 0; i < osm->n; i++) {
567: KSPSetOperators(osm->ksp[i], osm->pmat[i], osm->pmat[i]);
568: KSPGetOptionsPrefix(osm->ksp[i], &prefix);
569: MatSetOptionsPrefix(osm->pmat[i], prefix);
570: if (!pc->setupcalled) KSPSetFromOptions(osm->ksp[i]);
571: }
572: if (osm->pcmat) {
573: MatDestroy(&pc->pmat);
574: pc->pmat = osm->pcmat;
575: osm->pcmat = NULL;
576: }
577: return 0;
578: }
580: static PetscErrorCode PCSetUpOnBlocks_GASM(PC pc)
581: {
582: PC_GASM *osm = (PC_GASM *)pc->data;
583: PetscInt i;
585: for (i = 0; i < osm->n; i++) KSPSetUp(osm->ksp[i]);
586: return 0;
587: }
589: static PetscErrorCode PCApply_GASM(PC pc, Vec xin, Vec yout)
590: {
591: PC_GASM *osm = (PC_GASM *)pc->data;
592: PetscInt i;
593: Vec x, y;
594: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
596: if (osm->pctoouter) {
597: VecScatterBegin(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE);
598: VecScatterEnd(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE);
599: x = osm->pcx;
600: y = osm->pcy;
601: } else {
602: x = xin;
603: y = yout;
604: }
605: /*
606: support for limiting the restriction or interpolation only to the inner
607: subdomain values (leaving the other values 0).
608: */
609: if (!(osm->type & PC_GASM_RESTRICT)) {
610: /* have to zero the work RHS since scatter may leave some slots empty */
611: VecZeroEntries(osm->gx);
612: VecScatterBegin(osm->girestriction, x, osm->gx, INSERT_VALUES, forward);
613: } else {
614: VecScatterBegin(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward);
615: }
616: VecZeroEntries(osm->gy);
617: if (!(osm->type & PC_GASM_RESTRICT)) {
618: VecScatterEnd(osm->girestriction, x, osm->gx, INSERT_VALUES, forward);
619: } else {
620: VecScatterEnd(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward);
621: }
622: /* do the subdomain solves */
623: for (i = 0; i < osm->n; ++i) {
624: KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]);
625: KSPCheckSolve(osm->ksp[i], pc, osm->y[i]);
626: }
627: /* do we need to zero y? */
628: VecZeroEntries(y);
629: if (!(osm->type & PC_GASM_INTERPOLATE)) {
630: VecScatterBegin(osm->girestriction, osm->gy, y, ADD_VALUES, reverse);
631: VecScatterEnd(osm->girestriction, osm->gy, y, ADD_VALUES, reverse);
632: } else {
633: VecScatterBegin(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse);
634: VecScatterEnd(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse);
635: }
636: if (osm->pctoouter) {
637: VecScatterBegin(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD);
638: VecScatterEnd(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD);
639: }
640: return 0;
641: }
643: static PetscErrorCode PCMatApply_GASM(PC pc, Mat Xin, Mat Yout)
644: {
645: PC_GASM *osm = (PC_GASM *)pc->data;
646: Mat X, Y, O = NULL, Z, W;
647: Vec x, y;
648: PetscInt i, m, M, N;
649: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
652: MatGetSize(Xin, NULL, &N);
653: if (osm->pctoouter) {
654: VecGetLocalSize(osm->pcx, &m);
655: VecGetSize(osm->pcx, &M);
656: MatCreateDense(PetscObjectComm((PetscObject)osm->ois[0]), m, PETSC_DECIDE, M, N, NULL, &O);
657: for (i = 0; i < N; ++i) {
658: MatDenseGetColumnVecRead(Xin, i, &x);
659: MatDenseGetColumnVecWrite(O, i, &y);
660: VecScatterBegin(osm->pctoouter, x, y, INSERT_VALUES, SCATTER_REVERSE);
661: VecScatterEnd(osm->pctoouter, x, y, INSERT_VALUES, SCATTER_REVERSE);
662: MatDenseRestoreColumnVecWrite(O, i, &y);
663: MatDenseRestoreColumnVecRead(Xin, i, &x);
664: }
665: X = Y = O;
666: } else {
667: X = Xin;
668: Y = Yout;
669: }
670: /*
671: support for limiting the restriction or interpolation only to the inner
672: subdomain values (leaving the other values 0).
673: */
674: VecGetLocalSize(osm->x[0], &m);
675: VecGetSize(osm->x[0], &M);
676: MatCreateDense(PetscObjectComm((PetscObject)osm->ois[0]), m, PETSC_DECIDE, M, N, NULL, &Z);
677: for (i = 0; i < N; ++i) {
678: MatDenseGetColumnVecRead(X, i, &x);
679: MatDenseGetColumnVecWrite(Z, i, &y);
680: if (!(osm->type & PC_GASM_RESTRICT)) {
681: /* have to zero the work RHS since scatter may leave some slots empty */
682: VecZeroEntries(y);
683: VecScatterBegin(osm->girestriction, x, y, INSERT_VALUES, forward);
684: VecScatterEnd(osm->girestriction, x, y, INSERT_VALUES, forward);
685: } else {
686: VecScatterBegin(osm->gorestriction, x, y, INSERT_VALUES, forward);
687: VecScatterEnd(osm->gorestriction, x, y, INSERT_VALUES, forward);
688: }
689: MatDenseRestoreColumnVecWrite(Z, i, &y);
690: MatDenseRestoreColumnVecRead(X, i, &x);
691: }
692: MatCreateDense(PetscObjectComm((PetscObject)osm->ois[0]), m, PETSC_DECIDE, M, N, NULL, &W);
693: MatSetOption(Z, MAT_NO_OFF_PROC_ENTRIES, PETSC_TRUE);
694: MatAssemblyBegin(Z, MAT_FINAL_ASSEMBLY);
695: MatAssemblyEnd(Z, MAT_FINAL_ASSEMBLY);
696: /* do the subdomain solve */
697: KSPMatSolve(osm->ksp[0], Z, W);
698: KSPCheckSolve(osm->ksp[0], pc, NULL);
699: MatDestroy(&Z);
700: /* do we need to zero y? */
701: MatZeroEntries(Y);
702: for (i = 0; i < N; ++i) {
703: MatDenseGetColumnVecWrite(Y, i, &y);
704: MatDenseGetColumnVecRead(W, i, &x);
705: if (!(osm->type & PC_GASM_INTERPOLATE)) {
706: VecScatterBegin(osm->girestriction, x, y, ADD_VALUES, reverse);
707: VecScatterEnd(osm->girestriction, x, y, ADD_VALUES, reverse);
708: } else {
709: VecScatterBegin(osm->gorestriction, x, y, ADD_VALUES, reverse);
710: VecScatterEnd(osm->gorestriction, x, y, ADD_VALUES, reverse);
711: }
712: MatDenseRestoreColumnVecRead(W, i, &x);
713: if (osm->pctoouter) {
714: MatDenseGetColumnVecWrite(Yout, i, &x);
715: VecScatterBegin(osm->pctoouter, y, x, INSERT_VALUES, SCATTER_FORWARD);
716: VecScatterEnd(osm->pctoouter, y, x, INSERT_VALUES, SCATTER_FORWARD);
717: MatDenseRestoreColumnVecRead(Yout, i, &x);
718: }
719: MatDenseRestoreColumnVecWrite(Y, i, &y);
720: }
721: MatDestroy(&W);
722: MatDestroy(&O);
723: return 0;
724: }
726: static PetscErrorCode PCApplyTranspose_GASM(PC pc, Vec xin, Vec yout)
727: {
728: PC_GASM *osm = (PC_GASM *)pc->data;
729: PetscInt i;
730: Vec x, y;
731: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
733: if (osm->pctoouter) {
734: VecScatterBegin(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE);
735: VecScatterEnd(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE);
736: x = osm->pcx;
737: y = osm->pcy;
738: } else {
739: x = xin;
740: y = yout;
741: }
742: /*
743: Support for limiting the restriction or interpolation to only local
744: subdomain values (leaving the other values 0).
746: Note: these are reversed from the PCApply_GASM() because we are applying the
747: transpose of the three terms
748: */
749: if (!(osm->type & PC_GASM_INTERPOLATE)) {
750: /* have to zero the work RHS since scatter may leave some slots empty */
751: VecZeroEntries(osm->gx);
752: VecScatterBegin(osm->girestriction, x, osm->gx, INSERT_VALUES, forward);
753: } else {
754: VecScatterBegin(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward);
755: }
756: VecZeroEntries(osm->gy);
757: if (!(osm->type & PC_GASM_INTERPOLATE)) {
758: VecScatterEnd(osm->girestriction, x, osm->gx, INSERT_VALUES, forward);
759: } else {
760: VecScatterEnd(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward);
761: }
762: /* do the local solves */
763: for (i = 0; i < osm->n; ++i) { /* Note that the solves are local, so we can go to osm->n, rather than osm->nmax. */
764: KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i]);
765: KSPCheckSolve(osm->ksp[i], pc, osm->y[i]);
766: }
767: VecZeroEntries(y);
768: if (!(osm->type & PC_GASM_RESTRICT)) {
769: VecScatterBegin(osm->girestriction, osm->gy, y, ADD_VALUES, reverse);
770: VecScatterEnd(osm->girestriction, osm->gy, y, ADD_VALUES, reverse);
771: } else {
772: VecScatterBegin(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse);
773: VecScatterEnd(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse);
774: }
775: if (osm->pctoouter) {
776: VecScatterBegin(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD);
777: VecScatterEnd(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD);
778: }
779: return 0;
780: }
782: static PetscErrorCode PCReset_GASM(PC pc)
783: {
784: PC_GASM *osm = (PC_GASM *)pc->data;
785: PetscInt i;
787: if (osm->ksp) {
788: for (i = 0; i < osm->n; i++) KSPReset(osm->ksp[i]);
789: }
790: if (osm->pmat) {
791: if (osm->n > 0) {
792: PetscMPIInt size;
793: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
794: if (size > 1) {
795: /* osm->pmat is created by MatCreateSubMatricesMPI(), cannot use MatDestroySubMatrices() */
796: MatDestroyMatrices(osm->n, &osm->pmat);
797: } else {
798: MatDestroySubMatrices(osm->n, &osm->pmat);
799: }
800: }
801: }
802: if (osm->x) {
803: for (i = 0; i < osm->n; i++) {
804: VecDestroy(&osm->x[i]);
805: VecDestroy(&osm->y[i]);
806: }
807: }
808: VecDestroy(&osm->gx);
809: VecDestroy(&osm->gy);
811: VecScatterDestroy(&osm->gorestriction);
812: VecScatterDestroy(&osm->girestriction);
813: if (!osm->user_subdomains) {
814: PCGASMDestroySubdomains(osm->n, &osm->ois, &osm->iis);
815: osm->N = PETSC_DETERMINE;
816: osm->nmax = PETSC_DETERMINE;
817: }
818: if (osm->pctoouter) VecScatterDestroy(&(osm->pctoouter));
819: if (osm->permutationIS) ISDestroy(&(osm->permutationIS));
820: if (osm->pcx) VecDestroy(&(osm->pcx));
821: if (osm->pcy) VecDestroy(&(osm->pcy));
822: if (osm->permutationP) MatDestroy(&(osm->permutationP));
823: if (osm->pcmat) MatDestroy(&osm->pcmat);
824: return 0;
825: }
827: static PetscErrorCode PCDestroy_GASM(PC pc)
828: {
829: PC_GASM *osm = (PC_GASM *)pc->data;
830: PetscInt i;
832: PCReset_GASM(pc);
833: /* PCReset will not destroy subdomains, if user_subdomains is true. */
834: PCGASMDestroySubdomains(osm->n, &osm->ois, &osm->iis);
835: if (osm->ksp) {
836: for (i = 0; i < osm->n; i++) KSPDestroy(&osm->ksp[i]);
837: PetscFree(osm->ksp);
838: }
839: PetscFree(osm->x);
840: PetscFree(osm->y);
841: PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSubdomains_C", NULL);
842: PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetOverlap_C", NULL);
843: PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetType_C", NULL);
844: PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSortIndices_C", NULL);
845: PetscObjectComposeFunction((PetscObject)pc, "PCGASMGetSubKSP_C", NULL);
846: PetscFree(pc->data);
847: return 0;
848: }
850: static PetscErrorCode PCSetFromOptions_GASM(PC pc, PetscOptionItems *PetscOptionsObject)
851: {
852: PC_GASM *osm = (PC_GASM *)pc->data;
853: PetscInt blocks, ovl;
854: PetscBool flg;
855: PCGASMType gasmtype;
857: PetscOptionsHeadBegin(PetscOptionsObject, "Generalized additive Schwarz options");
858: PetscOptionsBool("-pc_gasm_use_dm_subdomains", "If subdomains aren't set, use DMCreateDomainDecomposition() to define subdomains.", "PCGASMSetUseDMSubdomains", osm->dm_subdomains, &osm->dm_subdomains, &flg);
859: PetscOptionsInt("-pc_gasm_total_subdomains", "Total number of subdomains across communicator", "PCGASMSetTotalSubdomains", osm->N, &blocks, &flg);
860: if (flg) PCGASMSetTotalSubdomains(pc, blocks);
861: PetscOptionsInt("-pc_gasm_overlap", "Number of overlapping degrees of freedom", "PCGASMSetOverlap", osm->overlap, &ovl, &flg);
862: if (flg) {
863: PCGASMSetOverlap(pc, ovl);
864: osm->dm_subdomains = PETSC_FALSE;
865: }
866: flg = PETSC_FALSE;
867: PetscOptionsEnum("-pc_gasm_type", "Type of restriction/extension", "PCGASMSetType", PCGASMTypes, (PetscEnum)osm->type, (PetscEnum *)&gasmtype, &flg);
868: if (flg) PCGASMSetType(pc, gasmtype);
869: PetscOptionsBool("-pc_gasm_use_hierachical_partitioning", "use hierarchical partitioning", NULL, osm->hierarchicalpartitioning, &osm->hierarchicalpartitioning, &flg);
870: PetscOptionsHeadEnd();
871: return 0;
872: }
874: /*------------------------------------------------------------------------------------*/
876: /*@
877: PCGASMSetTotalSubdomains - sets the total number of subdomains to use across the communicator for `PCGASM`
878: Logically collective on pc
880: Input Parameters:
881: + pc - the preconditioner
882: - N - total number of subdomains
884: Level: beginner
886: .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMSetOverlap()`
887: `PCGASMCreateSubdomains2D()`
888: @*/
889: PetscErrorCode PCGASMSetTotalSubdomains(PC pc, PetscInt N)
890: {
891: PC_GASM *osm = (PC_GASM *)pc->data;
892: PetscMPIInt size, rank;
897: PCGASMDestroySubdomains(osm->n, &osm->iis, &osm->ois);
898: osm->ois = osm->iis = NULL;
900: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
901: MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
902: osm->N = N;
903: osm->n = PETSC_DETERMINE;
904: osm->nmax = PETSC_DETERMINE;
905: osm->dm_subdomains = PETSC_FALSE;
906: return 0;
907: }
909: static PetscErrorCode PCGASMSetSubdomains_GASM(PC pc, PetscInt n, IS iis[], IS ois[])
910: {
911: PC_GASM *osm = (PC_GASM *)pc->data;
912: PetscInt i;
917: PCGASMDestroySubdomains(osm->n, &osm->iis, &osm->ois);
918: osm->iis = osm->ois = NULL;
919: osm->n = n;
920: osm->N = PETSC_DETERMINE;
921: osm->nmax = PETSC_DETERMINE;
922: if (ois) {
923: PetscMalloc1(n, &osm->ois);
924: for (i = 0; i < n; i++) {
925: PetscObjectReference((PetscObject)ois[i]);
926: osm->ois[i] = ois[i];
927: }
928: /*
929: Since the user set the outer subdomains, even if nontrivial overlap was requested via PCGASMSetOverlap(),
930: it will be ignored. To avoid confusion later on (e.g., when viewing the PC), the overlap size is set to -1.
931: */
932: osm->overlap = -1;
933: /* inner subdomains must be provided */
935: } /* end if */
936: if (iis) {
937: PetscMalloc1(n, &osm->iis);
938: for (i = 0; i < n; i++) {
939: PetscObjectReference((PetscObject)iis[i]);
940: osm->iis[i] = iis[i];
941: }
942: if (!ois) {
943: osm->ois = NULL;
944: /* if user does not provide outer indices, we will create the corresponding outer indices using osm->overlap =1 in PCSetUp_GASM */
945: }
946: }
947: if (PetscDefined(USE_DEBUG)) {
948: PetscInt j, rstart, rend, *covered, lsize;
949: const PetscInt *indices;
950: /* check if the inner indices cover and only cover the local portion of the preconditioning matrix */
951: MatGetOwnershipRange(pc->pmat, &rstart, &rend);
952: PetscCalloc1(rend - rstart, &covered);
953: /* check if the current MPI rank owns indices from others */
954: for (i = 0; i < n; i++) {
955: ISGetIndices(osm->iis[i], &indices);
956: ISGetLocalSize(osm->iis[i], &lsize);
957: for (j = 0; j < lsize; j++) {
960: covered[indices[j] - rstart] = 1;
961: }
962: ISRestoreIndices(osm->iis[i], &indices);
963: }
964: /* check if we miss any indices */
966: PetscFree(covered);
967: }
968: if (iis) osm->user_subdomains = PETSC_TRUE;
969: return 0;
970: }
972: static PetscErrorCode PCGASMSetOverlap_GASM(PC pc, PetscInt ovl)
973: {
974: PC_GASM *osm = (PC_GASM *)pc->data;
978: if (!pc->setupcalled) osm->overlap = ovl;
979: return 0;
980: }
982: static PetscErrorCode PCGASMSetType_GASM(PC pc, PCGASMType type)
983: {
984: PC_GASM *osm = (PC_GASM *)pc->data;
986: osm->type = type;
987: osm->type_set = PETSC_TRUE;
988: return 0;
989: }
991: static PetscErrorCode PCGASMSetSortIndices_GASM(PC pc, PetscBool doSort)
992: {
993: PC_GASM *osm = (PC_GASM *)pc->data;
995: osm->sort_indices = doSort;
996: return 0;
997: }
999: /*
1000: FIXME: This routine might need to be modified now that multiple ranks per subdomain are allowed.
1001: In particular, it would upset the global subdomain number calculation.
1002: */
1003: static PetscErrorCode PCGASMGetSubKSP_GASM(PC pc, PetscInt *n, PetscInt *first, KSP **ksp)
1004: {
1005: PC_GASM *osm = (PC_GASM *)pc->data;
1009: if (n) *n = osm->n;
1010: if (first) {
1011: MPI_Scan(&osm->n, first, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)pc));
1012: *first -= osm->n;
1013: }
1014: if (ksp) {
1015: /* Assume that local solves are now different; not necessarily
1016: true, though! This flag is used only for PCView_GASM() */
1017: *ksp = osm->ksp;
1018: osm->same_subdomain_solvers = PETSC_FALSE;
1019: }
1020: return 0;
1021: } /* PCGASMGetSubKSP_GASM() */
1023: /*@C
1024: PCGASMSetSubdomains - Sets the subdomains for this MPI rank
1025: for the additive Schwarz preconditioner with multiple MPI ranks per subdomain, `PCGASM`
1027: Collective on pc
1029: Input Parameters:
1030: + pc - the preconditioner object
1031: . n - the number of subdomains for this MPI rank
1032: . iis - the index sets that define the inner subdomains (or NULL for PETSc to determine subdomains)
1033: - 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)
1035: Notes:
1036: The `IS` indices use the parallel, global numbering of the vector entries.
1037: Inner subdomains are those where the correction is applied.
1038: Outer subdomains are those where the residual necessary to obtain the
1039: corrections is obtained (see `PCGASMType` for the use of inner/outer subdomains).
1040: Both inner and outer subdomains can extend over several MPI ranks.
1041: This rank's portion of a subdomain is known as a local subdomain.
1043: Inner subdomains can not overlap with each other, do not have any entities from remote ranks,
1044: and have to cover the entire local subdomain owned by the current rank. The index sets on each
1045: rank should be ordered such that the ith local subdomain is connected to the ith remote subdomain
1046: on another MPI rank.
1048: By default the `PGASM` preconditioner uses 1 (local) subdomain per MPI rank.
1050: Level: advanced
1052: .seealso: `PCGASM`, `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`,
1053: `PCGASMCreateSubdomains2D()`, `PCGASMGetSubdomains()`
1054: @*/
1055: PetscErrorCode PCGASMSetSubdomains(PC pc, PetscInt n, IS iis[], IS ois[])
1056: {
1057: PC_GASM *osm = (PC_GASM *)pc->data;
1060: PetscTryMethod(pc, "PCGASMSetSubdomains_C", (PC, PetscInt, IS[], IS[]), (pc, n, iis, ois));
1061: osm->dm_subdomains = PETSC_FALSE;
1062: return 0;
1063: }
1065: /*@
1066: PCGASMSetOverlap - Sets the overlap between a pair of subdomains for the
1067: additive Schwarz preconditioner `PCGASM`. Either all or no MPI ranks in the
1068: pc communicator must call this routine.
1070: Logically Collective on pc
1072: Input Parameters:
1073: + pc - the preconditioner context
1074: - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 0)
1076: Options Database Key:
1077: . -pc_gasm_overlap <overlap> - Sets overlap
1079: Notes:
1080: By default the `PCGASM` preconditioner uses 1 subdomain per rank. To use
1081: multiple subdomain per perocessor or "straddling" subdomains that intersect
1082: multiple ranks use `PCGASMSetSubdomains()` (or option -pc_gasm_total_subdomains <n>).
1084: The overlap defaults to 0, so if one desires that no additional
1085: overlap be computed beyond what may have been set with a call to
1086: `PCGASMSetSubdomains()`, then ovl must be set to be 0. In particular, if one does
1087: not explicitly set the subdomains in application code, then all overlap would be computed
1088: internally by PETSc, and using an overlap of 0 would result in an `PCGASM`
1089: variant that is equivalent to the block Jacobi preconditioner.
1091: One can define initial index sets with any overlap via
1092: `PCGASMSetSubdomains()`; the routine `PCGASMSetOverlap()` merely allows
1093: PETSc to extend that overlap further, if desired.
1095: Level: intermediate
1097: .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`,
1098: `PCGASMCreateSubdomains2D()`, `PCGASMGetSubdomains()`
1099: @*/
1100: PetscErrorCode PCGASMSetOverlap(PC pc, PetscInt ovl)
1101: {
1102: PC_GASM *osm = (PC_GASM *)pc->data;
1106: PetscTryMethod(pc, "PCGASMSetOverlap_C", (PC, PetscInt), (pc, ovl));
1107: osm->dm_subdomains = PETSC_FALSE;
1108: return 0;
1109: }
1111: /*@
1112: PCGASMSetType - Sets the type of restriction and interpolation used
1113: for local problems in the `PCGASM` additive Schwarz method.
1115: Logically Collective on pc
1117: Input Parameters:
1118: + pc - the preconditioner context
1119: - type - variant of `PCGASM`, one of
1120: .vb
1121: `PC_GASM_BASIC` - full interpolation and restriction
1122: `PC_GASM_RESTRICT` - full restriction, local MPI rank interpolation
1123: `PC_GASM_INTERPOLATE` - full interpolation, local MPI rank restriction
1124: `PC_GASM_NONE` - local MPI rank restriction and interpolation
1125: .ve
1127: Options Database Key:
1128: . -pc_gasm_type [basic,restrict,interpolate,none] - Sets GASM type
1130: Level: intermediate
1132: .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`,
1133: `PCGASMCreateSubdomains2D()`, `PCASM`, `PCASMSetType()`
1134: @*/
1135: PetscErrorCode PCGASMSetType(PC pc, PCGASMType type)
1136: {
1139: PetscTryMethod(pc, "PCGASMSetType_C", (PC, PCGASMType), (pc, type));
1140: return 0;
1141: }
1143: /*@
1144: PCGASMSetSortIndices - Determines whether subdomain indices are sorted.
1146: Logically Collective on pc
1148: Input Parameters:
1149: + pc - the preconditioner context
1150: - doSort - sort the subdomain indices
1152: Level: intermediate
1154: .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`,
1155: `PCGASMCreateSubdomains2D()`
1156: @*/
1157: PetscErrorCode PCGASMSetSortIndices(PC pc, PetscBool doSort)
1158: {
1161: PetscTryMethod(pc, "PCGASMSetSortIndices_C", (PC, PetscBool), (pc, doSort));
1162: return 0;
1163: }
1165: /*@C
1166: PCGASMGetSubKSP - Gets the local `KSP` contexts for all subdomains on
1167: this MPI rank.
1169: Collective on pc iff first_local is requested
1171: Input Parameter:
1172: . pc - the preconditioner context
1174: Output Parameters:
1175: + n_local - the number of blocks on this MPI rank or NULL
1176: . first_local - the global number of the first block on this rank or NULL,
1177: all ranks must request or all must pass NULL
1178: - ksp - the array of `KSP` contexts
1180: Note:
1181: After `PCGASMGetSubKSP()` the array of `KSP`es is not to be freed
1183: Currently for some matrix implementations only 1 block per MPI rank
1184: is supported.
1186: You must call `KSPSetUp()` before calling `PCGASMGetSubKSP()`.
1188: Level: advanced
1190: .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMSetOverlap()`,
1191: `PCGASMCreateSubdomains2D()`,
1192: @*/
1193: PetscErrorCode PCGASMGetSubKSP(PC pc, PetscInt *n_local, PetscInt *first_local, KSP *ksp[])
1194: {
1196: PetscUseMethod(pc, "PCGASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (pc, n_local, first_local, ksp));
1197: return 0;
1198: }
1200: /*MC
1201: PCGASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1202: its own `KSP` object on a subset of MPI ranks
1204: Options Database Keys:
1205: + -pc_gasm_total_subdomains <n> - Sets total number of local subdomains to be distributed among MPI ranks
1206: . -pc_gasm_view_subdomains - activates the printing of subdomain indices in `PCView()`, -ksp_view or -snes_view
1207: . -pc_gasm_print_subdomains - activates the printing of subdomain indices in `PCSetUp()`
1208: . -pc_gasm_overlap <ovl> - Sets overlap by which to (automatically) extend local subdomains
1209: - -pc_gasm_type [basic,restrict,interpolate,none] - Sets `PCGASMType`
1211: Notes:
1212: To set options on the solvers for each block append -sub_ to all the `KSP`, and `PC`
1213: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
1215: To set the options on the solvers separate for each block call `PCGASMGetSubKSP()`
1216: and set the options directly on the resulting `KSP` object (you can access its `PC`
1217: with `KSPGetPC())`
1219: Level: beginner
1221: References:
1222: + * - M Dryja, OB Widlund, An additive variant of the Schwarz alternating method for the case of many subregions
1223: Courant Institute, New York University Technical report
1224: - * - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations,
1225: Cambridge University Press.
1227: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCASM`, `PCGASMType`, `PCGASMSetType()`,
1228: `PCBJACOBI`, `PCGASMGetSubKSP()`, `PCGASMSetSubdomains()`,
1229: `PCSetModifySubMatrices()`, `PCGASMSetOverlap()`, `PCGASMSetType()`
1230: M*/
1232: PETSC_EXTERN PetscErrorCode PCCreate_GASM(PC pc)
1233: {
1234: PC_GASM *osm;
1236: PetscNew(&osm);
1238: osm->N = PETSC_DETERMINE;
1239: osm->n = PETSC_DECIDE;
1240: osm->nmax = PETSC_DETERMINE;
1241: osm->overlap = 0;
1242: osm->ksp = NULL;
1243: osm->gorestriction = NULL;
1244: osm->girestriction = NULL;
1245: osm->pctoouter = NULL;
1246: osm->gx = NULL;
1247: osm->gy = NULL;
1248: osm->x = NULL;
1249: osm->y = NULL;
1250: osm->pcx = NULL;
1251: osm->pcy = NULL;
1252: osm->permutationIS = NULL;
1253: osm->permutationP = NULL;
1254: osm->pcmat = NULL;
1255: osm->ois = NULL;
1256: osm->iis = NULL;
1257: osm->pmat = NULL;
1258: osm->type = PC_GASM_RESTRICT;
1259: osm->same_subdomain_solvers = PETSC_TRUE;
1260: osm->sort_indices = PETSC_TRUE;
1261: osm->dm_subdomains = PETSC_FALSE;
1262: osm->hierarchicalpartitioning = PETSC_FALSE;
1264: pc->data = (void *)osm;
1265: pc->ops->apply = PCApply_GASM;
1266: pc->ops->matapply = PCMatApply_GASM;
1267: pc->ops->applytranspose = PCApplyTranspose_GASM;
1268: pc->ops->setup = PCSetUp_GASM;
1269: pc->ops->reset = PCReset_GASM;
1270: pc->ops->destroy = PCDestroy_GASM;
1271: pc->ops->setfromoptions = PCSetFromOptions_GASM;
1272: pc->ops->setuponblocks = PCSetUpOnBlocks_GASM;
1273: pc->ops->view = PCView_GASM;
1274: pc->ops->applyrichardson = NULL;
1276: PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSubdomains_C", PCGASMSetSubdomains_GASM);
1277: PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetOverlap_C", PCGASMSetOverlap_GASM);
1278: PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetType_C", PCGASMSetType_GASM);
1279: PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSortIndices_C", PCGASMSetSortIndices_GASM);
1280: PetscObjectComposeFunction((PetscObject)pc, "PCGASMGetSubKSP_C", PCGASMGetSubKSP_GASM);
1281: return 0;
1282: }
1284: PetscErrorCode PCGASMCreateLocalSubdomains(Mat A, PetscInt nloc, IS *iis[])
1285: {
1286: MatPartitioning mpart;
1287: const char *prefix;
1288: PetscInt i, j, rstart, rend, bs;
1289: PetscBool hasop, isbaij = PETSC_FALSE, foundpart = PETSC_FALSE;
1290: Mat Ad = NULL, adj;
1291: IS ispart, isnumb, *is;
1295: /* Get prefix, row distribution, and block size */
1296: MatGetOptionsPrefix(A, &prefix);
1297: MatGetOwnershipRange(A, &rstart, &rend);
1298: MatGetBlockSize(A, &bs);
1301: /* Get diagonal block from matrix if possible */
1302: MatHasOperation(A, MATOP_GET_DIAGONAL_BLOCK, &hasop);
1303: if (hasop) MatGetDiagonalBlock(A, &Ad);
1304: if (Ad) {
1305: PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQBAIJ, &isbaij);
1306: if (!isbaij) PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQSBAIJ, &isbaij);
1307: }
1308: if (Ad && nloc > 1) {
1309: PetscBool match, done;
1310: /* Try to setup a good matrix partitioning if available */
1311: MatPartitioningCreate(PETSC_COMM_SELF, &mpart);
1312: PetscObjectSetOptionsPrefix((PetscObject)mpart, prefix);
1313: MatPartitioningSetFromOptions(mpart);
1314: PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGCURRENT, &match);
1315: if (!match) PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGSQUARE, &match);
1316: if (!match) { /* assume a "good" partitioner is available */
1317: PetscInt na;
1318: const PetscInt *ia, *ja;
1319: MatGetRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done);
1320: if (done) {
1321: /* Build adjacency matrix by hand. Unfortunately a call to
1322: MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1323: remove the block-aij structure and we cannot expect
1324: MatPartitioning to split vertices as we need */
1325: PetscInt i, j, len, nnz, cnt, *iia = NULL, *jja = NULL;
1326: const PetscInt *row;
1327: nnz = 0;
1328: for (i = 0; i < na; i++) { /* count number of nonzeros */
1329: len = ia[i + 1] - ia[i];
1330: row = ja + ia[i];
1331: for (j = 0; j < len; j++) {
1332: if (row[j] == i) { /* don't count diagonal */
1333: len--;
1334: break;
1335: }
1336: }
1337: nnz += len;
1338: }
1339: PetscMalloc1(na + 1, &iia);
1340: PetscMalloc1(nnz, &jja);
1341: nnz = 0;
1342: iia[0] = 0;
1343: for (i = 0; i < na; i++) { /* fill adjacency */
1344: cnt = 0;
1345: len = ia[i + 1] - ia[i];
1346: row = ja + ia[i];
1347: for (j = 0; j < len; j++) {
1348: if (row[j] != i) jja[nnz + cnt++] = row[j]; /* if not diagonal */
1349: }
1350: nnz += cnt;
1351: iia[i + 1] = nnz;
1352: }
1353: /* Partitioning of the adjacency matrix */
1354: MatCreateMPIAdj(PETSC_COMM_SELF, na, na, iia, jja, NULL, &adj);
1355: MatPartitioningSetAdjacency(mpart, adj);
1356: MatPartitioningSetNParts(mpart, nloc);
1357: MatPartitioningApply(mpart, &ispart);
1358: ISPartitioningToNumbering(ispart, &isnumb);
1359: MatDestroy(&adj);
1360: foundpart = PETSC_TRUE;
1361: }
1362: MatRestoreRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done);
1363: }
1364: MatPartitioningDestroy(&mpart);
1365: }
1366: PetscMalloc1(nloc, &is);
1367: if (!foundpart) {
1368: /* Partitioning by contiguous chunks of rows */
1370: PetscInt mbs = (rend - rstart) / bs;
1371: PetscInt start = rstart;
1372: for (i = 0; i < nloc; i++) {
1373: PetscInt count = (mbs / nloc + ((mbs % nloc) > i)) * bs;
1374: ISCreateStride(PETSC_COMM_SELF, count, start, 1, &is[i]);
1375: start += count;
1376: }
1378: } else {
1379: /* Partitioning by adjacency of diagonal block */
1381: const PetscInt *numbering;
1382: PetscInt *count, nidx, *indices, *newidx, start = 0;
1383: /* Get node count in each partition */
1384: PetscMalloc1(nloc, &count);
1385: ISPartitioningCount(ispart, nloc, count);
1386: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1387: for (i = 0; i < nloc; i++) count[i] *= bs;
1388: }
1389: /* Build indices from node numbering */
1390: ISGetLocalSize(isnumb, &nidx);
1391: PetscMalloc1(nidx, &indices);
1392: for (i = 0; i < nidx; i++) indices[i] = i; /* needs to be initialized */
1393: ISGetIndices(isnumb, &numbering);
1394: PetscSortIntWithPermutation(nidx, numbering, indices);
1395: ISRestoreIndices(isnumb, &numbering);
1396: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1397: PetscMalloc1(nidx * bs, &newidx);
1398: for (i = 0; i < nidx; i++) {
1399: for (j = 0; j < bs; j++) newidx[i * bs + j] = indices[i] * bs + j;
1400: }
1401: PetscFree(indices);
1402: nidx *= bs;
1403: indices = newidx;
1404: }
1405: /* Shift to get global indices */
1406: for (i = 0; i < nidx; i++) indices[i] += rstart;
1408: /* Build the index sets for each block */
1409: for (i = 0; i < nloc; i++) {
1410: ISCreateGeneral(PETSC_COMM_SELF, count[i], &indices[start], PETSC_COPY_VALUES, &is[i]);
1411: ISSort(is[i]);
1412: start += count[i];
1413: }
1415: PetscFree(count);
1416: PetscFree(indices);
1417: ISDestroy(&isnumb);
1418: ISDestroy(&ispart);
1419: }
1420: *iis = is;
1421: return 0;
1422: }
1424: PETSC_INTERN PetscErrorCode PCGASMCreateStraddlingSubdomains(Mat A, PetscInt N, PetscInt *n, IS *iis[])
1425: {
1426: MatSubdomainsCreateCoalesce(A, N, n, iis);
1427: return 0;
1428: }
1430: /*@C
1431: PCGASMCreateSubdomains - Creates n index sets defining n nonoverlapping subdomains for the `PCGASM` additive
1432: Schwarz preconditioner for a any problem based on its matrix.
1434: Collective
1436: Input Parameters:
1437: + A - The global matrix operator
1438: - N - the number of global subdomains requested
1440: Output Parameters:
1441: + n - the number of subdomains created on this MPI rank
1442: - iis - the array of index sets defining the local inner subdomains (on which the correction is applied)
1444: Level: advanced
1446: Note:
1447: When N >= A's communicator size, each subdomain is local -- contained within a single MPI rank.
1448: When N < size, the subdomains are 'straddling' (rank boundaries) and are no longer local.
1449: The resulting subdomains can be use in `PCGASMSetSubdomains`(pc,n,iss,NULL). The overlapping
1450: outer subdomains will be automatically generated from these according to the requested amount of
1451: overlap; this is currently supported only with local subdomains.
1453: .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMDestroySubdomains()`
1454: @*/
1455: PetscErrorCode PCGASMCreateSubdomains(Mat A, PetscInt N, PetscInt *n, IS *iis[])
1456: {
1457: PetscMPIInt size;
1463: MPI_Comm_size(PetscObjectComm((PetscObject)A), &size);
1464: if (N >= size) {
1465: *n = N / size + (N % size);
1466: PCGASMCreateLocalSubdomains(A, *n, iis);
1467: } else {
1468: PCGASMCreateStraddlingSubdomains(A, N, n, iis);
1469: }
1470: return 0;
1471: }
1473: /*@C
1474: PCGASMDestroySubdomains - Destroys the index sets created with
1475: `PCGASMCreateSubdomains()` or `PCGASMCreateSubdomains2D()`. Should be
1476: called after setting subdomains with `PCGASMSetSubdomains()`.
1478: Collective
1480: Input Parameters:
1481: + n - the number of index sets
1482: . iis - the array of inner subdomains,
1483: - ois - the array of outer subdomains, can be NULL
1485: Level: intermediate
1487: Note:
1488: This is a convenience subroutine that walks each list,
1489: destroys each `IS` on the list, and then frees the list. At the end the
1490: list pointers are set to NULL.
1492: .seealso: `PCGASM`, `PCGASMCreateSubdomains()`, `PCGASMSetSubdomains()`
1493: @*/
1494: PetscErrorCode PCGASMDestroySubdomains(PetscInt n, IS **iis, IS **ois)
1495: {
1496: PetscInt i;
1498: if (n <= 0) return 0;
1499: if (ois) {
1501: if (*ois) {
1503: for (i = 0; i < n; i++) ISDestroy(&(*ois)[i]);
1504: PetscFree((*ois));
1505: }
1506: }
1507: if (iis) {
1509: if (*iis) {
1511: for (i = 0; i < n; i++) ISDestroy(&(*iis)[i]);
1512: PetscFree((*iis));
1513: }
1514: }
1515: return 0;
1516: }
1518: #define PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, xleft_loc, ylow_loc, xright_loc, yhigh_loc, n) \
1519: { \
1520: PetscInt first_row = first / M, last_row = last / M + 1; \
1521: /* \
1522: Compute ylow_loc and yhigh_loc so that (ylow_loc,xleft) and (yhigh_loc,xright) are the corners \
1523: of the bounding box of the intersection of the subdomain with the local ownership range (local \
1524: subdomain). \
1525: Also compute xleft_loc and xright_loc as the lower and upper bounds on the first and last rows \
1526: of the intersection. \
1527: */ \
1528: /* ylow_loc is the grid row containing the first element of the local sumbdomain */ \
1529: *ylow_loc = PetscMax(first_row, ylow); \
1530: /* xleft_loc is the offset of first element of the local subdomain within its grid row (might actually be outside the local subdomain) */ \
1531: *xleft_loc = *ylow_loc == first_row ? PetscMax(first % M, xleft) : xleft; \
1532: /* yhigh_loc is the grid row above the last local subdomain element */ \
1533: *yhigh_loc = PetscMin(last_row, yhigh); \
1534: /* xright is the offset of the end of the local subdomain within its grid row (might actually be outside the local subdomain) */ \
1535: *xright_loc = *yhigh_loc == last_row ? PetscMin(xright, last % M) : xright; \
1536: /* Now compute the size of the local subdomain n. */ \
1537: *n = 0; \
1538: if (*ylow_loc < *yhigh_loc) { \
1539: PetscInt width = xright - xleft; \
1540: *n += width * (*yhigh_loc - *ylow_loc - 1); \
1541: *n += PetscMin(PetscMax(*xright_loc - xleft, 0), width); \
1542: *n -= PetscMin(PetscMax(*xleft_loc - xleft, 0), width); \
1543: } \
1544: }
1546: /*@
1547: PCGASMCreateSubdomains2D - Creates the index sets for the `PCGASM` overlapping Schwarz
1548: preconditioner for a two-dimensional problem on a regular grid.
1550: Collective
1552: Input Parameters:
1553: + pc - the preconditioner context
1554: . M - the global number of grid points in the x direction
1555: . N - the global number of grid points in the y direction
1556: . Mdomains - the global number of subdomains in the x direction
1557: . Ndomains - the global number of subdomains in the y direction
1558: . dof - degrees of freedom per node
1559: - overlap - overlap in mesh lines
1561: Output Parameters:
1562: + Nsub - the number of local subdomains created
1563: . iis - array of index sets defining inner (nonoverlapping) subdomains
1564: - ois - array of index sets defining outer (overlapping, if overlap > 0) subdomains
1566: Level: advanced
1568: .seealso: `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, `PCGASMSetOverlap()`
1569: @*/
1570: PetscErrorCode PCGASMCreateSubdomains2D(PC pc, PetscInt M, PetscInt N, PetscInt Mdomains, PetscInt Ndomains, PetscInt dof, PetscInt overlap, PetscInt *nsub, IS **iis, IS **ois)
1571: {
1572: PetscMPIInt size, rank;
1573: PetscInt i, j;
1574: PetscInt maxheight, maxwidth;
1575: PetscInt xstart, xleft, xright, xleft_loc, xright_loc;
1576: PetscInt ystart, ylow, yhigh, ylow_loc, yhigh_loc;
1577: PetscInt x[2][2], y[2][2], n[2];
1578: PetscInt first, last;
1579: PetscInt nidx, *idx;
1580: PetscInt ii, jj, s, q, d;
1581: PetscInt k, kk;
1582: PetscMPIInt color;
1583: MPI_Comm comm, subcomm;
1584: IS **xis = NULL, **is = ois, **is_local = iis;
1586: PetscObjectGetComm((PetscObject)pc, &comm);
1587: MPI_Comm_size(comm, &size);
1588: MPI_Comm_rank(comm, &rank);
1589: MatGetOwnershipRange(pc->pmat, &first, &last);
1591: "Matrix row partitioning unsuitable for domain decomposition: local row range (%" PetscInt_FMT ",%" PetscInt_FMT ") "
1592: "does not respect the number of degrees of freedom per grid point %" PetscInt_FMT,
1593: first, last, dof);
1595: /* Determine the number of domains with nonzero intersections with the local ownership range. */
1596: s = 0;
1597: ystart = 0;
1598: for (j = 0; j < Ndomains; ++j) {
1599: maxheight = N / Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */
1601: /* Vertical domain limits with an overlap. */
1602: ylow = PetscMax(ystart - overlap, 0);
1603: yhigh = PetscMin(ystart + maxheight + overlap, N);
1604: xstart = 0;
1605: for (i = 0; i < Mdomains; ++i) {
1606: maxwidth = M / Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */
1608: /* Horizontal domain limits with an overlap. */
1609: xleft = PetscMax(xstart - overlap, 0);
1610: xright = PetscMin(xstart + maxwidth + overlap, M);
1611: /*
1612: Determine whether this subdomain intersects this rank's ownership range of pc->pmat.
1613: */
1614: PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, (&xleft_loc), (&ylow_loc), (&xright_loc), (&yhigh_loc), (&nidx));
1615: if (nidx) ++s;
1616: xstart += maxwidth;
1617: } /* for (i = 0; i < Mdomains; ++i) */
1618: ystart += maxheight;
1619: } /* for (j = 0; j < Ndomains; ++j) */
1621: /* Now we can allocate the necessary number of ISs. */
1622: *nsub = s;
1623: PetscMalloc1(*nsub, is);
1624: PetscMalloc1(*nsub, is_local);
1625: s = 0;
1626: ystart = 0;
1627: for (j = 0; j < Ndomains; ++j) {
1628: maxheight = N / Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */
1630: /* Vertical domain limits with an overlap. */
1631: y[0][0] = PetscMax(ystart - overlap, 0);
1632: y[0][1] = PetscMin(ystart + maxheight + overlap, N);
1633: /* Vertical domain limits without an overlap. */
1634: y[1][0] = ystart;
1635: y[1][1] = PetscMin(ystart + maxheight, N);
1636: xstart = 0;
1637: for (i = 0; i < Mdomains; ++i) {
1638: maxwidth = M / Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */
1640: /* Horizontal domain limits with an overlap. */
1641: x[0][0] = PetscMax(xstart - overlap, 0);
1642: x[0][1] = PetscMin(xstart + maxwidth + overlap, M);
1643: /* Horizontal domain limits without an overlap. */
1644: x[1][0] = xstart;
1645: x[1][1] = PetscMin(xstart + maxwidth, M);
1646: /*
1647: Determine whether this domain intersects this rank's ownership range of pc->pmat.
1648: Do this twice: first for the domains with overlaps, and once without.
1649: During the first pass create the subcommunicators, and use them on the second pass as well.
1650: */
1651: for (q = 0; q < 2; ++q) {
1652: PetscBool split = PETSC_FALSE;
1653: /*
1654: domain limits, (xleft, xright) and (ylow, yheigh) are adjusted
1655: according to whether the domain with an overlap or without is considered.
1656: */
1657: xleft = x[q][0];
1658: xright = x[q][1];
1659: ylow = y[q][0];
1660: yhigh = y[q][1];
1661: PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, (&xleft_loc), (&ylow_loc), (&xright_loc), (&yhigh_loc), (&nidx));
1662: nidx *= dof;
1663: n[q] = nidx;
1664: /*
1665: Based on the counted number of indices in the local domain *with an overlap*,
1666: construct a subcommunicator of all the MPI ranks supporting this domain.
1667: Observe that a domain with an overlap might have nontrivial local support,
1668: while the domain without an overlap might not. Hence, the decision to participate
1669: in the subcommunicator must be based on the domain with an overlap.
1670: */
1671: if (q == 0) {
1672: if (nidx) color = 1;
1673: else color = MPI_UNDEFINED;
1674: MPI_Comm_split(comm, color, rank, &subcomm);
1675: split = PETSC_TRUE;
1676: }
1677: /*
1678: Proceed only if the number of local indices *with an overlap* is nonzero.
1679: */
1680: if (n[0]) {
1681: if (q == 0) xis = is;
1682: if (q == 1) {
1683: /*
1684: The IS for the no-overlap subdomain shares a communicator with the overlapping domain.
1685: Moreover, if the overlap is zero, the two ISs are identical.
1686: */
1687: if (overlap == 0) {
1688: (*is_local)[s] = (*is)[s];
1689: PetscObjectReference((PetscObject)(*is)[s]);
1690: continue;
1691: } else {
1692: xis = is_local;
1693: subcomm = ((PetscObject)(*is)[s])->comm;
1694: }
1695: } /* if (q == 1) */
1696: idx = NULL;
1697: PetscMalloc1(nidx, &idx);
1698: if (nidx) {
1699: k = 0;
1700: for (jj = ylow_loc; jj < yhigh_loc; ++jj) {
1701: PetscInt x0 = (jj == ylow_loc) ? xleft_loc : xleft;
1702: PetscInt x1 = (jj == yhigh_loc - 1) ? xright_loc : xright;
1703: kk = dof * (M * jj + x0);
1704: for (ii = x0; ii < x1; ++ii) {
1705: for (d = 0; d < dof; ++d) idx[k++] = kk++;
1706: }
1707: }
1708: }
1709: ISCreateGeneral(subcomm, nidx, idx, PETSC_OWN_POINTER, (*xis) + s);
1710: if (split) MPI_Comm_free(&subcomm);
1711: } /* if (n[0]) */
1712: } /* for (q = 0; q < 2; ++q) */
1713: if (n[0]) ++s;
1714: xstart += maxwidth;
1715: } /* for (i = 0; i < Mdomains; ++i) */
1716: ystart += maxheight;
1717: } /* for (j = 0; j < Ndomains; ++j) */
1718: return 0;
1719: }
1721: /*@C
1722: PCGASMGetSubdomains - Gets the subdomains supported on this MPI rank
1723: for the `PCGASM` additive Schwarz preconditioner.
1725: Not Collective
1727: Input Parameter:
1728: . pc - the preconditioner context
1730: Output Parameters:
1731: + n - the number of subdomains for this MPI rank (default value = 1)
1732: . iis - the index sets that define the inner subdomains (without overlap) supported on this rank (can be NULL)
1733: - ois - the index sets that define the outer subdomains (with overlap) supported on this rank (can be NULL)
1735: Note:
1736: The user is responsible for destroying the `IS`s and freeing the returned arrays.
1737: The `IS` numbering is in the parallel, global numbering of the vector.
1739: Level: advanced
1741: .seealso: `PCGASM`, `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`, `PCGASMCreateSubdomains2D()`,
1742: `PCGASMSetSubdomains()`, `PCGASMGetSubmatrices()`
1743: @*/
1744: PetscErrorCode PCGASMGetSubdomains(PC pc, PetscInt *n, IS *iis[], IS *ois[])
1745: {
1746: PC_GASM *osm;
1747: PetscBool match;
1748: PetscInt i;
1751: PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match);
1753: osm = (PC_GASM *)pc->data;
1754: if (n) *n = osm->n;
1755: if (iis) PetscMalloc1(osm->n, iis);
1756: if (ois) PetscMalloc1(osm->n, ois);
1757: if (iis || ois) {
1758: for (i = 0; i < osm->n; ++i) {
1759: if (iis) (*iis)[i] = osm->iis[i];
1760: if (ois) (*ois)[i] = osm->ois[i];
1761: }
1762: }
1763: return 0;
1764: }
1766: /*@C
1767: PCGASMGetSubmatrices - Gets the local submatrices (for this MPI rank
1768: only) for the `PCGASM` additive Schwarz preconditioner.
1770: Not Collective
1772: Input Parameter:
1773: . pc - the preconditioner context
1775: Output Parameters:
1776: + n - the number of matrices for this MPI rank (default value = 1)
1777: - mat - the matrices
1779: Note:
1780: Matrices returned by this routine have the same communicators as the index sets (IS)
1781: used to define subdomains in `PCGASMSetSubdomains()`
1783: Level: advanced
1785: .seealso: `PCGASM`, `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`,
1786: `PCGASMCreateSubdomains2D()`, `PCGASMSetSubdomains()`, `PCGASMGetSubdomains()`
1787: @*/
1788: PetscErrorCode PCGASMGetSubmatrices(PC pc, PetscInt *n, Mat *mat[])
1789: {
1790: PC_GASM *osm;
1791: PetscBool match;
1797: PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match);
1799: osm = (PC_GASM *)pc->data;
1800: if (n) *n = osm->n;
1801: if (mat) *mat = osm->pmat;
1802: return 0;
1803: }
1805: /*@
1806: PCGASMSetUseDMSubdomains - Indicates whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible for `PCGASM`
1808: Logically Collective
1810: Input Parameters:
1811: + pc - the preconditioner
1812: - flg - boolean indicating whether to use subdomains defined by the `DM`
1814: Options Database Key:
1815: . -pc_gasm_dm_subdomains -pc_gasm_overlap -pc_gasm_total_subdomains
1817: Level: intermediate
1819: Note:
1820: `PCGASMSetSubdomains()`, `PCGASMSetTotalSubdomains()` or `PCGASMSetOverlap()` take precedence over `PCGASMSetUseDMSubdomains()`,
1821: so setting `PCGASMSetSubdomains()` with nontrivial subdomain ISs or any of `PCGASMSetTotalSubdomains()` and `PCGASMSetOverlap()`
1822: automatically turns the latter off.
1824: .seealso: `PCGASM`, `PCGASMGetUseDMSubdomains()`, `PCGASMSetSubdomains()`, `PCGASMSetOverlap()`
1825: `PCGASMCreateSubdomains2D()`
1826: @*/
1827: PetscErrorCode PCGASMSetUseDMSubdomains(PC pc, PetscBool flg)
1828: {
1829: PC_GASM *osm = (PC_GASM *)pc->data;
1830: PetscBool match;
1835: PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match);
1836: if (match) {
1837: if (!osm->user_subdomains && osm->N == PETSC_DETERMINE && osm->overlap < 0) osm->dm_subdomains = flg;
1838: }
1839: return 0;
1840: }
1842: /*@
1843: PCGASMGetUseDMSubdomains - Returns flag indicating whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible with `PCGASM`
1845: Not Collective
1847: Input Parameter:
1848: . pc - the preconditioner
1850: Output Parameter:
1851: . flg - boolean indicating whether to use subdomains defined by the `DM`
1853: Level: intermediate
1855: .seealso: `PCGASM`, `PCGASMSetUseDMSubdomains()`, `PCGASMSetOverlap()`
1856: `PCGASMCreateSubdomains2D()`
1857: @*/
1858: PetscErrorCode PCGASMGetUseDMSubdomains(PC pc, PetscBool *flg)
1859: {
1860: PC_GASM *osm = (PC_GASM *)pc->data;
1861: PetscBool match;
1865: PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match);
1866: if (match) {
1867: if (flg) *flg = osm->dm_subdomains;
1868: }
1869: return 0;
1870: }