Actual source code: asm.c
petsc-3.12.5 2020-03-29
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: */
12: #include <petsc/private/pcimpl.h>
13: #include <petscdm.h>
15: typedef struct {
16: PetscInt n, n_local, n_local_true;
17: PetscInt overlap; /* overlap requested by user */
18: KSP *ksp; /* linear solvers for each block */
19: VecScatter restriction; /* mapping from global to overlapping (process) subdomain*/
20: VecScatter *lrestriction; /* mapping from subregion to overlapping (process) subdomain */
21: VecScatter *lprolongation; /* mapping from non-overlapping subregion to overlapping (process) subdomain; used for restrict additive version of algorithms */
22: Vec lx, ly; /* work vectors */
23: Vec *x,*y; /* work vectors */
24: IS lis; /* index set that defines each overlapping multiplicative (process) subdomain */
25: IS *is; /* index set that defines each overlapping subdomain */
26: IS *is_local; /* index set that defines each non-overlapping subdomain, may be NULL */
27: Mat *mat,*pmat; /* mat is not currently used */
28: PCASMType type; /* use reduced interpolation, restriction or both */
29: PetscBool type_set; /* if user set this value (so won't change it for symmetric problems) */
30: PetscBool same_local_solves; /* flag indicating whether all local solvers are same */
31: PetscBool sort_indices; /* flag to sort subdomain indices */
32: PetscBool dm_subdomains; /* whether DM is allowed to define subdomains */
33: PCCompositeType loctype; /* the type of composition for local solves */
34: MatType sub_mat_type; /* the type of Mat used for subdomain solves (can be MATSAME or NULL) */
35: /* For multiplicative solve */
36: Mat *lmats; /* submatrices for overlapping multiplicative (process) subdomain */
37: } PC_ASM;
39: static PetscErrorCode PCView_ASM(PC pc,PetscViewer viewer)
40: {
41: PC_ASM *osm = (PC_ASM*)pc->data;
43: PetscMPIInt rank;
44: PetscInt i,bsz;
45: PetscBool iascii,isstring;
46: PetscViewer sviewer;
49: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
50: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
51: if (iascii) {
52: char overlaps[256] = "user-defined overlap",blocks[256] = "total subdomain blocks not yet set";
53: if (osm->overlap >= 0) {PetscSNPrintf(overlaps,sizeof(overlaps),"amount of overlap = %D",osm->overlap);}
54: if (osm->n > 0) {PetscSNPrintf(blocks,sizeof(blocks),"total subdomain blocks = %D",osm->n);}
55: PetscViewerASCIIPrintf(viewer," %s, %s\n",blocks,overlaps);
56: PetscViewerASCIIPrintf(viewer," restriction/interpolation type - %s\n",PCASMTypes[osm->type]);
57: if (osm->dm_subdomains) {PetscViewerASCIIPrintf(viewer," Additive Schwarz: using DM to define subdomains\n");}
58: if (osm->loctype != PC_COMPOSITE_ADDITIVE) {PetscViewerASCIIPrintf(viewer," Additive Schwarz: local solve composition type - %s\n",PCCompositeTypes[osm->loctype]);}
59: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
60: if (osm->same_local_solves) {
61: if (osm->ksp) {
62: PetscViewerASCIIPrintf(viewer," Local solve is same for all blocks, in the following KSP and PC objects:\n");
63: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
64: if (!rank) {
65: PetscViewerASCIIPushTab(viewer);
66: KSPView(osm->ksp[0],sviewer);
67: PetscViewerASCIIPopTab(viewer);
68: }
69: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
70: }
71: } else {
72: PetscViewerASCIIPushSynchronized(viewer);
73: PetscViewerASCIISynchronizedPrintf(viewer," [%d] number of local blocks = %D\n",(int)rank,osm->n_local_true);
74: PetscViewerFlush(viewer);
75: PetscViewerASCIIPrintf(viewer," Local solve info for each block is in the following KSP and PC objects:\n");
76: PetscViewerASCIIPushTab(viewer);
77: PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
78: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
79: for (i=0; i<osm->n_local_true; i++) {
80: ISGetLocalSize(osm->is[i],&bsz);
81: PetscViewerASCIISynchronizedPrintf(sviewer,"[%d] local block number %D, size = %D\n",(int)rank,i,bsz);
82: KSPView(osm->ksp[i],sviewer);
83: PetscViewerASCIISynchronizedPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
84: }
85: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
86: PetscViewerASCIIPopTab(viewer);
87: PetscViewerFlush(viewer);
88: PetscViewerASCIIPopSynchronized(viewer);
89: }
90: } else if (isstring) {
91: PetscViewerStringSPrintf(viewer," blocks=%D, overlap=%D, type=%s",osm->n,osm->overlap,PCASMTypes[osm->type]);
92: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
93: if (osm->ksp) {KSPView(osm->ksp[0],sviewer);}
94: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
95: }
96: return(0);
97: }
99: static PetscErrorCode PCASMPrintSubdomains(PC pc)
100: {
101: PC_ASM *osm = (PC_ASM*)pc->data;
102: const char *prefix;
103: char fname[PETSC_MAX_PATH_LEN+1];
104: PetscViewer viewer, sviewer;
105: char *s;
106: PetscInt i,j,nidx;
107: const PetscInt *idx;
108: PetscMPIInt rank, size;
112: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
113: MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
114: PCGetOptionsPrefix(pc,&prefix);
115: PetscOptionsGetString(NULL,prefix,"-pc_asm_print_subdomains",fname,PETSC_MAX_PATH_LEN,NULL);
116: if (fname[0] == 0) { PetscStrcpy(fname,"stdout"); };
117: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),fname,&viewer);
118: for (i=0; i<osm->n_local; i++) {
119: if (i < osm->n_local_true) {
120: ISGetLocalSize(osm->is[i],&nidx);
121: ISGetIndices(osm->is[i],&idx);
122: /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
123: #define len 16*(nidx+1)+512
124: PetscMalloc1(len,&s);
125: PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer);
126: #undef len
127: PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D with overlap:\n", rank, size, i);
128: for (j=0; j<nidx; j++) {
129: PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
130: }
131: ISRestoreIndices(osm->is[i],&idx);
132: PetscViewerStringSPrintf(sviewer,"\n");
133: PetscViewerDestroy(&sviewer);
134: PetscViewerASCIIPushSynchronized(viewer);
135: PetscViewerASCIISynchronizedPrintf(viewer, s);
136: PetscViewerFlush(viewer);
137: PetscViewerASCIIPopSynchronized(viewer);
138: PetscFree(s);
139: if (osm->is_local) {
140: /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
141: #define len 16*(nidx+1)+512
142: PetscMalloc1(len, &s);
143: PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer);
144: #undef len
145: PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D without overlap:\n", rank, size, i);
146: ISGetLocalSize(osm->is_local[i],&nidx);
147: ISGetIndices(osm->is_local[i],&idx);
148: for (j=0; j<nidx; j++) {
149: PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
150: }
151: ISRestoreIndices(osm->is_local[i],&idx);
152: PetscViewerStringSPrintf(sviewer,"\n");
153: PetscViewerDestroy(&sviewer);
154: PetscViewerASCIIPushSynchronized(viewer);
155: PetscViewerASCIISynchronizedPrintf(viewer, s);
156: PetscViewerFlush(viewer);
157: PetscViewerASCIIPopSynchronized(viewer);
158: PetscFree(s);
159: }
160: } else {
161: /* Participate in collective viewer calls. */
162: PetscViewerASCIIPushSynchronized(viewer);
163: PetscViewerFlush(viewer);
164: PetscViewerASCIIPopSynchronized(viewer);
165: /* Assume either all ranks have is_local or none do. */
166: if (osm->is_local) {
167: PetscViewerASCIIPushSynchronized(viewer);
168: PetscViewerFlush(viewer);
169: PetscViewerASCIIPopSynchronized(viewer);
170: }
171: }
172: }
173: PetscViewerFlush(viewer);
174: PetscViewerDestroy(&viewer);
175: return(0);
176: }
178: static PetscErrorCode PCSetUp_ASM(PC pc)
179: {
180: PC_ASM *osm = (PC_ASM*)pc->data;
182: PetscBool symset,flg;
183: PetscInt i,m,m_local;
184: MatReuse scall = MAT_REUSE_MATRIX;
185: IS isl;
186: KSP ksp;
187: PC subpc;
188: const char *prefix,*pprefix;
189: Vec vec;
190: DM *domain_dm = NULL;
193: if (!pc->setupcalled) {
194: PetscInt m;
196: if (!osm->type_set) {
197: MatIsSymmetricKnown(pc->pmat,&symset,&flg);
198: if (symset && flg) osm->type = PC_ASM_BASIC;
199: }
201: /* Note: if subdomains have been set either via PCASMSetTotalSubdomains() or via PCASMSetLocalSubdomains(), osm->n_local_true will not be PETSC_DECIDE */
202: if (osm->n_local_true == PETSC_DECIDE) {
203: /* no subdomains given */
204: /* try pc->dm first, if allowed */
205: if (osm->dm_subdomains && pc->dm) {
206: PetscInt num_domains, d;
207: char **domain_names;
208: IS *inner_domain_is, *outer_domain_is;
209: DMCreateDomainDecomposition(pc->dm, &num_domains, &domain_names, &inner_domain_is, &outer_domain_is, &domain_dm);
210: osm->overlap = -1; /* We do not want to increase the overlap of the IS.
211: A future improvement of this code might allow one to use
212: DM-defined subdomains and also increase the overlap,
213: but that is not currently supported */
214: if (num_domains) {
215: PCASMSetLocalSubdomains(pc, num_domains, outer_domain_is, inner_domain_is);
216: }
217: for (d = 0; d < num_domains; ++d) {
218: if (domain_names) {PetscFree(domain_names[d]);}
219: if (inner_domain_is) {ISDestroy(&inner_domain_is[d]);}
220: if (outer_domain_is) {ISDestroy(&outer_domain_is[d]);}
221: }
222: PetscFree(domain_names);
223: PetscFree(inner_domain_is);
224: PetscFree(outer_domain_is);
225: }
226: if (osm->n_local_true == PETSC_DECIDE) {
227: /* still no subdomains; use one subdomain per processor */
228: osm->n_local_true = 1;
229: }
230: }
231: { /* determine the global and max number of subdomains */
232: struct {PetscInt max,sum;} inwork,outwork;
233: PetscMPIInt size;
235: inwork.max = osm->n_local_true;
236: inwork.sum = osm->n_local_true;
237: MPIU_Allreduce(&inwork,&outwork,1,MPIU_2INT,MPIU_MAXSUM_OP,PetscObjectComm((PetscObject)pc));
238: osm->n_local = outwork.max;
239: osm->n = outwork.sum;
241: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
242: if (outwork.max == 1 && outwork.sum == size) {
243: /* osm->n_local_true = 1 on all processes, set this option may enable use of optimized MatCreateSubMatrices() implementation */
244: MatSetOption(pc->pmat,MAT_SUBMAT_SINGLEIS,PETSC_TRUE);
245: }
246: }
247: if (!osm->is) { /* create the index sets */
248: PCASMCreateSubdomains(pc->pmat,osm->n_local_true,&osm->is);
249: }
250: if (osm->n_local_true > 1 && !osm->is_local) {
251: PetscMalloc1(osm->n_local_true,&osm->is_local);
252: for (i=0; i<osm->n_local_true; i++) {
253: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
254: ISDuplicate(osm->is[i],&osm->is_local[i]);
255: ISCopy(osm->is[i],osm->is_local[i]);
256: } else {
257: PetscObjectReference((PetscObject)osm->is[i]);
258: osm->is_local[i] = osm->is[i];
259: }
260: }
261: }
262: PCGetOptionsPrefix(pc,&prefix);
263: flg = PETSC_FALSE;
264: PetscOptionsGetBool(NULL,prefix,"-pc_asm_print_subdomains",&flg,NULL);
265: if (flg) { PCASMPrintSubdomains(pc); }
267: if (osm->overlap > 0) {
268: /* Extend the "overlapping" regions by a number of steps */
269: MatIncreaseOverlap(pc->pmat,osm->n_local_true,osm->is,osm->overlap);
270: }
271: if (osm->sort_indices) {
272: for (i=0; i<osm->n_local_true; i++) {
273: ISSort(osm->is[i]);
274: if (osm->is_local) {
275: ISSort(osm->is_local[i]);
276: }
277: }
278: }
280: if (!osm->ksp) {
281: /* Create the local solvers */
282: PetscMalloc1(osm->n_local_true,&osm->ksp);
283: if (domain_dm) {
284: PetscInfo(pc,"Setting up ASM subproblems using the embedded DM\n");
285: }
286: for (i=0; i<osm->n_local_true; i++) {
287: KSPCreate(PETSC_COMM_SELF,&ksp);
288: KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
289: PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
290: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
291: KSPSetType(ksp,KSPPREONLY);
292: KSPGetPC(ksp,&subpc);
293: PCGetOptionsPrefix(pc,&prefix);
294: KSPSetOptionsPrefix(ksp,prefix);
295: KSPAppendOptionsPrefix(ksp,"sub_");
296: if (domain_dm) {
297: KSPSetDM(ksp, domain_dm[i]);
298: KSPSetDMActive(ksp, PETSC_FALSE);
299: DMDestroy(&domain_dm[i]);
300: }
301: osm->ksp[i] = ksp;
302: }
303: if (domain_dm) {
304: PetscFree(domain_dm);
305: }
306: }
308: ISConcatenate(PETSC_COMM_SELF, osm->n_local_true, osm->is, &osm->lis);
309: ISSortRemoveDups(osm->lis);
310: ISGetLocalSize(osm->lis, &m);
311: VecCreateSeq(PETSC_COMM_SELF, m, &osm->lx);
312: VecDuplicate(osm->lx, &osm->ly);
314: scall = MAT_INITIAL_MATRIX;
315: } else {
316: /*
317: Destroy the blocks from the previous iteration
318: */
319: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
320: MatDestroyMatrices(osm->n_local_true,&osm->pmat);
321: scall = MAT_INITIAL_MATRIX;
322: }
323: }
325: /*
326: Extract out the submatrices
327: */
328: MatCreateSubMatrices(pc->pmat,osm->n_local_true,osm->is,osm->is,scall,&osm->pmat);
329: if (scall == MAT_INITIAL_MATRIX) {
330: PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
331: for (i=0; i<osm->n_local_true; i++) {
332: PetscLogObjectParent((PetscObject)pc,(PetscObject)osm->pmat[i]);
333: PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);
334: }
335: }
337: /* Convert the types of the submatrices (if needbe) */
338: if (osm->sub_mat_type) {
339: for (i=0; i<osm->n_local_true; i++) {
340: MatConvert(osm->pmat[i],osm->sub_mat_type,MAT_INPLACE_MATRIX,&(osm->pmat[i]));
341: }
342: }
344: if(!pc->setupcalled){
345: /* Create the local work vectors (from the local matrices) and scatter contexts */
346: MatCreateVecs(pc->pmat,&vec,0);
348: if (osm->is_local && (osm->type == PC_ASM_INTERPOLATE || osm->type == PC_ASM_NONE )) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Cannot use interpolate or none PCASMType if is_local was provided to PCASMSetLocalSubdomains()");
349: if (osm->is_local && osm->type == PC_ASM_RESTRICT && osm->loctype == PC_COMPOSITE_ADDITIVE) {
350: PetscMalloc1(osm->n_local_true,&osm->lprolongation);
351: }
352: PetscMalloc1(osm->n_local_true,&osm->lrestriction);
353: PetscMalloc1(osm->n_local_true,&osm->x);
354: PetscMalloc1(osm->n_local_true,&osm->y);
356: ISGetLocalSize(osm->lis,&m);
357: ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
358: VecScatterCreate(vec,osm->lis,osm->lx,isl,&osm->restriction);
359: ISDestroy(&isl);
362: for (i=0; i<osm->n_local_true; ++i) {
363: ISLocalToGlobalMapping ltog;
364: IS isll;
365: const PetscInt *idx_is;
366: PetscInt *idx_lis,nout;
368: ISGetLocalSize(osm->is[i],&m);
369: MatCreateVecs(osm->pmat[i],&osm->x[i],NULL);
370: VecDuplicate(osm->x[i],&osm->y[i]);
372: /* generate a scatter from ly to y[i] picking all the overlapping is[i] entries */
373: ISLocalToGlobalMappingCreateIS(osm->lis,<og);
374: ISGetLocalSize(osm->is[i],&m);
375: ISGetIndices(osm->is[i], &idx_is);
376: PetscMalloc1(m,&idx_lis);
377: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m,idx_is,&nout,idx_lis);
378: if (nout != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is not a subset of lis");
379: ISRestoreIndices(osm->is[i], &idx_is);
380: ISCreateGeneral(PETSC_COMM_SELF,m,idx_lis,PETSC_OWN_POINTER,&isll);
381: ISLocalToGlobalMappingDestroy(<og);
382: ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
383: VecScatterCreate(osm->ly,isll,osm->y[i],isl,&osm->lrestriction[i]);
384: ISDestroy(&isll);
385: ISDestroy(&isl);
386: if (osm->lprolongation) { /* generate a scatter from y[i] to ly picking only the the non-overalapping is_local[i] entries */
387: ISLocalToGlobalMapping ltog;
388: IS isll,isll_local;
389: const PetscInt *idx_local;
390: PetscInt *idx1, *idx2, nout;
392: ISGetLocalSize(osm->is_local[i],&m_local);
393: ISGetIndices(osm->is_local[i], &idx_local);
395: ISLocalToGlobalMappingCreateIS(osm->is[i],<og);
396: PetscMalloc1(m_local,&idx1);
397: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx1);
398: ISLocalToGlobalMappingDestroy(<og);
399: if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of is");
400: ISCreateGeneral(PETSC_COMM_SELF,m_local,idx1,PETSC_OWN_POINTER,&isll);
402: ISLocalToGlobalMappingCreateIS(osm->lis,<og);
403: PetscMalloc1(m_local,&idx2);
404: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx2);
405: ISLocalToGlobalMappingDestroy(<og);
406: if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of lis");
407: ISCreateGeneral(PETSC_COMM_SELF,m_local,idx2,PETSC_OWN_POINTER,&isll_local);
409: ISRestoreIndices(osm->is_local[i], &idx_local);
410: VecScatterCreate(osm->y[i],isll,osm->ly,isll_local,&osm->lprolongation[i]);
412: ISDestroy(&isll);
413: ISDestroy(&isll_local);
414: }
415: }
416: VecDestroy(&vec);
417: }
419: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
420: IS *cis;
421: PetscInt c;
423: PetscMalloc1(osm->n_local_true, &cis);
424: for (c = 0; c < osm->n_local_true; ++c) cis[c] = osm->lis;
425: MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, cis, scall, &osm->lmats);
426: PetscFree(cis);
427: }
429: /* Return control to the user so that the submatrices can be modified (e.g., to apply
430: different boundary conditions for the submatrices than for the global problem) */
431: PCModifySubMatrices(pc,osm->n_local_true,osm->is,osm->is,osm->pmat,pc->modifysubmatricesP);
433: /*
434: Loop over subdomains putting them into local ksp
435: */
436: for (i=0; i<osm->n_local_true; i++) {
437: KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i]);
438: if (!pc->setupcalled) {
439: KSPSetFromOptions(osm->ksp[i]);
440: }
441: }
442: return(0);
443: }
445: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
446: {
447: PC_ASM *osm = (PC_ASM*)pc->data;
448: PetscErrorCode ierr;
449: PetscInt i;
450: KSPConvergedReason reason;
453: for (i=0; i<osm->n_local_true; i++) {
454: KSPSetUp(osm->ksp[i]);
455: KSPGetConvergedReason(osm->ksp[i],&reason);
456: if (reason == KSP_DIVERGED_PC_FAILED) {
457: pc->failedreason = PC_SUBPC_ERROR;
458: }
459: }
460: return(0);
461: }
463: static PetscErrorCode PCApply_ASM(PC pc,Vec x,Vec y)
464: {
465: PC_ASM *osm = (PC_ASM*)pc->data;
467: PetscInt i,n_local_true = osm->n_local_true;
468: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
471: /*
472: Support for limiting the restriction or interpolation to only local
473: subdomain values (leaving the other values 0).
474: */
475: if (!(osm->type & PC_ASM_RESTRICT)) {
476: forward = SCATTER_FORWARD_LOCAL;
477: /* have to zero the work RHS since scatter may leave some slots empty */
478: VecSet(osm->lx, 0.0);
479: }
480: if (!(osm->type & PC_ASM_INTERPOLATE)) {
481: reverse = SCATTER_REVERSE_LOCAL;
482: }
484: if(osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE){
485: /* zero the global and the local solutions */
486: VecZeroEntries(y);
487: VecSet(osm->ly, 0.0);
489: /* Copy the global RHS to local RHS including the ghost nodes */
490: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
491: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
493: /* Restrict local RHS to the overlapping 0-block RHS */
494: VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
495: VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
497: /* do the local solves */
498: for (i = 0; i < n_local_true; ++i) {
500: /* solve the overlapping i-block */
501: PetscLogEventBegin(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
502: KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]);
503: KSPCheckSolve(osm->ksp[i],pc,osm->y[i]);
504: PetscLogEventEnd(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
506: if (osm->lprolongation) { /* interpolate the non-overalapping i-block solution to the local solution (only for restrictive additive) */
507: VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
508: VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
509: }
510: else{ /* interpolate the overalapping i-block solution to the local solution */
511: VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
512: VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
513: }
515: if (i < n_local_true-1) {
516: /* Restrict local RHS to the overlapping (i+1)-block RHS */
517: VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
518: VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
520: if ( osm->loctype == PC_COMPOSITE_MULTIPLICATIVE){
521: /* update the overlapping (i+1)-block RHS using the current local solution */
522: MatMult(osm->lmats[i+1], osm->ly, osm->y[i+1]);
523: VecAXPBY(osm->x[i+1],-1.,1., osm->y[i+1]);
524: }
525: }
526: }
527: /* Add the local solution to the global solution including the ghost nodes */
528: VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
529: VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
530: }else{
531: SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
532: }
533: return(0);
534: }
536: static PetscErrorCode PCApplyTranspose_ASM(PC pc,Vec x,Vec y)
537: {
538: PC_ASM *osm = (PC_ASM*)pc->data;
540: PetscInt i,n_local_true = osm->n_local_true;
541: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
544: /*
545: Support for limiting the restriction or interpolation to only local
546: subdomain values (leaving the other values 0).
548: Note: these are reversed from the PCApply_ASM() because we are applying the
549: transpose of the three terms
550: */
552: if (!(osm->type & PC_ASM_INTERPOLATE)) {
553: forward = SCATTER_FORWARD_LOCAL;
554: /* have to zero the work RHS since scatter may leave some slots empty */
555: VecSet(osm->lx, 0.0);
556: }
557: if (!(osm->type & PC_ASM_RESTRICT)) reverse = SCATTER_REVERSE_LOCAL;
559: /* zero the global and the local solutions */
560: VecZeroEntries(y);
561: VecSet(osm->ly, 0.0);
563: /* Copy the global RHS to local RHS including the ghost nodes */
564: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
565: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
567: /* Restrict local RHS to the overlapping 0-block RHS */
568: VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
569: VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
571: /* do the local solves */
572: for (i = 0; i < n_local_true; ++i) {
574: /* solve the overlapping i-block */
575: PetscLogEventBegin(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
576: KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i]);
577: KSPCheckSolve(osm->ksp[i],pc,osm->y[i]);
578: PetscLogEventEnd(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
580: if (osm->lprolongation) { /* interpolate the non-overalapping i-block solution to the local solution */
581: VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
582: VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
583: }
584: else{ /* interpolate the overalapping i-block solution to the local solution */
585: VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
586: VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
587: }
589: if (i < n_local_true-1) {
590: /* Restrict local RHS to the overlapping (i+1)-block RHS */
591: VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
592: VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
593: }
594: }
595: /* Add the local solution to the global solution including the ghost nodes */
596: VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
597: VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
599: return(0);
601: }
603: static PetscErrorCode PCReset_ASM(PC pc)
604: {
605: PC_ASM *osm = (PC_ASM*)pc->data;
607: PetscInt i;
610: if (osm->ksp) {
611: for (i=0; i<osm->n_local_true; i++) {
612: KSPReset(osm->ksp[i]);
613: }
614: }
615: if (osm->pmat) {
616: if (osm->n_local_true > 0) {
617: MatDestroySubMatrices(osm->n_local_true,&osm->pmat);
618: }
619: }
620: if (osm->lrestriction) {
621: VecScatterDestroy(&osm->restriction);
622: for (i=0; i<osm->n_local_true; i++) {
623: VecScatterDestroy(&osm->lrestriction[i]);
624: if (osm->lprolongation) {VecScatterDestroy(&osm->lprolongation[i]);}
625: VecDestroy(&osm->x[i]);
626: VecDestroy(&osm->y[i]);
627: }
628: PetscFree(osm->lrestriction);
629: if (osm->lprolongation) {PetscFree(osm->lprolongation);}
630: PetscFree(osm->x);
631: PetscFree(osm->y);
633: }
634: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
635: ISDestroy(&osm->lis);
636: VecDestroy(&osm->lx);
637: VecDestroy(&osm->ly);
638: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
639: MatDestroyMatrices(osm->n_local_true, &osm->lmats);
640: }
642: PetscFree(osm->sub_mat_type);
644: osm->is = 0;
645: osm->is_local = 0;
646: return(0);
647: }
649: static PetscErrorCode PCDestroy_ASM(PC pc)
650: {
651: PC_ASM *osm = (PC_ASM*)pc->data;
653: PetscInt i;
656: PCReset_ASM(pc);
657: if (osm->ksp) {
658: for (i=0; i<osm->n_local_true; i++) {
659: KSPDestroy(&osm->ksp[i]);
660: }
661: PetscFree(osm->ksp);
662: }
663: PetscFree(pc->data);
664: return(0);
665: }
667: static PetscErrorCode PCSetFromOptions_ASM(PetscOptionItems *PetscOptionsObject,PC pc)
668: {
669: PC_ASM *osm = (PC_ASM*)pc->data;
671: PetscInt blocks,ovl;
672: PetscBool symset,flg;
673: PCASMType asmtype;
674: PCCompositeType loctype;
675: char sub_mat_type[256];
678: /* set the type to symmetric if matrix is symmetric */
679: if (!osm->type_set && pc->pmat) {
680: MatIsSymmetricKnown(pc->pmat,&symset,&flg);
681: if (symset && flg) osm->type = PC_ASM_BASIC;
682: }
683: PetscOptionsHead(PetscOptionsObject,"Additive Schwarz options");
684: PetscOptionsBool("-pc_asm_dm_subdomains","Use DMCreateDomainDecomposition() to define subdomains","PCASMSetDMSubdomains",osm->dm_subdomains,&osm->dm_subdomains,&flg);
685: PetscOptionsInt("-pc_asm_blocks","Number of subdomains","PCASMSetTotalSubdomains",osm->n,&blocks,&flg);
686: if (flg) {
687: PCASMSetTotalSubdomains(pc,blocks,NULL,NULL);
688: osm->dm_subdomains = PETSC_FALSE;
689: }
690: PetscOptionsInt("-pc_asm_local_blocks","Number of local subdomains","PCASMSetLocalSubdomains",osm->n_local_true,&blocks,&flg);
691: if (flg) {
692: PCASMSetLocalSubdomains(pc,blocks,NULL,NULL);
693: osm->dm_subdomains = PETSC_FALSE;
694: }
695: PetscOptionsInt("-pc_asm_overlap","Number of grid points overlap","PCASMSetOverlap",osm->overlap,&ovl,&flg);
696: if (flg) {
697: PCASMSetOverlap(pc,ovl);
698: osm->dm_subdomains = PETSC_FALSE;
699: }
700: flg = PETSC_FALSE;
701: PetscOptionsEnum("-pc_asm_type","Type of restriction/extension","PCASMSetType",PCASMTypes,(PetscEnum)osm->type,(PetscEnum*)&asmtype,&flg);
702: if (flg) {PCASMSetType(pc,asmtype); }
703: flg = PETSC_FALSE;
704: PetscOptionsEnum("-pc_asm_local_type","Type of local solver composition","PCASMSetLocalType",PCCompositeTypes,(PetscEnum)osm->loctype,(PetscEnum*)&loctype,&flg);
705: if (flg) {PCASMSetLocalType(pc,loctype); }
706: PetscOptionsFList("-pc_asm_sub_mat_type","Subsolve Matrix Type","PCASMSetSubMatType",MatList,NULL,sub_mat_type,256,&flg);
707: if(flg){
708: PCASMSetSubMatType(pc,sub_mat_type);
709: }
710: PetscOptionsTail();
711: return(0);
712: }
714: /*------------------------------------------------------------------------------------*/
716: static PetscErrorCode PCASMSetLocalSubdomains_ASM(PC pc,PetscInt n,IS is[],IS is_local[])
717: {
718: PC_ASM *osm = (PC_ASM*)pc->data;
720: PetscInt i;
723: if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Each process must have 1 or more blocks, n = %D",n);
724: if (pc->setupcalled && (n != osm->n_local_true || is)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetLocalSubdomains() should be called before calling PCSetUp().");
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 = 0;
737: osm->is_local = 0;
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;
768: PetscMPIInt rank,size;
769: PetscInt n;
772: if (N < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of total blocks must be > 0, N = %D",N);
773: if (is || is_local) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Use PCASMSetLocalSubdomains() to set specific index sets\n\they cannot be set globally yet.");
775: /*
776: Split the subdomains equally among all processors
777: */
778: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
779: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
780: n = N/size + ((N % size) > rank);
781: if (!n) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Process %d must have at least one block: total processors %d total blocks %D",(int)rank,(int)size,N);
782: if (pc->setupcalled && n != osm->n_local_true) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCASMSetTotalSubdomains() should be called before PCSetUp().");
783: if (!pc->setupcalled) {
784: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
786: osm->n_local_true = n;
787: osm->is = 0;
788: osm->is_local = 0;
789: }
790: return(0);
791: }
793: static PetscErrorCode PCASMSetOverlap_ASM(PC pc,PetscInt ovl)
794: {
795: PC_ASM *osm = (PC_ASM*)pc->data;
798: if (ovl < 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested");
799: if (pc->setupcalled && ovl != osm->overlap) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetOverlap() should be called before PCSetUp().");
800: if (!pc->setupcalled) osm->overlap = ovl;
801: return(0);
802: }
804: static PetscErrorCode PCASMSetType_ASM(PC pc,PCASMType type)
805: {
806: PC_ASM *osm = (PC_ASM*)pc->data;
809: osm->type = type;
810: osm->type_set = PETSC_TRUE;
811: return(0);
812: }
814: static PetscErrorCode PCASMGetType_ASM(PC pc,PCASMType *type)
815: {
816: PC_ASM *osm = (PC_ASM*)pc->data;
819: *type = osm->type;
820: return(0);
821: }
823: static PetscErrorCode PCASMSetLocalType_ASM(PC pc, PCCompositeType type)
824: {
825: PC_ASM *osm = (PC_ASM *) pc->data;
828: if (type != PC_COMPOSITE_ADDITIVE && type != PC_COMPOSITE_MULTIPLICATIVE) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Only supports additive or multiplicative as the local type");
829: osm->loctype = type;
830: return(0);
831: }
833: static PetscErrorCode PCASMGetLocalType_ASM(PC pc, PCCompositeType *type)
834: {
835: PC_ASM *osm = (PC_ASM *) pc->data;
838: *type = osm->loctype;
839: return(0);
840: }
842: static PetscErrorCode PCASMSetSortIndices_ASM(PC pc,PetscBool doSort)
843: {
844: PC_ASM *osm = (PC_ASM*)pc->data;
847: osm->sort_indices = doSort;
848: return(0);
849: }
851: static PetscErrorCode PCASMGetSubKSP_ASM(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
852: {
853: PC_ASM *osm = (PC_ASM*)pc->data;
857: if (osm->n_local_true < 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Need to call PCSetUp() on PC (or KSPSetUp() on the outer KSP object) before calling here");
859: if (n_local) *n_local = osm->n_local_true;
860: if (first_local) {
861: MPI_Scan(&osm->n_local_true,first_local,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
862: *first_local -= osm->n_local_true;
863: }
864: if (ksp) {
865: /* Assume that local solves are now different; not necessarily
866: true though! This flag is used only for PCView_ASM() */
867: *ksp = osm->ksp;
868: osm->same_local_solves = PETSC_FALSE;
869: }
870: return(0);
871: }
873: static PetscErrorCode PCASMGetSubMatType_ASM(PC pc,MatType *sub_mat_type)
874: {
875: PC_ASM *osm = (PC_ASM*)pc->data;
880: *sub_mat_type = osm->sub_mat_type;
881: return(0);
882: }
884: static PetscErrorCode PCASMSetSubMatType_ASM(PC pc,MatType sub_mat_type)
885: {
886: PetscErrorCode ierr;
887: PC_ASM *osm = (PC_ASM*)pc->data;
891: PetscFree(osm->sub_mat_type);
892: PetscStrallocpy(sub_mat_type,(char**)&osm->sub_mat_type);
893: return(0);
894: }
896: /*@C
897: PCASMSetLocalSubdomains - Sets the local subdomains (for this processor only) for the additive Schwarz preconditioner.
899: Collective on pc
901: Input Parameters:
902: + pc - the preconditioner context
903: . n - the number of subdomains for this processor (default value = 1)
904: . is - the index set that defines the subdomains for this processor
905: (or NULL for PETSc to determine subdomains)
906: - is_local - the index sets that define the local part of the subdomains for this processor, not used unless PCASMType is PC_ASM_RESTRICT
907: (or NULL to not provide these)
909: Options Database Key:
910: To set the total number of subdomain blocks rather than specify the
911: index sets, use the option
912: . -pc_asm_local_blocks <blks> - Sets local blocks
914: Notes:
915: The IS numbering is in the parallel, global numbering of the vector for both is and is_local
917: By default the ASM preconditioner uses 1 block per processor.
919: Use PCASMSetTotalSubdomains() to set the subdomains for all processors.
921: If is_local is provided and PCASMType is PC_ASM_RESTRICT then the solution only over the is_local region is interpolated
922: back to form the global solution (this is the standard restricted additive Schwarz method)
924: If the is_local is provided and PCASMType is PC_ASM_INTERPOLATE or PC_ASM_NONE then an error is generated since there is
925: no code to handle that case.
927: Level: advanced
929: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
930: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), PCASMType, PCASMSetType()
931: @*/
932: PetscErrorCode PCASMSetLocalSubdomains(PC pc,PetscInt n,IS is[],IS is_local[])
933: {
938: PetscTryMethod(pc,"PCASMSetLocalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,is,is_local));
939: return(0);
940: }
942: /*@C
943: PCASMSetTotalSubdomains - Sets the subdomains for all processors for the
944: additive Schwarz preconditioner. Either all or no processors in the
945: PC communicator must call this routine, with the same index sets.
947: Collective on pc
949: Input Parameters:
950: + pc - the preconditioner context
951: . N - the number of subdomains for all processors
952: . is - the index sets that define the subdomains for all processors
953: (or NULL to ask PETSc to determine the subdomains)
954: - is_local - the index sets that define the local part of the subdomains for this processor
955: (or NULL to not provide this information)
957: Options Database Key:
958: To set the total number of subdomain blocks rather than specify the
959: index sets, use the option
960: . -pc_asm_blocks <blks> - Sets total blocks
962: Notes:
963: Currently you cannot use this to set the actual subdomains with the argument is or is_local.
965: By default the ASM preconditioner uses 1 block per processor.
967: These index sets cannot be destroyed until after completion of the
968: linear solves for which the ASM preconditioner is being used.
970: Use PCASMSetLocalSubdomains() to set local subdomains.
972: The IS numbering is in the parallel, global numbering of the vector for both is and is_local
974: Level: advanced
976: .seealso: PCASMSetLocalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
977: PCASMCreateSubdomains2D()
978: @*/
979: PetscErrorCode PCASMSetTotalSubdomains(PC pc,PetscInt N,IS is[],IS is_local[])
980: {
985: PetscTryMethod(pc,"PCASMSetTotalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,N,is,is_local));
986: return(0);
987: }
989: /*@
990: PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
991: additive Schwarz preconditioner. Either all or no processors in the
992: PC communicator must call this routine.
994: Logically Collective on pc
996: Input Parameters:
997: + pc - the preconditioner context
998: - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1)
1000: Options Database Key:
1001: . -pc_asm_overlap <ovl> - Sets overlap
1003: Notes:
1004: By default the ASM preconditioner uses 1 block per processor. To use
1005: multiple blocks per perocessor, see PCASMSetTotalSubdomains() and
1006: PCASMSetLocalSubdomains() (and the option -pc_asm_blocks <blks>).
1008: The overlap defaults to 1, so if one desires that no additional
1009: overlap be computed beyond what may have been set with a call to
1010: PCASMSetTotalSubdomains() or PCASMSetLocalSubdomains(), then ovl
1011: must be set to be 0. In particular, if one does not explicitly set
1012: the subdomains an application code, then all overlap would be computed
1013: internally by PETSc, and using an overlap of 0 would result in an ASM
1014: variant that is equivalent to the block Jacobi preconditioner.
1016: The default algorithm used by PETSc to increase overlap is fast, but not scalable,
1017: use the option -mat_increase_overlap_scalable when the problem and number of processes is large.
1019: Note that one can define initial index sets with any overlap via
1020: PCASMSetLocalSubdomains(); the routine
1021: PCASMSetOverlap() merely allows PETSc to extend that overlap further
1022: if desired.
1024: Level: intermediate
1026: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1027: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), MatIncreaseOverlap()
1028: @*/
1029: PetscErrorCode PCASMSetOverlap(PC pc,PetscInt ovl)
1030: {
1036: PetscTryMethod(pc,"PCASMSetOverlap_C",(PC,PetscInt),(pc,ovl));
1037: return(0);
1038: }
1040: /*@
1041: PCASMSetType - Sets the type of restriction and interpolation used
1042: for local problems in the additive Schwarz method.
1044: Logically Collective on pc
1046: Input Parameters:
1047: + pc - the preconditioner context
1048: - type - variant of ASM, one of
1049: .vb
1050: PC_ASM_BASIC - full interpolation and restriction
1051: PC_ASM_RESTRICT - full restriction, local processor interpolation
1052: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1053: PC_ASM_NONE - local processor restriction and interpolation
1054: .ve
1056: Options Database Key:
1057: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
1059: Notes:
1060: if the is_local arguments are passed to PCASMSetLocalSubdomains() then they are used when PC_ASM_RESTRICT has been selected
1061: to limit the local processor interpolation
1063: Level: intermediate
1065: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1066: PCASMCreateSubdomains2D(), PCASMType, PCASMSetLocalType(), PCASMGetLocalType()
1067: @*/
1068: PetscErrorCode PCASMSetType(PC pc,PCASMType type)
1069: {
1075: PetscTryMethod(pc,"PCASMSetType_C",(PC,PCASMType),(pc,type));
1076: return(0);
1077: }
1079: /*@
1080: PCASMGetType - Gets the type of restriction and interpolation used
1081: for local problems in the additive Schwarz method.
1083: Logically Collective on pc
1085: Input Parameter:
1086: . pc - the preconditioner context
1088: Output Parameter:
1089: . type - variant of ASM, one of
1091: .vb
1092: PC_ASM_BASIC - full interpolation and restriction
1093: PC_ASM_RESTRICT - full restriction, local processor interpolation
1094: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1095: PC_ASM_NONE - local processor restriction and interpolation
1096: .ve
1098: Options Database Key:
1099: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
1101: Level: intermediate
1103: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1104: PCASMCreateSubdomains2D(), PCASMType, PCASMSetType(), PCASMSetLocalType(), PCASMGetLocalType()
1105: @*/
1106: PetscErrorCode PCASMGetType(PC pc,PCASMType *type)
1107: {
1112: PetscUseMethod(pc,"PCASMGetType_C",(PC,PCASMType*),(pc,type));
1113: return(0);
1114: }
1116: /*@
1117: PCASMSetLocalType - Sets the type of composition used for local problems in the additive Schwarz method.
1119: Logically Collective on pc
1121: Input Parameters:
1122: + pc - the preconditioner context
1123: - type - type of composition, one of
1124: .vb
1125: PC_COMPOSITE_ADDITIVE - local additive combination
1126: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1127: .ve
1129: Options Database Key:
1130: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1132: Level: intermediate
1134: .seealso: PCASMSetType(), PCASMGetType(), PCASMGetLocalType(), PCASM, PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1135: @*/
1136: PetscErrorCode PCASMSetLocalType(PC pc, PCCompositeType type)
1137: {
1143: PetscTryMethod(pc, "PCASMSetLocalType_C", (PC, PCCompositeType), (pc, type));
1144: return(0);
1145: }
1147: /*@
1148: PCASMGetLocalType - Gets the type of composition used for local problems in the additive Schwarz method.
1150: Logically Collective on pc
1152: Input Parameter:
1153: . pc - the preconditioner context
1155: Output Parameter:
1156: . type - type of composition, one of
1157: .vb
1158: PC_COMPOSITE_ADDITIVE - local additive combination
1159: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1160: .ve
1162: Options Database Key:
1163: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1165: Level: intermediate
1167: .seealso: PCASMSetType(), PCASMGetType(), PCASMSetLocalType(), PCASMCreate(), PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1168: @*/
1169: PetscErrorCode PCASMGetLocalType(PC pc, PCCompositeType *type)
1170: {
1176: PetscUseMethod(pc, "PCASMGetLocalType_C", (PC, PCCompositeType *), (pc, type));
1177: return(0);
1178: }
1180: /*@
1181: PCASMSetSortIndices - Determines whether subdomain indices are sorted.
1183: Logically Collective on pc
1185: Input Parameters:
1186: + pc - the preconditioner context
1187: - doSort - sort the subdomain indices
1189: Level: intermediate
1191: .seealso: PCASMSetLocalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1192: PCASMCreateSubdomains2D()
1193: @*/
1194: PetscErrorCode PCASMSetSortIndices(PC pc,PetscBool doSort)
1195: {
1201: PetscTryMethod(pc,"PCASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));
1202: return(0);
1203: }
1205: /*@C
1206: PCASMGetSubKSP - Gets the local KSP contexts for all blocks on
1207: this processor.
1209: Collective on pc iff first_local is requested
1211: Input Parameter:
1212: . pc - the preconditioner context
1214: Output Parameters:
1215: + n_local - the number of blocks on this processor or NULL
1216: . first_local - the global number of the first block on this processor or NULL,
1217: all processors must request or all must pass NULL
1218: - ksp - the array of KSP contexts
1220: Note:
1221: After PCASMGetSubKSP() the array of KSPes is not to be freed.
1223: You must call KSPSetUp() before calling PCASMGetSubKSP().
1225: Fortran note:
1226: 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.
1228: Level: advanced
1230: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap(),
1231: PCASMCreateSubdomains2D(),
1232: @*/
1233: PetscErrorCode PCASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
1234: {
1239: PetscUseMethod(pc,"PCASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
1240: return(0);
1241: }
1243: /* -------------------------------------------------------------------------------------*/
1244: /*MC
1245: PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1246: its own KSP object.
1248: Options Database Keys:
1249: + -pc_asm_blocks <blks> - Sets total blocks
1250: . -pc_asm_overlap <ovl> - Sets overlap
1251: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type, default is restrict
1252: - -pc_asm_local_type [additive, multiplicative] - Sets ASM type, default is additive
1254: IMPORTANT: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you
1255: will get a different convergence rate due to the default option of -pc_asm_type restrict. Use
1256: -pc_asm_type basic to use the standard ASM.
1258: Notes:
1259: Each processor can have one or more blocks, but a block cannot be shared by more
1260: than one processor. Use PCGASM for subdomains shared by multiple processes. Defaults to one block per processor.
1262: To set options on the solvers for each block append -sub_ to all the KSP, and PC
1263: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
1265: To set the options on the solvers separate for each block call PCASMGetSubKSP()
1266: and set the options directly on the resulting KSP object (you can access its PC
1267: with KSPGetPC())
1269: Level: beginner
1271: References:
1272: + 1. - M Dryja, OB Widlund, An additive variant of the Schwarz alternating method for the case of many subregions
1273: Courant Institute, New York University Technical report
1274: - 2. - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations,
1275: Cambridge University Press.
1277: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
1278: PCBJACOBI, PCASMGetSubKSP(), PCASMSetLocalSubdomains(), PCASMType, PCASMGetType(), PCASMSetLocalType(), PCASMGetLocalType()
1279: PCASMSetTotalSubdomains(), PCSetModifySubMatrices(), PCASMSetOverlap(), PCASMSetType(), PCCompositeType
1281: M*/
1283: PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1284: {
1286: PC_ASM *osm;
1289: PetscNewLog(pc,&osm);
1291: osm->n = PETSC_DECIDE;
1292: osm->n_local = 0;
1293: osm->n_local_true = PETSC_DECIDE;
1294: osm->overlap = 1;
1295: osm->ksp = 0;
1296: osm->restriction = 0;
1297: osm->lprolongation = 0;
1298: osm->lrestriction = 0;
1299: osm->x = 0;
1300: osm->y = 0;
1301: osm->is = 0;
1302: osm->is_local = 0;
1303: osm->mat = 0;
1304: osm->pmat = 0;
1305: osm->type = PC_ASM_RESTRICT;
1306: osm->loctype = PC_COMPOSITE_ADDITIVE;
1307: osm->same_local_solves = PETSC_TRUE;
1308: osm->sort_indices = PETSC_TRUE;
1309: osm->dm_subdomains = PETSC_FALSE;
1310: osm->sub_mat_type = NULL;
1312: pc->data = (void*)osm;
1313: pc->ops->apply = PCApply_ASM;
1314: pc->ops->applytranspose = PCApplyTranspose_ASM;
1315: pc->ops->setup = PCSetUp_ASM;
1316: pc->ops->reset = PCReset_ASM;
1317: pc->ops->destroy = PCDestroy_ASM;
1318: pc->ops->setfromoptions = PCSetFromOptions_ASM;
1319: pc->ops->setuponblocks = PCSetUpOnBlocks_ASM;
1320: pc->ops->view = PCView_ASM;
1321: pc->ops->applyrichardson = 0;
1323: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",PCASMSetLocalSubdomains_ASM);
1324: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",PCASMSetTotalSubdomains_ASM);
1325: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetOverlap_C",PCASMSetOverlap_ASM);
1326: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetType_C",PCASMSetType_ASM);
1327: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetType_C",PCASMGetType_ASM);
1328: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalType_C",PCASMSetLocalType_ASM);
1329: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetLocalType_C",PCASMGetLocalType_ASM);
1330: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSortIndices_C",PCASMSetSortIndices_ASM);
1331: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubKSP_C",PCASMGetSubKSP_ASM);
1332: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubMatType_C",PCASMGetSubMatType_ASM);
1333: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSubMatType_C",PCASMSetSubMatType_ASM);
1334: return(0);
1335: }
1337: /*@C
1338: PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1339: preconditioner for a any problem on a general grid.
1341: Collective
1343: Input Parameters:
1344: + A - The global matrix operator
1345: - n - the number of local blocks
1347: Output Parameters:
1348: . outis - the array of index sets defining the subdomains
1350: Level: advanced
1352: Note: this generates nonoverlapping subdomains; the PCASM will generate the overlap
1353: from these if you use PCASMSetLocalSubdomains()
1355: In the Fortran version you must provide the array outis[] already allocated of length n.
1357: .seealso: PCASMSetLocalSubdomains(), PCASMDestroySubdomains()
1358: @*/
1359: PetscErrorCode PCASMCreateSubdomains(Mat A, PetscInt n, IS* outis[])
1360: {
1361: MatPartitioning mpart;
1362: const char *prefix;
1363: PetscInt i,j,rstart,rend,bs;
1364: PetscBool hasop, isbaij = PETSC_FALSE,foundpart = PETSC_FALSE;
1365: Mat Ad = NULL, adj;
1366: IS ispart,isnumb,*is;
1367: PetscErrorCode ierr;
1372: if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of local blocks must be > 0, n = %D",n);
1374: /* Get prefix, row distribution, and block size */
1375: MatGetOptionsPrefix(A,&prefix);
1376: MatGetOwnershipRange(A,&rstart,&rend);
1377: MatGetBlockSize(A,&bs);
1378: if (rstart/bs*bs != rstart || rend/bs*bs != rend) SETERRQ3(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"bad row distribution [%D,%D) for matrix block size %D",rstart,rend,bs);
1380: /* Get diagonal block from matrix if possible */
1381: MatHasOperation(A,MATOP_GET_DIAGONAL_BLOCK,&hasop);
1382: if (hasop) {
1383: MatGetDiagonalBlock(A,&Ad);
1384: }
1385: if (Ad) {
1386: PetscObjectBaseTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij);
1387: if (!isbaij) {PetscObjectBaseTypeCompare((PetscObject)Ad,MATSEQSBAIJ,&isbaij);}
1388: }
1389: if (Ad && n > 1) {
1390: PetscBool match,done;
1391: /* Try to setup a good matrix partitioning if available */
1392: MatPartitioningCreate(PETSC_COMM_SELF,&mpart);
1393: PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);
1394: MatPartitioningSetFromOptions(mpart);
1395: PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGCURRENT,&match);
1396: if (!match) {
1397: PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGSQUARE,&match);
1398: }
1399: if (!match) { /* assume a "good" partitioner is available */
1400: PetscInt na;
1401: const PetscInt *ia,*ja;
1402: MatGetRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1403: if (done) {
1404: /* Build adjacency matrix by hand. Unfortunately a call to
1405: MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1406: remove the block-aij structure and we cannot expect
1407: MatPartitioning to split vertices as we need */
1408: PetscInt i,j,len,nnz,cnt,*iia=0,*jja=0;
1409: const PetscInt *row;
1410: nnz = 0;
1411: for (i=0; i<na; i++) { /* count number of nonzeros */
1412: len = ia[i+1] - ia[i];
1413: row = ja + ia[i];
1414: for (j=0; j<len; j++) {
1415: if (row[j] == i) { /* don't count diagonal */
1416: len--; break;
1417: }
1418: }
1419: nnz += len;
1420: }
1421: PetscMalloc1(na+1,&iia);
1422: PetscMalloc1(nnz,&jja);
1423: nnz = 0;
1424: iia[0] = 0;
1425: for (i=0; i<na; i++) { /* fill adjacency */
1426: cnt = 0;
1427: len = ia[i+1] - ia[i];
1428: row = ja + ia[i];
1429: for (j=0; j<len; j++) {
1430: if (row[j] != i) { /* if not diagonal */
1431: jja[nnz+cnt++] = row[j];
1432: }
1433: }
1434: nnz += cnt;
1435: iia[i+1] = nnz;
1436: }
1437: /* Partitioning of the adjacency matrix */
1438: MatCreateMPIAdj(PETSC_COMM_SELF,na,na,iia,jja,NULL,&adj);
1439: MatPartitioningSetAdjacency(mpart,adj);
1440: MatPartitioningSetNParts(mpart,n);
1441: MatPartitioningApply(mpart,&ispart);
1442: ISPartitioningToNumbering(ispart,&isnumb);
1443: MatDestroy(&adj);
1444: foundpart = PETSC_TRUE;
1445: }
1446: MatRestoreRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1447: }
1448: MatPartitioningDestroy(&mpart);
1449: }
1451: PetscMalloc1(n,&is);
1452: *outis = is;
1454: if (!foundpart) {
1456: /* Partitioning by contiguous chunks of rows */
1458: PetscInt mbs = (rend-rstart)/bs;
1459: PetscInt start = rstart;
1460: for (i=0; i<n; i++) {
1461: PetscInt count = (mbs/n + ((mbs % n) > i)) * bs;
1462: ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);
1463: start += count;
1464: }
1466: } else {
1468: /* Partitioning by adjacency of diagonal block */
1470: const PetscInt *numbering;
1471: PetscInt *count,nidx,*indices,*newidx,start=0;
1472: /* Get node count in each partition */
1473: PetscMalloc1(n,&count);
1474: ISPartitioningCount(ispart,n,count);
1475: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1476: for (i=0; i<n; i++) count[i] *= bs;
1477: }
1478: /* Build indices from node numbering */
1479: ISGetLocalSize(isnumb,&nidx);
1480: PetscMalloc1(nidx,&indices);
1481: for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */
1482: ISGetIndices(isnumb,&numbering);
1483: PetscSortIntWithPermutation(nidx,numbering,indices);
1484: ISRestoreIndices(isnumb,&numbering);
1485: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1486: PetscMalloc1(nidx*bs,&newidx);
1487: for (i=0; i<nidx; i++) {
1488: for (j=0; j<bs; j++) newidx[i*bs+j] = indices[i]*bs + j;
1489: }
1490: PetscFree(indices);
1491: nidx *= bs;
1492: indices = newidx;
1493: }
1494: /* Shift to get global indices */
1495: for (i=0; i<nidx; i++) indices[i] += rstart;
1497: /* Build the index sets for each block */
1498: for (i=0; i<n; i++) {
1499: ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);
1500: ISSort(is[i]);
1501: start += count[i];
1502: }
1504: PetscFree(count);
1505: PetscFree(indices);
1506: ISDestroy(&isnumb);
1507: ISDestroy(&ispart);
1509: }
1510: return(0);
1511: }
1513: /*@C
1514: PCASMDestroySubdomains - Destroys the index sets created with
1515: PCASMCreateSubdomains(). Should be called after setting subdomains
1516: with PCASMSetLocalSubdomains().
1518: Collective
1520: Input Parameters:
1521: + n - the number of index sets
1522: . is - the array of index sets
1523: - is_local - the array of local index sets, can be NULL
1525: Level: advanced
1527: .seealso: PCASMCreateSubdomains(), PCASMSetLocalSubdomains()
1528: @*/
1529: PetscErrorCode PCASMDestroySubdomains(PetscInt n, IS is[], IS is_local[])
1530: {
1531: PetscInt i;
1535: if (n <= 0) return(0);
1536: if (is) {
1538: for (i=0; i<n; i++) { ISDestroy(&is[i]); }
1539: PetscFree(is);
1540: }
1541: if (is_local) {
1543: for (i=0; i<n; i++) { ISDestroy(&is_local[i]); }
1544: PetscFree(is_local);
1545: }
1546: return(0);
1547: }
1549: /*@
1550: PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1551: preconditioner for a two-dimensional problem on a regular grid.
1553: Not Collective
1555: Input Parameters:
1556: + m, n - the number of mesh points in the x and y directions
1557: . M, N - the number of subdomains in the x and y directions
1558: . dof - degrees of freedom per node
1559: - overlap - overlap in mesh lines
1561: Output Parameters:
1562: + Nsub - the number of subdomains created
1563: . is - array of index sets defining overlapping (if overlap > 0) subdomains
1564: - is_local - array of index sets defining non-overlapping subdomains
1566: Note:
1567: Presently PCAMSCreateSubdomains2d() is valid only for sequential
1568: preconditioners. More general related routines are
1569: PCASMSetTotalSubdomains() and PCASMSetLocalSubdomains().
1571: Level: advanced
1573: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1574: PCASMSetOverlap()
1575: @*/
1576: PetscErrorCode PCASMCreateSubdomains2D(PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt dof,PetscInt overlap,PetscInt *Nsub,IS **is,IS **is_local)
1577: {
1578: PetscInt i,j,height,width,ystart,xstart,yleft,yright,xleft,xright,loc_outer;
1580: PetscInt nidx,*idx,loc,ii,jj,count;
1583: if (dof != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP," ");
1585: *Nsub = N*M;
1586: PetscMalloc1(*Nsub,is);
1587: PetscMalloc1(*Nsub,is_local);
1588: ystart = 0;
1589: loc_outer = 0;
1590: for (i=0; i<N; i++) {
1591: height = n/N + ((n % N) > i); /* height of subdomain */
1592: if (height < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many N subdomains for mesh dimension n");
1593: yleft = ystart - overlap; if (yleft < 0) yleft = 0;
1594: yright = ystart + height + overlap; if (yright > n) yright = n;
1595: xstart = 0;
1596: for (j=0; j<M; j++) {
1597: width = m/M + ((m % M) > j); /* width of subdomain */
1598: if (width < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many M subdomains for mesh dimension m");
1599: xleft = xstart - overlap; if (xleft < 0) xleft = 0;
1600: xright = xstart + width + overlap; if (xright > m) xright = m;
1601: nidx = (xright - xleft)*(yright - yleft);
1602: PetscMalloc1(nidx,&idx);
1603: loc = 0;
1604: for (ii=yleft; ii<yright; ii++) {
1605: count = m*ii + xleft;
1606: for (jj=xleft; jj<xright; jj++) idx[loc++] = count++;
1607: }
1608: ISCreateGeneral(PETSC_COMM_SELF,nidx,idx,PETSC_COPY_VALUES,(*is)+loc_outer);
1609: if (overlap == 0) {
1610: PetscObjectReference((PetscObject)(*is)[loc_outer]);
1612: (*is_local)[loc_outer] = (*is)[loc_outer];
1613: } else {
1614: for (loc=0,ii=ystart; ii<ystart+height; ii++) {
1615: for (jj=xstart; jj<xstart+width; jj++) {
1616: idx[loc++] = m*ii + jj;
1617: }
1618: }
1619: ISCreateGeneral(PETSC_COMM_SELF,loc,idx,PETSC_COPY_VALUES,*is_local+loc_outer);
1620: }
1621: PetscFree(idx);
1622: xstart += width;
1623: loc_outer++;
1624: }
1625: ystart += height;
1626: }
1627: for (i=0; i<*Nsub; i++) { ISSort((*is)[i]); }
1628: return(0);
1629: }
1631: /*@C
1632: PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
1633: only) for the additive Schwarz preconditioner.
1635: Not Collective
1637: Input Parameter:
1638: . pc - the preconditioner context
1640: Output Parameters:
1641: + n - the number of subdomains for this processor (default value = 1)
1642: . is - the index sets that define the subdomains for this processor
1643: - is_local - the index sets that define the local part of the subdomains for this processor (can be NULL)
1646: Notes:
1647: The IS numbering is in the parallel, global numbering of the vector.
1649: Level: advanced
1651: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1652: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubmatrices()
1653: @*/
1654: PetscErrorCode PCASMGetLocalSubdomains(PC pc,PetscInt *n,IS *is[],IS *is_local[])
1655: {
1656: PC_ASM *osm = (PC_ASM*)pc->data;
1658: PetscBool match;
1664: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1665: if (!match) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PC is not a PCASM");
1666: if (n) *n = osm->n_local_true;
1667: if (is) *is = osm->is;
1668: if (is_local) *is_local = osm->is_local;
1669: return(0);
1670: }
1672: /*@C
1673: PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1674: only) for the additive Schwarz preconditioner.
1676: Not Collective
1678: Input Parameter:
1679: . pc - the preconditioner context
1681: Output Parameters:
1682: + n - the number of matrices for this processor (default value = 1)
1683: - mat - the matrices
1685: Level: advanced
1687: Notes:
1688: Call after PCSetUp() (or KSPSetUp()) but before PCApply() and before PCSetUpOnBlocks())
1690: Usually one would use PCSetModifySubMatrices() to change the submatrices in building the preconditioner.
1692: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1693: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains(), PCSetModifySubMatrices()
1694: @*/
1695: PetscErrorCode PCASMGetLocalSubmatrices(PC pc,PetscInt *n,Mat *mat[])
1696: {
1697: PC_ASM *osm;
1699: PetscBool match;
1705: if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUp() or PCSetUp().");
1706: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1707: if (!match) {
1708: if (n) *n = 0;
1709: if (mat) *mat = NULL;
1710: } else {
1711: osm = (PC_ASM*)pc->data;
1712: if (n) *n = osm->n_local_true;
1713: if (mat) *mat = osm->pmat;
1714: }
1715: return(0);
1716: }
1718: /*@
1719: PCASMSetDMSubdomains - Indicates whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1721: Logically Collective
1723: Input Parameter:
1724: + pc - the preconditioner
1725: - flg - boolean indicating whether to use subdomains defined by the DM
1727: Options Database Key:
1728: . -pc_asm_dm_subdomains
1730: Level: intermediate
1732: Notes:
1733: PCASMSetTotalSubdomains() and PCASMSetOverlap() take precedence over PCASMSetDMSubdomains(),
1734: so setting either of the first two effectively turns the latter off.
1736: .seealso: PCASMGetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1737: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1738: @*/
1739: PetscErrorCode PCASMSetDMSubdomains(PC pc,PetscBool flg)
1740: {
1741: PC_ASM *osm = (PC_ASM*)pc->data;
1743: PetscBool match;
1748: if (pc->setupcalled) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for a setup PC.");
1749: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1750: if (match) {
1751: osm->dm_subdomains = flg;
1752: }
1753: return(0);
1754: }
1756: /*@
1757: PCASMGetDMSubdomains - Returns flag indicating whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1758: Not Collective
1760: Input Parameter:
1761: . pc - the preconditioner
1763: Output Parameter:
1764: . flg - boolean indicating whether to use subdomains defined by the DM
1766: Level: intermediate
1768: .seealso: PCASMSetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1769: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1770: @*/
1771: PetscErrorCode PCASMGetDMSubdomains(PC pc,PetscBool* flg)
1772: {
1773: PC_ASM *osm = (PC_ASM*)pc->data;
1775: PetscBool match;
1780: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1781: if (match) {
1782: if (flg) *flg = osm->dm_subdomains;
1783: }
1784: return(0);
1785: }
1787: /*@
1788: PCASMGetSubMatType - Gets the matrix type used for ASM subsolves, as a string.
1790: Not Collective
1792: Input Parameter:
1793: . pc - the PC
1795: Output Parameter:
1796: . -pc_asm_sub_mat_type - name of matrix type
1798: Level: advanced
1800: .seealso: PCASMSetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1801: @*/
1802: PetscErrorCode PCASMGetSubMatType(PC pc,MatType *sub_mat_type){
1805: PetscTryMethod(pc,"PCASMGetSubMatType_C",(PC,MatType*),(pc,sub_mat_type));
1806: return(0);
1807: }
1809: /*@
1810: PCASMSetSubMatType - Set the type of matrix used for ASM subsolves
1812: Collective on Mat
1814: Input Parameters:
1815: + pc - the PC object
1816: - sub_mat_type - matrix type
1818: Options Database Key:
1819: . -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.
1821: Notes:
1822: See "${PETSC_DIR}/include/petscmat.h" for available types
1824: Level: advanced
1826: .seealso: PCASMGetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1827: @*/
1828: PetscErrorCode PCASMSetSubMatType(PC pc,MatType sub_mat_type)
1829: {
1832: PetscTryMethod(pc,"PCASMSetSubMatType_C",(PC,MatType),(pc,sub_mat_type));
1833: return(0);
1834: }