Actual source code: asm.c
1: /*
2: This file defines an additive Schwarz preconditioner for any Mat implementation.
4: Note that each processor may have any number of subdomains. But in order to
5: deal easily with the VecScatter(), we treat each processor as if it has the
6: same number of subdomains.
8: n - total number of true subdomains on all processors
9: n_local_true - actual number of subdomains on this processor
10: n_local = maximum over all processors of n_local_true
11: */
13: #include <petsc/private/pcasmimpl.h>
14: #include <petsc/private/matimpl.h>
16: static PetscErrorCode PCView_ASM(PC pc, PetscViewer viewer)
17: {
18: PC_ASM *osm = (PC_ASM *)pc->data;
19: PetscMPIInt rank;
20: PetscInt i, bsz;
21: PetscBool isascii, isstring;
22: PetscViewer sviewer;
23: PetscViewerFormat format;
24: const char *prefix;
26: PetscFunctionBegin;
27: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
28: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring));
29: if (isascii) {
30: char overlaps[256] = "user-defined overlap", blocks[256] = "total subdomain blocks not yet set";
31: if (osm->overlap >= 0) PetscCall(PetscSNPrintf(overlaps, sizeof(overlaps), "amount of overlap = %" PetscInt_FMT, osm->overlap));
32: if (osm->n > 0) PetscCall(PetscSNPrintf(blocks, sizeof(blocks), "total subdomain blocks = %" PetscInt_FMT, osm->n));
33: PetscCall(PetscViewerASCIIPrintf(viewer, " %s, %s\n", blocks, overlaps));
34: PetscCall(PetscViewerASCIIPrintf(viewer, " restriction/interpolation type - %s\n", PCASMTypes[osm->type]));
35: if (osm->dm_subdomains) PetscCall(PetscViewerASCIIPrintf(viewer, " Additive Schwarz: using DM to define subdomains\n"));
36: if (osm->loctype != PC_COMPOSITE_ADDITIVE) PetscCall(PetscViewerASCIIPrintf(viewer, " Additive Schwarz: local solve composition type - %s\n", PCCompositeTypes[osm->loctype]));
37: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
38: PetscCall(PetscViewerGetFormat(viewer, &format));
39: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
40: if (osm->ksp) {
41: PetscCall(PetscViewerASCIIPrintf(viewer, " Local solver information for first block is in the following KSP and PC objects on rank 0:\n"));
42: PetscCall(PCGetOptionsPrefix(pc, &prefix));
43: PetscCall(PetscViewerASCIIPrintf(viewer, " Use -%sksp_view ::ascii_info_detail to display information for all blocks\n", prefix ? prefix : ""));
44: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
45: if (rank == 0) {
46: PetscCall(PetscViewerASCIIPushTab(sviewer));
47: PetscCall(KSPView(osm->ksp[0], sviewer));
48: PetscCall(PetscViewerASCIIPopTab(sviewer));
49: }
50: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
51: }
52: } else {
53: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
54: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, " [%d] number of local blocks = %" PetscInt_FMT "\n", rank, osm->n_local_true));
55: PetscCall(PetscViewerFlush(viewer));
56: PetscCall(PetscViewerASCIIPrintf(viewer, " Local solver information for each block is in the following KSP and PC objects:\n"));
57: PetscCall(PetscViewerASCIIPushTab(viewer));
58: PetscCall(PetscViewerASCIIPrintf(viewer, "- - - - - - - - - - - - - - - - - -\n"));
59: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
60: for (i = 0; i < osm->n_local_true; i++) {
61: PetscCall(ISGetLocalSize(osm->is[i], &bsz));
62: PetscCall(PetscViewerASCIIPrintf(sviewer, "[%d] local block number %" PetscInt_FMT ", size = %" PetscInt_FMT "\n", rank, i, bsz));
63: PetscCall(KSPView(osm->ksp[i], sviewer));
64: PetscCall(PetscViewerASCIIPrintf(sviewer, "- - - - - - - - - - - - - - - - - -\n"));
65: }
66: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
67: PetscCall(PetscViewerASCIIPopTab(viewer));
68: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
69: }
70: } else if (isstring) {
71: PetscCall(PetscViewerStringSPrintf(viewer, " blocks=%" PetscInt_FMT ", overlap=%" PetscInt_FMT ", type=%s", osm->n, osm->overlap, PCASMTypes[osm->type]));
72: PetscCall(PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
73: if (osm->ksp) PetscCall(KSPView(osm->ksp[0], sviewer));
74: PetscCall(PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer));
75: }
76: PetscFunctionReturn(PETSC_SUCCESS);
77: }
79: static PetscErrorCode PCASMPrintSubdomains(PC pc)
80: {
81: PC_ASM *osm = (PC_ASM *)pc->data;
82: const char *prefix;
83: char fname[PETSC_MAX_PATH_LEN + 1];
84: PetscViewer viewer, sviewer;
85: char *s;
86: PetscInt i, j, nidx;
87: const PetscInt *idx;
88: PetscMPIInt rank, size;
90: PetscFunctionBegin;
91: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
92: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
93: PetscCall(PCGetOptionsPrefix(pc, &prefix));
94: PetscCall(PetscOptionsGetString(NULL, prefix, "-pc_asm_print_subdomains", fname, sizeof(fname), NULL));
95: if (fname[0] == 0) PetscCall(PetscStrncpy(fname, "stdout", sizeof(fname)));
96: PetscCall(PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc), fname, &viewer));
97: for (i = 0; i < osm->n_local; i++) {
98: if (i < osm->n_local_true) {
99: PetscCall(ISGetLocalSize(osm->is[i], &nidx));
100: PetscCall(ISGetIndices(osm->is[i], &idx));
101: /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
102: #define len 16 * (nidx + 1) + 512
103: PetscCall(PetscMalloc1(len, &s));
104: PetscCall(PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer));
105: #undef len
106: PetscCall(PetscViewerStringSPrintf(sviewer, "[%d:%d] Subdomain %" PetscInt_FMT " with overlap:\n", rank, size, i));
107: for (j = 0; j < nidx; j++) PetscCall(PetscViewerStringSPrintf(sviewer, "%" PetscInt_FMT " ", idx[j]));
108: PetscCall(ISRestoreIndices(osm->is[i], &idx));
109: PetscCall(PetscViewerStringSPrintf(sviewer, "\n"));
110: PetscCall(PetscViewerDestroy(&sviewer));
111: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
112: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%s", s));
113: PetscCall(PetscViewerFlush(viewer));
114: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
115: PetscCall(PetscFree(s));
116: if (osm->is_local) {
117: /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
118: #define len 16 * (nidx + 1) + 512
119: PetscCall(PetscMalloc1(len, &s));
120: PetscCall(PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer));
121: #undef len
122: PetscCall(PetscViewerStringSPrintf(sviewer, "[%d:%d] Subdomain %" PetscInt_FMT " without overlap:\n", rank, size, i));
123: PetscCall(ISGetLocalSize(osm->is_local[i], &nidx));
124: PetscCall(ISGetIndices(osm->is_local[i], &idx));
125: for (j = 0; j < nidx; j++) PetscCall(PetscViewerStringSPrintf(sviewer, "%" PetscInt_FMT " ", idx[j]));
126: PetscCall(ISRestoreIndices(osm->is_local[i], &idx));
127: PetscCall(PetscViewerStringSPrintf(sviewer, "\n"));
128: PetscCall(PetscViewerDestroy(&sviewer));
129: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
130: PetscCall(PetscViewerASCIISynchronizedPrintf(viewer, "%s", s));
131: PetscCall(PetscViewerFlush(viewer));
132: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
133: PetscCall(PetscFree(s));
134: }
135: } else {
136: /* Participate in collective viewer calls. */
137: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
138: PetscCall(PetscViewerFlush(viewer));
139: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
140: /* Assume either all ranks have is_local or none do. */
141: if (osm->is_local) {
142: PetscCall(PetscViewerASCIIPushSynchronized(viewer));
143: PetscCall(PetscViewerFlush(viewer));
144: PetscCall(PetscViewerASCIIPopSynchronized(viewer));
145: }
146: }
147: }
148: PetscCall(PetscViewerFlush(viewer));
149: PetscCall(PetscViewerDestroy(&viewer));
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: static PetscErrorCode PCSetUp_ASM(PC pc)
154: {
155: PC_ASM *osm = (PC_ASM *)pc->data;
156: PetscBool flg;
157: PetscInt i, m, m_local;
158: MatReuse scall = MAT_REUSE_MATRIX;
159: IS isl;
160: KSP ksp;
161: PC subpc;
162: const char *prefix, *pprefix;
163: Vec vec;
164: DM *domain_dm = NULL;
165: MatNullSpace *nullsp = NULL;
167: PetscFunctionBegin;
168: if (!pc->setupcalled) {
169: PetscInt m;
171: /* Note: if subdomains have been set either via PCASMSetTotalSubdomains() or via PCASMSetLocalSubdomains(), osm->n_local_true will not be PETSC_DECIDE */
172: if (osm->n_local_true == PETSC_DECIDE) {
173: /* no subdomains given */
174: /* try pc->dm first, if allowed */
175: if (osm->dm_subdomains && pc->dm) {
176: PetscInt num_domains;
177: char **domain_names;
178: IS *inner_domain_is, *outer_domain_is;
179: PetscCall(DMCreateDomainDecomposition(pc->dm, &num_domains, &domain_names, &inner_domain_is, &outer_domain_is, &domain_dm));
180: osm->overlap = -1; /* We do not want to increase the overlap of the IS.
181: A future improvement of this code might allow one to use
182: DM-defined subdomains and also increase the overlap,
183: but that is not currently supported */
184: if (num_domains) PetscCall(PCASMSetLocalSubdomains(pc, num_domains, outer_domain_is, inner_domain_is));
185: for (PetscInt d = 0; d < num_domains; ++d) {
186: if (domain_names) PetscCall(PetscFree(domain_names[d]));
187: if (inner_domain_is) PetscCall(ISDestroy(&inner_domain_is[d]));
188: if (outer_domain_is) PetscCall(ISDestroy(&outer_domain_is[d]));
189: }
190: PetscCall(PetscFree(domain_names));
191: PetscCall(PetscFree(inner_domain_is));
192: PetscCall(PetscFree(outer_domain_is));
193: }
194: if (osm->n_local_true == PETSC_DECIDE) {
195: /* still no subdomains; use one subdomain per processor */
196: osm->n_local_true = 1;
197: }
198: }
199: { /* determine the global and max number of subdomains */
200: struct {
201: PetscInt max, sum;
202: } outwork;
203: PetscMPIInt size;
205: outwork.max = osm->n_local_true;
206: outwork.sum = osm->n_local_true;
207: PetscCallMPI(MPIU_Allreduce(MPI_IN_PLACE, &outwork, 1, MPIU_2INT, MPIU_MAXSUM_OP, PetscObjectComm((PetscObject)pc)));
208: osm->n_local = outwork.max;
209: osm->n = outwork.sum;
211: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
212: if (outwork.max == 1 && outwork.sum == size) {
213: /* osm->n_local_true = 1 on all processes, set this option may enable use of optimized MatCreateSubMatrices() implementation */
214: PetscCall(MatSetOption(pc->pmat, MAT_SUBMAT_SINGLEIS, PETSC_TRUE));
215: }
216: }
217: if (!osm->is) { /* create the index sets */
218: PetscCall(PCASMCreateSubdomains(pc->pmat, osm->n_local_true, &osm->is));
219: }
220: if (osm->n_local_true > 1 && !osm->is_local) {
221: PetscCall(PetscMalloc1(osm->n_local_true, &osm->is_local));
222: for (i = 0; i < osm->n_local_true; i++) {
223: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
224: PetscCall(ISDuplicate(osm->is[i], &osm->is_local[i]));
225: PetscCall(ISCopy(osm->is[i], osm->is_local[i]));
226: } else {
227: PetscCall(PetscObjectReference((PetscObject)osm->is[i]));
228: osm->is_local[i] = osm->is[i];
229: }
230: }
231: }
232: PetscCall(PCGetOptionsPrefix(pc, &prefix));
233: if (osm->overlap > 0) {
234: /* Extend the "overlapping" regions by a number of steps */
235: PetscCall(MatIncreaseOverlap(pc->pmat, osm->n_local_true, osm->is, osm->overlap));
236: }
237: if (osm->sort_indices) {
238: for (i = 0; i < osm->n_local_true; i++) {
239: PetscCall(ISSort(osm->is[i]));
240: if (osm->is_local) PetscCall(ISSort(osm->is_local[i]));
241: }
242: }
243: flg = PETSC_FALSE;
244: PetscCall(PetscOptionsHasName(NULL, prefix, "-pc_asm_print_subdomains", &flg));
245: if (flg) PetscCall(PCASMPrintSubdomains(pc));
246: if (!osm->ksp) {
247: /* Create the local solvers */
248: PetscCall(PetscMalloc1(osm->n_local_true, &osm->ksp));
249: if (domain_dm) PetscCall(PetscInfo(pc, "Setting up ASM subproblems using the embedded DM\n"));
250: for (i = 0; i < osm->n_local_true; i++) {
251: PetscCall(KSPCreate(PETSC_COMM_SELF, &ksp));
252: PetscCall(KSPSetNestLevel(ksp, pc->kspnestlevel));
253: PetscCall(KSPSetErrorIfNotConverged(ksp, pc->erroriffailure));
254: PetscCall(PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1));
255: PetscCall(KSPSetType(ksp, KSPPREONLY));
256: PetscCall(KSPGetPC(ksp, &subpc));
257: PetscCall(PCGetOptionsPrefix(pc, &prefix));
258: PetscCall(KSPSetOptionsPrefix(ksp, prefix));
259: PetscCall(KSPAppendOptionsPrefix(ksp, "sub_"));
260: if (domain_dm) {
261: PetscCall(KSPSetDM(ksp, domain_dm[i]));
262: PetscCall(KSPSetDMActive(ksp, KSP_DMACTIVE_ALL, PETSC_FALSE));
263: PetscCall(DMDestroy(&domain_dm[i]));
264: }
265: osm->ksp[i] = ksp;
266: }
267: PetscCall(PetscFree(domain_dm));
268: }
270: PetscCall(ISConcatenate(PETSC_COMM_SELF, osm->n_local_true, osm->is, &osm->lis));
271: PetscCall(ISSortRemoveDups(osm->lis));
272: PetscCall(ISGetLocalSize(osm->lis, &m));
274: scall = MAT_INITIAL_MATRIX;
275: } else {
276: /*
277: Destroy the blocks from the previous iteration
278: */
279: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
280: PetscCall(MatGetNullSpaces(osm->n_local_true, osm->pmat, &nullsp));
281: PetscCall(MatDestroyMatrices(osm->n_local_true, &osm->pmat));
282: scall = MAT_INITIAL_MATRIX;
283: }
284: }
286: /* Destroy previous submatrices of a different type than pc->pmat since MAT_REUSE_MATRIX won't work in that case */
287: if (scall == MAT_REUSE_MATRIX && osm->sub_mat_type) {
288: PetscCall(MatGetNullSpaces(osm->n_local_true, osm->pmat, &nullsp));
289: if (osm->n_local_true > 0) PetscCall(MatDestroySubMatrices(osm->n_local_true, &osm->pmat));
290: scall = MAT_INITIAL_MATRIX;
291: }
293: /* A subsolver may have factored its submatrix in place, which leaves it
294: flagged as factored and so unfillable by MatCreateSubMatrices() below. */
295: if (scall == MAT_REUSE_MATRIX) {
296: for (i = 0; i < osm->n_local_true; i++) PetscCall(MatSetUnfactored(osm->pmat[i]));
297: }
299: /*
300: Extract out the submatrices
301: */
302: PetscCall(MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, osm->is, scall, &osm->pmat));
303: if (scall == MAT_INITIAL_MATRIX) {
304: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc->pmat, &pprefix));
305: for (i = 0; i < osm->n_local_true; i++) PetscCall(PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i], pprefix));
306: if (nullsp) PetscCall(MatRestoreNullSpaces(osm->n_local_true, osm->pmat, &nullsp));
307: }
309: /* Convert the types of the submatrices (if needbe) */
310: if (osm->sub_mat_type) {
311: for (i = 0; i < osm->n_local_true; i++) PetscCall(MatConvert(osm->pmat[i], osm->sub_mat_type, MAT_INPLACE_MATRIX, &osm->pmat[i]));
312: }
314: if (!pc->setupcalled) {
315: VecType vtype;
317: /* Create the local work vectors (from the local matrices) and scatter contexts */
318: PetscCall(MatCreateVecs(pc->pmat, &vec, NULL));
320: PetscCheck(!osm->is_local || osm->n_local_true == 1 || (osm->type != PC_ASM_INTERPOLATE && osm->type != PC_ASM_NONE), PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot use interpolate or none PCASMType if is_local was provided to PCASMSetLocalSubdomains() with more than a single subdomain");
321: if (osm->is_local && osm->type != PC_ASM_BASIC && osm->type != PC_ASM_WEIGHTED && osm->loctype == PC_COMPOSITE_ADDITIVE) PetscCall(PetscMalloc1(osm->n_local_true, &osm->lprolongation));
322: PetscCall(PetscMalloc1(osm->n_local_true, &osm->lrestriction));
323: PetscCall(PetscMalloc1(osm->n_local_true, &osm->x));
324: PetscCall(PetscMalloc1(osm->n_local_true, &osm->y));
326: PetscCall(ISGetLocalSize(osm->lis, &m));
327: PetscCall(ISCreateStride(PETSC_COMM_SELF, m, 0, 1, &isl));
328: PetscCall(MatGetVecType(osm->pmat[0], &vtype));
329: PetscCall(VecCreate(PETSC_COMM_SELF, &osm->lx));
330: PetscCall(VecSetSizes(osm->lx, m, m));
331: PetscCall(VecSetType(osm->lx, vtype));
332: PetscCall(VecDuplicate(osm->lx, &osm->ly));
333: PetscCall(VecScatterCreate(vec, osm->lis, osm->lx, isl, &osm->restriction));
334: PetscCall(ISDestroy(&isl));
336: for (i = 0; i < osm->n_local_true; ++i) {
337: ISLocalToGlobalMapping ltog;
338: IS isll;
339: const PetscInt *idx_is;
340: PetscInt *idx_lis, nout;
342: PetscCall(ISGetLocalSize(osm->is[i], &m));
343: PetscCall(MatCreateVecs(osm->pmat[i], &osm->x[i], NULL));
344: PetscCall(VecDuplicate(osm->x[i], &osm->y[i]));
346: /* generate a scatter from ly to y[i] picking all the overlapping is[i] entries */
347: PetscCall(ISLocalToGlobalMappingCreateIS(osm->lis, <og));
348: PetscCall(ISGetLocalSize(osm->is[i], &m));
349: PetscCall(ISGetIndices(osm->is[i], &idx_is));
350: PetscCall(PetscMalloc1(m, &idx_lis));
351: PetscCall(ISGlobalToLocalMappingApply(ltog, IS_GTOLM_DROP, m, idx_is, &nout, idx_lis));
352: PetscCheck(nout == m, PETSC_COMM_SELF, PETSC_ERR_PLIB, "is not a subset of lis");
353: PetscCall(ISRestoreIndices(osm->is[i], &idx_is));
354: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, m, idx_lis, PETSC_OWN_POINTER, &isll));
355: PetscCall(ISLocalToGlobalMappingDestroy(<og));
356: PetscCall(ISCreateStride(PETSC_COMM_SELF, m, 0, 1, &isl));
357: PetscCall(VecScatterCreate(osm->ly, isll, osm->y[i], isl, &osm->lrestriction[i]));
358: PetscCall(ISDestroy(&isll));
359: PetscCall(ISDestroy(&isl));
360: if (osm->lprolongation) { /* generate a scatter from y[i] to ly picking only the non-overlapping is_local[i] entries */
361: ISLocalToGlobalMapping ltog;
362: IS isll, isll_local;
363: const PetscInt *idx_local;
364: PetscInt *idx1, *idx2, nout;
366: PetscCall(ISGetLocalSize(osm->is_local[i], &m_local));
367: PetscCall(ISGetIndices(osm->is_local[i], &idx_local));
369: PetscCall(ISLocalToGlobalMappingCreateIS(osm->is[i], <og));
370: PetscCall(PetscMalloc1(m_local, &idx1));
371: PetscCall(ISGlobalToLocalMappingApply(ltog, IS_GTOLM_DROP, m_local, idx_local, &nout, idx1));
372: PetscCall(ISLocalToGlobalMappingDestroy(<og));
373: PetscCheck(nout == m_local, PETSC_COMM_SELF, PETSC_ERR_PLIB, "is_local not a subset of is");
374: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, m_local, idx1, PETSC_OWN_POINTER, &isll));
376: PetscCall(ISLocalToGlobalMappingCreateIS(osm->lis, <og));
377: PetscCall(PetscMalloc1(m_local, &idx2));
378: PetscCall(ISGlobalToLocalMappingApply(ltog, IS_GTOLM_DROP, m_local, idx_local, &nout, idx2));
379: PetscCall(ISLocalToGlobalMappingDestroy(<og));
380: PetscCheck(nout == m_local, PETSC_COMM_SELF, PETSC_ERR_PLIB, "is_local not a subset of lis");
381: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, m_local, idx2, PETSC_OWN_POINTER, &isll_local));
383: PetscCall(ISRestoreIndices(osm->is_local[i], &idx_local));
384: PetscCall(VecScatterCreate(osm->y[i], isll, osm->ly, isll_local, &osm->lprolongation[i]));
386: PetscCall(ISDestroy(&isll));
387: PetscCall(ISDestroy(&isll_local));
388: }
389: }
390: PetscCall(VecDestroy(&vec));
391: }
393: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
394: IS *cis;
396: PetscCall(PetscMalloc1(osm->n_local_true, &cis));
397: for (PetscInt c = 0; c < osm->n_local_true; ++c) cis[c] = osm->lis;
398: PetscCall(MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, cis, scall, &osm->lmats));
399: PetscCall(PetscFree(cis));
400: }
402: /* Return control to the user so that the submatrices can be modified (e.g., to apply
403: different boundary conditions for the submatrices than for the global problem) */
404: PetscCall(PCModifySubMatrices(pc, osm->n_local_true, osm->is, osm->is, osm->pmat, pc->modifysubmatricesP));
406: /*
407: Loop over subdomains putting them into local ksp
408: */
409: PetscCall(KSPGetOptionsPrefix(osm->ksp[0], &prefix));
410: for (i = 0; i < osm->n_local_true; i++) {
411: PetscCall(KSPSetOperators(osm->ksp[i], osm->pmat[i], osm->pmat[i]));
412: PetscCall(MatSetOptionsPrefix(osm->pmat[i], prefix));
413: if (!pc->setupcalled) PetscCall(KSPSetFromOptions(osm->ksp[i]));
414: }
415: if (osm->type == PC_ASM_WEIGHTED && osm->computescaling) {
416: if (!osm->scaling) PetscCall(PetscCalloc1(osm->n_local_true, &osm->scaling));
417: for (i = 0; i < osm->n_local_true; i++) {
418: if (!osm->scaling[i]) PetscCall(VecDuplicate(osm->x[i], &osm->scaling[i]));
419: PetscCallBack("PCASMWeightedComputeScalingFn", (*osm->computescaling)(pc, i, osm->scaling[i], osm->computescalingctx));
420: }
421: }
422: PetscFunctionReturn(PETSC_SUCCESS);
423: }
425: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
426: {
427: PC_ASM *osm = (PC_ASM *)pc->data;
428: PetscInt i;
429: KSPConvergedReason reason;
431: PetscFunctionBegin;
432: for (i = 0; i < osm->n_local_true; i++) {
433: PetscCall(KSPSetUp(osm->ksp[i]));
434: PetscCall(KSPGetConvergedReason(osm->ksp[i], &reason));
435: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
436: }
437: PetscFunctionReturn(PETSC_SUCCESS);
438: }
440: static PetscErrorCode PCApply_ASM(PC pc, Vec x, Vec y)
441: {
442: PC_ASM *osm = (PC_ASM *)pc->data;
443: PCASMType type = osm->type == PC_ASM_WEIGHTED ? PC_ASM_BASIC : osm->type; /* PC_ASM_WEIGHTED scatters like PC_ASM_BASIC, then applies the PCASMWeightedSetScaling() weights */
444: PetscInt i, n_local_true = osm->n_local_true;
445: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
447: PetscFunctionBegin;
448: PetscCheck(osm->type != PC_ASM_WEIGHTED || osm->scaling, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Call PCASMWeightedSetScaling() after PCSetUp() before applying PC_ASM_WEIGHTED");
449: /*
450: support for limiting the restriction or interpolation to only local
451: subdomain values (leaving the other values 0).
452: */
453: if (!(type & PC_ASM_RESTRICT)) {
454: forward = SCATTER_FORWARD_LOCAL;
455: /* have to zero the work RHS since scatter may leave some slots empty */
456: PetscCall(VecSet(osm->lx, 0.0));
457: }
458: if (!(type & PC_ASM_INTERPOLATE)) reverse = SCATTER_REVERSE_LOCAL;
460: PetscCheck(osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
461: /* zero the global and the local solutions */
462: PetscCall(VecSet(y, 0.0));
463: PetscCall(VecSet(osm->ly, 0.0));
465: /* copy the global RHS to local RHS including the ghost nodes */
466: PetscCall(VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward));
467: PetscCall(VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward));
469: /* restrict local RHS to the overlapping 0-block RHS */
470: PetscCall(VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward));
471: PetscCall(VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward));
473: /* do the local solves */
474: for (i = 0; i < n_local_true; ++i) {
475: /* solve the overlapping i-block */
476: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0));
477: PetscCall(KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]));
478: PetscCall(KSPCheckSolve(osm->ksp[i], pc, osm->y[i]));
479: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0));
480: if (osm->type == PC_ASM_WEIGHTED) PetscCall(VecPointwiseMult(osm->y[i], osm->scaling[i], osm->y[i]));
482: if (osm->lprolongation && !(type & PC_ASM_INTERPOLATE)) { /* interpolate the non-overlapping i-block solution to the local solution (only for restrictive additive) */
483: PetscCall(VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward));
484: PetscCall(VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward));
485: } else { /* interpolate the overlapping i-block solution to the local solution */
486: PetscCall(VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse));
487: PetscCall(VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse));
488: }
490: if (i < n_local_true - 1) {
491: /* restrict local RHS to the overlapping (i+1)-block RHS */
492: PetscCall(VecScatterBegin(osm->lrestriction[i + 1], osm->lx, osm->x[i + 1], INSERT_VALUES, forward));
493: PetscCall(VecScatterEnd(osm->lrestriction[i + 1], osm->lx, osm->x[i + 1], INSERT_VALUES, forward));
495: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
496: /* update the overlapping (i+1)-block RHS using the current local solution */
497: PetscCall(MatMult(osm->lmats[i + 1], osm->ly, osm->y[i + 1]));
498: PetscCall(VecAXPBY(osm->x[i + 1], -1., 1., osm->y[i + 1]));
499: }
500: }
501: }
502: /* add the local solution to the global solution including the ghost nodes */
503: PetscCall(VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse));
504: PetscCall(VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse));
505: PetscFunctionReturn(PETSC_SUCCESS);
506: }
508: static PetscErrorCode PCMatApply_ASM_Private(PC pc, Mat X, Mat Y, PetscBool transpose)
509: {
510: PC_ASM *osm = (PC_ASM *)pc->data;
511: PCASMType type = osm->type == PC_ASM_WEIGHTED ? PC_ASM_BASIC : osm->type; /* PC_ASM_WEIGHTED scatters like PC_ASM_BASIC, then applies the PCASMWeightedSetScaling() weights */
512: Mat Z, W;
513: Vec x;
514: PetscInt i, m, N;
515: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
517: PetscFunctionBegin;
518: PetscCheck(osm->n_local_true <= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not yet implemented");
519: PetscCheck(osm->type != PC_ASM_WEIGHTED || osm->scaling, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Call PCASMWeightedSetScaling() after PCSetUp() before applying PC_ASM_WEIGHTED");
520: /*
521: support for limiting the restriction or interpolation to only local
522: subdomain values (leaving the other values 0).
523: */
524: if ((!transpose && !(type & PC_ASM_RESTRICT)) || (transpose && !(type & PC_ASM_INTERPOLATE))) {
525: forward = SCATTER_FORWARD_LOCAL;
526: /* have to zero the work RHS since scatter may leave some slots empty */
527: PetscCall(VecSet(osm->lx, 0.0));
528: }
529: if ((!transpose && !(type & PC_ASM_INTERPOLATE)) || (transpose && !(type & PC_ASM_RESTRICT))) reverse = SCATTER_REVERSE_LOCAL;
530: PetscCall(VecGetLocalSize(osm->x[0], &m));
531: PetscCall(MatGetSize(X, NULL, &N));
532: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, m, N, NULL, &Z));
534: PetscCheck(osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
535: /* zero the global and the local solutions */
536: PetscCall(MatZeroEntries(Y));
537: PetscCall(VecSet(osm->ly, 0.0));
539: for (i = 0; i < N; ++i) {
540: PetscCall(MatDenseGetColumnVecRead(X, i, &x));
541: /* copy the global RHS to local RHS including the ghost nodes */
542: PetscCall(VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward));
543: PetscCall(VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward));
544: PetscCall(MatDenseRestoreColumnVecRead(X, i, &x));
546: PetscCall(MatDenseGetColumnVecWrite(Z, i, &x));
547: /* restrict local RHS to the overlapping 0-block RHS */
548: PetscCall(VecScatterBegin(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward));
549: PetscCall(VecScatterEnd(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward));
550: PetscCall(MatDenseRestoreColumnVecWrite(Z, i, &x));
551: }
552: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, m, N, NULL, &W));
553: /* solve the overlapping 0-block */
554: if (!transpose) {
555: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[0], Z, W, 0));
556: PetscCall(KSPMatSolve(osm->ksp[0], Z, W));
557: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[0], Z, W, 0));
558: } else {
559: if (osm->type == PC_ASM_WEIGHTED) PetscCall(MatDiagonalScale(Z, osm->scaling[0], NULL));
560: PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, osm->ksp[0], Z, W, 0));
561: PetscCall(KSPMatSolveTranspose(osm->ksp[0], Z, W));
562: PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, osm->ksp[0], Z, W, 0));
563: }
564: PetscCall(KSPCheckMatSolve(osm->ksp[0], pc, W));
565: if (!transpose && osm->type == PC_ASM_WEIGHTED) PetscCall(MatDiagonalScale(W, osm->scaling[0], NULL));
566: PetscCall(MatDestroy(&Z));
568: for (i = 0; i < N; ++i) {
569: PetscCall(VecSet(osm->ly, 0.0));
570: PetscCall(MatDenseGetColumnVecRead(W, i, &x));
571: if (osm->lprolongation && ((!transpose && !(type & PC_ASM_INTERPOLATE)) || (transpose && !(type & PC_ASM_RESTRICT)))) { /* interpolate the non-overlapping 0-block solution to the local solution (only for restrictive additive) */
572: PetscCall(VecScatterBegin(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward));
573: PetscCall(VecScatterEnd(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward));
574: } else { /* interpolate the overlapping 0-block solution to the local solution */
575: PetscCall(VecScatterBegin(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse));
576: PetscCall(VecScatterEnd(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse));
577: }
578: PetscCall(MatDenseRestoreColumnVecRead(W, i, &x));
580: PetscCall(MatDenseGetColumnVecWrite(Y, i, &x));
581: /* add the local solution to the global solution including the ghost nodes */
582: PetscCall(VecScatterBegin(osm->restriction, osm->ly, x, ADD_VALUES, reverse));
583: PetscCall(VecScatterEnd(osm->restriction, osm->ly, x, ADD_VALUES, reverse));
584: PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &x));
585: }
586: PetscCall(MatDestroy(&W));
587: PetscFunctionReturn(PETSC_SUCCESS);
588: }
590: static PetscErrorCode PCMatApply_ASM(PC pc, Mat X, Mat Y)
591: {
592: PetscFunctionBegin;
593: PetscCall(PCMatApply_ASM_Private(pc, X, Y, PETSC_FALSE));
594: PetscFunctionReturn(PETSC_SUCCESS);
595: }
597: static PetscErrorCode PCMatApplyTranspose_ASM(PC pc, Mat X, Mat Y)
598: {
599: PetscFunctionBegin;
600: PetscCall(PCMatApply_ASM_Private(pc, X, Y, PETSC_TRUE));
601: PetscFunctionReturn(PETSC_SUCCESS);
602: }
604: static PetscErrorCode PCApplyTranspose_ASM(PC pc, Vec x, Vec y)
605: {
606: PC_ASM *osm = (PC_ASM *)pc->data;
607: PCASMType type = osm->type == PC_ASM_WEIGHTED ? PC_ASM_BASIC : osm->type; /* PC_ASM_WEIGHTED scatters like PC_ASM_BASIC, then applies the PCASMWeightedSetScaling() weights */
608: PetscInt i, n_local_true = osm->n_local_true;
609: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
611: PetscFunctionBegin;
612: PetscCheck(osm->n_local_true <= 1 || osm->loctype == PC_COMPOSITE_ADDITIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not yet implemented");
613: PetscCheck(osm->type != PC_ASM_WEIGHTED || osm->scaling, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Call PCASMWeightedSetScaling() after PCSetUp() before applying PC_ASM_WEIGHTED");
614: /*
615: Support for limiting the restriction or interpolation to only local
616: subdomain values (leaving the other values 0).
618: Note: these are reversed from the PCApply_ASM() because we are applying the
619: transpose of the three terms
620: */
622: if (!(type & PC_ASM_INTERPOLATE)) {
623: forward = SCATTER_FORWARD_LOCAL;
624: /* have to zero the work RHS since scatter may leave some slots empty */
625: PetscCall(VecSet(osm->lx, 0.0));
626: }
627: if (!(type & PC_ASM_RESTRICT)) reverse = SCATTER_REVERSE_LOCAL;
629: /* zero the global and the local solutions */
630: PetscCall(VecSet(y, 0.0));
631: PetscCall(VecSet(osm->ly, 0.0));
633: /* Copy the global RHS to local RHS including the ghost nodes */
634: PetscCall(VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward));
635: PetscCall(VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward));
637: /* Restrict local RHS to the overlapping 0-block RHS */
638: PetscCall(VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward));
639: PetscCall(VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward));
641: /* do the local solves */
642: for (i = 0; i < n_local_true; ++i) {
643: /* solve the overlapping i-block */
644: if (osm->type == PC_ASM_WEIGHTED) PetscCall(VecPointwiseMult(osm->x[i], osm->scaling[i], osm->x[i]));
645: PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0));
646: PetscCall(KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i]));
647: PetscCall(KSPCheckSolve(osm->ksp[i], pc, osm->y[i]));
648: PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0));
650: if (osm->lprolongation && !(type & PC_ASM_RESTRICT)) { /* interpolate the non-overlapping i-block solution to the local solution */
651: PetscCall(VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward));
652: PetscCall(VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward));
653: } else { /* interpolate the overlapping i-block solution to the local solution */
654: PetscCall(VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse));
655: PetscCall(VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse));
656: }
658: if (i < n_local_true - 1) {
659: /* Restrict local RHS to the overlapping (i+1)-block RHS */
660: PetscCall(VecScatterBegin(osm->lrestriction[i + 1], osm->lx, osm->x[i + 1], INSERT_VALUES, forward));
661: PetscCall(VecScatterEnd(osm->lrestriction[i + 1], osm->lx, osm->x[i + 1], INSERT_VALUES, forward));
662: }
663: }
664: /* Add the local solution to the global solution including the ghost nodes */
665: PetscCall(VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse));
666: PetscCall(VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse));
667: PetscFunctionReturn(PETSC_SUCCESS);
668: }
670: static PetscErrorCode PCReset_ASM(PC pc)
671: {
672: PC_ASM *osm = (PC_ASM *)pc->data;
674: PetscFunctionBegin;
675: if (osm->scaling) {
676: for (PetscInt i = 0; i < osm->n_local_true; i++) PetscCall(VecDestroy(&osm->scaling[i]));
677: PetscCall(PetscFree(osm->scaling));
678: }
679: if (osm->ksp) {
680: for (PetscInt i = 0; i < osm->n_local_true; i++) PetscCall(KSPReset(osm->ksp[i]));
681: }
682: if (osm->pmat) {
683: if (osm->n_local_true > 0) PetscCall(MatDestroySubMatrices(osm->n_local_true, &osm->pmat));
684: }
685: if (osm->lrestriction) {
686: PetscCall(VecScatterDestroy(&osm->restriction));
687: for (PetscInt i = 0; i < osm->n_local_true; i++) {
688: PetscCall(VecScatterDestroy(&osm->lrestriction[i]));
689: if (osm->lprolongation) PetscCall(VecScatterDestroy(&osm->lprolongation[i]));
690: PetscCall(VecDestroy(&osm->x[i]));
691: PetscCall(VecDestroy(&osm->y[i]));
692: }
693: PetscCall(PetscFree(osm->lrestriction));
694: PetscCall(PetscFree(osm->lprolongation));
695: PetscCall(PetscFree(osm->x));
696: PetscCall(PetscFree(osm->y));
697: }
698: PetscCall(PCASMDestroySubdomains(osm->n_local_true, &osm->is, &osm->is_local));
699: PetscCall(ISDestroy(&osm->lis));
700: PetscCall(VecDestroy(&osm->lx));
701: PetscCall(VecDestroy(&osm->ly));
702: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(MatDestroyMatrices(osm->n_local_true, &osm->lmats));
704: PetscCall(PetscFree(osm->sub_mat_type));
706: osm->is = NULL;
707: osm->is_local = NULL;
708: PetscFunctionReturn(PETSC_SUCCESS);
709: }
711: static PetscErrorCode PCDestroy_ASM(PC pc)
712: {
713: PC_ASM *osm = (PC_ASM *)pc->data;
715: PetscFunctionBegin;
716: PetscCall(PCReset_ASM(pc));
717: if (osm->ksp) {
718: for (PetscInt i = 0; i < osm->n_local_true; i++) PetscCall(KSPDestroy(&osm->ksp[i]));
719: PetscCall(PetscFree(osm->ksp));
720: }
721: PetscCall(PetscFree(pc->data));
723: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetLocalSubdomains_C", NULL));
724: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetTotalSubdomains_C", NULL));
725: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetOverlap_C", NULL));
726: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetType_C", NULL));
727: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetType_C", NULL));
728: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMWeightedSetScaling_C", NULL));
729: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMWeightedSetComputeScaling_C", NULL));
730: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetLocalType_C", NULL));
731: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetLocalType_C", NULL));
732: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetSortIndices_C", NULL));
733: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetSubKSP_C", NULL));
734: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetSubMatType_C", NULL));
735: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetSubMatType_C", NULL));
736: PetscFunctionReturn(PETSC_SUCCESS);
737: }
739: static PetscErrorCode PCSetFromOptions_ASM(PC pc, PetscOptionItems PetscOptionsObject)
740: {
741: PC_ASM *osm = (PC_ASM *)pc->data;
742: PetscInt blocks, ovl;
743: PetscBool flg;
744: PCASMType asmtype;
745: PCCompositeType loctype;
746: char sub_mat_type[256];
748: PetscFunctionBegin;
749: PetscOptionsHeadBegin(PetscOptionsObject, "Additive Schwarz options");
750: PetscCall(PetscOptionsBool("-pc_asm_dm_subdomains", "Use DMCreateDomainDecomposition() to define subdomains", "PCASMSetDMSubdomains", osm->dm_subdomains, &osm->dm_subdomains, &flg));
751: PetscCall(PetscOptionsInt("-pc_asm_blocks", "Number of subdomains", "PCASMSetTotalSubdomains", osm->n, &blocks, &flg));
752: if (flg) {
753: PetscCall(PCASMSetTotalSubdomains(pc, blocks, NULL, NULL));
754: osm->dm_subdomains = PETSC_FALSE;
755: }
756: PetscCall(PetscOptionsInt("-pc_asm_local_blocks", "Number of local subdomains", "PCASMSetLocalSubdomains", osm->n_local_true, &blocks, &flg));
757: if (flg) {
758: PetscCall(PCASMSetLocalSubdomains(pc, blocks, NULL, NULL));
759: osm->dm_subdomains = PETSC_FALSE;
760: }
761: PetscCall(PetscOptionsInt("-pc_asm_overlap", "Number of grid points overlap", "PCASMSetOverlap", osm->overlap, &ovl, &flg));
762: if (flg) {
763: PetscCall(PCASMSetOverlap(pc, ovl));
764: osm->dm_subdomains = PETSC_FALSE;
765: }
766: flg = PETSC_FALSE;
767: PetscCall(PetscOptionsEnum("-pc_asm_type", "Type of restriction/extension", "PCASMSetType", PCASMTypes, (PetscEnum)osm->type, (PetscEnum *)&asmtype, &flg));
768: if (flg) PetscCall(PCASMSetType(pc, asmtype));
769: flg = PETSC_FALSE;
770: PetscCall(PetscOptionsEnum("-pc_asm_local_type", "Type of local solver composition", "PCASMSetLocalType", PCCompositeTypes, (PetscEnum)osm->loctype, (PetscEnum *)&loctype, &flg));
771: if (flg) PetscCall(PCASMSetLocalType(pc, loctype));
772: PetscCall(PetscOptionsFList("-pc_asm_sub_mat_type", "Subsolve Matrix Type", "PCASMSetSubMatType", MatList, NULL, sub_mat_type, sizeof(sub_mat_type), &flg));
773: if (flg) PetscCall(PCASMSetSubMatType(pc, sub_mat_type));
774: PetscOptionsHeadEnd();
775: PetscFunctionReturn(PETSC_SUCCESS);
776: }
778: static PetscErrorCode PCASMSetLocalSubdomains_ASM(PC pc, PetscInt n, IS is[], IS is_local[])
779: {
780: PC_ASM *osm = (PC_ASM *)pc->data;
782: PetscFunctionBegin;
783: PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Each process must have 1 or more blocks, n = %" PetscInt_FMT, n);
784: PetscCheck(!pc->setupcalled || (n == osm->n_local_true && !is), PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCASMSetLocalSubdomains() should be called before calling PCSetUp().");
786: if (!pc->setupcalled) {
787: if (is) {
788: for (PetscInt i = 0; i < n; i++) PetscCall(PetscObjectReference((PetscObject)is[i]));
789: }
790: if (is_local) {
791: for (PetscInt i = 0; i < n; i++) PetscCall(PetscObjectReference((PetscObject)is_local[i]));
792: }
793: PetscCall(PCASMDestroySubdomains(osm->n_local_true, &osm->is, &osm->is_local));
795: if (osm->ksp && osm->n_local_true != n) {
796: for (PetscInt i = 0; i < osm->n_local_true; i++) PetscCall(KSPDestroy(&osm->ksp[i]));
797: PetscCall(PetscFree(osm->ksp));
798: }
800: osm->n_local_true = n;
801: osm->is = NULL;
802: osm->is_local = NULL;
803: if (is) {
804: PetscCall(PetscMalloc1(n, &osm->is));
805: for (PetscInt i = 0; i < n; i++) osm->is[i] = is[i];
806: /* Flag indicating that the user has set overlapping subdomains so PCASM should not increase their size. */
807: osm->overlap = -1;
808: }
809: if (is_local) {
810: PetscCall(PetscMalloc1(n, &osm->is_local));
811: for (PetscInt i = 0; i < n; i++) osm->is_local[i] = is_local[i];
812: if (!is) {
813: PetscCall(PetscMalloc1(osm->n_local_true, &osm->is));
814: for (PetscInt i = 0; i < osm->n_local_true; i++) {
815: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
816: PetscCall(ISDuplicate(osm->is_local[i], &osm->is[i]));
817: PetscCall(ISCopy(osm->is_local[i], osm->is[i]));
818: } else {
819: PetscCall(PetscObjectReference((PetscObject)osm->is_local[i]));
820: osm->is[i] = osm->is_local[i];
821: }
822: }
823: }
824: }
825: }
826: PetscFunctionReturn(PETSC_SUCCESS);
827: }
829: static PetscErrorCode PCASMSetTotalSubdomains_ASM(PC pc, PetscInt N, IS *is, IS *is_local)
830: {
831: PC_ASM *osm = (PC_ASM *)pc->data;
832: PetscMPIInt rank, size;
833: PetscInt n;
835: PetscFunctionBegin;
836: PetscCheck(N >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Number of total blocks must be > 0, N = %" PetscInt_FMT, N);
837: PetscCheck(!is && !is_local, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Use PCASMSetLocalSubdomains() to set specific index sets, they cannot be set globally yet.");
839: /*
840: Split the subdomains equally among all processors
841: */
842: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
843: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
844: n = N / size + ((N % size) > rank);
845: PetscCheck(n, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Process %d must have at least one block: total processors %d total blocks %" PetscInt_FMT, rank, size, N);
846: PetscCheck(!pc->setupcalled || n == osm->n_local_true, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "PCASMSetTotalSubdomains() should be called before PCSetUp().");
847: if (!pc->setupcalled) {
848: PetscCall(PCASMDestroySubdomains(osm->n_local_true, &osm->is, &osm->is_local));
850: osm->n_local_true = n;
851: osm->is = NULL;
852: osm->is_local = NULL;
853: }
854: PetscFunctionReturn(PETSC_SUCCESS);
855: }
857: static PetscErrorCode PCASMSetOverlap_ASM(PC pc, PetscInt ovl)
858: {
859: PC_ASM *osm = (PC_ASM *)pc->data;
861: PetscFunctionBegin;
862: PetscCheck(ovl >= 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Negative overlap value requested");
863: PetscCheck(!pc->setupcalled || ovl == osm->overlap, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCASMSetOverlap() should be called before PCSetUp().");
864: if (!pc->setupcalled) osm->overlap = ovl;
865: PetscFunctionReturn(PETSC_SUCCESS);
866: }
868: static PetscErrorCode PCASMSetType_ASM(PC pc, PCASMType type)
869: {
870: PC_ASM *osm = (PC_ASM *)pc->data;
872: PetscFunctionBegin;
873: PetscCheck(type != PC_ASM_WEIGHTED || osm->loctype == PC_COMPOSITE_ADDITIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PC_ASM_WEIGHTED requires additive local composition");
874: osm->type = type;
875: osm->type_set = PETSC_TRUE;
876: PetscFunctionReturn(PETSC_SUCCESS);
877: }
879: static PetscErrorCode PCASMGetType_ASM(PC pc, PCASMType *type)
880: {
881: PC_ASM *osm = (PC_ASM *)pc->data;
883: PetscFunctionBegin;
884: *type = osm->type;
885: PetscFunctionReturn(PETSC_SUCCESS);
886: }
888: static PetscErrorCode PCASMWeightedSetComputeScaling_ASM(PC pc, PCASMWeightedComputeScalingFn *fn, PetscCtx ctx)
889: {
890: PC_ASM *osm = (PC_ASM *)pc->data;
892: PetscFunctionBegin;
893: osm->computescaling = fn;
894: osm->computescalingctx = ctx;
895: PetscFunctionReturn(PETSC_SUCCESS);
896: }
898: static PetscErrorCode PCASMWeightedSetScaling_ASM(PC pc, PetscInt n, Vec scaling[])
899: {
900: PC_ASM *osm = (PC_ASM *)pc->data;
901: Vec *newscaling;
902: VecType type;
903: PetscInt m, nvec;
904: PetscMPIInt size;
905: PetscBool match;
907: PetscFunctionBegin;
908: PetscCheck(osm->x, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Call PCSetUp() before PCASMWeightedSetScaling() so that the subdomain sizes and ordering are final");
909: PetscCheck(n == osm->n_local_true, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Number of scaling vectors %" PetscInt_FMT " must match number of local subdomains %" PetscInt_FMT, n, osm->n_local_true);
910: for (PetscInt i = 0; i < n; i++) {
911: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)scaling[i]), &size));
912: PetscCheck(size == 1, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Scaling vector %" PetscInt_FMT " must have a single-process communicator; create it with MatCreateVecs() from PCASMGetLocalSubmatrices()", i);
913: PetscCall(VecGetSize(scaling[i], &nvec));
914: PetscCall(VecGetSize(osm->x[i], &m));
915: PetscCheck(nvec == m, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "Scaling vector %" PetscInt_FMT " has size %" PetscInt_FMT ", expected %" PetscInt_FMT "; create it with MatCreateVecs() from PCASMGetLocalSubmatrices()", i, nvec, m);
916: PetscCall(VecGetType(osm->x[i], &type));
917: PetscCall(PetscObjectTypeCompare((PetscObject)scaling[i], type, &match));
918: PetscCheck(match, PETSC_COMM_SELF, PETSC_ERR_ARG_INCOMP, "Scaling vector %" PetscInt_FMT " must have local solver vector type %s; create it with MatCreateVecs() from PCASMGetLocalSubmatrices()", i, type);
919: }
920: PetscCall(PetscMalloc1(n, &newscaling));
921: for (PetscInt i = 0; i < n; i++) {
922: PetscCall(PetscObjectReference((PetscObject)scaling[i]));
923: newscaling[i] = scaling[i];
924: }
925: if (osm->scaling) {
926: for (PetscInt i = 0; i < n; i++) PetscCall(VecDestroy(&osm->scaling[i]));
927: PetscCall(PetscFree(osm->scaling));
928: }
929: osm->scaling = newscaling;
930: PetscFunctionReturn(PETSC_SUCCESS);
931: }
933: static PetscErrorCode PCASMSetLocalType_ASM(PC pc, PCCompositeType type)
934: {
935: PC_ASM *osm = (PC_ASM *)pc->data;
937: PetscFunctionBegin;
938: PetscCheck(type == PC_COMPOSITE_ADDITIVE || type == PC_COMPOSITE_MULTIPLICATIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Only supports additive or multiplicative as the local type");
939: PetscCheck(osm->type != PC_ASM_WEIGHTED || type == PC_COMPOSITE_ADDITIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PC_ASM_WEIGHTED requires additive local composition");
940: osm->loctype = type;
941: PetscFunctionReturn(PETSC_SUCCESS);
942: }
944: static PetscErrorCode PCASMGetLocalType_ASM(PC pc, PCCompositeType *type)
945: {
946: PC_ASM *osm = (PC_ASM *)pc->data;
948: PetscFunctionBegin;
949: *type = osm->loctype;
950: PetscFunctionReturn(PETSC_SUCCESS);
951: }
953: static PetscErrorCode PCASMSetSortIndices_ASM(PC pc, PetscBool doSort)
954: {
955: PC_ASM *osm = (PC_ASM *)pc->data;
957: PetscFunctionBegin;
958: osm->sort_indices = doSort;
959: PetscFunctionReturn(PETSC_SUCCESS);
960: }
962: static PetscErrorCode PCASMGetSubKSP_ASM(PC pc, PetscInt *n_local, PetscInt *first_local, KSP **ksp)
963: {
964: PC_ASM *osm = (PC_ASM *)pc->data;
966: PetscFunctionBegin;
967: PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Need to call PCSetUp() on PC (or KSPSetUp() on the outer KSP object) before calling here");
969: if (n_local) *n_local = osm->n_local_true;
970: if (first_local) {
971: PetscCallMPI(MPI_Scan(&osm->n_local_true, first_local, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)pc)));
972: *first_local -= osm->n_local_true;
973: }
974: if (ksp) *ksp = osm->ksp;
975: PetscFunctionReturn(PETSC_SUCCESS);
976: }
978: static PetscErrorCode PCASMGetSubMatType_ASM(PC pc, MatType *sub_mat_type)
979: {
980: PC_ASM *osm = (PC_ASM *)pc->data;
982: PetscFunctionBegin;
984: PetscAssertPointer(sub_mat_type, 2);
985: *sub_mat_type = osm->sub_mat_type;
986: PetscFunctionReturn(PETSC_SUCCESS);
987: }
989: static PetscErrorCode PCASMSetSubMatType_ASM(PC pc, MatType sub_mat_type)
990: {
991: PC_ASM *osm = (PC_ASM *)pc->data;
993: PetscFunctionBegin;
995: PetscCall(PetscFree(osm->sub_mat_type));
996: PetscCall(PetscStrallocpy(sub_mat_type, (char **)&osm->sub_mat_type));
997: PetscFunctionReturn(PETSC_SUCCESS);
998: }
1000: /*@
1001: PCASMSetLocalSubdomains - Sets the local subdomains (for this processor only) for the additive Schwarz preconditioner `PCASM`.
1003: Collective
1005: Input Parameters:
1006: + pc - the preconditioner context
1007: . n - the number of subdomains for this processor (default value = 1)
1008: . is - the index set that defines the subdomains for this processor (or `NULL` for PETSc to determine subdomains)
1009: the values of the `is` array are copied so you can free the array (not the `IS` in the array) after this call
1010: - is_local - the index sets that define the local part of the subdomains for this processor, not used unless `PCASMType` is `PC_ASM_RESTRICT`
1011: (or `NULL` to not provide these). The values of the `is_local` array are copied so you can free the array
1012: (not the `IS` in the array) after this call
1014: Options Database Key:
1015: . -pc_asm_local_blocks blks - Sets number of local blocks
1017: Level: advanced
1019: Notes:
1020: The `IS` numbering is in the parallel, global numbering of the vector for both `is` and `is_local`
1022: By default the `PCASM` preconditioner uses 1 block per processor.
1024: Use `PCASMSetTotalSubdomains()` to set the subdomains for all processors.
1026: If `is_local` is provided and `PCASMType` is `PC_ASM_RESTRICT` then the solution only over the `is_local` region is interpolated
1027: back to form the global solution (this is the standard restricted additive Schwarz method, RASM)
1029: If `is_local` is provided and `PCASMType` is `PC_ASM_INTERPOLATE` or `PC_ASM_NONE` then an error is generated since there is
1030: no code to handle that case.
1032: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`, `PCASMGetSubKSP()`,
1033: `PCASMCreateSubdomains2D()`, `PCASMGetLocalSubdomains()`, `PCASMType`, `PCASMSetType()`, `PCGASM`
1034: @*/
1035: PetscErrorCode PCASMSetLocalSubdomains(PC pc, PetscInt n, IS is[], IS is_local[])
1036: {
1037: PetscFunctionBegin;
1039: PetscTryMethod(pc, "PCASMSetLocalSubdomains_C", (PC, PetscInt, IS[], IS[]), (pc, n, is, is_local));
1040: PetscFunctionReturn(PETSC_SUCCESS);
1041: }
1043: /*@
1044: PCASMSetTotalSubdomains - Sets the subdomains for all processors for the
1045: additive Schwarz preconditioner, `PCASM`.
1047: Collective, all MPI ranks must pass in the same array of `IS`
1049: Input Parameters:
1050: + pc - the preconditioner context
1051: . N - the number of subdomains for all processors
1052: . is - the index sets that define the subdomains for all processors (or `NULL` to ask PETSc to determine the subdomains)
1053: the values of the `is` array are copied so you can free the array (not the `IS` in the array) after this call
1054: - is_local - the index sets that define the local part of the subdomains for this processor (or `NULL` to not provide this information)
1055: The values of the `is_local` array are copied so you can free the array (not the `IS` in the array) after this call
1057: Options Database Key:
1058: . -pc_asm_blocks blks - Sets total blocks
1060: Level: advanced
1062: Notes:
1063: Currently you cannot use this to set the actual subdomains with the argument `is` or `is_local`.
1065: By default the `PCASM` preconditioner uses 1 block per processor.
1067: These index sets cannot be destroyed until after completion of the
1068: linear solves for which the `PCASM` preconditioner is being used.
1070: Use `PCASMSetLocalSubdomains()` to set local subdomains.
1072: The `IS` numbering is in the parallel, global numbering of the vector for both is and is_local
1074: .seealso: [](ch_ksp), `PCASM`, `PCASMSetLocalSubdomains()`, `PCASMSetOverlap()`, `PCASMGetSubKSP()`,
1075: `PCASMCreateSubdomains2D()`, `PCGASM`
1076: @*/
1077: PetscErrorCode PCASMSetTotalSubdomains(PC pc, PetscInt N, IS is[], IS is_local[])
1078: {
1079: PetscFunctionBegin;
1081: PetscTryMethod(pc, "PCASMSetTotalSubdomains_C", (PC, PetscInt, IS[], IS[]), (pc, N, is, is_local));
1082: PetscFunctionReturn(PETSC_SUCCESS);
1083: }
1085: /*@
1086: PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
1087: additive Schwarz preconditioner, `PCASM`.
1089: Logically Collective
1091: Input Parameters:
1092: + pc - the preconditioner context
1093: - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1)
1095: Options Database Key:
1096: . -pc_asm_overlap ovl - Sets overlap
1098: Level: intermediate
1100: Notes:
1101: By default the `PCASM` preconditioner uses 1 block per processor. To use
1102: multiple blocks per perocessor, see `PCASMSetTotalSubdomains()` and
1103: `PCASMSetLocalSubdomains()` (and the option -pc_asm_blocks <blks>).
1105: The overlap defaults to 1, so if one desires that no additional
1106: overlap be computed beyond what may have been set with a call to
1107: `PCASMSetTotalSubdomains()` or `PCASMSetLocalSubdomains()`, then ovl
1108: must be set to be 0. In particular, if one does not explicitly set
1109: the subdomains an application code, then all overlap would be computed
1110: internally by PETSc, and using an overlap of 0 would result in an `PCASM`
1111: variant that is equivalent to the block Jacobi preconditioner.
1113: The default algorithm used by PETSc to increase overlap is fast, but not scalable,
1114: use the option -mat_increase_overlap_scalable when the problem and number of processes is large.
1116: One can define initial index sets with any overlap via
1117: `PCASMSetLocalSubdomains()`; the routine
1118: `PCASMSetOverlap()` merely allows PETSc to extend that overlap further
1119: if desired.
1121: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetLocalSubdomains()`, `PCASMGetSubKSP()`,
1122: `PCASMCreateSubdomains2D()`, `PCASMGetLocalSubdomains()`, `MatIncreaseOverlap()`, `PCGASM`
1123: @*/
1124: PetscErrorCode PCASMSetOverlap(PC pc, PetscInt ovl)
1125: {
1126: PetscFunctionBegin;
1129: PetscTryMethod(pc, "PCASMSetOverlap_C", (PC, PetscInt), (pc, ovl));
1130: PetscFunctionReturn(PETSC_SUCCESS);
1131: }
1133: /*@
1134: PCASMSetType - Sets the type of restriction and interpolation used
1135: for local problems in the additive Schwarz method, `PCASM`.
1137: Logically Collective
1139: Input Parameters:
1140: + pc - the preconditioner context
1141: - type - variant of `PCASM`, one of
1142: .vb
1143: PC_ASM_NONE - local processor restriction and interpolation
1144: PC_ASM_RESTRICT - full restriction, local processor interpolation (default)
1145: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1146: PC_ASM_BASIC - full interpolation and restriction
1147: PC_ASM_WEIGHTED - full restriction and interpolation with user-provided diagonal scaling
1148: .ve
1150: Options Database Key:
1151: . -pc_asm_type (none|restrict|interpolate|basic|weighted) - Sets `PCASMType`
1153: Level: intermediate
1155: Note:
1156: if the is_local arguments are passed to `PCASMSetLocalSubdomains()` then they are used when `PC_ASM_RESTRICT` has been selected
1157: to limit the local processor interpolation. `PC_ASM_WEIGHTED` ignores these inner index sets and
1158: uses the weights supplied with `PCASMWeightedSetScaling()` instead. Weighted ASM requires additive local composition.
1160: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMGetSubKSP()`,
1161: `PCASMCreateSubdomains2D()`, `PCASMType`, `PCASMWeightedSetScaling()`, `PCASMSetLocalType()`, `PCASMGetLocalType()`, `PCGASM`
1162: @*/
1163: PetscErrorCode PCASMSetType(PC pc, PCASMType type)
1164: {
1165: PetscFunctionBegin;
1168: PetscTryMethod(pc, "PCASMSetType_C", (PC, PCASMType), (pc, type));
1169: PetscFunctionReturn(PETSC_SUCCESS);
1170: }
1172: /*@
1173: PCASMWeightedSetComputeScaling - Sets a callback to compute `PC_ASM_WEIGHTED` scaling during `PCSetUp()`.
1175: Logically Collective
1177: Input Parameters:
1178: + pc - the `PCASM` preconditioner
1179: . fn - function to fill each local scaling vector, or `NULL` to disable the callback
1180: - ctx - function context passed to `fn`
1182: Level: intermediate
1184: Notes:
1185: Register before setup and select `PC_ASM_WEIGHTED`.
1187: Whenever `PCSetUp()` rebuilds weighted ASM,
1188: `fn` is called once per local subdomain, after overlap expansion and index sorting, with a vector
1189: of the correct size and type.
1191: No explicit setup or vector allocation is needed by the caller.
1193: The callback overwrites any existing weights, including those supplied directly by `PCASMWeightedSetScaling()`.
1195: Registration does not trigger setup. Passing `NULL` leaves the current weights in place.
1197: The callback and context survive `PCReset()`. The caller owns `ctx` and must keep it valid
1198: until the callback is replaced or disabled, or the preconditioner is destroyed.
1200: Fortran Note:
1201: `fn` is a subroutine with arguments `(pc, local, scaling, ctx, ierr)`. Pass `PETSC_NULL_FUNCTION` to disable the callback.
1203: .seealso: [](ch_ksp), `PCASM`, `PCASMWeightedComputeScalingFn`, `PCASMWeightedSetScaling()`, `PCASMWeightedGetScaling()`, `PCASMSetType()`
1204: @*/
1205: PetscErrorCode PCASMWeightedSetComputeScaling(PC pc, PCASMWeightedComputeScalingFn *fn, PetscCtx ctx)
1206: {
1207: PetscFunctionBegin;
1209: PetscTryMethod(pc, "PCASMWeightedSetComputeScaling_C", (PC, PCASMWeightedComputeScalingFn *, PetscCtx), (pc, fn, ctx));
1210: PetscFunctionReturn(PETSC_SUCCESS);
1211: }
1213: /*@
1214: PCASMWeightedSetScaling - Sets the diagonal weights for the overlapping local corrections in weighted additive Schwarz.
1216: Logically Collective
1218: Input Parameters:
1219: + pc - the `PCASM` preconditioner
1220: . n - the number of local subdomains
1221: - scaling - one local scaling `Vec` per overlapping subdomain
1223: Level: intermediate
1225: Notes:
1226: Call `PCSetUp()` before this routine, then use `PCASMGetLocalSubdomains()` to obtain the final
1227: overlapping index sets. Each vector must have a single-process communicator and the same length
1228: and ordering as its corresponding overlapping index set, after overlap expansion and sorting.
1229: Its vector type must match the local solver vectors. Compatible vectors can be created with
1230: `MatCreateVecs()` from the corresponding matrix returned by `PCASMGetLocalSubmatrices()`.
1231: These are local subdomain vectors, not vectors in the parallel global layout.
1232: Select `PC_ASM_WEIGHTED` with `PCASMSetType()` or `-pc_asm_type weighted` to use the weights.
1233: The weights are ignored by the other ASM types.
1235: With restriction operators $R_i$, local solvers $A_i^{-1}$, and $D_i = \text{diag}(scaling[i])$, the action is
1236: $B = \sum_i R_i^T D_i A_i^{-1} R_i$.
1237: PETSc uses the supplied weights as-is, without checking whether they are real, non-negative,
1238: or satisfy $\sum_i R_i^T D_i R_i = I$.
1239: The PC increments the reference count of the vectors but does not copy them.
1240: `PCReset()` discards the weights along with the subdomains.
1241: Alternatively, use `PCASMWeightedSetComputeScaling()` to fill internally-created vectors during `PCSetUp()`.
1243: Example Usage:
1244: .vb
1245: KSPGetPC(ksp, &pc);
1246: PCASMSetType(pc, PC_ASM_WEIGHTED);
1247: PCSetUp(pc);
1248: PCASMGetLocalSubmatrices(pc, &n, &submat);
1249: PCASMGetLocalSubdomains(pc, NULL, &is, NULL);
1250: for (i = 0; i < n; i++) {
1251: MatCreateVecs(submat[i], &scaling[i], NULL); // a Vec of the right size and type
1252: // fill scaling[i] in the local ordering of is[i]
1253: }
1254: PCASMWeightedSetScaling(pc, n, scaling);
1255: .ve
1257: .seealso: [](ch_ksp), `PCASM`, `PCASMType`, `PCASMSetType()`, `PCASMWeightedGetScaling()`, `PCASMGetLocalSubdomains()`, `PCASMGetLocalSubmatrices()`, `PCASMSetLocalSubdomains()`, `PCASMSetLocalType()`
1258: @*/
1259: PetscErrorCode PCASMWeightedSetScaling(PC pc, PetscInt n, Vec scaling[])
1260: {
1261: PetscFunctionBegin;
1263: PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of scaling vectors must be nonnegative");
1264: if (n) PetscAssertPointer(scaling, 3);
1266: PetscTryMethod(pc, "PCASMWeightedSetScaling_C", (PC, PetscInt, Vec[]), (pc, n, scaling));
1267: PetscFunctionReturn(PETSC_SUCCESS);
1268: }
1270: /*@
1271: PCASMWeightedGetScaling - Gets the diagonal weights supplied with `PCASMWeightedSetScaling()` or computed by the function provided with `PCASMWeightedSetComputeScaling()`.
1273: Not Collective
1275: Input Parameter:
1276: . pc - the `PCASM` preconditioner
1278: Output Parameters:
1279: + n - if requested, the number of local subdomains for this processor, or zero if no weights have been supplied
1280: - scaling - if requested, the local scaling `Vec`, or `NULL` if none have been supplied
1282: Level: intermediate
1284: Note:
1285: The returned array and its vectors are owned by `pc`; do not free or destroy them. They
1286: are released by `PCReset()` and `PCDestroy()`, and replaced by a further call to
1287: `PCASMWeightedSetScaling()`.
1289: Fortran Note:
1290: Declare `scaling` as `Vec, pointer :: scaling(:)`. It is always returned and is disassociated when no weights
1291: have been supplied; `n` may be `PETSC_NULL_INTEGER`. There is no restore routine, and the pointer must not be
1292: used after `PCReset()`, `PCDestroy()`, or a further call to `PCASMWeightedSetScaling()`.
1294: .seealso: [](ch_ksp), `PCASM`, `PCASMType`, `PCASMWeightedSetScaling()`, `PCASMSetType()`, `PCASMGetLocalSubdomains()`
1295: @*/
1296: PetscErrorCode PCASMWeightedGetScaling(PC pc, PetscInt *n, Vec *scaling[])
1297: {
1298: PC_ASM *osm = (PC_ASM *)pc->data;
1299: PetscBool match;
1301: PetscFunctionBegin;
1303: if (n) PetscAssertPointer(n, 2);
1304: if (scaling) PetscAssertPointer(scaling, 3);
1305: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCASM, &match));
1306: PetscCheck(match, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC is not a PCASM");
1307: if (n) *n = osm->scaling ? osm->n_local_true : 0;
1308: if (scaling) *scaling = osm->scaling;
1309: PetscFunctionReturn(PETSC_SUCCESS);
1310: }
1312: /*@
1313: PCASMGetType - Gets the type of restriction and interpolation used
1314: for local problems in the additive Schwarz method, `PCASM`.
1316: Logically Collective
1318: Input Parameter:
1319: . pc - the preconditioner context
1321: Output Parameter:
1322: . type - variant of `PCASM`, one of
1323: .vb
1324: PC_ASM_NONE - local processor restriction and interpolation
1325: PC_ASM_RESTRICT - full restriction, local processor interpolation
1326: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1327: PC_ASM_BASIC - full interpolation and restriction
1328: PC_ASM_WEIGHTED - full restriction and interpolation with user-provided diagonal scaling
1329: .ve
1331: Options Database Key:
1332: . -pc_asm_type (none|restrict|interpolate|basic|weighted) - Sets `PCASM` type
1334: Level: intermediate
1336: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMGetSubKSP()`, `PCGASM`,
1337: `PCASMCreateSubdomains2D()`, `PCASMType`, `PCASMSetType()`, `PCASMSetLocalType()`, `PCASMGetLocalType()`
1338: @*/
1339: PetscErrorCode PCASMGetType(PC pc, PCASMType *type)
1340: {
1341: PetscFunctionBegin;
1343: PetscUseMethod(pc, "PCASMGetType_C", (PC, PCASMType *), (pc, type));
1344: PetscFunctionReturn(PETSC_SUCCESS);
1345: }
1347: /*@
1348: PCASMSetLocalType - Sets the type of composition used for local problems in the additive Schwarz method, `PCASM`.
1350: Logically Collective
1352: Input Parameters:
1353: + pc - the preconditioner context
1354: - type - type of composition, one of
1355: .vb
1356: PC_COMPOSITE_ADDITIVE - local additive combination
1357: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1358: .ve
1360: Options Database Key:
1361: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1363: Level: intermediate
1365: .seealso: [](ch_ksp), `PCASM`, `PCASMSetType()`, `PCASMGetType()`, `PCASMGetLocalType()`, `PCASMType`, `PCCompositeType`
1366: @*/
1367: PetscErrorCode PCASMSetLocalType(PC pc, PCCompositeType type)
1368: {
1369: PetscFunctionBegin;
1372: PetscTryMethod(pc, "PCASMSetLocalType_C", (PC, PCCompositeType), (pc, type));
1373: PetscFunctionReturn(PETSC_SUCCESS);
1374: }
1376: /*@
1377: PCASMGetLocalType - Gets the type of composition used for local problems in the additive Schwarz method, `PCASM`.
1379: Logically Collective
1381: Input Parameter:
1382: . pc - the preconditioner context
1384: Output Parameter:
1385: . type - type of composition, one of
1386: .vb
1387: PC_COMPOSITE_ADDITIVE - local additive combination
1388: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1389: .ve
1391: Options Database Key:
1392: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1394: Level: intermediate
1396: .seealso: [](ch_ksp), `PCASM`, `PCASMSetType()`, `PCASMGetType()`, `PCASMSetLocalType()`, `PCASMType`, `PCCompositeType`
1397: @*/
1398: PetscErrorCode PCASMGetLocalType(PC pc, PCCompositeType *type)
1399: {
1400: PetscFunctionBegin;
1402: PetscAssertPointer(type, 2);
1403: PetscUseMethod(pc, "PCASMGetLocalType_C", (PC, PCCompositeType *), (pc, type));
1404: PetscFunctionReturn(PETSC_SUCCESS);
1405: }
1407: /*@
1408: PCASMSetSortIndices - Determines whether subdomain indices are sorted.
1410: Logically Collective
1412: Input Parameters:
1413: + pc - the preconditioner context
1414: - doSort - sort the subdomain indices
1416: Level: intermediate
1418: .seealso: [](ch_ksp), `PCASM`, `PCASMSetLocalSubdomains()`, `PCASMSetTotalSubdomains()`, `PCASMGetSubKSP()`,
1419: `PCASMCreateSubdomains2D()`
1420: @*/
1421: PetscErrorCode PCASMSetSortIndices(PC pc, PetscBool doSort)
1422: {
1423: PetscFunctionBegin;
1426: PetscTryMethod(pc, "PCASMSetSortIndices_C", (PC, PetscBool), (pc, doSort));
1427: PetscFunctionReturn(PETSC_SUCCESS);
1428: }
1430: /*@
1431: PCASMGetSubKSP - Gets the local `KSP` contexts for all blocks on
1432: this processor.
1434: Collective iff first_local is requested
1436: Input Parameter:
1437: . pc - the preconditioner context
1439: Output Parameters:
1440: + n_local - the number of blocks on this processor or `NULL`
1441: . first_local - the global number of the first block on this processor or `NULL`, all processors must request or all must pass `NULL`
1442: - ksp - the array of `KSP` contexts
1444: Level: advanced
1446: Notes:
1447: After `PCASMGetSubKSP()` the array of `KSP`s is not to be freed.
1449: You must call `KSPSetUp()` before calling `PCASMGetSubKSP()`.
1451: Fortran Note:
1452: Call `PCASMRestoreSubKSP()` when access to the array of `KSP` is no longer needed. Pass `PETSC_NULL_KSP_POINTER` for `ksp` if not needed.
1454: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`,
1455: `PCASMCreateSubdomains2D()`
1456: @*/
1457: PetscErrorCode PCASMGetSubKSP(PC pc, PetscInt *n_local, PetscInt *first_local, KSP *ksp[])
1458: {
1459: PetscFunctionBegin;
1461: PetscUseMethod(pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (pc, n_local, first_local, ksp));
1462: PetscFunctionReturn(PETSC_SUCCESS);
1463: }
1465: /*MC
1466: PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1467: its own `KSP` object, {cite}`dryja1987additive` and {cite}`1sbg`
1469: Options Database Keys:
1470: + -pc_asm_blocks blks - Sets total blocks. Defaults to one block per MPI process.
1471: . -pc_asm_overlap ovl - Sets overlap
1472: . -pc_asm_type (none|restrict|interpolate|basic|weighted) - Sets `PCASMType`, default is restrict. See `PCASMSetType()`
1473: . -pc_asm_dm_subdomains (true|false) - use subdomains defined by the `DM` with `DMCreateDomainDecomposition()`
1474: - -pc_asm_local_type (additive|multiplicative) - Sets `PCCompositeType`, default is additive. See `PCASMSetLocalType()`
1476: Level: beginner
1478: Notes:
1479: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you
1480: will get a different convergence rate due to the default option of `-pc_asm_type restrict`. Use
1481: `-pc_asm_type basic` to get the same convergence behavior
1483: Each processor can have one or more blocks, but a block cannot be shared by more
1484: than one processor. Use `PCGASM` for subdomains shared by multiple processes.
1486: To set options on the solvers for each block append `-sub_` to all the `KSP`, and `PC`
1487: options database keys. For example, `-sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly`
1489: To set the options on the solvers separate for each block call `PCASMGetSubKSP()`
1490: and set the options directly on the resulting `KSP` object (you can access its `PC` with `KSPGetPC()`)
1492: If the `PC` has an associated `DM`, then, by default, `DMCreateDomainDecomposition()` is used to create the subdomains
1494: Use `PCASMWeightedSetScaling()` with `PC_ASM_WEIGHTED` to supply a diagonal partition of unity on the overlapping subdomains.
1496: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCASMType`, `PCCompositeType`,
1497: `PCBJACOBI`, `PCASMGetSubKSP()`, `PCASMSetLocalSubdomains()`, `PCASMGetType()`, `PCASMSetLocalType()`, `PCASMGetLocalType()`,
1498: `PCASMSetTotalSubdomains()`, `PCSetModifySubMatrices()`, `PCASMSetOverlap()`, `PCASMSetType()`, `PCASMWeightedSetScaling()`
1499: M*/
1501: PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1502: {
1503: PC_ASM *osm;
1505: PetscFunctionBegin;
1506: PetscCall(PetscNew(&osm));
1508: osm->n = PETSC_DECIDE;
1509: osm->n_local = 0;
1510: osm->n_local_true = PETSC_DECIDE;
1511: osm->overlap = 1;
1512: osm->ksp = NULL;
1513: osm->restriction = NULL;
1514: osm->lprolongation = NULL;
1515: osm->lrestriction = NULL;
1516: osm->x = NULL;
1517: osm->y = NULL;
1518: osm->scaling = NULL;
1519: osm->is = NULL;
1520: osm->is_local = NULL;
1521: osm->mat = NULL;
1522: osm->pmat = NULL;
1523: osm->type = PC_ASM_RESTRICT;
1524: osm->loctype = PC_COMPOSITE_ADDITIVE;
1525: osm->sort_indices = PETSC_TRUE;
1526: osm->dm_subdomains = PETSC_FALSE;
1527: osm->sub_mat_type = NULL;
1529: pc->data = (void *)osm;
1530: pc->ops->apply = PCApply_ASM;
1531: pc->ops->matapply = PCMatApply_ASM;
1532: pc->ops->applytranspose = PCApplyTranspose_ASM;
1533: pc->ops->matapplytranspose = PCMatApplyTranspose_ASM;
1534: pc->ops->setup = PCSetUp_ASM;
1535: pc->ops->reset = PCReset_ASM;
1536: pc->ops->destroy = PCDestroy_ASM;
1537: pc->ops->setfromoptions = PCSetFromOptions_ASM;
1538: pc->ops->setuponblocks = PCSetUpOnBlocks_ASM;
1539: pc->ops->view = PCView_ASM;
1540: pc->ops->applyrichardson = NULL;
1542: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetLocalSubdomains_C", PCASMSetLocalSubdomains_ASM));
1543: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetTotalSubdomains_C", PCASMSetTotalSubdomains_ASM));
1544: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetOverlap_C", PCASMSetOverlap_ASM));
1545: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetType_C", PCASMSetType_ASM));
1546: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetType_C", PCASMGetType_ASM));
1547: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMWeightedSetScaling_C", PCASMWeightedSetScaling_ASM));
1548: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMWeightedSetComputeScaling_C", PCASMWeightedSetComputeScaling_ASM));
1549: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetLocalType_C", PCASMSetLocalType_ASM));
1550: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetLocalType_C", PCASMGetLocalType_ASM));
1551: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetSortIndices_C", PCASMSetSortIndices_ASM));
1552: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetSubKSP_C", PCASMGetSubKSP_ASM));
1553: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetSubMatType_C", PCASMGetSubMatType_ASM));
1554: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetSubMatType_C", PCASMSetSubMatType_ASM));
1555: PetscFunctionReturn(PETSC_SUCCESS);
1556: }
1558: /*@
1559: PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1560: preconditioner, `PCASM`, for any problem on a general grid.
1562: Collective
1564: Input Parameters:
1565: + A - The global matrix operator
1566: - n - the number of local blocks
1568: Output Parameter:
1569: . outis - the array of index sets defining the subdomains
1571: Level: advanced
1573: Note:
1574: This generates nonoverlapping subdomains; the `PCASM` will generate the overlap
1575: from these if you use `PCASMSetLocalSubdomains()`
1577: Fortran Note:
1578: `outis` cannot be `PETSC_NULL_IS_POINTER`. Destroy the returned array with `PCASMDestroySubdomains()`,
1579: passing `PETSC_NULL_IS_POINTER` for `is_local` because no local index sets are created.
1581: .seealso: [](ch_ksp), `PCASM`, `PCASMSetLocalSubdomains()`, `PCASMDestroySubdomains()`
1582: @*/
1583: PetscErrorCode PCASMCreateSubdomains(Mat A, PetscInt n, IS *outis[])
1584: {
1585: MatPartitioning mpart;
1586: const char *prefix;
1587: PetscInt i, j, rstart, rend, bs;
1588: PetscBool hasop, isbaij = PETSC_FALSE, foundpart = PETSC_FALSE;
1589: Mat Ad = NULL, adj;
1590: IS ispart, isnumb, *is;
1592: PetscFunctionBegin;
1594: PetscAssertPointer(outis, 3);
1595: PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "number of local blocks must be > 0, n = %" PetscInt_FMT, n);
1597: /* Get prefix, row distribution, and block size */
1598: PetscCall(MatGetOptionsPrefix(A, &prefix));
1599: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
1600: PetscCall(MatGetBlockSize(A, &bs));
1601: 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);
1603: /* Get diagonal block from matrix if possible */
1604: PetscCall(MatHasOperation(A, MATOP_GET_DIAGONAL_BLOCK, &hasop));
1605: if (hasop) PetscCall(MatGetDiagonalBlock(A, &Ad));
1606: if (Ad) {
1607: PetscCall(PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQBAIJ, &isbaij));
1608: if (!isbaij) PetscCall(PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQSBAIJ, &isbaij));
1609: }
1610: if (Ad && n > 1) {
1611: PetscBool match, done;
1612: /* Try to setup a good matrix partitioning if available */
1613: PetscCall(MatPartitioningCreate(PETSC_COMM_SELF, &mpart));
1614: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mpart, prefix));
1615: PetscCall(MatPartitioningSetFromOptions(mpart));
1616: PetscCall(PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGCURRENT, &match));
1617: if (!match) PetscCall(PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGSQUARE, &match));
1618: if (!match) { /* assume a "good" partitioner is available */
1619: PetscInt na;
1620: const PetscInt *ia, *ja;
1621: PetscCall(MatGetRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done));
1622: if (done) {
1623: /* Build adjacency matrix by hand. Unfortunately a call to
1624: MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1625: remove the block-aij structure and we cannot expect
1626: MatPartitioning to split vertices as we need */
1627: PetscInt i, j, len, nnz, cnt, *iia = NULL, *jja = NULL;
1628: const PetscInt *row;
1629: nnz = 0;
1630: for (i = 0; i < na; i++) { /* count number of nonzeros */
1631: len = ia[i + 1] - ia[i];
1632: row = ja + ia[i];
1633: for (j = 0; j < len; j++) {
1634: if (row[j] == i) { /* don't count diagonal */
1635: len--;
1636: break;
1637: }
1638: }
1639: nnz += len;
1640: }
1641: PetscCall(PetscMalloc1(na + 1, &iia));
1642: PetscCall(PetscMalloc1(nnz, &jja));
1643: nnz = 0;
1644: iia[0] = 0;
1645: for (i = 0; i < na; i++) { /* fill adjacency */
1646: cnt = 0;
1647: len = ia[i + 1] - ia[i];
1648: row = ja + ia[i];
1649: for (j = 0; j < len; j++) {
1650: if (row[j] != i) { /* if not diagonal */
1651: jja[nnz + cnt++] = row[j];
1652: }
1653: }
1654: nnz += cnt;
1655: iia[i + 1] = nnz;
1656: }
1657: /* Partitioning of the adjacency matrix */
1658: PetscCall(MatCreateMPIAdj(PETSC_COMM_SELF, na, na, iia, jja, NULL, &adj));
1659: PetscCall(MatPartitioningSetAdjacency(mpart, adj));
1660: PetscCall(MatPartitioningSetNParts(mpart, n));
1661: PetscCall(MatPartitioningApply(mpart, &ispart));
1662: PetscCall(ISPartitioningToNumbering(ispart, &isnumb));
1663: PetscCall(MatDestroy(&adj));
1664: foundpart = PETSC_TRUE;
1665: }
1666: PetscCall(MatRestoreRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done));
1667: }
1668: PetscCall(MatPartitioningDestroy(&mpart));
1669: }
1671: PetscCall(PetscMalloc1(n, &is));
1672: *outis = is;
1674: if (!foundpart) {
1675: /* Partitioning by contiguous chunks of rows */
1677: PetscInt mbs = (rend - rstart) / bs;
1678: PetscInt start = rstart;
1679: for (i = 0; i < n; i++) {
1680: PetscInt count = (mbs / n + ((mbs % n) > i)) * bs;
1681: PetscCall(ISCreateStride(PETSC_COMM_SELF, count, start, 1, &is[i]));
1682: start += count;
1683: }
1685: } else {
1686: /* Partitioning by adjacency of diagonal block */
1688: const PetscInt *numbering;
1689: PetscInt *count, nidx, *indices, *newidx, start = 0;
1690: /* Get node count in each partition */
1691: PetscCall(PetscMalloc1(n, &count));
1692: PetscCall(ISPartitioningCount(ispart, n, count));
1693: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1694: for (i = 0; i < n; i++) count[i] *= bs;
1695: }
1696: /* Build indices from node numbering */
1697: PetscCall(ISGetLocalSize(isnumb, &nidx));
1698: PetscCall(PetscMalloc1(nidx, &indices));
1699: for (i = 0; i < nidx; i++) indices[i] = i; /* needs to be initialized */
1700: PetscCall(ISGetIndices(isnumb, &numbering));
1701: PetscCall(PetscSortIntWithPermutation(nidx, numbering, indices));
1702: PetscCall(ISRestoreIndices(isnumb, &numbering));
1703: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1704: PetscCall(PetscMalloc1(nidx * bs, &newidx));
1705: for (i = 0; i < nidx; i++) {
1706: for (j = 0; j < bs; j++) newidx[i * bs + j] = indices[i] * bs + j;
1707: }
1708: PetscCall(PetscFree(indices));
1709: nidx *= bs;
1710: indices = newidx;
1711: }
1712: /* Shift to get global indices */
1713: for (i = 0; i < nidx; i++) indices[i] += rstart;
1715: /* Build the index sets for each block */
1716: for (i = 0; i < n; i++) {
1717: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, count[i], &indices[start], PETSC_COPY_VALUES, &is[i]));
1718: PetscCall(ISSort(is[i]));
1719: start += count[i];
1720: }
1722: PetscCall(PetscFree(count));
1723: PetscCall(PetscFree(indices));
1724: PetscCall(ISDestroy(&isnumb));
1725: PetscCall(ISDestroy(&ispart));
1726: }
1727: PetscFunctionReturn(PETSC_SUCCESS);
1728: }
1730: /*@
1731: PCASMDestroySubdomains - Destroys the index sets created with
1732: `PCASMCreateSubdomains()` or `PCASMCreateSubdomains2D()`. Should be called after setting subdomains with `PCASMSetLocalSubdomains()`.
1734: Collective
1736: Input Parameters:
1737: + n - the number of index sets
1738: . is - the array of index sets
1739: - is_local - the array of local index sets, can be `NULL`
1741: Level: advanced
1743: Note:
1744: `PCASMCreateSubdomains2D()` also creates an array of non-overlapping local index sets, which must be passed as `is_local`.
1745: `PCASMCreateSubdomains()` creates no local index sets, so pass `NULL` for `is_local`.
1747: Fortran Note:
1748: `is` cannot be `PETSC_NULL_IS_POINTER`. For arrays from `PCASMCreateSubdomains2D()`, pass the returned `is_local`;
1749: for arrays from `PCASMCreateSubdomains()`, pass `PETSC_NULL_IS_POINTER` for `is_local`.
1751: Developer Note:
1752: The `IS` arguments should be a *[]
1754: .seealso: [](ch_ksp), `PCASM`, `PCASMCreateSubdomains()`, `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`
1755: @*/
1756: PetscErrorCode PCASMDestroySubdomains(PetscInt n, IS *is[], IS *is_local[])
1757: {
1758: PetscInt i;
1760: PetscFunctionBegin;
1761: if (n <= 0) PetscFunctionReturn(PETSC_SUCCESS);
1762: if (*is) {
1763: PetscAssertPointer(*is, 2);
1764: for (i = 0; i < n; i++) PetscCall(ISDestroy(&(*is)[i]));
1765: PetscCall(PetscFree(*is));
1766: }
1767: if (is_local && *is_local) {
1768: PetscAssertPointer(*is_local, 3);
1769: for (i = 0; i < n; i++) PetscCall(ISDestroy(&(*is_local)[i]));
1770: PetscCall(PetscFree(*is_local));
1771: }
1772: PetscFunctionReturn(PETSC_SUCCESS);
1773: }
1775: /*@
1776: PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1777: preconditioner, `PCASM`, for a two-dimensional problem on a regular grid.
1779: Not Collective
1781: Input Parameters:
1782: + m - the number of mesh points in the x direction
1783: . n - the number of mesh points in the y direction
1784: . M - the number of subdomains in the x direction
1785: . N - the number of subdomains in the y direction
1786: . dof - degrees of freedom per node
1787: - overlap - overlap in mesh lines
1789: Output Parameters:
1790: + Nsub - the number of subdomains created
1791: . is - array of index sets defining overlapping (if overlap > 0) subdomains
1792: - is_local - array of index sets defining non-overlapping subdomains
1794: Level: advanced
1796: Note:
1797: Presently `PCAMSCreateSubdomains2d()` is valid only for sequential
1798: preconditioners. More general related routines are
1799: `PCASMSetTotalSubdomains()` and `PCASMSetLocalSubdomains()`.
1801: Fortran Note:
1802: Both `is` and `is_local` are created, so neither can be `PETSC_NULL_IS_POINTER`.
1804: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetLocalSubdomains()`, `PCASMGetSubKSP()`,
1805: `PCASMSetOverlap()`
1806: @*/
1807: PetscErrorCode PCASMCreateSubdomains2D(PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt dof, PetscInt overlap, PetscInt *Nsub, IS *is[], IS *is_local[])
1808: {
1809: PetscInt i, j, height, width, ystart, xstart, yleft, yright, xleft, xright, loc_outer;
1810: PetscInt nidx, *idx, loc, ii, jj, count;
1812: PetscFunctionBegin;
1813: PetscCheck(dof == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "dof must be 1");
1815: *Nsub = N * M;
1816: PetscCall(PetscMalloc1(*Nsub, is));
1817: PetscCall(PetscMalloc1(*Nsub, is_local));
1818: ystart = 0;
1819: loc_outer = 0;
1820: for (i = 0; i < N; i++) {
1821: height = n / N + ((n % N) > i); /* height of subdomain */
1822: PetscCheck(height >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many N subdomains for mesh dimension n");
1823: yleft = ystart - overlap;
1824: if (yleft < 0) yleft = 0;
1825: yright = ystart + height + overlap;
1826: if (yright > n) yright = n;
1827: xstart = 0;
1828: for (j = 0; j < M; j++) {
1829: width = m / M + ((m % M) > j); /* width of subdomain */
1830: PetscCheck(width >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many M subdomains for mesh dimension m");
1831: xleft = xstart - overlap;
1832: if (xleft < 0) xleft = 0;
1833: xright = xstart + width + overlap;
1834: if (xright > m) xright = m;
1835: nidx = (xright - xleft) * (yright - yleft);
1836: PetscCall(PetscMalloc1(nidx, &idx));
1837: loc = 0;
1838: for (ii = yleft; ii < yright; ii++) {
1839: count = m * ii + xleft;
1840: for (jj = xleft; jj < xright; jj++) idx[loc++] = count++;
1841: }
1842: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nidx, idx, PETSC_COPY_VALUES, (*is) + loc_outer));
1843: if (overlap == 0) {
1844: PetscCall(PetscObjectReference((PetscObject)(*is)[loc_outer]));
1846: (*is_local)[loc_outer] = (*is)[loc_outer];
1847: } else {
1848: for (loc = 0, ii = ystart; ii < ystart + height; ii++) {
1849: for (jj = xstart; jj < xstart + width; jj++) idx[loc++] = m * ii + jj;
1850: }
1851: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, loc, idx, PETSC_COPY_VALUES, *is_local + loc_outer));
1852: }
1853: PetscCall(PetscFree(idx));
1854: xstart += width;
1855: loc_outer++;
1856: }
1857: ystart += height;
1858: }
1859: for (i = 0; i < *Nsub; i++) PetscCall(ISSort((*is)[i]));
1860: PetscFunctionReturn(PETSC_SUCCESS);
1861: }
1863: /*@
1864: PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
1865: only) for the additive Schwarz preconditioner, `PCASM`.
1867: Not Collective
1869: Input Parameter:
1870: . pc - the preconditioner context
1872: Output Parameters:
1873: + n - if requested, the number of subdomains for this processor (default value = 1)
1874: . is - if requested, the index sets that define the subdomains for this processor
1875: - is_local - if requested, the index sets that define the local part of the subdomains for this processor (can be `NULL`)
1877: Level: advanced
1879: Note:
1880: The `IS` numbering is in the parallel, global numbering of the vector.
1882: Fortran Note:
1883: Pass `PETSC_NULL_IS_POINTER` for `is` or `is_local` if not needed. A requested array that does not exist is returned disassociated.
1885: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`, `PCASMGetSubKSP()`,
1886: `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`, `PCASMGetLocalSubmatrices()`
1887: @*/
1888: PetscErrorCode PCASMGetLocalSubdomains(PC pc, PetscInt *n, IS *is[], IS *is_local[])
1889: {
1890: PC_ASM *osm = (PC_ASM *)pc->data;
1891: PetscBool match;
1893: PetscFunctionBegin;
1895: if (n) PetscAssertPointer(n, 2);
1896: if (is) PetscAssertPointer(is, 3);
1897: if (is_local) PetscAssertPointer(is_local, 4);
1898: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCASM, &match));
1899: PetscCheck(match, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC is not a PCASM");
1900: if (n) *n = osm->n_local_true;
1901: if (is) *is = osm->is;
1902: if (is_local) *is_local = osm->is_local;
1903: PetscFunctionReturn(PETSC_SUCCESS);
1904: }
1906: /*@
1907: PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1908: only) for the additive Schwarz preconditioner, `PCASM`.
1910: Not Collective
1912: Input Parameter:
1913: . pc - the preconditioner context
1915: Output Parameters:
1916: + n - if requested, the number of matrices for this processor (default value = 1)
1917: - mat - if requested, the matrices
1919: Level: advanced
1921: Notes:
1922: Call after `PCSetUp()` (or `KSPSetUp()`) but before `PCApply()` and before `PCSetUpOnBlocks()`)
1924: Usually one would use `PCSetModifySubMatrices()` to change the submatrices in building the preconditioner.
1926: Fortran Note:
1927: Pass `PETSC_NULL_MAT_POINTER` for `mat` if not needed. If the `PC` is not a `PCASM`, `mat` is returned disassociated.
1929: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`, `PCASMGetSubKSP()`,
1930: `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`, `PCASMGetLocalSubdomains()`, `PCSetModifySubMatrices()`
1931: @*/
1932: PetscErrorCode PCASMGetLocalSubmatrices(PC pc, PetscInt *n, Mat *mat[])
1933: {
1934: PC_ASM *osm;
1935: PetscBool match;
1937: PetscFunctionBegin;
1939: if (n) PetscAssertPointer(n, 2);
1940: if (mat) PetscAssertPointer(mat, 3);
1941: PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call after KSPSetUp() or PCSetUp().");
1942: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCASM, &match));
1943: if (!match) {
1944: if (n) *n = 0;
1945: if (mat) *mat = NULL;
1946: } else {
1947: osm = (PC_ASM *)pc->data;
1948: if (n) *n = osm->n_local_true;
1949: if (mat) *mat = osm->pmat;
1950: }
1951: PetscFunctionReturn(PETSC_SUCCESS);
1952: }
1954: /*@
1955: PCASMSetDMSubdomains - Indicates whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible.
1957: Logically Collective
1959: Input Parameters:
1960: + pc - the preconditioner
1961: - flg - boolean indicating whether to use subdomains defined by the `DM`
1963: Options Database Key:
1964: . -pc_asm_dm_subdomains (true|false) - use subdomains defined by the `DM` with `DMCreateDomainDecomposition()`
1966: Level: intermediate
1968: Note:
1969: `PCASMSetTotalSubdomains()` and `PCASMSetOverlap()` take precedence over `PCASMSetDMSubdomains()`,
1970: so setting either of the first two effectively turns the latter off.
1972: Developer Note:
1973: This should be `PCASMSetUseDMSubdomains()`, similarly for the options database key
1975: .seealso: [](ch_ksp), `PCASM`, `PCASMGetDMSubdomains()`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`,
1976: `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`, `PCASMGetLocalSubdomains()`
1977: @*/
1978: PetscErrorCode PCASMSetDMSubdomains(PC pc, PetscBool flg)
1979: {
1980: PC_ASM *osm = (PC_ASM *)pc->data;
1981: PetscBool match;
1983: PetscFunctionBegin;
1986: PetscCheck(!pc->setupcalled, ((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONGSTATE, "Not for a setup PC.");
1987: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCASM, &match));
1988: if (match) osm->dm_subdomains = flg;
1989: PetscFunctionReturn(PETSC_SUCCESS);
1990: }
1992: /*@
1993: PCASMGetDMSubdomains - Returns flag indicating whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible.
1995: Not Collective
1997: Input Parameter:
1998: . pc - the preconditioner
2000: Output Parameter:
2001: . flg - boolean indicating whether to use subdomains defined by the `DM`
2003: Level: intermediate
2005: Developer Note:
2006: This should be `PCASMSetUseDMSubdomains()`
2008: .seealso: [](ch_ksp), `PCASM`, `PCASMSetDMSubdomains()`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`,
2009: `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`, `PCASMGetLocalSubdomains()`
2010: @*/
2011: PetscErrorCode PCASMGetDMSubdomains(PC pc, PetscBool *flg)
2012: {
2013: PC_ASM *osm = (PC_ASM *)pc->data;
2014: PetscBool match;
2016: PetscFunctionBegin;
2018: PetscAssertPointer(flg, 2);
2019: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCASM, &match));
2020: if (match) *flg = osm->dm_subdomains;
2021: else *flg = PETSC_FALSE;
2022: PetscFunctionReturn(PETSC_SUCCESS);
2023: }
2025: /*@
2026: PCASMGetSubMatType - Gets the matrix type used for `PCASM` subsolves, as a string.
2028: Not Collective
2030: Input Parameter:
2031: . pc - the `PC`
2033: Output Parameter:
2034: . sub_mat_type - name of matrix type
2036: Level: advanced
2038: .seealso: [](ch_ksp), `PCASM`, `PCASMSetSubMatType()`, `PCSetType()`, `VecSetType()`, `MatType`, `Mat`
2039: @*/
2040: PetscErrorCode PCASMGetSubMatType(PC pc, MatType *sub_mat_type)
2041: {
2042: PetscFunctionBegin;
2044: PetscTryMethod(pc, "PCASMGetSubMatType_C", (PC, MatType *), (pc, sub_mat_type));
2045: PetscFunctionReturn(PETSC_SUCCESS);
2046: }
2048: /*@
2049: PCASMSetSubMatType - Set the type of matrix used for `PCASM` subsolves
2051: Collective
2053: Input Parameters:
2054: + pc - the `PC` object
2055: - sub_mat_type - the `MatType`
2057: Options Database Key:
2058: . -pc_asm_sub_mat_type sub_mat_type - Sets the matrix type used for subsolves, for example, seqaijviennacl.
2059: If you specify a base name like aijviennacl, the corresponding sequential type is assumed.
2061: Note:
2062: See `MatType` for available types
2064: Level: advanced
2066: .seealso: [](ch_ksp), `PCASM`, `PCASMGetSubMatType()`, `PCSetType()`, `VecSetType()`, `MatType`, `Mat`
2067: @*/
2068: PetscErrorCode PCASMSetSubMatType(PC pc, MatType sub_mat_type)
2069: {
2070: PetscFunctionBegin;
2072: PetscTryMethod(pc, "PCASMSetSubMatType_C", (PC, MatType), (pc, sub_mat_type));
2073: PetscFunctionReturn(PETSC_SUCCESS);
2074: }