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, d;
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 (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: } inwork, outwork;
203: PetscMPIInt size;
205: inwork.max = osm->n_local_true;
206: inwork.sum = osm->n_local_true;
207: PetscCallMPI(MPIU_Allreduce(&inwork, &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, PETSC_FALSE));
263: PetscCall(DMDestroy(&domain_dm[i]));
264: }
265: osm->ksp[i] = ksp;
266: }
267: if (domain_dm) 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: /*
294: Extract out the submatrices
295: */
296: PetscCall(MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, osm->is, scall, &osm->pmat));
297: if (scall == MAT_INITIAL_MATRIX) {
298: PetscCall(PetscObjectGetOptionsPrefix((PetscObject)pc->pmat, &pprefix));
299: for (i = 0; i < osm->n_local_true; i++) PetscCall(PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i], pprefix));
300: if (nullsp) PetscCall(MatRestoreNullSpaces(osm->n_local_true, osm->pmat, &nullsp));
301: }
303: /* Convert the types of the submatrices (if needbe) */
304: if (osm->sub_mat_type) {
305: for (i = 0; i < osm->n_local_true; i++) PetscCall(MatConvert(osm->pmat[i], osm->sub_mat_type, MAT_INPLACE_MATRIX, &osm->pmat[i]));
306: }
308: if (!pc->setupcalled) {
309: VecType vtype;
311: /* Create the local work vectors (from the local matrices) and scatter contexts */
312: PetscCall(MatCreateVecs(pc->pmat, &vec, NULL));
314: 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");
315: if (osm->is_local && osm->type != PC_ASM_BASIC && osm->loctype == PC_COMPOSITE_ADDITIVE) PetscCall(PetscMalloc1(osm->n_local_true, &osm->lprolongation));
316: PetscCall(PetscMalloc1(osm->n_local_true, &osm->lrestriction));
317: PetscCall(PetscMalloc1(osm->n_local_true, &osm->x));
318: PetscCall(PetscMalloc1(osm->n_local_true, &osm->y));
320: PetscCall(ISGetLocalSize(osm->lis, &m));
321: PetscCall(ISCreateStride(PETSC_COMM_SELF, m, 0, 1, &isl));
322: PetscCall(MatGetVecType(osm->pmat[0], &vtype));
323: PetscCall(VecCreate(PETSC_COMM_SELF, &osm->lx));
324: PetscCall(VecSetSizes(osm->lx, m, m));
325: PetscCall(VecSetType(osm->lx, vtype));
326: PetscCall(VecDuplicate(osm->lx, &osm->ly));
327: PetscCall(VecScatterCreate(vec, osm->lis, osm->lx, isl, &osm->restriction));
328: PetscCall(ISDestroy(&isl));
330: for (i = 0; i < osm->n_local_true; ++i) {
331: ISLocalToGlobalMapping ltog;
332: IS isll;
333: const PetscInt *idx_is;
334: PetscInt *idx_lis, nout;
336: PetscCall(ISGetLocalSize(osm->is[i], &m));
337: PetscCall(MatCreateVecs(osm->pmat[i], &osm->x[i], NULL));
338: PetscCall(VecDuplicate(osm->x[i], &osm->y[i]));
340: /* generate a scatter from ly to y[i] picking all the overlapping is[i] entries */
341: PetscCall(ISLocalToGlobalMappingCreateIS(osm->lis, <og));
342: PetscCall(ISGetLocalSize(osm->is[i], &m));
343: PetscCall(ISGetIndices(osm->is[i], &idx_is));
344: PetscCall(PetscMalloc1(m, &idx_lis));
345: PetscCall(ISGlobalToLocalMappingApply(ltog, IS_GTOLM_DROP, m, idx_is, &nout, idx_lis));
346: PetscCheck(nout == m, PETSC_COMM_SELF, PETSC_ERR_PLIB, "is not a subset of lis");
347: PetscCall(ISRestoreIndices(osm->is[i], &idx_is));
348: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, m, idx_lis, PETSC_OWN_POINTER, &isll));
349: PetscCall(ISLocalToGlobalMappingDestroy(<og));
350: PetscCall(ISCreateStride(PETSC_COMM_SELF, m, 0, 1, &isl));
351: PetscCall(VecScatterCreate(osm->ly, isll, osm->y[i], isl, &osm->lrestriction[i]));
352: PetscCall(ISDestroy(&isll));
353: PetscCall(ISDestroy(&isl));
354: if (osm->lprolongation) { /* generate a scatter from y[i] to ly picking only the non-overlapping is_local[i] entries */
355: ISLocalToGlobalMapping ltog;
356: IS isll, isll_local;
357: const PetscInt *idx_local;
358: PetscInt *idx1, *idx2, nout;
360: PetscCall(ISGetLocalSize(osm->is_local[i], &m_local));
361: PetscCall(ISGetIndices(osm->is_local[i], &idx_local));
363: PetscCall(ISLocalToGlobalMappingCreateIS(osm->is[i], <og));
364: PetscCall(PetscMalloc1(m_local, &idx1));
365: PetscCall(ISGlobalToLocalMappingApply(ltog, IS_GTOLM_DROP, m_local, idx_local, &nout, idx1));
366: PetscCall(ISLocalToGlobalMappingDestroy(<og));
367: PetscCheck(nout == m_local, PETSC_COMM_SELF, PETSC_ERR_PLIB, "is_local not a subset of is");
368: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, m_local, idx1, PETSC_OWN_POINTER, &isll));
370: PetscCall(ISLocalToGlobalMappingCreateIS(osm->lis, <og));
371: PetscCall(PetscMalloc1(m_local, &idx2));
372: PetscCall(ISGlobalToLocalMappingApply(ltog, IS_GTOLM_DROP, m_local, idx_local, &nout, idx2));
373: PetscCall(ISLocalToGlobalMappingDestroy(<og));
374: PetscCheck(nout == m_local, PETSC_COMM_SELF, PETSC_ERR_PLIB, "is_local not a subset of lis");
375: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, m_local, idx2, PETSC_OWN_POINTER, &isll_local));
377: PetscCall(ISRestoreIndices(osm->is_local[i], &idx_local));
378: PetscCall(VecScatterCreate(osm->y[i], isll, osm->ly, isll_local, &osm->lprolongation[i]));
380: PetscCall(ISDestroy(&isll));
381: PetscCall(ISDestroy(&isll_local));
382: }
383: }
384: PetscCall(VecDestroy(&vec));
385: }
387: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
388: IS *cis;
389: PetscInt c;
391: PetscCall(PetscMalloc1(osm->n_local_true, &cis));
392: for (c = 0; c < osm->n_local_true; ++c) cis[c] = osm->lis;
393: PetscCall(MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, cis, scall, &osm->lmats));
394: PetscCall(PetscFree(cis));
395: }
397: /* Return control to the user so that the submatrices can be modified (e.g., to apply
398: different boundary conditions for the submatrices than for the global problem) */
399: PetscCall(PCModifySubMatrices(pc, osm->n_local_true, osm->is, osm->is, osm->pmat, pc->modifysubmatricesP));
401: /*
402: Loop over subdomains putting them into local ksp
403: */
404: PetscCall(KSPGetOptionsPrefix(osm->ksp[0], &prefix));
405: for (i = 0; i < osm->n_local_true; i++) {
406: PetscCall(KSPSetOperators(osm->ksp[i], osm->pmat[i], osm->pmat[i]));
407: PetscCall(MatSetOptionsPrefix(osm->pmat[i], prefix));
408: if (!pc->setupcalled) PetscCall(KSPSetFromOptions(osm->ksp[i]));
409: }
410: PetscFunctionReturn(PETSC_SUCCESS);
411: }
413: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
414: {
415: PC_ASM *osm = (PC_ASM *)pc->data;
416: PetscInt i;
417: KSPConvergedReason reason;
419: PetscFunctionBegin;
420: for (i = 0; i < osm->n_local_true; i++) {
421: PetscCall(KSPSetUp(osm->ksp[i]));
422: PetscCall(KSPGetConvergedReason(osm->ksp[i], &reason));
423: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
424: }
425: PetscFunctionReturn(PETSC_SUCCESS);
426: }
428: static PetscErrorCode PCApply_ASM(PC pc, Vec x, Vec y)
429: {
430: PC_ASM *osm = (PC_ASM *)pc->data;
431: PetscInt i, n_local_true = osm->n_local_true;
432: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
434: PetscFunctionBegin;
435: /*
436: support for limiting the restriction or interpolation to only local
437: subdomain values (leaving the other values 0).
438: */
439: if (!(osm->type & PC_ASM_RESTRICT)) {
440: forward = SCATTER_FORWARD_LOCAL;
441: /* have to zero the work RHS since scatter may leave some slots empty */
442: PetscCall(VecSet(osm->lx, 0.0));
443: }
444: if (!(osm->type & PC_ASM_INTERPOLATE)) reverse = SCATTER_REVERSE_LOCAL;
446: 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]);
447: /* zero the global and the local solutions */
448: PetscCall(VecSet(y, 0.0));
449: PetscCall(VecSet(osm->ly, 0.0));
451: /* copy the global RHS to local RHS including the ghost nodes */
452: PetscCall(VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward));
453: PetscCall(VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward));
455: /* restrict local RHS to the overlapping 0-block RHS */
456: PetscCall(VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward));
457: PetscCall(VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward));
459: /* do the local solves */
460: for (i = 0; i < n_local_true; ++i) {
461: /* solve the overlapping i-block */
462: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0));
463: PetscCall(KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]));
464: PetscCall(KSPCheckSolve(osm->ksp[i], pc, osm->y[i]));
465: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0));
467: if (osm->lprolongation && osm->type != PC_ASM_INTERPOLATE) { /* interpolate the non-overlapping i-block solution to the local solution (only for restrictive additive) */
468: PetscCall(VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward));
469: PetscCall(VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward));
470: } else { /* interpolate the overlapping i-block solution to the local solution */
471: PetscCall(VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse));
472: PetscCall(VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse));
473: }
475: if (i < n_local_true - 1) {
476: /* restrict local RHS to the overlapping (i+1)-block RHS */
477: PetscCall(VecScatterBegin(osm->lrestriction[i + 1], osm->lx, osm->x[i + 1], INSERT_VALUES, forward));
478: PetscCall(VecScatterEnd(osm->lrestriction[i + 1], osm->lx, osm->x[i + 1], INSERT_VALUES, forward));
480: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
481: /* update the overlapping (i+1)-block RHS using the current local solution */
482: PetscCall(MatMult(osm->lmats[i + 1], osm->ly, osm->y[i + 1]));
483: PetscCall(VecAXPBY(osm->x[i + 1], -1., 1., osm->y[i + 1]));
484: }
485: }
486: }
487: /* add the local solution to the global solution including the ghost nodes */
488: PetscCall(VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse));
489: PetscCall(VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse));
490: PetscFunctionReturn(PETSC_SUCCESS);
491: }
493: static PetscErrorCode PCMatApply_ASM_Private(PC pc, Mat X, Mat Y, PetscBool transpose)
494: {
495: PC_ASM *osm = (PC_ASM *)pc->data;
496: Mat Z, W;
497: Vec x;
498: PetscInt i, m, N;
499: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
501: PetscFunctionBegin;
502: PetscCheck(osm->n_local_true <= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not yet implemented");
503: /*
504: support for limiting the restriction or interpolation to only local
505: subdomain values (leaving the other values 0).
506: */
507: if ((!transpose && !(osm->type & PC_ASM_RESTRICT)) || (transpose && !(osm->type & PC_ASM_INTERPOLATE))) {
508: forward = SCATTER_FORWARD_LOCAL;
509: /* have to zero the work RHS since scatter may leave some slots empty */
510: PetscCall(VecSet(osm->lx, 0.0));
511: }
512: if ((!transpose && !(osm->type & PC_ASM_INTERPOLATE)) || (transpose && !(osm->type & PC_ASM_RESTRICT))) reverse = SCATTER_REVERSE_LOCAL;
513: PetscCall(VecGetLocalSize(osm->x[0], &m));
514: PetscCall(MatGetSize(X, NULL, &N));
515: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, m, N, NULL, &Z));
517: 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]);
518: /* zero the global and the local solutions */
519: PetscCall(MatZeroEntries(Y));
520: PetscCall(VecSet(osm->ly, 0.0));
522: for (i = 0; i < N; ++i) {
523: PetscCall(MatDenseGetColumnVecRead(X, i, &x));
524: /* copy the global RHS to local RHS including the ghost nodes */
525: PetscCall(VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward));
526: PetscCall(VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward));
527: PetscCall(MatDenseRestoreColumnVecRead(X, i, &x));
529: PetscCall(MatDenseGetColumnVecWrite(Z, i, &x));
530: /* restrict local RHS to the overlapping 0-block RHS */
531: PetscCall(VecScatterBegin(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward));
532: PetscCall(VecScatterEnd(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward));
533: PetscCall(MatDenseRestoreColumnVecWrite(Z, i, &x));
534: }
535: PetscCall(MatCreateSeqDense(PETSC_COMM_SELF, m, N, NULL, &W));
536: /* solve the overlapping 0-block */
537: if (!transpose) {
538: PetscCall(PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[0], Z, W, 0));
539: PetscCall(KSPMatSolve(osm->ksp[0], Z, W));
540: PetscCall(PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[0], Z, W, 0));
541: } else {
542: PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, osm->ksp[0], Z, W, 0));
543: PetscCall(KSPMatSolveTranspose(osm->ksp[0], Z, W));
544: PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, osm->ksp[0], Z, W, 0));
545: }
546: PetscCall(KSPCheckSolve(osm->ksp[0], pc, NULL));
547: PetscCall(MatDestroy(&Z));
549: for (i = 0; i < N; ++i) {
550: PetscCall(VecSet(osm->ly, 0.0));
551: PetscCall(MatDenseGetColumnVecRead(W, i, &x));
552: if (osm->lprolongation && ((!transpose && osm->type != PC_ASM_INTERPOLATE) || (transpose && osm->type != PC_ASM_RESTRICT))) { /* interpolate the non-overlapping 0-block solution to the local solution (only for restrictive additive) */
553: PetscCall(VecScatterBegin(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward));
554: PetscCall(VecScatterEnd(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward));
555: } else { /* interpolate the overlapping 0-block solution to the local solution */
556: PetscCall(VecScatterBegin(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse));
557: PetscCall(VecScatterEnd(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse));
558: }
559: PetscCall(MatDenseRestoreColumnVecRead(W, i, &x));
561: PetscCall(MatDenseGetColumnVecWrite(Y, i, &x));
562: /* add the local solution to the global solution including the ghost nodes */
563: PetscCall(VecScatterBegin(osm->restriction, osm->ly, x, ADD_VALUES, reverse));
564: PetscCall(VecScatterEnd(osm->restriction, osm->ly, x, ADD_VALUES, reverse));
565: PetscCall(MatDenseRestoreColumnVecWrite(Y, i, &x));
566: }
567: PetscCall(MatDestroy(&W));
568: PetscFunctionReturn(PETSC_SUCCESS);
569: }
571: static PetscErrorCode PCMatApply_ASM(PC pc, Mat X, Mat Y)
572: {
573: PetscFunctionBegin;
574: PetscCall(PCMatApply_ASM_Private(pc, X, Y, PETSC_FALSE));
575: PetscFunctionReturn(PETSC_SUCCESS);
576: }
578: static PetscErrorCode PCMatApplyTranspose_ASM(PC pc, Mat X, Mat Y)
579: {
580: PetscFunctionBegin;
581: PetscCall(PCMatApply_ASM_Private(pc, X, Y, PETSC_TRUE));
582: PetscFunctionReturn(PETSC_SUCCESS);
583: }
585: static PetscErrorCode PCApplyTranspose_ASM(PC pc, Vec x, Vec y)
586: {
587: PC_ASM *osm = (PC_ASM *)pc->data;
588: PetscInt i, n_local_true = osm->n_local_true;
589: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
591: PetscFunctionBegin;
592: PetscCheck(osm->n_local_true <= 1 || osm->loctype == PC_COMPOSITE_ADDITIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not yet implemented");
593: /*
594: Support for limiting the restriction or interpolation to only local
595: subdomain values (leaving the other values 0).
597: Note: these are reversed from the PCApply_ASM() because we are applying the
598: transpose of the three terms
599: */
601: if (!(osm->type & PC_ASM_INTERPOLATE)) {
602: forward = SCATTER_FORWARD_LOCAL;
603: /* have to zero the work RHS since scatter may leave some slots empty */
604: PetscCall(VecSet(osm->lx, 0.0));
605: }
606: if (!(osm->type & PC_ASM_RESTRICT)) reverse = SCATTER_REVERSE_LOCAL;
608: /* zero the global and the local solutions */
609: PetscCall(VecSet(y, 0.0));
610: PetscCall(VecSet(osm->ly, 0.0));
612: /* Copy the global RHS to local RHS including the ghost nodes */
613: PetscCall(VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward));
614: PetscCall(VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward));
616: /* Restrict local RHS to the overlapping 0-block RHS */
617: PetscCall(VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward));
618: PetscCall(VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward));
620: /* do the local solves */
621: for (i = 0; i < n_local_true; ++i) {
622: /* solve the overlapping i-block */
623: PetscCall(PetscLogEventBegin(PC_ApplyTransposeOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0));
624: PetscCall(KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i]));
625: PetscCall(KSPCheckSolve(osm->ksp[i], pc, osm->y[i]));
626: PetscCall(PetscLogEventEnd(PC_ApplyTransposeOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0));
628: if (osm->lprolongation && osm->type != PC_ASM_RESTRICT) { /* interpolate the non-overlapping i-block solution to the local solution */
629: PetscCall(VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward));
630: PetscCall(VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward));
631: } else { /* interpolate the overlapping i-block solution to the local solution */
632: PetscCall(VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse));
633: PetscCall(VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse));
634: }
636: if (i < n_local_true - 1) {
637: /* Restrict local RHS to the overlapping (i+1)-block RHS */
638: PetscCall(VecScatterBegin(osm->lrestriction[i + 1], osm->lx, osm->x[i + 1], INSERT_VALUES, forward));
639: PetscCall(VecScatterEnd(osm->lrestriction[i + 1], osm->lx, osm->x[i + 1], INSERT_VALUES, forward));
640: }
641: }
642: /* Add the local solution to the global solution including the ghost nodes */
643: PetscCall(VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse));
644: PetscCall(VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse));
645: PetscFunctionReturn(PETSC_SUCCESS);
646: }
648: static PetscErrorCode PCReset_ASM(PC pc)
649: {
650: PC_ASM *osm = (PC_ASM *)pc->data;
651: PetscInt i;
653: PetscFunctionBegin;
654: if (osm->ksp) {
655: for (i = 0; i < osm->n_local_true; i++) PetscCall(KSPReset(osm->ksp[i]));
656: }
657: if (osm->pmat) {
658: if (osm->n_local_true > 0) PetscCall(MatDestroySubMatrices(osm->n_local_true, &osm->pmat));
659: }
660: if (osm->lrestriction) {
661: PetscCall(VecScatterDestroy(&osm->restriction));
662: for (i = 0; i < osm->n_local_true; i++) {
663: PetscCall(VecScatterDestroy(&osm->lrestriction[i]));
664: if (osm->lprolongation) PetscCall(VecScatterDestroy(&osm->lprolongation[i]));
665: PetscCall(VecDestroy(&osm->x[i]));
666: PetscCall(VecDestroy(&osm->y[i]));
667: }
668: PetscCall(PetscFree(osm->lrestriction));
669: if (osm->lprolongation) PetscCall(PetscFree(osm->lprolongation));
670: PetscCall(PetscFree(osm->x));
671: PetscCall(PetscFree(osm->y));
672: }
673: PetscCall(PCASMDestroySubdomains(osm->n_local_true, &osm->is, &osm->is_local));
674: PetscCall(ISDestroy(&osm->lis));
675: PetscCall(VecDestroy(&osm->lx));
676: PetscCall(VecDestroy(&osm->ly));
677: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) PetscCall(MatDestroyMatrices(osm->n_local_true, &osm->lmats));
679: PetscCall(PetscFree(osm->sub_mat_type));
681: osm->is = NULL;
682: osm->is_local = NULL;
683: PetscFunctionReturn(PETSC_SUCCESS);
684: }
686: static PetscErrorCode PCDestroy_ASM(PC pc)
687: {
688: PC_ASM *osm = (PC_ASM *)pc->data;
689: PetscInt i;
691: PetscFunctionBegin;
692: PetscCall(PCReset_ASM(pc));
693: if (osm->ksp) {
694: for (i = 0; i < osm->n_local_true; i++) PetscCall(KSPDestroy(&osm->ksp[i]));
695: PetscCall(PetscFree(osm->ksp));
696: }
697: PetscCall(PetscFree(pc->data));
699: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetLocalSubdomains_C", NULL));
700: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetTotalSubdomains_C", NULL));
701: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetOverlap_C", NULL));
702: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetType_C", NULL));
703: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetType_C", NULL));
704: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetLocalType_C", NULL));
705: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetLocalType_C", NULL));
706: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetSortIndices_C", NULL));
707: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetSubKSP_C", NULL));
708: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetSubMatType_C", NULL));
709: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetSubMatType_C", NULL));
710: PetscFunctionReturn(PETSC_SUCCESS);
711: }
713: static PetscErrorCode PCSetFromOptions_ASM(PC pc, PetscOptionItems PetscOptionsObject)
714: {
715: PC_ASM *osm = (PC_ASM *)pc->data;
716: PetscInt blocks, ovl;
717: PetscBool flg;
718: PCASMType asmtype;
719: PCCompositeType loctype;
720: char sub_mat_type[256];
722: PetscFunctionBegin;
723: PetscOptionsHeadBegin(PetscOptionsObject, "Additive Schwarz options");
724: PetscCall(PetscOptionsBool("-pc_asm_dm_subdomains", "Use DMCreateDomainDecomposition() to define subdomains", "PCASMSetDMSubdomains", osm->dm_subdomains, &osm->dm_subdomains, &flg));
725: PetscCall(PetscOptionsInt("-pc_asm_blocks", "Number of subdomains", "PCASMSetTotalSubdomains", osm->n, &blocks, &flg));
726: if (flg) {
727: PetscCall(PCASMSetTotalSubdomains(pc, blocks, NULL, NULL));
728: osm->dm_subdomains = PETSC_FALSE;
729: }
730: PetscCall(PetscOptionsInt("-pc_asm_local_blocks", "Number of local subdomains", "PCASMSetLocalSubdomains", osm->n_local_true, &blocks, &flg));
731: if (flg) {
732: PetscCall(PCASMSetLocalSubdomains(pc, blocks, NULL, NULL));
733: osm->dm_subdomains = PETSC_FALSE;
734: }
735: PetscCall(PetscOptionsInt("-pc_asm_overlap", "Number of grid points overlap", "PCASMSetOverlap", osm->overlap, &ovl, &flg));
736: if (flg) {
737: PetscCall(PCASMSetOverlap(pc, ovl));
738: osm->dm_subdomains = PETSC_FALSE;
739: }
740: flg = PETSC_FALSE;
741: PetscCall(PetscOptionsEnum("-pc_asm_type", "Type of restriction/extension", "PCASMSetType", PCASMTypes, (PetscEnum)osm->type, (PetscEnum *)&asmtype, &flg));
742: if (flg) PetscCall(PCASMSetType(pc, asmtype));
743: flg = PETSC_FALSE;
744: PetscCall(PetscOptionsEnum("-pc_asm_local_type", "Type of local solver composition", "PCASMSetLocalType", PCCompositeTypes, (PetscEnum)osm->loctype, (PetscEnum *)&loctype, &flg));
745: if (flg) PetscCall(PCASMSetLocalType(pc, loctype));
746: PetscCall(PetscOptionsFList("-pc_asm_sub_mat_type", "Subsolve Matrix Type", "PCASMSetSubMatType", MatList, NULL, sub_mat_type, 256, &flg));
747: if (flg) PetscCall(PCASMSetSubMatType(pc, sub_mat_type));
748: PetscOptionsHeadEnd();
749: PetscFunctionReturn(PETSC_SUCCESS);
750: }
752: static PetscErrorCode PCASMSetLocalSubdomains_ASM(PC pc, PetscInt n, IS is[], IS is_local[])
753: {
754: PC_ASM *osm = (PC_ASM *)pc->data;
755: PetscInt i;
757: PetscFunctionBegin;
758: PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Each process must have 1 or more blocks, n = %" PetscInt_FMT, n);
759: PetscCheck(!pc->setupcalled || (n == osm->n_local_true && !is), PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCASMSetLocalSubdomains() should be called before calling PCSetUp().");
761: if (!pc->setupcalled) {
762: if (is) {
763: for (i = 0; i < n; i++) PetscCall(PetscObjectReference((PetscObject)is[i]));
764: }
765: if (is_local) {
766: for (i = 0; i < n; i++) PetscCall(PetscObjectReference((PetscObject)is_local[i]));
767: }
768: PetscCall(PCASMDestroySubdomains(osm->n_local_true, &osm->is, &osm->is_local));
770: if (osm->ksp && osm->n_local_true != n) {
771: for (i = 0; i < osm->n_local_true; i++) PetscCall(KSPDestroy(&osm->ksp[i]));
772: PetscCall(PetscFree(osm->ksp));
773: }
775: osm->n_local_true = n;
776: osm->is = NULL;
777: osm->is_local = NULL;
778: if (is) {
779: PetscCall(PetscMalloc1(n, &osm->is));
780: for (i = 0; i < n; i++) osm->is[i] = is[i];
781: /* Flag indicating that the user has set overlapping subdomains so PCASM should not increase their size. */
782: osm->overlap = -1;
783: }
784: if (is_local) {
785: PetscCall(PetscMalloc1(n, &osm->is_local));
786: for (i = 0; i < n; i++) osm->is_local[i] = is_local[i];
787: if (!is) {
788: PetscCall(PetscMalloc1(osm->n_local_true, &osm->is));
789: for (i = 0; i < osm->n_local_true; i++) {
790: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
791: PetscCall(ISDuplicate(osm->is_local[i], &osm->is[i]));
792: PetscCall(ISCopy(osm->is_local[i], osm->is[i]));
793: } else {
794: PetscCall(PetscObjectReference((PetscObject)osm->is_local[i]));
795: osm->is[i] = osm->is_local[i];
796: }
797: }
798: }
799: }
800: }
801: PetscFunctionReturn(PETSC_SUCCESS);
802: }
804: static PetscErrorCode PCASMSetTotalSubdomains_ASM(PC pc, PetscInt N, IS *is, IS *is_local)
805: {
806: PC_ASM *osm = (PC_ASM *)pc->data;
807: PetscMPIInt rank, size;
808: PetscInt n;
810: PetscFunctionBegin;
811: PetscCheck(N >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Number of total blocks must be > 0, N = %" PetscInt_FMT, N);
812: PetscCheck(!is && !is_local, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Use PCASMSetLocalSubdomains() to set specific index sets, they cannot be set globally yet.");
814: /*
815: Split the subdomains equally among all processors
816: */
817: PetscCallMPI(MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank));
818: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
819: n = N / size + ((N % size) > rank);
820: 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);
821: PetscCheck(!pc->setupcalled || n == osm->n_local_true, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "PCASMSetTotalSubdomains() should be called before PCSetUp().");
822: if (!pc->setupcalled) {
823: PetscCall(PCASMDestroySubdomains(osm->n_local_true, &osm->is, &osm->is_local));
825: osm->n_local_true = n;
826: osm->is = NULL;
827: osm->is_local = NULL;
828: }
829: PetscFunctionReturn(PETSC_SUCCESS);
830: }
832: static PetscErrorCode PCASMSetOverlap_ASM(PC pc, PetscInt ovl)
833: {
834: PC_ASM *osm = (PC_ASM *)pc->data;
836: PetscFunctionBegin;
837: PetscCheck(ovl >= 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Negative overlap value requested");
838: PetscCheck(!pc->setupcalled || ovl == osm->overlap, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "PCASMSetOverlap() should be called before PCSetUp().");
839: if (!pc->setupcalled) osm->overlap = ovl;
840: PetscFunctionReturn(PETSC_SUCCESS);
841: }
843: static PetscErrorCode PCASMSetType_ASM(PC pc, PCASMType type)
844: {
845: PC_ASM *osm = (PC_ASM *)pc->data;
847: PetscFunctionBegin;
848: osm->type = type;
849: osm->type_set = PETSC_TRUE;
850: PetscFunctionReturn(PETSC_SUCCESS);
851: }
853: static PetscErrorCode PCASMGetType_ASM(PC pc, PCASMType *type)
854: {
855: PC_ASM *osm = (PC_ASM *)pc->data;
857: PetscFunctionBegin;
858: *type = osm->type;
859: PetscFunctionReturn(PETSC_SUCCESS);
860: }
862: static PetscErrorCode PCASMSetLocalType_ASM(PC pc, PCCompositeType type)
863: {
864: PC_ASM *osm = (PC_ASM *)pc->data;
866: PetscFunctionBegin;
867: PetscCheck(type == PC_COMPOSITE_ADDITIVE || type == PC_COMPOSITE_MULTIPLICATIVE, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Only supports additive or multiplicative as the local type");
868: osm->loctype = type;
869: PetscFunctionReturn(PETSC_SUCCESS);
870: }
872: static PetscErrorCode PCASMGetLocalType_ASM(PC pc, PCCompositeType *type)
873: {
874: PC_ASM *osm = (PC_ASM *)pc->data;
876: PetscFunctionBegin;
877: *type = osm->loctype;
878: PetscFunctionReturn(PETSC_SUCCESS);
879: }
881: static PetscErrorCode PCASMSetSortIndices_ASM(PC pc, PetscBool doSort)
882: {
883: PC_ASM *osm = (PC_ASM *)pc->data;
885: PetscFunctionBegin;
886: osm->sort_indices = doSort;
887: PetscFunctionReturn(PETSC_SUCCESS);
888: }
890: static PetscErrorCode PCASMGetSubKSP_ASM(PC pc, PetscInt *n_local, PetscInt *first_local, KSP **ksp)
891: {
892: PC_ASM *osm = (PC_ASM *)pc->data;
894: PetscFunctionBegin;
895: PetscCheck(osm->n_local_true >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Need to call PCSetUp() on PC (or KSPSetUp() on the outer KSP object) before calling here");
897: if (n_local) *n_local = osm->n_local_true;
898: if (first_local) {
899: PetscCallMPI(MPI_Scan(&osm->n_local_true, first_local, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)pc)));
900: *first_local -= osm->n_local_true;
901: }
902: if (ksp) *ksp = osm->ksp;
903: PetscFunctionReturn(PETSC_SUCCESS);
904: }
906: static PetscErrorCode PCASMGetSubMatType_ASM(PC pc, MatType *sub_mat_type)
907: {
908: PC_ASM *osm = (PC_ASM *)pc->data;
910: PetscFunctionBegin;
912: PetscAssertPointer(sub_mat_type, 2);
913: *sub_mat_type = osm->sub_mat_type;
914: PetscFunctionReturn(PETSC_SUCCESS);
915: }
917: static PetscErrorCode PCASMSetSubMatType_ASM(PC pc, MatType sub_mat_type)
918: {
919: PC_ASM *osm = (PC_ASM *)pc->data;
921: PetscFunctionBegin;
923: PetscCall(PetscFree(osm->sub_mat_type));
924: PetscCall(PetscStrallocpy(sub_mat_type, (char **)&osm->sub_mat_type));
925: PetscFunctionReturn(PETSC_SUCCESS);
926: }
928: /*@
929: PCASMSetLocalSubdomains - Sets the local subdomains (for this processor only) for the additive Schwarz preconditioner `PCASM`.
931: Collective
933: Input Parameters:
934: + pc - the preconditioner context
935: . n - the number of subdomains for this processor (default value = 1)
936: . is - the index set that defines the subdomains for this processor (or `NULL` for PETSc to determine subdomains)
937: the values of the `is` array are copied so you can free the array (not the `IS` in the array) after this call
938: - is_local - the index sets that define the local part of the subdomains for this processor, not used unless `PCASMType` is `PC_ASM_RESTRICT`
939: (or `NULL` to not provide these). The values of the `is_local` array are copied so you can free the array
940: (not the `IS` in the array) after this call
942: Options Database Key:
943: . -pc_asm_local_blocks <blks> - Sets number of local blocks
945: Level: advanced
947: Notes:
948: The `IS` numbering is in the parallel, global numbering of the vector for both `is` and `is_local`
950: By default the `PCASM` preconditioner uses 1 block per processor.
952: Use `PCASMSetTotalSubdomains()` to set the subdomains for all processors.
954: If `is_local` is provided and `PCASMType` is `PC_ASM_RESTRICT` then the solution only over the `is_local` region is interpolated
955: back to form the global solution (this is the standard restricted additive Schwarz method, RASM)
957: If `is_local` is provided and `PCASMType` is `PC_ASM_INTERPOLATE` or `PC_ASM_NONE` then an error is generated since there is
958: no code to handle that case.
960: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`, `PCASMGetSubKSP()`,
961: `PCASMCreateSubdomains2D()`, `PCASMGetLocalSubdomains()`, `PCASMType`, `PCASMSetType()`, `PCGASM`
962: @*/
963: PetscErrorCode PCASMSetLocalSubdomains(PC pc, PetscInt n, IS is[], IS is_local[])
964: {
965: PetscFunctionBegin;
967: PetscTryMethod(pc, "PCASMSetLocalSubdomains_C", (PC, PetscInt, IS[], IS[]), (pc, n, is, is_local));
968: PetscFunctionReturn(PETSC_SUCCESS);
969: }
971: /*@
972: PCASMSetTotalSubdomains - Sets the subdomains for all processors for the
973: additive Schwarz preconditioner, `PCASM`.
975: Collective, all MPI ranks must pass in the same array of `IS`
977: Input Parameters:
978: + pc - the preconditioner context
979: . N - the number of subdomains for all processors
980: . is - the index sets that define the subdomains for all processors (or `NULL` to ask PETSc to determine the subdomains)
981: the values of the `is` array are copied so you can free the array (not the `IS` in the array) after this call
982: - is_local - the index sets that define the local part of the subdomains for this processor (or `NULL` to not provide this information)
983: The values of the `is_local` array are copied so you can free the array (not the `IS` in the array) after this call
985: Options Database Key:
986: . -pc_asm_blocks <blks> - Sets total blocks
988: Level: advanced
990: Notes:
991: Currently you cannot use this to set the actual subdomains with the argument `is` or `is_local`.
993: By default the `PCASM` preconditioner uses 1 block per processor.
995: These index sets cannot be destroyed until after completion of the
996: linear solves for which the `PCASM` preconditioner is being used.
998: Use `PCASMSetLocalSubdomains()` to set local subdomains.
1000: The `IS` numbering is in the parallel, global numbering of the vector for both is and is_local
1002: .seealso: [](ch_ksp), `PCASM`, `PCASMSetLocalSubdomains()`, `PCASMSetOverlap()`, `PCASMGetSubKSP()`,
1003: `PCASMCreateSubdomains2D()`, `PCGASM`
1004: @*/
1005: PetscErrorCode PCASMSetTotalSubdomains(PC pc, PetscInt N, IS is[], IS is_local[])
1006: {
1007: PetscFunctionBegin;
1009: PetscTryMethod(pc, "PCASMSetTotalSubdomains_C", (PC, PetscInt, IS[], IS[]), (pc, N, is, is_local));
1010: PetscFunctionReturn(PETSC_SUCCESS);
1011: }
1013: /*@
1014: PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
1015: additive Schwarz preconditioner, `PCASM`.
1017: Logically Collective
1019: Input Parameters:
1020: + pc - the preconditioner context
1021: - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1)
1023: Options Database Key:
1024: . -pc_asm_overlap <ovl> - Sets overlap
1026: Level: intermediate
1028: Notes:
1029: By default the `PCASM` preconditioner uses 1 block per processor. To use
1030: multiple blocks per perocessor, see `PCASMSetTotalSubdomains()` and
1031: `PCASMSetLocalSubdomains()` (and the option -pc_asm_blocks <blks>).
1033: The overlap defaults to 1, so if one desires that no additional
1034: overlap be computed beyond what may have been set with a call to
1035: `PCASMSetTotalSubdomains()` or `PCASMSetLocalSubdomains()`, then ovl
1036: must be set to be 0. In particular, if one does not explicitly set
1037: the subdomains an application code, then all overlap would be computed
1038: internally by PETSc, and using an overlap of 0 would result in an `PCASM`
1039: variant that is equivalent to the block Jacobi preconditioner.
1041: The default algorithm used by PETSc to increase overlap is fast, but not scalable,
1042: use the option -mat_increase_overlap_scalable when the problem and number of processes is large.
1044: One can define initial index sets with any overlap via
1045: `PCASMSetLocalSubdomains()`; the routine
1046: `PCASMSetOverlap()` merely allows PETSc to extend that overlap further
1047: if desired.
1049: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetLocalSubdomains()`, `PCASMGetSubKSP()`,
1050: `PCASMCreateSubdomains2D()`, `PCASMGetLocalSubdomains()`, `MatIncreaseOverlap()`, `PCGASM`
1051: @*/
1052: PetscErrorCode PCASMSetOverlap(PC pc, PetscInt ovl)
1053: {
1054: PetscFunctionBegin;
1057: PetscTryMethod(pc, "PCASMSetOverlap_C", (PC, PetscInt), (pc, ovl));
1058: PetscFunctionReturn(PETSC_SUCCESS);
1059: }
1061: /*@
1062: PCASMSetType - Sets the type of restriction and interpolation used
1063: for local problems in the additive Schwarz method, `PCASM`.
1065: Logically Collective
1067: Input Parameters:
1068: + pc - the preconditioner context
1069: - type - variant of `PCASM`, one of
1070: .vb
1071: PC_ASM_BASIC - full interpolation and restriction
1072: PC_ASM_RESTRICT - full restriction, local processor interpolation (default)
1073: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1074: PC_ASM_NONE - local processor restriction and interpolation
1075: .ve
1077: Options Database Key:
1078: . -pc_asm_type [basic,restrict,interpolate,none] - Sets `PCASMType`
1080: Level: intermediate
1082: Note:
1083: if the is_local arguments are passed to `PCASMSetLocalSubdomains()` then they are used when `PC_ASM_RESTRICT` has been selected
1084: to limit the local processor interpolation
1086: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMGetSubKSP()`,
1087: `PCASMCreateSubdomains2D()`, `PCASMType`, `PCASMSetLocalType()`, `PCASMGetLocalType()`, `PCGASM`
1088: @*/
1089: PetscErrorCode PCASMSetType(PC pc, PCASMType type)
1090: {
1091: PetscFunctionBegin;
1094: PetscTryMethod(pc, "PCASMSetType_C", (PC, PCASMType), (pc, type));
1095: PetscFunctionReturn(PETSC_SUCCESS);
1096: }
1098: /*@
1099: PCASMGetType - Gets the type of restriction and interpolation used
1100: for local problems in the additive Schwarz method, `PCASM`.
1102: Logically Collective
1104: Input Parameter:
1105: . pc - the preconditioner context
1107: Output Parameter:
1108: . type - variant of `PCASM`, one of
1109: .vb
1110: PC_ASM_BASIC - full interpolation and restriction
1111: PC_ASM_RESTRICT - full restriction, local processor interpolation
1112: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1113: PC_ASM_NONE - local processor restriction and interpolation
1114: .ve
1116: Options Database Key:
1117: . -pc_asm_type [basic,restrict,interpolate,none] - Sets `PCASM` type
1119: Level: intermediate
1121: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMGetSubKSP()`, `PCGASM`,
1122: `PCASMCreateSubdomains2D()`, `PCASMType`, `PCASMSetType()`, `PCASMSetLocalType()`, `PCASMGetLocalType()`
1123: @*/
1124: PetscErrorCode PCASMGetType(PC pc, PCASMType *type)
1125: {
1126: PetscFunctionBegin;
1128: PetscUseMethod(pc, "PCASMGetType_C", (PC, PCASMType *), (pc, type));
1129: PetscFunctionReturn(PETSC_SUCCESS);
1130: }
1132: /*@
1133: PCASMSetLocalType - Sets the type of composition used for local problems in the additive Schwarz method, `PCASM`.
1135: Logically Collective
1137: Input Parameters:
1138: + pc - the preconditioner context
1139: - type - type of composition, one of
1140: .vb
1141: PC_COMPOSITE_ADDITIVE - local additive combination
1142: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1143: .ve
1145: Options Database Key:
1146: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1148: Level: intermediate
1150: .seealso: [](ch_ksp), `PCASM`, `PCASMSetType()`, `PCASMGetType()`, `PCASMGetLocalType()`, `PCASMType`, `PCCompositeType`
1151: @*/
1152: PetscErrorCode PCASMSetLocalType(PC pc, PCCompositeType type)
1153: {
1154: PetscFunctionBegin;
1157: PetscTryMethod(pc, "PCASMSetLocalType_C", (PC, PCCompositeType), (pc, type));
1158: PetscFunctionReturn(PETSC_SUCCESS);
1159: }
1161: /*@
1162: PCASMGetLocalType - Gets the type of composition used for local problems in the additive Schwarz method, `PCASM`.
1164: Logically Collective
1166: Input Parameter:
1167: . pc - the preconditioner context
1169: Output Parameter:
1170: . type - type of composition, one of
1171: .vb
1172: PC_COMPOSITE_ADDITIVE - local additive combination
1173: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1174: .ve
1176: Options Database Key:
1177: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1179: Level: intermediate
1181: .seealso: [](ch_ksp), `PCASM`, `PCASMSetType()`, `PCASMGetType()`, `PCASMSetLocalType()`, `PCASMType`, `PCCompositeType`
1182: @*/
1183: PetscErrorCode PCASMGetLocalType(PC pc, PCCompositeType *type)
1184: {
1185: PetscFunctionBegin;
1187: PetscAssertPointer(type, 2);
1188: PetscUseMethod(pc, "PCASMGetLocalType_C", (PC, PCCompositeType *), (pc, type));
1189: PetscFunctionReturn(PETSC_SUCCESS);
1190: }
1192: /*@
1193: PCASMSetSortIndices - Determines whether subdomain indices are sorted.
1195: Logically Collective
1197: Input Parameters:
1198: + pc - the preconditioner context
1199: - doSort - sort the subdomain indices
1201: Level: intermediate
1203: .seealso: [](ch_ksp), `PCASM`, `PCASMSetLocalSubdomains()`, `PCASMSetTotalSubdomains()`, `PCASMGetSubKSP()`,
1204: `PCASMCreateSubdomains2D()`
1205: @*/
1206: PetscErrorCode PCASMSetSortIndices(PC pc, PetscBool doSort)
1207: {
1208: PetscFunctionBegin;
1211: PetscTryMethod(pc, "PCASMSetSortIndices_C", (PC, PetscBool), (pc, doSort));
1212: PetscFunctionReturn(PETSC_SUCCESS);
1213: }
1215: /*@C
1216: PCASMGetSubKSP - Gets the local `KSP` contexts for all blocks on
1217: this processor.
1219: Collective iff first_local is requested
1221: Input Parameter:
1222: . pc - the preconditioner context
1224: Output Parameters:
1225: + n_local - the number of blocks on this processor or `NULL`
1226: . first_local - the global number of the first block on this processor or `NULL`, all processors must request or all must pass `NULL`
1227: - ksp - the array of `KSP` contexts
1229: Level: advanced
1231: Notes:
1232: After `PCASMGetSubKSP()` the array of `KSP`s is not to be freed.
1234: You must call `KSPSetUp()` before calling `PCASMGetSubKSP()`.
1236: Fortran Note:
1237: Call `PCASMRestoreSubKSP()` when access to the array of `KSP` is no longer needed
1239: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`,
1240: `PCASMCreateSubdomains2D()`,
1241: @*/
1242: PetscErrorCode PCASMGetSubKSP(PC pc, PetscInt *n_local, PetscInt *first_local, KSP *ksp[])
1243: {
1244: PetscFunctionBegin;
1246: PetscUseMethod(pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (pc, n_local, first_local, ksp));
1247: PetscFunctionReturn(PETSC_SUCCESS);
1248: }
1250: /*MC
1251: PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1252: its own `KSP` object, {cite}`dryja1987additive` and {cite}`1sbg`
1254: Options Database Keys:
1255: + -pc_asm_blocks <blks> - Sets total blocks. Defaults to one block per MPI process.
1256: . -pc_asm_overlap <ovl> - Sets overlap
1257: . -pc_asm_type [basic,restrict,interpolate,none] - Sets `PCASMType`, default is restrict. See `PCASMSetType()`
1258: . -pc_asm_dm_subdomains <bool> - use subdomains defined by the `DM` with `DMCreateDomainDecomposition()`
1259: - -pc_asm_local_type [additive, multiplicative] - Sets `PCCompositeType`, default is additive. See `PCASMSetLocalType()`
1261: Level: beginner
1263: Notes:
1264: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you
1265: will get a different convergence rate due to the default option of `-pc_asm_type restrict`. Use
1266: `-pc_asm_type basic` to get the same convergence behavior
1268: Each processor can have one or more blocks, but a block cannot be shared by more
1269: than one processor. Use `PCGASM` for subdomains shared by multiple processes.
1271: To set options on the solvers for each block append `-sub_` to all the `KSP`, and `PC`
1272: options database keys. For example, `-sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly`
1274: To set the options on the solvers separate for each block call `PCASMGetSubKSP()`
1275: and set the options directly on the resulting `KSP` object (you can access its `PC` with `KSPGetPC()`)
1277: If the `PC` has an associated `DM`, then, by default, `DMCreateDomainDecomposition()` is used to create the subdomains
1279: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCASMType`, `PCCompositeType`,
1280: `PCBJACOBI`, `PCASMGetSubKSP()`, `PCASMSetLocalSubdomains()`, `PCASMType`, `PCASMGetType()`, `PCASMSetLocalType()`, `PCASMGetLocalType()`
1281: `PCASMSetTotalSubdomains()`, `PCSetModifySubMatrices()`, `PCASMSetOverlap()`, `PCASMSetType()`, `PCCompositeType`
1282: M*/
1284: PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1285: {
1286: PC_ASM *osm;
1288: PetscFunctionBegin;
1289: PetscCall(PetscNew(&osm));
1291: osm->n = PETSC_DECIDE;
1292: osm->n_local = 0;
1293: osm->n_local_true = PETSC_DECIDE;
1294: osm->overlap = 1;
1295: osm->ksp = NULL;
1296: osm->restriction = NULL;
1297: osm->lprolongation = NULL;
1298: osm->lrestriction = NULL;
1299: osm->x = NULL;
1300: osm->y = NULL;
1301: osm->is = NULL;
1302: osm->is_local = NULL;
1303: osm->mat = NULL;
1304: osm->pmat = NULL;
1305: osm->type = PC_ASM_RESTRICT;
1306: osm->loctype = PC_COMPOSITE_ADDITIVE;
1307: osm->sort_indices = PETSC_TRUE;
1308: osm->dm_subdomains = PETSC_FALSE;
1309: osm->sub_mat_type = NULL;
1311: pc->data = (void *)osm;
1312: pc->ops->apply = PCApply_ASM;
1313: pc->ops->matapply = PCMatApply_ASM;
1314: pc->ops->applytranspose = PCApplyTranspose_ASM;
1315: pc->ops->matapplytranspose = PCMatApplyTranspose_ASM;
1316: pc->ops->setup = PCSetUp_ASM;
1317: pc->ops->reset = PCReset_ASM;
1318: pc->ops->destroy = PCDestroy_ASM;
1319: pc->ops->setfromoptions = PCSetFromOptions_ASM;
1320: pc->ops->setuponblocks = PCSetUpOnBlocks_ASM;
1321: pc->ops->view = PCView_ASM;
1322: pc->ops->applyrichardson = NULL;
1324: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetLocalSubdomains_C", PCASMSetLocalSubdomains_ASM));
1325: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetTotalSubdomains_C", PCASMSetTotalSubdomains_ASM));
1326: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetOverlap_C", PCASMSetOverlap_ASM));
1327: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetType_C", PCASMSetType_ASM));
1328: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetType_C", PCASMGetType_ASM));
1329: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetLocalType_C", PCASMSetLocalType_ASM));
1330: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetLocalType_C", PCASMGetLocalType_ASM));
1331: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetSortIndices_C", PCASMSetSortIndices_ASM));
1332: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetSubKSP_C", PCASMGetSubKSP_ASM));
1333: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMGetSubMatType_C", PCASMGetSubMatType_ASM));
1334: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCASMSetSubMatType_C", PCASMSetSubMatType_ASM));
1335: PetscFunctionReturn(PETSC_SUCCESS);
1336: }
1338: /*@C
1339: PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1340: preconditioner, `PCASM`, for any problem on a general grid.
1342: Collective
1344: Input Parameters:
1345: + A - The global matrix operator
1346: - n - the number of local blocks
1348: Output Parameter:
1349: . outis - the array of index sets defining the subdomains
1351: Level: advanced
1353: Note:
1354: This generates nonoverlapping subdomains; the `PCASM` will generate the overlap
1355: from these if you use `PCASMSetLocalSubdomains()`
1357: .seealso: [](ch_ksp), `PCASM`, `PCASMSetLocalSubdomains()`, `PCASMDestroySubdomains()`
1358: @*/
1359: PetscErrorCode PCASMCreateSubdomains(Mat A, PetscInt n, IS *outis[])
1360: {
1361: MatPartitioning mpart;
1362: const char *prefix;
1363: PetscInt i, j, rstart, rend, bs;
1364: PetscBool hasop, isbaij = PETSC_FALSE, foundpart = PETSC_FALSE;
1365: Mat Ad = NULL, adj;
1366: IS ispart, isnumb, *is;
1368: PetscFunctionBegin;
1370: PetscAssertPointer(outis, 3);
1371: PetscCheck(n >= 1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "number of local blocks must be > 0, n = %" PetscInt_FMT, n);
1373: /* Get prefix, row distribution, and block size */
1374: PetscCall(MatGetOptionsPrefix(A, &prefix));
1375: PetscCall(MatGetOwnershipRange(A, &rstart, &rend));
1376: PetscCall(MatGetBlockSize(A, &bs));
1377: 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);
1379: /* Get diagonal block from matrix if possible */
1380: PetscCall(MatHasOperation(A, MATOP_GET_DIAGONAL_BLOCK, &hasop));
1381: if (hasop) PetscCall(MatGetDiagonalBlock(A, &Ad));
1382: if (Ad) {
1383: PetscCall(PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQBAIJ, &isbaij));
1384: if (!isbaij) PetscCall(PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQSBAIJ, &isbaij));
1385: }
1386: if (Ad && n > 1) {
1387: PetscBool match, done;
1388: /* Try to setup a good matrix partitioning if available */
1389: PetscCall(MatPartitioningCreate(PETSC_COMM_SELF, &mpart));
1390: PetscCall(PetscObjectSetOptionsPrefix((PetscObject)mpart, prefix));
1391: PetscCall(MatPartitioningSetFromOptions(mpart));
1392: PetscCall(PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGCURRENT, &match));
1393: if (!match) PetscCall(PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGSQUARE, &match));
1394: if (!match) { /* assume a "good" partitioner is available */
1395: PetscInt na;
1396: const PetscInt *ia, *ja;
1397: PetscCall(MatGetRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done));
1398: if (done) {
1399: /* Build adjacency matrix by hand. Unfortunately a call to
1400: MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1401: remove the block-aij structure and we cannot expect
1402: MatPartitioning to split vertices as we need */
1403: PetscInt i, j, len, nnz, cnt, *iia = NULL, *jja = NULL;
1404: const PetscInt *row;
1405: nnz = 0;
1406: for (i = 0; i < na; i++) { /* count number of nonzeros */
1407: len = ia[i + 1] - ia[i];
1408: row = ja + ia[i];
1409: for (j = 0; j < len; j++) {
1410: if (row[j] == i) { /* don't count diagonal */
1411: len--;
1412: break;
1413: }
1414: }
1415: nnz += len;
1416: }
1417: PetscCall(PetscMalloc1(na + 1, &iia));
1418: PetscCall(PetscMalloc1(nnz, &jja));
1419: nnz = 0;
1420: iia[0] = 0;
1421: for (i = 0; i < na; i++) { /* fill adjacency */
1422: cnt = 0;
1423: len = ia[i + 1] - ia[i];
1424: row = ja + ia[i];
1425: for (j = 0; j < len; j++) {
1426: if (row[j] != i) { /* if not diagonal */
1427: jja[nnz + cnt++] = row[j];
1428: }
1429: }
1430: nnz += cnt;
1431: iia[i + 1] = nnz;
1432: }
1433: /* Partitioning of the adjacency matrix */
1434: PetscCall(MatCreateMPIAdj(PETSC_COMM_SELF, na, na, iia, jja, NULL, &adj));
1435: PetscCall(MatPartitioningSetAdjacency(mpart, adj));
1436: PetscCall(MatPartitioningSetNParts(mpart, n));
1437: PetscCall(MatPartitioningApply(mpart, &ispart));
1438: PetscCall(ISPartitioningToNumbering(ispart, &isnumb));
1439: PetscCall(MatDestroy(&adj));
1440: foundpart = PETSC_TRUE;
1441: }
1442: PetscCall(MatRestoreRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done));
1443: }
1444: PetscCall(MatPartitioningDestroy(&mpart));
1445: }
1447: PetscCall(PetscMalloc1(n, &is));
1448: *outis = is;
1450: if (!foundpart) {
1451: /* Partitioning by contiguous chunks of rows */
1453: PetscInt mbs = (rend - rstart) / bs;
1454: PetscInt start = rstart;
1455: for (i = 0; i < n; i++) {
1456: PetscInt count = (mbs / n + ((mbs % n) > i)) * bs;
1457: PetscCall(ISCreateStride(PETSC_COMM_SELF, count, start, 1, &is[i]));
1458: start += count;
1459: }
1461: } else {
1462: /* Partitioning by adjacency of diagonal block */
1464: const PetscInt *numbering;
1465: PetscInt *count, nidx, *indices, *newidx, start = 0;
1466: /* Get node count in each partition */
1467: PetscCall(PetscMalloc1(n, &count));
1468: PetscCall(ISPartitioningCount(ispart, n, count));
1469: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1470: for (i = 0; i < n; i++) count[i] *= bs;
1471: }
1472: /* Build indices from node numbering */
1473: PetscCall(ISGetLocalSize(isnumb, &nidx));
1474: PetscCall(PetscMalloc1(nidx, &indices));
1475: for (i = 0; i < nidx; i++) indices[i] = i; /* needs to be initialized */
1476: PetscCall(ISGetIndices(isnumb, &numbering));
1477: PetscCall(PetscSortIntWithPermutation(nidx, numbering, indices));
1478: PetscCall(ISRestoreIndices(isnumb, &numbering));
1479: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1480: PetscCall(PetscMalloc1(nidx * bs, &newidx));
1481: for (i = 0; i < nidx; i++) {
1482: for (j = 0; j < bs; j++) newidx[i * bs + j] = indices[i] * bs + j;
1483: }
1484: PetscCall(PetscFree(indices));
1485: nidx *= bs;
1486: indices = newidx;
1487: }
1488: /* Shift to get global indices */
1489: for (i = 0; i < nidx; i++) indices[i] += rstart;
1491: /* Build the index sets for each block */
1492: for (i = 0; i < n; i++) {
1493: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, count[i], &indices[start], PETSC_COPY_VALUES, &is[i]));
1494: PetscCall(ISSort(is[i]));
1495: start += count[i];
1496: }
1498: PetscCall(PetscFree(count));
1499: PetscCall(PetscFree(indices));
1500: PetscCall(ISDestroy(&isnumb));
1501: PetscCall(ISDestroy(&ispart));
1502: }
1503: PetscFunctionReturn(PETSC_SUCCESS);
1504: }
1506: /*@C
1507: PCASMDestroySubdomains - Destroys the index sets created with
1508: `PCASMCreateSubdomains()`. Should be called after setting subdomains with `PCASMSetLocalSubdomains()`.
1510: Collective
1512: Input Parameters:
1513: + n - the number of index sets
1514: . is - the array of index sets
1515: - is_local - the array of local index sets, can be `NULL`
1517: Level: advanced
1519: Developer Note:
1520: The `IS` arguments should be a *[]
1522: .seealso: [](ch_ksp), `PCASM`, `PCASMCreateSubdomains()`, `PCASMSetLocalSubdomains()`
1523: @*/
1524: PetscErrorCode PCASMDestroySubdomains(PetscInt n, IS *is[], IS *is_local[])
1525: {
1526: PetscInt i;
1528: PetscFunctionBegin;
1529: if (n <= 0) PetscFunctionReturn(PETSC_SUCCESS);
1530: if (*is) {
1531: PetscAssertPointer(*is, 2);
1532: for (i = 0; i < n; i++) PetscCall(ISDestroy(&(*is)[i]));
1533: PetscCall(PetscFree(*is));
1534: }
1535: if (is_local && *is_local) {
1536: PetscAssertPointer(*is_local, 3);
1537: for (i = 0; i < n; i++) PetscCall(ISDestroy(&(*is_local)[i]));
1538: PetscCall(PetscFree(*is_local));
1539: }
1540: PetscFunctionReturn(PETSC_SUCCESS);
1541: }
1543: /*@C
1544: PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1545: preconditioner, `PCASM`, for a two-dimensional problem on a regular grid.
1547: Not Collective
1549: Input Parameters:
1550: + m - the number of mesh points in the x direction
1551: . n - the number of mesh points in the y direction
1552: . M - the number of subdomains in the x direction
1553: . N - the number of subdomains in the y direction
1554: . dof - degrees of freedom per node
1555: - overlap - overlap in mesh lines
1557: Output Parameters:
1558: + Nsub - the number of subdomains created
1559: . is - array of index sets defining overlapping (if overlap > 0) subdomains
1560: - is_local - array of index sets defining non-overlapping subdomains
1562: Level: advanced
1564: Note:
1565: Presently `PCAMSCreateSubdomains2d()` is valid only for sequential
1566: preconditioners. More general related routines are
1567: `PCASMSetTotalSubdomains()` and `PCASMSetLocalSubdomains()`.
1569: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetLocalSubdomains()`, `PCASMGetSubKSP()`,
1570: `PCASMSetOverlap()`
1571: @*/
1572: PetscErrorCode PCASMCreateSubdomains2D(PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt dof, PetscInt overlap, PetscInt *Nsub, IS *is[], IS *is_local[])
1573: {
1574: PetscInt i, j, height, width, ystart, xstart, yleft, yright, xleft, xright, loc_outer;
1575: PetscInt nidx, *idx, loc, ii, jj, count;
1577: PetscFunctionBegin;
1578: PetscCheck(dof == 1, PETSC_COMM_SELF, PETSC_ERR_SUP, "dof must be 1");
1580: *Nsub = N * M;
1581: PetscCall(PetscMalloc1(*Nsub, is));
1582: PetscCall(PetscMalloc1(*Nsub, is_local));
1583: ystart = 0;
1584: loc_outer = 0;
1585: for (i = 0; i < N; i++) {
1586: height = n / N + ((n % N) > i); /* height of subdomain */
1587: PetscCheck(height >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many N subdomains for mesh dimension n");
1588: yleft = ystart - overlap;
1589: if (yleft < 0) yleft = 0;
1590: yright = ystart + height + overlap;
1591: if (yright > n) yright = n;
1592: xstart = 0;
1593: for (j = 0; j < M; j++) {
1594: width = m / M + ((m % M) > j); /* width of subdomain */
1595: PetscCheck(width >= 2, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Too many M subdomains for mesh dimension m");
1596: xleft = xstart - overlap;
1597: if (xleft < 0) xleft = 0;
1598: xright = xstart + width + overlap;
1599: if (xright > m) xright = m;
1600: nidx = (xright - xleft) * (yright - yleft);
1601: PetscCall(PetscMalloc1(nidx, &idx));
1602: loc = 0;
1603: for (ii = yleft; ii < yright; ii++) {
1604: count = m * ii + xleft;
1605: for (jj = xleft; jj < xright; jj++) idx[loc++] = count++;
1606: }
1607: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, nidx, idx, PETSC_COPY_VALUES, (*is) + loc_outer));
1608: if (overlap == 0) {
1609: PetscCall(PetscObjectReference((PetscObject)(*is)[loc_outer]));
1611: (*is_local)[loc_outer] = (*is)[loc_outer];
1612: } else {
1613: for (loc = 0, ii = ystart; ii < ystart + height; ii++) {
1614: for (jj = xstart; jj < xstart + width; jj++) idx[loc++] = m * ii + jj;
1615: }
1616: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, loc, idx, PETSC_COPY_VALUES, *is_local + loc_outer));
1617: }
1618: PetscCall(PetscFree(idx));
1619: xstart += width;
1620: loc_outer++;
1621: }
1622: ystart += height;
1623: }
1624: for (i = 0; i < *Nsub; i++) PetscCall(ISSort((*is)[i]));
1625: PetscFunctionReturn(PETSC_SUCCESS);
1626: }
1628: /*@C
1629: PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
1630: only) for the additive Schwarz preconditioner, `PCASM`.
1632: Not Collective
1634: Input Parameter:
1635: . pc - the preconditioner context
1637: Output Parameters:
1638: + n - if requested, the number of subdomains for this processor (default value = 1)
1639: . is - if requested, the index sets that define the subdomains for this processor
1640: - is_local - if requested, the index sets that define the local part of the subdomains for this processor (can be `NULL`)
1642: Level: advanced
1644: Note:
1645: The `IS` numbering is in the parallel, global numbering of the vector.
1647: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`, `PCASMGetSubKSP()`,
1648: `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`, `PCASMGetLocalSubmatrices()`
1649: @*/
1650: PetscErrorCode PCASMGetLocalSubdomains(PC pc, PetscInt *n, IS *is[], IS *is_local[])
1651: {
1652: PC_ASM *osm = (PC_ASM *)pc->data;
1653: PetscBool match;
1655: PetscFunctionBegin;
1657: if (n) PetscAssertPointer(n, 2);
1658: if (is) PetscAssertPointer(is, 3);
1659: if (is_local) PetscAssertPointer(is_local, 4);
1660: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCASM, &match));
1661: PetscCheck(match, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "PC is not a PCASM");
1662: if (n) *n = osm->n_local_true;
1663: if (is) *is = osm->is;
1664: if (is_local) *is_local = osm->is_local;
1665: PetscFunctionReturn(PETSC_SUCCESS);
1666: }
1668: /*@C
1669: PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1670: only) for the additive Schwarz preconditioner, `PCASM`.
1672: Not Collective
1674: Input Parameter:
1675: . pc - the preconditioner context
1677: Output Parameters:
1678: + n - if requested, the number of matrices for this processor (default value = 1)
1679: - mat - if requested, the matrices
1681: Level: advanced
1683: Notes:
1684: Call after `PCSetUp()` (or `KSPSetUp()`) but before `PCApply()` and before `PCSetUpOnBlocks()`)
1686: Usually one would use `PCSetModifySubMatrices()` to change the submatrices in building the preconditioner.
1688: .seealso: [](ch_ksp), `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`, `PCASMGetSubKSP()`,
1689: `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`, `PCASMGetLocalSubdomains()`, `PCSetModifySubMatrices()`
1690: @*/
1691: PetscErrorCode PCASMGetLocalSubmatrices(PC pc, PetscInt *n, Mat *mat[])
1692: {
1693: PC_ASM *osm;
1694: PetscBool match;
1696: PetscFunctionBegin;
1698: if (n) PetscAssertPointer(n, 2);
1699: if (mat) PetscAssertPointer(mat, 3);
1700: PetscCheck(pc->setupcalled, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must call after KSPSetUp() or PCSetUp().");
1701: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCASM, &match));
1702: if (!match) {
1703: if (n) *n = 0;
1704: if (mat) *mat = NULL;
1705: } else {
1706: osm = (PC_ASM *)pc->data;
1707: if (n) *n = osm->n_local_true;
1708: if (mat) *mat = osm->pmat;
1709: }
1710: PetscFunctionReturn(PETSC_SUCCESS);
1711: }
1713: /*@
1714: PCASMSetDMSubdomains - Indicates whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible.
1716: Logically Collective
1718: Input Parameters:
1719: + pc - the preconditioner
1720: - flg - boolean indicating whether to use subdomains defined by the `DM`
1722: Options Database Key:
1723: . -pc_asm_dm_subdomains <bool> - use subdomains defined by the `DM` with `DMCreateDomainDecomposition()`
1725: Level: intermediate
1727: Note:
1728: `PCASMSetTotalSubdomains()` and `PCASMSetOverlap()` take precedence over `PCASMSetDMSubdomains()`,
1729: so setting either of the first two effectively turns the latter off.
1731: Developer Note:
1732: This should be `PCASMSetUseDMSubdomains()`, similarly for the options database key
1734: .seealso: [](ch_ksp), `PCASM`, `PCASMGetDMSubdomains()`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`
1735: `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`, `PCASMGetLocalSubdomains()`
1736: @*/
1737: PetscErrorCode PCASMSetDMSubdomains(PC pc, PetscBool flg)
1738: {
1739: PC_ASM *osm = (PC_ASM *)pc->data;
1740: PetscBool match;
1742: PetscFunctionBegin;
1745: PetscCheck(!pc->setupcalled, ((PetscObject)pc)->comm, PETSC_ERR_ARG_WRONGSTATE, "Not for a setup PC.");
1746: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCASM, &match));
1747: if (match) osm->dm_subdomains = flg;
1748: PetscFunctionReturn(PETSC_SUCCESS);
1749: }
1751: /*@
1752: PCASMGetDMSubdomains - Returns flag indicating whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible.
1754: Not Collective
1756: Input Parameter:
1757: . pc - the preconditioner
1759: Output Parameter:
1760: . flg - boolean indicating whether to use subdomains defined by the `DM`
1762: Level: intermediate
1764: Developer Note:
1765: This should be `PCASMSetUseDMSubdomains()`
1767: .seealso: [](ch_ksp), `PCASM`, `PCASMSetDMSubdomains()`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`
1768: `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`, `PCASMGetLocalSubdomains()`
1769: @*/
1770: PetscErrorCode PCASMGetDMSubdomains(PC pc, PetscBool *flg)
1771: {
1772: PC_ASM *osm = (PC_ASM *)pc->data;
1773: PetscBool match;
1775: PetscFunctionBegin;
1777: PetscAssertPointer(flg, 2);
1778: PetscCall(PetscObjectTypeCompare((PetscObject)pc, PCASM, &match));
1779: if (match) *flg = osm->dm_subdomains;
1780: else *flg = PETSC_FALSE;
1781: PetscFunctionReturn(PETSC_SUCCESS);
1782: }
1784: /*@
1785: PCASMGetSubMatType - Gets the matrix type used for `PCASM` subsolves, as a string.
1787: Not Collective
1789: Input Parameter:
1790: . pc - the `PC`
1792: Output Parameter:
1793: . sub_mat_type - name of matrix type
1795: Level: advanced
1797: .seealso: [](ch_ksp), `PCASM`, `PCASMSetSubMatType()`, `PCSetType()`, `VecSetType()`, `MatType`, `Mat`
1798: @*/
1799: PetscErrorCode PCASMGetSubMatType(PC pc, MatType *sub_mat_type)
1800: {
1801: PetscFunctionBegin;
1803: PetscTryMethod(pc, "PCASMGetSubMatType_C", (PC, MatType *), (pc, sub_mat_type));
1804: PetscFunctionReturn(PETSC_SUCCESS);
1805: }
1807: /*@
1808: PCASMSetSubMatType - Set the type of matrix used for `PCASM` subsolves
1810: Collective
1812: Input Parameters:
1813: + pc - the `PC` object
1814: - sub_mat_type - the `MatType`
1816: Options Database Key:
1817: . -pc_asm_sub_mat_type <sub_mat_type> - Sets the matrix type used for subsolves, for example, seqaijviennacl.
1818: If you specify a base name like aijviennacl, the corresponding sequential type is assumed.
1820: Note:
1821: See `MatType` for available types
1823: Level: advanced
1825: .seealso: [](ch_ksp), `PCASM`, `PCASMGetSubMatType()`, `PCSetType()`, `VecSetType()`, `MatType`, `Mat`
1826: @*/
1827: PetscErrorCode PCASMSetSubMatType(PC pc, MatType sub_mat_type)
1828: {
1829: PetscFunctionBegin;
1831: PetscTryMethod(pc, "PCASMSetSubMatType_C", (PC, MatType), (pc, sub_mat_type));
1832: PetscFunctionReturn(PETSC_SUCCESS);
1833: }