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 iascii, isstring;
22: PetscViewer sviewer;
23: PetscViewerFormat format;
24: const char *prefix;
26: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii);
27: PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERSTRING, &isstring);
28: if (iascii) {
29: char overlaps[256] = "user-defined overlap", blocks[256] = "total subdomain blocks not yet set";
30: if (osm->overlap >= 0) PetscSNPrintf(overlaps, sizeof(overlaps), "amount of overlap = %" PetscInt_FMT, osm->overlap);
31: if (osm->n > 0) PetscSNPrintf(blocks, sizeof(blocks), "total subdomain blocks = %" PetscInt_FMT, osm->n);
32: PetscViewerASCIIPrintf(viewer, " %s, %s\n", blocks, overlaps);
33: PetscViewerASCIIPrintf(viewer, " restriction/interpolation type - %s\n", PCASMTypes[osm->type]);
34: if (osm->dm_subdomains) PetscViewerASCIIPrintf(viewer, " Additive Schwarz: using DM to define subdomains\n");
35: if (osm->loctype != PC_COMPOSITE_ADDITIVE) PetscViewerASCIIPrintf(viewer, " Additive Schwarz: local solve composition type - %s\n", PCCompositeTypes[osm->loctype]);
36: MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
37: PetscViewerGetFormat(viewer, &format);
38: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
39: if (osm->ksp) {
40: PetscViewerASCIIPrintf(viewer, " Local solver information for first block is in the following KSP and PC objects on rank 0:\n");
41: PCGetOptionsPrefix(pc, &prefix);
42: PetscViewerASCIIPrintf(viewer, " Use -%sksp_view ::ascii_info_detail to display information for all blocks\n", prefix ? prefix : "");
43: PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
44: if (rank == 0) {
45: PetscViewerASCIIPushTab(viewer);
46: KSPView(osm->ksp[0], sviewer);
47: PetscViewerASCIIPopTab(viewer);
48: }
49: PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
50: }
51: } else {
52: PetscViewerASCIIPushSynchronized(viewer);
53: PetscViewerASCIISynchronizedPrintf(viewer, " [%d] number of local blocks = %" PetscInt_FMT "\n", (int)rank, osm->n_local_true);
54: PetscViewerFlush(viewer);
55: PetscViewerASCIIPrintf(viewer, " Local solver information for each block is in the following KSP and PC objects:\n");
56: PetscViewerASCIIPushTab(viewer);
57: PetscViewerASCIIPrintf(viewer, "- - - - - - - - - - - - - - - - - -\n");
58: PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
59: for (i = 0; i < osm->n_local_true; i++) {
60: ISGetLocalSize(osm->is[i], &bsz);
61: PetscViewerASCIISynchronizedPrintf(sviewer, "[%d] local block number %" PetscInt_FMT ", size = %" PetscInt_FMT "\n", (int)rank, i, bsz);
62: KSPView(osm->ksp[i], sviewer);
63: PetscViewerASCIISynchronizedPrintf(sviewer, "- - - - - - - - - - - - - - - - - -\n");
64: }
65: PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
66: PetscViewerASCIIPopTab(viewer);
67: PetscViewerFlush(viewer);
68: PetscViewerASCIIPopSynchronized(viewer);
69: }
70: } else if (isstring) {
71: PetscViewerStringSPrintf(viewer, " blocks=%" PetscInt_FMT ", overlap=%" PetscInt_FMT ", type=%s", osm->n, osm->overlap, PCASMTypes[osm->type]);
72: PetscViewerGetSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
73: if (osm->ksp) KSPView(osm->ksp[0], sviewer);
74: PetscViewerRestoreSubViewer(viewer, PETSC_COMM_SELF, &sviewer);
75: }
76: return 0;
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: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
91: MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
92: PCGetOptionsPrefix(pc, &prefix);
93: PetscOptionsGetString(NULL, prefix, "-pc_asm_print_subdomains", fname, sizeof(fname), NULL);
94: if (fname[0] == 0) PetscStrcpy(fname, "stdout");
95: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc), fname, &viewer);
96: for (i = 0; i < osm->n_local; i++) {
97: if (i < osm->n_local_true) {
98: ISGetLocalSize(osm->is[i], &nidx);
99: ISGetIndices(osm->is[i], &idx);
100: /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
101: #define len 16 * (nidx + 1) + 512
102: PetscMalloc1(len, &s);
103: PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer);
104: #undef len
105: PetscViewerStringSPrintf(sviewer, "[%d:%d] Subdomain %" PetscInt_FMT " with overlap:\n", rank, size, i);
106: for (j = 0; j < nidx; j++) PetscViewerStringSPrintf(sviewer, "%" PetscInt_FMT " ", idx[j]);
107: ISRestoreIndices(osm->is[i], &idx);
108: PetscViewerStringSPrintf(sviewer, "\n");
109: PetscViewerDestroy(&sviewer);
110: PetscViewerASCIIPushSynchronized(viewer);
111: PetscViewerASCIISynchronizedPrintf(viewer, "%s", s);
112: PetscViewerFlush(viewer);
113: PetscViewerASCIIPopSynchronized(viewer);
114: PetscFree(s);
115: if (osm->is_local) {
116: /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
117: #define len 16 * (nidx + 1) + 512
118: PetscMalloc1(len, &s);
119: PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer);
120: #undef len
121: PetscViewerStringSPrintf(sviewer, "[%d:%d] Subdomain %" PetscInt_FMT " without overlap:\n", rank, size, i);
122: ISGetLocalSize(osm->is_local[i], &nidx);
123: ISGetIndices(osm->is_local[i], &idx);
124: for (j = 0; j < nidx; j++) PetscViewerStringSPrintf(sviewer, "%" PetscInt_FMT " ", idx[j]);
125: ISRestoreIndices(osm->is_local[i], &idx);
126: PetscViewerStringSPrintf(sviewer, "\n");
127: PetscViewerDestroy(&sviewer);
128: PetscViewerASCIIPushSynchronized(viewer);
129: PetscViewerASCIISynchronizedPrintf(viewer, "%s", s);
130: PetscViewerFlush(viewer);
131: PetscViewerASCIIPopSynchronized(viewer);
132: PetscFree(s);
133: }
134: } else {
135: /* Participate in collective viewer calls. */
136: PetscViewerASCIIPushSynchronized(viewer);
137: PetscViewerFlush(viewer);
138: PetscViewerASCIIPopSynchronized(viewer);
139: /* Assume either all ranks have is_local or none do. */
140: if (osm->is_local) {
141: PetscViewerASCIIPushSynchronized(viewer);
142: PetscViewerFlush(viewer);
143: PetscViewerASCIIPopSynchronized(viewer);
144: }
145: }
146: }
147: PetscViewerFlush(viewer);
148: PetscViewerDestroy(&viewer);
149: return 0;
150: }
152: static PetscErrorCode PCSetUp_ASM(PC pc)
153: {
154: PC_ASM *osm = (PC_ASM *)pc->data;
155: PetscBool flg;
156: PetscInt i, m, m_local;
157: MatReuse scall = MAT_REUSE_MATRIX;
158: IS isl;
159: KSP ksp;
160: PC subpc;
161: const char *prefix, *pprefix;
162: Vec vec;
163: DM *domain_dm = NULL;
165: if (!pc->setupcalled) {
166: PetscInt m;
168: /* Note: if subdomains have been set either via PCASMSetTotalSubdomains() or via PCASMSetLocalSubdomains(), osm->n_local_true will not be PETSC_DECIDE */
169: if (osm->n_local_true == PETSC_DECIDE) {
170: /* no subdomains given */
171: /* try pc->dm first, if allowed */
172: if (osm->dm_subdomains && pc->dm) {
173: PetscInt num_domains, d;
174: char **domain_names;
175: IS *inner_domain_is, *outer_domain_is;
176: DMCreateDomainDecomposition(pc->dm, &num_domains, &domain_names, &inner_domain_is, &outer_domain_is, &domain_dm);
177: osm->overlap = -1; /* We do not want to increase the overlap of the IS.
178: A future improvement of this code might allow one to use
179: DM-defined subdomains and also increase the overlap,
180: but that is not currently supported */
181: if (num_domains) PCASMSetLocalSubdomains(pc, num_domains, outer_domain_is, inner_domain_is);
182: for (d = 0; d < num_domains; ++d) {
183: if (domain_names) PetscFree(domain_names[d]);
184: if (inner_domain_is) ISDestroy(&inner_domain_is[d]);
185: if (outer_domain_is) ISDestroy(&outer_domain_is[d]);
186: }
187: PetscFree(domain_names);
188: PetscFree(inner_domain_is);
189: PetscFree(outer_domain_is);
190: }
191: if (osm->n_local_true == PETSC_DECIDE) {
192: /* still no subdomains; use one subdomain per processor */
193: osm->n_local_true = 1;
194: }
195: }
196: { /* determine the global and max number of subdomains */
197: struct {
198: PetscInt max, sum;
199: } inwork, outwork;
200: PetscMPIInt size;
202: inwork.max = osm->n_local_true;
203: inwork.sum = osm->n_local_true;
204: MPIU_Allreduce(&inwork, &outwork, 1, MPIU_2INT, MPIU_MAXSUM_OP, PetscObjectComm((PetscObject)pc));
205: osm->n_local = outwork.max;
206: osm->n = outwork.sum;
208: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
209: if (outwork.max == 1 && outwork.sum == size) {
210: /* osm->n_local_true = 1 on all processes, set this option may enable use of optimized MatCreateSubMatrices() implementation */
211: MatSetOption(pc->pmat, MAT_SUBMAT_SINGLEIS, PETSC_TRUE);
212: }
213: }
214: if (!osm->is) { /* create the index sets */
215: PCASMCreateSubdomains(pc->pmat, osm->n_local_true, &osm->is);
216: }
217: if (osm->n_local_true > 1 && !osm->is_local) {
218: PetscMalloc1(osm->n_local_true, &osm->is_local);
219: for (i = 0; i < osm->n_local_true; i++) {
220: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
221: ISDuplicate(osm->is[i], &osm->is_local[i]);
222: ISCopy(osm->is[i], osm->is_local[i]);
223: } else {
224: PetscObjectReference((PetscObject)osm->is[i]);
225: osm->is_local[i] = osm->is[i];
226: }
227: }
228: }
229: PCGetOptionsPrefix(pc, &prefix);
230: if (osm->overlap > 0) {
231: /* Extend the "overlapping" regions by a number of steps */
232: MatIncreaseOverlap(pc->pmat, osm->n_local_true, osm->is, osm->overlap);
233: }
234: if (osm->sort_indices) {
235: for (i = 0; i < osm->n_local_true; i++) {
236: ISSort(osm->is[i]);
237: if (osm->is_local) ISSort(osm->is_local[i]);
238: }
239: }
240: flg = PETSC_FALSE;
241: PetscOptionsHasName(NULL, prefix, "-pc_asm_print_subdomains", &flg);
242: if (flg) PCASMPrintSubdomains(pc);
243: if (!osm->ksp) {
244: /* Create the local solvers */
245: PetscMalloc1(osm->n_local_true, &osm->ksp);
246: if (domain_dm) PetscInfo(pc, "Setting up ASM subproblems using the embedded DM\n");
247: for (i = 0; i < osm->n_local_true; i++) {
248: KSPCreate(PETSC_COMM_SELF, &ksp);
249: KSPSetErrorIfNotConverged(ksp, pc->erroriffailure);
250: PetscObjectIncrementTabLevel((PetscObject)ksp, (PetscObject)pc, 1);
251: KSPSetType(ksp, KSPPREONLY);
252: KSPGetPC(ksp, &subpc);
253: PCGetOptionsPrefix(pc, &prefix);
254: KSPSetOptionsPrefix(ksp, prefix);
255: KSPAppendOptionsPrefix(ksp, "sub_");
256: if (domain_dm) {
257: KSPSetDM(ksp, domain_dm[i]);
258: KSPSetDMActive(ksp, PETSC_FALSE);
259: DMDestroy(&domain_dm[i]);
260: }
261: osm->ksp[i] = ksp;
262: }
263: if (domain_dm) PetscFree(domain_dm);
264: }
266: ISConcatenate(PETSC_COMM_SELF, osm->n_local_true, osm->is, &osm->lis);
267: ISSortRemoveDups(osm->lis);
268: ISGetLocalSize(osm->lis, &m);
270: scall = MAT_INITIAL_MATRIX;
271: } else {
272: /*
273: Destroy the blocks from the previous iteration
274: */
275: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
276: MatDestroyMatrices(osm->n_local_true, &osm->pmat);
277: scall = MAT_INITIAL_MATRIX;
278: }
279: }
281: /* Destroy previous submatrices of a different type than pc->pmat since MAT_REUSE_MATRIX won't work in that case */
282: if (scall == MAT_REUSE_MATRIX && osm->sub_mat_type) {
283: if (osm->n_local_true > 0) MatDestroySubMatrices(osm->n_local_true, &osm->pmat);
284: scall = MAT_INITIAL_MATRIX;
285: }
287: /*
288: Extract out the submatrices
289: */
290: MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, osm->is, scall, &osm->pmat);
291: if (scall == MAT_INITIAL_MATRIX) {
292: PetscObjectGetOptionsPrefix((PetscObject)pc->pmat, &pprefix);
293: for (i = 0; i < osm->n_local_true; i++) { PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i], pprefix); }
294: }
296: /* Convert the types of the submatrices (if needbe) */
297: if (osm->sub_mat_type) {
298: for (i = 0; i < osm->n_local_true; i++) MatConvert(osm->pmat[i], osm->sub_mat_type, MAT_INPLACE_MATRIX, &(osm->pmat[i]));
299: }
301: if (!pc->setupcalled) {
302: VecType vtype;
304: /* Create the local work vectors (from the local matrices) and scatter contexts */
305: MatCreateVecs(pc->pmat, &vec, NULL);
308: if (osm->is_local && osm->type == PC_ASM_RESTRICT && osm->loctype == PC_COMPOSITE_ADDITIVE) PetscMalloc1(osm->n_local_true, &osm->lprolongation);
309: PetscMalloc1(osm->n_local_true, &osm->lrestriction);
310: PetscMalloc1(osm->n_local_true, &osm->x);
311: PetscMalloc1(osm->n_local_true, &osm->y);
313: ISGetLocalSize(osm->lis, &m);
314: ISCreateStride(PETSC_COMM_SELF, m, 0, 1, &isl);
315: MatGetVecType(osm->pmat[0], &vtype);
316: VecCreate(PETSC_COMM_SELF, &osm->lx);
317: VecSetSizes(osm->lx, m, m);
318: VecSetType(osm->lx, vtype);
319: VecDuplicate(osm->lx, &osm->ly);
320: VecScatterCreate(vec, osm->lis, osm->lx, isl, &osm->restriction);
321: ISDestroy(&isl);
323: for (i = 0; i < osm->n_local_true; ++i) {
324: ISLocalToGlobalMapping ltog;
325: IS isll;
326: const PetscInt *idx_is;
327: PetscInt *idx_lis, nout;
329: ISGetLocalSize(osm->is[i], &m);
330: MatCreateVecs(osm->pmat[i], &osm->x[i], NULL);
331: VecDuplicate(osm->x[i], &osm->y[i]);
333: /* generate a scatter from ly to y[i] picking all the overlapping is[i] entries */
334: ISLocalToGlobalMappingCreateIS(osm->lis, <og);
335: ISGetLocalSize(osm->is[i], &m);
336: ISGetIndices(osm->is[i], &idx_is);
337: PetscMalloc1(m, &idx_lis);
338: ISGlobalToLocalMappingApply(ltog, IS_GTOLM_DROP, m, idx_is, &nout, idx_lis);
340: ISRestoreIndices(osm->is[i], &idx_is);
341: ISCreateGeneral(PETSC_COMM_SELF, m, idx_lis, PETSC_OWN_POINTER, &isll);
342: ISLocalToGlobalMappingDestroy(<og);
343: ISCreateStride(PETSC_COMM_SELF, m, 0, 1, &isl);
344: VecScatterCreate(osm->ly, isll, osm->y[i], isl, &osm->lrestriction[i]);
345: ISDestroy(&isll);
346: ISDestroy(&isl);
347: if (osm->lprolongation) { /* generate a scatter from y[i] to ly picking only the the non-overlapping is_local[i] entries */
348: ISLocalToGlobalMapping ltog;
349: IS isll, isll_local;
350: const PetscInt *idx_local;
351: PetscInt *idx1, *idx2, nout;
353: ISGetLocalSize(osm->is_local[i], &m_local);
354: ISGetIndices(osm->is_local[i], &idx_local);
356: ISLocalToGlobalMappingCreateIS(osm->is[i], <og);
357: PetscMalloc1(m_local, &idx1);
358: ISGlobalToLocalMappingApply(ltog, IS_GTOLM_DROP, m_local, idx_local, &nout, idx1);
359: ISLocalToGlobalMappingDestroy(<og);
361: ISCreateGeneral(PETSC_COMM_SELF, m_local, idx1, PETSC_OWN_POINTER, &isll);
363: ISLocalToGlobalMappingCreateIS(osm->lis, <og);
364: PetscMalloc1(m_local, &idx2);
365: ISGlobalToLocalMappingApply(ltog, IS_GTOLM_DROP, m_local, idx_local, &nout, idx2);
366: ISLocalToGlobalMappingDestroy(<og);
368: ISCreateGeneral(PETSC_COMM_SELF, m_local, idx2, PETSC_OWN_POINTER, &isll_local);
370: ISRestoreIndices(osm->is_local[i], &idx_local);
371: VecScatterCreate(osm->y[i], isll, osm->ly, isll_local, &osm->lprolongation[i]);
373: ISDestroy(&isll);
374: ISDestroy(&isll_local);
375: }
376: }
377: VecDestroy(&vec);
378: }
380: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
381: IS *cis;
382: PetscInt c;
384: PetscMalloc1(osm->n_local_true, &cis);
385: for (c = 0; c < osm->n_local_true; ++c) cis[c] = osm->lis;
386: MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, cis, scall, &osm->lmats);
387: PetscFree(cis);
388: }
390: /* Return control to the user so that the submatrices can be modified (e.g., to apply
391: different boundary conditions for the submatrices than for the global problem) */
392: PCModifySubMatrices(pc, osm->n_local_true, osm->is, osm->is, osm->pmat, pc->modifysubmatricesP);
394: /*
395: Loop over subdomains putting them into local ksp
396: */
397: KSPGetOptionsPrefix(osm->ksp[0], &prefix);
398: for (i = 0; i < osm->n_local_true; i++) {
399: KSPSetOperators(osm->ksp[i], osm->pmat[i], osm->pmat[i]);
400: MatSetOptionsPrefix(osm->pmat[i], prefix);
401: if (!pc->setupcalled) KSPSetFromOptions(osm->ksp[i]);
402: }
403: return 0;
404: }
406: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
407: {
408: PC_ASM *osm = (PC_ASM *)pc->data;
409: PetscInt i;
410: KSPConvergedReason reason;
412: for (i = 0; i < osm->n_local_true; i++) {
413: KSPSetUp(osm->ksp[i]);
414: KSPGetConvergedReason(osm->ksp[i], &reason);
415: if (reason == KSP_DIVERGED_PC_FAILED) pc->failedreason = PC_SUBPC_ERROR;
416: }
417: return 0;
418: }
420: static PetscErrorCode PCApply_ASM(PC pc, Vec x, Vec y)
421: {
422: PC_ASM *osm = (PC_ASM *)pc->data;
423: PetscInt i, n_local_true = osm->n_local_true;
424: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
426: /*
427: support for limiting the restriction or interpolation to only local
428: subdomain values (leaving the other values 0).
429: */
430: if (!(osm->type & PC_ASM_RESTRICT)) {
431: forward = SCATTER_FORWARD_LOCAL;
432: /* have to zero the work RHS since scatter may leave some slots empty */
433: VecSet(osm->lx, 0.0);
434: }
435: if (!(osm->type & PC_ASM_INTERPOLATE)) reverse = SCATTER_REVERSE_LOCAL;
437: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE) {
438: /* zero the global and the local solutions */
439: VecSet(y, 0.0);
440: VecSet(osm->ly, 0.0);
442: /* copy the global RHS to local RHS including the ghost nodes */
443: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
444: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
446: /* restrict local RHS to the overlapping 0-block RHS */
447: VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
448: VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
450: /* do the local solves */
451: for (i = 0; i < n_local_true; ++i) {
452: /* solve the overlapping i-block */
453: PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0);
454: KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]);
455: KSPCheckSolve(osm->ksp[i], pc, osm->y[i]);
456: PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0);
458: if (osm->lprolongation) { /* interpolate the non-overlapping i-block solution to the local solution (only for restrictive additive) */
459: VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
460: VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
461: } else { /* interpolate the overlapping i-block solution to the local solution */
462: VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
463: VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
464: }
466: if (i < n_local_true - 1) {
467: /* restrict local RHS to the overlapping (i+1)-block RHS */
468: VecScatterBegin(osm->lrestriction[i + 1], osm->lx, osm->x[i + 1], INSERT_VALUES, forward);
469: VecScatterEnd(osm->lrestriction[i + 1], osm->lx, osm->x[i + 1], INSERT_VALUES, forward);
471: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
472: /* update the overlapping (i+1)-block RHS using the current local solution */
473: MatMult(osm->lmats[i + 1], osm->ly, osm->y[i + 1]);
474: VecAXPBY(osm->x[i + 1], -1., 1., osm->y[i + 1]);
475: }
476: }
477: }
478: /* add the local solution to the global solution including the ghost nodes */
479: VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
480: VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
481: } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
482: return 0;
483: }
485: static PetscErrorCode PCMatApply_ASM(PC pc, Mat X, Mat Y)
486: {
487: PC_ASM *osm = (PC_ASM *)pc->data;
488: Mat Z, W;
489: Vec x;
490: PetscInt i, m, N;
491: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
494: /*
495: support for limiting the restriction or interpolation to only local
496: subdomain values (leaving the other values 0).
497: */
498: if (!(osm->type & PC_ASM_RESTRICT)) {
499: forward = SCATTER_FORWARD_LOCAL;
500: /* have to zero the work RHS since scatter may leave some slots empty */
501: VecSet(osm->lx, 0.0);
502: }
503: if (!(osm->type & PC_ASM_INTERPOLATE)) reverse = SCATTER_REVERSE_LOCAL;
504: VecGetLocalSize(osm->x[0], &m);
505: MatGetSize(X, NULL, &N);
506: MatCreateSeqDense(PETSC_COMM_SELF, m, N, NULL, &Z);
507: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE) {
508: /* zero the global and the local solutions */
509: MatZeroEntries(Y);
510: VecSet(osm->ly, 0.0);
512: for (i = 0; i < N; ++i) {
513: MatDenseGetColumnVecRead(X, i, &x);
514: /* copy the global RHS to local RHS including the ghost nodes */
515: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
516: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
517: MatDenseRestoreColumnVecRead(X, i, &x);
519: MatDenseGetColumnVecWrite(Z, i, &x);
520: /* restrict local RHS to the overlapping 0-block RHS */
521: VecScatterBegin(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward);
522: VecScatterEnd(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward);
523: MatDenseRestoreColumnVecWrite(Z, i, &x);
524: }
525: MatCreateSeqDense(PETSC_COMM_SELF, m, N, NULL, &W);
526: /* solve the overlapping 0-block */
527: PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[0], Z, W, 0);
528: KSPMatSolve(osm->ksp[0], Z, W);
529: KSPCheckSolve(osm->ksp[0], pc, NULL);
530: PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[0], Z, W, 0);
531: MatDestroy(&Z);
533: for (i = 0; i < N; ++i) {
534: VecSet(osm->ly, 0.0);
535: MatDenseGetColumnVecRead(W, i, &x);
536: if (osm->lprolongation) { /* interpolate the non-overlapping 0-block solution to the local solution (only for restrictive additive) */
537: VecScatterBegin(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward);
538: VecScatterEnd(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward);
539: } else { /* interpolate the overlapping 0-block solution to the local solution */
540: VecScatterBegin(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse);
541: VecScatterEnd(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse);
542: }
543: MatDenseRestoreColumnVecRead(W, i, &x);
545: MatDenseGetColumnVecWrite(Y, i, &x);
546: /* add the local solution to the global solution including the ghost nodes */
547: VecScatterBegin(osm->restriction, osm->ly, x, ADD_VALUES, reverse);
548: VecScatterEnd(osm->restriction, osm->ly, x, ADD_VALUES, reverse);
549: MatDenseRestoreColumnVecWrite(Y, i, &x);
550: }
551: MatDestroy(&W);
552: } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
553: return 0;
554: }
556: static PetscErrorCode PCApplyTranspose_ASM(PC pc, Vec x, Vec y)
557: {
558: PC_ASM *osm = (PC_ASM *)pc->data;
559: PetscInt i, n_local_true = osm->n_local_true;
560: ScatterMode forward = SCATTER_FORWARD, reverse = SCATTER_REVERSE;
562: /*
563: Support for limiting the restriction or interpolation to only local
564: subdomain values (leaving the other values 0).
566: Note: these are reversed from the PCApply_ASM() because we are applying the
567: transpose of the three terms
568: */
570: if (!(osm->type & PC_ASM_INTERPOLATE)) {
571: forward = SCATTER_FORWARD_LOCAL;
572: /* have to zero the work RHS since scatter may leave some slots empty */
573: VecSet(osm->lx, 0.0);
574: }
575: if (!(osm->type & PC_ASM_RESTRICT)) reverse = SCATTER_REVERSE_LOCAL;
577: /* zero the global and the local solutions */
578: VecSet(y, 0.0);
579: VecSet(osm->ly, 0.0);
581: /* Copy the global RHS to local RHS including the ghost nodes */
582: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
583: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
585: /* Restrict local RHS to the overlapping 0-block RHS */
586: VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
587: VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
589: /* do the local solves */
590: for (i = 0; i < n_local_true; ++i) {
591: /* solve the overlapping i-block */
592: PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0);
593: KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i]);
594: KSPCheckSolve(osm->ksp[i], pc, osm->y[i]);
595: PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0);
597: if (osm->lprolongation) { /* interpolate the non-overlapping i-block solution to the local solution */
598: VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
599: VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
600: } else { /* interpolate the overlapping i-block solution to the local solution */
601: VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
602: VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
603: }
605: if (i < n_local_true - 1) {
606: /* Restrict local RHS to the overlapping (i+1)-block RHS */
607: VecScatterBegin(osm->lrestriction[i + 1], osm->lx, osm->x[i + 1], INSERT_VALUES, forward);
608: VecScatterEnd(osm->lrestriction[i + 1], osm->lx, osm->x[i + 1], INSERT_VALUES, forward);
609: }
610: }
611: /* Add the local solution to the global solution including the ghost nodes */
612: VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
613: VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
614: return 0;
615: }
617: static PetscErrorCode PCReset_ASM(PC pc)
618: {
619: PC_ASM *osm = (PC_ASM *)pc->data;
620: PetscInt i;
622: if (osm->ksp) {
623: for (i = 0; i < osm->n_local_true; i++) KSPReset(osm->ksp[i]);
624: }
625: if (osm->pmat) {
626: if (osm->n_local_true > 0) MatDestroySubMatrices(osm->n_local_true, &osm->pmat);
627: }
628: if (osm->lrestriction) {
629: VecScatterDestroy(&osm->restriction);
630: for (i = 0; i < osm->n_local_true; i++) {
631: VecScatterDestroy(&osm->lrestriction[i]);
632: if (osm->lprolongation) VecScatterDestroy(&osm->lprolongation[i]);
633: VecDestroy(&osm->x[i]);
634: VecDestroy(&osm->y[i]);
635: }
636: PetscFree(osm->lrestriction);
637: if (osm->lprolongation) PetscFree(osm->lprolongation);
638: PetscFree(osm->x);
639: PetscFree(osm->y);
640: }
641: PCASMDestroySubdomains(osm->n_local_true, osm->is, osm->is_local);
642: ISDestroy(&osm->lis);
643: VecDestroy(&osm->lx);
644: VecDestroy(&osm->ly);
645: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) MatDestroyMatrices(osm->n_local_true, &osm->lmats);
647: PetscFree(osm->sub_mat_type);
649: osm->is = NULL;
650: osm->is_local = NULL;
651: return 0;
652: }
654: static PetscErrorCode PCDestroy_ASM(PC pc)
655: {
656: PC_ASM *osm = (PC_ASM *)pc->data;
657: PetscInt i;
659: PCReset_ASM(pc);
660: if (osm->ksp) {
661: for (i = 0; i < osm->n_local_true; i++) KSPDestroy(&osm->ksp[i]);
662: PetscFree(osm->ksp);
663: }
664: PetscFree(pc->data);
666: PetscObjectComposeFunction((PetscObject)pc, "PCASMSetLocalSubdomains_C", NULL);
667: PetscObjectComposeFunction((PetscObject)pc, "PCASMSetTotalSubdomains_C", NULL);
668: PetscObjectComposeFunction((PetscObject)pc, "PCASMSetOverlap_C", NULL);
669: PetscObjectComposeFunction((PetscObject)pc, "PCASMSetType_C", NULL);
670: PetscObjectComposeFunction((PetscObject)pc, "PCASMGetType_C", NULL);
671: PetscObjectComposeFunction((PetscObject)pc, "PCASMSetLocalType_C", NULL);
672: PetscObjectComposeFunction((PetscObject)pc, "PCASMGetLocalType_C", NULL);
673: PetscObjectComposeFunction((PetscObject)pc, "PCASMSetSortIndices_C", NULL);
674: PetscObjectComposeFunction((PetscObject)pc, "PCASMGetSubKSP_C", NULL);
675: PetscObjectComposeFunction((PetscObject)pc, "PCASMGetSubMatType_C", NULL);
676: PetscObjectComposeFunction((PetscObject)pc, "PCASMSetSubMatType_C", NULL);
677: return 0;
678: }
680: static PetscErrorCode PCSetFromOptions_ASM(PC pc, PetscOptionItems *PetscOptionsObject)
681: {
682: PC_ASM *osm = (PC_ASM *)pc->data;
683: PetscInt blocks, ovl;
684: PetscBool flg;
685: PCASMType asmtype;
686: PCCompositeType loctype;
687: char sub_mat_type[256];
689: PetscOptionsHeadBegin(PetscOptionsObject, "Additive Schwarz options");
690: PetscOptionsBool("-pc_asm_dm_subdomains", "Use DMCreateDomainDecomposition() to define subdomains", "PCASMSetDMSubdomains", osm->dm_subdomains, &osm->dm_subdomains, &flg);
691: PetscOptionsInt("-pc_asm_blocks", "Number of subdomains", "PCASMSetTotalSubdomains", osm->n, &blocks, &flg);
692: if (flg) {
693: PCASMSetTotalSubdomains(pc, blocks, NULL, NULL);
694: osm->dm_subdomains = PETSC_FALSE;
695: }
696: PetscOptionsInt("-pc_asm_local_blocks", "Number of local subdomains", "PCASMSetLocalSubdomains", osm->n_local_true, &blocks, &flg);
697: if (flg) {
698: PCASMSetLocalSubdomains(pc, blocks, NULL, NULL);
699: osm->dm_subdomains = PETSC_FALSE;
700: }
701: PetscOptionsInt("-pc_asm_overlap", "Number of grid points overlap", "PCASMSetOverlap", osm->overlap, &ovl, &flg);
702: if (flg) {
703: PCASMSetOverlap(pc, ovl);
704: osm->dm_subdomains = PETSC_FALSE;
705: }
706: flg = PETSC_FALSE;
707: PetscOptionsEnum("-pc_asm_type", "Type of restriction/extension", "PCASMSetType", PCASMTypes, (PetscEnum)osm->type, (PetscEnum *)&asmtype, &flg);
708: if (flg) PCASMSetType(pc, asmtype);
709: flg = PETSC_FALSE;
710: PetscOptionsEnum("-pc_asm_local_type", "Type of local solver composition", "PCASMSetLocalType", PCCompositeTypes, (PetscEnum)osm->loctype, (PetscEnum *)&loctype, &flg);
711: if (flg) PCASMSetLocalType(pc, loctype);
712: PetscOptionsFList("-pc_asm_sub_mat_type", "Subsolve Matrix Type", "PCASMSetSubMatType", MatList, NULL, sub_mat_type, 256, &flg);
713: if (flg) PCASMSetSubMatType(pc, sub_mat_type);
714: PetscOptionsHeadEnd();
715: return 0;
716: }
718: static PetscErrorCode PCASMSetLocalSubdomains_ASM(PC pc, PetscInt n, IS is[], IS is_local[])
719: {
720: PC_ASM *osm = (PC_ASM *)pc->data;
721: PetscInt i;
726: if (!pc->setupcalled) {
727: if (is) {
728: for (i = 0; i < n; i++) PetscObjectReference((PetscObject)is[i]);
729: }
730: if (is_local) {
731: for (i = 0; i < n; i++) PetscObjectReference((PetscObject)is_local[i]);
732: }
733: PCASMDestroySubdomains(osm->n_local_true, osm->is, osm->is_local);
735: osm->n_local_true = n;
736: osm->is = NULL;
737: osm->is_local = NULL;
738: if (is) {
739: PetscMalloc1(n, &osm->is);
740: for (i = 0; i < n; i++) osm->is[i] = is[i];
741: /* Flag indicating that the user has set overlapping subdomains so PCASM should not increase their size. */
742: osm->overlap = -1;
743: }
744: if (is_local) {
745: PetscMalloc1(n, &osm->is_local);
746: for (i = 0; i < n; i++) osm->is_local[i] = is_local[i];
747: if (!is) {
748: PetscMalloc1(osm->n_local_true, &osm->is);
749: for (i = 0; i < osm->n_local_true; i++) {
750: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
751: ISDuplicate(osm->is_local[i], &osm->is[i]);
752: ISCopy(osm->is_local[i], osm->is[i]);
753: } else {
754: PetscObjectReference((PetscObject)osm->is_local[i]);
755: osm->is[i] = osm->is_local[i];
756: }
757: }
758: }
759: }
760: }
761: return 0;
762: }
764: static PetscErrorCode PCASMSetTotalSubdomains_ASM(PC pc, PetscInt N, IS *is, IS *is_local)
765: {
766: PC_ASM *osm = (PC_ASM *)pc->data;
767: PetscMPIInt rank, size;
768: PetscInt n;
773: /*
774: Split the subdomains equally among all processors
775: */
776: MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
777: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
778: n = N / size + ((N % size) > rank);
781: if (!pc->setupcalled) {
782: PCASMDestroySubdomains(osm->n_local_true, osm->is, osm->is_local);
784: osm->n_local_true = n;
785: osm->is = NULL;
786: osm->is_local = NULL;
787: }
788: return 0;
789: }
791: static PetscErrorCode PCASMSetOverlap_ASM(PC pc, PetscInt ovl)
792: {
793: PC_ASM *osm = (PC_ASM *)pc->data;
797: if (!pc->setupcalled) osm->overlap = ovl;
798: return 0;
799: }
801: static PetscErrorCode PCASMSetType_ASM(PC pc, PCASMType type)
802: {
803: PC_ASM *osm = (PC_ASM *)pc->data;
805: osm->type = type;
806: osm->type_set = PETSC_TRUE;
807: return 0;
808: }
810: static PetscErrorCode PCASMGetType_ASM(PC pc, PCASMType *type)
811: {
812: PC_ASM *osm = (PC_ASM *)pc->data;
814: *type = osm->type;
815: return 0;
816: }
818: static PetscErrorCode PCASMSetLocalType_ASM(PC pc, PCCompositeType type)
819: {
820: PC_ASM *osm = (PC_ASM *)pc->data;
823: osm->loctype = type;
824: return 0;
825: }
827: static PetscErrorCode PCASMGetLocalType_ASM(PC pc, PCCompositeType *type)
828: {
829: PC_ASM *osm = (PC_ASM *)pc->data;
831: *type = osm->loctype;
832: return 0;
833: }
835: static PetscErrorCode PCASMSetSortIndices_ASM(PC pc, PetscBool doSort)
836: {
837: PC_ASM *osm = (PC_ASM *)pc->data;
839: osm->sort_indices = doSort;
840: return 0;
841: }
843: static PetscErrorCode PCASMGetSubKSP_ASM(PC pc, PetscInt *n_local, PetscInt *first_local, KSP **ksp)
844: {
845: PC_ASM *osm = (PC_ASM *)pc->data;
849: if (n_local) *n_local = osm->n_local_true;
850: if (first_local) {
851: MPI_Scan(&osm->n_local_true, first_local, 1, MPIU_INT, MPI_SUM, PetscObjectComm((PetscObject)pc));
852: *first_local -= osm->n_local_true;
853: }
854: if (ksp) *ksp = osm->ksp;
855: return 0;
856: }
858: static PetscErrorCode PCASMGetSubMatType_ASM(PC pc, MatType *sub_mat_type)
859: {
860: PC_ASM *osm = (PC_ASM *)pc->data;
864: *sub_mat_type = osm->sub_mat_type;
865: return 0;
866: }
868: static PetscErrorCode PCASMSetSubMatType_ASM(PC pc, MatType sub_mat_type)
869: {
870: PC_ASM *osm = (PC_ASM *)pc->data;
873: PetscFree(osm->sub_mat_type);
874: PetscStrallocpy(sub_mat_type, (char **)&osm->sub_mat_type);
875: return 0;
876: }
878: /*@C
879: PCASMSetLocalSubdomains - Sets the local subdomains (for this processor only) for the additive Schwarz preconditioner `PCASM`.
881: Collective
883: Input Parameters:
884: + pc - the preconditioner context
885: . n - the number of subdomains for this processor (default value = 1)
886: . is - the index set that defines the subdomains for this processor
887: (or NULL for PETSc to determine subdomains)
888: - is_local - the index sets that define the local part of the subdomains for this processor, not used unless PCASMType is PC_ASM_RESTRICT
889: (or NULL to not provide these)
891: Options Database Key:
892: To set the total number of subdomain blocks rather than specify the
893: index sets, use the option
894: . -pc_asm_local_blocks <blks> - Sets local blocks
896: Notes:
897: The `IS` numbering is in the parallel, global numbering of the vector for both is and is_local
899: By default the `PCASM` preconditioner uses 1 block per processor.
901: Use `PCASMSetTotalSubdomains()` to set the subdomains for all processors.
903: If is_local is provided and `PCASMType` is `PC_ASM_RESTRICT` then the solution only over the is_local region is interpolated
904: back to form the global solution (this is the standard restricted additive Schwarz method)
906: If the is_local is provided and `PCASMType` is `PC_ASM_INTERPOLATE` or `PC_ASM_NONE` then an error is generated since there is
907: no code to handle that case.
909: Level: advanced
911: .seealso: `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`, `PCASMGetSubKSP()`,
912: `PCASMCreateSubdomains2D()`, `PCASMGetLocalSubdomains()`, `PCASMType`, `PCASMSetType()`, `PCGASM`
913: @*/
914: PetscErrorCode PCASMSetLocalSubdomains(PC pc, PetscInt n, IS is[], IS is_local[])
915: {
917: PetscTryMethod(pc, "PCASMSetLocalSubdomains_C", (PC, PetscInt, IS[], IS[]), (pc, n, is, is_local));
918: return 0;
919: }
921: /*@C
922: PCASMSetTotalSubdomains - Sets the subdomains for all processors for the
923: additive Schwarz preconditioner, `PCASM`.
925: Collective, all MPI ranks must pass in the same array of `IS`
927: Input Parameters:
928: + pc - the preconditioner context
929: . N - the number of subdomains for all processors
930: . is - the index sets that define the subdomains for all processors
931: (or NULL to ask PETSc to determine the subdomains)
932: - is_local - the index sets that define the local part of the subdomains for this processor
933: (or NULL to not provide this information)
935: Options Database Key:
936: To set the total number of subdomain blocks rather than specify the
937: index sets, use the option
938: . -pc_asm_blocks <blks> - Sets total blocks
940: Notes:
941: Currently you cannot use this to set the actual subdomains with the argument is or is_local.
943: By default the `PCASM` preconditioner uses 1 block per processor.
945: These index sets cannot be destroyed until after completion of the
946: linear solves for which the `PCASM` preconditioner is being used.
948: Use `PCASMSetLocalSubdomains()` to set local subdomains.
950: The `IS` numbering is in the parallel, global numbering of the vector for both is and is_local
952: Level: advanced
954: .seealso: `PCASM`, `PCASMSetLocalSubdomains()`, `PCASMSetOverlap()`, `PCASMGetSubKSP()`,
955: `PCASMCreateSubdomains2D()`, `PCGASM`
956: @*/
957: PetscErrorCode PCASMSetTotalSubdomains(PC pc, PetscInt N, IS is[], IS is_local[])
958: {
960: PetscTryMethod(pc, "PCASMSetTotalSubdomains_C", (PC, PetscInt, IS[], IS[]), (pc, N, is, is_local));
961: return 0;
962: }
964: /*@
965: PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
966: additive Schwarz preconditioner, `PCASM`.
968: Logically Collective
970: Input Parameters:
971: + pc - the preconditioner context
972: - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1)
974: Options Database Key:
975: . -pc_asm_overlap <ovl> - Sets overlap
977: Notes:
978: By default the `PCASM` preconditioner uses 1 block per processor. To use
979: multiple blocks per perocessor, see `PCASMSetTotalSubdomains()` and
980: `PCASMSetLocalSubdomains()` (and the option -pc_asm_blocks <blks>).
982: The overlap defaults to 1, so if one desires that no additional
983: overlap be computed beyond what may have been set with a call to
984: `PCASMSetTotalSubdomains()` or `PCASMSetLocalSubdomains()`, then ovl
985: must be set to be 0. In particular, if one does not explicitly set
986: the subdomains an application code, then all overlap would be computed
987: internally by PETSc, and using an overlap of 0 would result in an `PCASM`
988: variant that is equivalent to the block Jacobi preconditioner.
990: The default algorithm used by PETSc to increase overlap is fast, but not scalable,
991: use the option -mat_increase_overlap_scalable when the problem and number of processes is large.
993: Note that one can define initial index sets with any overlap via
994: `PCASMSetLocalSubdomains()`; the routine
995: `PCASMSetOverlap()` merely allows PETSc to extend that overlap further
996: if desired.
998: Level: intermediate
1000: .seealso: `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetLocalSubdomains()`, `PCASMGetSubKSP()`,
1001: `PCASMCreateSubdomains2D()`, `PCASMGetLocalSubdomains()`, `MatIncreaseOverlap()`, `PCGASM`
1002: @*/
1003: PetscErrorCode PCASMSetOverlap(PC pc, PetscInt ovl)
1004: {
1007: PetscTryMethod(pc, "PCASMSetOverlap_C", (PC, PetscInt), (pc, ovl));
1008: return 0;
1009: }
1011: /*@
1012: PCASMSetType - Sets the type of restriction and interpolation used
1013: for local problems in the additive Schwarz method, `PCASM`.
1015: Logically Collective
1017: Input Parameters:
1018: + pc - the preconditioner context
1019: - type - variant of `PCASM`, one of
1020: .vb
1021: PC_ASM_BASIC - full interpolation and restriction
1022: PC_ASM_RESTRICT - full restriction, local processor interpolation (default)
1023: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1024: PC_ASM_NONE - local processor restriction and interpolation
1025: .ve
1027: Options Database Key:
1028: . -pc_asm_type [basic,restrict,interpolate,none] - Sets `PCASMType`
1030: Note:
1031: if the is_local arguments are passed to `PCASMSetLocalSubdomains()` then they are used when `PC_ASM_RESTRICT` has been selected
1032: to limit the local processor interpolation
1034: Level: intermediate
1036: .seealso: `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetTotalSubdomains()`, `PCASMGetSubKSP()`,
1037: `PCASMCreateSubdomains2D()`, `PCASMType`, `PCASMSetLocalType()`, `PCASMGetLocalType()`, `PCGASM`
1038: @*/
1039: PetscErrorCode PCASMSetType(PC pc, PCASMType type)
1040: {
1043: PetscTryMethod(pc, "PCASMSetType_C", (PC, PCASMType), (pc, type));
1044: return 0;
1045: }
1047: /*@
1048: PCASMGetType - Gets the type of restriction and interpolation used
1049: for local problems in the additive Schwarz method, `PCASM`.
1051: Logically Collective
1053: Input Parameter:
1054: . pc - the preconditioner context
1056: Output Parameter:
1057: . type - variant of `PCASM`, one of
1059: .vb
1060: PC_ASM_BASIC - full interpolation and restriction
1061: PC_ASM_RESTRICT - full restriction, local processor interpolation
1062: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1063: PC_ASM_NONE - local processor restriction and interpolation
1064: .ve
1066: Options Database Key:
1067: . -pc_asm_type [basic,restrict,interpolate,none] - Sets `PCASM` type
1069: Level: intermediate
1071: .seealso: `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetTotalSubdomains()`, `PCASMGetSubKSP()`, `PCGASM`,
1072: `PCASMCreateSubdomains2D()`, `PCASMType`, `PCASMSetType()`, `PCASMSetLocalType()`, `PCASMGetLocalType()`
1073: @*/
1074: PetscErrorCode PCASMGetType(PC pc, PCASMType *type)
1075: {
1077: PetscUseMethod(pc, "PCASMGetType_C", (PC, PCASMType *), (pc, type));
1078: return 0;
1079: }
1081: /*@
1082: PCASMSetLocalType - Sets the type of composition used for local problems in the additive Schwarz method, `PCASM`.
1084: Logically Collective
1086: Input Parameters:
1087: + pc - the preconditioner context
1088: - type - type of composition, one of
1089: .vb
1090: PC_COMPOSITE_ADDITIVE - local additive combination
1091: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1092: .ve
1094: Options Database Key:
1095: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1097: Level: intermediate
1099: .seealso: `PCASM`, `PCASMSetType()`, `PCASMGetType()`, `PCASMGetLocalType()`, `PCASM`, `PCASMType`, `PCASMSetType()`, `PCASMGetType()`, `PCCompositeType`
1100: @*/
1101: PetscErrorCode PCASMSetLocalType(PC pc, PCCompositeType type)
1102: {
1105: PetscTryMethod(pc, "PCASMSetLocalType_C", (PC, PCCompositeType), (pc, type));
1106: return 0;
1107: }
1109: /*@
1110: PCASMGetLocalType - Gets the type of composition used for local problems in the additive Schwarz method, `PCASM`.
1112: Logically Collective
1114: Input Parameter:
1115: . pc - the preconditioner context
1117: Output Parameter:
1118: . type - type of composition, one of
1119: .vb
1120: PC_COMPOSITE_ADDITIVE - local additive combination
1121: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1122: .ve
1124: Options Database Key:
1125: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1127: Level: intermediate
1129: .seealso: `PCASM`, `PCASMSetType()`, `PCASMGetType()`, `PCASMSetLocalType()`, `PCASMCreate()`, `PCASMType`, `PCASMSetType()`, `PCASMGetType()`, `PCCompositeType`
1130: @*/
1131: PetscErrorCode PCASMGetLocalType(PC pc, PCCompositeType *type)
1132: {
1135: PetscUseMethod(pc, "PCASMGetLocalType_C", (PC, PCCompositeType *), (pc, type));
1136: return 0;
1137: }
1139: /*@
1140: PCASMSetSortIndices - Determines whether subdomain indices are sorted.
1142: Logically Collective
1144: Input Parameters:
1145: + pc - the preconditioner context
1146: - doSort - sort the subdomain indices
1148: Level: intermediate
1150: .seealso: `PCASM`, `PCASMSetLocalSubdomains()`, `PCASMSetTotalSubdomains()`, `PCASMGetSubKSP()`,
1151: `PCASMCreateSubdomains2D()`
1152: @*/
1153: PetscErrorCode PCASMSetSortIndices(PC pc, PetscBool doSort)
1154: {
1157: PetscTryMethod(pc, "PCASMSetSortIndices_C", (PC, PetscBool), (pc, doSort));
1158: return 0;
1159: }
1161: /*@C
1162: PCASMGetSubKSP - Gets the local `KSP` contexts for all blocks on
1163: this processor.
1165: Collective iff first_local is requested
1167: Input Parameter:
1168: . pc - the preconditioner context
1170: Output Parameters:
1171: + n_local - the number of blocks on this processor or NULL
1172: . first_local - the global number of the first block on this processor or NULL,
1173: all processors must request or all must pass NULL
1174: - ksp - the array of `KSP` contexts
1176: Notes:
1177: After `PCASMGetSubKSP()` the array of `KSP`s is not to be freed.
1179: You must call `KSPSetUp()` before calling `PCASMGetSubKSP()`.
1181: Fortran Note:
1182: The output argument 'ksp' must be an array of sufficient length or `PETSC_NULL_KSP`. The latter can be used to learn the necessary length.
1184: Level: advanced
1186: .seealso: `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`,
1187: `PCASMCreateSubdomains2D()`,
1188: @*/
1189: PetscErrorCode PCASMGetSubKSP(PC pc, PetscInt *n_local, PetscInt *first_local, KSP *ksp[])
1190: {
1192: PetscUseMethod(pc, "PCASMGetSubKSP_C", (PC, PetscInt *, PetscInt *, KSP **), (pc, n_local, first_local, ksp));
1193: return 0;
1194: }
1196: /*MC
1197: PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1198: its own `KSP` object.
1200: Options Database Keys:
1201: + -pc_asm_blocks <blks> - Sets total blocks. Defaults to one block per MPI rank.
1202: . -pc_asm_overlap <ovl> - Sets overlap
1203: . -pc_asm_type [basic,restrict,interpolate,none] - Sets `PCASMType`, default is restrict. See `PCASMSetType()`
1204: - -pc_asm_local_type [additive, multiplicative] - Sets `PCCompositeType`, default is additive. See `PCASMSetLocalType()`
1206: Level: beginner
1208: Notes:
1209: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you
1210: will get a different convergence rate due to the default option of -pc_asm_type restrict. Use
1211: -pc_asm_type basic to get the same convergence behavior
1213: Each processor can have one or more blocks, but a block cannot be shared by more
1214: than one processor. Use `PCGASM` for subdomains shared by multiple processes.
1216: To set options on the solvers for each block append -sub_ to all the `KSP`, and `PC`
1217: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
1219: To set the options on the solvers separate for each block call `PCASMGetSubKSP()`
1220: and set the options directly on the resulting `KSP` object (you can access its `PC` with `KSPGetPC()`)
1222: References:
1223: + * - M Dryja, OB Widlund, An additive variant of the Schwarz alternating method for the case of many subregions
1224: Courant Institute, New York University Technical report
1225: - * - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations,
1226: Cambridge University Press.
1228: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCASMType`, `PCCompositeType`,
1229: `PCBJACOBI`, `PCASMGetSubKSP()`, `PCASMSetLocalSubdomains()`, `PCASMType`, `PCASMGetType()`, `PCASMSetLocalType()`, `PCASMGetLocalType()`
1230: `PCASMSetTotalSubdomains()`, `PCSetModifySubMatrices()`, `PCASMSetOverlap()`, `PCASMSetType()`, `PCCompositeType`
1231: M*/
1233: PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1234: {
1235: PC_ASM *osm;
1237: PetscNew(&osm);
1239: osm->n = PETSC_DECIDE;
1240: osm->n_local = 0;
1241: osm->n_local_true = PETSC_DECIDE;
1242: osm->overlap = 1;
1243: osm->ksp = NULL;
1244: osm->restriction = NULL;
1245: osm->lprolongation = NULL;
1246: osm->lrestriction = NULL;
1247: osm->x = NULL;
1248: osm->y = NULL;
1249: osm->is = NULL;
1250: osm->is_local = NULL;
1251: osm->mat = NULL;
1252: osm->pmat = NULL;
1253: osm->type = PC_ASM_RESTRICT;
1254: osm->loctype = PC_COMPOSITE_ADDITIVE;
1255: osm->sort_indices = PETSC_TRUE;
1256: osm->dm_subdomains = PETSC_FALSE;
1257: osm->sub_mat_type = NULL;
1259: pc->data = (void *)osm;
1260: pc->ops->apply = PCApply_ASM;
1261: pc->ops->matapply = PCMatApply_ASM;
1262: pc->ops->applytranspose = PCApplyTranspose_ASM;
1263: pc->ops->setup = PCSetUp_ASM;
1264: pc->ops->reset = PCReset_ASM;
1265: pc->ops->destroy = PCDestroy_ASM;
1266: pc->ops->setfromoptions = PCSetFromOptions_ASM;
1267: pc->ops->setuponblocks = PCSetUpOnBlocks_ASM;
1268: pc->ops->view = PCView_ASM;
1269: pc->ops->applyrichardson = NULL;
1271: PetscObjectComposeFunction((PetscObject)pc, "PCASMSetLocalSubdomains_C", PCASMSetLocalSubdomains_ASM);
1272: PetscObjectComposeFunction((PetscObject)pc, "PCASMSetTotalSubdomains_C", PCASMSetTotalSubdomains_ASM);
1273: PetscObjectComposeFunction((PetscObject)pc, "PCASMSetOverlap_C", PCASMSetOverlap_ASM);
1274: PetscObjectComposeFunction((PetscObject)pc, "PCASMSetType_C", PCASMSetType_ASM);
1275: PetscObjectComposeFunction((PetscObject)pc, "PCASMGetType_C", PCASMGetType_ASM);
1276: PetscObjectComposeFunction((PetscObject)pc, "PCASMSetLocalType_C", PCASMSetLocalType_ASM);
1277: PetscObjectComposeFunction((PetscObject)pc, "PCASMGetLocalType_C", PCASMGetLocalType_ASM);
1278: PetscObjectComposeFunction((PetscObject)pc, "PCASMSetSortIndices_C", PCASMSetSortIndices_ASM);
1279: PetscObjectComposeFunction((PetscObject)pc, "PCASMGetSubKSP_C", PCASMGetSubKSP_ASM);
1280: PetscObjectComposeFunction((PetscObject)pc, "PCASMGetSubMatType_C", PCASMGetSubMatType_ASM);
1281: PetscObjectComposeFunction((PetscObject)pc, "PCASMSetSubMatType_C", PCASMSetSubMatType_ASM);
1282: return 0;
1283: }
1285: /*@C
1286: PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1287: preconditioner, `PCASM`, for any problem on a general grid.
1289: Collective
1291: Input Parameters:
1292: + A - The global matrix operator
1293: - n - the number of local blocks
1295: Output Parameters:
1296: . outis - the array of index sets defining the subdomains
1298: Level: advanced
1300: Note:
1301: This generates nonoverlapping subdomains; the `PCASM` will generate the overlap
1302: from these if you use `PCASMSetLocalSubdomains()`
1304: Fortran Note:
1305: You must provide the array outis[] already allocated of length n.
1307: .seealso: `PCASM`, `PCASMSetLocalSubdomains()`, `PCASMDestroySubdomains()`
1308: @*/
1309: PetscErrorCode PCASMCreateSubdomains(Mat A, PetscInt n, IS *outis[])
1310: {
1311: MatPartitioning mpart;
1312: const char *prefix;
1313: PetscInt i, j, rstart, rend, bs;
1314: PetscBool hasop, isbaij = PETSC_FALSE, foundpart = PETSC_FALSE;
1315: Mat Ad = NULL, adj;
1316: IS ispart, isnumb, *is;
1322: /* Get prefix, row distribution, and block size */
1323: MatGetOptionsPrefix(A, &prefix);
1324: MatGetOwnershipRange(A, &rstart, &rend);
1325: MatGetBlockSize(A, &bs);
1328: /* Get diagonal block from matrix if possible */
1329: MatHasOperation(A, MATOP_GET_DIAGONAL_BLOCK, &hasop);
1330: if (hasop) MatGetDiagonalBlock(A, &Ad);
1331: if (Ad) {
1332: PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQBAIJ, &isbaij);
1333: if (!isbaij) PetscObjectBaseTypeCompare((PetscObject)Ad, MATSEQSBAIJ, &isbaij);
1334: }
1335: if (Ad && n > 1) {
1336: PetscBool match, done;
1337: /* Try to setup a good matrix partitioning if available */
1338: MatPartitioningCreate(PETSC_COMM_SELF, &mpart);
1339: PetscObjectSetOptionsPrefix((PetscObject)mpart, prefix);
1340: MatPartitioningSetFromOptions(mpart);
1341: PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGCURRENT, &match);
1342: if (!match) PetscObjectTypeCompare((PetscObject)mpart, MATPARTITIONINGSQUARE, &match);
1343: if (!match) { /* assume a "good" partitioner is available */
1344: PetscInt na;
1345: const PetscInt *ia, *ja;
1346: MatGetRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done);
1347: if (done) {
1348: /* Build adjacency matrix by hand. Unfortunately a call to
1349: MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1350: remove the block-aij structure and we cannot expect
1351: MatPartitioning to split vertices as we need */
1352: PetscInt i, j, len, nnz, cnt, *iia = NULL, *jja = NULL;
1353: const PetscInt *row;
1354: nnz = 0;
1355: for (i = 0; i < na; i++) { /* count number of nonzeros */
1356: len = ia[i + 1] - ia[i];
1357: row = ja + ia[i];
1358: for (j = 0; j < len; j++) {
1359: if (row[j] == i) { /* don't count diagonal */
1360: len--;
1361: break;
1362: }
1363: }
1364: nnz += len;
1365: }
1366: PetscMalloc1(na + 1, &iia);
1367: PetscMalloc1(nnz, &jja);
1368: nnz = 0;
1369: iia[0] = 0;
1370: for (i = 0; i < na; i++) { /* fill adjacency */
1371: cnt = 0;
1372: len = ia[i + 1] - ia[i];
1373: row = ja + ia[i];
1374: for (j = 0; j < len; j++) {
1375: if (row[j] != i) { /* if not diagonal */
1376: jja[nnz + cnt++] = row[j];
1377: }
1378: }
1379: nnz += cnt;
1380: iia[i + 1] = nnz;
1381: }
1382: /* Partitioning of the adjacency matrix */
1383: MatCreateMPIAdj(PETSC_COMM_SELF, na, na, iia, jja, NULL, &adj);
1384: MatPartitioningSetAdjacency(mpart, adj);
1385: MatPartitioningSetNParts(mpart, n);
1386: MatPartitioningApply(mpart, &ispart);
1387: ISPartitioningToNumbering(ispart, &isnumb);
1388: MatDestroy(&adj);
1389: foundpart = PETSC_TRUE;
1390: }
1391: MatRestoreRowIJ(Ad, 0, PETSC_TRUE, isbaij, &na, &ia, &ja, &done);
1392: }
1393: MatPartitioningDestroy(&mpart);
1394: }
1396: PetscMalloc1(n, &is);
1397: *outis = is;
1399: if (!foundpart) {
1400: /* Partitioning by contiguous chunks of rows */
1402: PetscInt mbs = (rend - rstart) / bs;
1403: PetscInt start = rstart;
1404: for (i = 0; i < n; i++) {
1405: PetscInt count = (mbs / n + ((mbs % n) > i)) * bs;
1406: ISCreateStride(PETSC_COMM_SELF, count, start, 1, &is[i]);
1407: start += count;
1408: }
1410: } else {
1411: /* Partitioning by adjacency of diagonal block */
1413: const PetscInt *numbering;
1414: PetscInt *count, nidx, *indices, *newidx, start = 0;
1415: /* Get node count in each partition */
1416: PetscMalloc1(n, &count);
1417: ISPartitioningCount(ispart, n, count);
1418: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1419: for (i = 0; i < n; i++) count[i] *= bs;
1420: }
1421: /* Build indices from node numbering */
1422: ISGetLocalSize(isnumb, &nidx);
1423: PetscMalloc1(nidx, &indices);
1424: for (i = 0; i < nidx; i++) indices[i] = i; /* needs to be initialized */
1425: ISGetIndices(isnumb, &numbering);
1426: PetscSortIntWithPermutation(nidx, numbering, indices);
1427: ISRestoreIndices(isnumb, &numbering);
1428: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1429: PetscMalloc1(nidx * bs, &newidx);
1430: for (i = 0; i < nidx; i++) {
1431: for (j = 0; j < bs; j++) newidx[i * bs + j] = indices[i] * bs + j;
1432: }
1433: PetscFree(indices);
1434: nidx *= bs;
1435: indices = newidx;
1436: }
1437: /* Shift to get global indices */
1438: for (i = 0; i < nidx; i++) indices[i] += rstart;
1440: /* Build the index sets for each block */
1441: for (i = 0; i < n; i++) {
1442: ISCreateGeneral(PETSC_COMM_SELF, count[i], &indices[start], PETSC_COPY_VALUES, &is[i]);
1443: ISSort(is[i]);
1444: start += count[i];
1445: }
1447: PetscFree(count);
1448: PetscFree(indices);
1449: ISDestroy(&isnumb);
1450: ISDestroy(&ispart);
1451: }
1452: return 0;
1453: }
1455: /*@C
1456: PCASMDestroySubdomains - Destroys the index sets created with
1457: `PCASMCreateSubdomains()`. Should be called after setting subdomains with `PCASMSetLocalSubdomains()`.
1459: Collective
1461: Input Parameters:
1462: + n - the number of index sets
1463: . is - the array of index sets
1464: - is_local - the array of local index sets, can be NULL
1466: Level: advanced
1468: .seealso: `PCASM`, `PCASMCreateSubdomains()`, `PCASMSetLocalSubdomains()`
1469: @*/
1470: PetscErrorCode PCASMDestroySubdomains(PetscInt n, IS is[], IS is_local[])
1471: {
1472: PetscInt i;
1474: if (n <= 0) return 0;
1475: if (is) {
1477: for (i = 0; i < n; i++) ISDestroy(&is[i]);
1478: PetscFree(is);
1479: }
1480: if (is_local) {
1482: for (i = 0; i < n; i++) ISDestroy(&is_local[i]);
1483: PetscFree(is_local);
1484: }
1485: return 0;
1486: }
1488: /*@
1489: PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1490: preconditioner, `PCASM`, for a two-dimensional problem on a regular grid.
1492: Not Collective
1494: Input Parameters:
1495: + m - the number of mesh points in the x direction
1496: . n - the number of mesh points in the y direction
1497: . M - the number of subdomains in the x direction
1498: . N - the number of subdomains in the y direction
1499: . dof - degrees of freedom per node
1500: - overlap - overlap in mesh lines
1502: Output Parameters:
1503: + Nsub - the number of subdomains created
1504: . is - array of index sets defining overlapping (if overlap > 0) subdomains
1505: - is_local - array of index sets defining non-overlapping subdomains
1507: Note:
1508: Presently `PCAMSCreateSubdomains2d()` is valid only for sequential
1509: preconditioners. More general related routines are
1510: `PCASMSetTotalSubdomains()` and `PCASMSetLocalSubdomains()`.
1512: Level: advanced
1514: .seealso: `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetLocalSubdomains()`, `PCASMGetSubKSP()`,
1515: `PCASMSetOverlap()`
1516: @*/
1517: PetscErrorCode PCASMCreateSubdomains2D(PetscInt m, PetscInt n, PetscInt M, PetscInt N, PetscInt dof, PetscInt overlap, PetscInt *Nsub, IS **is, IS **is_local)
1518: {
1519: PetscInt i, j, height, width, ystart, xstart, yleft, yright, xleft, xright, loc_outer;
1520: PetscInt nidx, *idx, loc, ii, jj, count;
1524: *Nsub = N * M;
1525: PetscMalloc1(*Nsub, is);
1526: PetscMalloc1(*Nsub, is_local);
1527: ystart = 0;
1528: loc_outer = 0;
1529: for (i = 0; i < N; i++) {
1530: height = n / N + ((n % N) > i); /* height of subdomain */
1532: yleft = ystart - overlap;
1533: if (yleft < 0) yleft = 0;
1534: yright = ystart + height + overlap;
1535: if (yright > n) yright = n;
1536: xstart = 0;
1537: for (j = 0; j < M; j++) {
1538: width = m / M + ((m % M) > j); /* width of subdomain */
1540: xleft = xstart - overlap;
1541: if (xleft < 0) xleft = 0;
1542: xright = xstart + width + overlap;
1543: if (xright > m) xright = m;
1544: nidx = (xright - xleft) * (yright - yleft);
1545: PetscMalloc1(nidx, &idx);
1546: loc = 0;
1547: for (ii = yleft; ii < yright; ii++) {
1548: count = m * ii + xleft;
1549: for (jj = xleft; jj < xright; jj++) idx[loc++] = count++;
1550: }
1551: ISCreateGeneral(PETSC_COMM_SELF, nidx, idx, PETSC_COPY_VALUES, (*is) + loc_outer);
1552: if (overlap == 0) {
1553: PetscObjectReference((PetscObject)(*is)[loc_outer]);
1555: (*is_local)[loc_outer] = (*is)[loc_outer];
1556: } else {
1557: for (loc = 0, ii = ystart; ii < ystart + height; ii++) {
1558: for (jj = xstart; jj < xstart + width; jj++) idx[loc++] = m * ii + jj;
1559: }
1560: ISCreateGeneral(PETSC_COMM_SELF, loc, idx, PETSC_COPY_VALUES, *is_local + loc_outer);
1561: }
1562: PetscFree(idx);
1563: xstart += width;
1564: loc_outer++;
1565: }
1566: ystart += height;
1567: }
1568: for (i = 0; i < *Nsub; i++) ISSort((*is)[i]);
1569: return 0;
1570: }
1572: /*@C
1573: PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
1574: only) for the additive Schwarz preconditioner, `PCASM`.
1576: Not Collective
1578: Input Parameter:
1579: . pc - the preconditioner context
1581: Output Parameters:
1582: + n - if requested, the number of subdomains for this processor (default value = 1)
1583: . is - if requested, the index sets that define the subdomains for this processor
1584: - is_local - if requested, the index sets that define the local part of the subdomains for this processor (can be NULL)
1586: Note:
1587: The `IS` numbering is in the parallel, global numbering of the vector.
1589: Level: advanced
1591: .seealso: `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`, `PCASMGetSubKSP()`,
1592: `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`, `PCASMGetLocalSubmatrices()`
1593: @*/
1594: PetscErrorCode PCASMGetLocalSubdomains(PC pc, PetscInt *n, IS *is[], IS *is_local[])
1595: {
1596: PC_ASM *osm = (PC_ASM *)pc->data;
1597: PetscBool match;
1603: PetscObjectTypeCompare((PetscObject)pc, PCASM, &match);
1605: if (n) *n = osm->n_local_true;
1606: if (is) *is = osm->is;
1607: if (is_local) *is_local = osm->is_local;
1608: return 0;
1609: }
1611: /*@C
1612: PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1613: only) for the additive Schwarz preconditioner, `PCASM`.
1615: Not Collective
1617: Input Parameter:
1618: . pc - the preconditioner context
1620: Output Parameters:
1621: + n - if requested, the number of matrices for this processor (default value = 1)
1622: - mat - if requested, the matrices
1624: Level: advanced
1626: Notes:
1627: Call after `PCSetUp()` (or `KSPSetUp()`) but before `PCApply()` and before `PCSetUpOnBlocks()`)
1629: Usually one would use `PCSetModifySubMatrices()` to change the submatrices in building the preconditioner.
1631: .seealso: `PCASM`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`, `PCASMGetSubKSP()`,
1632: `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`, `PCASMGetLocalSubdomains()`, `PCSetModifySubMatrices()`
1633: @*/
1634: PetscErrorCode PCASMGetLocalSubmatrices(PC pc, PetscInt *n, Mat *mat[])
1635: {
1636: PC_ASM *osm;
1637: PetscBool match;
1643: PetscObjectTypeCompare((PetscObject)pc, PCASM, &match);
1644: if (!match) {
1645: if (n) *n = 0;
1646: if (mat) *mat = NULL;
1647: } else {
1648: osm = (PC_ASM *)pc->data;
1649: if (n) *n = osm->n_local_true;
1650: if (mat) *mat = osm->pmat;
1651: }
1652: return 0;
1653: }
1655: /*@
1656: PCASMSetDMSubdomains - Indicates whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible.
1658: Logically Collective
1660: Input Parameters:
1661: + pc - the preconditioner
1662: - flg - boolean indicating whether to use subdomains defined by the `DM`
1664: Options Database Key:
1665: . -pc_asm_dm_subdomains <bool> - use subdomains defined by the `DM`
1667: Level: intermediate
1669: Note:
1670: `PCASMSetTotalSubdomains()` and `PCASMSetOverlap()` take precedence over `PCASMSetDMSubdomains()`,
1671: so setting either of the first two effectively turns the latter off.
1673: .seealso: `PCASM`, `PCASMGetDMSubdomains()`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`
1674: `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`, `PCASMGetLocalSubdomains()`
1675: @*/
1676: PetscErrorCode PCASMSetDMSubdomains(PC pc, PetscBool flg)
1677: {
1678: PC_ASM *osm = (PC_ASM *)pc->data;
1679: PetscBool match;
1684: PetscObjectTypeCompare((PetscObject)pc, PCASM, &match);
1685: if (match) osm->dm_subdomains = flg;
1686: return 0;
1687: }
1689: /*@
1690: PCASMGetDMSubdomains - Returns flag indicating whether to use `DMCreateDomainDecomposition()` to define the subdomains, whenever possible.
1692: Not Collective
1694: Input Parameter:
1695: . pc - the preconditioner
1697: Output Parameter:
1698: . flg - boolean indicating whether to use subdomains defined by the DM
1700: Level: intermediate
1702: .seealso: `PCASM`, `PCASMSetDMSubdomains()`, `PCASMSetTotalSubdomains()`, `PCASMSetOverlap()`
1703: `PCASMCreateSubdomains2D()`, `PCASMSetLocalSubdomains()`, `PCASMGetLocalSubdomains()`
1704: @*/
1705: PetscErrorCode PCASMGetDMSubdomains(PC pc, PetscBool *flg)
1706: {
1707: PC_ASM *osm = (PC_ASM *)pc->data;
1708: PetscBool match;
1712: PetscObjectTypeCompare((PetscObject)pc, PCASM, &match);
1713: if (match) *flg = osm->dm_subdomains;
1714: else *flg = PETSC_FALSE;
1715: return 0;
1716: }
1718: /*@
1719: PCASMGetSubMatType - Gets the matrix type used for `PCASM` subsolves, as a string.
1721: Not Collective
1723: Input Parameter:
1724: . pc - the `PC`
1726: Output Parameter:
1727: . pc_asm_sub_mat_type - name of matrix type
1729: Level: advanced
1731: .seealso: `PCASM`, `PCASMSetSubMatType()`, `PCASM`, `PCSetType()`, `VecSetType()`, `MatType`, `Mat`
1732: @*/
1733: PetscErrorCode PCASMGetSubMatType(PC pc, MatType *sub_mat_type)
1734: {
1736: PetscTryMethod(pc, "PCASMGetSubMatType_C", (PC, MatType *), (pc, sub_mat_type));
1737: return 0;
1738: }
1740: /*@
1741: PCASMSetSubMatType - Set the type of matrix used for `PCASM` subsolves
1743: Collective
1745: Input Parameters:
1746: + pc - the `PC` object
1747: - sub_mat_type - the `MatType`
1749: Options Database Key:
1750: . -pc_asm_sub_mat_type <sub_mat_type> - Sets the matrix type used for subsolves, for example, seqaijviennacl.
1751: If you specify a base name like aijviennacl, the corresponding sequential type is assumed.
1753: Note:
1754: See "${PETSC_DIR}/include/petscmat.h" for available types
1756: Level: advanced
1758: .seealso: `PCASM`, `PCASMGetSubMatType()`, `PCASM`, `PCSetType()`, `VecSetType()`, `MatType`, `Mat`
1759: @*/
1760: PetscErrorCode PCASMSetSubMatType(PC pc, MatType sub_mat_type)
1761: {
1763: PetscTryMethod(pc, "PCASMSetSubMatType_C", (PC, MatType), (pc, sub_mat_type));
1764: return 0;
1765: }