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 process may intersect multiple subdomains and any subdomain may
4: intersect multiple MPI processes. Intersections of subdomains with MPI processes 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 process (set in `PCGASMSetSubdomains()` or calculated in `PCGASMSetTotalSubdomains()`)
9: nmax - maximum number of local subdomains per process (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: PetscFunctionBegin;
45: /* Determine the number of globally-distinct subdomains and compute a global numbering for them. */
46: PetscCall(PetscMalloc2(osm->n, numbering, osm->n, permutation));
47: PetscCall(PetscObjectsListGetGlobalNumbering(PetscObjectComm((PetscObject)pc), osm->n, (PetscObject *)osm->iis, NULL, *numbering));
48: for (i = 0; i < osm->n; ++i) (*permutation)[i] = i;
49: PetscCall(PetscSortIntWithPermutation(osm->n, *numbering, *permutation));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: static PetscErrorCode PCGASMSubdomainView_Private(PC pc, PetscInt i, PetscViewer viewer)
54: {
55: PC_GASM *osm = (PC_GASM *)pc->data;
56: PetscInt nidx;
57: const PetscInt *idx;
58: PetscViewer sviewer;
59: char *cidx;
61: PetscFunctionBegin;
62: PetscCheck(i >= -1 && i < osm->n, PetscObjectComm((PetscObject)viewer), PETSC_ERR_ARG_WRONG, "Invalid subdomain %" PetscInt_FMT ": must nonnegative and less than %" PetscInt_FMT, i, osm->n);
64: /* Inner subdomains. */
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: PetscCall(PetscViewerASCIIPrintf(viewer, "Inner subdomain:\n"));
72: PetscCall(PetscViewerFlush(viewer));
73: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
74: if (i > -1) {
75: PetscCall(ISGetLocalSize(osm->iis[i], &nidx));
76: PetscCall(PetscMalloc1(16 * (nidx + 1) + 1, &cidx));
77: PetscCall(PetscViewerStringOpen(PETSC_COMM_SELF, cidx, 16 * (nidx + 1) + 1, &sviewer));
78: PetscCall(ISGetIndices(osm->iis[i], &idx));
79: for (PetscInt j = 0; j < nidx; ++j) PetscCall(PetscViewerStringSPrintf(sviewer, "%" PetscInt_FMT " ", idx[j]));
80: PetscCall(ISRestoreIndices(osm->iis[i], &idx));
81: PetscCall(PetscViewerDestroy(&sviewer));
82: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%s", cidx));
83: PetscCall(PetscFree(cidx));
84: }
85: PetscCall(PetscViewerFlush(viewer));
86: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
87: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
88: PetscCall(PetscViewerFlush(viewer));
90: /* Outer subdomains. */
91: /*
92: No more than 15 characters per index plus a space.
93: PetscViewerStringSPrintf requires a string of size at least 2, so use (nidx+1) instead of nidx,
94: in case nidx == 0. That will take care of the space for the trailing '\0' as well.
95: For nidx == 0, the whole string 16 '\0'.
96: */
97: PetscCall(PetscViewerASCIIPrintf(viewer, "Outer subdomain:\n"));
98: PetscCall(PetscViewerFlush(viewer));
99: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
100: if (i > -1) {
101: PetscCall(ISGetLocalSize(osm->ois[i], &nidx));
102: PetscCall(PetscMalloc1(16 * (nidx + 1) + 1, &cidx));
103: PetscCall(PetscViewerStringOpen(PETSC_COMM_SELF, cidx, 16 * (nidx + 1) + 1, &sviewer));
104: PetscCall(ISGetIndices(osm->ois[i], &idx));
105: for (PetscInt j = 0; j < nidx; ++j) PetscCall(PetscViewerStringSPrintf(sviewer, "%" PetscInt_FMT " ", idx[j]));
106: PetscCall(PetscViewerDestroy(&sviewer));
107: PetscCall(ISRestoreIndices(osm->ois[i], &idx));
108: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%s", cidx));
109: PetscCall(PetscFree(cidx));
110: }
111: PetscCall(PetscViewerFlush(viewer));
112: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
113: PetscCall(PetscViewerASCIIPrintf(viewer, "\n"));
114: PetscCall(PetscViewerFlush(viewer));
115: PetscFunctionReturn(PETSC_SUCCESS);
116: }
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 found;
125: PetscViewer viewer;
126: PetscInt *numbering, *permutation; /* global numbering of locally-supported subdomains and the permutation from the local ordering */
128: PetscFunctionBegin;
129: PetscCall(PCGetOptionsPrefix(pc, &prefix));
130: PetscCall(PetscOptionsHasName(NULL, prefix, "-pc_gasm_print_subdomains", &found));
131: if (!found) PetscFunctionReturn(PETSC_SUCCESS);
132: PetscCall(PetscOptionsGetString(NULL, prefix, "-pc_gasm_print_subdomains", fname, sizeof(fname), &found));
133: if (!found) PetscCall(PetscStrncpy(fname, "stdout", sizeof(fname)));
134: PetscCall(PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc), fname, &viewer));
135: /*
136: Make sure the viewer has a name. Otherwise this may cause a deadlock or other weird errors when creating a subcomm viewer:
137: the subcomm viewer will attempt to inherit the viewer's name, which, if not set, will be constructed collectively on the comm.
138: */
139: PetscCall(PetscObjectName((PetscObject)viewer));
140: l = 0;
141: PetscCall(PCGASMComputeGlobalSubdomainNumbering_Private(pc, &numbering, &permutation));
142: for (count = 0; count < osm->N; ++count) {
143: /* Now let subdomains go one at a time in the global numbering order and print their subdomain/solver info. */
144: if (l < osm->n) {
145: d = permutation[l]; /* d is the local number of the l-th smallest (in the global ordering) among the locally supported subdomains */
146: if (numbering[d] == count) l++;
147: else d = -1;
148: } else d = -1;
149: PetscCall(PCGASMSubdomainView_Private(pc, d, viewer));
150: }
151: PetscCall(PetscFree2(numbering, permutation));
152: PetscCall(PetscViewerDestroy(&viewer));
153: PetscFunctionReturn(PETSC_SUCCESS);
154: }
156: static PetscErrorCode PCView_GASM(PC pc, PetscViewer viewer)
157: {
158: PC_GASM *osm = (PC_GASM *)pc->data;
159: const char *prefix;
160: PetscMPIInt rank, size;
161: PetscInt bsz;
162: PetscBool isascii, view_subdomains = PETSC_FALSE;
163: PetscViewer sviewer;
164: PetscInt l;
165: char overlap[256] = "user-defined overlap";
166: char gsubdomains[256] = "unknown total number of subdomains";
167: char msubdomains[256] = "unknown max number of local subdomains";
168: PetscInt *numbering, *permutation; /* global numbering of locally-supported subdomains and the permutation from the local ordering */
170: PetscFunctionBegin;
171: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
172: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
174: if (osm->overlap >= 0) PetscCall(PetscSNPrintf(overlap, sizeof(overlap), "requested amount of overlap = %" PetscInt_FMT, osm->overlap));
175: if (osm->N != PETSC_DETERMINE) PetscCall(PetscSNPrintf(gsubdomains, sizeof(gsubdomains), "total number of subdomains = %" PetscInt_FMT, osm->N));
176: if (osm->nmax != PETSC_DETERMINE) PetscCall(PetscSNPrintf(msubdomains, sizeof(msubdomains), "max number of local subdomains = %" PetscInt_FMT, osm->nmax));
178: PetscCall(PCGetOptionsPrefix(pc, &prefix));
179: PetscCall(PetscOptionsGetBool(NULL, prefix, "-pc_gasm_view_subdomains", &view_subdomains, NULL));
181: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
182: if (isascii) {
183: /*
184: Make sure the viewer has a name. Otherwise this may cause a deadlock when creating a subcomm viewer:
185: the subcomm viewer will attempt to inherit the viewer's name, which, if not set, will be constructed
186: collectively on the comm.
187: */
188: PetscCall(PetscObjectName((PetscObject)viewer));
189: PetscCall(PetscViewerASCIIPrintf(viewer, " Restriction/interpolation type: %s\n", PCGASMTypes[osm->type]));
190: PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", overlap));
191: PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", gsubdomains));
192: PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", msubdomains));
193: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
194: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " [%d|%d] number of locally-supported subdomains = %" PetscInt_FMT "\n", rank, size, osm->n));
195: PetscCall(PetscViewerFlush(viewer));
196: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
197: /* Cannot take advantage of osm->same_subdomain_solvers without a global numbering of subdomains. */
198: PetscCall(PetscViewerASCIIPrintf(viewer, " Subdomain solver info:\n"));
199: PetscCall(PetscViewerASCIIPushTab(viewer));
200: PetscCall(PetscViewerASCIIPrintf(viewer, " - - - - - - - - - - - - - - - - - -\n"));
201: /* Now let subdomains go one at a time in the global numbering order and print their subdomain/solver info. */
202: PetscCall(PCGASMComputeGlobalSubdomainNumbering_Private(pc, &numbering, &permutation));
203: l = 0;
204: for (PetscInt 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: PetscCallMPI(MPI_Comm_size(((PetscObject)osm->ois[d])->comm, &ssize));
210: PetscCallMPI(MPI_Comm_rank(((PetscObject)osm->ois[d])->comm, &srank));
211: PetscCall(PetscViewerGetSubViewer(viewer, ((PetscObject)osm->ois[d])->comm, &sviewer));
212: PetscCall(ISGetLocalSize(osm->ois[d], &bsz));
213: PetscCall(PetscViewerASCIISynchronizedPrintf(sviewer, " [%d|%d] (subcomm [%d|%d]) local subdomain number %" PetscInt_FMT ", local size = %" PetscInt_FMT "\n", rank, size, srank, ssize, d, bsz));
214: PetscCall(PetscViewerFlush(sviewer));
215: PetscCall(PetscViewerASCIIPushTab(sviewer));
216: if (view_subdomains) PetscCall(PCGASMSubdomainView_Private(pc, d, sviewer));
217: if (!pc->setupcalled) {
218: PetscCall(PetscViewerASCIISynchronizedPrintf(sviewer, " Solver not set up yet: PCSetUp() not yet called\n"));
219: } else {
220: PetscCall(KSPView(osm->ksp[d], sviewer));
221: }
222: PetscCall(PetscViewerASCIIPopTab(sviewer));
223: PetscCall(PetscViewerASCIIPrintf(sviewer, " - - - - - - - - - - - - - - - - - -\n"));
224: PetscCall(PetscViewerFlush(sviewer));
225: PetscCall(PetscViewerRestoreSubViewer(viewer, ((PetscObject)osm->ois[d])->comm, &sviewer));
226: ++l;
227: } else {
228: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
229: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
230: }
231: } else {
232: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
233: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
234: }
235: }
236: PetscCall(PetscFree2(numbering, permutation));
237: PetscCall(PetscViewerASCIIPopTab(viewer));
238: PetscCall(PetscViewerFlush(viewer));
239: /* this line is needed to match the extra PetscViewerASCIIPushSynchronized() in PetscViewerGetSubViewer() */
240: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
241: }
242: PetscFunctionReturn(PETSC_SUCCESS);
243: }
245: PETSC_INTERN PetscErrorCode PCGASMCreateLocalSubdomains(Mat A, PetscInt nloc, IS *iis[]);
247: static PetscErrorCode PCGASMSetHierarchicalPartitioning(PC pc)
248: {
249: PC_GASM *osm = (PC_GASM *)pc->data;
250: MatPartitioning part;
251: MPI_Comm comm;
252: PetscMPIInt size;
253: PetscInt nlocalsubdomains, fromrows_localsize;
254: IS partitioning, fromrows, isn;
255: Vec outervec;
257: PetscFunctionBegin;
258: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
259: PetscCallMPI(MPI_Comm_size(comm, &size));
260: /* we do not need a hierarchical partitioning when
261: * the total number of subdomains is consistent with
262: * the number of MPI tasks.
263: * For the following cases, we do not need to use HP
264: * */
265: if (osm->N == PETSC_DETERMINE || osm->N >= size || osm->N == 1) PetscFunctionReturn(PETSC_SUCCESS);
266: PetscCheck(size % osm->N == 0, PETSC_COMM_WORLD, PETSC_ERR_ARG_INCOMP, "have to specify the total number of subdomains %" PetscInt_FMT " to be a factor of the number of ranks %d ", osm->N, size);
267: nlocalsubdomains = size / osm->N;
268: osm->n = 1;
269: PetscCall(MatPartitioningCreate(comm, &part));
270: PetscCall(MatPartitioningSetAdjacency(part, pc->pmat));
271: PetscCall(MatPartitioningSetType(part, MATPARTITIONINGHIERARCH));
272: PetscCall(MatPartitioningHierarchicalSetNcoarseparts(part, osm->N));
273: PetscCall(MatPartitioningHierarchicalSetNfineparts(part, nlocalsubdomains));
274: PetscCall(MatPartitioningSetFromOptions(part));
275: /* get new rank owner number of each vertex */
276: PetscCall(MatPartitioningApply(part, &partitioning));
277: PetscCall(ISBuildTwoSided(partitioning, NULL, &fromrows));
278: PetscCall(ISPartitioningToNumbering(partitioning, &isn));
279: PetscCall(ISDestroy(&isn));
280: PetscCall(ISGetLocalSize(fromrows, &fromrows_localsize));
281: PetscCall(MatPartitioningDestroy(&part));
282: PetscCall(MatCreateVecs(pc->pmat, &outervec, NULL));
283: PetscCall(VecCreateMPI(comm, fromrows_localsize, PETSC_DETERMINE, &osm->pcx));
284: PetscCall(VecDuplicate(osm->pcx, &osm->pcy));
285: PetscCall(VecScatterCreate(osm->pcx, NULL, outervec, fromrows, &osm->pctoouter));
286: PetscCall(MatCreateSubMatrix(pc->pmat, fromrows, fromrows, MAT_INITIAL_MATRIX, &osm->permutationP));
287: PetscCall(PetscObjectReference((PetscObject)fromrows));
288: osm->permutationIS = fromrows;
289: osm->pcmat = pc->pmat;
290: PetscCall(PetscObjectReference((PetscObject)osm->permutationP));
291: pc->pmat = osm->permutationP;
292: PetscCall(VecDestroy(&outervec));
293: PetscCall(ISDestroy(&fromrows));
294: PetscCall(ISDestroy(&partitioning));
295: osm->n = PETSC_DETERMINE;
296: PetscFunctionReturn(PETSC_SUCCESS);
297: }
299: static PetscErrorCode PCSetUp_GASM(PC pc)
300: {
301: PC_GASM *osm = (PC_GASM *)pc->data;
302: PetscInt nInnerIndices, nTotalInnerIndices;
303: PetscMPIInt rank, size;
304: MatReuse scall = MAT_REUSE_MATRIX;
305: KSP ksp;
306: PC subpc;
307: const char *prefix, *pprefix;
308: Vec x, y;
309: PetscInt oni; /* Number of indices in the i-th local outer subdomain. */
310: const PetscInt *oidxi; /* Indices from the i-th subdomain local outer subdomain. */
311: PetscInt on; /* Number of indices in the disjoint union of local outer subdomains. */
312: PetscInt *oidx; /* Indices in the disjoint union of local outer subdomains. */
313: IS gois; /* Disjoint union the global indices of outer subdomains. */
314: IS goid; /* Identity IS of the size of the disjoint union of outer subdomains. */
315: PetscInt gostart; /* Start of locally-owned indices in the vectors -- osm->gx,osm->gy -- over the disjoint union of outer subdomains. */
316: PetscInt num_subdomains = 0;
317: DM *subdomain_dm = NULL;
318: char **subdomain_names = NULL;
319: PetscInt *numbering;
321: PetscFunctionBegin;
322: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
323: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
324: if (!pc->setupcalled) {
325: /* use a hierarchical partitioning */
326: if (osm->hierarchicalpartitioning) PetscCall(PCGASMSetHierarchicalPartitioning(pc));
327: if (osm->n == PETSC_DETERMINE) {
328: if (osm->N != PETSC_DETERMINE) {
329: /* No local subdomains given, but the desired number of total subdomains is known, so construct them accordingly. */
330: PetscCall(PCGASMCreateSubdomains(pc->pmat, osm->N, &osm->n, &osm->iis));
331: } else if (osm->dm_subdomains && pc->dm) {
332: /* try pc->dm next, if allowed */
333: IS *inner_subdomain_is, *outer_subdomain_is;
334: PetscCall(DMCreateDomainDecomposition(pc->dm, &num_subdomains, &subdomain_names, &inner_subdomain_is, &outer_subdomain_is, &subdomain_dm));
335: if (num_subdomains) PetscCall(PCGASMSetSubdomains(pc, num_subdomains, inner_subdomain_is, outer_subdomain_is));
336: for (PetscInt d = 0; d < num_subdomains; ++d) {
337: if (inner_subdomain_is) PetscCall(ISDestroy(&inner_subdomain_is[d]));
338: if (outer_subdomain_is) PetscCall(ISDestroy(&outer_subdomain_is[d]));
339: }
340: PetscCall(PetscFree(inner_subdomain_is));
341: PetscCall(PetscFree(outer_subdomain_is));
342: } else {
343: /* still no subdomains; use one per rank */
344: osm->nmax = osm->n = 1;
345: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
346: osm->N = size;
347: PetscCall(PCGASMCreateLocalSubdomains(pc->pmat, osm->n, &osm->iis));
348: }
349: }
350: if (!osm->iis) {
351: /*
352: osm->n was set in PCGASMSetSubdomains(), but the actual subdomains have not been supplied.
353: We create the requisite number of local inner subdomains and then expand them into
354: out subdomains, if necessary.
355: */
356: PetscCall(PCGASMCreateLocalSubdomains(pc->pmat, osm->n, &osm->iis));
357: }
358: if (!osm->ois) {
359: /*
360: Initially make outer subdomains the same as inner subdomains. If nonzero additional overlap
361: has been requested, copy the inner subdomains over so they can be modified.
362: */
363: PetscCall(PetscMalloc1(osm->n, &osm->ois));
364: for (PetscInt i = 0; i < osm->n; ++i) {
365: if (osm->overlap > 0 && osm->N > 1) { /* With positive overlap, osm->iis[i] will be modified */
366: PetscCall(ISDuplicate(osm->iis[i], (osm->ois) + i));
367: PetscCall(ISCopy(osm->iis[i], osm->ois[i]));
368: } else {
369: PetscCall(PetscObjectReference((PetscObject)osm->iis[i]));
370: osm->ois[i] = osm->iis[i];
371: }
372: }
373: if (osm->overlap > 0 && osm->N > 1) {
374: /* Extend the "overlapping" regions by a number of steps */
375: PetscCall(MatIncreaseOverlapSplit(pc->pmat, osm->n, osm->ois, osm->overlap));
376: }
377: }
379: /* Now the subdomains are defined. Determine their global and max local numbers, if necessary. */
380: if (osm->nmax == PETSC_DETERMINE) {
381: PetscInt outwork;
382: /* determine global number of subdomains and the max number of local subdomains */
383: outwork = osm->n;
384: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &outwork, 1, MPIU_INT, MPI_MAX, PetscObjectComm((PetscObject)pc)));
385: osm->nmax = outwork;
386: }
387: if (osm->N == PETSC_DETERMINE) {
388: /* Determine the number of globally-distinct subdomains and compute a global numbering for them. */
389: PetscCall(PetscObjectsListGetGlobalNumbering(PetscObjectComm((PetscObject)pc), osm->n, (PetscObject *)osm->ois, &osm->N, NULL));
390: }
392: if (osm->sort_indices) {
393: for (PetscInt i = 0; i < osm->n; i++) {
394: PetscCall(ISSort(osm->ois[i]));
395: PetscCall(ISSort(osm->iis[i]));
396: }
397: }
398: PetscCall(PCGetOptionsPrefix(pc, &prefix));
399: PetscCall(PCGASMPrintSubdomains(pc));
401: /*
402: Merge the ISs, create merged vectors and restrictions.
403: */
404: /* Merge outer subdomain ISs and construct a restriction onto the disjoint union of local outer subdomains. */
405: on = 0;
406: for (PetscInt i = 0; i < osm->n; i++) {
407: PetscCall(ISGetLocalSize(osm->ois[i], &oni));
408: on += oni;
409: }
410: PetscCall(PetscMalloc1(on, &oidx));
411: on = 0;
412: /* Merge local indices together */
413: for (PetscInt i = 0; i < osm->n; i++) {
414: PetscCall(ISGetLocalSize(osm->ois[i], &oni));
415: PetscCall(ISGetIndices(osm->ois[i], &oidxi));
416: PetscCall(PetscArraycpy(oidx + on, oidxi, oni));
417: PetscCall(ISRestoreIndices(osm->ois[i], &oidxi));
418: on += oni;
419: }
420: PetscCall(ISCreateGeneral(((PetscObject)pc)->comm, on, oidx, PETSC_OWN_POINTER, &gois));
421: nTotalInnerIndices = 0;
422: for (PetscInt i = 0; i < osm->n; i++) {
423: PetscCall(ISGetLocalSize(osm->iis[i], &nInnerIndices));
424: nTotalInnerIndices += nInnerIndices;
425: }
426: PetscCall(VecCreateMPI(((PetscObject)pc)->comm, nTotalInnerIndices, PETSC_DETERMINE, &x));
427: PetscCall(VecDuplicate(x, &y));
429: /* Keep the merged vectors on the host so that their arrays also work with host-only subdomain solvers. */
430: PetscCall(VecCreateMPI(PetscObjectComm((PetscObject)pc), on, PETSC_DECIDE, &osm->gx));
431: PetscCall(VecDuplicate(osm->gx, &osm->gy));
432: PetscCall(VecGetOwnershipRange(osm->gx, &gostart, NULL));
433: PetscCall(ISCreateStride(PetscObjectComm((PetscObject)pc), on, gostart, 1, &goid));
434: /* gois might indices not on local */
435: PetscCall(VecScatterCreate(x, gois, osm->gx, goid, &osm->gorestriction));
436: PetscCall(PetscMalloc1(osm->n, &numbering));
437: PetscCall(PetscObjectsListGetGlobalNumbering(PetscObjectComm((PetscObject)pc), osm->n, (PetscObject *)osm->ois, NULL, numbering));
438: PetscCall(VecDestroy(&x));
439: PetscCall(ISDestroy(&gois));
441: /* Merge inner subdomain ISs and construct a restriction onto the disjoint union of local inner subdomains. */
442: {
443: PetscInt ini; /* Number of indices the i-th a local inner subdomain. */
444: PetscInt in; /* Number of indices in the disjoint union of local inner subdomains. */
445: PetscInt *iidx; /* Global indices in the merged local inner subdomain. */
446: PetscInt *ioidx; /* Global indices of the disjoint union of inner subdomains within the disjoint union of outer subdomains. */
447: IS giis; /* IS for the disjoint union of inner subdomains. */
448: IS giois; /* IS for the disjoint union of inner subdomains within the disjoint union of outer subdomains. */
449: PetscScalar *array;
450: const PetscInt *indices;
451: on = 0;
452: for (PetscInt i = 0; i < osm->n; i++) {
453: PetscCall(ISGetLocalSize(osm->ois[i], &oni));
454: on += oni;
455: }
456: PetscCall(PetscMalloc1(on, &iidx));
457: PetscCall(PetscMalloc1(on, &ioidx));
458: PetscCall(VecGetArray(y, &array));
459: /* set communicator id to determine where overlap is */
460: in = 0;
461: for (PetscInt i = 0; i < osm->n; i++) {
462: PetscCall(ISGetLocalSize(osm->iis[i], &ini));
463: for (PetscInt k = 0; k < ini; ++k) array[in + k] = numbering[i];
464: in += ini;
465: }
466: PetscCall(VecRestoreArray(y, &array));
467: PetscCall(VecScatterBegin(osm->gorestriction, y, osm->gy, INSERT_VALUES, SCATTER_FORWARD));
468: PetscCall(VecScatterEnd(osm->gorestriction, y, osm->gy, INSERT_VALUES, SCATTER_FORWARD));
469: PetscCall(VecGetOwnershipRange(osm->gy, &gostart, NULL));
470: PetscCall(VecGetArray(osm->gy, &array));
471: on = 0;
472: in = 0;
473: for (PetscInt i = 0; i < osm->n; i++) {
474: PetscCall(ISGetLocalSize(osm->ois[i], &oni));
475: PetscCall(ISGetIndices(osm->ois[i], &indices));
476: for (PetscInt k = 0; k < oni; k++) {
477: /* skip overlapping indices to get inner domain */
478: if (PetscRealPart(array[on + k]) != numbering[i]) continue;
479: iidx[in] = indices[k];
480: ioidx[in++] = gostart + on + k;
481: }
482: PetscCall(ISRestoreIndices(osm->ois[i], &indices));
483: on += oni;
484: }
485: PetscCall(VecRestoreArray(osm->gy, &array));
486: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), in, iidx, PETSC_OWN_POINTER, &giis));
487: PetscCall(ISCreateGeneral(PetscObjectComm((PetscObject)pc), in, ioidx, PETSC_OWN_POINTER, &giois));
488: PetscCall(VecScatterCreate(y, giis, osm->gy, giois, &osm->girestriction));
489: PetscCall(VecDestroy(&y));
490: PetscCall(ISDestroy(&giis));
491: PetscCall(ISDestroy(&giois));
492: }
493: PetscCall(ISDestroy(&goid));
494: PetscCall(PetscFree(numbering));
496: /* Create the subdomain solvers */
497: PetscCall(PetscMalloc1(osm->n, &osm->ksp));
498: for (PetscInt i = 0; i < osm->n; i++) {
499: char subprefix[PETSC_MAX_PATH_LEN + 1];
500: PetscCall(KSPCreate(((PetscObject)osm->ois[i])->comm, &ksp));
501: PetscCall(KSPSetNestLevel(ksp, pc->kspnestlevel));
502: PetscCall(KSPSetErrorIfNotConverged(ksp, pc->erroriffailure));
503: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1));
504: PetscCall(KSPSetType(ksp, KSPPREONLY));
505: PetscCall(KSPGetPC(ksp, &subpc)); /* Why do we need this here? */
506: if (subdomain_dm) {
507: PetscCall(KSPSetDM(ksp, subdomain_dm[i]));
508: PetscCall(DMDestroy(subdomain_dm + i));
509: }
510: PetscCall(PCGetOptionsPrefix(pc, &prefix));
511: PetscCall(KSPSetOptionsPrefix(ksp, prefix));
512: if (subdomain_names && subdomain_names[i]) {
513: PetscCall(PetscSNPrintf(subprefix, PETSC_MAX_PATH_LEN, "sub_%s_", subdomain_names[i]));
514: PetscCall(KSPAppendOptionsPrefix(ksp, subprefix));
515: PetscCall(PetscFree(subdomain_names[i]));
516: }
517: PetscCall(KSPAppendOptionsPrefix(ksp, "sub_"));
518: osm->ksp[i] = ksp;
519: }
520: PetscCall(PetscFree(subdomain_dm));
521: PetscCall(PetscFree(subdomain_names));
522: scall = MAT_INITIAL_MATRIX;
523: } else { /* if (pc->setupcalled) */
524: /*
525: Destroy the submatrices from the previous iteration
526: */
527: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
528: PetscCall(MatDestroyMatrices(osm->n, &osm->pmat));
529: scall = MAT_INITIAL_MATRIX;
530: }
531: if (osm->permutationIS) {
532: PetscCall(MatCreateSubMatrix(pc->pmat, osm->permutationIS, osm->permutationIS, scall, &osm->permutationP));
533: PetscCall(PetscObjectReference((PetscObject)osm->permutationP));
534: osm->pcmat = pc->pmat;
535: pc->pmat = osm->permutationP;
536: }
537: }
539: /*
540: Extract the submatrices.
541: */
542: if (size > 1) PetscCall(MatCreateSubMatricesMPI(pc->pmat, osm->n, osm->ois, osm->ois, scall, &osm->pmat));
543: else PetscCall(MatCreateSubMatrices(pc->pmat, osm->n, osm->ois, osm->ois, scall, &osm->pmat));
544: if (scall == MAT_INITIAL_MATRIX) {
545: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc->pmat, &pprefix));
546: for (PetscInt i = 0; i < osm->n; i++) PetscCall(PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i], pprefix));
547: }
549: /* Return control to the user so that the submatrices can be modified (e.g., to apply
550: different boundary conditions for the submatrices than for the global problem) */
551: PetscCall(PCModifySubMatrices(pc, osm->n, osm->ois, osm->ois, osm->pmat, pc->modifysubmatricesP));
553: /*
554: Loop over submatrices putting them into local ksps
555: */
556: for (PetscInt i = 0; i < osm->n; i++) {
557: PetscCall(KSPSetOperators(osm->ksp[i], osm->pmat[i], osm->pmat[i]));
558: PetscCall(KSPGetOptionsPrefix(osm->ksp[i], &prefix));
559: PetscCall(MatSetOptionsPrefix(osm->pmat[i], prefix));
560: if (!pc->setupcalled) PetscCall(KSPSetFromOptions(osm->ksp[i]));
561: }
562: if (!pc->setupcalled) {
563: /* MatCreateVecs() would allocate arrays before they are needed. Create arrayless work vectors where supported because
564: PCApply_GASM() and PCApplyTranspose_GASM() place the merged vectors' arrays into them. */
565: PetscCall(PetscMalloc1(osm->n, &osm->x));
566: PetscCall(PetscMalloc1(osm->n, &osm->y));
567: for (PetscInt i = 0; i < osm->n; ++i) {
568: VecType vtype;
569: PetscInt oNi, m, n, M, N;
570: PetscMPIInt subsize;
571: PetscBool boundtocpu, bindingpropagates, iskok;
573: PetscCall(MatGetVecType(osm->pmat[i], &vtype));
574: PetscCall(PetscStrcmpAny(vtype, &iskok, VECKOKKOS, VECSEQKOKKOS, VECMPIKOKKOS, ""));
575: PetscCall(ISGetLocalSize(osm->ois[i], &oni));
576: PetscCall(ISGetSize(osm->ois[i], &oNi));
577: PetscCall(MatGetLocalSize(osm->pmat[i], &m, &n));
578: PetscCall(MatGetSize(osm->pmat[i], &M, &N));
579: PetscCheck(m == oni && n == oni && M == oNi && N == oNi, PetscObjectComm((PetscObject)osm->pmat[i]), PETSC_ERR_ARG_SIZ, "Modified submatrix %" PetscInt_FMT " has size %" PetscInt_FMT " x %" PetscInt_FMT " (local %" PetscInt_FMT " x %" PetscInt_FMT "), but its outer index set has size %" PetscInt_FMT " (local %" PetscInt_FMT ")", i, M, N, m, n, oNi, oni);
580: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)osm->pmat[i]), &subsize));
581: PetscCall(MatBoundToCPU(osm->pmat[i], &boundtocpu));
582: PetscCall(MatGetBindingPropagates(osm->pmat[i], &bindingpropagates));
583: for (PetscInt j = 0; j < 2; ++j) {
584: Vec *v = j ? &osm->y[i] : &osm->x[i];
586: if (iskok) {
587: /* Kokkos vectors require valid host storage even when VecPlaceArray() will replace it. */
588: PetscCall(VecCreate(PetscObjectComm((PetscObject)osm->pmat[i]), v));
589: PetscCall(VecSetSizes(*v, oni, oNi));
590: PetscCall(VecSetBlockSize(*v, 1));
591: } else if (subsize == 1) PetscCall(VecCreateSeqWithArray(PetscObjectComm((PetscObject)osm->pmat[i]), 1, oni, NULL, v));
592: else PetscCall(VecCreateMPIWithArray(PetscObjectComm((PetscObject)osm->pmat[i]), 1, oni, oNi, NULL, v));
593: PetscCall(VecSetType(*v, vtype));
594: if (boundtocpu && bindingpropagates) {
595: PetscCall(VecSetBindingPropagates(*v, PETSC_TRUE));
596: PetscCall(VecBindToCPU(*v, PETSC_TRUE));
597: }
598: }
599: }
600: }
601: if (osm->pcmat) {
602: PetscCall(MatDestroy(&pc->pmat));
603: pc->pmat = osm->pcmat;
604: osm->pcmat = NULL;
605: }
606: PetscFunctionReturn(PETSC_SUCCESS);
607: }
609: static PetscErrorCode PCSetUpOnBlocks_GASM(PC pc)
610: {
611: PC_GASM *osm = (PC_GASM *)pc->data;
612: PetscInt i;
614: PetscFunctionBegin;
615: for (i = 0; i < osm->n; i++) PetscCall(KSPSetUp(osm->ksp[i]));
616: PetscFunctionReturn(PETSC_SUCCESS);
617: }
619: static PetscErrorCode PCApply_GASM(PC pc, Vec xin, Vec yout)
620: {
621: PC_GASM *osm = (PC_GASM *)pc->data;
622: const PetscScalar *gxarray;
623: PetscScalar *gyarray;
624: PetscInt i, oni, on;
625: Vec x, y;
626: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
628: PetscFunctionBegin;
629: if (osm->pctoouter) {
630: PetscCall(VecScatterBegin(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE));
631: PetscCall(VecScatterEnd(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE));
632: x = osm->pcx;
633: y = osm->pcy;
634: } else {
635: x = xin;
636: y = yout;
637: }
638: /*
639: support for limiting the restriction or interpolation only to the inner
640: subdomain values (leaving the other values 0).
641: */
642: if (!(osm->type & PC_GASM_RESTRICT)) {
643: /* have to zero the work RHS since scatter may leave some slots empty */
644: PetscCall(VecZeroEntries(osm->gx));
645: PetscCall(VecScatterBegin(osm->girestriction, x, osm->gx, INSERT_VALUES, forward));
646: } else {
647: PetscCall(VecScatterBegin(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward));
648: }
649: PetscCall(VecZeroEntries(osm->gy));
650: if (!(osm->type & PC_GASM_RESTRICT)) {
651: PetscCall(VecScatterEnd(osm->girestriction, x, osm->gx, INSERT_VALUES, forward));
652: } else {
653: PetscCall(VecScatterEnd(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward));
654: }
655: /* do the subdomain solves */
656: PetscCall(VecGetArrayRead(osm->gx, &gxarray));
657: PetscCall(VecGetArray(osm->gy, &gyarray));
658: for (i = 0, on = 0; i < osm->n; ++i, on += oni) {
659: PetscCall(ISGetLocalSize(osm->ois[i], &oni));
660: PetscCall(VecPlaceArray(osm->x[i], PetscSafePointerPlusOffset(gxarray, on)));
661: PetscCall(VecPlaceArray(osm->y[i], PetscSafePointerPlusOffset(gyarray, on)));
662: PetscCall(KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]));
663: PetscCall(KSPCheckSolve(osm->ksp[i], pc, osm->y[i]));
664: PetscCall(VecResetArray(osm->x[i]));
665: PetscCall(VecResetArray(osm->y[i]));
666: }
667: PetscCall(VecRestoreArrayRead(osm->gx, &gxarray));
668: PetscCall(VecRestoreArray(osm->gy, &gyarray));
669: /* do we need to zero y? */
670: PetscCall(VecZeroEntries(y));
671: if (!(osm->type & PC_GASM_INTERPOLATE)) {
672: PetscCall(VecScatterBegin(osm->girestriction, osm->gy, y, ADD_VALUES, reverse));
673: PetscCall(VecScatterEnd(osm->girestriction, osm->gy, y, ADD_VALUES, reverse));
674: } else {
675: PetscCall(VecScatterBegin(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse));
676: PetscCall(VecScatterEnd(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse));
677: }
678: if (osm->pctoouter) {
679: PetscCall(VecScatterBegin(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD));
680: PetscCall(VecScatterEnd(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD));
681: }
682: PetscFunctionReturn(PETSC_SUCCESS);
683: }
685: static PetscErrorCode PCMatApply_GASM_Private(PC pc, Mat Xin, Mat Yout, PetscBool transpose)
686: {
687: PC_GASM *osm = (PC_GASM *)pc->data;
688: Mat X, Y, O = NULL, Z, W;
689: Vec x, y;
690: VecType vtype;
691: PetscInt i, m, M, N;
692: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
694: PetscFunctionBegin;
695: PetscCheck(osm->n == 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not yet implemented");
696: PetscCall(MatGetSize(Xin, NULL, &N));
697: if (osm->pctoouter) {
698: PetscCall(VecGetLocalSize(osm->pcx, &m));
699: PetscCall(VecGetSize(osm->pcx, &M));
700: PetscCall(MatGetVecType(Xin, &vtype));
701: PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)osm->pcx), vtype, m, PETSC_DECIDE, M, N, PETSC_DECIDE, NULL, &O));
702: for (i = 0; i < N; ++i) {
703: PetscCall(MatDenseGetColumnVecRead(Xin, i, &x));
704: PetscCall(MatDenseGetColumnVecWrite(O, i, &y));
705: PetscCall(VecScatterBegin(osm->pctoouter, x, y, INSERT_VALUES, SCATTER_REVERSE));
706: PetscCall(VecScatterEnd(osm->pctoouter, x, y, INSERT_VALUES, SCATTER_REVERSE));
707: PetscCall(MatDenseRestoreColumnVecWrite(O, i, &y));
708: PetscCall(MatDenseRestoreColumnVecRead(Xin, i, &x));
709: }
710: X = Y = O;
711: } else {
712: X = Xin;
713: Y = Yout;
714: }
715: /*
716: support for limiting the restriction or interpolation only to the inner
717: subdomain values (leaving the other values 0).
718: */
719: PetscCall(VecGetLocalSize(osm->x[0], &m));
720: PetscCall(VecGetSize(osm->x[0], &M));
721: PetscCall(VecGetType(osm->x[0], &vtype));
722: PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)osm->ois[0]), vtype, m, PETSC_DECIDE, M, N, PETSC_DECIDE, NULL, &Z));
723: for (i = 0; i < N; ++i) {
724: PetscCall(MatDenseGetColumnVecRead(X, i, &x));
725: PetscCall(MatDenseGetColumnVecWrite(Z, i, &y));
726: if ((!transpose && !(osm->type & PC_GASM_RESTRICT)) || (transpose && !(osm->type & PC_GASM_INTERPOLATE))) {
727: /* have to zero the work RHS since scatter may leave some slots empty */
728: PetscCall(VecZeroEntries(y));
729: PetscCall(VecScatterBegin(osm->girestriction, x, y, INSERT_VALUES, forward));
730: PetscCall(VecScatterEnd(osm->girestriction, x, y, INSERT_VALUES, forward));
731: } else {
732: PetscCall(VecScatterBegin(osm->gorestriction, x, y, INSERT_VALUES, forward));
733: PetscCall(VecScatterEnd(osm->gorestriction, x, y, INSERT_VALUES, forward));
734: }
735: PetscCall(MatDenseRestoreColumnVecWrite(Z, i, &y));
736: PetscCall(MatDenseRestoreColumnVecRead(X, i, &x));
737: }
738: PetscCall(MatCreateDenseFromVecType(PetscObjectComm((PetscObject)osm->ois[0]), vtype, m, PETSC_DECIDE, M, N, PETSC_DECIDE, NULL, &W));
739: /* do the subdomain solve */
740: if (!transpose) PetscCall(KSPMatSolve(osm->ksp[0], Z, W));
741: else PetscCall(KSPMatSolveTranspose(osm->ksp[0], Z, W));
742: PetscCall(KSPCheckMatSolve(osm->ksp[0], pc, W));
743: PetscCall(MatDestroy(&Z));
744: /* do we need to zero y? */
745: PetscCall(MatZeroEntries(Y));
746: for (i = 0; i < N; ++i) {
747: PetscCall(MatDenseGetColumnVecWrite(Y, i, &y));
748: PetscCall(MatDenseGetColumnVecRead(W, i, &x));
749: if ((!transpose && !(osm->type & PC_GASM_INTERPOLATE)) || (transpose && !(osm->type & PC_GASM_RESTRICT))) {
750: PetscCall(VecScatterBegin(osm->girestriction, x, y, ADD_VALUES, reverse));
751: PetscCall(VecScatterEnd(osm->girestriction, x, y, ADD_VALUES, reverse));
752: } else {
753: PetscCall(VecScatterBegin(osm->gorestriction, x, y, ADD_VALUES, reverse));
754: PetscCall(VecScatterEnd(osm->gorestriction, x, y, ADD_VALUES, reverse));
755: }
756: PetscCall(MatDenseRestoreColumnVecRead(W, i, &x));
757: if (osm->pctoouter) {
758: PetscCall(MatDenseGetColumnVecWrite(Yout, i, &x));
759: PetscCall(VecScatterBegin(osm->pctoouter, y, x, INSERT_VALUES, SCATTER_FORWARD));
760: PetscCall(VecScatterEnd(osm->pctoouter, y, x, INSERT_VALUES, SCATTER_FORWARD));
761: PetscCall(MatDenseRestoreColumnVecWrite(Yout, i, &x));
762: }
763: PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &y));
764: }
765: PetscCall(MatDestroy(&W));
766: PetscCall(MatDestroy(&O));
767: PetscFunctionReturn(PETSC_SUCCESS);
768: }
770: static PetscErrorCode PCMatApply_GASM(PC pc, Mat X, Mat Y)
771: {
772: PetscFunctionBegin;
773: PetscCall(PCMatApply_GASM_Private(pc, X, Y, PETSC_FALSE));
774: PetscFunctionReturn(PETSC_SUCCESS);
775: }
777: static PetscErrorCode PCMatApplyTranspose_GASM(PC pc, Mat X, Mat Y)
778: {
779: PetscFunctionBegin;
780: PetscCall(PCMatApply_GASM_Private(pc, X, Y, PETSC_TRUE));
781: PetscFunctionReturn(PETSC_SUCCESS);
782: }
784: static PetscErrorCode PCApplyTranspose_GASM(PC pc, Vec xin, Vec yout)
785: {
786: PC_GASM *osm = (PC_GASM *)pc->data;
787: const PetscScalar *gxarray;
788: PetscScalar *gyarray;
789: PetscInt i, oni, on;
790: Vec x, y;
791: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
793: PetscFunctionBegin;
794: if (osm->pctoouter) {
795: PetscCall(VecScatterBegin(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE));
796: PetscCall(VecScatterEnd(osm->pctoouter, xin, osm->pcx, INSERT_VALUES, SCATTER_REVERSE));
797: x = osm->pcx;
798: y = osm->pcy;
799: } else {
800: x = xin;
801: y = yout;
802: }
803: /*
804: Support for limiting the restriction or interpolation to only local
805: subdomain values (leaving the other values 0).
807: Note: these are reversed from the PCApply_GASM() because we are applying the
808: transpose of the three terms
809: */
810: if (!(osm->type & PC_GASM_INTERPOLATE)) {
811: /* have to zero the work RHS since scatter may leave some slots empty */
812: PetscCall(VecZeroEntries(osm->gx));
813: PetscCall(VecScatterBegin(osm->girestriction, x, osm->gx, INSERT_VALUES, forward));
814: } else {
815: PetscCall(VecScatterBegin(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward));
816: }
817: PetscCall(VecZeroEntries(osm->gy));
818: if (!(osm->type & PC_GASM_INTERPOLATE)) {
819: PetscCall(VecScatterEnd(osm->girestriction, x, osm->gx, INSERT_VALUES, forward));
820: } else {
821: PetscCall(VecScatterEnd(osm->gorestriction, x, osm->gx, INSERT_VALUES, forward));
822: }
823: /* do the local solves */
824: PetscCall(VecGetArrayRead(osm->gx, &gxarray));
825: PetscCall(VecGetArray(osm->gy, &gyarray));
826: for (i = 0, on = 0; i < osm->n; ++i, on += oni) { /* Note that the solves are local, so we can go to osm->n, rather than osm->nmax. */
827: PetscCall(ISGetLocalSize(osm->ois[i], &oni));
828: PetscCall(VecPlaceArray(osm->x[i], PetscSafePointerPlusOffset(gxarray, on)));
829: PetscCall(VecPlaceArray(osm->y[i], PetscSafePointerPlusOffset(gyarray, on)));
830: PetscCall(KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i]));
831: PetscCall(KSPCheckSolve(osm->ksp[i], pc, osm->y[i]));
832: PetscCall(VecResetArray(osm->x[i]));
833: PetscCall(VecResetArray(osm->y[i]));
834: }
835: PetscCall(VecRestoreArrayRead(osm->gx, &gxarray));
836: PetscCall(VecRestoreArray(osm->gy, &gyarray));
837: PetscCall(VecZeroEntries(y));
838: if (!(osm->type & PC_GASM_RESTRICT)) {
839: PetscCall(VecScatterBegin(osm->girestriction, osm->gy, y, ADD_VALUES, reverse));
840: PetscCall(VecScatterEnd(osm->girestriction, osm->gy, y, ADD_VALUES, reverse));
841: } else {
842: PetscCall(VecScatterBegin(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse));
843: PetscCall(VecScatterEnd(osm->gorestriction, osm->gy, y, ADD_VALUES, reverse));
844: }
845: if (osm->pctoouter) {
846: PetscCall(VecScatterBegin(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD));
847: PetscCall(VecScatterEnd(osm->pctoouter, y, yout, INSERT_VALUES, SCATTER_FORWARD));
848: }
849: PetscFunctionReturn(PETSC_SUCCESS);
850: }
852: static PetscErrorCode PCReset_GASM(PC pc)
853: {
854: PC_GASM *osm = (PC_GASM *)pc->data;
856: PetscFunctionBegin;
857: if (osm->ksp) {
858: for (PetscInt i = 0; i < osm->n; i++) PetscCall(KSPReset(osm->ksp[i]));
859: }
860: if (osm->pmat) {
861: if (osm->n > 0) {
862: PetscMPIInt size;
864: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
865: if (size > 1) {
866: /* osm->pmat is created by MatCreateSubMatricesMPI(), cannot use MatDestroySubMatrices() */
867: PetscCall(MatDestroyMatrices(osm->n, &osm->pmat));
868: } else {
869: PetscCall(MatDestroySubMatrices(osm->n, &osm->pmat));
870: }
871: }
872: }
873: if (osm->x) {
874: for (PetscInt i = 0; i < osm->n; i++) {
875: PetscCall(VecDestroy(&osm->x[i]));
876: PetscCall(VecDestroy(&osm->y[i]));
877: }
878: }
879: PetscCall(PetscFree(osm->x));
880: PetscCall(PetscFree(osm->y));
881: PetscCall(VecDestroy(&osm->gx));
882: PetscCall(VecDestroy(&osm->gy));
884: PetscCall(VecScatterDestroy(&osm->gorestriction));
885: PetscCall(VecScatterDestroy(&osm->girestriction));
886: if (!osm->user_subdomains) {
887: PetscCall(PCGASMDestroySubdomains(osm->n, &osm->ois, &osm->iis));
888: osm->N = PETSC_DETERMINE;
889: osm->nmax = PETSC_DETERMINE;
890: }
891: PetscCall(VecScatterDestroy(&osm->pctoouter));
892: PetscCall(ISDestroy(&osm->permutationIS));
893: PetscCall(VecDestroy(&osm->pcx));
894: PetscCall(VecDestroy(&osm->pcy));
895: PetscCall(MatDestroy(&osm->permutationP));
896: PetscCall(MatDestroy(&osm->pcmat));
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
900: static PetscErrorCode PCDestroy_GASM(PC pc)
901: {
902: PC_GASM *osm = (PC_GASM *)pc->data;
904: PetscFunctionBegin;
905: PetscCall(PCReset_GASM(pc));
906: /* PCReset will not destroy subdomains, if user_subdomains is true. */
907: PetscCall(PCGASMDestroySubdomains(osm->n, &osm->ois, &osm->iis));
908: if (osm->ksp) {
909: for (PetscInt i = 0; i < osm->n; i++) PetscCall(KSPDestroy(&osm->ksp[i]));
910: PetscCall(PetscFree(osm->ksp));
911: }
912: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSubdomains_C", NULL));
913: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetOverlap_C", NULL));
914: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetType_C", NULL));
915: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSortIndices_C", NULL));
916: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMGetSubKSP_C", NULL));
917: PetscCall(PetscFree(pc->data));
918: PetscFunctionReturn(PETSC_SUCCESS);
919: }
921: static PetscErrorCode PCSetFromOptions_GASM(PC pc, PetscOptionItems PetscOptionsObject)
922: {
923: PC_GASM *osm = (PC_GASM *)pc->data;
924: PetscInt blocks, ovl;
925: PetscBool flg;
926: PCGASMType gasmtype;
928: PetscFunctionBegin;
929: PetscOptionsHeadBegin(PetscOptionsObject, "Generalized additive Schwarz options");
930: PetscCall(PetscOptionsBool("-pc_gasm_use_dm_subdomains", "If subdomains aren't set, use DMCreateDomainDecomposition() to define subdomains.", "PCGASMSetUseDMSubdomains", osm->dm_subdomains, &osm->dm_subdomains, &flg));
931: PetscCall(PetscOptionsInt("-pc_gasm_total_subdomains", "Total number of subdomains across communicator", "PCGASMSetTotalSubdomains", osm->N, &blocks, &flg));
932: if (flg) PetscCall(PCGASMSetTotalSubdomains(pc, blocks));
933: PetscCall(PetscOptionsInt("-pc_gasm_overlap", "Number of overlapping degrees of freedom", "PCGASMSetOverlap", osm->overlap, &ovl, &flg));
934: if (flg) {
935: PetscCall(PCGASMSetOverlap(pc, ovl));
936: osm->dm_subdomains = PETSC_FALSE;
937: }
938: flg = PETSC_FALSE;
939: PetscCall(PetscOptionsEnum("-pc_gasm_type", "Type of restriction/extension", "PCGASMSetType", PCGASMTypes, (PetscEnum)osm->type, (PetscEnum *)&gasmtype, &flg));
940: if (flg) PetscCall(PCGASMSetType(pc, gasmtype));
941: PetscCall(PetscOptionsBool("-pc_gasm_use_hierachical_partitioning", "use hierarchical partitioning", NULL, osm->hierarchicalpartitioning, &osm->hierarchicalpartitioning, &flg));
942: PetscOptionsHeadEnd();
943: PetscFunctionReturn(PETSC_SUCCESS);
944: }
946: /*@
947: PCGASMSetTotalSubdomains - sets the total number of subdomains to use across the communicator for `PCGASM`
949: Logically Collective
951: Input Parameters:
952: + pc - the preconditioner
953: - N - total number of subdomains
955: Level: beginner
957: .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMSetOverlap()`,
958: `PCGASMCreateSubdomains2D()`
959: @*/
960: PetscErrorCode PCGASMSetTotalSubdomains(PC pc, PetscInt N)
961: {
962: PC_GASM *osm = (PC_GASM *)pc->data;
963: PetscMPIInt size, rank;
965: PetscFunctionBegin;
966: PetscCheck(N >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Total number of subdomains must be 1 or more, got N = %" PetscInt_FMT, N);
967: PetscCheck(!pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCGASMSetTotalSubdomains() should be called before calling PCSetUp().");
969: PetscCall(PCGASMDestroySubdomains(osm->n, &osm->iis, &osm->ois));
970: osm->ois = osm->iis = NULL;
972: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
973: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
974: osm->N = N;
975: osm->n = PETSC_DETERMINE;
976: osm->nmax = PETSC_DETERMINE;
977: osm->dm_subdomains = PETSC_FALSE;
978: PetscFunctionReturn(PETSC_SUCCESS);
979: }
981: static PetscErrorCode PCGASMSetSubdomains_GASM(PC pc, PetscInt n, IS iis[], IS ois[])
982: {
983: PC_GASM *osm = (PC_GASM *)pc->data;
985: PetscFunctionBegin;
986: PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Each MPI rank must have 1 or more subdomains, got n = %" PetscInt_FMT, n);
987: PetscCheck(!pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCGASMSetSubdomains() should be called before calling PCSetUp().");
989: PetscCall(PCGASMDestroySubdomains(osm->n, &osm->iis, &osm->ois));
990: osm->iis = osm->ois = NULL;
991: osm->n = n;
992: osm->N = PETSC_DETERMINE;
993: osm->nmax = PETSC_DETERMINE;
994: if (ois) {
995: PetscCall(PetscMalloc1(n, &osm->ois));
996: for (PetscInt i = 0; i < n; i++) {
997: PetscCall(PetscObjectReference((PetscObject)ois[i]));
998: osm->ois[i] = ois[i];
999: }
1000: /*
1001: Since the user set the outer subdomains, even if nontrivial overlap was requested via PCGASMSetOverlap(),
1002: it will be ignored. To avoid confusion later on (e.g., when viewing the PC), the overlap size is set to -1.
1003: */
1004: osm->overlap = -1;
1005: /* inner subdomains must be provided */
1006: PetscCheck(iis, PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "inner indices have to be provided ");
1007: } /* end if */
1008: if (iis) {
1009: PetscCall(PetscMalloc1(n, &osm->iis));
1010: for (PetscInt i = 0; i < n; i++) {
1011: PetscCall(PetscObjectReference((PetscObject)iis[i]));
1012: osm->iis[i] = iis[i];
1013: }
1014: if (!ois) {
1015: osm->ois = NULL;
1016: /* if user does not provide outer indices, we will create the corresponding outer indices using osm->overlap =1 in PCSetUp_GASM */
1017: }
1018: }
1019: if (PetscDefined(USE_DEBUG)) {
1020: PetscInt j, rstart, rend, *covered, lsize;
1021: const PetscInt *indices;
1023: if (osm->iis) {
1024: /* check if the inner indices cover and only cover the local portion of the matrix */
1025: PetscCall(MatGetOwnershipRange(pc->pmat, &rstart, &rend));
1026: PetscCall(PetscCalloc1(rend - rstart, &covered));
1027: /* check if the current MPI process owns indices from others */
1028: for (PetscInt i = 0; i < n; i++) {
1029: PetscCall(ISGetIndices(osm->iis[i], &indices));
1030: PetscCall(ISGetLocalSize(osm->iis[i], &lsize));
1031: for (j = 0; j < lsize; j++) {
1032: PetscCheck(indices[j] >= rstart && indices[j] < rend, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "inner subdomains can not own an index %" PetscInt_FMT " from other ranks", indices[j]);
1033: PetscCheck(covered[indices[j] - rstart] != 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "inner subdomains can not have an overlapping index %" PetscInt_FMT " ", indices[j]);
1034: covered[indices[j] - rstart] = 1;
1035: }
1036: PetscCall(ISRestoreIndices(osm->iis[i], &indices));
1037: }
1038: /* check if we miss any indices */
1039: for (PetscInt i = rstart; i < rend; i++) PetscCheck(covered[i - rstart], PETSC_COMM_SELF, PETSC_ERR_ARG_NULL, "local entity %" PetscInt_FMT " was not covered by inner subdomains", i);
1040: PetscCall(PetscFree(covered));
1041: }
1042: }
1043: if (iis) osm->user_subdomains = PETSC_TRUE;
1044: PetscFunctionReturn(PETSC_SUCCESS);
1045: }
1047: static PetscErrorCode PCGASMSetOverlap_GASM(PC pc, PetscInt ovl)
1048: {
1049: PC_GASM *osm = (PC_GASM *)pc->data;
1051: PetscFunctionBegin;
1052: PetscCheck(ovl >= 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Negative overlap value requested");
1053: PetscCheck(!pc->setupcalled || ovl == osm->overlap, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCGASMSetOverlap() should be called before PCSetUp().");
1054: if (!pc->setupcalled) osm->overlap = ovl;
1055: PetscFunctionReturn(PETSC_SUCCESS);
1056: }
1058: static PetscErrorCode PCGASMSetType_GASM(PC pc, PCGASMType type)
1059: {
1060: PC_GASM *osm = (PC_GASM *)pc->data;
1062: PetscFunctionBegin;
1063: osm->type = type;
1064: osm->type_set = PETSC_TRUE;
1065: PetscFunctionReturn(PETSC_SUCCESS);
1066: }
1068: static PetscErrorCode PCGASMSetSortIndices_GASM(PC pc, PetscBool doSort)
1069: {
1070: PC_GASM *osm = (PC_GASM *)pc->data;
1072: PetscFunctionBegin;
1073: osm->sort_indices = doSort;
1074: PetscFunctionReturn(PETSC_SUCCESS);
1075: }
1077: /*
1078: FIXME: This routine might need to be modified now that multiple processes per subdomain are allowed.
1079: In particular, it would upset the global subdomain number calculation.
1080: */
1081: static PetscErrorCode PCGASMGetSubKSP_GASM(PC pc, PetscInt *n, PetscInt *first, KSP **ksp)
1082: {
1083: PC_GASM *osm = (PC_GASM *)pc->data;
1085: PetscFunctionBegin;
1086: PetscCheck(osm->n >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Need to call PCSetUp() on PC (or KSPSetUp() on the outer KSP object) before calling here");
1088: if (n) *n = osm->n;
1089: if (first) {
1090: PetscCallMPI(MPI_Scan(&osm->n, first, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)pc)));
1091: *first -= osm->n;
1092: }
1093: if (ksp) {
1094: /* Assume that local solves are now different; not necessarily
1095: true, though! This flag is used only for PCView_GASM() */
1096: *ksp = osm->ksp;
1097: osm->same_subdomain_solvers = PETSC_FALSE;
1098: }
1099: PetscFunctionReturn(PETSC_SUCCESS);
1100: } /* PCGASMGetSubKSP_GASM() */
1102: /*@
1103: PCGASMSetSubdomains - Sets the subdomains for this MPI process
1104: for the additive Schwarz preconditioner with multiple MPI processes per subdomain, `PCGASM`
1106: Collective
1108: Input Parameters:
1109: + pc - the preconditioner object
1110: . n - the number of subdomains for this MPI process
1111: . iis - the index sets that define the inner subdomains (or `NULL` for PETSc to determine subdomains), the `iis` array is
1112: copied so may be freed after this call.
1113: - ois - the index sets that define the outer subdomains (or `NULL` to use the same as `iis`, or to construct by expanding `iis` by
1114: the requested overlap), the `ois` array is copied so may be freed after this call.
1116: Level: advanced
1118: Notes:
1119: The `IS` indices use the parallel, global numbering of the vector entries.
1121: Inner subdomains are those where the correction is applied.
1123: Outer subdomains are those where the residual necessary to obtain the
1124: corrections is obtained (see `PCGASMType` for the use of inner/outer subdomains).
1126: Both inner and outer subdomains can extend over several MPI processes.
1127: This process' portion of a subdomain is known as a local subdomain.
1129: Inner subdomains can not overlap with each other, do not have any entities from remote processes,
1130: and have to cover the entire local subdomain owned by the current process. The index sets on each
1131: process should be ordered such that the ith local subdomain is connected to the ith remote subdomain
1132: on another MPI process.
1134: By default the `PGASM` preconditioner uses 1 (local) subdomain per MPI process.
1136: The `iis` and `ois` arrays may be freed after this call using `PCGASMDestroySubdomains()`
1138: .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`, `PCGASMDestroySubdomains()`,
1139: `PCGASMCreateSubdomains2D()`, `PCGASMGetSubdomains()`
1140: @*/
1141: PetscErrorCode PCGASMSetSubdomains(PC pc, PetscInt n, IS iis[], IS ois[])
1142: {
1143: PC_GASM *osm = (PC_GASM *)pc->data;
1145: PetscFunctionBegin;
1147: PetscTryMethod(pc, "PCGASMSetSubdomains_C", (PC, PetscInt, IS[], IS[]), (pc, n, iis, ois));
1148: osm->dm_subdomains = PETSC_FALSE;
1149: PetscFunctionReturn(PETSC_SUCCESS);
1150: }
1152: /*@
1153: PCGASMSetOverlap - Sets the overlap between a pair of subdomains for the
1154: additive Schwarz preconditioner `PCGASM`. Either all or no MPI processes in the
1155: pc communicator must call this routine.
1157: Logically Collective
1159: Input Parameters:
1160: + pc - the preconditioner context
1161: - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 0)
1163: Options Database Key:
1164: . -pc_gasm_overlap overlap - Sets overlap
1166: Level: intermediate
1168: Notes:
1169: By default the `PCGASM` preconditioner uses 1 subdomain per process. To use
1170: multiple subdomain per perocessor or "straddling" subdomains that intersect
1171: multiple processes use `PCGASMSetSubdomains()` (or option `-pc_gasm_total_subdomains` <n>).
1173: The overlap defaults to 0, so if one desires that no additional
1174: overlap be computed beyond what may have been set with a call to
1175: `PCGASMSetSubdomains()`, then `ovl` must be set to be 0. In particular, if one does
1176: not explicitly set the subdomains in application code, then all overlap would be computed
1177: internally by PETSc, and using an overlap of 0 would result in an `PCGASM`
1178: variant that is equivalent to the block Jacobi preconditioner.
1180: One can define initial index sets with any overlap via
1181: `PCGASMSetSubdomains()`; the routine `PCGASMSetOverlap()` merely allows
1182: PETSc to extend that overlap further, if desired.
1184: .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`,
1185: `PCGASMCreateSubdomains2D()`, `PCGASMGetSubdomains()`
1186: @*/
1187: PetscErrorCode PCGASMSetOverlap(PC pc, PetscInt ovl)
1188: {
1189: PC_GASM *osm = (PC_GASM *)pc->data;
1191: PetscFunctionBegin;
1194: PetscTryMethod(pc, "PCGASMSetOverlap_C", (PC, PetscInt), (pc, ovl));
1195: osm->dm_subdomains = PETSC_FALSE;
1196: PetscFunctionReturn(PETSC_SUCCESS);
1197: }
1199: /*@
1200: PCGASMSetType - Sets the type of restriction and interpolation used
1201: for local problems in the `PCGASM` additive Schwarz method.
1203: Logically Collective
1205: Input Parameters:
1206: + pc - the preconditioner context
1207: - type - variant of `PCGASM`, one of
1208: .vb
1209: `PC_GASM_BASIC` - full interpolation and restriction
1210: `PC_GASM_RESTRICT` - full restriction, local MPI process interpolation
1211: `PC_GASM_INTERPOLATE` - full interpolation, local MPI process restriction
1212: `PC_GASM_NONE` - local MPI process restriction and interpolation
1213: .ve
1215: Options Database Key:
1216: . -pc_gasm_type [basic,restrict,interpolate,none] - Sets `PCGASM` type
1218: Level: intermediate
1220: .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`,
1221: `PCGASMCreateSubdomains2D()`, `PCASM`, `PCASMSetType()`
1222: @*/
1223: PetscErrorCode PCGASMSetType(PC pc, PCGASMType type)
1224: {
1225: PetscFunctionBegin;
1228: PetscTryMethod(pc, "PCGASMSetType_C", (PC, PCGASMType), (pc, type));
1229: PetscFunctionReturn(PETSC_SUCCESS);
1230: }
1232: /*@
1233: PCGASMSetSortIndices - Determines whether subdomain indices are sorted.
1235: Logically Collective
1237: Input Parameters:
1238: + pc - the preconditioner context
1239: - doSort - sort the subdomain indices
1241: Level: intermediate
1243: .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`,
1244: `PCGASMCreateSubdomains2D()`
1245: @*/
1246: PetscErrorCode PCGASMSetSortIndices(PC pc, PetscBool doSort)
1247: {
1248: PetscFunctionBegin;
1251: PetscTryMethod(pc, "PCGASMSetSortIndices_C", (PC, PetscBool), (pc, doSort));
1252: PetscFunctionReturn(PETSC_SUCCESS);
1253: }
1255: /*@
1256: PCGASMGetSubKSP - Gets the local `KSP` contexts for all subdomains on this MPI process.
1258: Collective iff first_local is requested
1260: Input Parameter:
1261: . pc - the preconditioner context
1263: Output Parameters:
1264: + n_local - the number of blocks on this MPI process or `NULL`
1265: . first_local - the global number of the first block on this process or `NULL`, all processes must request or all must pass `NULL`
1266: - ksp - the array of `KSP` contexts
1268: Level: advanced
1270: Note:
1271: After `PCGASMGetSubKSP()` the array of `KSP`es is not to be freed
1273: Currently for some matrix implementations only 1 block per MPI process
1274: is supported.
1276: You must call `KSPSetUp()` before calling `PCGASMGetSubKSP()`.
1278: Fortran Note:
1279: Call `PCGASMRestoreSubKSP()` when the array of `KSP` is no longer needed
1281: .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMSetOverlap()`,
1282: `PCGASMCreateSubdomains2D()`
1283: @*/
1284: PetscErrorCode PCGASMGetSubKSP(PC pc, PetscInt *n_local, PetscInt *first_local, KSP *ksp[])
1285: {
1286: PetscFunctionBegin;
1288: PetscUseMethod(pc, "PCGASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (pc, n_local, first_local, ksp));
1289: PetscFunctionReturn(PETSC_SUCCESS);
1290: }
1292: /*MC
1293: PCGASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1294: its own `KSP` object on a subset of MPI processes
1296: Options Database Keys:
1297: + -pc_gasm_total_subdomains n - Sets total number of local subdomains to be distributed among the MPI processes
1298: . -pc_gasm_view_subdomains - activates the printing of subdomain indices in `PCView()`, -ksp_view or -snes_view
1299: . -pc_gasm_print_subdomains - activates the printing of subdomain indices in `PCSetUp()`
1300: . -pc_gasm_overlap ovl - Sets overlap by which to (automatically) extend local subdomains
1301: - -pc_gasm_type (basic|restrict|interpolate|none) - Sets `PCGASMType`
1303: Level: beginner
1305: Notes:
1306: To set options on the solvers for each block append `-sub_` to all the `KSP`, and `PC`
1307: options database keys. For example, `-sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly`
1309: To set the options on the solvers separate for each block call `PCGASMGetSubKSP()`
1310: and set the options directly on the resulting `KSP` object (you can access its `PC`
1311: with `KSPGetPC()`)
1313: `PCGASM` uses host memory for the merged subdomain vectors so that host-only subdomain solvers work with device global vectors.
1314: Device subdomain solvers therefore copy each right-hand side to the device and each solution back to the host during every application.
1316: See {cite}`dryja1987additive` and {cite}`1sbg` for details on additive Schwarz algorithms
1318: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCASM`, `PCGASMType`, `PCGASMSetType()`,
1319: `PCBJACOBI`, `PCGASMGetSubKSP()`, `PCGASMSetSubdomains()`,
1320: `PCSetModifySubMatrices()`, `PCGASMSetOverlap()`, `PCASMSetType()`
1321: M*/
1323: PETSC_EXTERN PetscErrorCode PCCreate_GASM(PC pc)
1324: {
1325: PC_GASM *osm;
1327: PetscFunctionBegin;
1328: PetscCall(PetscNew(&osm));
1330: osm->N = PETSC_DETERMINE;
1331: osm->n = PETSC_DECIDE;
1332: osm->nmax = PETSC_DETERMINE;
1333: osm->overlap = 0;
1334: osm->ksp = NULL;
1335: osm->gorestriction = NULL;
1336: osm->girestriction = NULL;
1337: osm->pctoouter = NULL;
1338: osm->gx = NULL;
1339: osm->gy = NULL;
1340: osm->x = NULL;
1341: osm->y = NULL;
1342: osm->pcx = NULL;
1343: osm->pcy = NULL;
1344: osm->permutationIS = NULL;
1345: osm->permutationP = NULL;
1346: osm->pcmat = NULL;
1347: osm->ois = NULL;
1348: osm->iis = NULL;
1349: osm->pmat = NULL;
1350: osm->type = PC_GASM_RESTRICT;
1351: osm->same_subdomain_solvers = PETSC_TRUE;
1352: osm->sort_indices = PETSC_TRUE;
1353: osm->dm_subdomains = PETSC_FALSE;
1354: osm->hierarchicalpartitioning = PETSC_FALSE;
1356: pc->data = (void *)osm;
1357: pc->ops->apply = PCApply_GASM;
1358: pc->ops->matapply = PCMatApply_GASM;
1359: pc->ops->matapplytranspose = PCMatApplyTranspose_GASM;
1360: pc->ops->applytranspose = PCApplyTranspose_GASM;
1361: pc->ops->setup = PCSetUp_GASM;
1362: pc->ops->reset = PCReset_GASM;
1363: pc->ops->destroy = PCDestroy_GASM;
1364: pc->ops->setfromoptions = PCSetFromOptions_GASM;
1365: pc->ops->setuponblocks = PCSetUpOnBlocks_GASM;
1366: pc->ops->view = PCView_GASM;
1367: pc->ops->applyrichardson = NULL;
1369: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSubdomains_C", PCGASMSetSubdomains_GASM));
1370: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetOverlap_C", PCGASMSetOverlap_GASM));
1371: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetType_C", PCGASMSetType_GASM));
1372: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMSetSortIndices_C", PCGASMSetSortIndices_GASM));
1373: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGASMGetSubKSP_C", PCGASMGetSubKSP_GASM));
1374: PetscFunctionReturn(PETSC_SUCCESS);
1375: }
1377: PetscErrorCode PCGASMCreateLocalSubdomains(Mat A, PetscInt nloc, IS *iis[])
1378: {
1379: MatPartitioning mpart;
1380: const char *prefix;
1381: PetscInt i, j, rstart, rend, bs;
1382: PetscBool hasop, isbaij = PETSC_FALSE, foundpart = PETSC_FALSE;
1383: Mat Ad = NULL, adj;
1384: IS ispart, isnumb, *is;
1386: PetscFunctionBegin;
1387: PetscCheck(nloc >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "number of local subdomains must > 0, got nloc = %" PetscInt_FMT, nloc);
1389: /* Get prefix, row distribution, and block size */
1390: PetscCall(MatGetOptionsPrefix(A, &prefix));
1391: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
1392: PetscCall(MatGetBlockSize(A, &bs));
1393: PetscCheck(rstart / bs * bs == rstart && rend / bs * bs == rend, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "bad row distribution [%" PetscInt_FMT ",%" PetscInt_FMT ") for matrix block size %" PetscInt_FMT, rstart, rend, bs);
1395: /* Get diagonal block from matrix if possible */
1396: PetscCall(MatHasOperation(A, MATOP_GET_DIAGONAL_BLOCK, &hasop));
1397: if (hasop) PetscCall(MatGetDiagonalBlock(A, &Ad));
1398: if (Ad) {
1399: PetscCall(PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQBAIJ, &isbaij));
1400: if (!isbaij) PetscCall(PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQSBAIJ, &isbaij));
1401: }
1402: if (Ad && nloc > 1) {
1403: PetscBool match, done;
1404: /* Try to setup a good matrix partitioning if available */
1405: PetscCall(MatPartitioningCreate(PETSC_COMM_SELF, &mpart));
1406: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mpart, prefix));
1407: PetscCall(MatPartitioningSetFromOptions(mpart));
1408: PetscCall(PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGCURRENT, &match));
1409: if (!match) PetscCall(PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGSQUARE, &match));
1410: if (!match) { /* assume a "good" partitioner is available */
1411: PetscInt na;
1412: const PetscInt *ia, *ja;
1413: PetscCall(MatGetRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done));
1414: if (done) {
1415: /* Build adjacency matrix by hand. Unfortunately a call to
1416: MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1417: remove the block-aij structure and we cannot expect
1418: MatPartitioning to split vertices as we need */
1419: PetscInt i, j, len, nnz, cnt, *iia = NULL, *jja = NULL;
1420: const PetscInt *row;
1421: nnz = 0;
1422: for (i = 0; i < na; i++) { /* count number of nonzeros */
1423: len = ia[i + 1] - ia[i];
1424: row = ja + ia[i];
1425: for (j = 0; j < len; j++) {
1426: if (row[j] == i) { /* don't count diagonal */
1427: len--;
1428: break;
1429: }
1430: }
1431: nnz += len;
1432: }
1433: PetscCall(PetscMalloc1(na + 1, &iia));
1434: PetscCall(PetscMalloc1(nnz, &jja));
1435: nnz = 0;
1436: iia[0] = 0;
1437: for (i = 0; i < na; i++) { /* fill adjacency */
1438: cnt = 0;
1439: len = ia[i + 1] - ia[i];
1440: row = ja + ia[i];
1441: for (j = 0; j < len; j++) {
1442: if (row[j] != i) jja[nnz + cnt++] = row[j]; /* if not diagonal */
1443: }
1444: nnz += cnt;
1445: iia[i + 1] = nnz;
1446: }
1447: /* Partitioning of the adjacency matrix */
1448: PetscCall(MatCreateMPIAdj(PETSC_COMM_SELF, na, na, iia, jja, NULL, &adj));
1449: PetscCall(MatPartitioningSetAdjacency(mpart, adj));
1450: PetscCall(MatPartitioningSetNParts(mpart, nloc));
1451: PetscCall(MatPartitioningApply(mpart, &ispart));
1452: PetscCall(ISPartitioningToNumbering(ispart, &isnumb));
1453: PetscCall(MatDestroy(&adj));
1454: foundpart = PETSC_TRUE;
1455: }
1456: PetscCall(MatRestoreRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done));
1457: }
1458: PetscCall(MatPartitioningDestroy(&mpart));
1459: }
1460: PetscCall(PetscMalloc1(nloc, &is));
1461: if (!foundpart) {
1462: /* Partitioning by contiguous chunks of rows */
1464: PetscInt mbs = (rend - rstart) / bs;
1465: PetscInt start = rstart;
1466: for (i = 0; i < nloc; i++) {
1467: PetscInt count = (mbs / nloc + ((mbs % nloc) > i)) * bs;
1468: PetscCall(ISCreateStride(PETSC_COMM_SELF, count, start, 1, &is[i]));
1469: start += count;
1470: }
1472: } else {
1473: /* Partitioning by adjacency of diagonal block */
1475: const PetscInt *numbering;
1476: PetscInt *count, nidx, *indices, *newidx, start = 0;
1477: /* Get node count in each partition */
1478: PetscCall(PetscMalloc1(nloc, &count));
1479: PetscCall(ISPartitioningCount(ispart, nloc, count));
1480: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1481: for (i = 0; i < nloc; i++) count[i] *= bs;
1482: }
1483: /* Build indices from node numbering */
1484: PetscCall(ISGetLocalSize(isnumb, &nidx));
1485: PetscCall(PetscMalloc1(nidx, &indices));
1486: for (i = 0; i < nidx; i++) indices[i] = i; /* needs to be initialized */
1487: PetscCall(ISGetIndices(isnumb, &numbering));
1488: PetscCall(PetscSortIntWithPermutation(nidx, numbering, indices));
1489: PetscCall(ISRestoreIndices(isnumb, &numbering));
1490: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1491: PetscCall(PetscMalloc1(nidx * bs, &newidx));
1492: for (i = 0; i < nidx; i++) {
1493: for (j = 0; j < bs; j++) newidx[i * bs + j] = indices[i] * bs + j;
1494: }
1495: PetscCall(PetscFree(indices));
1496: nidx *= bs;
1497: indices = newidx;
1498: }
1499: /* Shift to get global indices */
1500: for (i = 0; i < nidx; i++) indices[i] += rstart;
1502: /* Build the index sets for each block */
1503: for (i = 0; i < nloc; i++) {
1504: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, count[i], &indices[start], PETSC_COPY_VALUES, &is[i]));
1505: PetscCall(ISSort(is[i]));
1506: start += count[i];
1507: }
1509: PetscCall(PetscFree(count));
1510: PetscCall(PetscFree(indices));
1511: PetscCall(ISDestroy(&isnumb));
1512: PetscCall(ISDestroy(&ispart));
1513: }
1514: *iis = is;
1515: PetscFunctionReturn(PETSC_SUCCESS);
1516: }
1518: PETSC_INTERN PetscErrorCode PCGASMCreateStraddlingSubdomains(Mat A, PetscInt N, PetscInt *n, IS *iis[])
1519: {
1520: PetscFunctionBegin;
1521: PetscCall(MatSubdomainsCreateCoalesce(A, N, n, iis));
1522: PetscFunctionReturn(PETSC_SUCCESS);
1523: }
1525: /*@
1526: PCGASMCreateSubdomains - Creates `n` index sets defining `n` nonoverlapping subdomains on this MPI process for the `PCGASM` additive
1527: Schwarz preconditioner for a any problem based on its matrix.
1529: Collective
1531: Input Parameters:
1532: + A - The global matrix operator
1533: - N - the number of global subdomains requested
1535: Output Parameters:
1536: + n - the number of subdomains created on this MPI process
1537: - iis - the array of index sets defining the local inner subdomains (on which the correction is applied)
1539: Level: advanced
1541: Notes:
1542: When `N` >= A's communicator size, each subdomain is local -- contained within a single MPI process.
1543: When `N` < size, the subdomains are 'straddling' (process boundaries) and are no longer local.
1544: The resulting subdomains can be use in `PCGASMSetSubdomains`(pc,n,iss,`NULL`). The overlapping
1545: outer subdomains will be automatically generated from these according to the requested amount of
1546: overlap; this is currently supported only with local subdomains.
1548: Use `PCGASMDestroySubdomains()` to free the array and the list of index sets.
1550: .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMDestroySubdomains()`
1551: @*/
1552: PetscErrorCode PCGASMCreateSubdomains(Mat A, PetscInt N, PetscInt *n, IS *iis[])
1553: {
1554: PetscMPIInt size;
1556: PetscFunctionBegin;
1558: PetscAssertPointer(iis, 4);
1560: PetscCheck(N >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Number of subdomains must be > 0, N = %" PetscInt_FMT, N);
1561: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)A), &size));
1562: if (N >= size) {
1563: *n = N / size + (N % size);
1564: PetscCall(PCGASMCreateLocalSubdomains(A, *n, iis));
1565: } else {
1566: PetscCall(PCGASMCreateStraddlingSubdomains(A, N, n, iis));
1567: }
1568: PetscFunctionReturn(PETSC_SUCCESS);
1569: }
1571: /*@
1572: PCGASMDestroySubdomains - Destroys the index sets created with
1573: `PCGASMCreateSubdomains()` or `PCGASMCreateSubdomains2D()`. Should be
1574: called after setting subdomains with `PCGASMSetSubdomains()`.
1576: Collective
1578: Input Parameters:
1579: + n - the number of index sets
1580: . iis - the array of inner subdomains
1581: - ois - the array of outer subdomains, can be `NULL`
1583: Level: intermediate
1585: Note:
1586: This is a convenience subroutine that walks each list,
1587: destroys each `IS` on the list, and then frees the list. At the end the
1588: list pointers are set to `NULL`.
1590: .seealso: [](ch_ksp), `PCGASM`, `PCGASMCreateSubdomains()`, `PCGASMSetSubdomains()`
1591: @*/
1592: PetscErrorCode PCGASMDestroySubdomains(PetscInt n, IS *iis[], IS *ois[])
1593: {
1594: PetscFunctionBegin;
1595: if (n <= 0) PetscFunctionReturn(PETSC_SUCCESS);
1596: if (ois) {
1597: PetscAssertPointer(ois, 3);
1598: if (*ois) {
1599: PetscAssertPointer(*ois, 3);
1600: for (PetscInt i = 0; i < n; i++) PetscCall(ISDestroy(&(*ois)[i]));
1601: PetscCall(PetscFree(*ois));
1602: }
1603: }
1604: if (iis) {
1605: PetscAssertPointer(iis, 2);
1606: if (*iis) {
1607: PetscAssertPointer(*iis, 2);
1608: for (PetscInt i = 0; i < n; i++) PetscCall(ISDestroy(&(*iis)[i]));
1609: PetscCall(PetscFree(*iis));
1610: }
1611: }
1612: PetscFunctionReturn(PETSC_SUCCESS);
1613: }
1615: #define PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, xleft_loc, ylow_loc, xright_loc, yhigh_loc, n) \
1616: do { \
1617: PetscInt first_row = (first) / (M), last_row = (last) / (M) + 1; \
1618: /* \
1619: Compute ylow_loc and yhigh_loc so that (ylow_loc,xleft) and (yhigh_loc,xright) are the corners \
1620: of the bounding box of the intersection of the subdomain with the local ownership range (local \
1621: subdomain). \
1622: Also compute xleft_loc and xright_loc as the lower and upper bounds on the first and last rows \
1623: of the intersection. \
1624: */ \
1625: /* ylow_loc is the grid row containing the first element of the local sumbdomain */ \
1626: *(ylow_loc) = PetscMax(first_row, ylow); \
1627: /* xleft_loc is the offset of first element of the local subdomain within its grid row (might actually be outside the local subdomain) */ \
1628: *(xleft_loc) = *(ylow_loc) == first_row ? PetscMax((first) % (M), xleft) : xleft; \
1629: /* yhigh_loc is the grid row above the last local subdomain element */ \
1630: *(yhigh_loc) = PetscMin(last_row, yhigh); \
1631: /* xright is the offset of the end of the local subdomain within its grid row (might actually be outside the local subdomain) */ \
1632: *(xright_loc) = *(yhigh_loc) == last_row ? PetscMin(xright, (last) % (M)) : xright; \
1633: /* Now compute the size of the local subdomain n. */ \
1634: *(n) = 0; \
1635: if (*(ylow_loc) < *(yhigh_loc)) { \
1636: PetscInt width = (xright) - (xleft); \
1637: *(n) += width * (*(yhigh_loc) - *(ylow_loc) - 1); \
1638: *(n) += PetscMin(PetscMax(*(xright_loc) - (xleft), 0), width); \
1639: *(n) -= PetscMin(PetscMax(*(xleft_loc) - (xleft), 0), width); \
1640: } \
1641: } while (0)
1643: /*@
1644: PCGASMCreateSubdomains2D - Creates the index sets for the `PCGASM` overlapping Schwarz
1645: preconditioner for a two-dimensional problem on a regular grid.
1647: Collective
1649: Input Parameters:
1650: + pc - the preconditioner context
1651: . M - the global number of grid points in the x direction
1652: . N - the global number of grid points in the y direction
1653: . Mdomains - the global number of subdomains in the x direction
1654: . Ndomains - the global number of subdomains in the y direction
1655: . dof - degrees of freedom per node
1656: - overlap - overlap in mesh lines
1658: Output Parameters:
1659: + nsub - the number of local subdomains created
1660: . iis - array of index sets defining inner (nonoverlapping) subdomains
1661: - ois - array of index sets defining outer (overlapping, if overlap > 0) subdomains
1663: Level: advanced
1665: Note:
1666: Use `PCGASMDestroySubdomains()` to free the index sets and the arrays
1668: .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetSubdomains()`, `PCGASMGetSubKSP()`, `PCGASMSetOverlap()`, `PCASMCreateSubdomains2D()`,
1669: `PCGASMDestroySubdomains()`
1670: @*/
1671: PetscErrorCode PCGASMCreateSubdomains2D(PC pc, PetscInt M, PetscInt N, PetscInt Mdomains, PetscInt Ndomains, PetscInt dof, PetscInt overlap, PetscInt *nsub, IS *iis[], IS *ois[])
1672: {
1673: PetscMPIInt size, rank;
1674: PetscInt maxheight, maxwidth;
1675: PetscInt xstart, xleft, xright, xleft_loc, xright_loc;
1676: PetscInt ystart, ylow, yhigh, ylow_loc, yhigh_loc;
1677: PetscInt x[2][2], y[2][2], n[2];
1678: PetscInt first, last;
1679: PetscInt nidx, *idx;
1680: PetscInt ii, jj, s, q, d;
1681: PetscInt k, kk;
1682: PetscMPIInt color;
1683: MPI_Comm comm, subcomm;
1684: IS **xis = NULL, **is = ois, **is_local = iis;
1686: PetscFunctionBegin;
1687: PetscCall(PetscObjectGetComm((PetscObject)pc, &comm));
1688: PetscCallMPI(MPI_Comm_size(comm, &size));
1689: PetscCallMPI(MPI_Comm_rank(comm, &rank));
1690: PetscCall(MatGetOwnershipRange(pc->pmat, &first, &last));
1691: PetscCheck((first % dof) == 0 && (last % dof) == 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE,
1692: "Matrix row partitioning unsuitable for domain decomposition: local row range (%" PetscInt_FMT ",%" PetscInt_FMT ") "
1693: "does not respect the number of degrees of freedom per grid point %" PetscInt_FMT,
1694: first, last, dof);
1696: /* Determine the number of domains with nonzero intersections with the local ownership range. */
1697: s = 0;
1698: ystart = 0;
1699: for (PetscInt j = 0; j < Ndomains; ++j) {
1700: maxheight = N / Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */
1701: PetscCheck(maxheight >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many %" PetscInt_FMT " subdomains in the vertical direction for mesh height %" PetscInt_FMT, Ndomains, N);
1702: /* Vertical domain limits with an overlap. */
1703: ylow = PetscMax(ystart - overlap, 0);
1704: yhigh = PetscMin(ystart + maxheight + overlap, N);
1705: xstart = 0;
1706: for (PetscInt i = 0; i < Mdomains; ++i) {
1707: maxwidth = M / Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */
1708: PetscCheck(maxwidth >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many %" PetscInt_FMT " subdomains in the horizontal direction for mesh width %" PetscInt_FMT, Mdomains, M);
1709: /* Horizontal domain limits with an overlap. */
1710: xleft = PetscMax(xstart - overlap, 0);
1711: xright = PetscMin(xstart + maxwidth + overlap, M);
1712: /*
1713: Determine whether this subdomain intersects this rank's ownership range of pc->pmat.
1714: */
1715: PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, (&xleft_loc), (&ylow_loc), (&xright_loc), (&yhigh_loc), (&nidx));
1716: if (nidx) ++s;
1717: xstart += maxwidth;
1718: } /* for (PetscInt i = 0; i < Mdomains; ++i) */
1719: ystart += maxheight;
1720: } /* for (PetscInt j = 0; j < Ndomains; ++j) */
1722: /* Now we can allocate the necessary number of ISs. */
1723: *nsub = s;
1724: PetscCall(PetscMalloc1(*nsub, is));
1725: PetscCall(PetscMalloc1(*nsub, is_local));
1726: s = 0;
1727: ystart = 0;
1728: for (PetscInt j = 0; j < Ndomains; ++j) {
1729: maxheight = N / Ndomains + ((N % Ndomains) > j); /* Maximal height of subdomain */
1730: PetscCheck(maxheight >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many %" PetscInt_FMT " subdomains in the vertical direction for mesh height %" PetscInt_FMT, Ndomains, N);
1731: /* Vertical domain limits with an overlap. */
1732: y[0][0] = PetscMax(ystart - overlap, 0);
1733: y[0][1] = PetscMin(ystart + maxheight + overlap, N);
1734: /* Vertical domain limits without an overlap. */
1735: y[1][0] = ystart;
1736: y[1][1] = PetscMin(ystart + maxheight, N);
1737: xstart = 0;
1738: for (PetscInt i = 0; i < Mdomains; ++i) {
1739: maxwidth = M / Mdomains + ((M % Mdomains) > i); /* Maximal width of subdomain */
1740: PetscCheck(maxwidth >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many %" PetscInt_FMT " subdomains in the horizontal direction for mesh width %" PetscInt_FMT, Mdomains, M);
1741: /* Horizontal domain limits with an overlap. */
1742: x[0][0] = PetscMax(xstart - overlap, 0);
1743: x[0][1] = PetscMin(xstart + maxwidth + overlap, M);
1744: /* Horizontal domain limits without an overlap. */
1745: x[1][0] = xstart;
1746: x[1][1] = PetscMin(xstart + maxwidth, M);
1747: /*
1748: Determine whether this domain intersects this rank's ownership range of pc->pmat.
1749: Do this twice: first for the domains with overlaps, and once without.
1750: During the first pass create the subcommunicators, and use them on the second pass as well.
1751: */
1752: for (q = 0; q < 2; ++q) {
1753: PetscBool split = PETSC_FALSE;
1754: /*
1755: domain limits, (xleft, xright) and (ylow, yheigh) are adjusted
1756: according to whether the domain with an overlap or without is considered.
1757: */
1758: xleft = x[q][0];
1759: xright = x[q][1];
1760: ylow = y[q][0];
1761: yhigh = y[q][1];
1762: PCGASMLocalSubdomainBounds2D(M, N, xleft, ylow, xright, yhigh, first, last, (&xleft_loc), (&ylow_loc), (&xright_loc), (&yhigh_loc), (&nidx));
1763: nidx *= dof;
1764: n[q] = nidx;
1765: /*
1766: Based on the counted number of indices in the local domain *with an overlap*,
1767: construct a subcommunicator of all the MPI ranks supporting this domain.
1768: Observe that a domain with an overlap might have nontrivial local support,
1769: while the domain without an overlap might not. Hence, the decision to participate
1770: in the subcommunicator must be based on the domain with an overlap.
1771: */
1772: if (q == 0) {
1773: if (nidx) color = 1;
1774: else color = MPI_UNDEFINED;
1775: PetscCallMPI(MPI_Comm_split(comm, color, rank, &subcomm));
1776: split = PETSC_TRUE;
1777: }
1778: /*
1779: Proceed only if the number of local indices *with an overlap* is nonzero.
1780: */
1781: if (n[0]) {
1782: if (q == 0) xis = is;
1783: if (q == 1) {
1784: /*
1785: The IS for the no-overlap subdomain shares a communicator with the overlapping domain.
1786: Moreover, if the overlap is zero, the two ISs are identical.
1787: */
1788: if (overlap == 0) {
1789: (*is_local)[s] = (*is)[s];
1790: PetscCall(PetscObjectReference((PetscObject)(*is)[s]));
1791: continue;
1792: } else {
1793: xis = is_local;
1794: subcomm = ((PetscObject)(*is)[s])->comm;
1795: }
1796: } /* if (q == 1) */
1797: idx = NULL;
1798: PetscCall(PetscMalloc1(nidx, &idx));
1799: if (nidx) {
1800: k = 0;
1801: for (jj = ylow_loc; jj < yhigh_loc; ++jj) {
1802: PetscInt x0 = (jj == ylow_loc) ? xleft_loc : xleft;
1803: PetscInt x1 = (jj == yhigh_loc - 1) ? xright_loc : xright;
1804: kk = dof * (M * jj + x0);
1805: for (ii = x0; ii < x1; ++ii) {
1806: for (d = 0; d < dof; ++d) idx[k++] = kk++;
1807: }
1808: }
1809: }
1810: PetscCall(ISCreateGeneral(subcomm, nidx, idx, PETSC_OWN_POINTER, (*xis) + s));
1811: if (split) PetscCallMPI(MPI_Comm_free(&subcomm));
1812: } /* if (n[0]) */
1813: } /* for (q = 0; q < 2; ++q) */
1814: if (n[0]) ++s;
1815: xstart += maxwidth;
1816: } /* for (PetscInt i = 0; i < Mdomains; ++i) */
1817: ystart += maxheight;
1818: } /* for (PetscInt j = 0; j < Ndomains; ++j) */
1819: PetscFunctionReturn(PETSC_SUCCESS);
1820: }
1822: /*@
1823: PCGASMGetSubdomains - Gets the subdomains supported on this MPI process
1824: for the `PCGASM` additive Schwarz preconditioner.
1826: Not Collective
1828: Input Parameter:
1829: . pc - the preconditioner context
1831: Output Parameters:
1832: + n - the number of subdomains for this MPI process (default value = 1)
1833: . iis - the index sets that define the inner subdomains (without overlap) supported on this process (can be `NULL`)
1834: - ois - the index sets that define the outer subdomains (with overlap) supported on this process (can be `NULL`)
1836: Level: advanced
1838: Notes:
1839: The user is responsible for destroying the `IS`s and freeing the returned arrays, this can be done with
1840: `PCGASMDestroySubdomains()`
1842: The `IS` numbering is in the parallel, global numbering of the vector.
1844: .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`, `PCGASMCreateSubdomains2D()`,
1845: `PCGASMSetSubdomains()`, `PCGASMGetSubmatrices()`, `PCGASMDestroySubdomains()`
1846: @*/
1847: PetscErrorCode PCGASMGetSubdomains(PC pc, PetscInt *n, IS *iis[], IS *ois[])
1848: {
1849: PC_GASM *osm;
1850: PetscBool match;
1852: PetscFunctionBegin;
1854: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match));
1855: PetscCheck(match, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Incorrect object type: expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name);
1856: osm = (PC_GASM *)pc->data;
1857: if (n) *n = osm->n;
1858: if (iis) PetscCall(PetscMalloc1(osm->n, iis));
1859: if (ois) PetscCall(PetscMalloc1(osm->n, ois));
1860: if (iis || ois) {
1861: for (PetscInt i = 0; i < osm->n; ++i) {
1862: if (iis) (*iis)[i] = osm->iis[i];
1863: if (ois) (*ois)[i] = osm->ois[i];
1864: }
1865: }
1866: PetscFunctionReturn(PETSC_SUCCESS);
1867: }
1869: /*@
1870: PCGASMGetSubmatrices - Gets the local submatrices (for this MPI process
1871: only) for the `PCGASM` additive Schwarz preconditioner.
1873: Not Collective
1875: Input Parameter:
1876: . pc - the preconditioner context
1878: Output Parameters:
1879: + n - the number of matrices for this MPI process (default value = 1)
1880: - mat - the matrices
1882: Level: advanced
1884: Note:
1885: Matrices returned by this routine have the same communicators as the index sets (`IS`)
1886: used to define subdomains in `PCGASMSetSubdomains()`
1888: .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetOverlap()`, `PCGASMGetSubKSP()`,
1889: `PCGASMCreateSubdomains2D()`, `PCGASMSetSubdomains()`, `PCGASMGetSubdomains()`
1890: @*/
1891: PetscErrorCode PCGASMGetSubmatrices(PC pc, PetscInt *n, Mat *mat[])
1892: {
1893: PC_GASM *osm;
1894: PetscBool match;
1896: PetscFunctionBegin;
1898: PetscAssertPointer(n, 2);
1899: if (mat) PetscAssertPointer(mat, 3);
1900: PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call after KSPSetUp() or PCSetUp().");
1901: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match));
1902: PetscCheck(match, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Expected %s, got %s instead", PCGASM, ((PetscObject)pc)->type_name);
1903: osm = (PC_GASM *)pc->data;
1904: if (n) *n = osm->n;
1905: if (mat) *mat = osm->pmat;
1906: PetscFunctionReturn(PETSC_SUCCESS);
1907: }
1909: /*@
1910: PCGASMSetUseDMSubdomains - Indicates whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible for `PCGASM`
1912: Logically Collective
1914: Input Parameters:
1915: + pc - the preconditioner
1916: - flg - boolean indicating whether to use subdomains defined by the `DM`
1918: Options Database Key:
1919: + -pc_gasm_dm_subdomains - configure subdomains
1920: . -pc_gasm_overlap - set overlap
1921: - -pc_gasm_total_subdomains - set number of subdomains
1923: Level: intermediate
1925: Note:
1926: `PCGASMSetSubdomains()`, `PCGASMSetTotalSubdomains()` or `PCGASMSetOverlap()` take precedence over `PCGASMSetUseDMSubdomains()`,
1927: so setting `PCGASMSetSubdomains()` with nontrivial subdomain ISs or any of `PCGASMSetTotalSubdomains()` and `PCGASMSetOverlap()`
1928: automatically turns the latter off.
1930: .seealso: [](ch_ksp), `PCGASM`, `PCGASMGetUseDMSubdomains()`, `PCGASMSetSubdomains()`, `PCGASMSetOverlap()`,
1931: `PCGASMCreateSubdomains2D()`
1932: @*/
1933: PetscErrorCode PCGASMSetUseDMSubdomains(PC pc, PetscBool flg)
1934: {
1935: PC_GASM *osm = (PC_GASM *)pc->data;
1936: PetscBool match;
1938: PetscFunctionBegin;
1941: PetscCheck(!pc->setupcalled, ((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONGSTATE, "Not for a setup PC.");
1942: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match));
1943: if (match && !osm->user_subdomains && osm->N == PETSC_DETERMINE && osm->overlap < 0) osm->dm_subdomains = flg;
1944: PetscFunctionReturn(PETSC_SUCCESS);
1945: }
1947: /*@
1948: PCGASMGetUseDMSubdomains - Returns flag indicating whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible with `PCGASM`
1950: Not Collective
1952: Input Parameter:
1953: . pc - the preconditioner
1955: Output Parameter:
1956: . flg - boolean indicating whether to use subdomains defined by the `DM`
1958: Level: intermediate
1960: .seealso: [](ch_ksp), `PCGASM`, `PCGASMSetUseDMSubdomains()`, `PCGASMSetOverlap()`,
1961: `PCGASMCreateSubdomains2D()`
1962: @*/
1963: PetscErrorCode PCGASMGetUseDMSubdomains(PC pc, PetscBool *flg)
1964: {
1965: PC_GASM *osm = (PC_GASM *)pc->data;
1966: PetscBool match;
1968: PetscFunctionBegin;
1970: PetscAssertPointer(flg, 2);
1971: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCGASM, &match));
1972: if (match) {
1973: if (flg) *flg = osm->dm_subdomains;
1974: }
1975: PetscFunctionReturn(PETSC_SUCCESS);
1976: }