Actual source code: asm.c
petsc-3.14.6 2021-03-30
1: /*
2: This file defines an additive Schwarz preconditioner for any Mat implementation.
4: Note that each processor may have any number of subdomains. But in order to
5: deal easily with the VecScatter(), we treat each processor as if it has the
6: same number of subdomains.
8: n - total number of true subdomains on all processors
9: n_local_true - actual number of subdomains on this processor
10: n_local = maximum over all processors of n_local_true
11: */
13: #include <../src/ksp/pc/impls/asm/asm.h>
15: static PetscErrorCode PCView_ASM(PC pc,PetscViewer viewer)
16: {
17: PC_ASM *osm = (PC_ASM*)pc->data;
19: PetscMPIInt rank;
20: PetscInt i,bsz;
21: PetscBool iascii,isstring;
22: PetscViewer sviewer;
25: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
26: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
27: if (iascii) {
28: char overlaps[256] = "user-defined overlap",blocks[256] = "total subdomain blocks not yet set";
29: if (osm->overlap >= 0) {PetscSNPrintf(overlaps,sizeof(overlaps),"amount of overlap = %D",osm->overlap);}
30: if (osm->n > 0) {PetscSNPrintf(blocks,sizeof(blocks),"total subdomain blocks = %D",osm->n);}
31: PetscViewerASCIIPrintf(viewer," %s, %s\n",blocks,overlaps);
32: PetscViewerASCIIPrintf(viewer," restriction/interpolation type - %s\n",PCASMTypes[osm->type]);
33: if (osm->dm_subdomains) {PetscViewerASCIIPrintf(viewer," Additive Schwarz: using DM to define subdomains\n");}
34: if (osm->loctype != PC_COMPOSITE_ADDITIVE) {PetscViewerASCIIPrintf(viewer," Additive Schwarz: local solve composition type - %s\n",PCCompositeTypes[osm->loctype]);}
35: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
36: if (osm->same_local_solves) {
37: if (osm->ksp) {
38: PetscViewerASCIIPrintf(viewer," Local solver is the same for all blocks, as in the following KSP and PC objects on rank 0:\n");
39: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
40: if (!rank) {
41: PetscViewerASCIIPushTab(viewer);
42: KSPView(osm->ksp[0],sviewer);
43: PetscViewerASCIIPopTab(viewer);
44: }
45: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
46: }
47: } else {
48: PetscViewerASCIIPushSynchronized(viewer);
49: PetscViewerASCIISynchronizedPrintf(viewer," [%d] number of local blocks = %D\n",(int)rank,osm->n_local_true);
50: PetscViewerFlush(viewer);
51: PetscViewerASCIIPrintf(viewer," Local solve info for each block is in the following KSP and PC objects:\n");
52: PetscViewerASCIIPushTab(viewer);
53: PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
54: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
55: for (i=0; i<osm->n_local_true; i++) {
56: ISGetLocalSize(osm->is[i],&bsz);
57: PetscViewerASCIISynchronizedPrintf(sviewer,"[%d] local block number %D, size = %D\n",(int)rank,i,bsz);
58: KSPView(osm->ksp[i],sviewer);
59: PetscViewerASCIISynchronizedPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
60: }
61: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
62: PetscViewerASCIIPopTab(viewer);
63: PetscViewerFlush(viewer);
64: PetscViewerASCIIPopSynchronized(viewer);
65: }
66: } else if (isstring) {
67: PetscViewerStringSPrintf(viewer," blocks=%D, overlap=%D, type=%s",osm->n,osm->overlap,PCASMTypes[osm->type]);
68: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
69: if (osm->ksp) {KSPView(osm->ksp[0],sviewer);}
70: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
71: }
72: return(0);
73: }
75: static PetscErrorCode PCASMPrintSubdomains(PC pc)
76: {
77: PC_ASM *osm = (PC_ASM*)pc->data;
78: const char *prefix;
79: char fname[PETSC_MAX_PATH_LEN+1];
80: PetscViewer viewer, sviewer;
81: char *s;
82: PetscInt i,j,nidx;
83: const PetscInt *idx;
84: PetscMPIInt rank, size;
88: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
89: MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
90: PCGetOptionsPrefix(pc,&prefix);
91: PetscOptionsGetString(NULL,prefix,"-pc_asm_print_subdomains",fname,sizeof(fname),NULL);
92: if (fname[0] == 0) { PetscStrcpy(fname,"stdout"); };
93: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),fname,&viewer);
94: for (i=0; i<osm->n_local; i++) {
95: if (i < osm->n_local_true) {
96: ISGetLocalSize(osm->is[i],&nidx);
97: ISGetIndices(osm->is[i],&idx);
98: /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
99: #define len 16*(nidx+1)+512
100: PetscMalloc1(len,&s);
101: PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer);
102: #undef len
103: PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D with overlap:\n", rank, size, i);
104: for (j=0; j<nidx; j++) {
105: PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
106: }
107: ISRestoreIndices(osm->is[i],&idx);
108: PetscViewerStringSPrintf(sviewer,"\n");
109: PetscViewerDestroy(&sviewer);
110: PetscViewerASCIIPushSynchronized(viewer);
111: PetscViewerASCIISynchronizedPrintf(viewer, s);
112: PetscViewerFlush(viewer);
113: PetscViewerASCIIPopSynchronized(viewer);
114: PetscFree(s);
115: if (osm->is_local) {
116: /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
117: #define len 16*(nidx+1)+512
118: PetscMalloc1(len, &s);
119: PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer);
120: #undef len
121: PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D without overlap:\n", rank, size, i);
122: ISGetLocalSize(osm->is_local[i],&nidx);
123: ISGetIndices(osm->is_local[i],&idx);
124: for (j=0; j<nidx; j++) {
125: PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
126: }
127: ISRestoreIndices(osm->is_local[i],&idx);
128: PetscViewerStringSPrintf(sviewer,"\n");
129: PetscViewerDestroy(&sviewer);
130: PetscViewerASCIIPushSynchronized(viewer);
131: PetscViewerASCIISynchronizedPrintf(viewer, s);
132: PetscViewerFlush(viewer);
133: PetscViewerASCIIPopSynchronized(viewer);
134: PetscFree(s);
135: }
136: } else {
137: /* Participate in collective viewer calls. */
138: PetscViewerASCIIPushSynchronized(viewer);
139: PetscViewerFlush(viewer);
140: PetscViewerASCIIPopSynchronized(viewer);
141: /* Assume either all ranks have is_local or none do. */
142: if (osm->is_local) {
143: PetscViewerASCIIPushSynchronized(viewer);
144: PetscViewerFlush(viewer);
145: PetscViewerASCIIPopSynchronized(viewer);
146: }
147: }
148: }
149: PetscViewerFlush(viewer);
150: PetscViewerDestroy(&viewer);
151: return(0);
152: }
154: static PetscErrorCode PCSetUp_ASM(PC pc)
155: {
156: PC_ASM *osm = (PC_ASM*)pc->data;
158: PetscBool flg;
159: PetscInt i,m,m_local;
160: MatReuse scall = MAT_REUSE_MATRIX;
161: IS isl;
162: KSP ksp;
163: PC subpc;
164: const char *prefix,*pprefix;
165: Vec vec;
166: DM *domain_dm = NULL;
169: if (!pc->setupcalled) {
170: PetscInt m;
172: /* Note: if subdomains have been set either via PCASMSetTotalSubdomains() or via PCASMSetLocalSubdomains(), osm->n_local_true will not be PETSC_DECIDE */
173: if (osm->n_local_true == PETSC_DECIDE) {
174: /* no subdomains given */
175: /* try pc->dm first, if allowed */
176: if (osm->dm_subdomains && pc->dm) {
177: PetscInt num_domains, d;
178: char **domain_names;
179: IS *inner_domain_is, *outer_domain_is;
180: DMCreateDomainDecomposition(pc->dm, &num_domains, &domain_names, &inner_domain_is, &outer_domain_is, &domain_dm);
181: osm->overlap = -1; /* We do not want to increase the overlap of the IS.
182: A future improvement of this code might allow one to use
183: DM-defined subdomains and also increase the overlap,
184: but that is not currently supported */
185: if (num_domains) {
186: PCASMSetLocalSubdomains(pc, num_domains, outer_domain_is, inner_domain_is);
187: }
188: for (d = 0; d < num_domains; ++d) {
189: if (domain_names) {PetscFree(domain_names[d]);}
190: if (inner_domain_is) {ISDestroy(&inner_domain_is[d]);}
191: if (outer_domain_is) {ISDestroy(&outer_domain_is[d]);}
192: }
193: PetscFree(domain_names);
194: PetscFree(inner_domain_is);
195: PetscFree(outer_domain_is);
196: }
197: if (osm->n_local_true == PETSC_DECIDE) {
198: /* still no subdomains; use one subdomain per processor */
199: osm->n_local_true = 1;
200: }
201: }
202: { /* determine the global and max number of subdomains */
203: struct {PetscInt max,sum;} inwork,outwork;
204: PetscMPIInt size;
206: inwork.max = osm->n_local_true;
207: inwork.sum = osm->n_local_true;
208: MPIU_Allreduce(&inwork,&outwork,1,MPIU_2INT,MPIU_MAXSUM_OP,PetscObjectComm((PetscObject)pc));
209: osm->n_local = outwork.max;
210: osm->n = outwork.sum;
212: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
213: if (outwork.max == 1 && outwork.sum == size) {
214: /* osm->n_local_true = 1 on all processes, set this option may enable use of optimized MatCreateSubMatrices() implementation */
215: MatSetOption(pc->pmat,MAT_SUBMAT_SINGLEIS,PETSC_TRUE);
216: }
217: }
218: if (!osm->is) { /* create the index sets */
219: PCASMCreateSubdomains(pc->pmat,osm->n_local_true,&osm->is);
220: }
221: if (osm->n_local_true > 1 && !osm->is_local) {
222: PetscMalloc1(osm->n_local_true,&osm->is_local);
223: for (i=0; i<osm->n_local_true; i++) {
224: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
225: ISDuplicate(osm->is[i],&osm->is_local[i]);
226: ISCopy(osm->is[i],osm->is_local[i]);
227: } else {
228: PetscObjectReference((PetscObject)osm->is[i]);
229: osm->is_local[i] = osm->is[i];
230: }
231: }
232: }
233: PCGetOptionsPrefix(pc,&prefix);
234: flg = PETSC_FALSE;
235: PetscOptionsGetBool(NULL,prefix,"-pc_asm_print_subdomains",&flg,NULL);
236: if (flg) { PCASMPrintSubdomains(pc); }
238: if (osm->overlap > 0) {
239: /* Extend the "overlapping" regions by a number of steps */
240: MatIncreaseOverlap(pc->pmat,osm->n_local_true,osm->is,osm->overlap);
241: }
242: if (osm->sort_indices) {
243: for (i=0; i<osm->n_local_true; i++) {
244: ISSort(osm->is[i]);
245: if (osm->is_local) {
246: ISSort(osm->is_local[i]);
247: }
248: }
249: }
251: if (!osm->ksp) {
252: /* Create the local solvers */
253: PetscMalloc1(osm->n_local_true,&osm->ksp);
254: if (domain_dm) {
255: PetscInfo(pc,"Setting up ASM subproblems using the embedded DM\n");
256: }
257: for (i=0; i<osm->n_local_true; i++) {
258: KSPCreate(PETSC_COMM_SELF,&ksp);
259: KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
260: PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
261: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
262: KSPSetType(ksp,KSPPREONLY);
263: KSPGetPC(ksp,&subpc);
264: PCGetOptionsPrefix(pc,&prefix);
265: KSPSetOptionsPrefix(ksp,prefix);
266: KSPAppendOptionsPrefix(ksp,"sub_");
267: if (domain_dm) {
268: KSPSetDM(ksp, domain_dm[i]);
269: KSPSetDMActive(ksp, PETSC_FALSE);
270: DMDestroy(&domain_dm[i]);
271: }
272: osm->ksp[i] = ksp;
273: }
274: if (domain_dm) {
275: PetscFree(domain_dm);
276: }
277: }
279: ISConcatenate(PETSC_COMM_SELF, osm->n_local_true, osm->is, &osm->lis);
280: ISSortRemoveDups(osm->lis);
281: ISGetLocalSize(osm->lis, &m);
283: scall = MAT_INITIAL_MATRIX;
284: } else {
285: /*
286: Destroy the blocks from the previous iteration
287: */
288: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
289: MatDestroyMatrices(osm->n_local_true,&osm->pmat);
290: scall = MAT_INITIAL_MATRIX;
291: }
292: }
294: /* Destroy previous submatrices of a different type than pc->pmat since MAT_REUSE_MATRIX won't work in that case */
295: if ((scall == MAT_REUSE_MATRIX) && osm->sub_mat_type) {
296: if (osm->n_local_true > 0) {
297: MatDestroySubMatrices(osm->n_local_true,&osm->pmat);
298: }
299: scall = MAT_INITIAL_MATRIX;
300: }
302: /*
303: Extract out the submatrices
304: */
305: MatCreateSubMatrices(pc->pmat,osm->n_local_true,osm->is,osm->is,scall,&osm->pmat);
306: if (scall == MAT_INITIAL_MATRIX) {
307: PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
308: for (i=0; i<osm->n_local_true; i++) {
309: PetscLogObjectParent((PetscObject)pc,(PetscObject)osm->pmat[i]);
310: PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);
311: }
312: }
314: /* Convert the types of the submatrices (if needbe) */
315: if (osm->sub_mat_type) {
316: for (i=0; i<osm->n_local_true; i++) {
317: MatConvert(osm->pmat[i],osm->sub_mat_type,MAT_INPLACE_MATRIX,&(osm->pmat[i]));
318: }
319: }
321: if (!pc->setupcalled) {
322: VecType vtype;
324: /* Create the local work vectors (from the local matrices) and scatter contexts */
325: MatCreateVecs(pc->pmat,&vec,NULL);
327: 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()");
328: if (osm->is_local && osm->type == PC_ASM_RESTRICT && osm->loctype == PC_COMPOSITE_ADDITIVE) {
329: PetscMalloc1(osm->n_local_true,&osm->lprolongation);
330: }
331: PetscMalloc1(osm->n_local_true,&osm->lrestriction);
332: PetscMalloc1(osm->n_local_true,&osm->x);
333: PetscMalloc1(osm->n_local_true,&osm->y);
335: ISGetLocalSize(osm->lis,&m);
336: ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
337: MatGetVecType(osm->pmat[0],&vtype);
338: VecCreate(PETSC_COMM_SELF,&osm->lx);
339: VecSetSizes(osm->lx,m,m);
340: VecSetType(osm->lx,vtype);
341: VecDuplicate(osm->lx, &osm->ly);
342: VecScatterCreate(vec,osm->lis,osm->lx,isl,&osm->restriction);
343: ISDestroy(&isl);
345: for (i=0; i<osm->n_local_true; ++i) {
346: ISLocalToGlobalMapping ltog;
347: IS isll;
348: const PetscInt *idx_is;
349: PetscInt *idx_lis,nout;
351: ISGetLocalSize(osm->is[i],&m);
352: MatCreateVecs(osm->pmat[i],&osm->x[i],NULL);
353: VecDuplicate(osm->x[i],&osm->y[i]);
355: /* generate a scatter from ly to y[i] picking all the overlapping is[i] entries */
356: ISLocalToGlobalMappingCreateIS(osm->lis,<og);
357: ISGetLocalSize(osm->is[i],&m);
358: ISGetIndices(osm->is[i], &idx_is);
359: PetscMalloc1(m,&idx_lis);
360: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m,idx_is,&nout,idx_lis);
361: if (nout != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is not a subset of lis");
362: ISRestoreIndices(osm->is[i], &idx_is);
363: ISCreateGeneral(PETSC_COMM_SELF,m,idx_lis,PETSC_OWN_POINTER,&isll);
364: ISLocalToGlobalMappingDestroy(<og);
365: ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
366: VecScatterCreate(osm->ly,isll,osm->y[i],isl,&osm->lrestriction[i]);
367: ISDestroy(&isll);
368: ISDestroy(&isl);
369: if (osm->lprolongation) { /* generate a scatter from y[i] to ly picking only the the non-overlapping is_local[i] entries */
370: ISLocalToGlobalMapping ltog;
371: IS isll,isll_local;
372: const PetscInt *idx_local;
373: PetscInt *idx1, *idx2, nout;
375: ISGetLocalSize(osm->is_local[i],&m_local);
376: ISGetIndices(osm->is_local[i], &idx_local);
378: ISLocalToGlobalMappingCreateIS(osm->is[i],<og);
379: PetscMalloc1(m_local,&idx1);
380: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx1);
381: ISLocalToGlobalMappingDestroy(<og);
382: if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of is");
383: ISCreateGeneral(PETSC_COMM_SELF,m_local,idx1,PETSC_OWN_POINTER,&isll);
385: ISLocalToGlobalMappingCreateIS(osm->lis,<og);
386: PetscMalloc1(m_local,&idx2);
387: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx2);
388: ISLocalToGlobalMappingDestroy(<og);
389: if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of lis");
390: ISCreateGeneral(PETSC_COMM_SELF,m_local,idx2,PETSC_OWN_POINTER,&isll_local);
392: ISRestoreIndices(osm->is_local[i], &idx_local);
393: VecScatterCreate(osm->y[i],isll,osm->ly,isll_local,&osm->lprolongation[i]);
395: ISDestroy(&isll);
396: ISDestroy(&isll_local);
397: }
398: }
399: VecDestroy(&vec);
400: }
402: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
403: IS *cis;
404: PetscInt c;
406: PetscMalloc1(osm->n_local_true, &cis);
407: for (c = 0; c < osm->n_local_true; ++c) cis[c] = osm->lis;
408: MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, cis, scall, &osm->lmats);
409: PetscFree(cis);
410: }
412: /* Return control to the user so that the submatrices can be modified (e.g., to apply
413: different boundary conditions for the submatrices than for the global problem) */
414: PCModifySubMatrices(pc,osm->n_local_true,osm->is,osm->is,osm->pmat,pc->modifysubmatricesP);
416: /*
417: Loop over subdomains putting them into local ksp
418: */
419: for (i=0; i<osm->n_local_true; i++) {
420: KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i]);
421: if (!pc->setupcalled) {
422: KSPSetFromOptions(osm->ksp[i]);
423: }
424: }
425: return(0);
426: }
428: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
429: {
430: PC_ASM *osm = (PC_ASM*)pc->data;
431: PetscErrorCode ierr;
432: PetscInt i;
433: KSPConvergedReason reason;
436: for (i=0; i<osm->n_local_true; i++) {
437: KSPSetUp(osm->ksp[i]);
438: KSPGetConvergedReason(osm->ksp[i],&reason);
439: if (reason == KSP_DIVERGED_PC_FAILED) {
440: pc->failedreason = PC_SUBPC_ERROR;
441: }
442: }
443: return(0);
444: }
446: static PetscErrorCode PCApply_ASM(PC pc,Vec x,Vec y)
447: {
448: PC_ASM *osm = (PC_ASM*)pc->data;
450: PetscInt i,n_local_true = osm->n_local_true;
451: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
454: /*
455: support for limiting the restriction or interpolation to only local
456: subdomain values (leaving the other values 0).
457: */
458: if (!(osm->type & PC_ASM_RESTRICT)) {
459: forward = SCATTER_FORWARD_LOCAL;
460: /* have to zero the work RHS since scatter may leave some slots empty */
461: VecSet(osm->lx, 0.0);
462: }
463: if (!(osm->type & PC_ASM_INTERPOLATE)) {
464: reverse = SCATTER_REVERSE_LOCAL;
465: }
467: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE) {
468: /* zero the global and the local solutions */
469: VecSet(y, 0.0);
470: VecSet(osm->ly, 0.0);
472: /* copy the global RHS to local RHS including the ghost nodes */
473: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
474: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
476: /* restrict local RHS to the overlapping 0-block RHS */
477: VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
478: VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
480: /* do the local solves */
481: for (i = 0; i < n_local_true; ++i) {
483: /* solve the overlapping i-block */
484: PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i],0);
485: KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]);
486: KSPCheckSolve(osm->ksp[i], pc, osm->y[i]);
487: PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0);
489: if (osm->lprolongation) { /* interpolate the non-overlapping i-block solution to the local solution (only for restrictive additive) */
490: VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
491: VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
492: } else { /* interpolate the overlapping i-block solution to the local solution */
493: VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
494: VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
495: }
497: if (i < n_local_true-1) {
498: /* restrict local RHS to the overlapping (i+1)-block RHS */
499: VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
500: VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
502: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
503: /* update the overlapping (i+1)-block RHS using the current local solution */
504: MatMult(osm->lmats[i+1], osm->ly, osm->y[i+1]);
505: VecAXPBY(osm->x[i+1],-1.,1., osm->y[i+1]);
506: }
507: }
508: }
509: /* add the local solution to the global solution including the ghost nodes */
510: VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
511: VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
512: } else SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
513: return(0);
514: }
516: static PetscErrorCode PCMatApply_ASM(PC pc,Mat X,Mat Y)
517: {
518: PC_ASM *osm = (PC_ASM*)pc->data;
519: Mat Z,W;
520: Vec x;
521: PetscInt i,m,N;
522: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
526: if (osm->n_local_true > 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Not yet implemented");
527: /*
528: support for limiting the restriction or interpolation to only local
529: subdomain values (leaving the other values 0).
530: */
531: if (!(osm->type & PC_ASM_RESTRICT)) {
532: forward = SCATTER_FORWARD_LOCAL;
533: /* have to zero the work RHS since scatter may leave some slots empty */
534: VecSet(osm->lx, 0.0);
535: }
536: if (!(osm->type & PC_ASM_INTERPOLATE)) {
537: reverse = SCATTER_REVERSE_LOCAL;
538: }
539: VecGetLocalSize(osm->x[0], &m);
540: MatGetSize(X, NULL, &N);
541: MatCreateSeqDense(PETSC_COMM_SELF, m, N, NULL, &Z);
542: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE) {
543: /* zero the global and the local solutions */
544: MatZeroEntries(Y);
545: VecSet(osm->ly, 0.0);
547: for (i = 0; i < N; ++i) {
548: MatDenseGetColumnVecRead(X, i, &x);
549: /* copy the global RHS to local RHS including the ghost nodes */
550: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
551: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
552: MatDenseRestoreColumnVecRead(X, i, &x);
554: MatDenseGetColumnVecWrite(Z, i, &x);
555: /* restrict local RHS to the overlapping 0-block RHS */
556: VecScatterBegin(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward);
557: VecScatterEnd(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward);
558: MatDenseRestoreColumnVecWrite(Z, i, &x);
559: }
560: MatCreateSeqDense(PETSC_COMM_SELF, m, N, NULL, &W);
561: /* solve the overlapping 0-block */
562: PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[0], Z, W, 0);
563: KSPMatSolve(osm->ksp[0], Z, W);
564: KSPCheckSolve(osm->ksp[0], pc, NULL);
565: PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[0], Z, W,0);
566: MatDestroy(&Z);
568: for (i = 0; i < N; ++i) {
569: VecSet(osm->ly, 0.0);
570: MatDenseGetColumnVecRead(W, i, &x);
571: if (osm->lprolongation) { /* interpolate the non-overlapping 0-block solution to the local solution (only for restrictive additive) */
572: VecScatterBegin(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward);
573: VecScatterEnd(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward);
574: } else { /* interpolate the overlapping 0-block solution to the local solution */
575: VecScatterBegin(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse);
576: VecScatterEnd(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse);
577: }
578: MatDenseRestoreColumnVecRead(W, i, &x);
580: MatDenseGetColumnVecWrite(Y, i, &x);
581: /* add the local solution to the global solution including the ghost nodes */
582: VecScatterBegin(osm->restriction, osm->ly, x, ADD_VALUES, reverse);
583: VecScatterEnd(osm->restriction, osm->ly, x, ADD_VALUES, reverse);
584: MatDenseRestoreColumnVecWrite(Y, i, &x);
585: }
586: MatDestroy(&W);
587: } else SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
588: return(0);
589: }
591: static PetscErrorCode PCApplyTranspose_ASM(PC pc,Vec x,Vec y)
592: {
593: PC_ASM *osm = (PC_ASM*)pc->data;
595: PetscInt i,n_local_true = osm->n_local_true;
596: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
599: /*
600: Support for limiting the restriction or interpolation to only local
601: subdomain values (leaving the other values 0).
603: Note: these are reversed from the PCApply_ASM() because we are applying the
604: transpose of the three terms
605: */
607: if (!(osm->type & PC_ASM_INTERPOLATE)) {
608: forward = SCATTER_FORWARD_LOCAL;
609: /* have to zero the work RHS since scatter may leave some slots empty */
610: VecSet(osm->lx, 0.0);
611: }
612: if (!(osm->type & PC_ASM_RESTRICT)) reverse = SCATTER_REVERSE_LOCAL;
614: /* zero the global and the local solutions */
615: VecSet(y, 0.0);
616: VecSet(osm->ly, 0.0);
618: /* Copy the global RHS to local RHS including the ghost nodes */
619: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
620: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
622: /* Restrict local RHS to the overlapping 0-block RHS */
623: VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
624: VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
626: /* do the local solves */
627: for (i = 0; i < n_local_true; ++i) {
629: /* solve the overlapping i-block */
630: PetscLogEventBegin(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
631: KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i]);
632: KSPCheckSolve(osm->ksp[i],pc,osm->y[i]);
633: PetscLogEventEnd(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
635: if (osm->lprolongation) { /* interpolate the non-overlapping i-block solution to the local solution */
636: VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
637: VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
638: } else { /* interpolate the overlapping i-block solution to the local solution */
639: VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
640: VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
641: }
643: if (i < n_local_true-1) {
644: /* Restrict local RHS to the overlapping (i+1)-block RHS */
645: VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
646: VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
647: }
648: }
649: /* Add the local solution to the global solution including the ghost nodes */
650: VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
651: VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
653: return(0);
655: }
657: static PetscErrorCode PCReset_ASM(PC pc)
658: {
659: PC_ASM *osm = (PC_ASM*)pc->data;
661: PetscInt i;
664: if (osm->ksp) {
665: for (i=0; i<osm->n_local_true; i++) {
666: KSPReset(osm->ksp[i]);
667: }
668: }
669: if (osm->pmat) {
670: if (osm->n_local_true > 0) {
671: MatDestroySubMatrices(osm->n_local_true,&osm->pmat);
672: }
673: }
674: if (osm->lrestriction) {
675: VecScatterDestroy(&osm->restriction);
676: for (i=0; i<osm->n_local_true; i++) {
677: VecScatterDestroy(&osm->lrestriction[i]);
678: if (osm->lprolongation) {VecScatterDestroy(&osm->lprolongation[i]);}
679: VecDestroy(&osm->x[i]);
680: VecDestroy(&osm->y[i]);
681: }
682: PetscFree(osm->lrestriction);
683: if (osm->lprolongation) {PetscFree(osm->lprolongation);}
684: PetscFree(osm->x);
685: PetscFree(osm->y);
687: }
688: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
689: ISDestroy(&osm->lis);
690: VecDestroy(&osm->lx);
691: VecDestroy(&osm->ly);
692: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
693: MatDestroyMatrices(osm->n_local_true, &osm->lmats);
694: }
696: PetscFree(osm->sub_mat_type);
698: osm->is = NULL;
699: osm->is_local = NULL;
700: return(0);
701: }
703: static PetscErrorCode PCDestroy_ASM(PC pc)
704: {
705: PC_ASM *osm = (PC_ASM*)pc->data;
707: PetscInt i;
710: PCReset_ASM(pc);
711: if (osm->ksp) {
712: for (i=0; i<osm->n_local_true; i++) {
713: KSPDestroy(&osm->ksp[i]);
714: }
715: PetscFree(osm->ksp);
716: }
717: PetscFree(pc->data);
718: return(0);
719: }
721: static PetscErrorCode PCSetFromOptions_ASM(PetscOptionItems *PetscOptionsObject,PC pc)
722: {
723: PC_ASM *osm = (PC_ASM*)pc->data;
725: PetscInt blocks,ovl;
726: PetscBool flg;
727: PCASMType asmtype;
728: PCCompositeType loctype;
729: char sub_mat_type[256];
732: PetscOptionsHead(PetscOptionsObject,"Additive Schwarz options");
733: PetscOptionsBool("-pc_asm_dm_subdomains","Use DMCreateDomainDecomposition() to define subdomains","PCASMSetDMSubdomains",osm->dm_subdomains,&osm->dm_subdomains,&flg);
734: PetscOptionsInt("-pc_asm_blocks","Number of subdomains","PCASMSetTotalSubdomains",osm->n,&blocks,&flg);
735: if (flg) {
736: PCASMSetTotalSubdomains(pc,blocks,NULL,NULL);
737: osm->dm_subdomains = PETSC_FALSE;
738: }
739: PetscOptionsInt("-pc_asm_local_blocks","Number of local subdomains","PCASMSetLocalSubdomains",osm->n_local_true,&blocks,&flg);
740: if (flg) {
741: PCASMSetLocalSubdomains(pc,blocks,NULL,NULL);
742: osm->dm_subdomains = PETSC_FALSE;
743: }
744: PetscOptionsInt("-pc_asm_overlap","Number of grid points overlap","PCASMSetOverlap",osm->overlap,&ovl,&flg);
745: if (flg) {
746: PCASMSetOverlap(pc,ovl);
747: osm->dm_subdomains = PETSC_FALSE;
748: }
749: flg = PETSC_FALSE;
750: PetscOptionsEnum("-pc_asm_type","Type of restriction/extension","PCASMSetType",PCASMTypes,(PetscEnum)osm->type,(PetscEnum*)&asmtype,&flg);
751: if (flg) {PCASMSetType(pc,asmtype); }
752: flg = PETSC_FALSE;
753: PetscOptionsEnum("-pc_asm_local_type","Type of local solver composition","PCASMSetLocalType",PCCompositeTypes,(PetscEnum)osm->loctype,(PetscEnum*)&loctype,&flg);
754: if (flg) {PCASMSetLocalType(pc,loctype); }
755: PetscOptionsFList("-pc_asm_sub_mat_type","Subsolve Matrix Type","PCASMSetSubMatType",MatList,NULL,sub_mat_type,256,&flg);
756: if (flg) {
757: PCASMSetSubMatType(pc,sub_mat_type);
758: }
759: PetscOptionsTail();
760: return(0);
761: }
763: /*------------------------------------------------------------------------------------*/
765: static PetscErrorCode PCASMSetLocalSubdomains_ASM(PC pc,PetscInt n,IS is[],IS is_local[])
766: {
767: PC_ASM *osm = (PC_ASM*)pc->data;
769: PetscInt i;
772: if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Each process must have 1 or more blocks, n = %D",n);
773: if (pc->setupcalled && (n != osm->n_local_true || is)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetLocalSubdomains() should be called before calling PCSetUp().");
775: if (!pc->setupcalled) {
776: if (is) {
777: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is[i]);}
778: }
779: if (is_local) {
780: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is_local[i]);}
781: }
782: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
784: osm->n_local_true = n;
785: osm->is = NULL;
786: osm->is_local = NULL;
787: if (is) {
788: PetscMalloc1(n,&osm->is);
789: for (i=0; i<n; i++) osm->is[i] = is[i];
790: /* Flag indicating that the user has set overlapping subdomains so PCASM should not increase their size. */
791: osm->overlap = -1;
792: }
793: if (is_local) {
794: PetscMalloc1(n,&osm->is_local);
795: for (i=0; i<n; i++) osm->is_local[i] = is_local[i];
796: if (!is) {
797: PetscMalloc1(osm->n_local_true,&osm->is);
798: for (i=0; i<osm->n_local_true; i++) {
799: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
800: ISDuplicate(osm->is_local[i],&osm->is[i]);
801: ISCopy(osm->is_local[i],osm->is[i]);
802: } else {
803: PetscObjectReference((PetscObject)osm->is_local[i]);
804: osm->is[i] = osm->is_local[i];
805: }
806: }
807: }
808: }
809: }
810: return(0);
811: }
813: static PetscErrorCode PCASMSetTotalSubdomains_ASM(PC pc,PetscInt N,IS *is,IS *is_local)
814: {
815: PC_ASM *osm = (PC_ASM*)pc->data;
817: PetscMPIInt rank,size;
818: PetscInt n;
821: if (N < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of total blocks must be > 0, N = %D",N);
822: 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.");
824: /*
825: Split the subdomains equally among all processors
826: */
827: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
828: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
829: n = N/size + ((N % size) > rank);
830: 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);
831: if (pc->setupcalled && n != osm->n_local_true) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCASMSetTotalSubdomains() should be called before PCSetUp().");
832: if (!pc->setupcalled) {
833: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
835: osm->n_local_true = n;
836: osm->is = NULL;
837: osm->is_local = NULL;
838: }
839: return(0);
840: }
842: static PetscErrorCode PCASMSetOverlap_ASM(PC pc,PetscInt ovl)
843: {
844: PC_ASM *osm = (PC_ASM*)pc->data;
847: if (ovl < 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested");
848: if (pc->setupcalled && ovl != osm->overlap) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetOverlap() should be called before PCSetUp().");
849: if (!pc->setupcalled) osm->overlap = ovl;
850: return(0);
851: }
853: static PetscErrorCode PCASMSetType_ASM(PC pc,PCASMType type)
854: {
855: PC_ASM *osm = (PC_ASM*)pc->data;
858: osm->type = type;
859: osm->type_set = PETSC_TRUE;
860: return(0);
861: }
863: static PetscErrorCode PCASMGetType_ASM(PC pc,PCASMType *type)
864: {
865: PC_ASM *osm = (PC_ASM*)pc->data;
868: *type = osm->type;
869: return(0);
870: }
872: static PetscErrorCode PCASMSetLocalType_ASM(PC pc, PCCompositeType type)
873: {
874: PC_ASM *osm = (PC_ASM *) pc->data;
877: 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");
878: osm->loctype = type;
879: return(0);
880: }
882: static PetscErrorCode PCASMGetLocalType_ASM(PC pc, PCCompositeType *type)
883: {
884: PC_ASM *osm = (PC_ASM *) pc->data;
887: *type = osm->loctype;
888: return(0);
889: }
891: static PetscErrorCode PCASMSetSortIndices_ASM(PC pc,PetscBool doSort)
892: {
893: PC_ASM *osm = (PC_ASM*)pc->data;
896: osm->sort_indices = doSort;
897: return(0);
898: }
900: static PetscErrorCode PCASMGetSubKSP_ASM(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
901: {
902: PC_ASM *osm = (PC_ASM*)pc->data;
906: 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");
908: if (n_local) *n_local = osm->n_local_true;
909: if (first_local) {
910: MPI_Scan(&osm->n_local_true,first_local,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
911: *first_local -= osm->n_local_true;
912: }
913: if (ksp) {
914: /* Assume that local solves are now different; not necessarily
915: true though! This flag is used only for PCView_ASM() */
916: *ksp = osm->ksp;
917: osm->same_local_solves = PETSC_FALSE;
918: }
919: return(0);
920: }
922: static PetscErrorCode PCASMGetSubMatType_ASM(PC pc,MatType *sub_mat_type)
923: {
924: PC_ASM *osm = (PC_ASM*)pc->data;
929: *sub_mat_type = osm->sub_mat_type;
930: return(0);
931: }
933: static PetscErrorCode PCASMSetSubMatType_ASM(PC pc,MatType sub_mat_type)
934: {
935: PetscErrorCode ierr;
936: PC_ASM *osm = (PC_ASM*)pc->data;
940: PetscFree(osm->sub_mat_type);
941: PetscStrallocpy(sub_mat_type,(char**)&osm->sub_mat_type);
942: return(0);
943: }
945: /*@C
946: PCASMSetLocalSubdomains - Sets the local subdomains (for this processor only) for the additive Schwarz preconditioner.
948: Collective on pc
950: Input Parameters:
951: + pc - the preconditioner context
952: . n - the number of subdomains for this processor (default value = 1)
953: . is - the index set that defines the subdomains for this processor
954: (or NULL for PETSc to determine subdomains)
955: - is_local - the index sets that define the local part of the subdomains for this processor, not used unless PCASMType is PC_ASM_RESTRICT
956: (or NULL to not provide these)
958: Options Database Key:
959: To set the total number of subdomain blocks rather than specify the
960: index sets, use the option
961: . -pc_asm_local_blocks <blks> - Sets local blocks
963: Notes:
964: The IS numbering is in the parallel, global numbering of the vector for both is and is_local
966: By default the ASM preconditioner uses 1 block per processor.
968: Use PCASMSetTotalSubdomains() to set the subdomains for all processors.
970: If is_local is provided and PCASMType is PC_ASM_RESTRICT then the solution only over the is_local region is interpolated
971: back to form the global solution (this is the standard restricted additive Schwarz method)
973: If the is_local is provided and PCASMType is PC_ASM_INTERPOLATE or PC_ASM_NONE then an error is generated since there is
974: no code to handle that case.
976: Level: advanced
978: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
979: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), PCASMType, PCASMSetType()
980: @*/
981: PetscErrorCode PCASMSetLocalSubdomains(PC pc,PetscInt n,IS is[],IS is_local[])
982: {
987: PetscTryMethod(pc,"PCASMSetLocalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,is,is_local));
988: return(0);
989: }
991: /*@C
992: PCASMSetTotalSubdomains - Sets the subdomains for all processors for the
993: additive Schwarz preconditioner. Either all or no processors in the
994: PC communicator must call this routine, with the same index sets.
996: Collective on pc
998: Input Parameters:
999: + pc - the preconditioner context
1000: . N - the number of subdomains for all processors
1001: . is - the index sets that define the subdomains for all processors
1002: (or NULL to ask PETSc to determine the subdomains)
1003: - is_local - the index sets that define the local part of the subdomains for this processor
1004: (or NULL to not provide this information)
1006: Options Database Key:
1007: To set the total number of subdomain blocks rather than specify the
1008: index sets, use the option
1009: . -pc_asm_blocks <blks> - Sets total blocks
1011: Notes:
1012: Currently you cannot use this to set the actual subdomains with the argument is or is_local.
1014: By default the ASM preconditioner uses 1 block per processor.
1016: These index sets cannot be destroyed until after completion of the
1017: linear solves for which the ASM preconditioner is being used.
1019: Use PCASMSetLocalSubdomains() to set local subdomains.
1021: The IS numbering is in the parallel, global numbering of the vector for both is and is_local
1023: Level: advanced
1025: .seealso: PCASMSetLocalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1026: PCASMCreateSubdomains2D()
1027: @*/
1028: PetscErrorCode PCASMSetTotalSubdomains(PC pc,PetscInt N,IS is[],IS is_local[])
1029: {
1034: PetscTryMethod(pc,"PCASMSetTotalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,N,is,is_local));
1035: return(0);
1036: }
1038: /*@
1039: PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
1040: additive Schwarz preconditioner. Either all or no processors in the
1041: PC communicator must call this routine.
1043: Logically Collective on pc
1045: Input Parameters:
1046: + pc - the preconditioner context
1047: - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1)
1049: Options Database Key:
1050: . -pc_asm_overlap <ovl> - Sets overlap
1052: Notes:
1053: By default the ASM preconditioner uses 1 block per processor. To use
1054: multiple blocks per perocessor, see PCASMSetTotalSubdomains() and
1055: PCASMSetLocalSubdomains() (and the option -pc_asm_blocks <blks>).
1057: The overlap defaults to 1, so if one desires that no additional
1058: overlap be computed beyond what may have been set with a call to
1059: PCASMSetTotalSubdomains() or PCASMSetLocalSubdomains(), then ovl
1060: must be set to be 0. In particular, if one does not explicitly set
1061: the subdomains an application code, then all overlap would be computed
1062: internally by PETSc, and using an overlap of 0 would result in an ASM
1063: variant that is equivalent to the block Jacobi preconditioner.
1065: The default algorithm used by PETSc to increase overlap is fast, but not scalable,
1066: use the option -mat_increase_overlap_scalable when the problem and number of processes is large.
1068: Note that one can define initial index sets with any overlap via
1069: PCASMSetLocalSubdomains(); the routine
1070: PCASMSetOverlap() merely allows PETSc to extend that overlap further
1071: if desired.
1073: Level: intermediate
1075: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1076: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), MatIncreaseOverlap()
1077: @*/
1078: PetscErrorCode PCASMSetOverlap(PC pc,PetscInt ovl)
1079: {
1085: PetscTryMethod(pc,"PCASMSetOverlap_C",(PC,PetscInt),(pc,ovl));
1086: return(0);
1087: }
1089: /*@
1090: PCASMSetType - Sets the type of restriction and interpolation used
1091: for local problems in the additive Schwarz method.
1093: Logically Collective on pc
1095: Input Parameters:
1096: + pc - the preconditioner context
1097: - type - variant of ASM, one of
1098: .vb
1099: PC_ASM_BASIC - full interpolation and restriction
1100: PC_ASM_RESTRICT - full restriction, local processor interpolation
1101: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1102: PC_ASM_NONE - local processor restriction and interpolation
1103: .ve
1105: Options Database Key:
1106: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
1108: Notes:
1109: if the is_local arguments are passed to PCASMSetLocalSubdomains() then they are used when PC_ASM_RESTRICT has been selected
1110: to limit the local processor interpolation
1112: Level: intermediate
1114: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1115: PCASMCreateSubdomains2D(), PCASMType, PCASMSetLocalType(), PCASMGetLocalType()
1116: @*/
1117: PetscErrorCode PCASMSetType(PC pc,PCASMType type)
1118: {
1124: PetscTryMethod(pc,"PCASMSetType_C",(PC,PCASMType),(pc,type));
1125: return(0);
1126: }
1128: /*@
1129: PCASMGetType - Gets the type of restriction and interpolation used
1130: for local problems in the additive Schwarz method.
1132: Logically Collective on pc
1134: Input Parameter:
1135: . pc - the preconditioner context
1137: Output Parameter:
1138: . type - variant of ASM, one of
1140: .vb
1141: PC_ASM_BASIC - full interpolation and restriction
1142: PC_ASM_RESTRICT - full restriction, local processor interpolation
1143: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1144: PC_ASM_NONE - local processor restriction and interpolation
1145: .ve
1147: Options Database Key:
1148: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
1150: Level: intermediate
1152: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1153: PCASMCreateSubdomains2D(), PCASMType, PCASMSetType(), PCASMSetLocalType(), PCASMGetLocalType()
1154: @*/
1155: PetscErrorCode PCASMGetType(PC pc,PCASMType *type)
1156: {
1161: PetscUseMethod(pc,"PCASMGetType_C",(PC,PCASMType*),(pc,type));
1162: return(0);
1163: }
1165: /*@
1166: PCASMSetLocalType - Sets the type of composition used for local problems in the additive Schwarz method.
1168: Logically Collective on pc
1170: Input Parameters:
1171: + pc - the preconditioner context
1172: - type - type of composition, one of
1173: .vb
1174: PC_COMPOSITE_ADDITIVE - local additive combination
1175: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1176: .ve
1178: Options Database Key:
1179: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1181: Level: intermediate
1183: .seealso: PCASMSetType(), PCASMGetType(), PCASMGetLocalType(), PCASM, PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1184: @*/
1185: PetscErrorCode PCASMSetLocalType(PC pc, PCCompositeType type)
1186: {
1192: PetscTryMethod(pc, "PCASMSetLocalType_C", (PC, PCCompositeType), (pc, type));
1193: return(0);
1194: }
1196: /*@
1197: PCASMGetLocalType - Gets the type of composition used for local problems in the additive Schwarz method.
1199: Logically Collective on pc
1201: Input Parameter:
1202: . pc - the preconditioner context
1204: Output Parameter:
1205: . type - type of composition, one of
1206: .vb
1207: PC_COMPOSITE_ADDITIVE - local additive combination
1208: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1209: .ve
1211: Options Database Key:
1212: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1214: Level: intermediate
1216: .seealso: PCASMSetType(), PCASMGetType(), PCASMSetLocalType(), PCASMCreate(), PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1217: @*/
1218: PetscErrorCode PCASMGetLocalType(PC pc, PCCompositeType *type)
1219: {
1225: PetscUseMethod(pc, "PCASMGetLocalType_C", (PC, PCCompositeType *), (pc, type));
1226: return(0);
1227: }
1229: /*@
1230: PCASMSetSortIndices - Determines whether subdomain indices are sorted.
1232: Logically Collective on pc
1234: Input Parameters:
1235: + pc - the preconditioner context
1236: - doSort - sort the subdomain indices
1238: Level: intermediate
1240: .seealso: PCASMSetLocalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1241: PCASMCreateSubdomains2D()
1242: @*/
1243: PetscErrorCode PCASMSetSortIndices(PC pc,PetscBool doSort)
1244: {
1250: PetscTryMethod(pc,"PCASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));
1251: return(0);
1252: }
1254: /*@C
1255: PCASMGetSubKSP - Gets the local KSP contexts for all blocks on
1256: this processor.
1258: Collective on pc iff first_local is requested
1260: Input Parameter:
1261: . pc - the preconditioner context
1263: Output Parameters:
1264: + n_local - the number of blocks on this processor or NULL
1265: . first_local - the global number of the first block on this processor or NULL,
1266: all processors must request or all must pass NULL
1267: - ksp - the array of KSP contexts
1269: Note:
1270: After PCASMGetSubKSP() the array of KSPes is not to be freed.
1272: You must call KSPSetUp() before calling PCASMGetSubKSP().
1274: Fortran note:
1275: 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.
1277: Level: advanced
1279: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap(),
1280: PCASMCreateSubdomains2D(),
1281: @*/
1282: PetscErrorCode PCASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
1283: {
1288: PetscUseMethod(pc,"PCASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
1289: return(0);
1290: }
1292: /* -------------------------------------------------------------------------------------*/
1293: /*MC
1294: PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1295: its own KSP object.
1297: Options Database Keys:
1298: + -pc_asm_blocks <blks> - Sets total blocks
1299: . -pc_asm_overlap <ovl> - Sets overlap
1300: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type, default is restrict
1301: - -pc_asm_local_type [additive, multiplicative] - Sets ASM type, default is additive
1303: IMPORTANT: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you
1304: will get a different convergence rate due to the default option of -pc_asm_type restrict. Use
1305: -pc_asm_type basic to use the standard ASM.
1307: Notes:
1308: Each processor can have one or more blocks, but a block cannot be shared by more
1309: than one processor. Use PCGASM for subdomains shared by multiple processes. Defaults to one block per processor.
1311: To set options on the solvers for each block append -sub_ to all the KSP, and PC
1312: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
1314: To set the options on the solvers separate for each block call PCASMGetSubKSP()
1315: and set the options directly on the resulting KSP object (you can access its PC
1316: with KSPGetPC())
1318: Level: beginner
1320: References:
1321: + 1. - M Dryja, OB Widlund, An additive variant of the Schwarz alternating method for the case of many subregions
1322: Courant Institute, New York University Technical report
1323: - 2. - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations,
1324: Cambridge University Press.
1326: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
1327: PCBJACOBI, PCASMGetSubKSP(), PCASMSetLocalSubdomains(), PCASMType, PCASMGetType(), PCASMSetLocalType(), PCASMGetLocalType()
1328: PCASMSetTotalSubdomains(), PCSetModifySubMatrices(), PCASMSetOverlap(), PCASMSetType(), PCCompositeType
1330: M*/
1332: PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1333: {
1335: PC_ASM *osm;
1338: PetscNewLog(pc,&osm);
1340: osm->n = PETSC_DECIDE;
1341: osm->n_local = 0;
1342: osm->n_local_true = PETSC_DECIDE;
1343: osm->overlap = 1;
1344: osm->ksp = NULL;
1345: osm->restriction = NULL;
1346: osm->lprolongation = NULL;
1347: osm->lrestriction = NULL;
1348: osm->x = NULL;
1349: osm->y = NULL;
1350: osm->is = NULL;
1351: osm->is_local = NULL;
1352: osm->mat = NULL;
1353: osm->pmat = NULL;
1354: osm->type = PC_ASM_RESTRICT;
1355: osm->loctype = PC_COMPOSITE_ADDITIVE;
1356: osm->same_local_solves = PETSC_TRUE;
1357: osm->sort_indices = PETSC_TRUE;
1358: osm->dm_subdomains = PETSC_FALSE;
1359: osm->sub_mat_type = NULL;
1361: pc->data = (void*)osm;
1362: pc->ops->apply = PCApply_ASM;
1363: pc->ops->matapply = PCMatApply_ASM;
1364: pc->ops->applytranspose = PCApplyTranspose_ASM;
1365: pc->ops->setup = PCSetUp_ASM;
1366: pc->ops->reset = PCReset_ASM;
1367: pc->ops->destroy = PCDestroy_ASM;
1368: pc->ops->setfromoptions = PCSetFromOptions_ASM;
1369: pc->ops->setuponblocks = PCSetUpOnBlocks_ASM;
1370: pc->ops->view = PCView_ASM;
1371: pc->ops->applyrichardson = NULL;
1373: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",PCASMSetLocalSubdomains_ASM);
1374: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",PCASMSetTotalSubdomains_ASM);
1375: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetOverlap_C",PCASMSetOverlap_ASM);
1376: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetType_C",PCASMSetType_ASM);
1377: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetType_C",PCASMGetType_ASM);
1378: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalType_C",PCASMSetLocalType_ASM);
1379: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetLocalType_C",PCASMGetLocalType_ASM);
1380: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSortIndices_C",PCASMSetSortIndices_ASM);
1381: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubKSP_C",PCASMGetSubKSP_ASM);
1382: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubMatType_C",PCASMGetSubMatType_ASM);
1383: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSubMatType_C",PCASMSetSubMatType_ASM);
1384: return(0);
1385: }
1387: /*@C
1388: PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1389: preconditioner for a any problem on a general grid.
1391: Collective
1393: Input Parameters:
1394: + A - The global matrix operator
1395: - n - the number of local blocks
1397: Output Parameters:
1398: . outis - the array of index sets defining the subdomains
1400: Level: advanced
1402: Note: this generates nonoverlapping subdomains; the PCASM will generate the overlap
1403: from these if you use PCASMSetLocalSubdomains()
1405: In the Fortran version you must provide the array outis[] already allocated of length n.
1407: .seealso: PCASMSetLocalSubdomains(), PCASMDestroySubdomains()
1408: @*/
1409: PetscErrorCode PCASMCreateSubdomains(Mat A, PetscInt n, IS* outis[])
1410: {
1411: MatPartitioning mpart;
1412: const char *prefix;
1413: PetscInt i,j,rstart,rend,bs;
1414: PetscBool hasop, isbaij = PETSC_FALSE,foundpart = PETSC_FALSE;
1415: Mat Ad = NULL, adj;
1416: IS ispart,isnumb,*is;
1417: PetscErrorCode ierr;
1422: if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of local blocks must be > 0, n = %D",n);
1424: /* Get prefix, row distribution, and block size */
1425: MatGetOptionsPrefix(A,&prefix);
1426: MatGetOwnershipRange(A,&rstart,&rend);
1427: MatGetBlockSize(A,&bs);
1428: 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);
1430: /* Get diagonal block from matrix if possible */
1431: MatHasOperation(A,MATOP_GET_DIAGONAL_BLOCK,&hasop);
1432: if (hasop) {
1433: MatGetDiagonalBlock(A,&Ad);
1434: }
1435: if (Ad) {
1436: PetscObjectBaseTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij);
1437: if (!isbaij) {PetscObjectBaseTypeCompare((PetscObject)Ad,MATSEQSBAIJ,&isbaij);}
1438: }
1439: if (Ad && n > 1) {
1440: PetscBool match,done;
1441: /* Try to setup a good matrix partitioning if available */
1442: MatPartitioningCreate(PETSC_COMM_SELF,&mpart);
1443: PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);
1444: MatPartitioningSetFromOptions(mpart);
1445: PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGCURRENT,&match);
1446: if (!match) {
1447: PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGSQUARE,&match);
1448: }
1449: if (!match) { /* assume a "good" partitioner is available */
1450: PetscInt na;
1451: const PetscInt *ia,*ja;
1452: MatGetRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1453: if (done) {
1454: /* Build adjacency matrix by hand. Unfortunately a call to
1455: MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1456: remove the block-aij structure and we cannot expect
1457: MatPartitioning to split vertices as we need */
1458: PetscInt i,j,len,nnz,cnt,*iia=NULL,*jja=NULL;
1459: const PetscInt *row;
1460: nnz = 0;
1461: for (i=0; i<na; i++) { /* count number of nonzeros */
1462: len = ia[i+1] - ia[i];
1463: row = ja + ia[i];
1464: for (j=0; j<len; j++) {
1465: if (row[j] == i) { /* don't count diagonal */
1466: len--; break;
1467: }
1468: }
1469: nnz += len;
1470: }
1471: PetscMalloc1(na+1,&iia);
1472: PetscMalloc1(nnz,&jja);
1473: nnz = 0;
1474: iia[0] = 0;
1475: for (i=0; i<na; i++) { /* fill adjacency */
1476: cnt = 0;
1477: len = ia[i+1] - ia[i];
1478: row = ja + ia[i];
1479: for (j=0; j<len; j++) {
1480: if (row[j] != i) { /* if not diagonal */
1481: jja[nnz+cnt++] = row[j];
1482: }
1483: }
1484: nnz += cnt;
1485: iia[i+1] = nnz;
1486: }
1487: /* Partitioning of the adjacency matrix */
1488: MatCreateMPIAdj(PETSC_COMM_SELF,na,na,iia,jja,NULL,&adj);
1489: MatPartitioningSetAdjacency(mpart,adj);
1490: MatPartitioningSetNParts(mpart,n);
1491: MatPartitioningApply(mpart,&ispart);
1492: ISPartitioningToNumbering(ispart,&isnumb);
1493: MatDestroy(&adj);
1494: foundpart = PETSC_TRUE;
1495: }
1496: MatRestoreRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1497: }
1498: MatPartitioningDestroy(&mpart);
1499: }
1501: PetscMalloc1(n,&is);
1502: *outis = is;
1504: if (!foundpart) {
1506: /* Partitioning by contiguous chunks of rows */
1508: PetscInt mbs = (rend-rstart)/bs;
1509: PetscInt start = rstart;
1510: for (i=0; i<n; i++) {
1511: PetscInt count = (mbs/n + ((mbs % n) > i)) * bs;
1512: ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);
1513: start += count;
1514: }
1516: } else {
1518: /* Partitioning by adjacency of diagonal block */
1520: const PetscInt *numbering;
1521: PetscInt *count,nidx,*indices,*newidx,start=0;
1522: /* Get node count in each partition */
1523: PetscMalloc1(n,&count);
1524: ISPartitioningCount(ispart,n,count);
1525: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1526: for (i=0; i<n; i++) count[i] *= bs;
1527: }
1528: /* Build indices from node numbering */
1529: ISGetLocalSize(isnumb,&nidx);
1530: PetscMalloc1(nidx,&indices);
1531: for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */
1532: ISGetIndices(isnumb,&numbering);
1533: PetscSortIntWithPermutation(nidx,numbering,indices);
1534: ISRestoreIndices(isnumb,&numbering);
1535: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1536: PetscMalloc1(nidx*bs,&newidx);
1537: for (i=0; i<nidx; i++) {
1538: for (j=0; j<bs; j++) newidx[i*bs+j] = indices[i]*bs + j;
1539: }
1540: PetscFree(indices);
1541: nidx *= bs;
1542: indices = newidx;
1543: }
1544: /* Shift to get global indices */
1545: for (i=0; i<nidx; i++) indices[i] += rstart;
1547: /* Build the index sets for each block */
1548: for (i=0; i<n; i++) {
1549: ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);
1550: ISSort(is[i]);
1551: start += count[i];
1552: }
1554: PetscFree(count);
1555: PetscFree(indices);
1556: ISDestroy(&isnumb);
1557: ISDestroy(&ispart);
1559: }
1560: return(0);
1561: }
1563: /*@C
1564: PCASMDestroySubdomains - Destroys the index sets created with
1565: PCASMCreateSubdomains(). Should be called after setting subdomains
1566: with PCASMSetLocalSubdomains().
1568: Collective
1570: Input Parameters:
1571: + n - the number of index sets
1572: . is - the array of index sets
1573: - is_local - the array of local index sets, can be NULL
1575: Level: advanced
1577: .seealso: PCASMCreateSubdomains(), PCASMSetLocalSubdomains()
1578: @*/
1579: PetscErrorCode PCASMDestroySubdomains(PetscInt n, IS is[], IS is_local[])
1580: {
1581: PetscInt i;
1585: if (n <= 0) return(0);
1586: if (is) {
1588: for (i=0; i<n; i++) { ISDestroy(&is[i]); }
1589: PetscFree(is);
1590: }
1591: if (is_local) {
1593: for (i=0; i<n; i++) { ISDestroy(&is_local[i]); }
1594: PetscFree(is_local);
1595: }
1596: return(0);
1597: }
1599: /*@
1600: PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1601: preconditioner for a two-dimensional problem on a regular grid.
1603: Not Collective
1605: Input Parameters:
1606: + m, n - the number of mesh points in the x and y directions
1607: . M, N - the number of subdomains in the x and y directions
1608: . dof - degrees of freedom per node
1609: - overlap - overlap in mesh lines
1611: Output Parameters:
1612: + Nsub - the number of subdomains created
1613: . is - array of index sets defining overlapping (if overlap > 0) subdomains
1614: - is_local - array of index sets defining non-overlapping subdomains
1616: Note:
1617: Presently PCAMSCreateSubdomains2d() is valid only for sequential
1618: preconditioners. More general related routines are
1619: PCASMSetTotalSubdomains() and PCASMSetLocalSubdomains().
1621: Level: advanced
1623: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1624: PCASMSetOverlap()
1625: @*/
1626: PetscErrorCode PCASMCreateSubdomains2D(PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt dof,PetscInt overlap,PetscInt *Nsub,IS **is,IS **is_local)
1627: {
1628: PetscInt i,j,height,width,ystart,xstart,yleft,yright,xleft,xright,loc_outer;
1630: PetscInt nidx,*idx,loc,ii,jj,count;
1633: if (dof != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP," ");
1635: *Nsub = N*M;
1636: PetscMalloc1(*Nsub,is);
1637: PetscMalloc1(*Nsub,is_local);
1638: ystart = 0;
1639: loc_outer = 0;
1640: for (i=0; i<N; i++) {
1641: height = n/N + ((n % N) > i); /* height of subdomain */
1642: if (height < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many N subdomains for mesh dimension n");
1643: yleft = ystart - overlap; if (yleft < 0) yleft = 0;
1644: yright = ystart + height + overlap; if (yright > n) yright = n;
1645: xstart = 0;
1646: for (j=0; j<M; j++) {
1647: width = m/M + ((m % M) > j); /* width of subdomain */
1648: if (width < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many M subdomains for mesh dimension m");
1649: xleft = xstart - overlap; if (xleft < 0) xleft = 0;
1650: xright = xstart + width + overlap; if (xright > m) xright = m;
1651: nidx = (xright - xleft)*(yright - yleft);
1652: PetscMalloc1(nidx,&idx);
1653: loc = 0;
1654: for (ii=yleft; ii<yright; ii++) {
1655: count = m*ii + xleft;
1656: for (jj=xleft; jj<xright; jj++) idx[loc++] = count++;
1657: }
1658: ISCreateGeneral(PETSC_COMM_SELF,nidx,idx,PETSC_COPY_VALUES,(*is)+loc_outer);
1659: if (overlap == 0) {
1660: PetscObjectReference((PetscObject)(*is)[loc_outer]);
1662: (*is_local)[loc_outer] = (*is)[loc_outer];
1663: } else {
1664: for (loc=0,ii=ystart; ii<ystart+height; ii++) {
1665: for (jj=xstart; jj<xstart+width; jj++) {
1666: idx[loc++] = m*ii + jj;
1667: }
1668: }
1669: ISCreateGeneral(PETSC_COMM_SELF,loc,idx,PETSC_COPY_VALUES,*is_local+loc_outer);
1670: }
1671: PetscFree(idx);
1672: xstart += width;
1673: loc_outer++;
1674: }
1675: ystart += height;
1676: }
1677: for (i=0; i<*Nsub; i++) { ISSort((*is)[i]); }
1678: return(0);
1679: }
1681: /*@C
1682: PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
1683: only) for the additive Schwarz preconditioner.
1685: Not Collective
1687: Input Parameter:
1688: . pc - the preconditioner context
1690: Output Parameters:
1691: + n - if requested, the number of subdomains for this processor (default value = 1)
1692: . is - if requested, the index sets that define the subdomains for this processor
1693: - is_local - if requested, the index sets that define the local part of the subdomains for this processor (can be NULL)
1696: Notes:
1697: The IS numbering is in the parallel, global numbering of the vector.
1699: Level: advanced
1701: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1702: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubmatrices()
1703: @*/
1704: PetscErrorCode PCASMGetLocalSubdomains(PC pc,PetscInt *n,IS *is[],IS *is_local[])
1705: {
1706: PC_ASM *osm = (PC_ASM*)pc->data;
1708: PetscBool match;
1715: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1716: if (!match) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PC is not a PCASM");
1717: if (n) *n = osm->n_local_true;
1718: if (is) *is = osm->is;
1719: if (is_local) *is_local = osm->is_local;
1720: return(0);
1721: }
1723: /*@C
1724: PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1725: only) for the additive Schwarz preconditioner.
1727: Not Collective
1729: Input Parameter:
1730: . pc - the preconditioner context
1732: Output Parameters:
1733: + n - if requested, the number of matrices for this processor (default value = 1)
1734: - mat - if requested, the matrices
1736: Level: advanced
1738: Notes:
1739: Call after PCSetUp() (or KSPSetUp()) but before PCApply() and before PCSetUpOnBlocks())
1741: Usually one would use PCSetModifySubMatrices() to change the submatrices in building the preconditioner.
1743: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1744: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains(), PCSetModifySubMatrices()
1745: @*/
1746: PetscErrorCode PCASMGetLocalSubmatrices(PC pc,PetscInt *n,Mat *mat[])
1747: {
1748: PC_ASM *osm;
1750: PetscBool match;
1756: if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUp() or PCSetUp().");
1757: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1758: if (!match) {
1759: if (n) *n = 0;
1760: if (mat) *mat = NULL;
1761: } else {
1762: osm = (PC_ASM*)pc->data;
1763: if (n) *n = osm->n_local_true;
1764: if (mat) *mat = osm->pmat;
1765: }
1766: return(0);
1767: }
1769: /*@
1770: PCASMSetDMSubdomains - Indicates whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1772: Logically Collective
1774: Input Parameter:
1775: + pc - the preconditioner
1776: - flg - boolean indicating whether to use subdomains defined by the DM
1778: Options Database Key:
1779: . -pc_asm_dm_subdomains
1781: Level: intermediate
1783: Notes:
1784: PCASMSetTotalSubdomains() and PCASMSetOverlap() take precedence over PCASMSetDMSubdomains(),
1785: so setting either of the first two effectively turns the latter off.
1787: .seealso: PCASMGetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1788: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1789: @*/
1790: PetscErrorCode PCASMSetDMSubdomains(PC pc,PetscBool flg)
1791: {
1792: PC_ASM *osm = (PC_ASM*)pc->data;
1794: PetscBool match;
1799: if (pc->setupcalled) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for a setup PC.");
1800: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1801: if (match) {
1802: osm->dm_subdomains = flg;
1803: }
1804: return(0);
1805: }
1807: /*@
1808: PCASMGetDMSubdomains - Returns flag indicating whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1809: Not Collective
1811: Input Parameter:
1812: . pc - the preconditioner
1814: Output Parameter:
1815: . flg - boolean indicating whether to use subdomains defined by the DM
1817: Level: intermediate
1819: .seealso: PCASMSetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1820: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1821: @*/
1822: PetscErrorCode PCASMGetDMSubdomains(PC pc,PetscBool* flg)
1823: {
1824: PC_ASM *osm = (PC_ASM*)pc->data;
1826: PetscBool match;
1831: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1832: if (match) *flg = osm->dm_subdomains;
1833: else *flg = PETSC_FALSE;
1834: return(0);
1835: }
1837: /*@
1838: PCASMGetSubMatType - Gets the matrix type used for ASM subsolves, as a string.
1840: Not Collective
1842: Input Parameter:
1843: . pc - the PC
1845: Output Parameter:
1846: . -pc_asm_sub_mat_type - name of matrix type
1848: Level: advanced
1850: .seealso: PCASMSetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1851: @*/
1852: PetscErrorCode PCASMGetSubMatType(PC pc,MatType *sub_mat_type)
1853: {
1858: PetscTryMethod(pc,"PCASMGetSubMatType_C",(PC,MatType*),(pc,sub_mat_type));
1859: return(0);
1860: }
1862: /*@
1863: PCASMSetSubMatType - Set the type of matrix used for ASM subsolves
1865: Collective on Mat
1867: Input Parameters:
1868: + pc - the PC object
1869: - sub_mat_type - matrix type
1871: Options Database Key:
1872: . -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.
1874: Notes:
1875: See "${PETSC_DIR}/include/petscmat.h" for available types
1877: Level: advanced
1879: .seealso: PCASMGetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1880: @*/
1881: PetscErrorCode PCASMSetSubMatType(PC pc,MatType sub_mat_type)
1882: {
1887: PetscTryMethod(pc,"PCASMSetSubMatType_C",(PC,MatType),(pc,sub_mat_type));
1888: return(0);
1889: }