Actual source code: asm.c
petsc-3.11.4 2019-09-28
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: PetscMalloc1(16*(nidx+1)+512, &s);
124: PetscViewerStringOpen(PETSC_COMM_SELF, s, 16*(nidx+1)+512, &sviewer);
125: PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D with overlap:\n", rank, size, i);
126: for (j=0; j<nidx; j++) {
127: PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
128: }
129: ISRestoreIndices(osm->is[i],&idx);
130: PetscViewerStringSPrintf(sviewer,"\n");
131: PetscViewerDestroy(&sviewer);
132: PetscViewerASCIIPushSynchronized(viewer);
133: PetscViewerASCIISynchronizedPrintf(viewer, s);
134: PetscViewerFlush(viewer);
135: PetscViewerASCIIPopSynchronized(viewer);
136: PetscFree(s);
137: if (osm->is_local) {
138: /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
139: PetscMalloc1(16*(nidx+1)+512, &s);
140: PetscViewerStringOpen(PETSC_COMM_SELF, s, 16*(nidx+1)+512, &sviewer);
141: PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D without overlap:\n", rank, size, i);
142: ISGetLocalSize(osm->is_local[i],&nidx);
143: ISGetIndices(osm->is_local[i],&idx);
144: for (j=0; j<nidx; j++) {
145: PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
146: }
147: ISRestoreIndices(osm->is_local[i],&idx);
148: PetscViewerStringSPrintf(sviewer,"\n");
149: PetscViewerDestroy(&sviewer);
150: PetscViewerASCIIPushSynchronized(viewer);
151: PetscViewerASCIISynchronizedPrintf(viewer, s);
152: PetscViewerFlush(viewer);
153: PetscViewerASCIIPopSynchronized(viewer);
154: PetscFree(s);
155: }
156: } else {
157: /* Participate in collective viewer calls. */
158: PetscViewerASCIIPushSynchronized(viewer);
159: PetscViewerFlush(viewer);
160: PetscViewerASCIIPopSynchronized(viewer);
161: /* Assume either all ranks have is_local or none do. */
162: if (osm->is_local) {
163: PetscViewerASCIIPushSynchronized(viewer);
164: PetscViewerFlush(viewer);
165: PetscViewerASCIIPopSynchronized(viewer);
166: }
167: }
168: }
169: PetscViewerFlush(viewer);
170: PetscViewerDestroy(&viewer);
171: return(0);
172: }
174: static PetscErrorCode PCSetUp_ASM(PC pc)
175: {
176: PC_ASM *osm = (PC_ASM*)pc->data;
178: PetscBool symset,flg;
179: PetscInt i,m,m_local;
180: MatReuse scall = MAT_REUSE_MATRIX;
181: IS isl;
182: KSP ksp;
183: PC subpc;
184: const char *prefix,*pprefix;
185: Vec vec;
186: DM *domain_dm = NULL;
189: if (!pc->setupcalled) {
190: PetscInt m;
192: if (!osm->type_set) {
193: MatIsSymmetricKnown(pc->pmat,&symset,&flg);
194: if (symset && flg) osm->type = PC_ASM_BASIC;
195: }
197: /* Note: if subdomains have been set either via PCASMSetTotalSubdomains() or via PCASMSetLocalSubdomains(), osm->n_local_true will not be PETSC_DECIDE */
198: if (osm->n_local_true == PETSC_DECIDE) {
199: /* no subdomains given */
200: /* try pc->dm first, if allowed */
201: if (osm->dm_subdomains && pc->dm) {
202: PetscInt num_domains, d;
203: char **domain_names;
204: IS *inner_domain_is, *outer_domain_is;
205: DMCreateDomainDecomposition(pc->dm, &num_domains, &domain_names, &inner_domain_is, &outer_domain_is, &domain_dm);
206: osm->overlap = -1; /* We do not want to increase the overlap of the IS.
207: A future improvement of this code might allow one to use
208: DM-defined subdomains and also increase the overlap,
209: but that is not currently supported */
210: if (num_domains) {
211: PCASMSetLocalSubdomains(pc, num_domains, outer_domain_is, inner_domain_is);
212: }
213: for (d = 0; d < num_domains; ++d) {
214: if (domain_names) {PetscFree(domain_names[d]);}
215: if (inner_domain_is) {ISDestroy(&inner_domain_is[d]);}
216: if (outer_domain_is) {ISDestroy(&outer_domain_is[d]);}
217: }
218: PetscFree(domain_names);
219: PetscFree(inner_domain_is);
220: PetscFree(outer_domain_is);
221: }
222: if (osm->n_local_true == PETSC_DECIDE) {
223: /* still no subdomains; use one subdomain per processor */
224: osm->n_local_true = 1;
225: }
226: }
227: { /* determine the global and max number of subdomains */
228: struct {PetscInt max,sum;} inwork,outwork;
229: PetscMPIInt size;
231: inwork.max = osm->n_local_true;
232: inwork.sum = osm->n_local_true;
233: MPIU_Allreduce(&inwork,&outwork,1,MPIU_2INT,MPIU_MAXSUM_OP,PetscObjectComm((PetscObject)pc));
234: osm->n_local = outwork.max;
235: osm->n = outwork.sum;
237: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
238: if (outwork.max == 1 && outwork.sum == size) {
239: /* osm->n_local_true = 1 on all processes, set this option may enable use of optimized MatCreateSubMatrices() implementation */
240: MatSetOption(pc->pmat,MAT_SUBMAT_SINGLEIS,PETSC_TRUE);
241: }
242: }
243: if (!osm->is) { /* create the index sets */
244: PCASMCreateSubdomains(pc->pmat,osm->n_local_true,&osm->is);
245: }
246: if (osm->n_local_true > 1 && !osm->is_local) {
247: PetscMalloc1(osm->n_local_true,&osm->is_local);
248: for (i=0; i<osm->n_local_true; i++) {
249: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
250: ISDuplicate(osm->is[i],&osm->is_local[i]);
251: ISCopy(osm->is[i],osm->is_local[i]);
252: } else {
253: PetscObjectReference((PetscObject)osm->is[i]);
254: osm->is_local[i] = osm->is[i];
255: }
256: }
257: }
258: PCGetOptionsPrefix(pc,&prefix);
259: flg = PETSC_FALSE;
260: PetscOptionsGetBool(NULL,prefix,"-pc_asm_print_subdomains",&flg,NULL);
261: if (flg) { PCASMPrintSubdomains(pc); }
263: if (osm->overlap > 0) {
264: /* Extend the "overlapping" regions by a number of steps */
265: MatIncreaseOverlap(pc->pmat,osm->n_local_true,osm->is,osm->overlap);
266: }
267: if (osm->sort_indices) {
268: for (i=0; i<osm->n_local_true; i++) {
269: ISSort(osm->is[i]);
270: if (osm->is_local) {
271: ISSort(osm->is_local[i]);
272: }
273: }
274: }
276: if (!osm->ksp) {
277: /* Create the local solvers */
278: PetscMalloc1(osm->n_local_true,&osm->ksp);
279: if (domain_dm) {
280: PetscInfo(pc,"Setting up ASM subproblems using the embedded DM\n");
281: }
282: for (i=0; i<osm->n_local_true; i++) {
283: KSPCreate(PETSC_COMM_SELF,&ksp);
284: KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
285: PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
286: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
287: KSPSetType(ksp,KSPPREONLY);
288: KSPGetPC(ksp,&subpc);
289: PCGetOptionsPrefix(pc,&prefix);
290: KSPSetOptionsPrefix(ksp,prefix);
291: KSPAppendOptionsPrefix(ksp,"sub_");
292: if (domain_dm) {
293: KSPSetDM(ksp, domain_dm[i]);
294: KSPSetDMActive(ksp, PETSC_FALSE);
295: DMDestroy(&domain_dm[i]);
296: }
297: osm->ksp[i] = ksp;
298: }
299: if (domain_dm) {
300: PetscFree(domain_dm);
301: }
302: }
304: ISConcatenate(PETSC_COMM_SELF, osm->n_local_true, osm->is, &osm->lis);
305: ISSortRemoveDups(osm->lis);
306: ISGetLocalSize(osm->lis, &m);
307: VecCreateSeq(PETSC_COMM_SELF, m, &osm->lx);
308: VecDuplicate(osm->lx, &osm->ly);
310: scall = MAT_INITIAL_MATRIX;
311: } else {
312: /*
313: Destroy the blocks from the previous iteration
314: */
315: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
316: MatDestroyMatrices(osm->n_local_true,&osm->pmat);
317: scall = MAT_INITIAL_MATRIX;
318: }
319: }
321: /*
322: Extract out the submatrices
323: */
324: MatCreateSubMatrices(pc->pmat,osm->n_local_true,osm->is,osm->is,scall,&osm->pmat);
325: if (scall == MAT_INITIAL_MATRIX) {
326: PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
327: for (i=0; i<osm->n_local_true; i++) {
328: PetscLogObjectParent((PetscObject)pc,(PetscObject)osm->pmat[i]);
329: PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);
330: }
331: }
333: /* Convert the types of the submatrices (if needbe) */
334: if (osm->sub_mat_type) {
335: for (i=0; i<osm->n_local_true; i++) {
336: MatConvert(osm->pmat[i],osm->sub_mat_type,MAT_INPLACE_MATRIX,&(osm->pmat[i]));
337: }
338: }
340: if(!pc->setupcalled){
341: /* Create the local work vectors (from the local matrices) and scatter contexts */
342: MatCreateVecs(pc->pmat,&vec,0);
344: 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()");
345: if (osm->is_local && osm->type == PC_ASM_RESTRICT && osm->loctype == PC_COMPOSITE_ADDITIVE) {
346: PetscMalloc1(osm->n_local_true,&osm->lprolongation);
347: }
348: PetscMalloc1(osm->n_local_true,&osm->lrestriction);
349: PetscMalloc1(osm->n_local_true,&osm->x);
350: PetscMalloc1(osm->n_local_true,&osm->y);
352: ISGetLocalSize(osm->lis,&m);
353: ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
354: VecScatterCreate(vec,osm->lis,osm->lx,isl,&osm->restriction);
355: ISDestroy(&isl);
358: for (i=0; i<osm->n_local_true; ++i) {
359: ISLocalToGlobalMapping ltog;
360: IS isll;
361: const PetscInt *idx_is;
362: PetscInt *idx_lis,nout;
364: ISGetLocalSize(osm->is[i],&m);
365: MatCreateVecs(osm->pmat[i],&osm->x[i],NULL);
366: VecDuplicate(osm->x[i],&osm->y[i]);
368: /* generate a scatter from ly to y[i] picking all the overlapping is[i] entries */
369: ISLocalToGlobalMappingCreateIS(osm->lis,<og);
370: ISGetLocalSize(osm->is[i],&m);
371: ISGetIndices(osm->is[i], &idx_is);
372: PetscMalloc1(m,&idx_lis);
373: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m,idx_is,&nout,idx_lis);
374: if (nout != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is not a subset of lis");
375: ISRestoreIndices(osm->is[i], &idx_is);
376: ISCreateGeneral(PETSC_COMM_SELF,m,idx_lis,PETSC_OWN_POINTER,&isll);
377: ISLocalToGlobalMappingDestroy(<og);
378: ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
379: VecScatterCreate(osm->ly,isll,osm->y[i],isl,&osm->lrestriction[i]);
380: ISDestroy(&isll);
381: ISDestroy(&isl);
382: if (osm->lprolongation) { /* generate a scatter from y[i] to ly picking only the the non-overalapping is_local[i] entries */
383: ISLocalToGlobalMapping ltog;
384: IS isll,isll_local;
385: const PetscInt *idx_local;
386: PetscInt *idx1, *idx2, nout;
388: ISGetLocalSize(osm->is_local[i],&m_local);
389: ISGetIndices(osm->is_local[i], &idx_local);
391: ISLocalToGlobalMappingCreateIS(osm->is[i],<og);
392: PetscMalloc1(m_local,&idx1);
393: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx1);
394: ISLocalToGlobalMappingDestroy(<og);
395: if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of is");
396: ISCreateGeneral(PETSC_COMM_SELF,m_local,idx1,PETSC_OWN_POINTER,&isll);
398: ISLocalToGlobalMappingCreateIS(osm->lis,<og);
399: PetscMalloc1(m_local,&idx2);
400: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx2);
401: ISLocalToGlobalMappingDestroy(<og);
402: if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of lis");
403: ISCreateGeneral(PETSC_COMM_SELF,m_local,idx2,PETSC_OWN_POINTER,&isll_local);
405: ISRestoreIndices(osm->is_local[i], &idx_local);
406: VecScatterCreate(osm->y[i],isll,osm->ly,isll_local,&osm->lprolongation[i]);
408: ISDestroy(&isll);
409: ISDestroy(&isll_local);
410: }
411: }
412: VecDestroy(&vec);
413: }
415: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
416: IS *cis;
417: PetscInt c;
419: PetscMalloc1(osm->n_local_true, &cis);
420: for (c = 0; c < osm->n_local_true; ++c) cis[c] = osm->lis;
421: MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, cis, scall, &osm->lmats);
422: PetscFree(cis);
423: }
425: /* Return control to the user so that the submatrices can be modified (e.g., to apply
426: different boundary conditions for the submatrices than for the global problem) */
427: PCModifySubMatrices(pc,osm->n_local_true,osm->is,osm->is,osm->pmat,pc->modifysubmatricesP);
429: /*
430: Loop over subdomains putting them into local ksp
431: */
432: for (i=0; i<osm->n_local_true; i++) {
433: KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i]);
434: if (!pc->setupcalled) {
435: KSPSetFromOptions(osm->ksp[i]);
436: }
437: }
438: return(0);
439: }
441: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
442: {
443: PC_ASM *osm = (PC_ASM*)pc->data;
444: PetscErrorCode ierr;
445: PetscInt i;
446: KSPConvergedReason reason;
449: for (i=0; i<osm->n_local_true; i++) {
450: KSPSetUp(osm->ksp[i]);
451: KSPGetConvergedReason(osm->ksp[i],&reason);
452: if (reason == KSP_DIVERGED_PC_FAILED) {
453: pc->failedreason = PC_SUBPC_ERROR;
454: }
455: }
456: return(0);
457: }
459: static PetscErrorCode PCApply_ASM(PC pc,Vec x,Vec y)
460: {
461: PC_ASM *osm = (PC_ASM*)pc->data;
463: PetscInt i,n_local_true = osm->n_local_true;
464: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
467: /*
468: Support for limiting the restriction or interpolation to only local
469: subdomain values (leaving the other values 0).
470: */
471: if (!(osm->type & PC_ASM_RESTRICT)) {
472: forward = SCATTER_FORWARD_LOCAL;
473: /* have to zero the work RHS since scatter may leave some slots empty */
474: VecSet(osm->lx, 0.0);
475: }
476: if (!(osm->type & PC_ASM_INTERPOLATE)) {
477: reverse = SCATTER_REVERSE_LOCAL;
478: }
480: if(osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE){
481: /* zero the global and the local solutions */
482: VecZeroEntries(y);
483: VecSet(osm->ly, 0.0);
485: /* Copy the global RHS to local RHS including the ghost nodes */
486: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
487: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
489: /* Restrict local RHS to the overlapping 0-block RHS */
490: VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
491: VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
493: /* do the local solves */
494: for (i = 0; i < n_local_true; ++i) {
496: /* solve the overlapping i-block */
497: PetscLogEventBegin(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
498: KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]);
499: KSPCheckSolve(osm->ksp[i],pc,osm->y[i]);
500: PetscLogEventEnd(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
502: if (osm->lprolongation) { /* interpolate the non-overalapping i-block solution to the local solution (only for restrictive additive) */
503: VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
504: VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
505: }
506: else{ /* interpolate the overalapping i-block solution to the local solution */
507: VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
508: VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
509: }
511: if (i < n_local_true-1) {
512: /* Restrict local RHS to the overlapping (i+1)-block RHS */
513: VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
514: VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
516: if ( osm->loctype == PC_COMPOSITE_MULTIPLICATIVE){
517: /* udpdate the overlapping (i+1)-block RHS using the current local solution */
518: MatMult(osm->lmats[i+1], osm->ly, osm->y[i+1]);
519: VecAXPBY(osm->x[i+1],-1.,1., osm->y[i+1]);
520: }
521: }
522: }
523: /* Add the local solution to the global solution including the ghost nodes */
524: VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
525: VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
526: }else{
527: SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
528: }
529: return(0);
530: }
532: static PetscErrorCode PCApplyTranspose_ASM(PC pc,Vec x,Vec y)
533: {
534: PC_ASM *osm = (PC_ASM*)pc->data;
536: PetscInt i,n_local_true = osm->n_local_true;
537: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
540: /*
541: Support for limiting the restriction or interpolation to only local
542: subdomain values (leaving the other values 0).
544: Note: these are reversed from the PCApply_ASM() because we are applying the
545: transpose of the three terms
546: */
548: if (!(osm->type & PC_ASM_INTERPOLATE)) {
549: forward = SCATTER_FORWARD_LOCAL;
550: /* have to zero the work RHS since scatter may leave some slots empty */
551: VecSet(osm->lx, 0.0);
552: }
553: if (!(osm->type & PC_ASM_RESTRICT)) reverse = SCATTER_REVERSE_LOCAL;
555: /* zero the global and the local solutions */
556: VecZeroEntries(y);
557: VecSet(osm->ly, 0.0);
559: /* Copy the global RHS to local RHS including the ghost nodes */
560: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
561: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
563: /* Restrict local RHS to the overlapping 0-block RHS */
564: VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
565: VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
567: /* do the local solves */
568: for (i = 0; i < n_local_true; ++i) {
570: /* solve the overlapping i-block */
571: PetscLogEventBegin(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
572: KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i]);
573: KSPCheckSolve(osm->ksp[i],pc,osm->y[i]);
574: PetscLogEventEnd(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
576: if (osm->lprolongation) { /* interpolate the non-overalapping i-block solution to the local solution */
577: VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
578: VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
579: }
580: else{ /* interpolate the overalapping i-block solution to the local solution */
581: VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
582: VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
583: }
585: if (i < n_local_true-1) {
586: /* Restrict local RHS to the overlapping (i+1)-block RHS */
587: VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
588: VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
589: }
590: }
591: /* Add the local solution to the global solution including the ghost nodes */
592: VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
593: VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
595: return(0);
597: }
599: static PetscErrorCode PCReset_ASM(PC pc)
600: {
601: PC_ASM *osm = (PC_ASM*)pc->data;
603: PetscInt i;
606: if (osm->ksp) {
607: for (i=0; i<osm->n_local_true; i++) {
608: KSPReset(osm->ksp[i]);
609: }
610: }
611: if (osm->pmat) {
612: if (osm->n_local_true > 0) {
613: MatDestroySubMatrices(osm->n_local_true,&osm->pmat);
614: }
615: }
616: if (osm->lrestriction) {
617: VecScatterDestroy(&osm->restriction);
618: for (i=0; i<osm->n_local_true; i++) {
619: VecScatterDestroy(&osm->lrestriction[i]);
620: if (osm->lprolongation) {VecScatterDestroy(&osm->lprolongation[i]);}
621: VecDestroy(&osm->x[i]);
622: VecDestroy(&osm->y[i]);
623: }
624: PetscFree(osm->lrestriction);
625: if (osm->lprolongation) {PetscFree(osm->lprolongation);}
626: PetscFree(osm->x);
627: PetscFree(osm->y);
629: }
630: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
631: ISDestroy(&osm->lis);
632: VecDestroy(&osm->lx);
633: VecDestroy(&osm->ly);
634: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
635: MatDestroyMatrices(osm->n_local_true, &osm->lmats);
636: }
638: PetscFree(osm->sub_mat_type);
640: osm->is = 0;
641: osm->is_local = 0;
642: return(0);
643: }
645: static PetscErrorCode PCDestroy_ASM(PC pc)
646: {
647: PC_ASM *osm = (PC_ASM*)pc->data;
649: PetscInt i;
652: PCReset_ASM(pc);
653: if (osm->ksp) {
654: for (i=0; i<osm->n_local_true; i++) {
655: KSPDestroy(&osm->ksp[i]);
656: }
657: PetscFree(osm->ksp);
658: }
659: PetscFree(pc->data);
660: return(0);
661: }
663: static PetscErrorCode PCSetFromOptions_ASM(PetscOptionItems *PetscOptionsObject,PC pc)
664: {
665: PC_ASM *osm = (PC_ASM*)pc->data;
667: PetscInt blocks,ovl;
668: PetscBool symset,flg;
669: PCASMType asmtype;
670: PCCompositeType loctype;
671: char sub_mat_type[256];
674: /* set the type to symmetric if matrix is symmetric */
675: if (!osm->type_set && pc->pmat) {
676: MatIsSymmetricKnown(pc->pmat,&symset,&flg);
677: if (symset && flg) osm->type = PC_ASM_BASIC;
678: }
679: PetscOptionsHead(PetscOptionsObject,"Additive Schwarz options");
680: PetscOptionsBool("-pc_asm_dm_subdomains","Use DMCreateDomainDecomposition() to define subdomains","PCASMSetDMSubdomains",osm->dm_subdomains,&osm->dm_subdomains,&flg);
681: PetscOptionsInt("-pc_asm_blocks","Number of subdomains","PCASMSetTotalSubdomains",osm->n,&blocks,&flg);
682: if (flg) {
683: PCASMSetTotalSubdomains(pc,blocks,NULL,NULL);
684: osm->dm_subdomains = PETSC_FALSE;
685: }
686: PetscOptionsInt("-pc_asm_local_blocks","Number of local subdomains","PCASMSetLocalSubdomains",osm->n_local_true,&blocks,&flg);
687: if (flg) {
688: PCASMSetLocalSubdomains(pc,blocks,NULL,NULL);
689: osm->dm_subdomains = PETSC_FALSE;
690: }
691: PetscOptionsInt("-pc_asm_overlap","Number of grid points overlap","PCASMSetOverlap",osm->overlap,&ovl,&flg);
692: if (flg) {
693: PCASMSetOverlap(pc,ovl);
694: osm->dm_subdomains = PETSC_FALSE;
695: }
696: flg = PETSC_FALSE;
697: PetscOptionsEnum("-pc_asm_type","Type of restriction/extension","PCASMSetType",PCASMTypes,(PetscEnum)osm->type,(PetscEnum*)&asmtype,&flg);
698: if (flg) {PCASMSetType(pc,asmtype); }
699: flg = PETSC_FALSE;
700: PetscOptionsEnum("-pc_asm_local_type","Type of local solver composition","PCASMSetLocalType",PCCompositeTypes,(PetscEnum)osm->loctype,(PetscEnum*)&loctype,&flg);
701: if (flg) {PCASMSetLocalType(pc,loctype); }
702: PetscOptionsFList("-pc_asm_sub_mat_type","Subsolve Matrix Type","PCASMSetSubMatType",MatList,NULL,sub_mat_type,256,&flg);
703: if(flg){
704: PCASMSetSubMatType(pc,sub_mat_type);
705: }
706: PetscOptionsTail();
707: return(0);
708: }
710: /*------------------------------------------------------------------------------------*/
712: static PetscErrorCode PCASMSetLocalSubdomains_ASM(PC pc,PetscInt n,IS is[],IS is_local[])
713: {
714: PC_ASM *osm = (PC_ASM*)pc->data;
716: PetscInt i;
719: if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Each process must have 1 or more blocks, n = %D",n);
720: if (pc->setupcalled && (n != osm->n_local_true || is)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetLocalSubdomains() should be called before calling PCSetUp().");
722: if (!pc->setupcalled) {
723: if (is) {
724: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is[i]);}
725: }
726: if (is_local) {
727: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is_local[i]);}
728: }
729: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
731: osm->n_local_true = n;
732: osm->is = 0;
733: osm->is_local = 0;
734: if (is) {
735: PetscMalloc1(n,&osm->is);
736: for (i=0; i<n; i++) osm->is[i] = is[i];
737: /* Flag indicating that the user has set overlapping subdomains so PCASM should not increase their size. */
738: osm->overlap = -1;
739: }
740: if (is_local) {
741: PetscMalloc1(n,&osm->is_local);
742: for (i=0; i<n; i++) osm->is_local[i] = is_local[i];
743: if (!is) {
744: PetscMalloc1(osm->n_local_true,&osm->is);
745: for (i=0; i<osm->n_local_true; i++) {
746: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
747: ISDuplicate(osm->is_local[i],&osm->is[i]);
748: ISCopy(osm->is_local[i],osm->is[i]);
749: } else {
750: PetscObjectReference((PetscObject)osm->is_local[i]);
751: osm->is[i] = osm->is_local[i];
752: }
753: }
754: }
755: }
756: }
757: return(0);
758: }
760: static PetscErrorCode PCASMSetTotalSubdomains_ASM(PC pc,PetscInt N,IS *is,IS *is_local)
761: {
762: PC_ASM *osm = (PC_ASM*)pc->data;
764: PetscMPIInt rank,size;
765: PetscInt n;
768: if (N < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of total blocks must be > 0, N = %D",N);
769: 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.");
771: /*
772: Split the subdomains equally among all processors
773: */
774: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
775: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
776: n = N/size + ((N % size) > rank);
777: 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);
778: if (pc->setupcalled && n != osm->n_local_true) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCASMSetTotalSubdomains() should be called before PCSetUp().");
779: if (!pc->setupcalled) {
780: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
782: osm->n_local_true = n;
783: osm->is = 0;
784: osm->is_local = 0;
785: }
786: return(0);
787: }
789: static PetscErrorCode PCASMSetOverlap_ASM(PC pc,PetscInt ovl)
790: {
791: PC_ASM *osm = (PC_ASM*)pc->data;
794: if (ovl < 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested");
795: if (pc->setupcalled && ovl != osm->overlap) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetOverlap() should be called before PCSetUp().");
796: if (!pc->setupcalled) osm->overlap = ovl;
797: return(0);
798: }
800: static PetscErrorCode PCASMSetType_ASM(PC pc,PCASMType type)
801: {
802: PC_ASM *osm = (PC_ASM*)pc->data;
805: osm->type = type;
806: osm->type_set = PETSC_TRUE;
807: return(0);
808: }
810: static PetscErrorCode PCASMGetType_ASM(PC pc,PCASMType *type)
811: {
812: PC_ASM *osm = (PC_ASM*)pc->data;
815: *type = osm->type;
816: return(0);
817: }
819: static PetscErrorCode PCASMSetLocalType_ASM(PC pc, PCCompositeType type)
820: {
821: PC_ASM *osm = (PC_ASM *) pc->data;
824: 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");
825: osm->loctype = type;
826: return(0);
827: }
829: static PetscErrorCode PCASMGetLocalType_ASM(PC pc, PCCompositeType *type)
830: {
831: PC_ASM *osm = (PC_ASM *) pc->data;
834: *type = osm->loctype;
835: return(0);
836: }
838: static PetscErrorCode PCASMSetSortIndices_ASM(PC pc,PetscBool doSort)
839: {
840: PC_ASM *osm = (PC_ASM*)pc->data;
843: osm->sort_indices = doSort;
844: return(0);
845: }
847: static PetscErrorCode PCASMGetSubKSP_ASM(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
848: {
849: PC_ASM *osm = (PC_ASM*)pc->data;
853: 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");
855: if (n_local) *n_local = osm->n_local_true;
856: if (first_local) {
857: MPI_Scan(&osm->n_local_true,first_local,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
858: *first_local -= osm->n_local_true;
859: }
860: if (ksp) {
861: /* Assume that local solves are now different; not necessarily
862: true though! This flag is used only for PCView_ASM() */
863: *ksp = osm->ksp;
864: osm->same_local_solves = PETSC_FALSE;
865: }
866: return(0);
867: }
869: static PetscErrorCode PCASMGetSubMatType_ASM(PC pc,MatType *sub_mat_type)
870: {
871: PC_ASM *osm = (PC_ASM*)pc->data;
876: *sub_mat_type = osm->sub_mat_type;
877: return(0);
878: }
880: static PetscErrorCode PCASMSetSubMatType_ASM(PC pc,MatType sub_mat_type)
881: {
882: PetscErrorCode ierr;
883: PC_ASM *osm = (PC_ASM*)pc->data;
887: PetscFree(osm->sub_mat_type);
888: PetscStrallocpy(sub_mat_type,(char**)&osm->sub_mat_type);
889: return(0);
890: }
892: /*@C
893: PCASMSetLocalSubdomains - Sets the local subdomains (for this processor only) for the additive Schwarz preconditioner.
895: Collective on PC
897: Input Parameters:
898: + pc - the preconditioner context
899: . n - the number of subdomains for this processor (default value = 1)
900: . is - the index set that defines the subdomains for this processor
901: (or NULL for PETSc to determine subdomains)
902: - is_local - the index sets that define the local part of the subdomains for this processor, not used unless PCASMType is PC_ASM_RESTRICT
903: (or NULL to not provide these)
905: Options Database Key:
906: To set the total number of subdomain blocks rather than specify the
907: index sets, use the option
908: . -pc_asm_local_blocks <blks> - Sets local blocks
910: Notes:
911: The IS numbering is in the parallel, global numbering of the vector for both is and is_local
913: By default the ASM preconditioner uses 1 block per processor.
915: Use PCASMSetTotalSubdomains() to set the subdomains for all processors.
917: If is_local is provided and PCASMType is PC_ASM_RESTRICT then the solution only over the is_local region is interpolated
918: back to form the global solution (this is the standard restricted additive Schwarz method)
920: If the is_local is provided and PCASMType is PC_ASM_INTERPOLATE or PC_ASM_NONE then an error is generated since there is
921: no code to handle that case.
923: Level: advanced
925: .keywords: PC, ASM, set, local, subdomains, additive Schwarz
927: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
928: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), PCASMType, PCASMSetType()
929: @*/
930: PetscErrorCode PCASMSetLocalSubdomains(PC pc,PetscInt n,IS is[],IS is_local[])
931: {
936: PetscTryMethod(pc,"PCASMSetLocalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,is,is_local));
937: return(0);
938: }
940: /*@C
941: PCASMSetTotalSubdomains - Sets the subdomains for all processors for the
942: additive Schwarz preconditioner. Either all or no processors in the
943: PC communicator must call this routine, with the same index sets.
945: Collective on PC
947: Input Parameters:
948: + pc - the preconditioner context
949: . N - the number of subdomains for all processors
950: . is - the index sets that define the subdomains for all processors
951: (or NULL to ask PETSc to determine the subdomains)
952: - is_local - the index sets that define the local part of the subdomains for this processor
953: (or NULL to not provide this information)
955: Options Database Key:
956: To set the total number of subdomain blocks rather than specify the
957: index sets, use the option
958: . -pc_asm_blocks <blks> - Sets total blocks
960: Notes:
961: Currently you cannot use this to set the actual subdomains with the argument is or is_local.
963: By default the ASM preconditioner uses 1 block per processor.
965: These index sets cannot be destroyed until after completion of the
966: linear solves for which the ASM preconditioner is being used.
968: Use PCASMSetLocalSubdomains() to set local subdomains.
970: The IS numbering is in the parallel, global numbering of the vector for both is and is_local
972: Level: advanced
974: .keywords: PC, ASM, set, total, global, subdomains, additive Schwarz
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: .keywords: PC, ASM, set, overlap
1028: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1029: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), MatIncreaseOverlap()
1030: @*/
1031: PetscErrorCode PCASMSetOverlap(PC pc,PetscInt ovl)
1032: {
1038: PetscTryMethod(pc,"PCASMSetOverlap_C",(PC,PetscInt),(pc,ovl));
1039: return(0);
1040: }
1042: /*@
1043: PCASMSetType - Sets the type of restriction and interpolation used
1044: for local problems in the additive Schwarz method.
1046: Logically Collective on PC
1048: Input Parameters:
1049: + pc - the preconditioner context
1050: - type - variant of ASM, one of
1051: .vb
1052: PC_ASM_BASIC - full interpolation and restriction
1053: PC_ASM_RESTRICT - full restriction, local processor interpolation
1054: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1055: PC_ASM_NONE - local processor restriction and interpolation
1056: .ve
1058: Options Database Key:
1059: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
1061: Notes:
1062: if the is_local arguments are passed to PCASMSetLocalSubdomains() then they are used when PC_ASM_RESTRICT has been selected
1063: to limit the local processor interpolation
1065: Level: intermediate
1067: .keywords: PC, ASM, set, type
1069: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1070: PCASMCreateSubdomains2D(), PCASMType, PCASMSetLocalType(), PCASMGetLocalType()
1071: @*/
1072: PetscErrorCode PCASMSetType(PC pc,PCASMType type)
1073: {
1079: PetscTryMethod(pc,"PCASMSetType_C",(PC,PCASMType),(pc,type));
1080: return(0);
1081: }
1083: /*@
1084: PCASMGetType - Gets the type of restriction and interpolation used
1085: for local problems in the additive Schwarz method.
1087: Logically Collective on PC
1089: Input Parameter:
1090: . pc - the preconditioner context
1092: Output Parameter:
1093: . type - variant of ASM, one of
1095: .vb
1096: PC_ASM_BASIC - full interpolation and restriction
1097: PC_ASM_RESTRICT - full restriction, local processor interpolation
1098: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1099: PC_ASM_NONE - local processor restriction and interpolation
1100: .ve
1102: Options Database Key:
1103: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
1105: Level: intermediate
1107: .keywords: PC, ASM, set, type
1109: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1110: PCASMCreateSubdomains2D(), PCASMType, PCASMSetType(), PCASMSetLocalType(), PCASMGetLocalType()
1111: @*/
1112: PetscErrorCode PCASMGetType(PC pc,PCASMType *type)
1113: {
1118: PetscUseMethod(pc,"PCASMGetType_C",(PC,PCASMType*),(pc,type));
1119: return(0);
1120: }
1122: /*@
1123: PCASMSetLocalType - Sets the type of composition used for local problems in the additive Schwarz method.
1125: Logically Collective on PC
1127: Input Parameters:
1128: + pc - the preconditioner context
1129: - type - type of composition, one of
1130: .vb
1131: PC_COMPOSITE_ADDITIVE - local additive combination
1132: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1133: .ve
1135: Options Database Key:
1136: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1138: Level: intermediate
1140: .seealso: PCASMSetType(), PCASMGetType(), PCASMGetLocalType(), PCASM, PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1141: @*/
1142: PetscErrorCode PCASMSetLocalType(PC pc, PCCompositeType type)
1143: {
1149: PetscTryMethod(pc, "PCASMSetLocalType_C", (PC, PCCompositeType), (pc, type));
1150: return(0);
1151: }
1153: /*@
1154: PCASMGetLocalType - Gets the type of composition used for local problems in the additive Schwarz method.
1156: Logically Collective on PC
1158: Input Parameter:
1159: . pc - the preconditioner context
1161: Output Parameter:
1162: . type - type of composition, one of
1163: .vb
1164: PC_COMPOSITE_ADDITIVE - local additive combination
1165: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1166: .ve
1168: Options Database Key:
1169: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1171: Level: intermediate
1173: .seealso: PCASMSetType(), PCASMGetType(), PCASMSetLocalType(), PCASMCreate(), PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1174: @*/
1175: PetscErrorCode PCASMGetLocalType(PC pc, PCCompositeType *type)
1176: {
1182: PetscUseMethod(pc, "PCASMGetLocalType_C", (PC, PCCompositeType *), (pc, type));
1183: return(0);
1184: }
1186: /*@
1187: PCASMSetSortIndices - Determines whether subdomain indices are sorted.
1189: Logically Collective on PC
1191: Input Parameters:
1192: + pc - the preconditioner context
1193: - doSort - sort the subdomain indices
1195: Level: intermediate
1197: .keywords: PC, ASM, set, type
1199: .seealso: PCASMSetLocalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1200: PCASMCreateSubdomains2D()
1201: @*/
1202: PetscErrorCode PCASMSetSortIndices(PC pc,PetscBool doSort)
1203: {
1209: PetscTryMethod(pc,"PCASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));
1210: return(0);
1211: }
1213: /*@C
1214: PCASMGetSubKSP - Gets the local KSP contexts for all blocks on
1215: this processor.
1217: Collective on PC iff first_local is requested
1219: Input Parameter:
1220: . pc - the preconditioner context
1222: Output Parameters:
1223: + n_local - the number of blocks on this processor or NULL
1224: . first_local - the global number of the first block on this processor or NULL,
1225: all processors must request or all must pass NULL
1226: - ksp - the array of KSP contexts
1228: Note:
1229: After PCASMGetSubKSP() the array of KSPes is not to be freed.
1231: You must call KSPSetUp() before calling PCASMGetSubKSP().
1233: Fortran note:
1234: 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.
1236: Level: advanced
1238: .keywords: PC, ASM, additive Schwarz, get, sub, KSP, context
1240: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap(),
1241: PCASMCreateSubdomains2D(),
1242: @*/
1243: PetscErrorCode PCASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
1244: {
1249: PetscUseMethod(pc,"PCASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
1250: return(0);
1251: }
1253: /* -------------------------------------------------------------------------------------*/
1254: /*MC
1255: PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1256: its own KSP object.
1258: Options Database Keys:
1259: + -pc_asm_blocks <blks> - Sets total blocks
1260: . -pc_asm_overlap <ovl> - Sets overlap
1261: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type, default is restrict
1262: - -pc_asm_local_type [additive, multiplicative] - Sets ASM type, default is additive
1264: IMPORTANT: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you
1265: will get a different convergence rate due to the default option of -pc_asm_type restrict. Use
1266: -pc_asm_type basic to use the standard ASM.
1268: Notes:
1269: Each processor can have one or more blocks, but a block cannot be shared by more
1270: than one processor. Use PCGASM for subdomains shared by multiple processes. Defaults to one block per processor.
1272: To set options on the solvers for each block append -sub_ to all the KSP, and PC
1273: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
1275: To set the options on the solvers separate for each block call PCASMGetSubKSP()
1276: and set the options directly on the resulting KSP object (you can access its PC
1277: with KSPGetPC())
1279: Level: beginner
1281: Concepts: additive Schwarz method
1283: References:
1284: + 1. - M Dryja, OB Widlund, An additive variant of the Schwarz alternating method for the case of many subregions
1285: Courant Institute, New York University Technical report
1286: - 1. - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations,
1287: Cambridge University Press.
1289: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
1290: PCBJACOBI, PCASMGetSubKSP(), PCASMSetLocalSubdomains(), PCASMType, PCASMGetType(), PCASMSetLocalType(), PCASMGetLocalType()
1291: PCASMSetTotalSubdomains(), PCSetModifySubmatrices(), PCASMSetOverlap(), PCASMSetType(), PCCompositeType
1293: M*/
1295: PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1296: {
1298: PC_ASM *osm;
1301: PetscNewLog(pc,&osm);
1303: osm->n = PETSC_DECIDE;
1304: osm->n_local = 0;
1305: osm->n_local_true = PETSC_DECIDE;
1306: osm->overlap = 1;
1307: osm->ksp = 0;
1308: osm->restriction = 0;
1309: osm->lprolongation = 0;
1310: osm->lrestriction = 0;
1311: osm->x = 0;
1312: osm->y = 0;
1313: osm->is = 0;
1314: osm->is_local = 0;
1315: osm->mat = 0;
1316: osm->pmat = 0;
1317: osm->type = PC_ASM_RESTRICT;
1318: osm->loctype = PC_COMPOSITE_ADDITIVE;
1319: osm->same_local_solves = PETSC_TRUE;
1320: osm->sort_indices = PETSC_TRUE;
1321: osm->dm_subdomains = PETSC_FALSE;
1322: osm->sub_mat_type = NULL;
1324: pc->data = (void*)osm;
1325: pc->ops->apply = PCApply_ASM;
1326: pc->ops->applytranspose = PCApplyTranspose_ASM;
1327: pc->ops->setup = PCSetUp_ASM;
1328: pc->ops->reset = PCReset_ASM;
1329: pc->ops->destroy = PCDestroy_ASM;
1330: pc->ops->setfromoptions = PCSetFromOptions_ASM;
1331: pc->ops->setuponblocks = PCSetUpOnBlocks_ASM;
1332: pc->ops->view = PCView_ASM;
1333: pc->ops->applyrichardson = 0;
1335: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",PCASMSetLocalSubdomains_ASM);
1336: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",PCASMSetTotalSubdomains_ASM);
1337: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetOverlap_C",PCASMSetOverlap_ASM);
1338: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetType_C",PCASMSetType_ASM);
1339: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetType_C",PCASMGetType_ASM);
1340: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalType_C",PCASMSetLocalType_ASM);
1341: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetLocalType_C",PCASMGetLocalType_ASM);
1342: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSortIndices_C",PCASMSetSortIndices_ASM);
1343: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubKSP_C",PCASMGetSubKSP_ASM);
1344: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubMatType_C",PCASMGetSubMatType_ASM);
1345: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSubMatType_C",PCASMSetSubMatType_ASM);
1346: return(0);
1347: }
1349: /*@C
1350: PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1351: preconditioner for a any problem on a general grid.
1353: Collective
1355: Input Parameters:
1356: + A - The global matrix operator
1357: - n - the number of local blocks
1359: Output Parameters:
1360: . outis - the array of index sets defining the subdomains
1362: Level: advanced
1364: Note: this generates nonoverlapping subdomains; the PCASM will generate the overlap
1365: from these if you use PCASMSetLocalSubdomains()
1367: In the Fortran version you must provide the array outis[] already allocated of length n.
1369: .keywords: PC, ASM, additive Schwarz, create, subdomains, unstructured grid
1371: .seealso: PCASMSetLocalSubdomains(), PCASMDestroySubdomains()
1372: @*/
1373: PetscErrorCode PCASMCreateSubdomains(Mat A, PetscInt n, IS* outis[])
1374: {
1375: MatPartitioning mpart;
1376: const char *prefix;
1377: PetscInt i,j,rstart,rend,bs;
1378: PetscBool hasop, isbaij = PETSC_FALSE,foundpart = PETSC_FALSE;
1379: Mat Ad = NULL, adj;
1380: IS ispart,isnumb,*is;
1381: PetscErrorCode ierr;
1386: if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of local blocks must be > 0, n = %D",n);
1388: /* Get prefix, row distribution, and block size */
1389: MatGetOptionsPrefix(A,&prefix);
1390: MatGetOwnershipRange(A,&rstart,&rend);
1391: MatGetBlockSize(A,&bs);
1392: 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);
1394: /* Get diagonal block from matrix if possible */
1395: MatHasOperation(A,MATOP_GET_DIAGONAL_BLOCK,&hasop);
1396: if (hasop) {
1397: MatGetDiagonalBlock(A,&Ad);
1398: }
1399: if (Ad) {
1400: PetscObjectTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij);
1401: if (!isbaij) {PetscObjectTypeCompare((PetscObject)Ad,MATSEQSBAIJ,&isbaij);}
1402: }
1403: if (Ad && n > 1) {
1404: PetscBool match,done;
1405: /* Try to setup a good matrix partitioning if available */
1406: MatPartitioningCreate(PETSC_COMM_SELF,&mpart);
1407: PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);
1408: MatPartitioningSetFromOptions(mpart);
1409: PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGCURRENT,&match);
1410: if (!match) {
1411: PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGSQUARE,&match);
1412: }
1413: if (!match) { /* assume a "good" partitioner is available */
1414: PetscInt na;
1415: const PetscInt *ia,*ja;
1416: MatGetRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1417: if (done) {
1418: /* Build adjacency matrix by hand. Unfortunately a call to
1419: MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1420: remove the block-aij structure and we cannot expect
1421: MatPartitioning to split vertices as we need */
1422: PetscInt i,j,len,nnz,cnt,*iia=0,*jja=0;
1423: const PetscInt *row;
1424: nnz = 0;
1425: for (i=0; i<na; i++) { /* count number of nonzeros */
1426: len = ia[i+1] - ia[i];
1427: row = ja + ia[i];
1428: for (j=0; j<len; j++) {
1429: if (row[j] == i) { /* don't count diagonal */
1430: len--; break;
1431: }
1432: }
1433: nnz += len;
1434: }
1435: PetscMalloc1(na+1,&iia);
1436: PetscMalloc1(nnz,&jja);
1437: nnz = 0;
1438: iia[0] = 0;
1439: for (i=0; i<na; i++) { /* fill adjacency */
1440: cnt = 0;
1441: len = ia[i+1] - ia[i];
1442: row = ja + ia[i];
1443: for (j=0; j<len; j++) {
1444: if (row[j] != i) { /* if not diagonal */
1445: jja[nnz+cnt++] = row[j];
1446: }
1447: }
1448: nnz += cnt;
1449: iia[i+1] = nnz;
1450: }
1451: /* Partitioning of the adjacency matrix */
1452: MatCreateMPIAdj(PETSC_COMM_SELF,na,na,iia,jja,NULL,&adj);
1453: MatPartitioningSetAdjacency(mpart,adj);
1454: MatPartitioningSetNParts(mpart,n);
1455: MatPartitioningApply(mpart,&ispart);
1456: ISPartitioningToNumbering(ispart,&isnumb);
1457: MatDestroy(&adj);
1458: foundpart = PETSC_TRUE;
1459: }
1460: MatRestoreRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1461: }
1462: MatPartitioningDestroy(&mpart);
1463: }
1465: PetscMalloc1(n,&is);
1466: *outis = is;
1468: if (!foundpart) {
1470: /* Partitioning by contiguous chunks of rows */
1472: PetscInt mbs = (rend-rstart)/bs;
1473: PetscInt start = rstart;
1474: for (i=0; i<n; i++) {
1475: PetscInt count = (mbs/n + ((mbs % n) > i)) * bs;
1476: ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);
1477: start += count;
1478: }
1480: } else {
1482: /* Partitioning by adjacency of diagonal block */
1484: const PetscInt *numbering;
1485: PetscInt *count,nidx,*indices,*newidx,start=0;
1486: /* Get node count in each partition */
1487: PetscMalloc1(n,&count);
1488: ISPartitioningCount(ispart,n,count);
1489: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1490: for (i=0; i<n; i++) count[i] *= bs;
1491: }
1492: /* Build indices from node numbering */
1493: ISGetLocalSize(isnumb,&nidx);
1494: PetscMalloc1(nidx,&indices);
1495: for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */
1496: ISGetIndices(isnumb,&numbering);
1497: PetscSortIntWithPermutation(nidx,numbering,indices);
1498: ISRestoreIndices(isnumb,&numbering);
1499: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1500: PetscMalloc1(nidx*bs,&newidx);
1501: for (i=0; i<nidx; i++) {
1502: for (j=0; j<bs; j++) newidx[i*bs+j] = indices[i]*bs + j;
1503: }
1504: PetscFree(indices);
1505: nidx *= bs;
1506: indices = newidx;
1507: }
1508: /* Shift to get global indices */
1509: for (i=0; i<nidx; i++) indices[i] += rstart;
1511: /* Build the index sets for each block */
1512: for (i=0; i<n; i++) {
1513: ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);
1514: ISSort(is[i]);
1515: start += count[i];
1516: }
1518: PetscFree(count);
1519: PetscFree(indices);
1520: ISDestroy(&isnumb);
1521: ISDestroy(&ispart);
1523: }
1524: return(0);
1525: }
1527: /*@C
1528: PCASMDestroySubdomains - Destroys the index sets created with
1529: PCASMCreateSubdomains(). Should be called after setting subdomains
1530: with PCASMSetLocalSubdomains().
1532: Collective
1534: Input Parameters:
1535: + n - the number of index sets
1536: . is - the array of index sets
1537: - is_local - the array of local index sets, can be NULL
1539: Level: advanced
1541: .keywords: PC, ASM, additive Schwarz, create, subdomains, unstructured grid
1543: .seealso: PCASMCreateSubdomains(), PCASMSetLocalSubdomains()
1544: @*/
1545: PetscErrorCode PCASMDestroySubdomains(PetscInt n, IS is[], IS is_local[])
1546: {
1547: PetscInt i;
1551: if (n <= 0) return(0);
1552: if (is) {
1554: for (i=0; i<n; i++) { ISDestroy(&is[i]); }
1555: PetscFree(is);
1556: }
1557: if (is_local) {
1559: for (i=0; i<n; i++) { ISDestroy(&is_local[i]); }
1560: PetscFree(is_local);
1561: }
1562: return(0);
1563: }
1565: /*@
1566: PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1567: preconditioner for a two-dimensional problem on a regular grid.
1569: Not Collective
1571: Input Parameters:
1572: + m, n - the number of mesh points in the x and y directions
1573: . M, N - the number of subdomains in the x and y directions
1574: . dof - degrees of freedom per node
1575: - overlap - overlap in mesh lines
1577: Output Parameters:
1578: + Nsub - the number of subdomains created
1579: . is - array of index sets defining overlapping (if overlap > 0) subdomains
1580: - is_local - array of index sets defining non-overlapping subdomains
1582: Note:
1583: Presently PCAMSCreateSubdomains2d() is valid only for sequential
1584: preconditioners. More general related routines are
1585: PCASMSetTotalSubdomains() and PCASMSetLocalSubdomains().
1587: Level: advanced
1589: .keywords: PC, ASM, additive Schwarz, create, subdomains, 2D, regular grid
1591: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1592: PCASMSetOverlap()
1593: @*/
1594: PetscErrorCode PCASMCreateSubdomains2D(PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt dof,PetscInt overlap,PetscInt *Nsub,IS **is,IS **is_local)
1595: {
1596: PetscInt i,j,height,width,ystart,xstart,yleft,yright,xleft,xright,loc_outer;
1598: PetscInt nidx,*idx,loc,ii,jj,count;
1601: if (dof != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP," ");
1603: *Nsub = N*M;
1604: PetscMalloc1(*Nsub,is);
1605: PetscMalloc1(*Nsub,is_local);
1606: ystart = 0;
1607: loc_outer = 0;
1608: for (i=0; i<N; i++) {
1609: height = n/N + ((n % N) > i); /* height of subdomain */
1610: if (height < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many N subdomains for mesh dimension n");
1611: yleft = ystart - overlap; if (yleft < 0) yleft = 0;
1612: yright = ystart + height + overlap; if (yright > n) yright = n;
1613: xstart = 0;
1614: for (j=0; j<M; j++) {
1615: width = m/M + ((m % M) > j); /* width of subdomain */
1616: if (width < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many M subdomains for mesh dimension m");
1617: xleft = xstart - overlap; if (xleft < 0) xleft = 0;
1618: xright = xstart + width + overlap; if (xright > m) xright = m;
1619: nidx = (xright - xleft)*(yright - yleft);
1620: PetscMalloc1(nidx,&idx);
1621: loc = 0;
1622: for (ii=yleft; ii<yright; ii++) {
1623: count = m*ii + xleft;
1624: for (jj=xleft; jj<xright; jj++) idx[loc++] = count++;
1625: }
1626: ISCreateGeneral(PETSC_COMM_SELF,nidx,idx,PETSC_COPY_VALUES,(*is)+loc_outer);
1627: if (overlap == 0) {
1628: PetscObjectReference((PetscObject)(*is)[loc_outer]);
1630: (*is_local)[loc_outer] = (*is)[loc_outer];
1631: } else {
1632: for (loc=0,ii=ystart; ii<ystart+height; ii++) {
1633: for (jj=xstart; jj<xstart+width; jj++) {
1634: idx[loc++] = m*ii + jj;
1635: }
1636: }
1637: ISCreateGeneral(PETSC_COMM_SELF,loc,idx,PETSC_COPY_VALUES,*is_local+loc_outer);
1638: }
1639: PetscFree(idx);
1640: xstart += width;
1641: loc_outer++;
1642: }
1643: ystart += height;
1644: }
1645: for (i=0; i<*Nsub; i++) { ISSort((*is)[i]); }
1646: return(0);
1647: }
1649: /*@C
1650: PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
1651: only) for the additive Schwarz preconditioner.
1653: Not Collective
1655: Input Parameter:
1656: . pc - the preconditioner context
1658: Output Parameters:
1659: + n - the number of subdomains for this processor (default value = 1)
1660: . is - the index sets that define the subdomains for this processor
1661: - is_local - the index sets that define the local part of the subdomains for this processor (can be NULL)
1664: Notes:
1665: The IS numbering is in the parallel, global numbering of the vector.
1667: Level: advanced
1669: .keywords: PC, ASM, set, local, subdomains, additive Schwarz
1671: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1672: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubmatrices()
1673: @*/
1674: PetscErrorCode PCASMGetLocalSubdomains(PC pc,PetscInt *n,IS *is[],IS *is_local[])
1675: {
1676: PC_ASM *osm = (PC_ASM*)pc->data;
1678: PetscBool match;
1684: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1685: if (!match) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PC is not a PCASM");
1686: if (n) *n = osm->n_local_true;
1687: if (is) *is = osm->is;
1688: if (is_local) *is_local = osm->is_local;
1689: return(0);
1690: }
1692: /*@C
1693: PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1694: only) for the additive Schwarz preconditioner.
1696: Not Collective
1698: Input Parameter:
1699: . pc - the preconditioner context
1701: Output Parameters:
1702: + n - the number of matrices for this processor (default value = 1)
1703: - mat - the matrices
1705: Level: advanced
1707: Notes:
1708: Call after PCSetUp() (or KSPSetUp()) but before PCApply() (or KSPApply()) and before PCSetUpOnBlocks())
1710: Usually one would use PCSetModifySubmatrices() to change the submatrices in building the preconditioner.
1712: .keywords: PC, ASM, set, local, subdomains, additive Schwarz, block Jacobi
1714: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1715: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains(), PCSetModifySubmatrices()
1716: @*/
1717: PetscErrorCode PCASMGetLocalSubmatrices(PC pc,PetscInt *n,Mat *mat[])
1718: {
1719: PC_ASM *osm;
1721: PetscBool match;
1727: if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUP() or PCSetUp().");
1728: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1729: if (!match) {
1730: if (n) *n = 0;
1731: if (mat) *mat = NULL;
1732: } else {
1733: osm = (PC_ASM*)pc->data;
1734: if (n) *n = osm->n_local_true;
1735: if (mat) *mat = osm->pmat;
1736: }
1737: return(0);
1738: }
1740: /*@
1741: PCASMSetDMSubdomains - Indicates whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1743: Logically Collective
1745: Input Parameter:
1746: + pc - the preconditioner
1747: - flg - boolean indicating whether to use subdomains defined by the DM
1749: Options Database Key:
1750: . -pc_asm_dm_subdomains
1752: Level: intermediate
1754: Notes:
1755: PCASMSetTotalSubdomains() and PCASMSetOverlap() take precedence over PCASMSetDMSubdomains(),
1756: so setting either of the first two effectively turns the latter off.
1758: .keywords: PC, ASM, DM, set, subdomains, additive Schwarz
1760: .seealso: PCASMGetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1761: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1762: @*/
1763: PetscErrorCode PCASMSetDMSubdomains(PC pc,PetscBool flg)
1764: {
1765: PC_ASM *osm = (PC_ASM*)pc->data;
1767: PetscBool match;
1772: if (pc->setupcalled) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for a setup PC.");
1773: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1774: if (match) {
1775: osm->dm_subdomains = flg;
1776: }
1777: return(0);
1778: }
1780: /*@
1781: PCASMGetDMSubdomains - Returns flag indicating whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1782: Not Collective
1784: Input Parameter:
1785: . pc - the preconditioner
1787: Output Parameter:
1788: . flg - boolean indicating whether to use subdomains defined by the DM
1790: Level: intermediate
1792: .keywords: PC, ASM, DM, set, subdomains, additive Schwarz
1794: .seealso: PCASMSetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1795: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1796: @*/
1797: PetscErrorCode PCASMGetDMSubdomains(PC pc,PetscBool* flg)
1798: {
1799: PC_ASM *osm = (PC_ASM*)pc->data;
1801: PetscBool match;
1806: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1807: if (match) {
1808: if (flg) *flg = osm->dm_subdomains;
1809: }
1810: return(0);
1811: }
1813: /*@
1814: PCASMGetSubMatType - Gets the matrix type used for ASM subsolves, as a string.
1816: Not Collective
1818: Input Parameter:
1819: . pc - the PC
1821: Output Parameter:
1822: . -pc_asm_sub_mat_type - name of matrix type
1824: Level: advanced
1826: .keywords: PC, PCASM, MatType, set
1828: .seealso: PCASMSetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1829: @*/
1830: PetscErrorCode PCASMGetSubMatType(PC pc,MatType *sub_mat_type){
1833: PetscTryMethod(pc,"PCASMGetSubMatType_C",(PC,MatType*),(pc,sub_mat_type));
1834: return(0);
1835: }
1837: /*@
1838: PCASMSetSubMatType - Set the type of matrix used for ASM subsolves
1840: Collective on Mat
1842: Input Parameters:
1843: + pc - the PC object
1844: - sub_mat_type - matrix type
1846: Options Database Key:
1847: . -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.
1849: Notes:
1850: See "${PETSC_DIR}/include/petscmat.h" for available types
1852: Level: advanced
1854: .keywords: PC, PCASM, MatType, set
1856: .seealso: PCASMGetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1857: @*/
1858: PetscErrorCode PCASMSetSubMatType(PC pc,MatType sub_mat_type)
1859: {
1862: PetscTryMethod(pc,"PCASMSetSubMatType_C",(PC,MatType),(pc,sub_mat_type));
1863: return(0);
1864: }