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>
15: static PetscErrorCode PCView_ASM(PC pc,PetscViewer viewer)
16: {
17: PC_ASM *osm = (PC_ASM*)pc->data;
18: PetscMPIInt rank;
19: PetscInt i,bsz;
20: PetscBool iascii,isstring;
21: PetscViewer sviewer;
22: PetscViewerFormat format;
23: const char *prefix;
25: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
26: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
27: if (iascii) {
28: char overlaps[256] = "user-defined overlap",blocks[256] = "total subdomain blocks not yet set";
29: if (osm->overlap >= 0) PetscSNPrintf(overlaps,sizeof(overlaps),"amount of overlap = %D",osm->overlap);
30: if (osm->n > 0) PetscSNPrintf(blocks,sizeof(blocks),"total subdomain blocks = %D",osm->n);
31: PetscViewerASCIIPrintf(viewer," %s, %s\n",blocks,overlaps);
32: PetscViewerASCIIPrintf(viewer," restriction/interpolation type - %s\n",PCASMTypes[osm->type]);
33: if (osm->dm_subdomains) PetscViewerASCIIPrintf(viewer," Additive Schwarz: using DM to define subdomains\n");
34: if (osm->loctype != PC_COMPOSITE_ADDITIVE) PetscViewerASCIIPrintf(viewer," Additive Schwarz: local solve composition type - %s\n",PCCompositeTypes[osm->loctype]);
35: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
36: PetscViewerGetFormat(viewer,&format);
37: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
38: if (osm->ksp) {
39: PetscViewerASCIIPrintf(viewer," Local solver information for first block is in the following KSP and PC objects on rank 0:\n");
40: PCGetOptionsPrefix(pc,&prefix);
41: PetscViewerASCIIPrintf(viewer," Use -%sksp_view ::ascii_info_detail to display information for all blocks\n",prefix?prefix:"");
42: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
43: if (rank == 0) {
44: PetscViewerASCIIPushTab(viewer);
45: KSPView(osm->ksp[0],sviewer);
46: PetscViewerASCIIPopTab(viewer);
47: }
48: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
49: }
50: } else {
51: PetscViewerASCIIPushSynchronized(viewer);
52: PetscViewerASCIISynchronizedPrintf(viewer," [%d] number of local blocks = %D\n",(int)rank,osm->n_local_true);
53: PetscViewerFlush(viewer);
54: PetscViewerASCIIPrintf(viewer," Local solver information for each block is in the following KSP and PC objects:\n");
55: PetscViewerASCIIPushTab(viewer);
56: PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
57: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
58: for (i=0; i<osm->n_local_true; i++) {
59: ISGetLocalSize(osm->is[i],&bsz);
60: PetscViewerASCIISynchronizedPrintf(sviewer,"[%d] local block number %D, size = %D\n",(int)rank,i,bsz);
61: KSPView(osm->ksp[i],sviewer);
62: PetscViewerASCIISynchronizedPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
63: }
64: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
65: PetscViewerASCIIPopTab(viewer);
66: PetscViewerFlush(viewer);
67: PetscViewerASCIIPopSynchronized(viewer);
68: }
69: } else if (isstring) {
70: PetscViewerStringSPrintf(viewer," blocks=%D, overlap=%D, type=%s",osm->n,osm->overlap,PCASMTypes[osm->type]);
71: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
72: if (osm->ksp) KSPView(osm->ksp[0],sviewer);
73: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
74: }
75: return 0;
76: }
78: static PetscErrorCode PCASMPrintSubdomains(PC pc)
79: {
80: PC_ASM *osm = (PC_ASM*)pc->data;
81: const char *prefix;
82: char fname[PETSC_MAX_PATH_LEN+1];
83: PetscViewer viewer, sviewer;
84: char *s;
85: PetscInt i,j,nidx;
86: const PetscInt *idx;
87: PetscMPIInt rank, size;
89: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
90: MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
91: PCGetOptionsPrefix(pc,&prefix);
92: PetscOptionsGetString(NULL,prefix,"-pc_asm_print_subdomains",fname,sizeof(fname),NULL);
93: if (fname[0] == 0) PetscStrcpy(fname,"stdout");
94: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),fname,&viewer);
95: for (i=0; i<osm->n_local; i++) {
96: if (i < osm->n_local_true) {
97: ISGetLocalSize(osm->is[i],&nidx);
98: ISGetIndices(osm->is[i],&idx);
99: /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
100: #define len 16*(nidx+1)+512
101: PetscMalloc1(len,&s);
102: PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer);
103: #undef len
104: PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D with overlap:\n", rank, size, i);
105: for (j=0; j<nidx; j++) {
106: PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
107: }
108: ISRestoreIndices(osm->is[i],&idx);
109: PetscViewerStringSPrintf(sviewer,"\n");
110: PetscViewerDestroy(&sviewer);
111: PetscViewerASCIIPushSynchronized(viewer);
112: PetscViewerASCIISynchronizedPrintf(viewer, s);
113: PetscViewerFlush(viewer);
114: PetscViewerASCIIPopSynchronized(viewer);
115: 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: PetscMalloc1(len, &s);
120: PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer);
121: #undef len
122: PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D without overlap:\n", rank, size, i);
123: ISGetLocalSize(osm->is_local[i],&nidx);
124: ISGetIndices(osm->is_local[i],&idx);
125: for (j=0; j<nidx; j++) {
126: PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
127: }
128: ISRestoreIndices(osm->is_local[i],&idx);
129: PetscViewerStringSPrintf(sviewer,"\n");
130: PetscViewerDestroy(&sviewer);
131: PetscViewerASCIIPushSynchronized(viewer);
132: PetscViewerASCIISynchronizedPrintf(viewer, s);
133: PetscViewerFlush(viewer);
134: PetscViewerASCIIPopSynchronized(viewer);
135: PetscFree(s);
136: }
137: } else {
138: /* Participate in collective viewer calls. */
139: PetscViewerASCIIPushSynchronized(viewer);
140: PetscViewerFlush(viewer);
141: PetscViewerASCIIPopSynchronized(viewer);
142: /* Assume either all ranks have is_local or none do. */
143: if (osm->is_local) {
144: PetscViewerASCIIPushSynchronized(viewer);
145: PetscViewerFlush(viewer);
146: PetscViewerASCIIPopSynchronized(viewer);
147: }
148: }
149: }
150: PetscViewerFlush(viewer);
151: PetscViewerDestroy(&viewer);
152: return 0;
153: }
155: static PetscErrorCode PCSetUp_ASM(PC pc)
156: {
157: PC_ASM *osm = (PC_ASM*)pc->data;
158: PetscBool flg;
159: PetscInt i,m,m_local;
160: MatReuse scall = MAT_REUSE_MATRIX;
161: IS isl;
162: KSP ksp;
163: PC subpc;
164: const char *prefix,*pprefix;
165: Vec vec;
166: DM *domain_dm = NULL;
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: 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) {
185: PCASMSetLocalSubdomains(pc, num_domains, outer_domain_is, inner_domain_is);
186: }
187: for (d = 0; d < num_domains; ++d) {
188: if (domain_names) PetscFree(domain_names[d]);
189: if (inner_domain_is) ISDestroy(&inner_domain_is[d]);
190: if (outer_domain_is) ISDestroy(&outer_domain_is[d]);
191: }
192: PetscFree(domain_names);
193: PetscFree(inner_domain_is);
194: PetscFree(outer_domain_is);
195: }
196: if (osm->n_local_true == PETSC_DECIDE) {
197: /* still no subdomains; use one subdomain per processor */
198: osm->n_local_true = 1;
199: }
200: }
201: { /* determine the global and max number of subdomains */
202: struct {PetscInt max,sum;} inwork,outwork;
203: PetscMPIInt size;
205: inwork.max = osm->n_local_true;
206: inwork.sum = osm->n_local_true;
207: 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: 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: MatSetOption(pc->pmat,MAT_SUBMAT_SINGLEIS,PETSC_TRUE);
215: }
216: }
217: if (!osm->is) { /* create the index sets */
218: PCASMCreateSubdomains(pc->pmat,osm->n_local_true,&osm->is);
219: }
220: if (osm->n_local_true > 1 && !osm->is_local) {
221: 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: ISDuplicate(osm->is[i],&osm->is_local[i]);
225: ISCopy(osm->is[i],osm->is_local[i]);
226: } else {
227: PetscObjectReference((PetscObject)osm->is[i]);
228: osm->is_local[i] = osm->is[i];
229: }
230: }
231: }
232: PCGetOptionsPrefix(pc,&prefix);
233: flg = PETSC_FALSE;
234: PetscOptionsHasName(NULL,prefix,"-pc_asm_print_subdomains",&flg);
235: if (flg) PCASMPrintSubdomains(pc);
237: if (osm->overlap > 0) {
238: /* Extend the "overlapping" regions by a number of steps */
239: MatIncreaseOverlap(pc->pmat,osm->n_local_true,osm->is,osm->overlap);
240: }
241: if (osm->sort_indices) {
242: for (i=0; i<osm->n_local_true; i++) {
243: ISSort(osm->is[i]);
244: if (osm->is_local) {
245: ISSort(osm->is_local[i]);
246: }
247: }
248: }
250: if (!osm->ksp) {
251: /* Create the local solvers */
252: PetscMalloc1(osm->n_local_true,&osm->ksp);
253: if (domain_dm) {
254: PetscInfo(pc,"Setting up ASM subproblems using the embedded DM\n");
255: }
256: for (i=0; i<osm->n_local_true; i++) {
257: KSPCreate(PETSC_COMM_SELF,&ksp);
258: KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
259: PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
260: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
261: KSPSetType(ksp,KSPPREONLY);
262: KSPGetPC(ksp,&subpc);
263: PCGetOptionsPrefix(pc,&prefix);
264: KSPSetOptionsPrefix(ksp,prefix);
265: KSPAppendOptionsPrefix(ksp,"sub_");
266: if (domain_dm) {
267: KSPSetDM(ksp, domain_dm[i]);
268: KSPSetDMActive(ksp, PETSC_FALSE);
269: DMDestroy(&domain_dm[i]);
270: }
271: osm->ksp[i] = ksp;
272: }
273: if (domain_dm) {
274: PetscFree(domain_dm);
275: }
276: }
278: ISConcatenate(PETSC_COMM_SELF, osm->n_local_true, osm->is, &osm->lis);
279: ISSortRemoveDups(osm->lis);
280: ISGetLocalSize(osm->lis, &m);
282: scall = MAT_INITIAL_MATRIX;
283: } else {
284: /*
285: Destroy the blocks from the previous iteration
286: */
287: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
288: MatDestroyMatrices(osm->n_local_true,&osm->pmat);
289: scall = MAT_INITIAL_MATRIX;
290: }
291: }
293: /* Destroy previous submatrices of a different type than pc->pmat since MAT_REUSE_MATRIX won't work in that case */
294: if ((scall == MAT_REUSE_MATRIX) && osm->sub_mat_type) {
295: if (osm->n_local_true > 0) {
296: MatDestroySubMatrices(osm->n_local_true,&osm->pmat);
297: }
298: scall = MAT_INITIAL_MATRIX;
299: }
301: /*
302: Extract out the submatrices
303: */
304: MatCreateSubMatrices(pc->pmat,osm->n_local_true,osm->is,osm->is,scall,&osm->pmat);
305: if (scall == MAT_INITIAL_MATRIX) {
306: PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
307: for (i=0; i<osm->n_local_true; i++) {
308: PetscLogObjectParent((PetscObject)pc,(PetscObject)osm->pmat[i]);
309: PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);
310: }
311: }
313: /* Convert the types of the submatrices (if needbe) */
314: if (osm->sub_mat_type) {
315: for (i=0; i<osm->n_local_true; i++) {
316: MatConvert(osm->pmat[i],osm->sub_mat_type,MAT_INPLACE_MATRIX,&(osm->pmat[i]));
317: }
318: }
320: if (!pc->setupcalled) {
321: VecType vtype;
323: /* Create the local work vectors (from the local matrices) and scatter contexts */
324: MatCreateVecs(pc->pmat,&vec,NULL);
327: if (osm->is_local && osm->type == PC_ASM_RESTRICT && osm->loctype == PC_COMPOSITE_ADDITIVE) {
328: PetscMalloc1(osm->n_local_true,&osm->lprolongation);
329: }
330: PetscMalloc1(osm->n_local_true,&osm->lrestriction);
331: PetscMalloc1(osm->n_local_true,&osm->x);
332: PetscMalloc1(osm->n_local_true,&osm->y);
334: ISGetLocalSize(osm->lis,&m);
335: ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
336: MatGetVecType(osm->pmat[0],&vtype);
337: VecCreate(PETSC_COMM_SELF,&osm->lx);
338: VecSetSizes(osm->lx,m,m);
339: VecSetType(osm->lx,vtype);
340: VecDuplicate(osm->lx, &osm->ly);
341: VecScatterCreate(vec,osm->lis,osm->lx,isl,&osm->restriction);
342: ISDestroy(&isl);
344: for (i=0; i<osm->n_local_true; ++i) {
345: ISLocalToGlobalMapping ltog;
346: IS isll;
347: const PetscInt *idx_is;
348: PetscInt *idx_lis,nout;
350: ISGetLocalSize(osm->is[i],&m);
351: MatCreateVecs(osm->pmat[i],&osm->x[i],NULL);
352: VecDuplicate(osm->x[i],&osm->y[i]);
354: /* generate a scatter from ly to y[i] picking all the overlapping is[i] entries */
355: ISLocalToGlobalMappingCreateIS(osm->lis,<og);
356: ISGetLocalSize(osm->is[i],&m);
357: ISGetIndices(osm->is[i], &idx_is);
358: PetscMalloc1(m,&idx_lis);
359: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m,idx_is,&nout,idx_lis);
361: ISRestoreIndices(osm->is[i], &idx_is);
362: ISCreateGeneral(PETSC_COMM_SELF,m,idx_lis,PETSC_OWN_POINTER,&isll);
363: ISLocalToGlobalMappingDestroy(<og);
364: ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
365: VecScatterCreate(osm->ly,isll,osm->y[i],isl,&osm->lrestriction[i]);
366: ISDestroy(&isll);
367: ISDestroy(&isl);
368: if (osm->lprolongation) { /* generate a scatter from y[i] to ly picking only the the non-overlapping is_local[i] entries */
369: ISLocalToGlobalMapping ltog;
370: IS isll,isll_local;
371: const PetscInt *idx_local;
372: PetscInt *idx1, *idx2, nout;
374: ISGetLocalSize(osm->is_local[i],&m_local);
375: ISGetIndices(osm->is_local[i], &idx_local);
377: ISLocalToGlobalMappingCreateIS(osm->is[i],<og);
378: PetscMalloc1(m_local,&idx1);
379: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx1);
380: ISLocalToGlobalMappingDestroy(<og);
382: ISCreateGeneral(PETSC_COMM_SELF,m_local,idx1,PETSC_OWN_POINTER,&isll);
384: ISLocalToGlobalMappingCreateIS(osm->lis,<og);
385: PetscMalloc1(m_local,&idx2);
386: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx2);
387: ISLocalToGlobalMappingDestroy(<og);
389: ISCreateGeneral(PETSC_COMM_SELF,m_local,idx2,PETSC_OWN_POINTER,&isll_local);
391: ISRestoreIndices(osm->is_local[i], &idx_local);
392: VecScatterCreate(osm->y[i],isll,osm->ly,isll_local,&osm->lprolongation[i]);
394: ISDestroy(&isll);
395: ISDestroy(&isll_local);
396: }
397: }
398: VecDestroy(&vec);
399: }
401: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
402: IS *cis;
403: PetscInt c;
405: PetscMalloc1(osm->n_local_true, &cis);
406: for (c = 0; c < osm->n_local_true; ++c) cis[c] = osm->lis;
407: MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, cis, scall, &osm->lmats);
408: PetscFree(cis);
409: }
411: /* Return control to the user so that the submatrices can be modified (e.g., to apply
412: different boundary conditions for the submatrices than for the global problem) */
413: PCModifySubMatrices(pc,osm->n_local_true,osm->is,osm->is,osm->pmat,pc->modifysubmatricesP);
415: /*
416: Loop over subdomains putting them into local ksp
417: */
418: KSPGetOptionsPrefix(osm->ksp[0],&prefix);
419: for (i=0; i<osm->n_local_true; i++) {
420: KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i]);
421: MatSetOptionsPrefix(osm->pmat[i],prefix);
422: if (!pc->setupcalled) {
423: KSPSetFromOptions(osm->ksp[i]);
424: }
425: }
426: return 0;
427: }
429: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
430: {
431: PC_ASM *osm = (PC_ASM*)pc->data;
432: PetscInt i;
433: KSPConvergedReason reason;
435: for (i=0; i<osm->n_local_true; i++) {
436: KSPSetUp(osm->ksp[i]);
437: KSPGetConvergedReason(osm->ksp[i],&reason);
438: if (reason == KSP_DIVERGED_PC_FAILED) {
439: pc->failedreason = PC_SUBPC_ERROR;
440: }
441: }
442: return 0;
443: }
445: static PetscErrorCode PCApply_ASM(PC pc,Vec x,Vec y)
446: {
447: PC_ASM *osm = (PC_ASM*)pc->data;
448: PetscInt i,n_local_true = osm->n_local_true;
449: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
451: /*
452: support for limiting the restriction or interpolation to only local
453: subdomain values (leaving the other values 0).
454: */
455: if (!(osm->type & PC_ASM_RESTRICT)) {
456: forward = SCATTER_FORWARD_LOCAL;
457: /* have to zero the work RHS since scatter may leave some slots empty */
458: VecSet(osm->lx, 0.0);
459: }
460: if (!(osm->type & PC_ASM_INTERPOLATE)) {
461: reverse = SCATTER_REVERSE_LOCAL;
462: }
464: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE) {
465: /* zero the global and the local solutions */
466: VecSet(y, 0.0);
467: VecSet(osm->ly, 0.0);
469: /* copy the global RHS to local RHS including the ghost nodes */
470: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
471: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
473: /* restrict local RHS to the overlapping 0-block RHS */
474: VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
475: VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
477: /* do the local solves */
478: for (i = 0; i < n_local_true; ++i) {
480: /* solve the overlapping i-block */
481: PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i],0);
482: KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]);
483: KSPCheckSolve(osm->ksp[i], pc, osm->y[i]);
484: PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0);
486: if (osm->lprolongation) { /* interpolate the non-overlapping i-block solution to the local solution (only for restrictive additive) */
487: VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
488: VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
489: } else { /* interpolate the overlapping i-block solution to the local solution */
490: VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
491: VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
492: }
494: if (i < n_local_true-1) {
495: /* restrict local RHS to the overlapping (i+1)-block RHS */
496: VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
497: VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
499: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
500: /* update the overlapping (i+1)-block RHS using the current local solution */
501: MatMult(osm->lmats[i+1], osm->ly, osm->y[i+1]);
502: VecAXPBY(osm->x[i+1],-1.,1., osm->y[i+1]);
503: }
504: }
505: }
506: /* add the local solution to the global solution including the ghost nodes */
507: VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
508: VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
509: } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
510: return 0;
511: }
513: static PetscErrorCode PCMatApply_ASM(PC pc,Mat X,Mat Y)
514: {
515: PC_ASM *osm = (PC_ASM*)pc->data;
516: Mat Z,W;
517: Vec x;
518: PetscInt i,m,N;
519: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
522: /*
523: support for limiting the restriction or interpolation to only local
524: subdomain values (leaving the other values 0).
525: */
526: if (!(osm->type & PC_ASM_RESTRICT)) {
527: forward = SCATTER_FORWARD_LOCAL;
528: /* have to zero the work RHS since scatter may leave some slots empty */
529: VecSet(osm->lx, 0.0);
530: }
531: if (!(osm->type & PC_ASM_INTERPOLATE)) {
532: reverse = SCATTER_REVERSE_LOCAL;
533: }
534: VecGetLocalSize(osm->x[0], &m);
535: MatGetSize(X, NULL, &N);
536: MatCreateSeqDense(PETSC_COMM_SELF, m, N, NULL, &Z);
537: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE) {
538: /* zero the global and the local solutions */
539: MatZeroEntries(Y);
540: VecSet(osm->ly, 0.0);
542: for (i = 0; i < N; ++i) {
543: MatDenseGetColumnVecRead(X, i, &x);
544: /* copy the global RHS to local RHS including the ghost nodes */
545: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
546: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
547: MatDenseRestoreColumnVecRead(X, i, &x);
549: MatDenseGetColumnVecWrite(Z, i, &x);
550: /* restrict local RHS to the overlapping 0-block RHS */
551: VecScatterBegin(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward);
552: VecScatterEnd(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward);
553: MatDenseRestoreColumnVecWrite(Z, i, &x);
554: }
555: MatCreateSeqDense(PETSC_COMM_SELF, m, N, NULL, &W);
556: /* solve the overlapping 0-block */
557: PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[0], Z, W, 0);
558: KSPMatSolve(osm->ksp[0], Z, W);
559: KSPCheckSolve(osm->ksp[0], pc, NULL);
560: PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[0], Z, W,0);
561: MatDestroy(&Z);
563: for (i = 0; i < N; ++i) {
564: VecSet(osm->ly, 0.0);
565: MatDenseGetColumnVecRead(W, i, &x);
566: if (osm->lprolongation) { /* interpolate the non-overlapping 0-block solution to the local solution (only for restrictive additive) */
567: VecScatterBegin(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward);
568: VecScatterEnd(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward);
569: } else { /* interpolate the overlapping 0-block solution to the local solution */
570: VecScatterBegin(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse);
571: VecScatterEnd(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse);
572: }
573: MatDenseRestoreColumnVecRead(W, i, &x);
575: MatDenseGetColumnVecWrite(Y, i, &x);
576: /* add the local solution to the global solution including the ghost nodes */
577: VecScatterBegin(osm->restriction, osm->ly, x, ADD_VALUES, reverse);
578: VecScatterEnd(osm->restriction, osm->ly, x, ADD_VALUES, reverse);
579: MatDenseRestoreColumnVecWrite(Y, i, &x);
580: }
581: MatDestroy(&W);
582: } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
583: return 0;
584: }
586: static PetscErrorCode PCApplyTranspose_ASM(PC pc,Vec x,Vec y)
587: {
588: PC_ASM *osm = (PC_ASM*)pc->data;
589: PetscInt i,n_local_true = osm->n_local_true;
590: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
592: /*
593: Support for limiting the restriction or interpolation to only local
594: subdomain values (leaving the other values 0).
596: Note: these are reversed from the PCApply_ASM() because we are applying the
597: transpose of the three terms
598: */
600: if (!(osm->type & PC_ASM_INTERPOLATE)) {
601: forward = SCATTER_FORWARD_LOCAL;
602: /* have to zero the work RHS since scatter may leave some slots empty */
603: VecSet(osm->lx, 0.0);
604: }
605: if (!(osm->type & PC_ASM_RESTRICT)) reverse = SCATTER_REVERSE_LOCAL;
607: /* zero the global and the local solutions */
608: VecSet(y, 0.0);
609: VecSet(osm->ly, 0.0);
611: /* Copy the global RHS to local RHS including the ghost nodes */
612: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
613: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
615: /* Restrict local RHS to the overlapping 0-block RHS */
616: VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
617: VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
619: /* do the local solves */
620: for (i = 0; i < n_local_true; ++i) {
622: /* solve the overlapping i-block */
623: PetscLogEventBegin(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
624: KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i]);
625: KSPCheckSolve(osm->ksp[i],pc,osm->y[i]);
626: PetscLogEventEnd(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
628: if (osm->lprolongation) { /* interpolate the non-overlapping i-block solution to the local solution */
629: VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
630: 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: VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
633: 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: VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
639: 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: VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
644: VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
646: return 0;
648: }
650: static PetscErrorCode PCReset_ASM(PC pc)
651: {
652: PC_ASM *osm = (PC_ASM*)pc->data;
653: PetscInt i;
655: if (osm->ksp) {
656: for (i=0; i<osm->n_local_true; i++) {
657: KSPReset(osm->ksp[i]);
658: }
659: }
660: if (osm->pmat) {
661: if (osm->n_local_true > 0) {
662: MatDestroySubMatrices(osm->n_local_true,&osm->pmat);
663: }
664: }
665: if (osm->lrestriction) {
666: VecScatterDestroy(&osm->restriction);
667: for (i=0; i<osm->n_local_true; i++) {
668: VecScatterDestroy(&osm->lrestriction[i]);
669: if (osm->lprolongation) VecScatterDestroy(&osm->lprolongation[i]);
670: VecDestroy(&osm->x[i]);
671: VecDestroy(&osm->y[i]);
672: }
673: PetscFree(osm->lrestriction);
674: if (osm->lprolongation) PetscFree(osm->lprolongation);
675: PetscFree(osm->x);
676: PetscFree(osm->y);
678: }
679: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
680: ISDestroy(&osm->lis);
681: VecDestroy(&osm->lx);
682: VecDestroy(&osm->ly);
683: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
684: MatDestroyMatrices(osm->n_local_true, &osm->lmats);
685: }
687: PetscFree(osm->sub_mat_type);
689: osm->is = NULL;
690: osm->is_local = NULL;
691: return 0;
692: }
694: static PetscErrorCode PCDestroy_ASM(PC pc)
695: {
696: PC_ASM *osm = (PC_ASM*)pc->data;
697: PetscInt i;
699: PCReset_ASM(pc);
700: if (osm->ksp) {
701: for (i=0; i<osm->n_local_true; i++) {
702: KSPDestroy(&osm->ksp[i]);
703: }
704: PetscFree(osm->ksp);
705: }
706: PetscFree(pc->data);
708: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",NULL);
709: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",NULL);
710: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetOverlap_C",NULL);
711: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetType_C",NULL);
712: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetType_C",NULL);
713: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalType_C",NULL);
714: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetLocalType_C",NULL);
715: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSortIndices_C",NULL);
716: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubKSP_C",NULL);
717: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubMatType_C",NULL);
718: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSubMatType_C",NULL);
719: return 0;
720: }
722: static PetscErrorCode PCSetFromOptions_ASM(PetscOptionItems *PetscOptionsObject,PC pc)
723: {
724: PC_ASM *osm = (PC_ASM*)pc->data;
725: PetscInt blocks,ovl;
726: PetscBool flg;
727: PCASMType asmtype;
728: PCCompositeType loctype;
729: char sub_mat_type[256];
731: PetscOptionsHead(PetscOptionsObject,"Additive Schwarz options");
732: PetscOptionsBool("-pc_asm_dm_subdomains","Use DMCreateDomainDecomposition() to define subdomains","PCASMSetDMSubdomains",osm->dm_subdomains,&osm->dm_subdomains,&flg);
733: PetscOptionsInt("-pc_asm_blocks","Number of subdomains","PCASMSetTotalSubdomains",osm->n,&blocks,&flg);
734: if (flg) {
735: PCASMSetTotalSubdomains(pc,blocks,NULL,NULL);
736: osm->dm_subdomains = PETSC_FALSE;
737: }
738: PetscOptionsInt("-pc_asm_local_blocks","Number of local subdomains","PCASMSetLocalSubdomains",osm->n_local_true,&blocks,&flg);
739: if (flg) {
740: PCASMSetLocalSubdomains(pc,blocks,NULL,NULL);
741: osm->dm_subdomains = PETSC_FALSE;
742: }
743: PetscOptionsInt("-pc_asm_overlap","Number of grid points overlap","PCASMSetOverlap",osm->overlap,&ovl,&flg);
744: if (flg) {
745: PCASMSetOverlap(pc,ovl);
746: osm->dm_subdomains = PETSC_FALSE;
747: }
748: flg = PETSC_FALSE;
749: PetscOptionsEnum("-pc_asm_type","Type of restriction/extension","PCASMSetType",PCASMTypes,(PetscEnum)osm->type,(PetscEnum*)&asmtype,&flg);
750: if (flg) PCASMSetType(pc,asmtype);
751: flg = PETSC_FALSE;
752: PetscOptionsEnum("-pc_asm_local_type","Type of local solver composition","PCASMSetLocalType",PCCompositeTypes,(PetscEnum)osm->loctype,(PetscEnum*)&loctype,&flg);
753: if (flg) PCASMSetLocalType(pc,loctype);
754: PetscOptionsFList("-pc_asm_sub_mat_type","Subsolve Matrix Type","PCASMSetSubMatType",MatList,NULL,sub_mat_type,256,&flg);
755: if (flg) {
756: PCASMSetSubMatType(pc,sub_mat_type);
757: }
758: PetscOptionsTail();
759: return 0;
760: }
762: /*------------------------------------------------------------------------------------*/
764: static PetscErrorCode PCASMSetLocalSubdomains_ASM(PC pc,PetscInt n,IS is[],IS is_local[])
765: {
766: PC_ASM *osm = (PC_ASM*)pc->data;
767: PetscInt i;
772: if (!pc->setupcalled) {
773: if (is) {
774: for (i=0; i<n; i++) PetscObjectReference((PetscObject)is[i]);
775: }
776: if (is_local) {
777: for (i=0; i<n; i++) PetscObjectReference((PetscObject)is_local[i]);
778: }
779: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
781: osm->n_local_true = n;
782: osm->is = NULL;
783: osm->is_local = NULL;
784: if (is) {
785: PetscMalloc1(n,&osm->is);
786: for (i=0; i<n; i++) osm->is[i] = is[i];
787: /* Flag indicating that the user has set overlapping subdomains so PCASM should not increase their size. */
788: osm->overlap = -1;
789: }
790: if (is_local) {
791: PetscMalloc1(n,&osm->is_local);
792: for (i=0; i<n; i++) osm->is_local[i] = is_local[i];
793: if (!is) {
794: PetscMalloc1(osm->n_local_true,&osm->is);
795: for (i=0; i<osm->n_local_true; i++) {
796: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
797: ISDuplicate(osm->is_local[i],&osm->is[i]);
798: ISCopy(osm->is_local[i],osm->is[i]);
799: } else {
800: PetscObjectReference((PetscObject)osm->is_local[i]);
801: osm->is[i] = osm->is_local[i];
802: }
803: }
804: }
805: }
806: }
807: return 0;
808: }
810: static PetscErrorCode PCASMSetTotalSubdomains_ASM(PC pc,PetscInt N,IS *is,IS *is_local)
811: {
812: PC_ASM *osm = (PC_ASM*)pc->data;
813: PetscMPIInt rank,size;
814: PetscInt n;
819: /*
820: Split the subdomains equally among all processors
821: */
822: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
823: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
824: n = N/size + ((N % size) > rank);
827: if (!pc->setupcalled) {
828: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
830: osm->n_local_true = n;
831: osm->is = NULL;
832: osm->is_local = NULL;
833: }
834: return 0;
835: }
837: static PetscErrorCode PCASMSetOverlap_ASM(PC pc,PetscInt ovl)
838: {
839: PC_ASM *osm = (PC_ASM*)pc->data;
843: if (!pc->setupcalled) osm->overlap = ovl;
844: return 0;
845: }
847: static PetscErrorCode PCASMSetType_ASM(PC pc,PCASMType type)
848: {
849: PC_ASM *osm = (PC_ASM*)pc->data;
851: osm->type = type;
852: osm->type_set = PETSC_TRUE;
853: return 0;
854: }
856: static PetscErrorCode PCASMGetType_ASM(PC pc,PCASMType *type)
857: {
858: PC_ASM *osm = (PC_ASM*)pc->data;
860: *type = osm->type;
861: return 0;
862: }
864: static PetscErrorCode PCASMSetLocalType_ASM(PC pc, PCCompositeType type)
865: {
866: PC_ASM *osm = (PC_ASM *) pc->data;
869: osm->loctype = type;
870: return 0;
871: }
873: static PetscErrorCode PCASMGetLocalType_ASM(PC pc, PCCompositeType *type)
874: {
875: PC_ASM *osm = (PC_ASM *) pc->data;
877: *type = osm->loctype;
878: return 0;
879: }
881: static PetscErrorCode PCASMSetSortIndices_ASM(PC pc,PetscBool doSort)
882: {
883: PC_ASM *osm = (PC_ASM*)pc->data;
885: osm->sort_indices = doSort;
886: return 0;
887: }
889: static PetscErrorCode PCASMGetSubKSP_ASM(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
890: {
891: PC_ASM *osm = (PC_ASM*)pc->data;
895: if (n_local) *n_local = osm->n_local_true;
896: if (first_local) {
897: MPI_Scan(&osm->n_local_true,first_local,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
898: *first_local -= osm->n_local_true;
899: }
900: if (ksp) *ksp = osm->ksp;
901: return 0;
902: }
904: static PetscErrorCode PCASMGetSubMatType_ASM(PC pc,MatType *sub_mat_type)
905: {
906: PC_ASM *osm = (PC_ASM*)pc->data;
910: *sub_mat_type = osm->sub_mat_type;
911: return 0;
912: }
914: static PetscErrorCode PCASMSetSubMatType_ASM(PC pc,MatType sub_mat_type)
915: {
916: PC_ASM *osm = (PC_ASM*)pc->data;
919: PetscFree(osm->sub_mat_type);
920: PetscStrallocpy(sub_mat_type,(char**)&osm->sub_mat_type);
921: return 0;
922: }
924: /*@C
925: PCASMSetLocalSubdomains - Sets the local subdomains (for this processor only) for the additive Schwarz preconditioner.
927: Collective on pc
929: Input Parameters:
930: + pc - the preconditioner context
931: . n - the number of subdomains for this processor (default value = 1)
932: . is - the index set that defines the subdomains for this processor
933: (or NULL for PETSc to determine subdomains)
934: - is_local - the index sets that define the local part of the subdomains for this processor, not used unless PCASMType is PC_ASM_RESTRICT
935: (or NULL to not provide these)
937: Options Database Key:
938: To set the total number of subdomain blocks rather than specify the
939: index sets, use the option
940: . -pc_asm_local_blocks <blks> - Sets local blocks
942: Notes:
943: The IS numbering is in the parallel, global numbering of the vector for both is and is_local
945: By default the ASM preconditioner uses 1 block per processor.
947: Use PCASMSetTotalSubdomains() to set the subdomains for all processors.
949: If is_local is provided and PCASMType is PC_ASM_RESTRICT then the solution only over the is_local region is interpolated
950: back to form the global solution (this is the standard restricted additive Schwarz method)
952: If the is_local is provided and PCASMType is PC_ASM_INTERPOLATE or PC_ASM_NONE then an error is generated since there is
953: no code to handle that case.
955: Level: advanced
957: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
958: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), PCASMType, PCASMSetType()
959: @*/
960: PetscErrorCode PCASMSetLocalSubdomains(PC pc,PetscInt n,IS is[],IS is_local[])
961: {
963: PetscTryMethod(pc,"PCASMSetLocalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,is,is_local));
964: return 0;
965: }
967: /*@C
968: PCASMSetTotalSubdomains - Sets the subdomains for all processors for the
969: additive Schwarz preconditioner. Either all or no processors in the
970: PC communicator must call this routine, with the same index sets.
972: Collective on pc
974: Input Parameters:
975: + pc - the preconditioner context
976: . N - the number of subdomains for all processors
977: . is - the index sets that define the subdomains for all processors
978: (or NULL to ask PETSc to determine the subdomains)
979: - is_local - the index sets that define the local part of the subdomains for this processor
980: (or NULL to not provide this information)
982: Options Database Key:
983: To set the total number of subdomain blocks rather than specify the
984: index sets, use the option
985: . -pc_asm_blocks <blks> - Sets total blocks
987: Notes:
988: Currently you cannot use this to set the actual subdomains with the argument is or is_local.
990: By default the ASM preconditioner uses 1 block per processor.
992: These index sets cannot be destroyed until after completion of the
993: linear solves for which the ASM preconditioner is being used.
995: Use PCASMSetLocalSubdomains() to set local subdomains.
997: The IS numbering is in the parallel, global numbering of the vector for both is and is_local
999: Level: advanced
1001: .seealso: PCASMSetLocalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1002: PCASMCreateSubdomains2D()
1003: @*/
1004: PetscErrorCode PCASMSetTotalSubdomains(PC pc,PetscInt N,IS is[],IS is_local[])
1005: {
1007: PetscTryMethod(pc,"PCASMSetTotalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,N,is,is_local));
1008: return 0;
1009: }
1011: /*@
1012: PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
1013: additive Schwarz preconditioner. Either all or no processors in the
1014: PC communicator must call this routine.
1016: Logically Collective on pc
1018: Input Parameters:
1019: + pc - the preconditioner context
1020: - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1)
1022: Options Database Key:
1023: . -pc_asm_overlap <ovl> - Sets overlap
1025: Notes:
1026: By default the ASM preconditioner uses 1 block per processor. To use
1027: multiple blocks per perocessor, see PCASMSetTotalSubdomains() and
1028: PCASMSetLocalSubdomains() (and the option -pc_asm_blocks <blks>).
1030: The overlap defaults to 1, so if one desires that no additional
1031: overlap be computed beyond what may have been set with a call to
1032: PCASMSetTotalSubdomains() or PCASMSetLocalSubdomains(), then ovl
1033: must be set to be 0. In particular, if one does not explicitly set
1034: the subdomains an application code, then all overlap would be computed
1035: internally by PETSc, and using an overlap of 0 would result in an ASM
1036: variant that is equivalent to the block Jacobi preconditioner.
1038: The default algorithm used by PETSc to increase overlap is fast, but not scalable,
1039: use the option -mat_increase_overlap_scalable when the problem and number of processes is large.
1041: Note that one can define initial index sets with any overlap via
1042: PCASMSetLocalSubdomains(); the routine
1043: PCASMSetOverlap() merely allows PETSc to extend that overlap further
1044: if desired.
1046: Level: intermediate
1048: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1049: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), MatIncreaseOverlap()
1050: @*/
1051: PetscErrorCode PCASMSetOverlap(PC pc,PetscInt ovl)
1052: {
1055: PetscTryMethod(pc,"PCASMSetOverlap_C",(PC,PetscInt),(pc,ovl));
1056: return 0;
1057: }
1059: /*@
1060: PCASMSetType - Sets the type of restriction and interpolation used
1061: for local problems in the additive Schwarz method.
1063: Logically Collective on pc
1065: Input Parameters:
1066: + pc - the preconditioner context
1067: - type - variant of ASM, one of
1068: .vb
1069: PC_ASM_BASIC - full interpolation and restriction
1070: PC_ASM_RESTRICT - full restriction, local processor interpolation
1071: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1072: PC_ASM_NONE - local processor restriction and interpolation
1073: .ve
1075: Options Database Key:
1076: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
1078: Notes:
1079: if the is_local arguments are passed to PCASMSetLocalSubdomains() then they are used when PC_ASM_RESTRICT has been selected
1080: to limit the local processor interpolation
1082: Level: intermediate
1084: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1085: PCASMCreateSubdomains2D(), PCASMType, PCASMSetLocalType(), PCASMGetLocalType()
1086: @*/
1087: PetscErrorCode PCASMSetType(PC pc,PCASMType type)
1088: {
1091: PetscTryMethod(pc,"PCASMSetType_C",(PC,PCASMType),(pc,type));
1092: return 0;
1093: }
1095: /*@
1096: PCASMGetType - Gets the type of restriction and interpolation used
1097: for local problems in the additive Schwarz method.
1099: Logically Collective on pc
1101: Input Parameter:
1102: . pc - the preconditioner context
1104: Output Parameter:
1105: . type - variant of ASM, one of
1107: .vb
1108: PC_ASM_BASIC - full interpolation and restriction
1109: PC_ASM_RESTRICT - full restriction, local processor interpolation
1110: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1111: PC_ASM_NONE - local processor restriction and interpolation
1112: .ve
1114: Options Database Key:
1115: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
1117: Level: intermediate
1119: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1120: PCASMCreateSubdomains2D(), PCASMType, PCASMSetType(), PCASMSetLocalType(), PCASMGetLocalType()
1121: @*/
1122: PetscErrorCode PCASMGetType(PC pc,PCASMType *type)
1123: {
1125: PetscUseMethod(pc,"PCASMGetType_C",(PC,PCASMType*),(pc,type));
1126: return 0;
1127: }
1129: /*@
1130: PCASMSetLocalType - Sets the type of composition used for local problems in the additive Schwarz method.
1132: Logically Collective on pc
1134: Input Parameters:
1135: + pc - the preconditioner context
1136: - type - type of composition, one of
1137: .vb
1138: PC_COMPOSITE_ADDITIVE - local additive combination
1139: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1140: .ve
1142: Options Database Key:
1143: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1145: Level: intermediate
1147: .seealso: PCASMSetType(), PCASMGetType(), PCASMGetLocalType(), PCASM, PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1148: @*/
1149: PetscErrorCode PCASMSetLocalType(PC pc, PCCompositeType type)
1150: {
1153: PetscTryMethod(pc, "PCASMSetLocalType_C", (PC, PCCompositeType), (pc, type));
1154: return 0;
1155: }
1157: /*@
1158: PCASMGetLocalType - Gets the type of composition used for local problems in the additive Schwarz method.
1160: Logically Collective on pc
1162: Input Parameter:
1163: . pc - the preconditioner context
1165: Output Parameter:
1166: . type - type of composition, one of
1167: .vb
1168: PC_COMPOSITE_ADDITIVE - local additive combination
1169: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1170: .ve
1172: Options Database Key:
1173: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1175: Level: intermediate
1177: .seealso: PCASMSetType(), PCASMGetType(), PCASMSetLocalType(), PCASMCreate(), PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1178: @*/
1179: PetscErrorCode PCASMGetLocalType(PC pc, PCCompositeType *type)
1180: {
1183: PetscUseMethod(pc, "PCASMGetLocalType_C", (PC, PCCompositeType *), (pc, type));
1184: return 0;
1185: }
1187: /*@
1188: PCASMSetSortIndices - Determines whether subdomain indices are sorted.
1190: Logically Collective on pc
1192: Input Parameters:
1193: + pc - the preconditioner context
1194: - doSort - sort the subdomain indices
1196: Level: intermediate
1198: .seealso: PCASMSetLocalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1199: PCASMCreateSubdomains2D()
1200: @*/
1201: PetscErrorCode PCASMSetSortIndices(PC pc,PetscBool doSort)
1202: {
1205: PetscTryMethod(pc,"PCASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));
1206: return 0;
1207: }
1209: /*@C
1210: PCASMGetSubKSP - Gets the local KSP contexts for all blocks on
1211: this processor.
1213: Collective on pc iff first_local is requested
1215: Input Parameter:
1216: . pc - the preconditioner context
1218: Output Parameters:
1219: + n_local - the number of blocks on this processor or NULL
1220: . first_local - the global number of the first block on this processor or NULL,
1221: all processors must request or all must pass NULL
1222: - ksp - the array of KSP contexts
1224: Note:
1225: After PCASMGetSubKSP() the array of KSPes is not to be freed.
1227: You must call KSPSetUp() before calling PCASMGetSubKSP().
1229: Fortran note:
1230: 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.
1232: Level: advanced
1234: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap(),
1235: PCASMCreateSubdomains2D(),
1236: @*/
1237: PetscErrorCode PCASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
1238: {
1240: PetscUseMethod(pc,"PCASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
1241: return 0;
1242: }
1244: /* -------------------------------------------------------------------------------------*/
1245: /*MC
1246: PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1247: its own KSP object.
1249: Options Database Keys:
1250: + -pc_asm_blocks <blks> - Sets total blocks
1251: . -pc_asm_overlap <ovl> - Sets overlap
1252: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type, default is restrict
1253: - -pc_asm_local_type [additive, multiplicative] - Sets ASM type, default is additive
1255: IMPORTANT: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you
1256: will get a different convergence rate due to the default option of -pc_asm_type restrict. Use
1257: -pc_asm_type basic to use the standard ASM.
1259: Notes:
1260: Each processor can have one or more blocks, but a block cannot be shared by more
1261: than one processor. Use PCGASM for subdomains shared by multiple processes. Defaults to one block per processor.
1263: To set options on the solvers for each block append -sub_ to all the KSP, and PC
1264: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
1266: To set the options on the solvers separate for each block call PCASMGetSubKSP()
1267: and set the options directly on the resulting KSP object (you can access its PC
1268: with KSPGetPC())
1270: Level: beginner
1272: References:
1273: + * - M Dryja, OB Widlund, An additive variant of the Schwarz alternating method for the case of many subregions
1274: Courant Institute, New York University Technical report
1275: - * - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations,
1276: Cambridge University Press.
1278: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
1279: PCBJACOBI, PCASMGetSubKSP(), PCASMSetLocalSubdomains(), PCASMType, PCASMGetType(), PCASMSetLocalType(), PCASMGetLocalType()
1280: PCASMSetTotalSubdomains(), PCSetModifySubMatrices(), PCASMSetOverlap(), PCASMSetType(), PCCompositeType
1282: M*/
1284: PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1285: {
1286: PC_ASM *osm;
1288: PetscNewLog(pc,&osm);
1290: osm->n = PETSC_DECIDE;
1291: osm->n_local = 0;
1292: osm->n_local_true = PETSC_DECIDE;
1293: osm->overlap = 1;
1294: osm->ksp = NULL;
1295: osm->restriction = NULL;
1296: osm->lprolongation = NULL;
1297: osm->lrestriction = NULL;
1298: osm->x = NULL;
1299: osm->y = NULL;
1300: osm->is = NULL;
1301: osm->is_local = NULL;
1302: osm->mat = NULL;
1303: osm->pmat = NULL;
1304: osm->type = PC_ASM_RESTRICT;
1305: osm->loctype = PC_COMPOSITE_ADDITIVE;
1306: osm->sort_indices = PETSC_TRUE;
1307: osm->dm_subdomains = PETSC_FALSE;
1308: osm->sub_mat_type = NULL;
1310: pc->data = (void*)osm;
1311: pc->ops->apply = PCApply_ASM;
1312: pc->ops->matapply = PCMatApply_ASM;
1313: pc->ops->applytranspose = PCApplyTranspose_ASM;
1314: pc->ops->setup = PCSetUp_ASM;
1315: pc->ops->reset = PCReset_ASM;
1316: pc->ops->destroy = PCDestroy_ASM;
1317: pc->ops->setfromoptions = PCSetFromOptions_ASM;
1318: pc->ops->setuponblocks = PCSetUpOnBlocks_ASM;
1319: pc->ops->view = PCView_ASM;
1320: pc->ops->applyrichardson = NULL;
1322: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",PCASMSetLocalSubdomains_ASM);
1323: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",PCASMSetTotalSubdomains_ASM);
1324: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetOverlap_C",PCASMSetOverlap_ASM);
1325: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetType_C",PCASMSetType_ASM);
1326: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetType_C",PCASMGetType_ASM);
1327: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalType_C",PCASMSetLocalType_ASM);
1328: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetLocalType_C",PCASMGetLocalType_ASM);
1329: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSortIndices_C",PCASMSetSortIndices_ASM);
1330: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubKSP_C",PCASMGetSubKSP_ASM);
1331: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubMatType_C",PCASMGetSubMatType_ASM);
1332: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSubMatType_C",PCASMSetSubMatType_ASM);
1333: return 0;
1334: }
1336: /*@C
1337: PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1338: preconditioner for any problem on a general grid.
1340: Collective
1342: Input Parameters:
1343: + A - The global matrix operator
1344: - n - the number of local blocks
1346: Output Parameters:
1347: . outis - the array of index sets defining the subdomains
1349: Level: advanced
1351: Note: this generates nonoverlapping subdomains; the PCASM will generate the overlap
1352: from these if you use PCASMSetLocalSubdomains()
1354: In the Fortran version you must provide the array outis[] already allocated of length n.
1356: .seealso: PCASMSetLocalSubdomains(), PCASMDestroySubdomains()
1357: @*/
1358: PetscErrorCode PCASMCreateSubdomains(Mat A, PetscInt n, IS* outis[])
1359: {
1360: MatPartitioning mpart;
1361: const char *prefix;
1362: PetscInt i,j,rstart,rend,bs;
1363: PetscBool hasop, isbaij = PETSC_FALSE,foundpart = PETSC_FALSE;
1364: Mat Ad = NULL, adj;
1365: IS ispart,isnumb,*is;
1371: /* Get prefix, row distribution, and block size */
1372: MatGetOptionsPrefix(A,&prefix);
1373: MatGetOwnershipRange(A,&rstart,&rend);
1374: MatGetBlockSize(A,&bs);
1377: /* Get diagonal block from matrix if possible */
1378: MatHasOperation(A,MATOP_GET_DIAGONAL_BLOCK,&hasop);
1379: if (hasop) {
1380: MatGetDiagonalBlock(A,&Ad);
1381: }
1382: if (Ad) {
1383: PetscObjectBaseTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij);
1384: if (!isbaij) 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: MatPartitioningCreate(PETSC_COMM_SELF,&mpart);
1390: PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);
1391: MatPartitioningSetFromOptions(mpart);
1392: PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGCURRENT,&match);
1393: if (!match) {
1394: PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGSQUARE,&match);
1395: }
1396: if (!match) { /* assume a "good" partitioner is available */
1397: PetscInt na;
1398: const PetscInt *ia,*ja;
1399: MatGetRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1400: if (done) {
1401: /* Build adjacency matrix by hand. Unfortunately a call to
1402: MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1403: remove the block-aij structure and we cannot expect
1404: MatPartitioning to split vertices as we need */
1405: PetscInt i,j,len,nnz,cnt,*iia=NULL,*jja=NULL;
1406: const PetscInt *row;
1407: nnz = 0;
1408: for (i=0; i<na; i++) { /* count number of nonzeros */
1409: len = ia[i+1] - ia[i];
1410: row = ja + ia[i];
1411: for (j=0; j<len; j++) {
1412: if (row[j] == i) { /* don't count diagonal */
1413: len--; break;
1414: }
1415: }
1416: nnz += len;
1417: }
1418: PetscMalloc1(na+1,&iia);
1419: PetscMalloc1(nnz,&jja);
1420: nnz = 0;
1421: iia[0] = 0;
1422: for (i=0; i<na; i++) { /* fill adjacency */
1423: cnt = 0;
1424: len = ia[i+1] - ia[i];
1425: row = ja + ia[i];
1426: for (j=0; j<len; j++) {
1427: if (row[j] != i) { /* if not diagonal */
1428: jja[nnz+cnt++] = row[j];
1429: }
1430: }
1431: nnz += cnt;
1432: iia[i+1] = nnz;
1433: }
1434: /* Partitioning of the adjacency matrix */
1435: MatCreateMPIAdj(PETSC_COMM_SELF,na,na,iia,jja,NULL,&adj);
1436: MatPartitioningSetAdjacency(mpart,adj);
1437: MatPartitioningSetNParts(mpart,n);
1438: MatPartitioningApply(mpart,&ispart);
1439: ISPartitioningToNumbering(ispart,&isnumb);
1440: MatDestroy(&adj);
1441: foundpart = PETSC_TRUE;
1442: }
1443: MatRestoreRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1444: }
1445: MatPartitioningDestroy(&mpart);
1446: }
1448: PetscMalloc1(n,&is);
1449: *outis = is;
1451: if (!foundpart) {
1453: /* Partitioning by contiguous chunks of rows */
1455: PetscInt mbs = (rend-rstart)/bs;
1456: PetscInt start = rstart;
1457: for (i=0; i<n; i++) {
1458: PetscInt count = (mbs/n + ((mbs % n) > i)) * bs;
1459: ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);
1460: start += count;
1461: }
1463: } else {
1465: /* Partitioning by adjacency of diagonal block */
1467: const PetscInt *numbering;
1468: PetscInt *count,nidx,*indices,*newidx,start=0;
1469: /* Get node count in each partition */
1470: PetscMalloc1(n,&count);
1471: ISPartitioningCount(ispart,n,count);
1472: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1473: for (i=0; i<n; i++) count[i] *= bs;
1474: }
1475: /* Build indices from node numbering */
1476: ISGetLocalSize(isnumb,&nidx);
1477: PetscMalloc1(nidx,&indices);
1478: for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */
1479: ISGetIndices(isnumb,&numbering);
1480: PetscSortIntWithPermutation(nidx,numbering,indices);
1481: ISRestoreIndices(isnumb,&numbering);
1482: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1483: PetscMalloc1(nidx*bs,&newidx);
1484: for (i=0; i<nidx; i++) {
1485: for (j=0; j<bs; j++) newidx[i*bs+j] = indices[i]*bs + j;
1486: }
1487: PetscFree(indices);
1488: nidx *= bs;
1489: indices = newidx;
1490: }
1491: /* Shift to get global indices */
1492: for (i=0; i<nidx; i++) indices[i] += rstart;
1494: /* Build the index sets for each block */
1495: for (i=0; i<n; i++) {
1496: ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);
1497: ISSort(is[i]);
1498: start += count[i];
1499: }
1501: PetscFree(count);
1502: PetscFree(indices);
1503: ISDestroy(&isnumb);
1504: ISDestroy(&ispart);
1506: }
1507: return 0;
1508: }
1510: /*@C
1511: PCASMDestroySubdomains - Destroys the index sets created with
1512: PCASMCreateSubdomains(). Should be called after setting subdomains
1513: with PCASMSetLocalSubdomains().
1515: Collective
1517: Input Parameters:
1518: + n - the number of index sets
1519: . is - the array of index sets
1520: - is_local - the array of local index sets, can be NULL
1522: Level: advanced
1524: .seealso: PCASMCreateSubdomains(), PCASMSetLocalSubdomains()
1525: @*/
1526: PetscErrorCode PCASMDestroySubdomains(PetscInt n, IS is[], IS is_local[])
1527: {
1528: PetscInt i;
1530: if (n <= 0) return 0;
1531: if (is) {
1533: for (i=0; i<n; i++) ISDestroy(&is[i]);
1534: PetscFree(is);
1535: }
1536: if (is_local) {
1538: for (i=0; i<n; i++) ISDestroy(&is_local[i]);
1539: PetscFree(is_local);
1540: }
1541: return 0;
1542: }
1544: /*@
1545: PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1546: preconditioner for a two-dimensional problem on a regular grid.
1548: Not Collective
1550: Input Parameters:
1551: + m - the number of mesh points in the x direction
1552: . n - the number of mesh points in the y direction
1553: . M - the number of subdomains in the x direction
1554: . N - the number of subdomains in the y direction
1555: . dof - degrees of freedom per node
1556: - overlap - overlap in mesh lines
1558: Output Parameters:
1559: + Nsub - the number of subdomains created
1560: . is - array of index sets defining overlapping (if overlap > 0) subdomains
1561: - is_local - array of index sets defining non-overlapping subdomains
1563: Note:
1564: Presently PCAMSCreateSubdomains2d() is valid only for sequential
1565: preconditioners. More general related routines are
1566: PCASMSetTotalSubdomains() and PCASMSetLocalSubdomains().
1568: Level: advanced
1570: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1571: PCASMSetOverlap()
1572: @*/
1573: PetscErrorCode PCASMCreateSubdomains2D(PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt dof,PetscInt overlap,PetscInt *Nsub,IS **is,IS **is_local)
1574: {
1575: PetscInt i,j,height,width,ystart,xstart,yleft,yright,xleft,xright,loc_outer;
1576: PetscInt nidx,*idx,loc,ii,jj,count;
1580: *Nsub = N*M;
1581: PetscMalloc1(*Nsub,is);
1582: 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 */
1588: yleft = ystart - overlap; if (yleft < 0) yleft = 0;
1589: yright = ystart + height + overlap; if (yright > n) yright = n;
1590: xstart = 0;
1591: for (j=0; j<M; j++) {
1592: width = m/M + ((m % M) > j); /* width of subdomain */
1594: xleft = xstart - overlap; if (xleft < 0) xleft = 0;
1595: xright = xstart + width + overlap; if (xright > m) xright = m;
1596: nidx = (xright - xleft)*(yright - yleft);
1597: PetscMalloc1(nidx,&idx);
1598: loc = 0;
1599: for (ii=yleft; ii<yright; ii++) {
1600: count = m*ii + xleft;
1601: for (jj=xleft; jj<xright; jj++) idx[loc++] = count++;
1602: }
1603: ISCreateGeneral(PETSC_COMM_SELF,nidx,idx,PETSC_COPY_VALUES,(*is)+loc_outer);
1604: if (overlap == 0) {
1605: PetscObjectReference((PetscObject)(*is)[loc_outer]);
1607: (*is_local)[loc_outer] = (*is)[loc_outer];
1608: } else {
1609: for (loc=0,ii=ystart; ii<ystart+height; ii++) {
1610: for (jj=xstart; jj<xstart+width; jj++) {
1611: idx[loc++] = m*ii + jj;
1612: }
1613: }
1614: ISCreateGeneral(PETSC_COMM_SELF,loc,idx,PETSC_COPY_VALUES,*is_local+loc_outer);
1615: }
1616: PetscFree(idx);
1617: xstart += width;
1618: loc_outer++;
1619: }
1620: ystart += height;
1621: }
1622: for (i=0; i<*Nsub; i++) ISSort((*is)[i]);
1623: return 0;
1624: }
1626: /*@C
1627: PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
1628: only) for the additive Schwarz preconditioner.
1630: Not Collective
1632: Input Parameter:
1633: . pc - the preconditioner context
1635: Output Parameters:
1636: + n - if requested, the number of subdomains for this processor (default value = 1)
1637: . is - if requested, the index sets that define the subdomains for this processor
1638: - is_local - if requested, the index sets that define the local part of the subdomains for this processor (can be NULL)
1640: Notes:
1641: The IS numbering is in the parallel, global numbering of the vector.
1643: Level: advanced
1645: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1646: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubmatrices()
1647: @*/
1648: PetscErrorCode PCASMGetLocalSubdomains(PC pc,PetscInt *n,IS *is[],IS *is_local[])
1649: {
1650: PC_ASM *osm = (PC_ASM*)pc->data;
1651: PetscBool match;
1657: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1659: if (n) *n = osm->n_local_true;
1660: if (is) *is = osm->is;
1661: if (is_local) *is_local = osm->is_local;
1662: return 0;
1663: }
1665: /*@C
1666: PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1667: only) for the additive Schwarz preconditioner.
1669: Not Collective
1671: Input Parameter:
1672: . pc - the preconditioner context
1674: Output Parameters:
1675: + n - if requested, the number of matrices for this processor (default value = 1)
1676: - mat - if requested, the matrices
1678: Level: advanced
1680: Notes:
1681: Call after PCSetUp() (or KSPSetUp()) but before PCApply() and before PCSetUpOnBlocks())
1683: Usually one would use PCSetModifySubMatrices() to change the submatrices in building the preconditioner.
1685: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1686: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains(), PCSetModifySubMatrices()
1687: @*/
1688: PetscErrorCode PCASMGetLocalSubmatrices(PC pc,PetscInt *n,Mat *mat[])
1689: {
1690: PC_ASM *osm;
1691: PetscBool match;
1697: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1698: if (!match) {
1699: if (n) *n = 0;
1700: if (mat) *mat = NULL;
1701: } else {
1702: osm = (PC_ASM*)pc->data;
1703: if (n) *n = osm->n_local_true;
1704: if (mat) *mat = osm->pmat;
1705: }
1706: return 0;
1707: }
1709: /*@
1710: PCASMSetDMSubdomains - Indicates whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1712: Logically Collective
1714: Input Parameters:
1715: + pc - the preconditioner
1716: - flg - boolean indicating whether to use subdomains defined by the DM
1718: Options Database Key:
1719: . -pc_asm_dm_subdomains <bool> - use subdomains defined by the DM
1721: Level: intermediate
1723: Notes:
1724: PCASMSetTotalSubdomains() and PCASMSetOverlap() take precedence over PCASMSetDMSubdomains(),
1725: so setting either of the first two effectively turns the latter off.
1727: .seealso: PCASMGetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1728: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1729: @*/
1730: PetscErrorCode PCASMSetDMSubdomains(PC pc,PetscBool flg)
1731: {
1732: PC_ASM *osm = (PC_ASM*)pc->data;
1733: PetscBool match;
1738: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1739: if (match) {
1740: osm->dm_subdomains = flg;
1741: }
1742: return 0;
1743: }
1745: /*@
1746: PCASMGetDMSubdomains - Returns flag indicating whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1747: Not Collective
1749: Input Parameter:
1750: . pc - the preconditioner
1752: Output Parameter:
1753: . flg - boolean indicating whether to use subdomains defined by the DM
1755: Level: intermediate
1757: .seealso: PCASMSetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1758: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1759: @*/
1760: PetscErrorCode PCASMGetDMSubdomains(PC pc,PetscBool* flg)
1761: {
1762: PC_ASM *osm = (PC_ASM*)pc->data;
1763: PetscBool match;
1767: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1768: if (match) *flg = osm->dm_subdomains;
1769: else *flg = PETSC_FALSE;
1770: return 0;
1771: }
1773: /*@
1774: PCASMGetSubMatType - Gets the matrix type used for ASM subsolves, as a string.
1776: Not Collective
1778: Input Parameter:
1779: . pc - the PC
1781: Output Parameter:
1782: . -pc_asm_sub_mat_type - name of matrix type
1784: Level: advanced
1786: .seealso: PCASMSetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1787: @*/
1788: PetscErrorCode PCASMGetSubMatType(PC pc,MatType *sub_mat_type)
1789: {
1791: PetscTryMethod(pc,"PCASMGetSubMatType_C",(PC,MatType*),(pc,sub_mat_type));
1792: return 0;
1793: }
1795: /*@
1796: PCASMSetSubMatType - Set the type of matrix used for ASM subsolves
1798: Collective on Mat
1800: Input Parameters:
1801: + pc - the PC object
1802: - sub_mat_type - matrix type
1804: Options Database Key:
1805: . -pc_asm_sub_mat_type <sub_mat_type> - Sets the matrix type used for subsolves, for example, seqaijviennacl. If you specify a base name like aijviennacl, the corresponding sequential type is assumed.
1807: Notes:
1808: See "${PETSC_DIR}/include/petscmat.h" for available types
1810: Level: advanced
1812: .seealso: PCASMGetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1813: @*/
1814: PetscErrorCode PCASMSetSubMatType(PC pc,MatType sub_mat_type)
1815: {
1817: PetscTryMethod(pc,"PCASMSetSubMatType_C",(PC,MatType),(pc,sub_mat_type));
1818: return 0;
1819: }