Actual source code: asm.c
1: /*
2: This file defines an additive Schwarz preconditioner for any Mat implementation.
4: Note that each processor may have any number of subdomains. But in order to
5: deal easily with the VecScatter(), we treat each processor as if it has the
6: same number of subdomains.
8: n - total number of true subdomains on all processors
9: n_local_true - actual number of subdomains on this processor
10: n_local = maximum over all processors of n_local_true
11: */
13: #include <petsc/private/pcasmimpl.h>
15: static PetscErrorCode PCView_ASM(PC pc,PetscViewer viewer)
16: {
17: PC_ASM *osm = (PC_ASM*)pc->data;
18: PetscErrorCode ierr;
19: PetscMPIInt rank;
20: PetscInt i,bsz;
21: PetscBool iascii,isstring;
22: PetscViewer sviewer;
23: PetscViewerFormat format;
24: const char *prefix;
27: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
28: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
29: if (iascii) {
30: char overlaps[256] = "user-defined overlap",blocks[256] = "total subdomain blocks not yet set";
31: if (osm->overlap >= 0) {PetscSNPrintf(overlaps,sizeof(overlaps),"amount of overlap = %D",osm->overlap);}
32: if (osm->n > 0) {PetscSNPrintf(blocks,sizeof(blocks),"total subdomain blocks = %D",osm->n);}
33: PetscViewerASCIIPrintf(viewer," %s, %s\n",blocks,overlaps);
34: PetscViewerASCIIPrintf(viewer," restriction/interpolation type - %s\n",PCASMTypes[osm->type]);
35: if (osm->dm_subdomains) {PetscViewerASCIIPrintf(viewer," Additive Schwarz: using DM to define subdomains\n");}
36: if (osm->loctype != PC_COMPOSITE_ADDITIVE) {PetscViewerASCIIPrintf(viewer," Additive Schwarz: local solve composition type - %s\n",PCCompositeTypes[osm->loctype]);}
37: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
38: PetscViewerGetFormat(viewer,&format);
39: if (format != PETSC_VIEWER_ASCII_INFO_DETAIL) {
40: if (osm->ksp) {
41: PetscViewerASCIIPrintf(viewer," Local solver information for first block is in the following KSP and PC objects on rank 0:\n");
42: PCGetOptionsPrefix(pc,&prefix);
43: PetscViewerASCIIPrintf(viewer," Use -%sksp_view ::ascii_info_detail to display information for all blocks\n",prefix?prefix:"");
44: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
45: if (!rank) {
46: PetscViewerASCIIPushTab(viewer);
47: KSPView(osm->ksp[0],sviewer);
48: PetscViewerASCIIPopTab(viewer);
49: }
50: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
51: }
52: } else {
53: PetscViewerASCIIPushSynchronized(viewer);
54: PetscViewerASCIISynchronizedPrintf(viewer," [%d] number of local blocks = %D\n",(int)rank,osm->n_local_true);
55: PetscViewerFlush(viewer);
56: PetscViewerASCIIPrintf(viewer," Local solver information for each block is in the following KSP and PC objects:\n");
57: PetscViewerASCIIPushTab(viewer);
58: PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
59: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
60: for (i=0; i<osm->n_local_true; i++) {
61: ISGetLocalSize(osm->is[i],&bsz);
62: PetscViewerASCIISynchronizedPrintf(sviewer,"[%d] local block number %D, size = %D\n",(int)rank,i,bsz);
63: KSPView(osm->ksp[i],sviewer);
64: PetscViewerASCIISynchronizedPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
65: }
66: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
67: PetscViewerASCIIPopTab(viewer);
68: PetscViewerFlush(viewer);
69: PetscViewerASCIIPopSynchronized(viewer);
70: }
71: } else if (isstring) {
72: PetscViewerStringSPrintf(viewer," blocks=%D, overlap=%D, type=%s",osm->n,osm->overlap,PCASMTypes[osm->type]);
73: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
74: if (osm->ksp) {KSPView(osm->ksp[0],sviewer);}
75: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
76: }
77: return(0);
78: }
80: static PetscErrorCode PCASMPrintSubdomains(PC pc)
81: {
82: PC_ASM *osm = (PC_ASM*)pc->data;
83: const char *prefix;
84: char fname[PETSC_MAX_PATH_LEN+1];
85: PetscViewer viewer, sviewer;
86: char *s;
87: PetscInt i,j,nidx;
88: const PetscInt *idx;
89: PetscMPIInt rank, size;
93: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
94: MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
95: PCGetOptionsPrefix(pc,&prefix);
96: PetscOptionsGetString(NULL,prefix,"-pc_asm_print_subdomains",fname,sizeof(fname),NULL);
97: if (fname[0] == 0) { PetscStrcpy(fname,"stdout"); };
98: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),fname,&viewer);
99: for (i=0; i<osm->n_local; i++) {
100: if (i < osm->n_local_true) {
101: ISGetLocalSize(osm->is[i],&nidx);
102: ISGetIndices(osm->is[i],&idx);
103: /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
104: #define len 16*(nidx+1)+512
105: PetscMalloc1(len,&s);
106: PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer);
107: #undef len
108: PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D with overlap:\n", rank, size, i);
109: for (j=0; j<nidx; j++) {
110: PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
111: }
112: ISRestoreIndices(osm->is[i],&idx);
113: PetscViewerStringSPrintf(sviewer,"\n");
114: PetscViewerDestroy(&sviewer);
115: PetscViewerASCIIPushSynchronized(viewer);
116: PetscViewerASCIISynchronizedPrintf(viewer, s);
117: PetscViewerFlush(viewer);
118: PetscViewerASCIIPopSynchronized(viewer);
119: PetscFree(s);
120: if (osm->is_local) {
121: /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
122: #define len 16*(nidx+1)+512
123: PetscMalloc1(len, &s);
124: PetscViewerStringOpen(PETSC_COMM_SELF, s, len, &sviewer);
125: #undef len
126: PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D without overlap:\n", rank, size, i);
127: ISGetLocalSize(osm->is_local[i],&nidx);
128: ISGetIndices(osm->is_local[i],&idx);
129: for (j=0; j<nidx; j++) {
130: PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
131: }
132: ISRestoreIndices(osm->is_local[i],&idx);
133: PetscViewerStringSPrintf(sviewer,"\n");
134: PetscViewerDestroy(&sviewer);
135: PetscViewerASCIIPushSynchronized(viewer);
136: PetscViewerASCIISynchronizedPrintf(viewer, s);
137: PetscViewerFlush(viewer);
138: PetscViewerASCIIPopSynchronized(viewer);
139: PetscFree(s);
140: }
141: } else {
142: /* Participate in collective viewer calls. */
143: PetscViewerASCIIPushSynchronized(viewer);
144: PetscViewerFlush(viewer);
145: PetscViewerASCIIPopSynchronized(viewer);
146: /* Assume either all ranks have is_local or none do. */
147: if (osm->is_local) {
148: PetscViewerASCIIPushSynchronized(viewer);
149: PetscViewerFlush(viewer);
150: PetscViewerASCIIPopSynchronized(viewer);
151: }
152: }
153: }
154: PetscViewerFlush(viewer);
155: PetscViewerDestroy(&viewer);
156: return(0);
157: }
159: static PetscErrorCode PCSetUp_ASM(PC pc)
160: {
161: PC_ASM *osm = (PC_ASM*)pc->data;
163: PetscBool flg;
164: PetscInt i,m,m_local;
165: MatReuse scall = MAT_REUSE_MATRIX;
166: IS isl;
167: KSP ksp;
168: PC subpc;
169: const char *prefix,*pprefix;
170: Vec vec;
171: DM *domain_dm = NULL;
174: if (!pc->setupcalled) {
175: PetscInt m;
177: /* Note: if subdomains have been set either via PCASMSetTotalSubdomains() or via PCASMSetLocalSubdomains(), osm->n_local_true will not be PETSC_DECIDE */
178: if (osm->n_local_true == PETSC_DECIDE) {
179: /* no subdomains given */
180: /* try pc->dm first, if allowed */
181: if (osm->dm_subdomains && pc->dm) {
182: PetscInt num_domains, d;
183: char **domain_names;
184: IS *inner_domain_is, *outer_domain_is;
185: DMCreateDomainDecomposition(pc->dm, &num_domains, &domain_names, &inner_domain_is, &outer_domain_is, &domain_dm);
186: osm->overlap = -1; /* We do not want to increase the overlap of the IS.
187: A future improvement of this code might allow one to use
188: DM-defined subdomains and also increase the overlap,
189: but that is not currently supported */
190: if (num_domains) {
191: PCASMSetLocalSubdomains(pc, num_domains, outer_domain_is, inner_domain_is);
192: }
193: for (d = 0; d < num_domains; ++d) {
194: if (domain_names) {PetscFree(domain_names[d]);}
195: if (inner_domain_is) {ISDestroy(&inner_domain_is[d]);}
196: if (outer_domain_is) {ISDestroy(&outer_domain_is[d]);}
197: }
198: PetscFree(domain_names);
199: PetscFree(inner_domain_is);
200: PetscFree(outer_domain_is);
201: }
202: if (osm->n_local_true == PETSC_DECIDE) {
203: /* still no subdomains; use one subdomain per processor */
204: osm->n_local_true = 1;
205: }
206: }
207: { /* determine the global and max number of subdomains */
208: struct {PetscInt max,sum;} inwork,outwork;
209: PetscMPIInt size;
211: inwork.max = osm->n_local_true;
212: inwork.sum = osm->n_local_true;
213: MPIU_Allreduce(&inwork,&outwork,1,MPIU_2INT,MPIU_MAXSUM_OP,PetscObjectComm((PetscObject)pc));
214: osm->n_local = outwork.max;
215: osm->n = outwork.sum;
217: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
218: if (outwork.max == 1 && outwork.sum == size) {
219: /* osm->n_local_true = 1 on all processes, set this option may enable use of optimized MatCreateSubMatrices() implementation */
220: MatSetOption(pc->pmat,MAT_SUBMAT_SINGLEIS,PETSC_TRUE);
221: }
222: }
223: if (!osm->is) { /* create the index sets */
224: PCASMCreateSubdomains(pc->pmat,osm->n_local_true,&osm->is);
225: }
226: if (osm->n_local_true > 1 && !osm->is_local) {
227: PetscMalloc1(osm->n_local_true,&osm->is_local);
228: for (i=0; i<osm->n_local_true; i++) {
229: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
230: ISDuplicate(osm->is[i],&osm->is_local[i]);
231: ISCopy(osm->is[i],osm->is_local[i]);
232: } else {
233: PetscObjectReference((PetscObject)osm->is[i]);
234: osm->is_local[i] = osm->is[i];
235: }
236: }
237: }
238: PCGetOptionsPrefix(pc,&prefix);
239: flg = PETSC_FALSE;
240: PetscOptionsGetBool(NULL,prefix,"-pc_asm_print_subdomains",&flg,NULL);
241: if (flg) { PCASMPrintSubdomains(pc); }
243: if (osm->overlap > 0) {
244: /* Extend the "overlapping" regions by a number of steps */
245: MatIncreaseOverlap(pc->pmat,osm->n_local_true,osm->is,osm->overlap);
246: }
247: if (osm->sort_indices) {
248: for (i=0; i<osm->n_local_true; i++) {
249: ISSort(osm->is[i]);
250: if (osm->is_local) {
251: ISSort(osm->is_local[i]);
252: }
253: }
254: }
256: if (!osm->ksp) {
257: /* Create the local solvers */
258: PetscMalloc1(osm->n_local_true,&osm->ksp);
259: if (domain_dm) {
260: PetscInfo(pc,"Setting up ASM subproblems using the embedded DM\n");
261: }
262: for (i=0; i<osm->n_local_true; i++) {
263: KSPCreate(PETSC_COMM_SELF,&ksp);
264: KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
265: PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
266: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
267: KSPSetType(ksp,KSPPREONLY);
268: KSPGetPC(ksp,&subpc);
269: PCGetOptionsPrefix(pc,&prefix);
270: KSPSetOptionsPrefix(ksp,prefix);
271: KSPAppendOptionsPrefix(ksp,"sub_");
272: if (domain_dm) {
273: KSPSetDM(ksp, domain_dm[i]);
274: KSPSetDMActive(ksp, PETSC_FALSE);
275: DMDestroy(&domain_dm[i]);
276: }
277: osm->ksp[i] = ksp;
278: }
279: if (domain_dm) {
280: PetscFree(domain_dm);
281: }
282: }
284: ISConcatenate(PETSC_COMM_SELF, osm->n_local_true, osm->is, &osm->lis);
285: ISSortRemoveDups(osm->lis);
286: ISGetLocalSize(osm->lis, &m);
288: scall = MAT_INITIAL_MATRIX;
289: } else {
290: /*
291: Destroy the blocks from the previous iteration
292: */
293: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
294: MatDestroyMatrices(osm->n_local_true,&osm->pmat);
295: scall = MAT_INITIAL_MATRIX;
296: }
297: }
299: /* Destroy previous submatrices of a different type than pc->pmat since MAT_REUSE_MATRIX won't work in that case */
300: if ((scall == MAT_REUSE_MATRIX) && osm->sub_mat_type) {
301: if (osm->n_local_true > 0) {
302: MatDestroySubMatrices(osm->n_local_true,&osm->pmat);
303: }
304: scall = MAT_INITIAL_MATRIX;
305: }
307: /*
308: Extract out the submatrices
309: */
310: MatCreateSubMatrices(pc->pmat,osm->n_local_true,osm->is,osm->is,scall,&osm->pmat);
311: if (scall == MAT_INITIAL_MATRIX) {
312: PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
313: for (i=0; i<osm->n_local_true; i++) {
314: PetscLogObjectParent((PetscObject)pc,(PetscObject)osm->pmat[i]);
315: PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);
316: }
317: }
319: /* Convert the types of the submatrices (if needbe) */
320: if (osm->sub_mat_type) {
321: for (i=0; i<osm->n_local_true; i++) {
322: MatConvert(osm->pmat[i],osm->sub_mat_type,MAT_INPLACE_MATRIX,&(osm->pmat[i]));
323: }
324: }
326: if (!pc->setupcalled) {
327: VecType vtype;
329: /* Create the local work vectors (from the local matrices) and scatter contexts */
330: MatCreateVecs(pc->pmat,&vec,NULL);
332: 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()");
333: if (osm->is_local && osm->type == PC_ASM_RESTRICT && osm->loctype == PC_COMPOSITE_ADDITIVE) {
334: PetscMalloc1(osm->n_local_true,&osm->lprolongation);
335: }
336: PetscMalloc1(osm->n_local_true,&osm->lrestriction);
337: PetscMalloc1(osm->n_local_true,&osm->x);
338: PetscMalloc1(osm->n_local_true,&osm->y);
340: ISGetLocalSize(osm->lis,&m);
341: ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
342: MatGetVecType(osm->pmat[0],&vtype);
343: VecCreate(PETSC_COMM_SELF,&osm->lx);
344: VecSetSizes(osm->lx,m,m);
345: VecSetType(osm->lx,vtype);
346: VecDuplicate(osm->lx, &osm->ly);
347: VecScatterCreate(vec,osm->lis,osm->lx,isl,&osm->restriction);
348: ISDestroy(&isl);
350: for (i=0; i<osm->n_local_true; ++i) {
351: ISLocalToGlobalMapping ltog;
352: IS isll;
353: const PetscInt *idx_is;
354: PetscInt *idx_lis,nout;
356: ISGetLocalSize(osm->is[i],&m);
357: MatCreateVecs(osm->pmat[i],&osm->x[i],NULL);
358: VecDuplicate(osm->x[i],&osm->y[i]);
360: /* generate a scatter from ly to y[i] picking all the overlapping is[i] entries */
361: ISLocalToGlobalMappingCreateIS(osm->lis,<og);
362: ISGetLocalSize(osm->is[i],&m);
363: ISGetIndices(osm->is[i], &idx_is);
364: PetscMalloc1(m,&idx_lis);
365: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m,idx_is,&nout,idx_lis);
366: if (nout != m) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is not a subset of lis");
367: ISRestoreIndices(osm->is[i], &idx_is);
368: ISCreateGeneral(PETSC_COMM_SELF,m,idx_lis,PETSC_OWN_POINTER,&isll);
369: ISLocalToGlobalMappingDestroy(<og);
370: ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
371: VecScatterCreate(osm->ly,isll,osm->y[i],isl,&osm->lrestriction[i]);
372: ISDestroy(&isll);
373: ISDestroy(&isl);
374: if (osm->lprolongation) { /* generate a scatter from y[i] to ly picking only the the non-overlapping is_local[i] entries */
375: ISLocalToGlobalMapping ltog;
376: IS isll,isll_local;
377: const PetscInt *idx_local;
378: PetscInt *idx1, *idx2, nout;
380: ISGetLocalSize(osm->is_local[i],&m_local);
381: ISGetIndices(osm->is_local[i], &idx_local);
383: ISLocalToGlobalMappingCreateIS(osm->is[i],<og);
384: PetscMalloc1(m_local,&idx1);
385: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx1);
386: ISLocalToGlobalMappingDestroy(<og);
387: if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of is");
388: ISCreateGeneral(PETSC_COMM_SELF,m_local,idx1,PETSC_OWN_POINTER,&isll);
390: ISLocalToGlobalMappingCreateIS(osm->lis,<og);
391: PetscMalloc1(m_local,&idx2);
392: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx2);
393: ISLocalToGlobalMappingDestroy(<og);
394: if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of lis");
395: ISCreateGeneral(PETSC_COMM_SELF,m_local,idx2,PETSC_OWN_POINTER,&isll_local);
397: ISRestoreIndices(osm->is_local[i], &idx_local);
398: VecScatterCreate(osm->y[i],isll,osm->ly,isll_local,&osm->lprolongation[i]);
400: ISDestroy(&isll);
401: ISDestroy(&isll_local);
402: }
403: }
404: VecDestroy(&vec);
405: }
407: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
408: IS *cis;
409: PetscInt c;
411: PetscMalloc1(osm->n_local_true, &cis);
412: for (c = 0; c < osm->n_local_true; ++c) cis[c] = osm->lis;
413: MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, cis, scall, &osm->lmats);
414: PetscFree(cis);
415: }
417: /* Return control to the user so that the submatrices can be modified (e.g., to apply
418: different boundary conditions for the submatrices than for the global problem) */
419: PCModifySubMatrices(pc,osm->n_local_true,osm->is,osm->is,osm->pmat,pc->modifysubmatricesP);
421: /*
422: Loop over subdomains putting them into local ksp
423: */
424: KSPGetOptionsPrefix(osm->ksp[0],&prefix);
425: for (i=0; i<osm->n_local_true; i++) {
426: KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i]);
427: MatSetOptionsPrefix(osm->pmat[i],prefix);
428: if (!pc->setupcalled) {
429: KSPSetFromOptions(osm->ksp[i]);
430: }
431: }
432: return(0);
433: }
435: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
436: {
437: PC_ASM *osm = (PC_ASM*)pc->data;
438: PetscErrorCode ierr;
439: PetscInt i;
440: KSPConvergedReason reason;
443: for (i=0; i<osm->n_local_true; i++) {
444: KSPSetUp(osm->ksp[i]);
445: KSPGetConvergedReason(osm->ksp[i],&reason);
446: if (reason == KSP_DIVERGED_PC_FAILED) {
447: pc->failedreason = PC_SUBPC_ERROR;
448: }
449: }
450: return(0);
451: }
453: static PetscErrorCode PCApply_ASM(PC pc,Vec x,Vec y)
454: {
455: PC_ASM *osm = (PC_ASM*)pc->data;
457: PetscInt i,n_local_true = osm->n_local_true;
458: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
461: /*
462: support for limiting the restriction or interpolation to only local
463: subdomain values (leaving the other values 0).
464: */
465: if (!(osm->type & PC_ASM_RESTRICT)) {
466: forward = SCATTER_FORWARD_LOCAL;
467: /* have to zero the work RHS since scatter may leave some slots empty */
468: VecSet(osm->lx, 0.0);
469: }
470: if (!(osm->type & PC_ASM_INTERPOLATE)) {
471: reverse = SCATTER_REVERSE_LOCAL;
472: }
474: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE) {
475: /* zero the global and the local solutions */
476: VecSet(y, 0.0);
477: VecSet(osm->ly, 0.0);
479: /* copy the global RHS to local RHS including the ghost nodes */
480: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
481: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
483: /* restrict local RHS to the overlapping 0-block RHS */
484: VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
485: VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
487: /* do the local solves */
488: for (i = 0; i < n_local_true; ++i) {
490: /* solve the overlapping i-block */
491: PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i],0);
492: KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]);
493: KSPCheckSolve(osm->ksp[i], pc, osm->y[i]);
494: PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[i], osm->x[i], osm->y[i], 0);
496: if (osm->lprolongation) { /* interpolate the non-overlapping i-block solution to the local solution (only for restrictive additive) */
497: VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
498: VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
499: } else { /* interpolate the overlapping i-block solution to the local solution */
500: VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
501: VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
502: }
504: if (i < n_local_true-1) {
505: /* restrict local RHS to the overlapping (i+1)-block RHS */
506: VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
507: VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
509: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
510: /* update the overlapping (i+1)-block RHS using the current local solution */
511: MatMult(osm->lmats[i+1], osm->ly, osm->y[i+1]);
512: VecAXPBY(osm->x[i+1],-1.,1., osm->y[i+1]);
513: }
514: }
515: }
516: /* add the local solution to the global solution including the ghost nodes */
517: VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
518: VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
519: } else SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
520: return(0);
521: }
523: static PetscErrorCode PCMatApply_ASM(PC pc,Mat X,Mat Y)
524: {
525: PC_ASM *osm = (PC_ASM*)pc->data;
526: Mat Z,W;
527: Vec x;
528: PetscInt i,m,N;
529: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
533: if (osm->n_local_true > 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Not yet implemented");
534: /*
535: support for limiting the restriction or interpolation to only local
536: subdomain values (leaving the other values 0).
537: */
538: if (!(osm->type & PC_ASM_RESTRICT)) {
539: forward = SCATTER_FORWARD_LOCAL;
540: /* have to zero the work RHS since scatter may leave some slots empty */
541: VecSet(osm->lx, 0.0);
542: }
543: if (!(osm->type & PC_ASM_INTERPOLATE)) {
544: reverse = SCATTER_REVERSE_LOCAL;
545: }
546: VecGetLocalSize(osm->x[0], &m);
547: MatGetSize(X, NULL, &N);
548: MatCreateSeqDense(PETSC_COMM_SELF, m, N, NULL, &Z);
549: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE || osm->loctype == PC_COMPOSITE_ADDITIVE) {
550: /* zero the global and the local solutions */
551: MatZeroEntries(Y);
552: VecSet(osm->ly, 0.0);
554: for (i = 0; i < N; ++i) {
555: MatDenseGetColumnVecRead(X, i, &x);
556: /* copy the global RHS to local RHS including the ghost nodes */
557: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
558: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
559: MatDenseRestoreColumnVecRead(X, i, &x);
561: MatDenseGetColumnVecWrite(Z, i, &x);
562: /* restrict local RHS to the overlapping 0-block RHS */
563: VecScatterBegin(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward);
564: VecScatterEnd(osm->lrestriction[0], osm->lx, x, INSERT_VALUES, forward);
565: MatDenseRestoreColumnVecWrite(Z, i, &x);
566: }
567: MatCreateSeqDense(PETSC_COMM_SELF, m, N, NULL, &W);
568: /* solve the overlapping 0-block */
569: PetscLogEventBegin(PC_ApplyOnBlocks, osm->ksp[0], Z, W, 0);
570: KSPMatSolve(osm->ksp[0], Z, W);
571: KSPCheckSolve(osm->ksp[0], pc, NULL);
572: PetscLogEventEnd(PC_ApplyOnBlocks, osm->ksp[0], Z, W,0);
573: MatDestroy(&Z);
575: for (i = 0; i < N; ++i) {
576: VecSet(osm->ly, 0.0);
577: MatDenseGetColumnVecRead(W, i, &x);
578: if (osm->lprolongation) { /* interpolate the non-overlapping 0-block solution to the local solution (only for restrictive additive) */
579: VecScatterBegin(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward);
580: VecScatterEnd(osm->lprolongation[0], x, osm->ly, ADD_VALUES, forward);
581: } else { /* interpolate the overlapping 0-block solution to the local solution */
582: VecScatterBegin(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse);
583: VecScatterEnd(osm->lrestriction[0], x, osm->ly, ADD_VALUES, reverse);
584: }
585: MatDenseRestoreColumnVecRead(W, i, &x);
587: MatDenseGetColumnVecWrite(Y, i, &x);
588: /* add the local solution to the global solution including the ghost nodes */
589: VecScatterBegin(osm->restriction, osm->ly, x, ADD_VALUES, reverse);
590: VecScatterEnd(osm->restriction, osm->ly, x, ADD_VALUES, reverse);
591: MatDenseRestoreColumnVecWrite(Y, i, &x);
592: }
593: MatDestroy(&W);
594: } else SETERRQ1(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
595: return(0);
596: }
598: static PetscErrorCode PCApplyTranspose_ASM(PC pc,Vec x,Vec y)
599: {
600: PC_ASM *osm = (PC_ASM*)pc->data;
602: PetscInt i,n_local_true = osm->n_local_true;
603: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
606: /*
607: Support for limiting the restriction or interpolation to only local
608: subdomain values (leaving the other values 0).
610: Note: these are reversed from the PCApply_ASM() because we are applying the
611: transpose of the three terms
612: */
614: if (!(osm->type & PC_ASM_INTERPOLATE)) {
615: forward = SCATTER_FORWARD_LOCAL;
616: /* have to zero the work RHS since scatter may leave some slots empty */
617: VecSet(osm->lx, 0.0);
618: }
619: if (!(osm->type & PC_ASM_RESTRICT)) reverse = SCATTER_REVERSE_LOCAL;
621: /* zero the global and the local solutions */
622: VecSet(y, 0.0);
623: VecSet(osm->ly, 0.0);
625: /* Copy the global RHS to local RHS including the ghost nodes */
626: VecScatterBegin(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
627: VecScatterEnd(osm->restriction, x, osm->lx, INSERT_VALUES, forward);
629: /* Restrict local RHS to the overlapping 0-block RHS */
630: VecScatterBegin(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
631: VecScatterEnd(osm->lrestriction[0], osm->lx, osm->x[0], INSERT_VALUES, forward);
633: /* do the local solves */
634: for (i = 0; i < n_local_true; ++i) {
636: /* solve the overlapping i-block */
637: PetscLogEventBegin(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
638: KSPSolveTranspose(osm->ksp[i], osm->x[i], osm->y[i]);
639: KSPCheckSolve(osm->ksp[i],pc,osm->y[i]);
640: PetscLogEventEnd(PC_ApplyOnBlocks,osm->ksp[i],osm->x[i],osm->y[i],0);
642: if (osm->lprolongation) { /* interpolate the non-overlapping i-block solution to the local solution */
643: VecScatterBegin(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
644: VecScatterEnd(osm->lprolongation[i], osm->y[i], osm->ly, ADD_VALUES, forward);
645: } else { /* interpolate the overlapping i-block solution to the local solution */
646: VecScatterBegin(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
647: VecScatterEnd(osm->lrestriction[i], osm->y[i], osm->ly, ADD_VALUES, reverse);
648: }
650: if (i < n_local_true-1) {
651: /* Restrict local RHS to the overlapping (i+1)-block RHS */
652: VecScatterBegin(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
653: VecScatterEnd(osm->lrestriction[i+1], osm->lx, osm->x[i+1], INSERT_VALUES, forward);
654: }
655: }
656: /* Add the local solution to the global solution including the ghost nodes */
657: VecScatterBegin(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
658: VecScatterEnd(osm->restriction, osm->ly, y, ADD_VALUES, reverse);
660: return(0);
662: }
664: static PetscErrorCode PCReset_ASM(PC pc)
665: {
666: PC_ASM *osm = (PC_ASM*)pc->data;
668: PetscInt i;
671: if (osm->ksp) {
672: for (i=0; i<osm->n_local_true; i++) {
673: KSPReset(osm->ksp[i]);
674: }
675: }
676: if (osm->pmat) {
677: if (osm->n_local_true > 0) {
678: MatDestroySubMatrices(osm->n_local_true,&osm->pmat);
679: }
680: }
681: if (osm->lrestriction) {
682: VecScatterDestroy(&osm->restriction);
683: for (i=0; i<osm->n_local_true; i++) {
684: VecScatterDestroy(&osm->lrestriction[i]);
685: if (osm->lprolongation) {VecScatterDestroy(&osm->lprolongation[i]);}
686: VecDestroy(&osm->x[i]);
687: VecDestroy(&osm->y[i]);
688: }
689: PetscFree(osm->lrestriction);
690: if (osm->lprolongation) {PetscFree(osm->lprolongation);}
691: PetscFree(osm->x);
692: PetscFree(osm->y);
694: }
695: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
696: ISDestroy(&osm->lis);
697: VecDestroy(&osm->lx);
698: VecDestroy(&osm->ly);
699: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
700: MatDestroyMatrices(osm->n_local_true, &osm->lmats);
701: }
703: PetscFree(osm->sub_mat_type);
705: osm->is = NULL;
706: osm->is_local = NULL;
707: return(0);
708: }
710: static PetscErrorCode PCDestroy_ASM(PC pc)
711: {
712: PC_ASM *osm = (PC_ASM*)pc->data;
714: PetscInt i;
717: PCReset_ASM(pc);
718: if (osm->ksp) {
719: for (i=0; i<osm->n_local_true; i++) {
720: KSPDestroy(&osm->ksp[i]);
721: }
722: PetscFree(osm->ksp);
723: }
724: PetscFree(pc->data);
726: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",NULL);
727: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",NULL);
728: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetOverlap_C",NULL);
729: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetType_C",NULL);
730: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetType_C",NULL);
731: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalType_C",NULL);
732: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetLocalType_C",NULL);
733: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSortIndices_C",NULL);
734: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubKSP_C",NULL);
735: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubMatType_C",NULL);
736: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSubMatType_C",NULL);
737: return(0);
738: }
740: static PetscErrorCode PCSetFromOptions_ASM(PetscOptionItems *PetscOptionsObject,PC pc)
741: {
742: PC_ASM *osm = (PC_ASM*)pc->data;
744: PetscInt blocks,ovl;
745: PetscBool flg;
746: PCASMType asmtype;
747: PCCompositeType loctype;
748: char sub_mat_type[256];
751: PetscOptionsHead(PetscOptionsObject,"Additive Schwarz options");
752: PetscOptionsBool("-pc_asm_dm_subdomains","Use DMCreateDomainDecomposition() to define subdomains","PCASMSetDMSubdomains",osm->dm_subdomains,&osm->dm_subdomains,&flg);
753: PetscOptionsInt("-pc_asm_blocks","Number of subdomains","PCASMSetTotalSubdomains",osm->n,&blocks,&flg);
754: if (flg) {
755: PCASMSetTotalSubdomains(pc,blocks,NULL,NULL);
756: osm->dm_subdomains = PETSC_FALSE;
757: }
758: PetscOptionsInt("-pc_asm_local_blocks","Number of local subdomains","PCASMSetLocalSubdomains",osm->n_local_true,&blocks,&flg);
759: if (flg) {
760: PCASMSetLocalSubdomains(pc,blocks,NULL,NULL);
761: osm->dm_subdomains = PETSC_FALSE;
762: }
763: PetscOptionsInt("-pc_asm_overlap","Number of grid points overlap","PCASMSetOverlap",osm->overlap,&ovl,&flg);
764: if (flg) {
765: PCASMSetOverlap(pc,ovl);
766: osm->dm_subdomains = PETSC_FALSE;
767: }
768: flg = PETSC_FALSE;
769: PetscOptionsEnum("-pc_asm_type","Type of restriction/extension","PCASMSetType",PCASMTypes,(PetscEnum)osm->type,(PetscEnum*)&asmtype,&flg);
770: if (flg) {PCASMSetType(pc,asmtype); }
771: flg = PETSC_FALSE;
772: PetscOptionsEnum("-pc_asm_local_type","Type of local solver composition","PCASMSetLocalType",PCCompositeTypes,(PetscEnum)osm->loctype,(PetscEnum*)&loctype,&flg);
773: if (flg) {PCASMSetLocalType(pc,loctype); }
774: PetscOptionsFList("-pc_asm_sub_mat_type","Subsolve Matrix Type","PCASMSetSubMatType",MatList,NULL,sub_mat_type,256,&flg);
775: if (flg) {
776: PCASMSetSubMatType(pc,sub_mat_type);
777: }
778: PetscOptionsTail();
779: return(0);
780: }
782: /*------------------------------------------------------------------------------------*/
784: static PetscErrorCode PCASMSetLocalSubdomains_ASM(PC pc,PetscInt n,IS is[],IS is_local[])
785: {
786: PC_ASM *osm = (PC_ASM*)pc->data;
788: PetscInt i;
791: if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Each process must have 1 or more blocks, n = %D",n);
792: if (pc->setupcalled && (n != osm->n_local_true || is)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetLocalSubdomains() should be called before calling PCSetUp().");
794: if (!pc->setupcalled) {
795: if (is) {
796: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is[i]);}
797: }
798: if (is_local) {
799: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is_local[i]);}
800: }
801: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
803: osm->n_local_true = n;
804: osm->is = NULL;
805: osm->is_local = NULL;
806: if (is) {
807: PetscMalloc1(n,&osm->is);
808: for (i=0; i<n; i++) osm->is[i] = is[i];
809: /* Flag indicating that the user has set overlapping subdomains so PCASM should not increase their size. */
810: osm->overlap = -1;
811: }
812: if (is_local) {
813: PetscMalloc1(n,&osm->is_local);
814: for (i=0; i<n; i++) osm->is_local[i] = is_local[i];
815: if (!is) {
816: PetscMalloc1(osm->n_local_true,&osm->is);
817: for (i=0; i<osm->n_local_true; i++) {
818: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
819: ISDuplicate(osm->is_local[i],&osm->is[i]);
820: ISCopy(osm->is_local[i],osm->is[i]);
821: } else {
822: PetscObjectReference((PetscObject)osm->is_local[i]);
823: osm->is[i] = osm->is_local[i];
824: }
825: }
826: }
827: }
828: }
829: return(0);
830: }
832: static PetscErrorCode PCASMSetTotalSubdomains_ASM(PC pc,PetscInt N,IS *is,IS *is_local)
833: {
834: PC_ASM *osm = (PC_ASM*)pc->data;
836: PetscMPIInt rank,size;
837: PetscInt n;
840: if (N < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of total blocks must be > 0, N = %D",N);
841: 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.");
843: /*
844: Split the subdomains equally among all processors
845: */
846: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
847: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
848: n = N/size + ((N % size) > rank);
849: 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);
850: if (pc->setupcalled && n != osm->n_local_true) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCASMSetTotalSubdomains() should be called before PCSetUp().");
851: if (!pc->setupcalled) {
852: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
854: osm->n_local_true = n;
855: osm->is = NULL;
856: osm->is_local = NULL;
857: }
858: return(0);
859: }
861: static PetscErrorCode PCASMSetOverlap_ASM(PC pc,PetscInt ovl)
862: {
863: PC_ASM *osm = (PC_ASM*)pc->data;
866: if (ovl < 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested");
867: if (pc->setupcalled && ovl != osm->overlap) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetOverlap() should be called before PCSetUp().");
868: if (!pc->setupcalled) osm->overlap = ovl;
869: return(0);
870: }
872: static PetscErrorCode PCASMSetType_ASM(PC pc,PCASMType type)
873: {
874: PC_ASM *osm = (PC_ASM*)pc->data;
877: osm->type = type;
878: osm->type_set = PETSC_TRUE;
879: return(0);
880: }
882: static PetscErrorCode PCASMGetType_ASM(PC pc,PCASMType *type)
883: {
884: PC_ASM *osm = (PC_ASM*)pc->data;
887: *type = osm->type;
888: return(0);
889: }
891: static PetscErrorCode PCASMSetLocalType_ASM(PC pc, PCCompositeType type)
892: {
893: PC_ASM *osm = (PC_ASM *) pc->data;
896: 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");
897: osm->loctype = type;
898: return(0);
899: }
901: static PetscErrorCode PCASMGetLocalType_ASM(PC pc, PCCompositeType *type)
902: {
903: PC_ASM *osm = (PC_ASM *) pc->data;
906: *type = osm->loctype;
907: return(0);
908: }
910: static PetscErrorCode PCASMSetSortIndices_ASM(PC pc,PetscBool doSort)
911: {
912: PC_ASM *osm = (PC_ASM*)pc->data;
915: osm->sort_indices = doSort;
916: return(0);
917: }
919: static PetscErrorCode PCASMGetSubKSP_ASM(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
920: {
921: PC_ASM *osm = (PC_ASM*)pc->data;
925: 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");
927: if (n_local) *n_local = osm->n_local_true;
928: if (first_local) {
929: MPI_Scan(&osm->n_local_true,first_local,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
930: *first_local -= osm->n_local_true;
931: }
932: if (ksp) *ksp = osm->ksp;
933: return(0);
934: }
936: static PetscErrorCode PCASMGetSubMatType_ASM(PC pc,MatType *sub_mat_type)
937: {
938: PC_ASM *osm = (PC_ASM*)pc->data;
943: *sub_mat_type = osm->sub_mat_type;
944: return(0);
945: }
947: static PetscErrorCode PCASMSetSubMatType_ASM(PC pc,MatType sub_mat_type)
948: {
949: PetscErrorCode ierr;
950: PC_ASM *osm = (PC_ASM*)pc->data;
954: PetscFree(osm->sub_mat_type);
955: PetscStrallocpy(sub_mat_type,(char**)&osm->sub_mat_type);
956: return(0);
957: }
959: /*@C
960: PCASMSetLocalSubdomains - Sets the local subdomains (for this processor only) for the additive Schwarz preconditioner.
962: Collective on pc
964: Input Parameters:
965: + pc - the preconditioner context
966: . n - the number of subdomains for this processor (default value = 1)
967: . is - the index set that defines the subdomains for this processor
968: (or NULL for PETSc to determine subdomains)
969: - is_local - the index sets that define the local part of the subdomains for this processor, not used unless PCASMType is PC_ASM_RESTRICT
970: (or NULL to not provide these)
972: Options Database Key:
973: To set the total number of subdomain blocks rather than specify the
974: index sets, use the option
975: . -pc_asm_local_blocks <blks> - Sets local blocks
977: Notes:
978: The IS numbering is in the parallel, global numbering of the vector for both is and is_local
980: By default the ASM preconditioner uses 1 block per processor.
982: Use PCASMSetTotalSubdomains() to set the subdomains for all processors.
984: If is_local is provided and PCASMType is PC_ASM_RESTRICT then the solution only over the is_local region is interpolated
985: back to form the global solution (this is the standard restricted additive Schwarz method)
987: If the is_local is provided and PCASMType is PC_ASM_INTERPOLATE or PC_ASM_NONE then an error is generated since there is
988: no code to handle that case.
990: Level: advanced
992: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
993: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), PCASMType, PCASMSetType()
994: @*/
995: PetscErrorCode PCASMSetLocalSubdomains(PC pc,PetscInt n,IS is[],IS is_local[])
996: {
1001: PetscTryMethod(pc,"PCASMSetLocalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,is,is_local));
1002: return(0);
1003: }
1005: /*@C
1006: PCASMSetTotalSubdomains - Sets the subdomains for all processors for the
1007: additive Schwarz preconditioner. Either all or no processors in the
1008: PC communicator must call this routine, with the same index sets.
1010: Collective on pc
1012: Input Parameters:
1013: + pc - the preconditioner context
1014: . N - the number of subdomains for all processors
1015: . is - the index sets that define the subdomains for all processors
1016: (or NULL to ask PETSc to determine the subdomains)
1017: - is_local - the index sets that define the local part of the subdomains for this processor
1018: (or NULL to not provide this information)
1020: Options Database Key:
1021: To set the total number of subdomain blocks rather than specify the
1022: index sets, use the option
1023: . -pc_asm_blocks <blks> - Sets total blocks
1025: Notes:
1026: Currently you cannot use this to set the actual subdomains with the argument is or is_local.
1028: By default the ASM preconditioner uses 1 block per processor.
1030: These index sets cannot be destroyed until after completion of the
1031: linear solves for which the ASM preconditioner is being used.
1033: Use PCASMSetLocalSubdomains() to set local subdomains.
1035: The IS numbering is in the parallel, global numbering of the vector for both is and is_local
1037: Level: advanced
1039: .seealso: PCASMSetLocalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1040: PCASMCreateSubdomains2D()
1041: @*/
1042: PetscErrorCode PCASMSetTotalSubdomains(PC pc,PetscInt N,IS is[],IS is_local[])
1043: {
1048: PetscTryMethod(pc,"PCASMSetTotalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,N,is,is_local));
1049: return(0);
1050: }
1052: /*@
1053: PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
1054: additive Schwarz preconditioner. Either all or no processors in the
1055: PC communicator must call this routine.
1057: Logically Collective on pc
1059: Input Parameters:
1060: + pc - the preconditioner context
1061: - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1)
1063: Options Database Key:
1064: . -pc_asm_overlap <ovl> - Sets overlap
1066: Notes:
1067: By default the ASM preconditioner uses 1 block per processor. To use
1068: multiple blocks per perocessor, see PCASMSetTotalSubdomains() and
1069: PCASMSetLocalSubdomains() (and the option -pc_asm_blocks <blks>).
1071: The overlap defaults to 1, so if one desires that no additional
1072: overlap be computed beyond what may have been set with a call to
1073: PCASMSetTotalSubdomains() or PCASMSetLocalSubdomains(), then ovl
1074: must be set to be 0. In particular, if one does not explicitly set
1075: the subdomains an application code, then all overlap would be computed
1076: internally by PETSc, and using an overlap of 0 would result in an ASM
1077: variant that is equivalent to the block Jacobi preconditioner.
1079: The default algorithm used by PETSc to increase overlap is fast, but not scalable,
1080: use the option -mat_increase_overlap_scalable when the problem and number of processes is large.
1082: Note that one can define initial index sets with any overlap via
1083: PCASMSetLocalSubdomains(); the routine
1084: PCASMSetOverlap() merely allows PETSc to extend that overlap further
1085: if desired.
1087: Level: intermediate
1089: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1090: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), MatIncreaseOverlap()
1091: @*/
1092: PetscErrorCode PCASMSetOverlap(PC pc,PetscInt ovl)
1093: {
1099: PetscTryMethod(pc,"PCASMSetOverlap_C",(PC,PetscInt),(pc,ovl));
1100: return(0);
1101: }
1103: /*@
1104: PCASMSetType - Sets the type of restriction and interpolation used
1105: for local problems in the additive Schwarz method.
1107: Logically Collective on pc
1109: Input Parameters:
1110: + pc - the preconditioner context
1111: - type - variant of ASM, one of
1112: .vb
1113: PC_ASM_BASIC - full interpolation and restriction
1114: PC_ASM_RESTRICT - full restriction, local processor interpolation
1115: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1116: PC_ASM_NONE - local processor restriction and interpolation
1117: .ve
1119: Options Database Key:
1120: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
1122: Notes:
1123: if the is_local arguments are passed to PCASMSetLocalSubdomains() then they are used when PC_ASM_RESTRICT has been selected
1124: to limit the local processor interpolation
1126: Level: intermediate
1128: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1129: PCASMCreateSubdomains2D(), PCASMType, PCASMSetLocalType(), PCASMGetLocalType()
1130: @*/
1131: PetscErrorCode PCASMSetType(PC pc,PCASMType type)
1132: {
1138: PetscTryMethod(pc,"PCASMSetType_C",(PC,PCASMType),(pc,type));
1139: return(0);
1140: }
1142: /*@
1143: PCASMGetType - Gets the type of restriction and interpolation used
1144: for local problems in the additive Schwarz method.
1146: Logically Collective on pc
1148: Input Parameter:
1149: . pc - the preconditioner context
1151: Output Parameter:
1152: . type - variant of ASM, one of
1154: .vb
1155: PC_ASM_BASIC - full interpolation and restriction
1156: PC_ASM_RESTRICT - full restriction, local processor interpolation
1157: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1158: PC_ASM_NONE - local processor restriction and interpolation
1159: .ve
1161: Options Database Key:
1162: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
1164: Level: intermediate
1166: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1167: PCASMCreateSubdomains2D(), PCASMType, PCASMSetType(), PCASMSetLocalType(), PCASMGetLocalType()
1168: @*/
1169: PetscErrorCode PCASMGetType(PC pc,PCASMType *type)
1170: {
1175: PetscUseMethod(pc,"PCASMGetType_C",(PC,PCASMType*),(pc,type));
1176: return(0);
1177: }
1179: /*@
1180: PCASMSetLocalType - Sets the type of composition used for local problems in the additive Schwarz method.
1182: Logically Collective on pc
1184: Input Parameters:
1185: + pc - the preconditioner context
1186: - type - type of composition, one of
1187: .vb
1188: PC_COMPOSITE_ADDITIVE - local additive combination
1189: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1190: .ve
1192: Options Database Key:
1193: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1195: Level: intermediate
1197: .seealso: PCASMSetType(), PCASMGetType(), PCASMGetLocalType(), PCASM, PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1198: @*/
1199: PetscErrorCode PCASMSetLocalType(PC pc, PCCompositeType type)
1200: {
1206: PetscTryMethod(pc, "PCASMSetLocalType_C", (PC, PCCompositeType), (pc, type));
1207: return(0);
1208: }
1210: /*@
1211: PCASMGetLocalType - Gets the type of composition used for local problems in the additive Schwarz method.
1213: Logically Collective on pc
1215: Input Parameter:
1216: . pc - the preconditioner context
1218: Output Parameter:
1219: . type - type of composition, one of
1220: .vb
1221: PC_COMPOSITE_ADDITIVE - local additive combination
1222: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1223: .ve
1225: Options Database Key:
1226: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1228: Level: intermediate
1230: .seealso: PCASMSetType(), PCASMGetType(), PCASMSetLocalType(), PCASMCreate(), PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1231: @*/
1232: PetscErrorCode PCASMGetLocalType(PC pc, PCCompositeType *type)
1233: {
1239: PetscUseMethod(pc, "PCASMGetLocalType_C", (PC, PCCompositeType *), (pc, type));
1240: return(0);
1241: }
1243: /*@
1244: PCASMSetSortIndices - Determines whether subdomain indices are sorted.
1246: Logically Collective on pc
1248: Input Parameters:
1249: + pc - the preconditioner context
1250: - doSort - sort the subdomain indices
1252: Level: intermediate
1254: .seealso: PCASMSetLocalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1255: PCASMCreateSubdomains2D()
1256: @*/
1257: PetscErrorCode PCASMSetSortIndices(PC pc,PetscBool doSort)
1258: {
1264: PetscTryMethod(pc,"PCASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));
1265: return(0);
1266: }
1268: /*@C
1269: PCASMGetSubKSP - Gets the local KSP contexts for all blocks on
1270: this processor.
1272: Collective on pc iff first_local is requested
1274: Input Parameter:
1275: . pc - the preconditioner context
1277: Output Parameters:
1278: + n_local - the number of blocks on this processor or NULL
1279: . first_local - the global number of the first block on this processor or NULL,
1280: all processors must request or all must pass NULL
1281: - ksp - the array of KSP contexts
1283: Note:
1284: After PCASMGetSubKSP() the array of KSPes is not to be freed.
1286: You must call KSPSetUp() before calling PCASMGetSubKSP().
1288: Fortran note:
1289: 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.
1291: Level: advanced
1293: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap(),
1294: PCASMCreateSubdomains2D(),
1295: @*/
1296: PetscErrorCode PCASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
1297: {
1302: PetscUseMethod(pc,"PCASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
1303: return(0);
1304: }
1306: /* -------------------------------------------------------------------------------------*/
1307: /*MC
1308: PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1309: its own KSP object.
1311: Options Database Keys:
1312: + -pc_asm_blocks <blks> - Sets total blocks
1313: . -pc_asm_overlap <ovl> - Sets overlap
1314: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type, default is restrict
1315: - -pc_asm_local_type [additive, multiplicative] - Sets ASM type, default is additive
1317: IMPORTANT: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you
1318: will get a different convergence rate due to the default option of -pc_asm_type restrict. Use
1319: -pc_asm_type basic to use the standard ASM.
1321: Notes:
1322: Each processor can have one or more blocks, but a block cannot be shared by more
1323: than one processor. Use PCGASM for subdomains shared by multiple processes. Defaults to one block per processor.
1325: To set options on the solvers for each block append -sub_ to all the KSP, and PC
1326: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
1328: To set the options on the solvers separate for each block call PCASMGetSubKSP()
1329: and set the options directly on the resulting KSP object (you can access its PC
1330: with KSPGetPC())
1332: Level: beginner
1334: References:
1335: + 1. - M Dryja, OB Widlund, An additive variant of the Schwarz alternating method for the case of many subregions
1336: Courant Institute, New York University Technical report
1337: - 2. - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations,
1338: Cambridge University Press.
1340: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
1341: PCBJACOBI, PCASMGetSubKSP(), PCASMSetLocalSubdomains(), PCASMType, PCASMGetType(), PCASMSetLocalType(), PCASMGetLocalType()
1342: PCASMSetTotalSubdomains(), PCSetModifySubMatrices(), PCASMSetOverlap(), PCASMSetType(), PCCompositeType
1344: M*/
1346: PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1347: {
1349: PC_ASM *osm;
1352: PetscNewLog(pc,&osm);
1354: osm->n = PETSC_DECIDE;
1355: osm->n_local = 0;
1356: osm->n_local_true = PETSC_DECIDE;
1357: osm->overlap = 1;
1358: osm->ksp = NULL;
1359: osm->restriction = NULL;
1360: osm->lprolongation = NULL;
1361: osm->lrestriction = NULL;
1362: osm->x = NULL;
1363: osm->y = NULL;
1364: osm->is = NULL;
1365: osm->is_local = NULL;
1366: osm->mat = NULL;
1367: osm->pmat = NULL;
1368: osm->type = PC_ASM_RESTRICT;
1369: osm->loctype = PC_COMPOSITE_ADDITIVE;
1370: osm->sort_indices = PETSC_TRUE;
1371: osm->dm_subdomains = PETSC_FALSE;
1372: osm->sub_mat_type = NULL;
1374: pc->data = (void*)osm;
1375: pc->ops->apply = PCApply_ASM;
1376: pc->ops->matapply = PCMatApply_ASM;
1377: pc->ops->applytranspose = PCApplyTranspose_ASM;
1378: pc->ops->setup = PCSetUp_ASM;
1379: pc->ops->reset = PCReset_ASM;
1380: pc->ops->destroy = PCDestroy_ASM;
1381: pc->ops->setfromoptions = PCSetFromOptions_ASM;
1382: pc->ops->setuponblocks = PCSetUpOnBlocks_ASM;
1383: pc->ops->view = PCView_ASM;
1384: pc->ops->applyrichardson = NULL;
1386: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",PCASMSetLocalSubdomains_ASM);
1387: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",PCASMSetTotalSubdomains_ASM);
1388: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetOverlap_C",PCASMSetOverlap_ASM);
1389: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetType_C",PCASMSetType_ASM);
1390: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetType_C",PCASMGetType_ASM);
1391: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalType_C",PCASMSetLocalType_ASM);
1392: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetLocalType_C",PCASMGetLocalType_ASM);
1393: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSortIndices_C",PCASMSetSortIndices_ASM);
1394: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubKSP_C",PCASMGetSubKSP_ASM);
1395: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubMatType_C",PCASMGetSubMatType_ASM);
1396: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSubMatType_C",PCASMSetSubMatType_ASM);
1397: return(0);
1398: }
1400: /*@C
1401: PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1402: preconditioner for any problem on a general grid.
1404: Collective
1406: Input Parameters:
1407: + A - The global matrix operator
1408: - n - the number of local blocks
1410: Output Parameters:
1411: . outis - the array of index sets defining the subdomains
1413: Level: advanced
1415: Note: this generates nonoverlapping subdomains; the PCASM will generate the overlap
1416: from these if you use PCASMSetLocalSubdomains()
1418: In the Fortran version you must provide the array outis[] already allocated of length n.
1420: .seealso: PCASMSetLocalSubdomains(), PCASMDestroySubdomains()
1421: @*/
1422: PetscErrorCode PCASMCreateSubdomains(Mat A, PetscInt n, IS* outis[])
1423: {
1424: MatPartitioning mpart;
1425: const char *prefix;
1426: PetscInt i,j,rstart,rend,bs;
1427: PetscBool hasop, isbaij = PETSC_FALSE,foundpart = PETSC_FALSE;
1428: Mat Ad = NULL, adj;
1429: IS ispart,isnumb,*is;
1430: PetscErrorCode ierr;
1435: if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of local blocks must be > 0, n = %D",n);
1437: /* Get prefix, row distribution, and block size */
1438: MatGetOptionsPrefix(A,&prefix);
1439: MatGetOwnershipRange(A,&rstart,&rend);
1440: MatGetBlockSize(A,&bs);
1441: 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);
1443: /* Get diagonal block from matrix if possible */
1444: MatHasOperation(A,MATOP_GET_DIAGONAL_BLOCK,&hasop);
1445: if (hasop) {
1446: MatGetDiagonalBlock(A,&Ad);
1447: }
1448: if (Ad) {
1449: PetscObjectBaseTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij);
1450: if (!isbaij) {PetscObjectBaseTypeCompare((PetscObject)Ad,MATSEQSBAIJ,&isbaij);}
1451: }
1452: if (Ad && n > 1) {
1453: PetscBool match,done;
1454: /* Try to setup a good matrix partitioning if available */
1455: MatPartitioningCreate(PETSC_COMM_SELF,&mpart);
1456: PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);
1457: MatPartitioningSetFromOptions(mpart);
1458: PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGCURRENT,&match);
1459: if (!match) {
1460: PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGSQUARE,&match);
1461: }
1462: if (!match) { /* assume a "good" partitioner is available */
1463: PetscInt na;
1464: const PetscInt *ia,*ja;
1465: MatGetRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1466: if (done) {
1467: /* Build adjacency matrix by hand. Unfortunately a call to
1468: MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1469: remove the block-aij structure and we cannot expect
1470: MatPartitioning to split vertices as we need */
1471: PetscInt i,j,len,nnz,cnt,*iia=NULL,*jja=NULL;
1472: const PetscInt *row;
1473: nnz = 0;
1474: for (i=0; i<na; i++) { /* count number of nonzeros */
1475: len = ia[i+1] - ia[i];
1476: row = ja + ia[i];
1477: for (j=0; j<len; j++) {
1478: if (row[j] == i) { /* don't count diagonal */
1479: len--; break;
1480: }
1481: }
1482: nnz += len;
1483: }
1484: PetscMalloc1(na+1,&iia);
1485: PetscMalloc1(nnz,&jja);
1486: nnz = 0;
1487: iia[0] = 0;
1488: for (i=0; i<na; i++) { /* fill adjacency */
1489: cnt = 0;
1490: len = ia[i+1] - ia[i];
1491: row = ja + ia[i];
1492: for (j=0; j<len; j++) {
1493: if (row[j] != i) { /* if not diagonal */
1494: jja[nnz+cnt++] = row[j];
1495: }
1496: }
1497: nnz += cnt;
1498: iia[i+1] = nnz;
1499: }
1500: /* Partitioning of the adjacency matrix */
1501: MatCreateMPIAdj(PETSC_COMM_SELF,na,na,iia,jja,NULL,&adj);
1502: MatPartitioningSetAdjacency(mpart,adj);
1503: MatPartitioningSetNParts(mpart,n);
1504: MatPartitioningApply(mpart,&ispart);
1505: ISPartitioningToNumbering(ispart,&isnumb);
1506: MatDestroy(&adj);
1507: foundpart = PETSC_TRUE;
1508: }
1509: MatRestoreRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1510: }
1511: MatPartitioningDestroy(&mpart);
1512: }
1514: PetscMalloc1(n,&is);
1515: *outis = is;
1517: if (!foundpart) {
1519: /* Partitioning by contiguous chunks of rows */
1521: PetscInt mbs = (rend-rstart)/bs;
1522: PetscInt start = rstart;
1523: for (i=0; i<n; i++) {
1524: PetscInt count = (mbs/n + ((mbs % n) > i)) * bs;
1525: ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);
1526: start += count;
1527: }
1529: } else {
1531: /* Partitioning by adjacency of diagonal block */
1533: const PetscInt *numbering;
1534: PetscInt *count,nidx,*indices,*newidx,start=0;
1535: /* Get node count in each partition */
1536: PetscMalloc1(n,&count);
1537: ISPartitioningCount(ispart,n,count);
1538: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1539: for (i=0; i<n; i++) count[i] *= bs;
1540: }
1541: /* Build indices from node numbering */
1542: ISGetLocalSize(isnumb,&nidx);
1543: PetscMalloc1(nidx,&indices);
1544: for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */
1545: ISGetIndices(isnumb,&numbering);
1546: PetscSortIntWithPermutation(nidx,numbering,indices);
1547: ISRestoreIndices(isnumb,&numbering);
1548: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1549: PetscMalloc1(nidx*bs,&newidx);
1550: for (i=0; i<nidx; i++) {
1551: for (j=0; j<bs; j++) newidx[i*bs+j] = indices[i]*bs + j;
1552: }
1553: PetscFree(indices);
1554: nidx *= bs;
1555: indices = newidx;
1556: }
1557: /* Shift to get global indices */
1558: for (i=0; i<nidx; i++) indices[i] += rstart;
1560: /* Build the index sets for each block */
1561: for (i=0; i<n; i++) {
1562: ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);
1563: ISSort(is[i]);
1564: start += count[i];
1565: }
1567: PetscFree(count);
1568: PetscFree(indices);
1569: ISDestroy(&isnumb);
1570: ISDestroy(&ispart);
1572: }
1573: return(0);
1574: }
1576: /*@C
1577: PCASMDestroySubdomains - Destroys the index sets created with
1578: PCASMCreateSubdomains(). Should be called after setting subdomains
1579: with PCASMSetLocalSubdomains().
1581: Collective
1583: Input Parameters:
1584: + n - the number of index sets
1585: . is - the array of index sets
1586: - is_local - the array of local index sets, can be NULL
1588: Level: advanced
1590: .seealso: PCASMCreateSubdomains(), PCASMSetLocalSubdomains()
1591: @*/
1592: PetscErrorCode PCASMDestroySubdomains(PetscInt n, IS is[], IS is_local[])
1593: {
1594: PetscInt i;
1598: if (n <= 0) return(0);
1599: if (is) {
1601: for (i=0; i<n; i++) { ISDestroy(&is[i]); }
1602: PetscFree(is);
1603: }
1604: if (is_local) {
1606: for (i=0; i<n; i++) { ISDestroy(&is_local[i]); }
1607: PetscFree(is_local);
1608: }
1609: return(0);
1610: }
1612: /*@
1613: PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1614: preconditioner for a two-dimensional problem on a regular grid.
1616: Not Collective
1618: Input Parameters:
1619: + m, n - the number of mesh points in the x and y directions
1620: . M, N - the number of subdomains in the x and y directions
1621: . dof - degrees of freedom per node
1622: - overlap - overlap in mesh lines
1624: Output Parameters:
1625: + Nsub - the number of subdomains created
1626: . is - array of index sets defining overlapping (if overlap > 0) subdomains
1627: - is_local - array of index sets defining non-overlapping subdomains
1629: Note:
1630: Presently PCAMSCreateSubdomains2d() is valid only for sequential
1631: preconditioners. More general related routines are
1632: PCASMSetTotalSubdomains() and PCASMSetLocalSubdomains().
1634: Level: advanced
1636: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1637: PCASMSetOverlap()
1638: @*/
1639: PetscErrorCode PCASMCreateSubdomains2D(PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt dof,PetscInt overlap,PetscInt *Nsub,IS **is,IS **is_local)
1640: {
1641: PetscInt i,j,height,width,ystart,xstart,yleft,yright,xleft,xright,loc_outer;
1643: PetscInt nidx,*idx,loc,ii,jj,count;
1646: if (dof != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP," ");
1648: *Nsub = N*M;
1649: PetscMalloc1(*Nsub,is);
1650: PetscMalloc1(*Nsub,is_local);
1651: ystart = 0;
1652: loc_outer = 0;
1653: for (i=0; i<N; i++) {
1654: height = n/N + ((n % N) > i); /* height of subdomain */
1655: if (height < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many N subdomains for mesh dimension n");
1656: yleft = ystart - overlap; if (yleft < 0) yleft = 0;
1657: yright = ystart + height + overlap; if (yright > n) yright = n;
1658: xstart = 0;
1659: for (j=0; j<M; j++) {
1660: width = m/M + ((m % M) > j); /* width of subdomain */
1661: if (width < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many M subdomains for mesh dimension m");
1662: xleft = xstart - overlap; if (xleft < 0) xleft = 0;
1663: xright = xstart + width + overlap; if (xright > m) xright = m;
1664: nidx = (xright - xleft)*(yright - yleft);
1665: PetscMalloc1(nidx,&idx);
1666: loc = 0;
1667: for (ii=yleft; ii<yright; ii++) {
1668: count = m*ii + xleft;
1669: for (jj=xleft; jj<xright; jj++) idx[loc++] = count++;
1670: }
1671: ISCreateGeneral(PETSC_COMM_SELF,nidx,idx,PETSC_COPY_VALUES,(*is)+loc_outer);
1672: if (overlap == 0) {
1673: PetscObjectReference((PetscObject)(*is)[loc_outer]);
1675: (*is_local)[loc_outer] = (*is)[loc_outer];
1676: } else {
1677: for (loc=0,ii=ystart; ii<ystart+height; ii++) {
1678: for (jj=xstart; jj<xstart+width; jj++) {
1679: idx[loc++] = m*ii + jj;
1680: }
1681: }
1682: ISCreateGeneral(PETSC_COMM_SELF,loc,idx,PETSC_COPY_VALUES,*is_local+loc_outer);
1683: }
1684: PetscFree(idx);
1685: xstart += width;
1686: loc_outer++;
1687: }
1688: ystart += height;
1689: }
1690: for (i=0; i<*Nsub; i++) { ISSort((*is)[i]); }
1691: return(0);
1692: }
1694: /*@C
1695: PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
1696: only) for the additive Schwarz preconditioner.
1698: Not Collective
1700: Input Parameter:
1701: . pc - the preconditioner context
1703: Output Parameters:
1704: + n - if requested, the number of subdomains for this processor (default value = 1)
1705: . is - if requested, the index sets that define the subdomains for this processor
1706: - is_local - if requested, the index sets that define the local part of the subdomains for this processor (can be NULL)
1709: Notes:
1710: The IS numbering is in the parallel, global numbering of the vector.
1712: Level: advanced
1714: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1715: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubmatrices()
1716: @*/
1717: PetscErrorCode PCASMGetLocalSubdomains(PC pc,PetscInt *n,IS *is[],IS *is_local[])
1718: {
1719: PC_ASM *osm = (PC_ASM*)pc->data;
1721: PetscBool match;
1728: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1729: if (!match) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PC is not a PCASM");
1730: if (n) *n = osm->n_local_true;
1731: if (is) *is = osm->is;
1732: if (is_local) *is_local = osm->is_local;
1733: return(0);
1734: }
1736: /*@C
1737: PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1738: only) for the additive Schwarz preconditioner.
1740: Not Collective
1742: Input Parameter:
1743: . pc - the preconditioner context
1745: Output Parameters:
1746: + n - if requested, the number of matrices for this processor (default value = 1)
1747: - mat - if requested, the matrices
1749: Level: advanced
1751: Notes:
1752: Call after PCSetUp() (or KSPSetUp()) but before PCApply() and before PCSetUpOnBlocks())
1754: Usually one would use PCSetModifySubMatrices() to change the submatrices in building the preconditioner.
1756: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1757: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains(), PCSetModifySubMatrices()
1758: @*/
1759: PetscErrorCode PCASMGetLocalSubmatrices(PC pc,PetscInt *n,Mat *mat[])
1760: {
1761: PC_ASM *osm;
1763: PetscBool match;
1769: if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUp() or PCSetUp().");
1770: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1771: if (!match) {
1772: if (n) *n = 0;
1773: if (mat) *mat = NULL;
1774: } else {
1775: osm = (PC_ASM*)pc->data;
1776: if (n) *n = osm->n_local_true;
1777: if (mat) *mat = osm->pmat;
1778: }
1779: return(0);
1780: }
1782: /*@
1783: PCASMSetDMSubdomains - Indicates whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1785: Logically Collective
1787: Input Parameter:
1788: + pc - the preconditioner
1789: - flg - boolean indicating whether to use subdomains defined by the DM
1791: Options Database Key:
1792: . -pc_asm_dm_subdomains
1794: Level: intermediate
1796: Notes:
1797: PCASMSetTotalSubdomains() and PCASMSetOverlap() take precedence over PCASMSetDMSubdomains(),
1798: so setting either of the first two effectively turns the latter off.
1800: .seealso: PCASMGetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1801: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1802: @*/
1803: PetscErrorCode PCASMSetDMSubdomains(PC pc,PetscBool flg)
1804: {
1805: PC_ASM *osm = (PC_ASM*)pc->data;
1807: PetscBool match;
1812: if (pc->setupcalled) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for a setup PC.");
1813: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1814: if (match) {
1815: osm->dm_subdomains = flg;
1816: }
1817: return(0);
1818: }
1820: /*@
1821: PCASMGetDMSubdomains - Returns flag indicating whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1822: Not Collective
1824: Input Parameter:
1825: . pc - the preconditioner
1827: Output Parameter:
1828: . flg - boolean indicating whether to use subdomains defined by the DM
1830: Level: intermediate
1832: .seealso: PCASMSetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1833: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1834: @*/
1835: PetscErrorCode PCASMGetDMSubdomains(PC pc,PetscBool* flg)
1836: {
1837: PC_ASM *osm = (PC_ASM*)pc->data;
1839: PetscBool match;
1844: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1845: if (match) *flg = osm->dm_subdomains;
1846: else *flg = PETSC_FALSE;
1847: return(0);
1848: }
1850: /*@
1851: PCASMGetSubMatType - Gets the matrix type used for ASM subsolves, as a string.
1853: Not Collective
1855: Input Parameter:
1856: . pc - the PC
1858: Output Parameter:
1859: . -pc_asm_sub_mat_type - name of matrix type
1861: Level: advanced
1863: .seealso: PCASMSetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1864: @*/
1865: PetscErrorCode PCASMGetSubMatType(PC pc,MatType *sub_mat_type)
1866: {
1871: PetscTryMethod(pc,"PCASMGetSubMatType_C",(PC,MatType*),(pc,sub_mat_type));
1872: return(0);
1873: }
1875: /*@
1876: PCASMSetSubMatType - Set the type of matrix used for ASM subsolves
1878: Collective on Mat
1880: Input Parameters:
1881: + pc - the PC object
1882: - sub_mat_type - matrix type
1884: Options Database Key:
1885: . -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.
1887: Notes:
1888: See "${PETSC_DIR}/include/petscmat.h" for available types
1890: Level: advanced
1892: .seealso: PCASMGetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1893: @*/
1894: PetscErrorCode PCASMSetSubMatType(PC pc,MatType sub_mat_type)
1895: {
1900: PetscTryMethod(pc,"PCASMSetSubMatType_C",(PC,MatType),(pc,sub_mat_type));
1901: return(0);
1902: }