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: VecType vtype;
515: PetscInt i, m, N;
516: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
518: PetscFunctionBegin;
519: PetscCheck(osm->n_local_true <= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not yet implemented");
520: PetscCheck(osm->type != PC_ASM_WEIGHTED || osm->scaling, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Call PCASMWeightedSetScaling() after PCSetUp() before applying PC_ASM_WEIGHTED");
521: /*
522: support for limiting the restriction or interpolation to only local
523: subdomain values (leaving the other values 0).
524: */
525: if ((!transpose && !(type & PC_ASM_RESTRICT)) || (transpose && !(type & PC_ASM_INTERPOLATE))) {
526: forward = SCATTER_FORWARD_LOCAL;
527: /* have to zero the work RHS since scatter may leave some slots empty */
528: PetscCall(VecSet(osm->lx, 0.0));
529: }
530: if ((!transpose && !(type & PC_ASM_INTERPOLATE)) || (transpose && !(type & PC_ASM_RESTRICT))) reverse = SCATTER_REVERSE_LOCAL;
531: PetscCall(VecGetLocalSize(osm->x[0], &m));
532: PetscCall(MatGetSize(X, NULL, &N));
533: PetscCall(VecGetType(osm->x[0], &vtype));
534: PetscCall(MatCreateDenseFromVecType(PETSC_COMM_SELF, vtype, m, N, m, N, PETSC_DECIDE, NULL, &Z));
536: 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]);
537: /* zero the global and the local solutions */
538: PetscCall(MatZeroEntries(Y));
539: PetscCall(VecSet(osm->ly, 0.0));
541: for (i = 0; i < N; ++i) {
542: PetscCall(MatDenseGetColumnVecRead(X, i, &x));
543: /* copy the global RHS to local RHS including the ghost nodes */
544: PetscCall(VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward));
545: PetscCall(VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward));
546: PetscCall(MatDenseRestoreColumnVecRead(X, i, &x));
548: PetscCall(MatDenseGetColumnVecWrite(Z, i, &x));
549: /* restrict local RHS to the overlapping 0-block RHS */
550: PetscCall(VecScatterBegin(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward));
551: PetscCall(VecScatterEnd(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward));
552: PetscCall(MatDenseRestoreColumnVecWrite(Z, i, &x));
553: }
554: PetscCall(MatCreateDenseFromVecType(PETSC_COMM_SELF, vtype, m, N, m, N, PETSC_DECIDE, NULL, &W));
555: /* solve the overlapping 0-block */
556: if (!transpose) {
557: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[0], Z, W, 0));
558: PetscCall(KSPMatSolve(osm->ksp[0], Z, W));
559: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[0], Z, W, 0));
560: } else {
561: if (osm->type == PC_ASM_WEIGHTED) PetscCall(MatDiagonalScale(Z, osm->scaling[0], NULL));
562: PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, osm->ksp[0], Z, W, 0));
563: PetscCall(KSPMatSolveTranspose(osm->ksp[0], Z, W));
564: PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, osm->ksp[0], Z, W, 0));
565: }
566: PetscCall(KSPCheckMatSolve(osm->ksp[0], pc, W));
567: if (!transpose && osm->type == PC_ASM_WEIGHTED) PetscCall(MatDiagonalScale(W, osm->scaling[0], NULL));
568: PetscCall(MatDestroy(&Z));
570: for (i = 0; i < N; ++i) {
571: PetscCall(VecSet(osm->ly, 0.0));
572: PetscCall(MatDenseGetColumnVecRead(W, i, &x));
573: 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) */
574: PetscCall(VecScatterBegin(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward));
575: PetscCall(VecScatterEnd(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward));
576: } else { /* interpolate the overlapping 0-block solution to the local solution */
577: PetscCall(VecScatterBegin(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse));
578: PetscCall(VecScatterEnd(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse));
579: }
580: PetscCall(MatDenseRestoreColumnVecRead(W, i, &x));
582: PetscCall(MatDenseGetColumnVecWrite(Y, i, &x));
583: /* add the local solution to the global solution including the ghost nodes */
584: PetscCall(VecScatterBegin(osm->restriction, osm->ly, x, ADD_VALUES, reverse));
585: PetscCall(VecScatterEnd(osm->restriction, osm->ly, x, ADD_VALUES, reverse));
586: PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &x));
587: }
588: PetscCall(MatDestroy(&W));
589: PetscFunctionReturn(PETSC_SUCCESS);
590: }
592: static PetscErrorCode PCMatApply_ASM(PC pc, Mat X, Mat Y)
593: {
594: PetscFunctionBegin;
595: PetscCall(PCMatApply_ASM_Private(pc, X, Y, PETSC_FALSE));
596: PetscFunctionReturn(PETSC_SUCCESS);
597: }
599: static PetscErrorCode PCMatApplyTranspose_ASM(PC pc, Mat X, Mat Y)
600: {
601: PetscFunctionBegin;
602: PetscCall(PCMatApply_ASM_Private(pc, X, Y, PETSC_TRUE));
603: PetscFunctionReturn(PETSC_SUCCESS);
604: }
606: static PetscErrorCode PCApplyTranspose_ASM(PC pc, Vec x, Vec y)
607: {
608: PC_ASM *osm = (PC_ASM *)pc->data;
609: 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 */
610: PetscInt i, n_local_true = osm->n_local_true;
611: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
613: PetscFunctionBegin;
614: PetscCheck(osm->n_local_true <= 1 || osm->loctype == PC_COMPOSITE_ADDITIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not yet implemented");
615: PetscCheck(osm->type != PC_ASM_WEIGHTED || osm->scaling, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Call PCASMWeightedSetScaling() after PCSetUp() before applying PC_ASM_WEIGHTED");
616: /*
617: Support for limiting the restriction or interpolation to only local
618: subdomain values (leaving the other values 0).
620: Note: these are reversed from the PCApply_ASM() because we are applying the
621: transpose of the three terms
622: */
624: if (!(type & PC_ASM_INTERPOLATE)) {
625: forward = SCATTER_FORWARD_LOCAL;
626: /* have to zero the work RHS since scatter may leave some slots empty */
627: PetscCall(VecSet(osm->lx, 0.0));
628: }
629: if (!(type & PC_ASM_RESTRICT)) reverse = SCATTER_REVERSE_LOCAL;
631: /* zero the global and the local solutions */
632: PetscCall(VecSet(y, 0.0));
633: PetscCall(VecSet(osm->ly, 0.0));
635: /* Copy the global RHS to local RHS including the ghost nodes */
636: PetscCall(VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward));
637: PetscCall(VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward));
639: /* Restrict local RHS to the overlapping 0-block RHS */
640: PetscCall(VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward));
641: PetscCall(VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward));
643: /* do the local solves */
644: for (i = 0; i < n_local_true; ++i) {
645: /* solve the overlapping i-block */
646: if (osm->type == PC_ASM_WEIGHTED) PetscCall(VecPointwiseMult(osm->x[i], osm->scaling[i], osm->x[i]));
647: PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0));
648: PetscCall(KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i]));
649: PetscCall(KSPCheckSolve(osm->ksp[i], pc, osm->y[i]));
650: PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0));
652: if (osm->lprolongation && !(type & PC_ASM_RESTRICT)) { /* interpolate the non-overlapping i-block solution to the local solution */
653: PetscCall(VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward));
654: PetscCall(VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward));
655: } else { /* interpolate the overlapping i-block solution to the local solution */
656: PetscCall(VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse));
657: PetscCall(VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse));
658: }
660: if (i < n_local_true - 1) {
661: /* Restrict local RHS to the overlapping (i+1)-block RHS */
662: PetscCall(VecScatterBegin(osm->lrestriction[i + 1], osm->lx, osm->x[i + 1], INSERT_VALUES, forward));
663: PetscCall(VecScatterEnd(osm->lrestriction[i + 1], osm->lx, osm->x[i + 1], INSERT_VALUES, forward));
664: }
665: }
666: /* Add the local solution to the global solution including the ghost nodes */
667: PetscCall(VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse));
668: PetscCall(VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse));
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: static PetscErrorCode PCReset_ASM(PC pc)
673: {
674: PC_ASM *osm = (PC_ASM *)pc->data;
676: PetscFunctionBegin;
677: if (osm->scaling) {
678: for (PetscInt i = 0; i < osm->n_local_true; i++) PetscCall(VecDestroy(&osm->scaling[i]));
679: PetscCall(PetscFree(osm->scaling));
680: }
681: if (osm->ksp) {
682: for (PetscInt i = 0; i < osm->n_local_true; i++) PetscCall(KSPReset(osm->ksp[i]));
683: }
684: if (osm->pmat) {
685: if (osm->n_local_true > 0) PetscCall(MatDestroySubMatrices(osm->n_local_true, &osm->pmat));
686: }
687: if (osm->lrestriction) {
688: PetscCall(VecScatterDestroy(&osm->restriction));
689: for (PetscInt i = 0; i < osm->n_local_true; i++) {
690: PetscCall(VecScatterDestroy(&osm->lrestriction[i]));
691: if (osm->lprolongation) PetscCall(VecScatterDestroy(&osm->lprolongation[i]));
692: PetscCall(VecDestroy(&osm->x[i]));
693: PetscCall(VecDestroy(&osm->y[i]));
694: }
695: PetscCall(PetscFree(osm->lrestriction));
696: PetscCall(PetscFree(osm->lprolongation));
697: PetscCall(PetscFree(osm->x));
698: PetscCall(PetscFree(osm->y));
699: }
700: PetscCall(PCASMDestroySubdomains(osm->n_local_true, &osm->is, &osm->is_local));
701: PetscCall(ISDestroy(&osm->lis));
702: PetscCall(VecDestroy(&osm->lx));
703: PetscCall(VecDestroy(&osm->ly));
704: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(MatDestroyMatrices(osm->n_local_true, &osm->lmats));
706: PetscCall(PetscFree(osm->sub_mat_type));
708: osm->is = NULL;
709: osm->is_local = NULL;
710: PetscFunctionReturn(PETSC_SUCCESS);
711: }
713: static PetscErrorCode PCDestroy_ASM(PC pc)
714: {
715: PC_ASM *osm = (PC_ASM *)pc->data;
717: PetscFunctionBegin;
718: PetscCall(PCReset_ASM(pc));
719: if (osm->ksp) {
720: for (PetscInt i = 0; i < osm->n_local_true; i++) PetscCall(KSPDestroy(&osm->ksp[i]));
721: PetscCall(PetscFree(osm->ksp));
722: }
723: PetscCall(PetscFree(pc->data));
725: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetLocalSubdomains_C", NULL));
726: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetTotalSubdomains_C", NULL));
727: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetOverlap_C", NULL));
728: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetType_C", NULL));
729: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetType_C", NULL));
730: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMWeightedSetScaling_C", NULL));
731: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMWeightedSetComputeScaling_C", NULL));
732: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetLocalType_C", NULL));
733: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetLocalType_C", NULL));
734: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetSortIndices_C", NULL));
735: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetSubKSP_C", NULL));
736: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetSubMatType_C", NULL));
737: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetSubMatType_C", NULL));
738: PetscFunctionReturn(PETSC_SUCCESS);
739: }
741: static PetscErrorCode PCSetFromOptions_ASM(PC pc, PetscOptionItems PetscOptionsObject)
742: {
743: PC_ASM *osm = (PC_ASM *)pc->data;
744: PetscInt blocks, ovl;
745: PetscBool flg;
746: PCASMType asmtype;
747: PCCompositeType loctype;
748: char sub_mat_type[256];
750: PetscFunctionBegin;
751: PetscOptionsHeadBegin(PetscOptionsObject, "Additive Schwarz options");
752: PetscCall(PetscOptionsBool("-pc_asm_dm_subdomains", "Use DMCreateDomainDecomposition() to define subdomains", "PCASMSetDMSubdomains", osm->dm_subdomains, &osm->dm_subdomains, &flg));
753: PetscCall(PetscOptionsInt("-pc_asm_blocks", "Number of subdomains", "PCASMSetTotalSubdomains", osm->n, &blocks, &flg));
754: if (flg) {
755: PetscCall(PCASMSetTotalSubdomains(pc, blocks, NULL, NULL));
756: osm->dm_subdomains = PETSC_FALSE;
757: }
758: PetscCall(PetscOptionsInt("-pc_asm_local_blocks", "Number of local subdomains", "PCASMSetLocalSubdomains", osm->n_local_true, &blocks, &flg));
759: if (flg) {
760: PetscCall(PCASMSetLocalSubdomains(pc, blocks, NULL, NULL));
761: osm->dm_subdomains = PETSC_FALSE;
762: }
763: PetscCall(PetscOptionsInt("-pc_asm_overlap", "Number of grid points overlap", "PCASMSetOverlap", osm->overlap, &ovl, &flg));
764: if (flg) {
765: PetscCall(PCASMSetOverlap(pc, ovl));
766: osm->dm_subdomains = PETSC_FALSE;
767: }
768: flg = PETSC_FALSE;
769: PetscCall(PetscOptionsEnum("-pc_asm_type", "Type of restriction/extension", "PCASMSetType", PCASMTypes, (PetscEnum)osm->type, (PetscEnum *)&asmtype, &flg));
770: if (flg) PetscCall(PCASMSetType(pc, asmtype));
771: flg = PETSC_FALSE;
772: PetscCall(PetscOptionsEnum("-pc_asm_local_type", "Type of local solver composition", "PCASMSetLocalType", PCCompositeTypes, (PetscEnum)osm->loctype, (PetscEnum *)&loctype, &flg));
773: if (flg) PetscCall(PCASMSetLocalType(pc, loctype));
774: PetscCall(PetscOptionsFList("-pc_asm_sub_mat_type", "Subsolve Matrix Type", "PCASMSetSubMatType", MatList, NULL, sub_mat_type, sizeof(sub_mat_type), &flg));
775: if (flg) PetscCall(PCASMSetSubMatType(pc, sub_mat_type));
776: PetscOptionsHeadEnd();
777: PetscFunctionReturn(PETSC_SUCCESS);
778: }
780: static PetscErrorCode PCASMSetLocalSubdomains_ASM(PC pc, PetscInt n, IS is[], IS is_local[])
781: {
782: PC_ASM *osm = (PC_ASM *)pc->data;
784: PetscFunctionBegin;
785: PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Each process must have 1 or more blocks, n = %" PetscInt_FMT, n);
786: PetscCheck(!pc->setupcalled || (n == osm->n_local_true && !is), PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCASMSetLocalSubdomains() should be called before calling PCSetUp().");
788: if (!pc->setupcalled) {
789: if (is) {
790: for (PetscInt i = 0; i < n; i++) PetscCall(PetscObjectReference((PetscObject)is[i]));
791: }
792: if (is_local) {
793: for (PetscInt i = 0; i < n; i++) PetscCall(PetscObjectReference((PetscObject)is_local[i]));
794: }
795: PetscCall(PCASMDestroySubdomains(osm->n_local_true, &osm->is, &osm->is_local));
797: if (osm->ksp && osm->n_local_true != n) {
798: for (PetscInt i = 0; i < osm->n_local_true; i++) PetscCall(KSPDestroy(&osm->ksp[i]));
799: PetscCall(PetscFree(osm->ksp));
800: }
802: osm->n_local_true = n;
803: osm->is = NULL;
804: osm->is_local = NULL;
805: if (is) {
806: PetscCall(PetscMalloc1(n, &osm->is));
807: for (PetscInt i = 0; i < n; i++) osm->is[i] = is[i];
808: /* Flag indicating that the user has set overlapping subdomains so PCASM should not increase their size. */
809: osm->overlap = -1;
810: }
811: if (is_local) {
812: PetscCall(PetscMalloc1(n, &osm->is_local));
813: for (PetscInt i = 0; i < n; i++) osm->is_local[i] = is_local[i];
814: if (!is) {
815: PetscCall(PetscMalloc1(osm->n_local_true, &osm->is));
816: for (PetscInt i = 0; i < osm->n_local_true; i++) {
817: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
818: PetscCall(ISDuplicate(osm->is_local[i], &osm->is[i]));
819: PetscCall(ISCopy(osm->is_local[i], osm->is[i]));
820: } else {
821: PetscCall(PetscObjectReference((PetscObject)osm->is_local[i]));
822: osm->is[i] = osm->is_local[i];
823: }
824: }
825: }
826: }
827: }
828: PetscFunctionReturn(PETSC_SUCCESS);
829: }
831: static PetscErrorCode PCASMSetTotalSubdomains_ASM(PC pc, PetscInt N, IS *is, IS *is_local)
832: {
833: PC_ASM *osm = (PC_ASM *)pc->data;
834: PetscMPIInt rank, size;
835: PetscInt n;
837: PetscFunctionBegin;
838: PetscCheck(N >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Number of total blocks must be > 0, N = %" PetscInt_FMT, N);
839: PetscCheck(!is && !is_local, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Use PCASMSetLocalSubdomains() to set specific index sets, they cannot be set globally yet.");
841: /*
842: Split the subdomains equally among all processors
843: */
844: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
845: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
846: n = N / size + ((N % size) > rank);
847: 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);
848: PetscCheck(!pc->setupcalled || n == osm->n_local_true, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "PCASMSetTotalSubdomains() should be called before PCSetUp().");
849: if (!pc->setupcalled) {
850: PetscCall(PCASMDestroySubdomains(osm->n_local_true, &osm->is, &osm->is_local));
852: osm->n_local_true = n;
853: osm->is = NULL;
854: osm->is_local = NULL;
855: }
856: PetscFunctionReturn(PETSC_SUCCESS);
857: }
859: static PetscErrorCode PCASMSetOverlap_ASM(PC pc, PetscInt ovl)
860: {
861: PC_ASM *osm = (PC_ASM *)pc->data;
863: PetscFunctionBegin;
864: PetscCheck(ovl >= 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Negative overlap value requested");
865: PetscCheck(!pc->setupcalled || ovl == osm->overlap, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCASMSetOverlap() should be called before PCSetUp().");
866: if (!pc->setupcalled) osm->overlap = ovl;
867: PetscFunctionReturn(PETSC_SUCCESS);
868: }
870: static PetscErrorCode PCASMSetType_ASM(PC pc, PCASMType type)
871: {
872: PC_ASM *osm = (PC_ASM *)pc->data;
874: PetscFunctionBegin;
875: PetscCheck(type != PC_ASM_WEIGHTED || osm->loctype == PC_COMPOSITE_ADDITIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PC_ASM_WEIGHTED requires additive local composition");
876: osm->type = type;
877: osm->type_set = PETSC_TRUE;
878: PetscFunctionReturn(PETSC_SUCCESS);
879: }
881: static PetscErrorCode PCASMGetType_ASM(PC pc, PCASMType *type)
882: {
883: PC_ASM *osm = (PC_ASM *)pc->data;
885: PetscFunctionBegin;
886: *type = osm->type;
887: PetscFunctionReturn(PETSC_SUCCESS);
888: }
890: static PetscErrorCode PCASMWeightedSetComputeScaling_ASM(PC pc, PCASMWeightedComputeScalingFn *fn, PetscCtx ctx)
891: {
892: PC_ASM *osm = (PC_ASM *)pc->data;
894: PetscFunctionBegin;
895: osm->computescaling = fn;
896: osm->computescalingctx = ctx;
897: PetscFunctionReturn(PETSC_SUCCESS);
898: }
900: static PetscErrorCode PCASMWeightedSetScaling_ASM(PC pc, PetscInt n, Vec scaling[])
901: {
902: PC_ASM *osm = (PC_ASM *)pc->data;
903: Vec *newscaling;
904: VecType type;
905: PetscInt m, nvec;
906: PetscMPIInt size;
907: PetscBool match;
909: PetscFunctionBegin;
910: PetscCheck(osm->x, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Call PCSetUp() before PCASMWeightedSetScaling() so that the subdomain sizes and ordering are final");
911: 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);
912: for (PetscInt i = 0; i < n; i++) {
913: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)scaling[i]), &size));
914: 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);
915: PetscCall(VecGetSize(scaling[i], &nvec));
916: PetscCall(VecGetSize(osm->x[i], &m));
917: 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);
918: PetscCall(VecGetType(osm->x[i], &type));
919: PetscCall(PetscObjectTypeCompare((PetscObject)scaling[i], type, &match));
920: 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);
921: }
922: PetscCall(PetscMalloc1(n, &newscaling));
923: for (PetscInt i = 0; i < n; i++) {
924: PetscCall(PetscObjectReference((PetscObject)scaling[i]));
925: newscaling[i] = scaling[i];
926: }
927: if (osm->scaling) {
928: for (PetscInt i = 0; i < n; i++) PetscCall(VecDestroy(&osm->scaling[i]));
929: PetscCall(PetscFree(osm->scaling));
930: }
931: osm->scaling = newscaling;
932: PetscFunctionReturn(PETSC_SUCCESS);
933: }
935: static PetscErrorCode PCASMSetLocalType_ASM(PC pc, PCCompositeType type)
936: {
937: PC_ASM *osm = (PC_ASM *)pc->data;
939: PetscFunctionBegin;
940: PetscCheck(type == PC_COMPOSITE_ADDITIVE || type == PC_COMPOSITE_MULTIPLICATIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Only supports additive or multiplicative as the local type");
941: PetscCheck(osm->type != PC_ASM_WEIGHTED || type == PC_COMPOSITE_ADDITIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "PC_ASM_WEIGHTED requires additive local composition");
942: osm->loctype = type;
943: PetscFunctionReturn(PETSC_SUCCESS);
944: }
946: static PetscErrorCode PCASMGetLocalType_ASM(PC pc, PCCompositeType *type)
947: {
948: PC_ASM *osm = (PC_ASM *)pc->data;
950: PetscFunctionBegin;
951: *type = osm->loctype;
952: PetscFunctionReturn(PETSC_SUCCESS);
953: }
955: static PetscErrorCode PCASMSetSortIndices_ASM(PC pc, PetscBool doSort)
956: {
957: PC_ASM *osm = (PC_ASM *)pc->data;
959: PetscFunctionBegin;
960: osm->sort_indices = doSort;
961: PetscFunctionReturn(PETSC_SUCCESS);
962: }
964: static PetscErrorCode PCASMGetSubKSP_ASM(PC pc, PetscInt *n_local, PetscInt *first_local, KSP **ksp)
965: {
966: PC_ASM *osm = (PC_ASM *)pc->data;
968: PetscFunctionBegin;
969: 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");
971: if (n_local) *n_local = osm->n_local_true;
972: if (first_local) {
973: PetscCallMPI(MPI_Scan(&osm->n_local_true, first_local, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)pc)));
974: *first_local -= osm->n_local_true;
975: }
976: if (ksp) *ksp = osm->ksp;
977: PetscFunctionReturn(PETSC_SUCCESS);
978: }
980: static PetscErrorCode PCASMGetSubMatType_ASM(PC pc, MatType *sub_mat_type)
981: {
982: PC_ASM *osm = (PC_ASM *)pc->data;
984: PetscFunctionBegin;
986: PetscAssertPointer(sub_mat_type, 2);
987: *sub_mat_type = osm->sub_mat_type;
988: PetscFunctionReturn(PETSC_SUCCESS);
989: }
991: static PetscErrorCode PCASMSetSubMatType_ASM(PC pc, MatType sub_mat_type)
992: {
993: PC_ASM *osm = (PC_ASM *)pc->data;
995: PetscFunctionBegin;
997: PetscCall(PetscFree(osm->sub_mat_type));
998: PetscCall(PetscStrallocpy(sub_mat_type, (char **)&osm->sub_mat_type));
999: PetscFunctionReturn(PETSC_SUCCESS);
1000: }
1002: /*@
1003: PCASMSetLocalSubdomains - Sets the local subdomains (for this processor only) for the additive Schwarz preconditioner `PCASM`.
1005: Collective
1007: Input Parameters:
1008: + pc - the preconditioner context
1009: . n - the number of subdomains for this processor (default value = 1)
1010: . is - the index set that defines the subdomains for this processor (or `NULL` for PETSc to determine subdomains)
1011: the values of the `is` array are copied so you can free the array (not the `IS` in the array) after this call
1012: - is_local - the index sets that define the local part of the subdomains for this processor, not used unless `PCASMType` is `PC_ASM_RESTRICT`
1013: (or `NULL` to not provide these). The values of the `is_local` array are copied so you can free the array
1014: (not the `IS` in the array) after this call
1016: Options Database Key:
1017: . -pc_asm_local_blocks blks - Sets number of local blocks
1019: Level: advanced
1021: Notes:
1022: The `IS` numbering is in the parallel, global numbering of the vector for both `is` and `is_local`
1024: By default the `PCASM` preconditioner uses 1 block per processor.
1026: Use `PCASMSetTotalSubdomains()` to set the subdomains for all processors.
1028: If `is_local` is provided and `PCASMType` is `PC_ASM_RESTRICT` then the solution only over the `is_local` region is interpolated
1029: back to form the global solution (this is the standard restricted additive Schwarz method, RASM)
1031: If `is_local` is provided and `PCASMType` is `PC_ASM_INTERPOLATE` or `PC_ASM_NONE` then an error is generated since there is
1032: no code to handle that case.
1034: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`, `PCASMGetSubKSP()`,
1035: `PCASMCreateSubdomains2D()`, `PCASMGetLocalSubdomains()`, `PCASMType`, `PCASMSetType()`, `PCGASM`
1036: @*/
1037: PetscErrorCode PCASMSetLocalSubdomains(PC pc, PetscInt n, IS is[], IS is_local[])
1038: {
1039: PetscFunctionBegin;
1041: PetscTryMethod(pc, "PCASMSetLocalSubdomains_C", (PC, PetscInt, IS[], IS[]), (pc, n, is, is_local));
1042: PetscFunctionReturn(PETSC_SUCCESS);
1043: }
1045: /*@
1046: PCASMSetTotalSubdomains - Sets the subdomains for all processors for the
1047: additive Schwarz preconditioner, `PCASM`.
1049: Collective, all MPI ranks must pass in the same array of `IS`
1051: Input Parameters:
1052: + pc - the preconditioner context
1053: . N - the number of subdomains for all processors
1054: . is - the index sets that define the subdomains for all processors (or `NULL` to ask PETSc to determine the subdomains)
1055: the values of the `is` array are copied so you can free the array (not the `IS` in the array) after this call
1056: - is_local - the index sets that define the local part of the subdomains for this processor (or `NULL` to not provide this information)
1057: The values of the `is_local` array are copied so you can free the array (not the `IS` in the array) after this call
1059: Options Database Key:
1060: . -pc_asm_blocks blks - Sets total blocks
1062: Level: advanced
1064: Notes:
1065: Currently you cannot use this to set the actual subdomains with the argument `is` or `is_local`.
1067: By default the `PCASM` preconditioner uses 1 block per processor.
1069: These index sets cannot be destroyed until after completion of the
1070: linear solves for which the `PCASM` preconditioner is being used.
1072: Use `PCASMSetLocalSubdomains()` to set local subdomains.
1074: The `IS` numbering is in the parallel, global numbering of the vector for both is and is_local
1076: .seealso: [](ch_ksp), `PCASM`, `PCASMSetLocalSubdomains()`, `PCASMSetOverlap()`, `PCASMGetSubKSP()`,
1077: `PCASMCreateSubdomains2D()`, `PCGASM`
1078: @*/
1079: PetscErrorCode PCASMSetTotalSubdomains(PC pc, PetscInt N, IS is[], IS is_local[])
1080: {
1081: PetscFunctionBegin;
1083: PetscTryMethod(pc, "PCASMSetTotalSubdomains_C", (PC, PetscInt, IS[], IS[]), (pc, N, is, is_local));
1084: PetscFunctionReturn(PETSC_SUCCESS);
1085: }
1087: /*@
1088: PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
1089: additive Schwarz preconditioner, `PCASM`.
1091: Logically Collective
1093: Input Parameters:
1094: + pc - the preconditioner context
1095: - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1)
1097: Options Database Key:
1098: . -pc_asm_overlap ovl - Sets overlap
1100: Level: intermediate
1102: Notes:
1103: By default the `PCASM` preconditioner uses 1 block per processor. To use
1104: multiple blocks per perocessor, see `PCASMSetTotalSubdomains()` and
1105: `PCASMSetLocalSubdomains()` (and the option -pc_asm_blocks <blks>).
1107: The overlap defaults to 1, so if one desires that no additional
1108: overlap be computed beyond what may have been set with a call to
1109: `PCASMSetTotalSubdomains()` or `PCASMSetLocalSubdomains()`, then ovl
1110: must be set to be 0. In particular, if one does not explicitly set
1111: the subdomains an application code, then all overlap would be computed
1112: internally by PETSc, and using an overlap of 0 would result in an `PCASM`
1113: variant that is equivalent to the block Jacobi preconditioner.
1115: The default algorithm used by PETSc to increase overlap is fast, but not scalable,
1116: use the option -mat_increase_overlap_scalable when the problem and number of processes is large.
1118: One can define initial index sets with any overlap via
1119: `PCASMSetLocalSubdomains()`; the routine
1120: `PCASMSetOverlap()` merely allows PETSc to extend that overlap further
1121: if desired.
1123: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetLocalSubdomains()`, `PCASMGetSubKSP()`,
1124: `PCASMCreateSubdomains2D()`, `PCASMGetLocalSubdomains()`, `MatIncreaseOverlap()`, `PCGASM`
1125: @*/
1126: PetscErrorCode PCASMSetOverlap(PC pc, PetscInt ovl)
1127: {
1128: PetscFunctionBegin;
1131: PetscTryMethod(pc, "PCASMSetOverlap_C", (PC, PetscInt), (pc, ovl));
1132: PetscFunctionReturn(PETSC_SUCCESS);
1133: }
1135: /*@
1136: PCASMSetType - Sets the type of restriction and interpolation used
1137: for local problems in the additive Schwarz method, `PCASM`.
1139: Logically Collective
1141: Input Parameters:
1142: + pc - the preconditioner context
1143: - type - variant of `PCASM`, one of
1144: .vb
1145: PC_ASM_NONE - local processor restriction and interpolation
1146: PC_ASM_RESTRICT - full restriction, local processor interpolation (default)
1147: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1148: PC_ASM_BASIC - full interpolation and restriction
1149: PC_ASM_WEIGHTED - full restriction and interpolation with user-provided diagonal scaling
1150: .ve
1152: Options Database Key:
1153: . -pc_asm_type (none|restrict|interpolate|basic|weighted) - Sets `PCASMType`
1155: Level: intermediate
1157: Note:
1158: if the is_local arguments are passed to `PCASMSetLocalSubdomains()` then they are used when `PC_ASM_RESTRICT` has been selected
1159: to limit the local processor interpolation. `PC_ASM_WEIGHTED` ignores these inner index sets and
1160: uses the weights supplied with `PCASMWeightedSetScaling()` instead. Weighted ASM requires additive local composition.
1162: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMGetSubKSP()`,
1163: `PCASMCreateSubdomains2D()`, `PCASMType`, `PCASMWeightedSetScaling()`, `PCASMSetLocalType()`, `PCASMGetLocalType()`, `PCGASM`
1164: @*/
1165: PetscErrorCode PCASMSetType(PC pc, PCASMType type)
1166: {
1167: PetscFunctionBegin;
1170: PetscTryMethod(pc, "PCASMSetType_C", (PC, PCASMType), (pc, type));
1171: PetscFunctionReturn(PETSC_SUCCESS);
1172: }
1174: /*@
1175: PCASMWeightedSetComputeScaling - Sets a callback to compute `PC_ASM_WEIGHTED` scaling during `PCSetUp()`.
1177: Logically Collective
1179: Input Parameters:
1180: + pc - the `PCASM` preconditioner
1181: . fn - function to fill each local scaling vector, or `NULL` to disable the callback
1182: - ctx - function context passed to `fn`
1184: Level: intermediate
1186: Notes:
1187: Register before setup and select `PC_ASM_WEIGHTED`.
1189: Whenever `PCSetUp()` rebuilds weighted ASM,
1190: `fn` is called once per local subdomain, after overlap expansion and index sorting, with a vector
1191: of the correct size and type.
1193: No explicit setup or vector allocation is needed by the caller.
1195: The callback overwrites any existing weights, including those supplied directly by `PCASMWeightedSetScaling()`.
1197: Registration does not trigger setup. Passing `NULL` leaves the current weights in place.
1199: The callback and context survive `PCReset()`. The caller owns `ctx` and must keep it valid
1200: until the callback is replaced or disabled, or the preconditioner is destroyed.
1202: Fortran Note:
1203: `fn` is a subroutine with arguments `(pc, local, scaling, ctx, ierr)`. Pass `PETSC_NULL_FUNCTION` to disable the callback.
1205: .seealso: [](ch_ksp), `PCASM`, `PCASMWeightedComputeScalingFn`, `PCASMWeightedSetScaling()`, `PCASMWeightedGetScaling()`, `PCASMSetType()`
1206: @*/
1207: PetscErrorCode PCASMWeightedSetComputeScaling(PC pc, PCASMWeightedComputeScalingFn *fn, PetscCtx ctx)
1208: {
1209: PetscFunctionBegin;
1211: PetscTryMethod(pc, "PCASMWeightedSetComputeScaling_C", (PC, PCASMWeightedComputeScalingFn *, PetscCtx), (pc, fn, ctx));
1212: PetscFunctionReturn(PETSC_SUCCESS);
1213: }
1215: /*@
1216: PCASMWeightedSetScaling - Sets the diagonal weights for the overlapping local corrections in weighted additive Schwarz.
1218: Logically Collective
1220: Input Parameters:
1221: + pc - the `PCASM` preconditioner
1222: . n - the number of local subdomains
1223: - scaling - one local scaling `Vec` per overlapping subdomain
1225: Level: intermediate
1227: Notes:
1228: Call `PCSetUp()` before this routine, then use `PCASMGetLocalSubdomains()` to obtain the final
1229: overlapping index sets. Each vector must have a single-process communicator and the same length
1230: and ordering as its corresponding overlapping index set, after overlap expansion and sorting.
1231: Its vector type must match the local solver vectors. Compatible vectors can be created with
1232: `MatCreateVecs()` from the corresponding matrix returned by `PCASMGetLocalSubmatrices()`.
1233: These are local subdomain vectors, not vectors in the parallel global layout.
1234: Select `PC_ASM_WEIGHTED` with `PCASMSetType()` or `-pc_asm_type weighted` to use the weights.
1235: The weights are ignored by the other ASM types.
1237: With restriction operators $R_i$, local solvers $A_i^{-1}$, and $D_i = \text{diag}(scaling[i])$, the action is
1238: $B = \sum_i R_i^T D_i A_i^{-1} R_i$.
1239: PETSc uses the supplied weights as-is, without checking whether they are real, non-negative,
1240: or satisfy $\sum_i R_i^T D_i R_i = I$.
1241: The PC increments the reference count of the vectors but does not copy them.
1242: `PCReset()` discards the weights along with the subdomains.
1243: Alternatively, use `PCASMWeightedSetComputeScaling()` to fill internally-created vectors during `PCSetUp()`.
1245: Example Usage:
1246: .vb
1247: KSPGetPC(ksp, &pc);
1248: PCASMSetType(pc, PC_ASM_WEIGHTED);
1249: PCSetUp(pc);
1250: PCASMGetLocalSubmatrices(pc, &n, &submat);
1251: PCASMGetLocalSubdomains(pc, NULL, &is, NULL);
1252: for (i = 0; i < n; i++) {
1253: MatCreateVecs(submat[i], &scaling[i], NULL); // a Vec of the right size and type
1254: // fill scaling[i] in the local ordering of is[i]
1255: }
1256: PCASMWeightedSetScaling(pc, n, scaling);
1257: .ve
1259: .seealso: [](ch_ksp), `PCASM`, `PCASMType`, `PCASMSetType()`, `PCASMWeightedGetScaling()`, `PCASMGetLocalSubdomains()`, `PCASMGetLocalSubmatrices()`, `PCASMSetLocalSubdomains()`, `PCASMSetLocalType()`
1260: @*/
1261: PetscErrorCode PCASMWeightedSetScaling(PC pc, PetscInt n, Vec scaling[])
1262: {
1263: PetscFunctionBegin;
1265: PetscCheck(n >= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Number of scaling vectors must be nonnegative");
1266: if (n) PetscAssertPointer(scaling, 3);
1268: PetscTryMethod(pc, "PCASMWeightedSetScaling_C", (PC, PetscInt, Vec[]), (pc, n, scaling));
1269: PetscFunctionReturn(PETSC_SUCCESS);
1270: }
1272: /*@
1273: PCASMWeightedGetScaling - Gets the diagonal weights supplied with `PCASMWeightedSetScaling()` or computed by the function provided with `PCASMWeightedSetComputeScaling()`.
1275: Not Collective
1277: Input Parameter:
1278: . pc - the `PCASM` preconditioner
1280: Output Parameters:
1281: + n - if requested, the number of local subdomains for this processor, or zero if no weights have been supplied
1282: - scaling - if requested, the local scaling `Vec`, or `NULL` if none have been supplied
1284: Level: intermediate
1286: Note:
1287: The returned array and its vectors are owned by `pc`; do not free or destroy them. They
1288: are released by `PCReset()` and `PCDestroy()`, and replaced by a further call to
1289: `PCASMWeightedSetScaling()`.
1291: Fortran Note:
1292: Declare `scaling` as `Vec, pointer :: scaling(:)`. It is always returned and is disassociated when no weights
1293: have been supplied; `n` may be `PETSC_NULL_INTEGER`. There is no restore routine, and the pointer must not be
1294: used after `PCReset()`, `PCDestroy()`, or a further call to `PCASMWeightedSetScaling()`.
1296: .seealso: [](ch_ksp), `PCASM`, `PCASMType`, `PCASMWeightedSetScaling()`, `PCASMSetType()`, `PCASMGetLocalSubdomains()`
1297: @*/
1298: PetscErrorCode PCASMWeightedGetScaling(PC pc, PetscInt *n, Vec *scaling[])
1299: {
1300: PC_ASM *osm = (PC_ASM *)pc->data;
1301: PetscBool match;
1303: PetscFunctionBegin;
1305: if (n) PetscAssertPointer(n, 2);
1306: if (scaling) PetscAssertPointer(scaling, 3);
1307: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCASM, &match));
1308: PetscCheck(match, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC is not a PCASM");
1309: if (n) *n = osm->scaling ? osm->n_local_true : 0;
1310: if (scaling) *scaling = osm->scaling;
1311: PetscFunctionReturn(PETSC_SUCCESS);
1312: }
1314: /*@
1315: PCASMGetType - Gets the type of restriction and interpolation used
1316: for local problems in the additive Schwarz method, `PCASM`.
1318: Logically Collective
1320: Input Parameter:
1321: . pc - the preconditioner context
1323: Output Parameter:
1324: . type - variant of `PCASM`, one of
1325: .vb
1326: PC_ASM_NONE - local processor restriction and interpolation
1327: PC_ASM_RESTRICT - full restriction, local processor interpolation
1328: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1329: PC_ASM_BASIC - full interpolation and restriction
1330: PC_ASM_WEIGHTED - full restriction and interpolation with user-provided diagonal scaling
1331: .ve
1333: Options Database Key:
1334: . -pc_asm_type (none|restrict|interpolate|basic|weighted) - Sets `PCASM` type
1336: Level: intermediate
1338: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMGetSubKSP()`, `PCGASM`,
1339: `PCASMCreateSubdomains2D()`, `PCASMType`, `PCASMSetType()`, `PCASMSetLocalType()`, `PCASMGetLocalType()`
1340: @*/
1341: PetscErrorCode PCASMGetType(PC pc, PCASMType *type)
1342: {
1343: PetscFunctionBegin;
1345: PetscUseMethod(pc, "PCASMGetType_C", (PC, PCASMType *), (pc, type));
1346: PetscFunctionReturn(PETSC_SUCCESS);
1347: }
1349: /*@
1350: PCASMSetLocalType - Sets the type of composition used for local problems in the additive Schwarz method, `PCASM`.
1352: Logically Collective
1354: Input Parameters:
1355: + pc - the preconditioner context
1356: - type - type of composition, one of
1357: .vb
1358: PC_COMPOSITE_ADDITIVE - local additive combination
1359: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1360: .ve
1362: Options Database Key:
1363: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1365: Level: intermediate
1367: .seealso: [](ch_ksp), `PCASM`, `PCASMSetType()`, `PCASMGetType()`, `PCASMGetLocalType()`, `PCASMType`, `PCCompositeType`
1368: @*/
1369: PetscErrorCode PCASMSetLocalType(PC pc, PCCompositeType type)
1370: {
1371: PetscFunctionBegin;
1374: PetscTryMethod(pc, "PCASMSetLocalType_C", (PC, PCCompositeType), (pc, type));
1375: PetscFunctionReturn(PETSC_SUCCESS);
1376: }
1378: /*@
1379: PCASMGetLocalType - Gets the type of composition used for local problems in the additive Schwarz method, `PCASM`.
1381: Logically Collective
1383: Input Parameter:
1384: . pc - the preconditioner context
1386: Output Parameter:
1387: . type - type of composition, one of
1388: .vb
1389: PC_COMPOSITE_ADDITIVE - local additive combination
1390: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1391: .ve
1393: Options Database Key:
1394: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1396: Level: intermediate
1398: .seealso: [](ch_ksp), `PCASM`, `PCASMSetType()`, `PCASMGetType()`, `PCASMSetLocalType()`, `PCASMType`, `PCCompositeType`
1399: @*/
1400: PetscErrorCode PCASMGetLocalType(PC pc, PCCompositeType *type)
1401: {
1402: PetscFunctionBegin;
1404: PetscAssertPointer(type, 2);
1405: PetscUseMethod(pc, "PCASMGetLocalType_C", (PC, PCCompositeType *), (pc, type));
1406: PetscFunctionReturn(PETSC_SUCCESS);
1407: }
1409: /*@
1410: PCASMSetSortIndices - Determines whether subdomain indices are sorted.
1412: Logically Collective
1414: Input Parameters:
1415: + pc - the preconditioner context
1416: - doSort - sort the subdomain indices
1418: Level: intermediate
1420: .seealso: [](ch_ksp), `PCASM`, `PCASMSetLocalSubdomains()`, `PCASMSetTotalSubdomains()`, `PCASMGetSubKSP()`,
1421: `PCASMCreateSubdomains2D()`
1422: @*/
1423: PetscErrorCode PCASMSetSortIndices(PC pc, PetscBool doSort)
1424: {
1425: PetscFunctionBegin;
1428: PetscTryMethod(pc, "PCASMSetSortIndices_C", (PC, PetscBool), (pc, doSort));
1429: PetscFunctionReturn(PETSC_SUCCESS);
1430: }
1432: /*@
1433: PCASMGetSubKSP - Gets the local `KSP` contexts for all blocks on
1434: this processor.
1436: Collective iff first_local is requested
1438: Input Parameter:
1439: . pc - the preconditioner context
1441: Output Parameters:
1442: + n_local - the number of blocks on this processor or `NULL`
1443: . first_local - the global number of the first block on this processor or `NULL`, all processors must request or all must pass `NULL`
1444: - ksp - the array of `KSP` contexts
1446: Level: advanced
1448: Notes:
1449: After `PCASMGetSubKSP()` the array of `KSP`s is not to be freed.
1451: You must call `KSPSetUp()` before calling `PCASMGetSubKSP()`.
1453: Fortran Note:
1454: Call `PCASMRestoreSubKSP()` when access to the array of `KSP` is no longer needed. Pass `PETSC_NULL_KSP_POINTER` for `ksp` if not needed.
1456: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`,
1457: `PCASMCreateSubdomains2D()`
1458: @*/
1459: PetscErrorCode PCASMGetSubKSP(PC pc, PetscInt *n_local, PetscInt *first_local, KSP *ksp[])
1460: {
1461: PetscFunctionBegin;
1463: PetscUseMethod(pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (pc, n_local, first_local, ksp));
1464: PetscFunctionReturn(PETSC_SUCCESS);
1465: }
1467: /*MC
1468: PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1469: its own `KSP` object, {cite}`dryja1987additive` and {cite}`1sbg`
1471: Options Database Keys:
1472: + -pc_asm_blocks blks - Sets total blocks. Defaults to one block per MPI process.
1473: . -pc_asm_overlap ovl - Sets overlap
1474: . -pc_asm_type (none|restrict|interpolate|basic|weighted) - Sets `PCASMType`, default is restrict. See `PCASMSetType()`
1475: . -pc_asm_dm_subdomains (true|false) - use subdomains defined by the `DM` with `DMCreateDomainDecomposition()`
1476: - -pc_asm_local_type (additive|multiplicative) - Sets `PCCompositeType`, default is additive. See `PCASMSetLocalType()`
1478: Level: beginner
1480: Notes:
1481: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you
1482: will get a different convergence rate due to the default option of `-pc_asm_type restrict`. Use
1483: `-pc_asm_type basic` to get the same convergence behavior
1485: Each processor can have one or more blocks, but a block cannot be shared by more
1486: than one processor. Use `PCGASM` for subdomains shared by multiple processes.
1488: To set options on the solvers for each block append `-sub_` to all the `KSP`, and `PC`
1489: options database keys. For example, `-sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly`
1491: To set the options on the solvers separate for each block call `PCASMGetSubKSP()`
1492: and set the options directly on the resulting `KSP` object (you can access its `PC` with `KSPGetPC()`)
1494: If the `PC` has an associated `DM`, then, by default, `DMCreateDomainDecomposition()` is used to create the subdomains
1496: Use `PCASMWeightedSetScaling()` with `PC_ASM_WEIGHTED` to supply a diagonal partition of unity on the overlapping subdomains.
1498: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCASMType`, `PCCompositeType`,
1499: `PCBJACOBI`, `PCASMGetSubKSP()`, `PCASMSetLocalSubdomains()`, `PCASMGetType()`, `PCASMSetLocalType()`, `PCASMGetLocalType()`,
1500: `PCASMSetTotalSubdomains()`, `PCSetModifySubMatrices()`, `PCASMSetOverlap()`, `PCASMSetType()`, `PCASMWeightedSetScaling()`
1501: M*/
1503: PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1504: {
1505: PC_ASM *osm;
1507: PetscFunctionBegin;
1508: PetscCall(PetscNew(&osm));
1510: osm->n = PETSC_DECIDE;
1511: osm->n_local = 0;
1512: osm->n_local_true = PETSC_DECIDE;
1513: osm->overlap = 1;
1514: osm->ksp = NULL;
1515: osm->restriction = NULL;
1516: osm->lprolongation = NULL;
1517: osm->lrestriction = NULL;
1518: osm->x = NULL;
1519: osm->y = NULL;
1520: osm->scaling = NULL;
1521: osm->is = NULL;
1522: osm->is_local = NULL;
1523: osm->mat = NULL;
1524: osm->pmat = NULL;
1525: osm->type = PC_ASM_RESTRICT;
1526: osm->loctype = PC_COMPOSITE_ADDITIVE;
1527: osm->sort_indices = PETSC_TRUE;
1528: osm->dm_subdomains = PETSC_FALSE;
1529: osm->sub_mat_type = NULL;
1531: pc->data = (void *)osm;
1532: pc->ops->apply = PCApply_ASM;
1533: pc->ops->matapply = PCMatApply_ASM;
1534: pc->ops->applytranspose = PCApplyTranspose_ASM;
1535: pc->ops->matapplytranspose = PCMatApplyTranspose_ASM;
1536: pc->ops->setup = PCSetUp_ASM;
1537: pc->ops->reset = PCReset_ASM;
1538: pc->ops->destroy = PCDestroy_ASM;
1539: pc->ops->setfromoptions = PCSetFromOptions_ASM;
1540: pc->ops->setuponblocks = PCSetUpOnBlocks_ASM;
1541: pc->ops->view = PCView_ASM;
1542: pc->ops->applyrichardson = NULL;
1544: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetLocalSubdomains_C", PCASMSetLocalSubdomains_ASM));
1545: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetTotalSubdomains_C", PCASMSetTotalSubdomains_ASM));
1546: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetOverlap_C", PCASMSetOverlap_ASM));
1547: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetType_C", PCASMSetType_ASM));
1548: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetType_C", PCASMGetType_ASM));
1549: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMWeightedSetScaling_C", PCASMWeightedSetScaling_ASM));
1550: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMWeightedSetComputeScaling_C", PCASMWeightedSetComputeScaling_ASM));
1551: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetLocalType_C", PCASMSetLocalType_ASM));
1552: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetLocalType_C", PCASMGetLocalType_ASM));
1553: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetSortIndices_C", PCASMSetSortIndices_ASM));
1554: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetSubKSP_C", PCASMGetSubKSP_ASM));
1555: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetSubMatType_C", PCASMGetSubMatType_ASM));
1556: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetSubMatType_C", PCASMSetSubMatType_ASM));
1557: PetscFunctionReturn(PETSC_SUCCESS);
1558: }
1560: /*@
1561: PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1562: preconditioner, `PCASM`, for any problem on a general grid.
1564: Collective
1566: Input Parameters:
1567: + A - The global matrix operator
1568: - n - the number of local blocks
1570: Output Parameter:
1571: . outis - the array of index sets defining the subdomains
1573: Level: advanced
1575: Note:
1576: This generates nonoverlapping subdomains; the `PCASM` will generate the overlap
1577: from these if you use `PCASMSetLocalSubdomains()`
1579: Fortran Note:
1580: `outis` cannot be `PETSC_NULL_IS_POINTER`. Destroy the returned array with `PCASMDestroySubdomains()`,
1581: passing `PETSC_NULL_IS_POINTER` for `is_local` because no local index sets are created.
1583: .seealso: [](ch_ksp), `PCASM`, `PCASMSetLocalSubdomains()`, `PCASMDestroySubdomains()`
1584: @*/
1585: PetscErrorCode PCASMCreateSubdomains(Mat A, PetscInt n, IS *outis[])
1586: {
1587: MatPartitioning mpart;
1588: const char *prefix;
1589: PetscInt i, j, rstart, rend, bs;
1590: PetscBool hasop, isbaij = PETSC_FALSE, foundpart = PETSC_FALSE;
1591: Mat Ad = NULL, adj;
1592: IS ispart, isnumb, *is;
1594: PetscFunctionBegin;
1596: PetscAssertPointer(outis, 3);
1597: PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "number of local blocks must be > 0, n = %" PetscInt_FMT, n);
1599: /* Get prefix, row distribution, and block size */
1600: PetscCall(MatGetOptionsPrefix(A, &prefix));
1601: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
1602: PetscCall(MatGetBlockSize(A, &bs));
1603: 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);
1605: /* Get diagonal block from matrix if possible */
1606: PetscCall(MatHasOperation(A, MATOP_GET_DIAGONAL_BLOCK, &hasop));
1607: if (hasop) PetscCall(MatGetDiagonalBlock(A, &Ad));
1608: if (Ad) {
1609: PetscCall(PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQBAIJ, &isbaij));
1610: if (!isbaij) PetscCall(PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQSBAIJ, &isbaij));
1611: }
1612: if (Ad && n > 1) {
1613: PetscBool match, done;
1614: /* Try to setup a good matrix partitioning if available */
1615: PetscCall(MatPartitioningCreate(PETSC_COMM_SELF, &mpart));
1616: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mpart, prefix));
1617: PetscCall(MatPartitioningSetFromOptions(mpart));
1618: PetscCall(PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGCURRENT, &match));
1619: if (!match) PetscCall(PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGSQUARE, &match));
1620: if (!match) { /* assume a "good" partitioner is available */
1621: PetscInt na;
1622: const PetscInt *ia, *ja;
1623: PetscCall(MatGetRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done));
1624: if (done) {
1625: /* Build adjacency matrix by hand. Unfortunately a call to
1626: MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1627: remove the block-aij structure and we cannot expect
1628: MatPartitioning to split vertices as we need */
1629: PetscInt i, j, len, nnz, cnt, *iia = NULL, *jja = NULL;
1630: const PetscInt *row;
1631: nnz = 0;
1632: for (i = 0; i < na; i++) { /* count number of nonzeros */
1633: len = ia[i + 1] - ia[i];
1634: row = ja + ia[i];
1635: for (j = 0; j < len; j++) {
1636: if (row[j] == i) { /* don't count diagonal */
1637: len--;
1638: break;
1639: }
1640: }
1641: nnz += len;
1642: }
1643: PetscCall(PetscMalloc1(na + 1, &iia));
1644: PetscCall(PetscMalloc1(nnz, &jja));
1645: nnz = 0;
1646: iia[0] = 0;
1647: for (i = 0; i < na; i++) { /* fill adjacency */
1648: cnt = 0;
1649: len = ia[i + 1] - ia[i];
1650: row = ja + ia[i];
1651: for (j = 0; j < len; j++) {
1652: if (row[j] != i) { /* if not diagonal */
1653: jja[nnz + cnt++] = row[j];
1654: }
1655: }
1656: nnz += cnt;
1657: iia[i + 1] = nnz;
1658: }
1659: /* Partitioning of the adjacency matrix */
1660: PetscCall(MatCreateMPIAdj(PETSC_COMM_SELF, na, na, iia, jja, NULL, &adj));
1661: PetscCall(MatPartitioningSetAdjacency(mpart, adj));
1662: PetscCall(MatPartitioningSetNParts(mpart, n));
1663: PetscCall(MatPartitioningApply(mpart, &ispart));
1664: PetscCall(ISPartitioningToNumbering(ispart, &isnumb));
1665: PetscCall(MatDestroy(&adj));
1666: foundpart = PETSC_TRUE;
1667: }
1668: PetscCall(MatRestoreRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done));
1669: }
1670: PetscCall(MatPartitioningDestroy(&mpart));
1671: }
1673: PetscCall(PetscMalloc1(n, &is));
1674: *outis = is;
1676: if (!foundpart) {
1677: /* Partitioning by contiguous chunks of rows */
1679: PetscInt mbs = (rend - rstart) / bs;
1680: PetscInt start = rstart;
1681: for (i = 0; i < n; i++) {
1682: PetscInt count = (mbs / n + ((mbs % n) > i)) * bs;
1683: PetscCall(ISCreateStride(PETSC_COMM_SELF, count, start, 1, &is[i]));
1684: start += count;
1685: }
1687: } else {
1688: /* Partitioning by adjacency of diagonal block */
1690: const PetscInt *numbering;
1691: PetscInt *count, nidx, *indices, *newidx, start = 0;
1692: /* Get node count in each partition */
1693: PetscCall(PetscMalloc1(n, &count));
1694: PetscCall(ISPartitioningCount(ispart, n, count));
1695: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1696: for (i = 0; i < n; i++) count[i] *= bs;
1697: }
1698: /* Build indices from node numbering */
1699: PetscCall(ISGetLocalSize(isnumb, &nidx));
1700: PetscCall(PetscMalloc1(nidx, &indices));
1701: for (i = 0; i < nidx; i++) indices[i] = i; /* needs to be initialized */
1702: PetscCall(ISGetIndices(isnumb, &numbering));
1703: PetscCall(PetscSortIntWithPermutation(nidx, numbering, indices));
1704: PetscCall(ISRestoreIndices(isnumb, &numbering));
1705: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1706: PetscCall(PetscMalloc1(nidx * bs, &newidx));
1707: for (i = 0; i < nidx; i++) {
1708: for (j = 0; j < bs; j++) newidx[i * bs + j] = indices[i] * bs + j;
1709: }
1710: PetscCall(PetscFree(indices));
1711: nidx *= bs;
1712: indices = newidx;
1713: }
1714: /* Shift to get global indices */
1715: for (i = 0; i < nidx; i++) indices[i] += rstart;
1717: /* Build the index sets for each block */
1718: for (i = 0; i < n; i++) {
1719: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, count[i], &indices[start], PETSC_COPY_VALUES, &is[i]));
1720: PetscCall(ISSort(is[i]));
1721: start += count[i];
1722: }
1724: PetscCall(PetscFree(count));
1725: PetscCall(PetscFree(indices));
1726: PetscCall(ISDestroy(&isnumb));
1727: PetscCall(ISDestroy(&ispart));
1728: }
1729: PetscFunctionReturn(PETSC_SUCCESS);
1730: }
1732: /*@
1733: PCASMDestroySubdomains - Destroys the index sets created with
1734: `PCASMCreateSubdomains()` or `PCASMCreateSubdomains2D()`. Should be called after setting subdomains with `PCASMSetLocalSubdomains()`.
1736: Collective
1738: Input Parameters:
1739: + n - the number of index sets
1740: . is - the array of index sets
1741: - is_local - the array of local index sets, can be `NULL`
1743: Level: advanced
1745: Note:
1746: `PCASMCreateSubdomains2D()` also creates an array of non-overlapping local index sets, which must be passed as `is_local`.
1747: `PCASMCreateSubdomains()` creates no local index sets, so pass `NULL` for `is_local`.
1749: Fortran Note:
1750: `is` cannot be `PETSC_NULL_IS_POINTER`. For arrays from `PCASMCreateSubdomains2D()`, pass the returned `is_local`;
1751: for arrays from `PCASMCreateSubdomains()`, pass `PETSC_NULL_IS_POINTER` for `is_local`.
1753: Developer Note:
1754: The `IS` arguments should be a *[]
1756: .seealso: [](ch_ksp), `PCASM`, `PCASMCreateSubdomains()`, `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`
1757: @*/
1758: PetscErrorCode PCASMDestroySubdomains(PetscInt n, IS *is[], IS *is_local[])
1759: {
1760: PetscInt i;
1762: PetscFunctionBegin;
1763: if (n <= 0) PetscFunctionReturn(PETSC_SUCCESS);
1764: if (*is) {
1765: PetscAssertPointer(*is, 2);
1766: for (i = 0; i < n; i++) PetscCall(ISDestroy(&(*is)[i]));
1767: PetscCall(PetscFree(*is));
1768: }
1769: if (is_local && *is_local) {
1770: PetscAssertPointer(*is_local, 3);
1771: for (i = 0; i < n; i++) PetscCall(ISDestroy(&(*is_local)[i]));
1772: PetscCall(PetscFree(*is_local));
1773: }
1774: PetscFunctionReturn(PETSC_SUCCESS);
1775: }
1777: /*@
1778: PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1779: preconditioner, `PCASM`, for a two-dimensional problem on a regular grid.
1781: Not Collective
1783: Input Parameters:
1784: + m - the number of mesh points in the x direction
1785: . n - the number of mesh points in the y direction
1786: . M - the number of subdomains in the x direction
1787: . N - the number of subdomains in the y direction
1788: . dof - degrees of freedom per node
1789: - overlap - overlap in mesh lines
1791: Output Parameters:
1792: + Nsub - the number of subdomains created
1793: . is - array of index sets defining overlapping (if overlap > 0) subdomains
1794: - is_local - array of index sets defining non-overlapping subdomains
1796: Level: advanced
1798: Note:
1799: Presently `PCAMSCreateSubdomains2d()` is valid only for sequential
1800: preconditioners. More general related routines are
1801: `PCASMSetTotalSubdomains()` and `PCASMSetLocalSubdomains()`.
1803: Fortran Note:
1804: Both `is` and `is_local` are created, so neither can be `PETSC_NULL_IS_POINTER`.
1806: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetLocalSubdomains()`, `PCASMGetSubKSP()`,
1807: `PCASMSetOverlap()`
1808: @*/
1809: PetscErrorCode PCASMCreateSubdomains2D(PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt dof, PetscInt overlap, PetscInt *Nsub, IS *is[], IS *is_local[])
1810: {
1811: PetscInt i, j, height, width, ystart, xstart, yleft, yright, xleft, xright, loc_outer;
1812: PetscInt nidx, *idx, loc, ii, jj, count;
1814: PetscFunctionBegin;
1815: PetscCheck(dof == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "dof must be 1");
1817: *Nsub = N * M;
1818: PetscCall(PetscMalloc1(*Nsub, is));
1819: PetscCall(PetscMalloc1(*Nsub, is_local));
1820: ystart = 0;
1821: loc_outer = 0;
1822: for (i = 0; i < N; i++) {
1823: height = n / N + ((n % N) > i); /* height of subdomain */
1824: PetscCheck(height >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many N subdomains for mesh dimension n");
1825: yleft = ystart - overlap;
1826: if (yleft < 0) yleft = 0;
1827: yright = ystart + height + overlap;
1828: if (yright > n) yright = n;
1829: xstart = 0;
1830: for (j = 0; j < M; j++) {
1831: width = m / M + ((m % M) > j); /* width of subdomain */
1832: PetscCheck(width >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many M subdomains for mesh dimension m");
1833: xleft = xstart - overlap;
1834: if (xleft < 0) xleft = 0;
1835: xright = xstart + width + overlap;
1836: if (xright > m) xright = m;
1837: nidx = (xright - xleft) * (yright - yleft);
1838: PetscCall(PetscMalloc1(nidx, &idx));
1839: loc = 0;
1840: for (ii = yleft; ii < yright; ii++) {
1841: count = m * ii + xleft;
1842: for (jj = xleft; jj < xright; jj++) idx[loc++] = count++;
1843: }
1844: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nidx, idx, PETSC_COPY_VALUES, (*is) + loc_outer));
1845: if (overlap == 0) {
1846: PetscCall(PetscObjectReference((PetscObject)(*is)[loc_outer]));
1848: (*is_local)[loc_outer] = (*is)[loc_outer];
1849: } else {
1850: for (loc = 0, ii = ystart; ii < ystart + height; ii++) {
1851: for (jj = xstart; jj < xstart + width; jj++) idx[loc++] = m * ii + jj;
1852: }
1853: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, loc, idx, PETSC_COPY_VALUES, *is_local + loc_outer));
1854: }
1855: PetscCall(PetscFree(idx));
1856: xstart += width;
1857: loc_outer++;
1858: }
1859: ystart += height;
1860: }
1861: for (i = 0; i < *Nsub; i++) PetscCall(ISSort((*is)[i]));
1862: PetscFunctionReturn(PETSC_SUCCESS);
1863: }
1865: /*@
1866: PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
1867: only) for the additive Schwarz preconditioner, `PCASM`.
1869: Not Collective
1871: Input Parameter:
1872: . pc - the preconditioner context
1874: Output Parameters:
1875: + n - if requested, the number of subdomains for this processor (default value = 1)
1876: . is - if requested, the index sets that define the subdomains for this processor
1877: - is_local - if requested, the index sets that define the local part of the subdomains for this processor (can be `NULL`)
1879: Level: advanced
1881: Note:
1882: The `IS` numbering is in the parallel, global numbering of the vector.
1884: Fortran Note:
1885: Pass `PETSC_NULL_IS_POINTER` for `is` or `is_local` if not needed. A requested array that does not exist is returned disassociated.
1887: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`, `PCASMGetSubKSP()`,
1888: `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`, `PCASMGetLocalSubmatrices()`
1889: @*/
1890: PetscErrorCode PCASMGetLocalSubdomains(PC pc, PetscInt *n, IS *is[], IS *is_local[])
1891: {
1892: PC_ASM *osm = (PC_ASM *)pc->data;
1893: PetscBool match;
1895: PetscFunctionBegin;
1897: if (n) PetscAssertPointer(n, 2);
1898: if (is) PetscAssertPointer(is, 3);
1899: if (is_local) PetscAssertPointer(is_local, 4);
1900: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCASM, &match));
1901: PetscCheck(match, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC is not a PCASM");
1902: if (n) *n = osm->n_local_true;
1903: if (is) *is = osm->is;
1904: if (is_local) *is_local = osm->is_local;
1905: PetscFunctionReturn(PETSC_SUCCESS);
1906: }
1908: /*@
1909: PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1910: only) for the additive Schwarz preconditioner, `PCASM`.
1912: Not Collective
1914: Input Parameter:
1915: . pc - the preconditioner context
1917: Output Parameters:
1918: + n - if requested, the number of matrices for this processor (default value = 1)
1919: - mat - if requested, the matrices
1921: Level: advanced
1923: Notes:
1924: Call after `PCSetUp()` (or `KSPSetUp()`) but before `PCApply()` and before `PCSetUpOnBlocks()`)
1926: Usually one would use `PCSetModifySubMatrices()` to change the submatrices in building the preconditioner.
1928: Fortran Note:
1929: Pass `PETSC_NULL_MAT_POINTER` for `mat` if not needed. If the `PC` is not a `PCASM`, `mat` is returned disassociated.
1931: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`, `PCASMGetSubKSP()`,
1932: `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`, `PCASMGetLocalSubdomains()`, `PCSetModifySubMatrices()`
1933: @*/
1934: PetscErrorCode PCASMGetLocalSubmatrices(PC pc, PetscInt *n, Mat *mat[])
1935: {
1936: PC_ASM *osm;
1937: PetscBool match;
1939: PetscFunctionBegin;
1941: if (n) PetscAssertPointer(n, 2);
1942: if (mat) PetscAssertPointer(mat, 3);
1943: PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call after KSPSetUp() or PCSetUp().");
1944: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCASM, &match));
1945: if (!match) {
1946: if (n) *n = 0;
1947: if (mat) *mat = NULL;
1948: } else {
1949: osm = (PC_ASM *)pc->data;
1950: if (n) *n = osm->n_local_true;
1951: if (mat) *mat = osm->pmat;
1952: }
1953: PetscFunctionReturn(PETSC_SUCCESS);
1954: }
1956: /*@
1957: PCASMSetDMSubdomains - Indicates whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible.
1959: Logically Collective
1961: Input Parameters:
1962: + pc - the preconditioner
1963: - flg - boolean indicating whether to use subdomains defined by the `DM`
1965: Options Database Key:
1966: . -pc_asm_dm_subdomains (true|false) - use subdomains defined by the `DM` with `DMCreateDomainDecomposition()`
1968: Level: intermediate
1970: Note:
1971: `PCASMSetTotalSubdomains()` and `PCASMSetOverlap()` take precedence over `PCASMSetDMSubdomains()`,
1972: so setting either of the first two effectively turns the latter off.
1974: Developer Note:
1975: This should be `PCASMSetUseDMSubdomains()`, similarly for the options database key
1977: .seealso: [](ch_ksp), `PCASM`, `PCASMGetDMSubdomains()`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`,
1978: `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`, `PCASMGetLocalSubdomains()`
1979: @*/
1980: PetscErrorCode PCASMSetDMSubdomains(PC pc, PetscBool flg)
1981: {
1982: PC_ASM *osm = (PC_ASM *)pc->data;
1983: PetscBool match;
1985: PetscFunctionBegin;
1988: PetscCheck(!pc->setupcalled, ((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONGSTATE, "Not for a setup PC.");
1989: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCASM, &match));
1990: if (match) osm->dm_subdomains = flg;
1991: PetscFunctionReturn(PETSC_SUCCESS);
1992: }
1994: /*@
1995: PCASMGetDMSubdomains - Returns flag indicating whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible.
1997: Not Collective
1999: Input Parameter:
2000: . pc - the preconditioner
2002: Output Parameter:
2003: . flg - boolean indicating whether to use subdomains defined by the `DM`
2005: Level: intermediate
2007: Developer Note:
2008: This should be `PCASMSetUseDMSubdomains()`
2010: .seealso: [](ch_ksp), `PCASM`, `PCASMSetDMSubdomains()`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`,
2011: `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`, `PCASMGetLocalSubdomains()`
2012: @*/
2013: PetscErrorCode PCASMGetDMSubdomains(PC pc, PetscBool *flg)
2014: {
2015: PC_ASM *osm = (PC_ASM *)pc->data;
2016: PetscBool match;
2018: PetscFunctionBegin;
2020: PetscAssertPointer(flg, 2);
2021: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCASM, &match));
2022: if (match) *flg = osm->dm_subdomains;
2023: else *flg = PETSC_FALSE;
2024: PetscFunctionReturn(PETSC_SUCCESS);
2025: }
2027: /*@
2028: PCASMGetSubMatType - Gets the matrix type used for `PCASM` subsolves, as a string.
2030: Not Collective
2032: Input Parameter:
2033: . pc - the `PC`
2035: Output Parameter:
2036: . sub_mat_type - name of matrix type
2038: Level: advanced
2040: .seealso: [](ch_ksp), `PCASM`, `PCASMSetSubMatType()`, `PCSetType()`, `VecSetType()`, `MatType`, `Mat`
2041: @*/
2042: PetscErrorCode PCASMGetSubMatType(PC pc, MatType *sub_mat_type)
2043: {
2044: PetscFunctionBegin;
2046: PetscTryMethod(pc, "PCASMGetSubMatType_C", (PC, MatType *), (pc, sub_mat_type));
2047: PetscFunctionReturn(PETSC_SUCCESS);
2048: }
2050: /*@
2051: PCASMSetSubMatType - Set the type of matrix used for `PCASM` subsolves
2053: Collective
2055: Input Parameters:
2056: + pc - the `PC` object
2057: - sub_mat_type - the `MatType`
2059: Options Database Key:
2060: . -pc_asm_sub_mat_type sub_mat_type - Sets the matrix type used for subsolves, for example, seqaijviennacl.
2061: If you specify a base name like aijviennacl, the corresponding sequential type is assumed.
2063: Note:
2064: See `MatType` for available types
2066: Level: advanced
2068: .seealso: [](ch_ksp), `PCASM`, `PCASMGetSubMatType()`, `PCSetType()`, `VecSetType()`, `MatType`, `Mat`
2069: @*/
2070: PetscErrorCode PCASMSetSubMatType(PC pc, MatType sub_mat_type)
2071: {
2072: PetscFunctionBegin;
2074: PetscTryMethod(pc, "PCASMSetSubMatType_C", (PC, MatType), (pc, sub_mat_type));
2075: PetscFunctionReturn(PETSC_SUCCESS);
2076: }