Actual source code: asm.c
petsc-3.9.4 2018-09-11
2: /*
3: This file defines an additive Schwarz preconditioner for any Mat implementation.
5: Note that each processor may have any number of subdomains. But in order to
6: deal easily with the VecScatter(), we treat each processor as if it has the
7: same number of subdomains.
9: n - total number of true subdomains on all processors
10: n_local_true - actual number of subdomains on this processor
11: n_local = maximum over all processors of n_local_true
12: */
13: #include <petsc/private/pcimpl.h>
14: #include <petscdm.h>
16: typedef struct {
17: PetscInt n, n_local, n_local_true;
18: PetscInt overlap; /* overlap requested by user */
19: KSP *ksp; /* linear solvers for each block */
20: VecScatter restriction; /* mapping from global to overlapping (process) subdomain*/
21: VecScatter *lrestriction; /* mapping from subregion to overlapping (process) subdomain */
22: VecScatter *lprolongation; /* mapping from non-overlapping subregion to overlapping (process) subdomain; used for restrict additive version of algorithms */
23: Vec lx, ly; /* work vectors */
24: Vec *x,*y; /* work vectors */
25: IS lis; /* index set that defines each overlapping multiplicative (process) subdomain */
26: IS *is; /* index set that defines each overlapping subdomain */
27: IS *is_local; /* index set that defines each non-overlapping subdomain, may be NULL */
28: Mat *mat,*pmat; /* mat is not currently used */
29: PCASMType type; /* use reduced interpolation, restriction or both */
30: PetscBool type_set; /* if user set this value (so won't change it for symmetric problems) */
31: PetscBool same_local_solves; /* flag indicating whether all local solvers are same */
32: PetscBool sort_indices; /* flag to sort subdomain indices */
33: PetscBool dm_subdomains; /* whether DM is allowed to define subdomains */
34: PCCompositeType loctype; /* the type of composition for local solves */
35: MatType sub_mat_type; /* the type of Mat used for subdomain solves (can be MATSAME or NULL) */
36: /* For multiplicative solve */
37: Mat *lmats; /* submatrices for overlapping multiplicative (process) subdomain */
38: } PC_ASM;
40: static PetscErrorCode PCView_ASM(PC pc,PetscViewer viewer)
41: {
42: PC_ASM *osm = (PC_ASM*)pc->data;
44: PetscMPIInt rank;
45: PetscInt i,bsz;
46: PetscBool iascii,isstring;
47: PetscViewer sviewer;
50: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
51: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
52: if (iascii) {
53: char overlaps[256] = "user-defined overlap",blocks[256] = "total subdomain blocks not yet set";
54: if (osm->overlap >= 0) {PetscSNPrintf(overlaps,sizeof(overlaps),"amount of overlap = %D",osm->overlap);}
55: if (osm->n > 0) {PetscSNPrintf(blocks,sizeof(blocks),"total subdomain blocks = %D",osm->n);}
56: PetscViewerASCIIPrintf(viewer," %s, %s\n",blocks,overlaps);
57: PetscViewerASCIIPrintf(viewer," restriction/interpolation type - %s\n",PCASMTypes[osm->type]);
58: if (osm->dm_subdomains) {PetscViewerASCIIPrintf(viewer," Additive Schwarz: using DM to define subdomains\n");}
59: if (osm->loctype != PC_COMPOSITE_ADDITIVE) {PetscViewerASCIIPrintf(viewer," Additive Schwarz: local solve composition type - %s\n",PCCompositeTypes[osm->loctype]);}
60: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
61: if (osm->same_local_solves) {
62: if (osm->ksp) {
63: PetscViewerASCIIPrintf(viewer," Local solve is same for all blocks, in the following KSP and PC objects:\n");
64: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
65: if (!rank) {
66: PetscViewerASCIIPushTab(viewer);
67: KSPView(osm->ksp[0],sviewer);
68: PetscViewerASCIIPopTab(viewer);
69: }
70: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
71: }
72: } else {
73: PetscViewerASCIIPushSynchronized(viewer);
74: PetscViewerASCIISynchronizedPrintf(viewer," [%d] number of local blocks = %D\n",(int)rank,osm->n_local_true);
75: PetscViewerFlush(viewer);
76: PetscViewerASCIIPrintf(viewer," Local solve info for each block is in the following KSP and PC objects:\n");
77: PetscViewerASCIIPushTab(viewer);
78: PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
79: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
80: for (i=0; i<osm->n_local_true; i++) {
81: ISGetLocalSize(osm->is[i],&bsz);
82: PetscViewerASCIISynchronizedPrintf(sviewer,"[%d] local block number %D, size = %D\n",(int)rank,i,bsz);
83: KSPView(osm->ksp[i],sviewer);
84: PetscViewerASCIISynchronizedPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
85: }
86: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
87: PetscViewerASCIIPopTab(viewer);
88: PetscViewerFlush(viewer);
89: PetscViewerASCIIPopSynchronized(viewer);
90: }
91: } else if (isstring) {
92: PetscViewerStringSPrintf(viewer," blocks=%D, overlap=%D, type=%s",osm->n,osm->overlap,PCASMTypes[osm->type]);
93: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
94: if (osm->ksp) {KSPView(osm->ksp[0],sviewer);}
95: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
96: }
97: return(0);
98: }
100: static PetscErrorCode PCASMPrintSubdomains(PC pc)
101: {
102: PC_ASM *osm = (PC_ASM*)pc->data;
103: const char *prefix;
104: char fname[PETSC_MAX_PATH_LEN+1];
105: PetscViewer viewer, sviewer;
106: char *s;
107: PetscInt i,j,nidx;
108: const PetscInt *idx;
109: PetscMPIInt rank, size;
113: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
114: MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
115: PCGetOptionsPrefix(pc,&prefix);
116: PetscOptionsGetString(NULL,prefix,"-pc_asm_print_subdomains",fname,PETSC_MAX_PATH_LEN,NULL);
117: if (fname[0] == 0) { PetscStrcpy(fname,"stdout"); };
118: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),fname,&viewer);
119: for (i=0; i<osm->n_local; i++) {
120: if (i < osm->n_local_true) {
121: ISGetLocalSize(osm->is[i],&nidx);
122: ISGetIndices(osm->is[i],&idx);
123: /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
124: PetscMalloc1(16*(nidx+1)+512, &s);
125: PetscViewerStringOpen(PETSC_COMM_SELF, s, 16*(nidx+1)+512, &sviewer);
126: PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D with overlap:\n", rank, size, i);
127: for (j=0; j<nidx; j++) {
128: PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
129: }
130: ISRestoreIndices(osm->is[i],&idx);
131: PetscViewerStringSPrintf(sviewer,"\n");
132: PetscViewerDestroy(&sviewer);
133: PetscViewerASCIIPushSynchronized(viewer);
134: PetscViewerASCIISynchronizedPrintf(viewer, s);
135: PetscViewerFlush(viewer);
136: PetscViewerASCIIPopSynchronized(viewer);
137: PetscFree(s);
138: if (osm->is_local) {
139: /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
140: PetscMalloc1(16*(nidx+1)+512, &s);
141: PetscViewerStringOpen(PETSC_COMM_SELF, s, 16*(nidx+1)+512, &sviewer);
142: PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D without overlap:\n", rank, size, i);
143: ISGetLocalSize(osm->is_local[i],&nidx);
144: ISGetIndices(osm->is_local[i],&idx);
145: for (j=0; j<nidx; j++) {
146: PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
147: }
148: ISRestoreIndices(osm->is_local[i],&idx);
149: PetscViewerStringSPrintf(sviewer,"\n");
150: PetscViewerDestroy(&sviewer);
151: PetscViewerASCIIPushSynchronized(viewer);
152: PetscViewerASCIISynchronizedPrintf(viewer, s);
153: PetscViewerFlush(viewer);
154: PetscViewerASCIIPopSynchronized(viewer);
155: PetscFree(s);
156: }
157: } else {
158: /* Participate in collective viewer calls. */
159: PetscViewerASCIIPushSynchronized(viewer);
160: PetscViewerFlush(viewer);
161: PetscViewerASCIIPopSynchronized(viewer);
162: /* Assume either all ranks have is_local or none do. */
163: if (osm->is_local) {
164: PetscViewerASCIIPushSynchronized(viewer);
165: PetscViewerFlush(viewer);
166: PetscViewerASCIIPopSynchronized(viewer);
167: }
168: }
169: }
170: PetscViewerFlush(viewer);
171: PetscViewerDestroy(&viewer);
172: return(0);
173: }
175: static PetscErrorCode PCSetUp_ASM(PC pc)
176: {
177: PC_ASM *osm = (PC_ASM*)pc->data;
179: PetscBool symset,flg;
180: PetscInt i,m,m_local;
181: MatReuse scall = MAT_REUSE_MATRIX;
182: IS isl;
183: KSP ksp;
184: PC subpc;
185: const char *prefix,*pprefix;
186: Vec vec;
187: DM *domain_dm = NULL;
190: if (!pc->setupcalled) {
191: PetscInt m;
193: if (!osm->type_set) {
194: MatIsSymmetricKnown(pc->pmat,&symset,&flg);
195: if (symset && flg) osm->type = PC_ASM_BASIC;
196: }
198: /* Note: if subdomains have been set either via PCASMSetTotalSubdomains() or via PCASMSetLocalSubdomains(), osm->n_local_true will not be PETSC_DECIDE */
199: if (osm->n_local_true == PETSC_DECIDE) {
200: /* no subdomains given */
201: /* try pc->dm first, if allowed */
202: if (osm->dm_subdomains && pc->dm) {
203: PetscInt num_domains, d;
204: char **domain_names;
205: IS *inner_domain_is, *outer_domain_is;
206: DMCreateDomainDecomposition(pc->dm, &num_domains, &domain_names, &inner_domain_is, &outer_domain_is, &domain_dm);
207: osm->overlap = -1; /* We do not want to increase the overlap of the IS.
208: A future improvement of this code might allow one to use
209: DM-defined subdomains and also increase the overlap,
210: but that is not currently supported */
211: if (num_domains) {
212: PCASMSetLocalSubdomains(pc, num_domains, outer_domain_is, inner_domain_is);
213: }
214: for (d = 0; d < num_domains; ++d) {
215: if (domain_names) {PetscFree(domain_names[d]);}
216: if (inner_domain_is) {ISDestroy(&inner_domain_is[d]);}
217: if (outer_domain_is) {ISDestroy(&outer_domain_is[d]);}
218: }
219: PetscFree(domain_names);
220: PetscFree(inner_domain_is);
221: PetscFree(outer_domain_is);
222: }
223: if (osm->n_local_true == PETSC_DECIDE) {
224: /* still no subdomains; use one subdomain per processor */
225: osm->n_local_true = 1;
226: }
227: }
228: { /* determine the global and max number of subdomains */
229: struct {PetscInt max,sum;} inwork,outwork;
230: PetscMPIInt size;
232: inwork.max = osm->n_local_true;
233: inwork.sum = osm->n_local_true;
234: MPIU_Allreduce(&inwork,&outwork,1,MPIU_2INT,MPIU_MAXSUM_OP,PetscObjectComm((PetscObject)pc));
235: osm->n_local = outwork.max;
236: osm->n = outwork.sum;
238: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
239: if (outwork.max == 1 && outwork.sum == size) {
240: /* osm->n_local_true = 1 on all processes, set this option may enable use of optimized MatCreateSubMatrices() implementation */
241: MatSetOption(pc->pmat,MAT_SUBMAT_SINGLEIS,PETSC_TRUE);
242: }
243: }
244: if (!osm->is) { /* create the index sets */
245: PCASMCreateSubdomains(pc->pmat,osm->n_local_true,&osm->is);
246: }
247: if (osm->n_local_true > 1 && !osm->is_local) {
248: PetscMalloc1(osm->n_local_true,&osm->is_local);
249: for (i=0; i<osm->n_local_true; i++) {
250: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
251: ISDuplicate(osm->is[i],&osm->is_local[i]);
252: ISCopy(osm->is[i],osm->is_local[i]);
253: } else {
254: PetscObjectReference((PetscObject)osm->is[i]);
255: osm->is_local[i] = osm->is[i];
256: }
257: }
258: }
259: PCGetOptionsPrefix(pc,&prefix);
260: flg = PETSC_FALSE;
261: PetscOptionsGetBool(NULL,prefix,"-pc_asm_print_subdomains",&flg,NULL);
262: if (flg) { PCASMPrintSubdomains(pc); }
264: if (osm->overlap > 0) {
265: /* Extend the "overlapping" regions by a number of steps */
266: MatIncreaseOverlap(pc->pmat,osm->n_local_true,osm->is,osm->overlap);
267: }
268: if (osm->sort_indices) {
269: for (i=0; i<osm->n_local_true; i++) {
270: ISSort(osm->is[i]);
271: if (osm->is_local) {
272: ISSort(osm->is_local[i]);
273: }
274: }
275: }
277: if (!osm->ksp) {
278: /* Create the local solvers */
279: PetscMalloc1(osm->n_local_true,&osm->ksp);
280: if (domain_dm) {
281: PetscInfo(pc,"Setting up ASM subproblems using the embedded DM\n");
282: }
283: for (i=0; i<osm->n_local_true; i++) {
284: KSPCreate(PETSC_COMM_SELF,&ksp);
285: KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
286: PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
287: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
288: KSPSetType(ksp,KSPPREONLY);
289: KSPGetPC(ksp,&subpc);
290: PCGetOptionsPrefix(pc,&prefix);
291: KSPSetOptionsPrefix(ksp,prefix);
292: KSPAppendOptionsPrefix(ksp,"sub_");
293: if (domain_dm) {
294: KSPSetDM(ksp, domain_dm[i]);
295: KSPSetDMActive(ksp, PETSC_FALSE);
296: DMDestroy(&domain_dm[i]);
297: }
298: osm->ksp[i] = ksp;
299: }
300: if (domain_dm) {
301: PetscFree(domain_dm);
302: }
303: }
304:
305: ISConcatenate(PETSC_COMM_SELF, osm->n_local_true, osm->is, &osm->lis);
306: ISSortRemoveDups(osm->lis);
307: ISGetLocalSize(osm->lis, &m);
308: VecCreateSeq(PETSC_COMM_SELF, m, &osm->lx);
309: VecDuplicate(osm->lx, &osm->ly);
310:
311: scall = MAT_INITIAL_MATRIX;
312: } else {
313: /*
314: Destroy the blocks from the previous iteration
315: */
316: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
317: MatDestroyMatrices(osm->n_local_true,&osm->pmat);
318: scall = MAT_INITIAL_MATRIX;
319: }
320: }
322: /*
323: Extract out the submatrices
324: */
325: MatCreateSubMatrices(pc->pmat,osm->n_local_true,osm->is,osm->is,scall,&osm->pmat);
326: if (scall == MAT_INITIAL_MATRIX) {
327: PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
328: for (i=0; i<osm->n_local_true; i++) {
329: PetscLogObjectParent((PetscObject)pc,(PetscObject)osm->pmat[i]);
330: PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);
331: }
332: }
334: /* Convert the types of the submatrices (if needbe) */
335: if (osm->sub_mat_type) {
336: for (i=0; i<osm->n_local_true; i++) {
337: MatConvert(osm->pmat[i],osm->sub_mat_type,MAT_INPLACE_MATRIX,&(osm->pmat[i]));
338: }
339: }
341: if(!pc->setupcalled){
342: /* Create the local work vectors (from the local matrices) and scatter contexts */
343: MatCreateVecs(pc->pmat,&vec,0);
344:
345: 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()");
346: if (osm->is_local && osm->type == PC_ASM_RESTRICT && osm->loctype == PC_COMPOSITE_ADDITIVE) {
347: PetscMalloc1(osm->n_local_true,&osm->lprolongation);
348: }
349: PetscMalloc1(osm->n_local_true,&osm->lrestriction);
350: PetscMalloc1(osm->n_local_true,&osm->x);
351: PetscMalloc1(osm->n_local_true,&osm->y);
352:
353: ISGetLocalSize(osm->lis,&m);
354: ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
355: VecScatterCreate(vec,osm->lis,osm->lx,isl,&osm->restriction);
356: ISDestroy(&isl);
357:
358:
359: for (i=0; i<osm->n_local_true; ++i) {
360: ISLocalToGlobalMapping ltog;
361: IS isll;
362: const PetscInt *idx_is;
363: PetscInt *idx_lis,nout;
365: ISGetLocalSize(osm->is[i],&m);
366: MatCreateVecs(osm->pmat[i],&osm->x[i],NULL);
367: VecDuplicate(osm->x[i],&osm->y[i]);
368:
369: /* generate a scatter from ly to y[i] picking all the overlapping is[i] entries */
370: ISLocalToGlobalMappingCreateIS(osm->lis,<og);
371: ISGetLocalSize(osm->is[i],&m);
372: ISGetIndices(osm->is[i], &idx_is);
373: PetscMalloc1(m,&idx_lis);
374: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m,idx_is,&nout,idx_lis);
375: if (nout != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is not a subset of lis");
376: ISRestoreIndices(osm->is[i], &idx_is);
377: ISCreateGeneral(PETSC_COMM_SELF,m,idx_lis,PETSC_OWN_POINTER,&isll);
378: ISLocalToGlobalMappingDestroy(<og);
379: ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
380: VecScatterCreate(osm->ly,isll,osm->y[i],isl,&osm->lrestriction[i]);
381: ISDestroy(&isll);
382: ISDestroy(&isl);
383: if (osm->lprolongation) { /* generate a scatter from y[i] to ly picking only the the non-overalapping is_local[i] entries */
384: ISLocalToGlobalMapping ltog;
385: IS isll,isll_local;
386: const PetscInt *idx_local;
387: PetscInt *idx1, *idx2, nout;
388:
389: ISGetLocalSize(osm->is_local[i],&m_local);
390: ISGetIndices(osm->is_local[i], &idx_local);
391:
392: ISLocalToGlobalMappingCreateIS(osm->is[i],<og);
393: PetscMalloc1(m_local,&idx1);
394: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx1);
395: ISLocalToGlobalMappingDestroy(<og);
396: if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of is");
397: ISCreateGeneral(PETSC_COMM_SELF,m_local,idx1,PETSC_OWN_POINTER,&isll);
398:
399: ISLocalToGlobalMappingCreateIS(osm->lis,<og);
400: PetscMalloc1(m_local,&idx2);
401: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx2);
402: ISLocalToGlobalMappingDestroy(<og);
403: if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of lis");
404: ISCreateGeneral(PETSC_COMM_SELF,m_local,idx2,PETSC_OWN_POINTER,&isll_local);
405:
406: ISRestoreIndices(osm->is_local[i], &idx_local);
407: VecScatterCreate(osm->y[i],isll,osm->ly,isll_local,&osm->lprolongation[i]);
408:
409: ISDestroy(&isll);
410: ISDestroy(&isll_local);
411: }
412: }
413: VecDestroy(&vec);
414: }
416: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
417: IS *cis;
418: PetscInt c;
420: PetscMalloc1(osm->n_local_true, &cis);
421: for (c = 0; c < osm->n_local_true; ++c) cis[c] = osm->lis;
422: MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, cis, scall, &osm->lmats);
423: PetscFree(cis);
424: }
426: /* Return control to the user so that the submatrices can be modified (e.g., to apply
427: different boundary conditions for the submatrices than for the global problem) */
428: PCModifySubMatrices(pc,osm->n_local_true,osm->is,osm->is,osm->pmat,pc->modifysubmatricesP);
430: /*
431: Loop over subdomains putting them into local ksp
432: */
433: for (i=0; i<osm->n_local_true; i++) {
434: KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i]);
435: if (!pc->setupcalled) {
436: KSPSetFromOptions(osm->ksp[i]);
437: }
438: }
439: return(0);
440: }
442: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
443: {
444: PC_ASM *osm = (PC_ASM*)pc->data;
445: PetscErrorCode ierr;
446: PetscInt i;
447: KSPConvergedReason reason;
450: for (i=0; i<osm->n_local_true; i++) {
451: KSPSetUp(osm->ksp[i]);
452: KSPGetConvergedReason(osm->ksp[i],&reason);
453: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
454: pc->failedreason = PC_SUBPC_ERROR;
455: }
456: }
457: return(0);
458: }
460: static PetscErrorCode PCApply_ASM(PC pc,Vec x,Vec y)
461: {
462: PC_ASM *osm = (PC_ASM*)pc->data;
464: PetscInt i,n_local_true = osm->n_local_true;
465: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
468: /*
469: Support for limiting the restriction or interpolation to only local
470: subdomain values (leaving the other values 0).
471: */
472: if (!(osm->type & PC_ASM_RESTRICT)) {
473: forward = SCATTER_FORWARD_LOCAL;
474: /* have to zero the work RHS since scatter may leave some slots empty */
475: VecSet(osm->lx, 0.0);
476: }
477: if (!(osm->type & PC_ASM_INTERPOLATE)) {
478: reverse = SCATTER_REVERSE_LOCAL;
479: }
480:
481: if(osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE){
482: /* zero the global and the local solutions */
483: VecZeroEntries(y);
484: VecSet(osm->ly, 0.0);
485:
486: /* Copy the global RHS to local RHS including the ghost nodes */
487: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
488: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
489:
490: /* Restrict local RHS to the overlapping 0-block RHS */
491: VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
492: VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
493:
494: /* do the local solves */
495: for (i = 0; i < n_local_true; ++i) {
496:
497: /* solve the overlapping i-block */
498: KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]);
500: if (osm->lprolongation) { /* interpolate the non-overalapping i-block solution to the local solution (only for restrictive additive) */
501: VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
502: VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
503: }
504: else{ /* interpolate the overalapping i-block solution to the local solution */
505: VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
506: VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
507: }
508:
509: if (i < n_local_true-1) {
510: /* Restrict local RHS to the overlapping (i+1)-block RHS */
511: VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
512: VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
513:
514: if ( osm->loctype == PC_COMPOSITE_MULTIPLICATIVE){
515: /* udpdate the overlapping (i+1)-block RHS using the current local solution */
516: MatMult(osm->lmats[i+1], osm->ly, osm->y[i+1]);
517: VecAXPBY(osm->x[i+1],-1.,1., osm->y[i+1]);
518: }
519: }
520: }
521: /* Add the local solution to the global solution including the ghost nodes */
522: VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
523: VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
524: }else{
525: SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
526: }
527: return(0);
528: }
530: static PetscErrorCode PCApplyTranspose_ASM(PC pc,Vec x,Vec y)
531: {
532: PC_ASM *osm = (PC_ASM*)pc->data;
534: PetscInt i,n_local_true = osm->n_local_true;
535: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
538: /*
539: Support for limiting the restriction or interpolation to only local
540: subdomain values (leaving the other values 0).
542: Note: these are reversed from the PCApply_ASM() because we are applying the
543: transpose of the three terms
544: */
545:
546: if (!(osm->type & PC_ASM_INTERPOLATE)) {
547: forward = SCATTER_FORWARD_LOCAL;
548: /* have to zero the work RHS since scatter may leave some slots empty */
549: VecSet(osm->lx, 0.0);
550: }
551: if (!(osm->type & PC_ASM_RESTRICT)) reverse = SCATTER_REVERSE_LOCAL;
553: /* zero the global and the local solutions */
554: VecZeroEntries(y);
555: VecSet(osm->ly, 0.0);
556:
557: /* Copy the global RHS to local RHS including the ghost nodes */
558: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
559: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
560:
561: /* Restrict local RHS to the overlapping 0-block RHS */
562: VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
563: VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
564:
565: /* do the local solves */
566: for (i = 0; i < n_local_true; ++i) {
567:
568: /* solve the overlapping i-block */
569: KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i]);
571: if (osm->lprolongation) { /* interpolate the non-overalapping i-block solution to the local solution */
572: VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
573: VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
574: }
575: else{ /* interpolate the overalapping i-block solution to the local solution */
576: VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
577: VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
578: }
579:
580: if (i < n_local_true-1) {
581: /* Restrict local RHS to the overlapping (i+1)-block RHS */
582: VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
583: VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
584: }
585: }
586: /* Add the local solution to the global solution including the ghost nodes */
587: VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
588: VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
589:
590: return(0);
591:
592: }
594: static PetscErrorCode PCReset_ASM(PC pc)
595: {
596: PC_ASM *osm = (PC_ASM*)pc->data;
598: PetscInt i;
601: if (osm->ksp) {
602: for (i=0; i<osm->n_local_true; i++) {
603: KSPReset(osm->ksp[i]);
604: }
605: }
606: if (osm->pmat) {
607: if (osm->n_local_true > 0) {
608: MatDestroySubMatrices(osm->n_local_true,&osm->pmat);
609: }
610: }
611: if (osm->lrestriction) {
612: VecScatterDestroy(&osm->restriction);
613: for (i=0; i<osm->n_local_true; i++) {
614: VecScatterDestroy(&osm->lrestriction[i]);
615: if (osm->lprolongation) {VecScatterDestroy(&osm->lprolongation[i]);}
616: VecDestroy(&osm->x[i]);
617: VecDestroy(&osm->y[i]);
618: }
619: PetscFree(osm->lrestriction);
620: if (osm->lprolongation) {PetscFree(osm->lprolongation);}
621: PetscFree(osm->x);
622: PetscFree(osm->y);
623:
624: }
625: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
626: ISDestroy(&osm->lis);
627: VecDestroy(&osm->lx);
628: VecDestroy(&osm->ly);
629: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
630: MatDestroyMatrices(osm->n_local_true, &osm->lmats);
631: }
633: PetscFree(osm->sub_mat_type);
635: osm->is = 0;
636: osm->is_local = 0;
637: return(0);
638: }
640: static PetscErrorCode PCDestroy_ASM(PC pc)
641: {
642: PC_ASM *osm = (PC_ASM*)pc->data;
644: PetscInt i;
647: PCReset_ASM(pc);
648: if (osm->ksp) {
649: for (i=0; i<osm->n_local_true; i++) {
650: KSPDestroy(&osm->ksp[i]);
651: }
652: PetscFree(osm->ksp);
653: }
654: PetscFree(pc->data);
655: return(0);
656: }
658: static PetscErrorCode PCSetFromOptions_ASM(PetscOptionItems *PetscOptionsObject,PC pc)
659: {
660: PC_ASM *osm = (PC_ASM*)pc->data;
662: PetscInt blocks,ovl;
663: PetscBool symset,flg;
664: PCASMType asmtype;
665: PCCompositeType loctype;
666: char sub_mat_type[256];
669: /* set the type to symmetric if matrix is symmetric */
670: if (!osm->type_set && pc->pmat) {
671: MatIsSymmetricKnown(pc->pmat,&symset,&flg);
672: if (symset && flg) osm->type = PC_ASM_BASIC;
673: }
674: PetscOptionsHead(PetscOptionsObject,"Additive Schwarz options");
675: PetscOptionsBool("-pc_asm_dm_subdomains","Use DMCreateDomainDecomposition() to define subdomains","PCASMSetDMSubdomains",osm->dm_subdomains,&osm->dm_subdomains,&flg);
676: PetscOptionsInt("-pc_asm_blocks","Number of subdomains","PCASMSetTotalSubdomains",osm->n,&blocks,&flg);
677: if (flg) {
678: PCASMSetTotalSubdomains(pc,blocks,NULL,NULL);
679: osm->dm_subdomains = PETSC_FALSE;
680: }
681: PetscOptionsInt("-pc_asm_overlap","Number of grid points overlap","PCASMSetOverlap",osm->overlap,&ovl,&flg);
682: if (flg) {
683: PCASMSetOverlap(pc,ovl);
684: osm->dm_subdomains = PETSC_FALSE;
685: }
686: flg = PETSC_FALSE;
687: PetscOptionsEnum("-pc_asm_type","Type of restriction/extension","PCASMSetType",PCASMTypes,(PetscEnum)osm->type,(PetscEnum*)&asmtype,&flg);
688: if (flg) {PCASMSetType(pc,asmtype); }
689: flg = PETSC_FALSE;
690: PetscOptionsEnum("-pc_asm_local_type","Type of local solver composition","PCASMSetLocalType",PCCompositeTypes,(PetscEnum)osm->loctype,(PetscEnum*)&loctype,&flg);
691: if (flg) {PCASMSetLocalType(pc,loctype); }
692: PetscOptionsFList("-pc_asm_sub_mat_type","Subsolve Matrix Type","PCASMSetSubMatType",MatList,NULL,sub_mat_type,256,&flg);
693: if(flg){
694: PCASMSetSubMatType(pc,sub_mat_type);
695: }
696: PetscOptionsTail();
697: return(0);
698: }
700: /*------------------------------------------------------------------------------------*/
702: static PetscErrorCode PCASMSetLocalSubdomains_ASM(PC pc,PetscInt n,IS is[],IS is_local[])
703: {
704: PC_ASM *osm = (PC_ASM*)pc->data;
706: PetscInt i;
709: if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Each process must have 1 or more blocks, n = %D",n);
710: if (pc->setupcalled && (n != osm->n_local_true || is)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetLocalSubdomains() should be called before calling PCSetUp().");
712: if (!pc->setupcalled) {
713: if (is) {
714: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is[i]);}
715: }
716: if (is_local) {
717: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is_local[i]);}
718: }
719: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
721: osm->n_local_true = n;
722: osm->is = 0;
723: osm->is_local = 0;
724: if (is) {
725: PetscMalloc1(n,&osm->is);
726: for (i=0; i<n; i++) osm->is[i] = is[i];
727: /* Flag indicating that the user has set overlapping subdomains so PCASM should not increase their size. */
728: osm->overlap = -1;
729: }
730: if (is_local) {
731: PetscMalloc1(n,&osm->is_local);
732: for (i=0; i<n; i++) osm->is_local[i] = is_local[i];
733: if (!is) {
734: PetscMalloc1(osm->n_local_true,&osm->is);
735: for (i=0; i<osm->n_local_true; i++) {
736: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
737: ISDuplicate(osm->is_local[i],&osm->is[i]);
738: ISCopy(osm->is_local[i],osm->is[i]);
739: } else {
740: PetscObjectReference((PetscObject)osm->is_local[i]);
741: osm->is[i] = osm->is_local[i];
742: }
743: }
744: }
745: }
746: }
747: return(0);
748: }
750: static PetscErrorCode PCASMSetTotalSubdomains_ASM(PC pc,PetscInt N,IS *is,IS *is_local)
751: {
752: PC_ASM *osm = (PC_ASM*)pc->data;
754: PetscMPIInt rank,size;
755: PetscInt n;
758: if (N < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of total blocks must be > 0, N = %D",N);
759: 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.");
761: /*
762: Split the subdomains equally among all processors
763: */
764: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
765: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
766: n = N/size + ((N % size) > rank);
767: 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);
768: if (pc->setupcalled && n != osm->n_local_true) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCASMSetTotalSubdomains() should be called before PCSetUp().");
769: if (!pc->setupcalled) {
770: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
772: osm->n_local_true = n;
773: osm->is = 0;
774: osm->is_local = 0;
775: }
776: return(0);
777: }
779: static PetscErrorCode PCASMSetOverlap_ASM(PC pc,PetscInt ovl)
780: {
781: PC_ASM *osm = (PC_ASM*)pc->data;
784: if (ovl < 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested");
785: if (pc->setupcalled && ovl != osm->overlap) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetOverlap() should be called before PCSetUp().");
786: if (!pc->setupcalled) osm->overlap = ovl;
787: return(0);
788: }
790: static PetscErrorCode PCASMSetType_ASM(PC pc,PCASMType type)
791: {
792: PC_ASM *osm = (PC_ASM*)pc->data;
795: osm->type = type;
796: osm->type_set = PETSC_TRUE;
797: return(0);
798: }
800: static PetscErrorCode PCASMGetType_ASM(PC pc,PCASMType *type)
801: {
802: PC_ASM *osm = (PC_ASM*)pc->data;
805: *type = osm->type;
806: return(0);
807: }
809: static PetscErrorCode PCASMSetLocalType_ASM(PC pc, PCCompositeType type)
810: {
811: PC_ASM *osm = (PC_ASM *) pc->data;
814: 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");
815: osm->loctype = type;
816: return(0);
817: }
819: static PetscErrorCode PCASMGetLocalType_ASM(PC pc, PCCompositeType *type)
820: {
821: PC_ASM *osm = (PC_ASM *) pc->data;
824: *type = osm->loctype;
825: return(0);
826: }
828: static PetscErrorCode PCASMSetSortIndices_ASM(PC pc,PetscBool doSort)
829: {
830: PC_ASM *osm = (PC_ASM*)pc->data;
833: osm->sort_indices = doSort;
834: return(0);
835: }
837: static PetscErrorCode PCASMGetSubKSP_ASM(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
838: {
839: PC_ASM *osm = (PC_ASM*)pc->data;
843: 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");
845: if (n_local) *n_local = osm->n_local_true;
846: if (first_local) {
847: MPI_Scan(&osm->n_local_true,first_local,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
848: *first_local -= osm->n_local_true;
849: }
850: if (ksp) {
851: /* Assume that local solves are now different; not necessarily
852: true though! This flag is used only for PCView_ASM() */
853: *ksp = osm->ksp;
854: osm->same_local_solves = PETSC_FALSE;
855: }
856: return(0);
857: }
859: static PetscErrorCode PCASMGetSubMatType_ASM(PC pc,MatType *sub_mat_type)
860: {
861: PC_ASM *osm = (PC_ASM*)pc->data;
866: *sub_mat_type = osm->sub_mat_type;
867: return(0);
868: }
870: static PetscErrorCode PCASMSetSubMatType_ASM(PC pc,MatType sub_mat_type)
871: {
872: PetscErrorCode ierr;
873: PC_ASM *osm = (PC_ASM*)pc->data;
877: PetscFree(osm->sub_mat_type);
878: PetscStrallocpy(sub_mat_type,(char**)&osm->sub_mat_type);
879: return(0);
880: }
882: /*@C
883: PCASMSetLocalSubdomains - Sets the local subdomains (for this processor only) for the additive Schwarz preconditioner.
885: Collective on PC
887: Input Parameters:
888: + pc - the preconditioner context
889: . n - the number of subdomains for this processor (default value = 1)
890: . is - the index set that defines the subdomains for this processor
891: (or NULL for PETSc to determine subdomains)
892: - is_local - the index sets that define the local part of the subdomains for this processor, not used unless PCASMType is PC_ASM_RESTRICT
893: (or NULL to not provide these)
895: Notes:
896: The IS numbering is in the parallel, global numbering of the vector for both is and is_local
898: By default the ASM preconditioner uses 1 block per processor.
900: Use PCASMSetTotalSubdomains() to set the subdomains for all processors.
902: If is_local is provided and PCASMType is PC_ASM_RESTRICT then the solution only over the is_local region is interpolated
903: back to form the global solution (this is the standard restricted additive Schwarz method)
905: If the is_local is provided and PCASMType is PC_ASM_INTERPOLATE or PC_ASM_NONE then an error is generated since there is
906: no code to handle that case.
908: Level: advanced
910: .keywords: PC, ASM, set, local, subdomains, additive Schwarz
912: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
913: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), PCASMType, PCASMSetType()
914: @*/
915: PetscErrorCode PCASMSetLocalSubdomains(PC pc,PetscInt n,IS is[],IS is_local[])
916: {
921: PetscTryMethod(pc,"PCASMSetLocalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,is,is_local));
922: return(0);
923: }
925: /*@C
926: PCASMSetTotalSubdomains - Sets the subdomains for all processors for the
927: additive Schwarz preconditioner. Either all or no processors in the
928: PC communicator must call this routine, with the same index sets.
930: Collective on PC
932: Input Parameters:
933: + pc - the preconditioner context
934: . N - the number of subdomains for all processors
935: . is - the index sets that define the subdomains for all processors
936: (or NULL to ask PETSc to determine the subdomains)
937: - is_local - the index sets that define the local part of the subdomains for this processor
938: (or NULL to not provide this information)
940: Options Database Key:
941: To set the total number of subdomain blocks rather than specify the
942: index sets, use the option
943: . -pc_asm_blocks <blks> - Sets total blocks
945: Notes:
946: Currently you cannot use this to set the actual subdomains with the argument is or is_local.
948: By default the ASM preconditioner uses 1 block per processor.
950: These index sets cannot be destroyed until after completion of the
951: linear solves for which the ASM preconditioner is being used.
953: Use PCASMSetLocalSubdomains() to set local subdomains.
955: The IS numbering is in the parallel, global numbering of the vector for both is and is_local
957: Level: advanced
959: .keywords: PC, ASM, set, total, global, subdomains, additive Schwarz
961: .seealso: PCASMSetLocalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
962: PCASMCreateSubdomains2D()
963: @*/
964: PetscErrorCode PCASMSetTotalSubdomains(PC pc,PetscInt N,IS is[],IS is_local[])
965: {
970: PetscTryMethod(pc,"PCASMSetTotalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,N,is,is_local));
971: return(0);
972: }
974: /*@
975: PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
976: additive Schwarz preconditioner. Either all or no processors in the
977: PC communicator must call this routine.
979: Logically Collective on PC
981: Input Parameters:
982: + pc - the preconditioner context
983: - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1)
985: Options Database Key:
986: . -pc_asm_overlap <ovl> - Sets overlap
988: Notes:
989: By default the ASM preconditioner uses 1 block per processor. To use
990: multiple blocks per perocessor, see PCASMSetTotalSubdomains() and
991: PCASMSetLocalSubdomains() (and the option -pc_asm_blocks <blks>).
993: The overlap defaults to 1, so if one desires that no additional
994: overlap be computed beyond what may have been set with a call to
995: PCASMSetTotalSubdomains() or PCASMSetLocalSubdomains(), then ovl
996: must be set to be 0. In particular, if one does not explicitly set
997: the subdomains an application code, then all overlap would be computed
998: internally by PETSc, and using an overlap of 0 would result in an ASM
999: variant that is equivalent to the block Jacobi preconditioner.
1001: The default algorithm used by PETSc to increase overlap is fast, but not scalable,
1002: use the option -mat_increase_overlap_scalable when the problem and number of processes is large.
1004: Note that one can define initial index sets with any overlap via
1005: PCASMSetLocalSubdomains(); the routine
1006: PCASMSetOverlap() merely allows PETSc to extend that overlap further
1007: if desired.
1009: Level: intermediate
1011: .keywords: PC, ASM, set, overlap
1013: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1014: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), MatIncreaseOverlap()
1015: @*/
1016: PetscErrorCode PCASMSetOverlap(PC pc,PetscInt ovl)
1017: {
1023: PetscTryMethod(pc,"PCASMSetOverlap_C",(PC,PetscInt),(pc,ovl));
1024: return(0);
1025: }
1027: /*@
1028: PCASMSetType - Sets the type of restriction and interpolation used
1029: for local problems in the additive Schwarz method.
1031: Logically Collective on PC
1033: Input Parameters:
1034: + pc - the preconditioner context
1035: - type - variant of ASM, one of
1036: .vb
1037: PC_ASM_BASIC - full interpolation and restriction
1038: PC_ASM_RESTRICT - full restriction, local processor interpolation
1039: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1040: PC_ASM_NONE - local processor restriction and interpolation
1041: .ve
1043: Options Database Key:
1044: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
1046: Notes: if the is_local arguments are passed to PCASMSetLocalSubdomains() then they are used when PC_ASM_RESTRICT has been selected
1047: to limit the local processor interpolation
1049: Level: intermediate
1051: .keywords: PC, ASM, set, type
1053: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1054: PCASMCreateSubdomains2D(), PCASMType, PCASMSetLocalType(), PCASMGetLocalType()
1055: @*/
1056: PetscErrorCode PCASMSetType(PC pc,PCASMType type)
1057: {
1063: PetscTryMethod(pc,"PCASMSetType_C",(PC,PCASMType),(pc,type));
1064: return(0);
1065: }
1067: /*@
1068: PCASMGetType - Gets the type of restriction and interpolation used
1069: for local problems in the additive Schwarz method.
1071: Logically Collective on PC
1073: Input Parameter:
1074: . pc - the preconditioner context
1076: Output Parameter:
1077: . type - variant of ASM, one of
1079: .vb
1080: PC_ASM_BASIC - full interpolation and restriction
1081: PC_ASM_RESTRICT - full restriction, local processor interpolation
1082: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1083: PC_ASM_NONE - local processor restriction and interpolation
1084: .ve
1086: Options Database Key:
1087: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
1089: Level: intermediate
1091: .keywords: PC, ASM, set, type
1093: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1094: PCASMCreateSubdomains2D(), PCASMType, PCASMSetType(), PCASMSetLocalType(), PCASMGetLocalType()
1095: @*/
1096: PetscErrorCode PCASMGetType(PC pc,PCASMType *type)
1097: {
1102: PetscUseMethod(pc,"PCASMGetType_C",(PC,PCASMType*),(pc,type));
1103: return(0);
1104: }
1106: /*@
1107: PCASMSetLocalType - Sets the type of composition used for local problems in the additive Schwarz method.
1109: Logically Collective on PC
1111: Input Parameters:
1112: + pc - the preconditioner context
1113: - type - type of composition, one of
1114: .vb
1115: PC_COMPOSITE_ADDITIVE - local additive combination
1116: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1117: .ve
1119: Options Database Key:
1120: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1122: Level: intermediate
1124: .seealso: PCASMSetType(), PCASMGetType(), PCASMGetLocalType(), PCASM, PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1125: @*/
1126: PetscErrorCode PCASMSetLocalType(PC pc, PCCompositeType type)
1127: {
1133: PetscTryMethod(pc, "PCASMSetLocalType_C", (PC, PCCompositeType), (pc, type));
1134: return(0);
1135: }
1137: /*@
1138: PCASMGetLocalType - Gets the type of composition used for local problems in the additive Schwarz method.
1140: Logically Collective on PC
1142: Input Parameter:
1143: . pc - the preconditioner context
1145: Output Parameter:
1146: . type - type of composition, one of
1147: .vb
1148: PC_COMPOSITE_ADDITIVE - local additive combination
1149: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1150: .ve
1152: Options Database Key:
1153: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1155: Level: intermediate
1157: .seealso: PCASMSetType(), PCASMGetType(), PCASMSetLocalType(), PCASMCreate(), PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1158: @*/
1159: PetscErrorCode PCASMGetLocalType(PC pc, PCCompositeType *type)
1160: {
1166: PetscUseMethod(pc, "PCASMGetLocalType_C", (PC, PCCompositeType *), (pc, type));
1167: return(0);
1168: }
1170: /*@
1171: PCASMSetSortIndices - Determines whether subdomain indices are sorted.
1173: Logically Collective on PC
1175: Input Parameters:
1176: + pc - the preconditioner context
1177: - doSort - sort the subdomain indices
1179: Level: intermediate
1181: .keywords: PC, ASM, set, type
1183: .seealso: PCASMSetLocalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1184: PCASMCreateSubdomains2D()
1185: @*/
1186: PetscErrorCode PCASMSetSortIndices(PC pc,PetscBool doSort)
1187: {
1193: PetscTryMethod(pc,"PCASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));
1194: return(0);
1195: }
1197: /*@C
1198: PCASMGetSubKSP - Gets the local KSP contexts for all blocks on
1199: this processor.
1201: Collective on PC iff first_local is requested
1203: Input Parameter:
1204: . pc - the preconditioner context
1206: Output Parameters:
1207: + n_local - the number of blocks on this processor or NULL
1208: . first_local - the global number of the first block on this processor or NULL,
1209: all processors must request or all must pass NULL
1210: - ksp - the array of KSP contexts
1212: Note:
1213: After PCASMGetSubKSP() the array of KSPes is not to be freed.
1215: You must call KSPSetUp() before calling PCASMGetSubKSP().
1217: Fortran note:
1218: 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.
1220: Level: advanced
1222: .keywords: PC, ASM, additive Schwarz, get, sub, KSP, context
1224: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap(),
1225: PCASMCreateSubdomains2D(),
1226: @*/
1227: PetscErrorCode PCASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
1228: {
1233: PetscUseMethod(pc,"PCASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
1234: return(0);
1235: }
1237: /* -------------------------------------------------------------------------------------*/
1238: /*MC
1239: PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1240: its own KSP object.
1242: Options Database Keys:
1243: + -pc_asm_blocks <blks> - Sets total blocks
1244: . -pc_asm_overlap <ovl> - Sets overlap
1245: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type, default is restrict
1246: - -pc_asm_local_type [additive, multiplicative] - Sets ASM type, default is additive
1248: IMPORTANT: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you
1249: will get a different convergence rate due to the default option of -pc_asm_type restrict. Use
1250: -pc_asm_type basic to use the standard ASM.
1252: Notes: Each processor can have one or more blocks, but a block cannot be shared by more
1253: than one processor. Use PCGASM for subdomains shared by multiple processes. Defaults to one block per processor.
1255: To set options on the solvers for each block append -sub_ to all the KSP, and PC
1256: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
1258: To set the options on the solvers separate for each block call PCASMGetSubKSP()
1259: and set the options directly on the resulting KSP object (you can access its PC
1260: with KSPGetPC())
1262: Level: beginner
1264: Concepts: additive Schwarz method
1266: References:
1267: + 1. - M Dryja, OB Widlund, An additive variant of the Schwarz alternating method for the case of many subregions
1268: Courant Institute, New York University Technical report
1269: - 1. - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations,
1270: Cambridge University Press.
1272: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
1273: PCBJACOBI, PCASMGetSubKSP(), PCASMSetLocalSubdomains(), PCASMType, PCASMGetType(), PCASMSetLocalType(), PCASMGetLocalType()
1274: PCASMSetTotalSubdomains(), PCSetModifySubmatrices(), PCASMSetOverlap(), PCASMSetType(), PCCompositeType
1276: M*/
1278: PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1279: {
1281: PC_ASM *osm;
1284: PetscNewLog(pc,&osm);
1286: osm->n = PETSC_DECIDE;
1287: osm->n_local = 0;
1288: osm->n_local_true = PETSC_DECIDE;
1289: osm->overlap = 1;
1290: osm->ksp = 0;
1291: osm->restriction = 0;
1292: osm->lprolongation = 0;
1293: osm->lrestriction = 0;
1294: osm->x = 0;
1295: osm->y = 0;
1296: osm->is = 0;
1297: osm->is_local = 0;
1298: osm->mat = 0;
1299: osm->pmat = 0;
1300: osm->type = PC_ASM_RESTRICT;
1301: osm->loctype = PC_COMPOSITE_ADDITIVE;
1302: osm->same_local_solves = PETSC_TRUE;
1303: osm->sort_indices = PETSC_TRUE;
1304: osm->dm_subdomains = PETSC_FALSE;
1305: osm->sub_mat_type = NULL;
1307: pc->data = (void*)osm;
1308: pc->ops->apply = PCApply_ASM;
1309: pc->ops->applytranspose = PCApplyTranspose_ASM;
1310: pc->ops->setup = PCSetUp_ASM;
1311: pc->ops->reset = PCReset_ASM;
1312: pc->ops->destroy = PCDestroy_ASM;
1313: pc->ops->setfromoptions = PCSetFromOptions_ASM;
1314: pc->ops->setuponblocks = PCSetUpOnBlocks_ASM;
1315: pc->ops->view = PCView_ASM;
1316: pc->ops->applyrichardson = 0;
1318: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",PCASMSetLocalSubdomains_ASM);
1319: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",PCASMSetTotalSubdomains_ASM);
1320: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetOverlap_C",PCASMSetOverlap_ASM);
1321: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetType_C",PCASMSetType_ASM);
1322: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetType_C",PCASMGetType_ASM);
1323: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalType_C",PCASMSetLocalType_ASM);
1324: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetLocalType_C",PCASMGetLocalType_ASM);
1325: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSortIndices_C",PCASMSetSortIndices_ASM);
1326: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubKSP_C",PCASMGetSubKSP_ASM);
1327: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubMatType_C",PCASMGetSubMatType_ASM);
1328: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSubMatType_C",PCASMSetSubMatType_ASM);
1329: return(0);
1330: }
1332: /*@C
1333: PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1334: preconditioner for a any problem on a general grid.
1336: Collective
1338: Input Parameters:
1339: + A - The global matrix operator
1340: - n - the number of local blocks
1342: Output Parameters:
1343: . outis - the array of index sets defining the subdomains
1345: Level: advanced
1347: Note: this generates nonoverlapping subdomains; the PCASM will generate the overlap
1348: from these if you use PCASMSetLocalSubdomains()
1350: In the Fortran version you must provide the array outis[] already allocated of length n.
1352: .keywords: PC, ASM, additive Schwarz, create, subdomains, unstructured grid
1354: .seealso: PCASMSetLocalSubdomains(), PCASMDestroySubdomains()
1355: @*/
1356: PetscErrorCode PCASMCreateSubdomains(Mat A, PetscInt n, IS* outis[])
1357: {
1358: MatPartitioning mpart;
1359: const char *prefix;
1360: PetscInt i,j,rstart,rend,bs;
1361: PetscBool hasop, isbaij = PETSC_FALSE,foundpart = PETSC_FALSE;
1362: Mat Ad = NULL, adj;
1363: IS ispart,isnumb,*is;
1364: PetscErrorCode ierr;
1369: if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of local blocks must be > 0, n = %D",n);
1371: /* Get prefix, row distribution, and block size */
1372: MatGetOptionsPrefix(A,&prefix);
1373: MatGetOwnershipRange(A,&rstart,&rend);
1374: MatGetBlockSize(A,&bs);
1375: 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);
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: PetscObjectTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij);
1384: if (!isbaij) {PetscObjectTypeCompare((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=0,*jja=0;
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: .keywords: PC, ASM, additive Schwarz, create, subdomains, unstructured grid
1526: .seealso: PCASMCreateSubdomains(), PCASMSetLocalSubdomains()
1527: @*/
1528: PetscErrorCode PCASMDestroySubdomains(PetscInt n, IS is[], IS is_local[])
1529: {
1530: PetscInt i;
1534: if (n <= 0) return(0);
1535: if (is) {
1537: for (i=0; i<n; i++) { ISDestroy(&is[i]); }
1538: PetscFree(is);
1539: }
1540: if (is_local) {
1542: for (i=0; i<n; i++) { ISDestroy(&is_local[i]); }
1543: PetscFree(is_local);
1544: }
1545: return(0);
1546: }
1548: /*@
1549: PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1550: preconditioner for a two-dimensional problem on a regular grid.
1552: Not Collective
1554: Input Parameters:
1555: + m, n - the number of mesh points in the x and y directions
1556: . M, N - the number of subdomains in the x and y directions
1557: . dof - degrees of freedom per node
1558: - overlap - overlap in mesh lines
1560: Output Parameters:
1561: + Nsub - the number of subdomains created
1562: . is - array of index sets defining overlapping (if overlap > 0) subdomains
1563: - is_local - array of index sets defining non-overlapping subdomains
1565: Note:
1566: Presently PCAMSCreateSubdomains2d() is valid only for sequential
1567: preconditioners. More general related routines are
1568: PCASMSetTotalSubdomains() and PCASMSetLocalSubdomains().
1570: Level: advanced
1572: .keywords: PC, ASM, additive Schwarz, create, subdomains, 2D, regular grid
1574: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1575: PCASMSetOverlap()
1576: @*/
1577: PetscErrorCode PCASMCreateSubdomains2D(PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt dof,PetscInt overlap,PetscInt *Nsub,IS **is,IS **is_local)
1578: {
1579: PetscInt i,j,height,width,ystart,xstart,yleft,yright,xleft,xright,loc_outer;
1581: PetscInt nidx,*idx,loc,ii,jj,count;
1584: if (dof != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP," ");
1586: *Nsub = N*M;
1587: PetscMalloc1(*Nsub,is);
1588: PetscMalloc1(*Nsub,is_local);
1589: ystart = 0;
1590: loc_outer = 0;
1591: for (i=0; i<N; i++) {
1592: height = n/N + ((n % N) > i); /* height of subdomain */
1593: if (height < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many N subdomains for mesh dimension n");
1594: yleft = ystart - overlap; if (yleft < 0) yleft = 0;
1595: yright = ystart + height + overlap; if (yright > n) yright = n;
1596: xstart = 0;
1597: for (j=0; j<M; j++) {
1598: width = m/M + ((m % M) > j); /* width of subdomain */
1599: if (width < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many M subdomains for mesh dimension m");
1600: xleft = xstart - overlap; if (xleft < 0) xleft = 0;
1601: xright = xstart + width + overlap; if (xright > m) xright = m;
1602: nidx = (xright - xleft)*(yright - yleft);
1603: PetscMalloc1(nidx,&idx);
1604: loc = 0;
1605: for (ii=yleft; ii<yright; ii++) {
1606: count = m*ii + xleft;
1607: for (jj=xleft; jj<xright; jj++) idx[loc++] = count++;
1608: }
1609: ISCreateGeneral(PETSC_COMM_SELF,nidx,idx,PETSC_COPY_VALUES,(*is)+loc_outer);
1610: if (overlap == 0) {
1611: PetscObjectReference((PetscObject)(*is)[loc_outer]);
1613: (*is_local)[loc_outer] = (*is)[loc_outer];
1614: } else {
1615: for (loc=0,ii=ystart; ii<ystart+height; ii++) {
1616: for (jj=xstart; jj<xstart+width; jj++) {
1617: idx[loc++] = m*ii + jj;
1618: }
1619: }
1620: ISCreateGeneral(PETSC_COMM_SELF,loc,idx,PETSC_COPY_VALUES,*is_local+loc_outer);
1621: }
1622: PetscFree(idx);
1623: xstart += width;
1624: loc_outer++;
1625: }
1626: ystart += height;
1627: }
1628: for (i=0; i<*Nsub; i++) { ISSort((*is)[i]); }
1629: return(0);
1630: }
1632: /*@C
1633: PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
1634: only) for the additive Schwarz preconditioner.
1636: Not Collective
1638: Input Parameter:
1639: . pc - the preconditioner context
1641: Output Parameters:
1642: + n - the number of subdomains for this processor (default value = 1)
1643: . is - the index sets that define the subdomains for this processor
1644: - is_local - the index sets that define the local part of the subdomains for this processor (can be NULL)
1647: Notes:
1648: The IS numbering is in the parallel, global numbering of the vector.
1650: Level: advanced
1652: .keywords: PC, ASM, set, local, subdomains, additive Schwarz
1654: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1655: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubmatrices()
1656: @*/
1657: PetscErrorCode PCASMGetLocalSubdomains(PC pc,PetscInt *n,IS *is[],IS *is_local[])
1658: {
1659: PC_ASM *osm = (PC_ASM*)pc->data;
1661: PetscBool match;
1667: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1668: if (!match) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PC is not a PCASM");
1669: if (n) *n = osm->n_local_true;
1670: if (is) *is = osm->is;
1671: if (is_local) *is_local = osm->is_local;
1672: return(0);
1673: }
1675: /*@C
1676: PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1677: only) for the additive Schwarz preconditioner.
1679: Not Collective
1681: Input Parameter:
1682: . pc - the preconditioner context
1684: Output Parameters:
1685: + n - the number of matrices for this processor (default value = 1)
1686: - mat - the matrices
1688: Level: advanced
1690: Notes: Call after PCSetUp() (or KSPSetUp()) but before PCApply() (or KSPApply()) and before PCSetUpOnBlocks())
1692: Usually one would use PCSetModifySubmatrices() to change the submatrices in building the preconditioner.
1694: .keywords: PC, ASM, set, local, subdomains, additive Schwarz, block Jacobi
1696: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1697: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains(), PCSetModifySubmatrices()
1698: @*/
1699: PetscErrorCode PCASMGetLocalSubmatrices(PC pc,PetscInt *n,Mat *mat[])
1700: {
1701: PC_ASM *osm;
1703: PetscBool match;
1709: if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUP() or PCSetUp().");
1710: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1711: if (!match) {
1712: if (n) *n = 0;
1713: if (mat) *mat = NULL;
1714: } else {
1715: osm = (PC_ASM*)pc->data;
1716: if (n) *n = osm->n_local_true;
1717: if (mat) *mat = osm->pmat;
1718: }
1719: return(0);
1720: }
1722: /*@
1723: PCASMSetDMSubdomains - Indicates whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1725: Logically Collective
1727: Input Parameter:
1728: + pc - the preconditioner
1729: - flg - boolean indicating whether to use subdomains defined by the DM
1731: Options Database Key:
1732: . -pc_asm_dm_subdomains
1734: Level: intermediate
1736: Notes:
1737: PCASMSetTotalSubdomains() and PCASMSetOverlap() take precedence over PCASMSetDMSubdomains(),
1738: so setting either of the first two effectively turns the latter off.
1740: .keywords: PC, ASM, DM, set, subdomains, additive Schwarz
1742: .seealso: PCASMGetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1743: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1744: @*/
1745: PetscErrorCode PCASMSetDMSubdomains(PC pc,PetscBool flg)
1746: {
1747: PC_ASM *osm = (PC_ASM*)pc->data;
1749: PetscBool match;
1754: if (pc->setupcalled) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for a setup PC.");
1755: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1756: if (match) {
1757: osm->dm_subdomains = flg;
1758: }
1759: return(0);
1760: }
1762: /*@
1763: PCASMGetDMSubdomains - Returns flag indicating whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1764: Not Collective
1766: Input Parameter:
1767: . pc - the preconditioner
1769: Output Parameter:
1770: . flg - boolean indicating whether to use subdomains defined by the DM
1772: Level: intermediate
1774: .keywords: PC, ASM, DM, set, subdomains, additive Schwarz
1776: .seealso: PCASMSetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1777: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1778: @*/
1779: PetscErrorCode PCASMGetDMSubdomains(PC pc,PetscBool* flg)
1780: {
1781: PC_ASM *osm = (PC_ASM*)pc->data;
1783: PetscBool match;
1788: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1789: if (match) {
1790: if (flg) *flg = osm->dm_subdomains;
1791: }
1792: return(0);
1793: }
1795: /*@
1796: PCASMGetSubMatType - Gets the matrix type used for ASM subsolves, as a string.
1798: Not Collective
1800: Input Parameter:
1801: . pc - the PC
1803: Output Parameter:
1804: . -pc_asm_sub_mat_type - name of matrix type
1806: Level: advanced
1808: .keywords: PC, PCASM, MatType, set
1810: .seealso: PCASMSetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1811: @*/
1812: PetscErrorCode PCASMGetSubMatType(PC pc,MatType *sub_mat_type){
1815: PetscTryMethod(pc,"PCASMGetSubMatType_C",(PC,MatType*),(pc,sub_mat_type));
1816: return(0);
1817: }
1819: /*@
1820: PCASMSetSubMatType - Set the type of matrix used for ASM subsolves
1822: Collective on Mat
1824: Input Parameters:
1825: + pc - the PC object
1826: - sub_mat_type - matrix type
1828: Options Database Key:
1829: . -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.
1831: Notes:
1832: See "${PETSC_DIR}/include/petscmat.h" for available types
1834: Level: advanced
1836: .keywords: PC, PCASM, MatType, set
1838: .seealso: PCASMGetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1839: @*/
1840: PetscErrorCode PCASMSetSubMatType(PC pc,MatType sub_mat_type)
1841: {
1844: PetscTryMethod(pc,"PCASMSetSubMatType_C",(PC,MatType),(pc,sub_mat_type));
1845: return(0);
1846: }