Actual source code: asm.c
petsc-3.8.4 2018-03-24
2: /*
3: This file defines an additive Schwarz preconditioner for any Mat implementation.
5: Note that each processor may have any number of subdomains. But in order to
6: deal easily with the VecScatter(), we treat each processor as if it has the
7: same number of subdomains.
9: n - total number of true subdomains on all processors
10: n_local_true - actual number of subdomains on this processor
11: n_local = maximum over all processors of n_local_true
12: */
13: #include <petsc/private/pcimpl.h>
14: #include <petscdm.h>
16: typedef struct {
17: PetscInt n, n_local, n_local_true;
18: PetscInt overlap; /* overlap requested by user */
19: KSP *ksp; /* linear solvers for each block */
20: VecScatter *restriction; /* mapping from global to subregion */
21: VecScatter *localization; /* mapping from overlapping to non-overlapping subregion; used for restrict version of algorithms */
22: VecScatter *prolongation; /* mapping from subregion to global */
23: Vec *x,*y,*y_local; /* work vectors */
24: IS *is; /* index set that defines each overlapping subdomain */
25: IS *is_local; /* index set that defines each non-overlapping subdomain, may be NULL */
26: Mat *mat,*pmat; /* mat is not currently used */
27: PCASMType type; /* use reduced interpolation, restriction or both */
28: PetscBool type_set; /* if user set this value (so won't change it for symmetric problems) */
29: PetscBool same_local_solves; /* flag indicating whether all local solvers are same */
30: PetscBool sort_indices; /* flag to sort subdomain indices */
31: PetscBool dm_subdomains; /* whether DM is allowed to define subdomains */
32: PCCompositeType loctype; /* the type of composition for local solves */
33: MatType sub_mat_type; /* the type of Mat used for subdomain solves (can be MATSAME or NULL) */
34: /* For multiplicative solve */
35: Mat *lmats; /* submatrices for overlapping multiplicative (process) subdomain */
36: Vec lx, ly; /* work vectors */
37: IS lis; /* index set that defines each overlapping multiplicative (process) subdomain */
38: } PC_ASM;
40: static PetscErrorCode PCView_ASM(PC pc,PetscViewer viewer)
41: {
42: PC_ASM *osm = (PC_ASM*)pc->data;
44: PetscMPIInt rank;
45: PetscInt i,bsz;
46: PetscBool iascii,isstring;
47: PetscViewer sviewer;
50: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
51: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERSTRING,&isstring);
52: if (iascii) {
53: char overlaps[256] = "user-defined overlap",blocks[256] = "total subdomain blocks not yet set";
54: if (osm->overlap >= 0) {PetscSNPrintf(overlaps,sizeof(overlaps),"amount of overlap = %D",osm->overlap);}
55: if (osm->n > 0) {PetscSNPrintf(blocks,sizeof(blocks),"total subdomain blocks = %D",osm->n);}
56: PetscViewerASCIIPrintf(viewer," %s, %s\n",blocks,overlaps);
57: PetscViewerASCIIPrintf(viewer," restriction/interpolation type - %s\n",PCASMTypes[osm->type]);
58: if (osm->dm_subdomains) {PetscViewerASCIIPrintf(viewer," Additive Schwarz: using DM to define subdomains\n");}
59: if (osm->loctype != PC_COMPOSITE_ADDITIVE) {PetscViewerASCIIPrintf(viewer," Additive Schwarz: local solve composition type - %s\n",PCCompositeTypes[osm->loctype]);}
60: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
61: if (osm->same_local_solves) {
62: if (osm->ksp) {
63: PetscViewerASCIIPrintf(viewer," Local solve is same for all blocks, in the following KSP and PC objects:\n");
64: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
65: if (!rank) {
66: PetscViewerASCIIPushTab(viewer);
67: KSPView(osm->ksp[0],sviewer);
68: PetscViewerASCIIPopTab(viewer);
69: }
70: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
71: }
72: } else {
73: PetscViewerASCIIPushSynchronized(viewer);
74: PetscViewerASCIISynchronizedPrintf(viewer," [%d] number of local blocks = %D\n",(int)rank,osm->n_local_true);
75: PetscViewerFlush(viewer);
76: PetscViewerASCIIPrintf(viewer," Local solve info for each block is in the following KSP and PC objects:\n");
77: PetscViewerASCIIPushTab(viewer);
78: PetscViewerASCIIPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
79: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
80: for (i=0; i<osm->n_local_true; i++) {
81: ISGetLocalSize(osm->is[i],&bsz);
82: PetscViewerASCIISynchronizedPrintf(sviewer,"[%d] local block number %D, size = %D\n",(int)rank,i,bsz);
83: KSPView(osm->ksp[i],sviewer);
84: PetscViewerASCIISynchronizedPrintf(sviewer,"- - - - - - - - - - - - - - - - - -\n");
85: }
86: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
87: PetscViewerASCIIPopTab(viewer);
88: PetscViewerFlush(viewer);
89: PetscViewerASCIIPopSynchronized(viewer);
90: }
91: } else if (isstring) {
92: PetscViewerStringSPrintf(viewer," blocks=%D, overlap=%D, type=%s",osm->n,osm->overlap,PCASMTypes[osm->type]);
93: PetscViewerGetSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
94: if (osm->ksp) {KSPView(osm->ksp[0],sviewer);}
95: PetscViewerRestoreSubViewer(viewer,PETSC_COMM_SELF,&sviewer);
96: }
97: return(0);
98: }
100: static PetscErrorCode PCASMPrintSubdomains(PC pc)
101: {
102: PC_ASM *osm = (PC_ASM*)pc->data;
103: const char *prefix;
104: char fname[PETSC_MAX_PATH_LEN+1];
105: PetscViewer viewer, sviewer;
106: char *s;
107: PetscInt i,j,nidx;
108: const PetscInt *idx;
109: PetscMPIInt rank, size;
113: MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size);
114: MPI_Comm_rank(PetscObjectComm((PetscObject)pc), &rank);
115: PCGetOptionsPrefix(pc,&prefix);
116: PetscOptionsGetString(NULL,prefix,"-pc_asm_print_subdomains",fname,PETSC_MAX_PATH_LEN,NULL);
117: if (fname[0] == 0) { PetscStrcpy(fname,"stdout"); };
118: PetscViewerASCIIOpen(PetscObjectComm((PetscObject)pc),fname,&viewer);
119: for (i=0; i<osm->n_local; i++) {
120: if (i < osm->n_local_true) {
121: ISGetLocalSize(osm->is[i],&nidx);
122: ISGetIndices(osm->is[i],&idx);
123: /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
124: PetscMalloc1(16*(nidx+1)+512, &s);
125: PetscViewerStringOpen(PETSC_COMM_SELF, s, 16*(nidx+1)+512, &sviewer);
126: PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D with overlap:\n", rank, size, i);
127: for (j=0; j<nidx; j++) {
128: PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
129: }
130: ISRestoreIndices(osm->is[i],&idx);
131: PetscViewerStringSPrintf(sviewer,"\n");
132: PetscViewerDestroy(&sviewer);
133: PetscViewerASCIIPushSynchronized(viewer);
134: PetscViewerASCIISynchronizedPrintf(viewer, s);
135: PetscViewerFlush(viewer);
136: PetscViewerASCIIPopSynchronized(viewer);
137: PetscFree(s);
138: if (osm->is_local) {
139: /* Print to a string viewer; no more than 15 characters per index plus 512 char for the header.*/
140: PetscMalloc1(16*(nidx+1)+512, &s);
141: PetscViewerStringOpen(PETSC_COMM_SELF, s, 16*(nidx+1)+512, &sviewer);
142: PetscViewerStringSPrintf(sviewer, "[%D:%D] Subdomain %D without overlap:\n", rank, size, i);
143: ISGetLocalSize(osm->is_local[i],&nidx);
144: ISGetIndices(osm->is_local[i],&idx);
145: for (j=0; j<nidx; j++) {
146: PetscViewerStringSPrintf(sviewer,"%D ",idx[j]);
147: }
148: ISRestoreIndices(osm->is_local[i],&idx);
149: PetscViewerStringSPrintf(sviewer,"\n");
150: PetscViewerDestroy(&sviewer);
151: PetscViewerASCIIPushSynchronized(viewer);
152: PetscViewerASCIISynchronizedPrintf(viewer, s);
153: PetscViewerFlush(viewer);
154: PetscViewerASCIIPopSynchronized(viewer);
155: PetscFree(s);
156: }
157: } else {
158: /* Participate in collective viewer calls. */
159: PetscViewerASCIIPushSynchronized(viewer);
160: PetscViewerFlush(viewer);
161: PetscViewerASCIIPopSynchronized(viewer);
162: /* Assume either all ranks have is_local or none do. */
163: if (osm->is_local) {
164: PetscViewerASCIIPushSynchronized(viewer);
165: PetscViewerFlush(viewer);
166: PetscViewerASCIIPopSynchronized(viewer);
167: }
168: }
169: }
170: PetscViewerFlush(viewer);
171: PetscViewerDestroy(&viewer);
172: return(0);
173: }
175: static PetscErrorCode PCSetUp_ASM(PC pc)
176: {
177: PC_ASM *osm = (PC_ASM*)pc->data;
179: PetscBool symset,flg;
180: PetscInt i,m,m_local;
181: MatReuse scall = MAT_REUSE_MATRIX;
182: IS isl;
183: KSP ksp;
184: PC subpc;
185: const char *prefix,*pprefix;
186: Vec vec;
187: DM *domain_dm = NULL;
190: if (!pc->setupcalled) {
192: if (!osm->type_set) {
193: MatIsSymmetricKnown(pc->pmat,&symset,&flg);
194: if (symset && flg) osm->type = PC_ASM_BASIC;
195: }
197: /* Note: if subdomains have been set either via PCASMSetTotalSubdomains() or via PCASMSetLocalSubdomains(), osm->n_local_true will not be PETSC_DECIDE */
198: if (osm->n_local_true == PETSC_DECIDE) {
199: /* no subdomains given */
200: /* try pc->dm first, if allowed */
201: if (osm->dm_subdomains && pc->dm) {
202: PetscInt num_domains, d;
203: char **domain_names;
204: IS *inner_domain_is, *outer_domain_is;
205: DMCreateDomainDecomposition(pc->dm, &num_domains, &domain_names, &inner_domain_is, &outer_domain_is, &domain_dm);
206: osm->overlap = -1; /* We do not want to increase the overlap of the IS.
207: A future improvement of this code might allow one to use
208: DM-defined subdomains and also increase the overlap,
209: but that is not currently supported */
210: if (num_domains) {
211: PCASMSetLocalSubdomains(pc, num_domains, outer_domain_is, inner_domain_is);
212: }
213: for (d = 0; d < num_domains; ++d) {
214: if (domain_names) {PetscFree(domain_names[d]);}
215: if (inner_domain_is) {ISDestroy(&inner_domain_is[d]);}
216: if (outer_domain_is) {ISDestroy(&outer_domain_is[d]);}
217: }
218: PetscFree(domain_names);
219: PetscFree(inner_domain_is);
220: PetscFree(outer_domain_is);
221: }
222: if (osm->n_local_true == PETSC_DECIDE) {
223: /* still no subdomains; use one subdomain per processor */
224: osm->n_local_true = 1;
225: }
226: }
227: { /* determine the global and max number of subdomains */
228: struct {PetscInt max,sum;} inwork,outwork;
229: PetscMPIInt size;
231: inwork.max = osm->n_local_true;
232: inwork.sum = osm->n_local_true;
233: MPIU_Allreduce(&inwork,&outwork,1,MPIU_2INT,MPIU_MAXSUM_OP,PetscObjectComm((PetscObject)pc));
234: osm->n_local = outwork.max;
235: osm->n = outwork.sum;
237: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
238: if (outwork.max == 1 && outwork.sum == size) {
239: /* osm->n_local_true = 1 on all processes, set this option may enable use of optimized MatCreateSubMatrices() implementation */
240: MatSetOption(pc->pmat,MAT_SUBMAT_SINGLEIS,PETSC_TRUE);
241: }
242: }
243: if (!osm->is) { /* create the index sets */
244: PCASMCreateSubdomains(pc->pmat,osm->n_local_true,&osm->is);
245: }
246: if (osm->n_local_true > 1 && !osm->is_local) {
247: PetscMalloc1(osm->n_local_true,&osm->is_local);
248: for (i=0; i<osm->n_local_true; i++) {
249: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
250: ISDuplicate(osm->is[i],&osm->is_local[i]);
251: ISCopy(osm->is[i],osm->is_local[i]);
252: } else {
253: PetscObjectReference((PetscObject)osm->is[i]);
254: osm->is_local[i] = osm->is[i];
255: }
256: }
257: }
258: PCGetOptionsPrefix(pc,&prefix);
259: flg = PETSC_FALSE;
260: PetscOptionsGetBool(NULL,prefix,"-pc_asm_print_subdomains",&flg,NULL);
261: if (flg) { PCASMPrintSubdomains(pc); }
263: if (osm->overlap > 0) {
264: /* Extend the "overlapping" regions by a number of steps */
265: MatIncreaseOverlap(pc->pmat,osm->n_local_true,osm->is,osm->overlap);
266: }
267: if (osm->sort_indices) {
268: for (i=0; i<osm->n_local_true; i++) {
269: ISSort(osm->is[i]);
270: if (osm->is_local) {
271: ISSort(osm->is_local[i]);
272: }
273: }
274: }
276: if (!osm->ksp) {
277: /* Create the local solvers */
278: PetscMalloc1(osm->n_local_true,&osm->ksp);
279: if (domain_dm) {
280: PetscInfo(pc,"Setting up ASM subproblems using the embedded DM\n");
281: }
282: for (i=0; i<osm->n_local_true; i++) {
283: KSPCreate(PETSC_COMM_SELF,&ksp);
284: KSPSetErrorIfNotConverged(ksp,pc->erroriffailure);
285: PetscLogObjectParent((PetscObject)pc,(PetscObject)ksp);
286: PetscObjectIncrementTabLevel((PetscObject)ksp,(PetscObject)pc,1);
287: KSPSetType(ksp,KSPPREONLY);
288: KSPGetPC(ksp,&subpc);
289: PCGetOptionsPrefix(pc,&prefix);
290: KSPSetOptionsPrefix(ksp,prefix);
291: KSPAppendOptionsPrefix(ksp,"sub_");
292: if (domain_dm) {
293: KSPSetDM(ksp, domain_dm[i]);
294: KSPSetDMActive(ksp, PETSC_FALSE);
295: DMDestroy(&domain_dm[i]);
296: }
297: osm->ksp[i] = ksp;
298: }
299: if (domain_dm) {
300: PetscFree(domain_dm);
301: }
302: }
303: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
304: PetscInt m;
306: ISConcatenate(PETSC_COMM_SELF, osm->n_local_true, osm->is, &osm->lis);
307: ISSortRemoveDups(osm->lis);
308: ISGetLocalSize(osm->lis, &m);
309: VecCreateSeq(PETSC_COMM_SELF, m, &osm->lx);
310: VecDuplicate(osm->lx, &osm->ly);
311: }
312: scall = MAT_INITIAL_MATRIX;
313: } else {
314: /*
315: Destroy the blocks from the previous iteration
316: */
317: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
318: MatDestroyMatrices(osm->n_local_true,&osm->pmat);
319: scall = MAT_INITIAL_MATRIX;
320: }
321: }
323: /*
324: Extract out the submatrices
325: */
326: MatCreateSubMatrices(pc->pmat,osm->n_local_true,osm->is,osm->is,scall,&osm->pmat);
327: if (scall == MAT_INITIAL_MATRIX) {
328: PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
329: for (i=0; i<osm->n_local_true; i++) {
330: PetscLogObjectParent((PetscObject)pc,(PetscObject)osm->pmat[i]);
331: PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);
332: }
333: }
335: /* Convert the types of the submatrices (if needbe) */
336: if (osm->sub_mat_type) {
337: for (i=0; i<osm->n_local_true; i++) {
338: MatConvert(osm->pmat[i],osm->sub_mat_type,MAT_INPLACE_MATRIX,&(osm->pmat[i]));
339: }
340: }
342: if(!pc->setupcalled){
343: /* Create the local work vectors (from the local matrices) and scatter contexts */
344: MatCreateVecs(pc->pmat,&vec,0);
345: PetscMalloc1(osm->n_local,&osm->restriction);
346: 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()");
347: if (osm->is_local && osm->type == PC_ASM_RESTRICT) {PetscMalloc1(osm->n_local,&osm->localization);}
348: PetscMalloc1(osm->n_local,&osm->prolongation);
349: PetscMalloc1(osm->n_local,&osm->x);
350: PetscMalloc1(osm->n_local,&osm->y);
351: PetscMalloc1(osm->n_local,&osm->y_local);
352: for (i=0; i<osm->n_local_true; ++i) {
353: ISGetLocalSize(osm->is[i],&m);
354: MatCreateVecs(osm->pmat[i],&osm->x[i],NULL);
355: ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
356: VecScatterCreate(vec,osm->is[i],osm->x[i],isl,&osm->restriction[i]);
357: ISDestroy(&isl);
358: VecDuplicate(osm->x[i],&osm->y[i]);
359: if (osm->localization) {
360: ISLocalToGlobalMapping ltog;
361: IS isll;
362: const PetscInt *idx_local;
363: PetscInt *idx,nout;
365: ISLocalToGlobalMappingCreateIS(osm->is[i],<og);
366: ISGetLocalSize(osm->is_local[i],&m_local);
367: ISGetIndices(osm->is_local[i], &idx_local);
368: PetscMalloc1(m_local,&idx);
369: ISGlobalToLocalMappingApply(ltog,IS_GTOLM_DROP,m_local,idx_local,&nout,idx);
370: ISLocalToGlobalMappingDestroy(<og);
371: if (nout != m_local) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"is_local not a subset of is");
372: ISRestoreIndices(osm->is_local[i], &idx_local);
373: ISCreateGeneral(PETSC_COMM_SELF,m_local,idx,PETSC_OWN_POINTER,&isll);
374: ISCreateStride(PETSC_COMM_SELF,m_local,0,1,&isl);
375: VecCreateSeq(PETSC_COMM_SELF,m_local,&osm->y_local[i]);
376: VecScatterCreate(osm->y[i],isll,osm->y_local[i],isl,&osm->localization[i]);
377: ISDestroy(&isll);
379: VecScatterCreate(vec,osm->is_local[i],osm->y_local[i],isl,&osm->prolongation[i]);
380: ISDestroy(&isl);
381: } else {
382: osm->y_local[i] = osm->y[i];
383: PetscObjectReference((PetscObject) osm->y[i]);
384: osm->prolongation[i] = osm->restriction[i];
385: PetscObjectReference((PetscObject) osm->restriction[i]);
386: }
387: }
388: for (i=osm->n_local_true; i<osm->n_local; i++) {
389: VecCreateSeq(PETSC_COMM_SELF,0,&osm->x[i]);
390: VecDuplicate(osm->x[i],&osm->y[i]);
391: VecDuplicate(osm->x[i],&osm->y_local[i]);
392: ISCreateStride(PETSC_COMM_SELF,0,0,1,&isl);
393: VecScatterCreate(vec,isl,osm->x[i],isl,&osm->restriction[i]);
394: if (osm->localization) {
395: VecScatterCreate(osm->y[i],isl,osm->y_local[i],isl,&osm->localization[i]);
396: VecScatterCreate(vec,isl,osm->x[i],isl,&osm->prolongation[i]);
397: } else {
398: osm->prolongation[i] = osm->restriction[i];
399: PetscObjectReference((PetscObject) osm->restriction[i]);
400: }
401: ISDestroy(&isl);
402: }
403: VecDestroy(&vec);
404: }
406: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
407: IS *cis;
408: PetscInt c;
410: PetscMalloc1(osm->n_local_true, &cis);
411: for (c = 0; c < osm->n_local_true; ++c) cis[c] = osm->lis;
412: MatCreateSubMatrices(pc->pmat, osm->n_local_true, osm->is, cis, scall, &osm->lmats);
413: PetscFree(cis);
414: }
416: /* Return control to the user so that the submatrices can be modified (e.g., to apply
417: different boundary conditions for the submatrices than for the global problem) */
418: PCModifySubMatrices(pc,osm->n_local_true,osm->is,osm->is,osm->pmat,pc->modifysubmatricesP);
420: /*
421: Loop over subdomains putting them into local ksp
422: */
423: for (i=0; i<osm->n_local_true; i++) {
424: KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i]);
425: if (!pc->setupcalled) {
426: KSPSetFromOptions(osm->ksp[i]);
427: }
428: }
429: return(0);
430: }
432: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
433: {
434: PC_ASM *osm = (PC_ASM*)pc->data;
435: PetscErrorCode ierr;
436: PetscInt i;
437: KSPConvergedReason reason;
440: for (i=0; i<osm->n_local_true; i++) {
441: KSPSetUp(osm->ksp[i]);
442: KSPGetConvergedReason(osm->ksp[i],&reason);
443: if (reason == KSP_DIVERGED_PCSETUP_FAILED) {
444: pc->failedreason = PC_SUBPC_ERROR;
445: }
446: }
447: return(0);
448: }
450: static PetscErrorCode PCApply_ASM(PC pc,Vec x,Vec y)
451: {
452: PC_ASM *osm = (PC_ASM*)pc->data;
454: PetscInt i,n_local = osm->n_local,n_local_true = osm->n_local_true;
455: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
458: /*
459: Support for limiting the restriction or interpolation to only local
460: subdomain values (leaving the other values 0).
461: */
462: if (!(osm->type & PC_ASM_RESTRICT)) {
463: forward = SCATTER_FORWARD_LOCAL;
464: /* have to zero the work RHS since scatter may leave some slots empty */
465: for (i=0; i<n_local_true; i++) {
466: VecZeroEntries(osm->x[i]);
467: }
468: }
469: if (!(osm->type & PC_ASM_INTERPOLATE)) reverse = SCATTER_REVERSE_LOCAL;
471: switch (osm->loctype)
472: {
473: case PC_COMPOSITE_ADDITIVE:
474: for (i=0; i<n_local; i++) {
475: VecScatterBegin(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
476: }
477: VecZeroEntries(y);
478: /* do the local solves */
479: for (i=0; i<n_local_true; i++) {
480: VecScatterEnd(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
481: KSPSolve(osm->ksp[i],osm->x[i],osm->y[i]);
482: if (osm->localization) {
483: VecScatterBegin(osm->localization[i],osm->y[i],osm->y_local[i],INSERT_VALUES,forward);
484: VecScatterEnd(osm->localization[i],osm->y[i],osm->y_local[i],INSERT_VALUES,forward);
485: }
486: VecScatterBegin(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
487: }
488: /* handle the rest of the scatters that do not have local solves */
489: for (i=n_local_true; i<n_local; i++) {
490: VecScatterEnd(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
491: VecScatterBegin(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
492: }
493: for (i=0; i<n_local; i++) {
494: VecScatterEnd(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
495: }
496: break;
497: case PC_COMPOSITE_MULTIPLICATIVE:
498: VecZeroEntries(y);
499: /* do the local solves */
500: for (i = 0; i < n_local_true; ++i) {
501: if (i > 0) {
502: /* Update rhs */
503: VecScatterBegin(osm->restriction[i], osm->lx, osm->x[i], INSERT_VALUES, forward);
504: VecScatterEnd(osm->restriction[i], osm->lx, osm->x[i], INSERT_VALUES, forward);
505: } else {
506: VecZeroEntries(osm->x[i]);
507: }
508: VecScatterBegin(osm->restriction[i], x, osm->x[i], ADD_VALUES, forward);
509: VecScatterEnd(osm->restriction[i], x, osm->x[i], ADD_VALUES, forward);
510: KSPSolve(osm->ksp[i], osm->x[i], osm->y[i]);
511: if (osm->localization) {
512: VecScatterBegin(osm->localization[i], osm->y[i], osm->y_local[i], INSERT_VALUES, forward);
513: VecScatterEnd(osm->localization[i], osm->y[i], osm->y_local[i], INSERT_VALUES, forward);
514: }
515: VecScatterBegin(osm->prolongation[i], osm->y_local[i], y, ADD_VALUES, reverse);
516: VecScatterEnd(osm->prolongation[i], osm->y_local[i], y, ADD_VALUES, reverse);
517: if (i < n_local_true-1) {
518: VecSet(osm->ly, 0.0);
519: VecScatterBegin(osm->prolongation[i], osm->y_local[i], osm->ly, INSERT_VALUES, reverse);
520: VecScatterEnd(osm->prolongation[i], osm->y_local[i], osm->ly, INSERT_VALUES, reverse);
521: VecScale(osm->ly, -1.0);
522: MatMult(osm->lmats[i+1], osm->ly, osm->y[i+1]);
523: VecScatterBegin(osm->restriction[i+1], osm->y[i+1], osm->lx, INSERT_VALUES, reverse);
524: VecScatterEnd(osm->restriction[i+1], osm->y[i+1], osm->lx, INSERT_VALUES, reverse);
525: }
526: }
527: /* handle the rest of the scatters that do not have local solves */
528: for (i = n_local_true; i < n_local; ++i) {
529: VecScatterBegin(osm->restriction[i], x, osm->x[i], INSERT_VALUES, forward);
530: VecScatterEnd(osm->restriction[i], x, osm->x[i], INSERT_VALUES, forward);
531: VecScatterBegin(osm->prolongation[i], osm->y_local[i], y, ADD_VALUES, reverse);
532: VecScatterEnd(osm->prolongation[i], osm->y_local[i], y, ADD_VALUES, reverse);
533: }
534: break;
535: default: SETERRQ1(PetscObjectComm((PetscObject) pc), PETSC_ERR_ARG_WRONG, "Invalid local composition type: %s", PCCompositeTypes[osm->loctype]);
536: }
537: return(0);
538: }
540: static PetscErrorCode PCApplyTranspose_ASM(PC pc,Vec x,Vec y)
541: {
542: PC_ASM *osm = (PC_ASM*)pc->data;
544: PetscInt i,n_local = osm->n_local,n_local_true = osm->n_local_true;
545: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
548: /*
549: Support for limiting the restriction or interpolation to only local
550: subdomain values (leaving the other values 0).
552: Note: these are reversed from the PCApply_ASM() because we are applying the
553: transpose of the three terms
554: */
555: if (!(osm->type & PC_ASM_INTERPOLATE)) {
556: forward = SCATTER_FORWARD_LOCAL;
557: /* have to zero the work RHS since scatter may leave some slots empty */
558: for (i=0; i<n_local_true; i++) {
559: VecZeroEntries(osm->x[i]);
560: }
561: }
562: if (!(osm->type & PC_ASM_RESTRICT)) reverse = SCATTER_REVERSE_LOCAL;
564: for (i=0; i<n_local; i++) {
565: VecScatterBegin(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
566: }
567: VecZeroEntries(y);
568: /* do the local solves */
569: for (i=0; i<n_local_true; i++) {
570: VecScatterEnd(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
571: KSPSolveTranspose(osm->ksp[i],osm->x[i],osm->y[i]);
572: if (osm->localization) {
573: VecScatterBegin(osm->localization[i],osm->y[i],osm->y_local[i],INSERT_VALUES,forward);
574: VecScatterEnd(osm->localization[i],osm->y[i],osm->y_local[i],INSERT_VALUES,forward);
575: }
576: VecScatterBegin(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
577: }
578: /* handle the rest of the scatters that do not have local solves */
579: for (i=n_local_true; i<n_local; i++) {
580: VecScatterEnd(osm->restriction[i],x,osm->x[i],INSERT_VALUES,forward);
581: VecScatterBegin(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
582: }
583: for (i=0; i<n_local; i++) {
584: VecScatterEnd(osm->prolongation[i],osm->y_local[i],y,ADD_VALUES,reverse);
585: }
586: return(0);
587: }
589: static PetscErrorCode PCReset_ASM(PC pc)
590: {
591: PC_ASM *osm = (PC_ASM*)pc->data;
593: PetscInt i;
596: if (osm->ksp) {
597: for (i=0; i<osm->n_local_true; i++) {
598: KSPReset(osm->ksp[i]);
599: }
600: }
601: if (osm->pmat) {
602: if (osm->n_local_true > 0) {
603: MatDestroySubMatrices(osm->n_local_true,&osm->pmat);
604: }
605: }
606: if (osm->restriction) {
607: for (i=0; i<osm->n_local; i++) {
608: VecScatterDestroy(&osm->restriction[i]);
609: if (osm->localization) {VecScatterDestroy(&osm->localization[i]);}
610: VecScatterDestroy(&osm->prolongation[i]);
611: VecDestroy(&osm->x[i]);
612: VecDestroy(&osm->y[i]);
613: VecDestroy(&osm->y_local[i]);
614: }
615: PetscFree(osm->restriction);
616: if (osm->localization) {PetscFree(osm->localization);}
617: PetscFree(osm->prolongation);
618: PetscFree(osm->x);
619: PetscFree(osm->y);
620: PetscFree(osm->y_local);
621: }
622: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
623: if (osm->loctype == PC_COMPOSITE_MULTIPLICATIVE) {
624: ISDestroy(&osm->lis);
625: MatDestroyMatrices(osm->n_local_true, &osm->lmats);
626: VecDestroy(&osm->lx);
627: VecDestroy(&osm->ly);
628: }
630: PetscFree(osm->sub_mat_type);
632: osm->is = 0;
633: osm->is_local = 0;
634: return(0);
635: }
637: static PetscErrorCode PCDestroy_ASM(PC pc)
638: {
639: PC_ASM *osm = (PC_ASM*)pc->data;
641: PetscInt i;
644: PCReset_ASM(pc);
645: if (osm->ksp) {
646: for (i=0; i<osm->n_local_true; i++) {
647: KSPDestroy(&osm->ksp[i]);
648: }
649: PetscFree(osm->ksp);
650: }
651: PetscFree(pc->data);
652: return(0);
653: }
655: static PetscErrorCode PCSetFromOptions_ASM(PetscOptionItems *PetscOptionsObject,PC pc)
656: {
657: PC_ASM *osm = (PC_ASM*)pc->data;
659: PetscInt blocks,ovl;
660: PetscBool symset,flg;
661: PCASMType asmtype;
662: PCCompositeType loctype;
663: char sub_mat_type[256];
666: /* set the type to symmetric if matrix is symmetric */
667: if (!osm->type_set && pc->pmat) {
668: MatIsSymmetricKnown(pc->pmat,&symset,&flg);
669: if (symset && flg) osm->type = PC_ASM_BASIC;
670: }
671: PetscOptionsHead(PetscOptionsObject,"Additive Schwarz options");
672: PetscOptionsBool("-pc_asm_dm_subdomains","Use DMCreateDomainDecomposition() to define subdomains","PCASMSetDMSubdomains",osm->dm_subdomains,&osm->dm_subdomains,&flg);
673: PetscOptionsInt("-pc_asm_blocks","Number of subdomains","PCASMSetTotalSubdomains",osm->n,&blocks,&flg);
674: if (flg) {
675: PCASMSetTotalSubdomains(pc,blocks,NULL,NULL);
676: osm->dm_subdomains = PETSC_FALSE;
677: }
678: PetscOptionsInt("-pc_asm_overlap","Number of grid points overlap","PCASMSetOverlap",osm->overlap,&ovl,&flg);
679: if (flg) {
680: PCASMSetOverlap(pc,ovl);
681: osm->dm_subdomains = PETSC_FALSE;
682: }
683: flg = PETSC_FALSE;
684: PetscOptionsEnum("-pc_asm_type","Type of restriction/extension","PCASMSetType",PCASMTypes,(PetscEnum)osm->type,(PetscEnum*)&asmtype,&flg);
685: if (flg) {PCASMSetType(pc,asmtype); }
686: flg = PETSC_FALSE;
687: PetscOptionsEnum("-pc_asm_local_type","Type of local solver composition","PCASMSetLocalType",PCCompositeTypes,(PetscEnum)osm->loctype,(PetscEnum*)&loctype,&flg);
688: if (flg) {PCASMSetLocalType(pc,loctype); }
689: PetscOptionsFList("-pc_asm_sub_mat_type","Subsolve Matrix Type","PCASMSetSubMatType",MatList,NULL,sub_mat_type,256,&flg);
690: if(flg){
691: PCASMSetSubMatType(pc,sub_mat_type);
692: }
693: PetscOptionsTail();
694: return(0);
695: }
697: /*------------------------------------------------------------------------------------*/
699: static PetscErrorCode PCASMSetLocalSubdomains_ASM(PC pc,PetscInt n,IS is[],IS is_local[])
700: {
701: PC_ASM *osm = (PC_ASM*)pc->data;
703: PetscInt i;
706: if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Each process must have 1 or more blocks, n = %D",n);
707: if (pc->setupcalled && (n != osm->n_local_true || is)) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetLocalSubdomains() should be called before calling PCSetUp().");
709: if (!pc->setupcalled) {
710: if (is) {
711: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is[i]);}
712: }
713: if (is_local) {
714: for (i=0; i<n; i++) {PetscObjectReference((PetscObject)is_local[i]);}
715: }
716: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
718: osm->n_local_true = n;
719: osm->is = 0;
720: osm->is_local = 0;
721: if (is) {
722: PetscMalloc1(n,&osm->is);
723: for (i=0; i<n; i++) osm->is[i] = is[i];
724: /* Flag indicating that the user has set overlapping subdomains so PCASM should not increase their size. */
725: osm->overlap = -1;
726: }
727: if (is_local) {
728: PetscMalloc1(n,&osm->is_local);
729: for (i=0; i<n; i++) osm->is_local[i] = is_local[i];
730: if (!is) {
731: PetscMalloc1(osm->n_local_true,&osm->is);
732: for (i=0; i<osm->n_local_true; i++) {
733: if (osm->overlap > 0) { /* With positive overlap, osm->is[i] will be modified */
734: ISDuplicate(osm->is_local[i],&osm->is[i]);
735: ISCopy(osm->is_local[i],osm->is[i]);
736: } else {
737: PetscObjectReference((PetscObject)osm->is_local[i]);
738: osm->is[i] = osm->is_local[i];
739: }
740: }
741: }
742: }
743: }
744: return(0);
745: }
747: static PetscErrorCode PCASMSetTotalSubdomains_ASM(PC pc,PetscInt N,IS *is,IS *is_local)
748: {
749: PC_ASM *osm = (PC_ASM*)pc->data;
751: PetscMPIInt rank,size;
752: PetscInt n;
755: if (N < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of total blocks must be > 0, N = %D",N);
756: 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.");
758: /*
759: Split the subdomains equally among all processors
760: */
761: MPI_Comm_rank(PetscObjectComm((PetscObject)pc),&rank);
762: MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
763: n = N/size + ((N % size) > rank);
764: 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);
765: if (pc->setupcalled && n != osm->n_local_true) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONGSTATE,"PCASMSetTotalSubdomains() should be called before PCSetUp().");
766: if (!pc->setupcalled) {
767: PCASMDestroySubdomains(osm->n_local_true,osm->is,osm->is_local);
769: osm->n_local_true = n;
770: osm->is = 0;
771: osm->is_local = 0;
772: }
773: return(0);
774: }
776: static PetscErrorCode PCASMSetOverlap_ASM(PC pc,PetscInt ovl)
777: {
778: PC_ASM *osm = (PC_ASM*)pc->data;
781: if (ovl < 0) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested");
782: if (pc->setupcalled && ovl != osm->overlap) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"PCASMSetOverlap() should be called before PCSetUp().");
783: if (!pc->setupcalled) osm->overlap = ovl;
784: return(0);
785: }
787: static PetscErrorCode PCASMSetType_ASM(PC pc,PCASMType type)
788: {
789: PC_ASM *osm = (PC_ASM*)pc->data;
792: osm->type = type;
793: osm->type_set = PETSC_TRUE;
794: return(0);
795: }
797: static PetscErrorCode PCASMGetType_ASM(PC pc,PCASMType *type)
798: {
799: PC_ASM *osm = (PC_ASM*)pc->data;
802: *type = osm->type;
803: return(0);
804: }
806: static PetscErrorCode PCASMSetLocalType_ASM(PC pc, PCCompositeType type)
807: {
808: PC_ASM *osm = (PC_ASM *) pc->data;
811: 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");
812: osm->loctype = type;
813: return(0);
814: }
816: static PetscErrorCode PCASMGetLocalType_ASM(PC pc, PCCompositeType *type)
817: {
818: PC_ASM *osm = (PC_ASM *) pc->data;
821: *type = osm->loctype;
822: return(0);
823: }
825: static PetscErrorCode PCASMSetSortIndices_ASM(PC pc,PetscBool doSort)
826: {
827: PC_ASM *osm = (PC_ASM*)pc->data;
830: osm->sort_indices = doSort;
831: return(0);
832: }
834: static PetscErrorCode PCASMGetSubKSP_ASM(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
835: {
836: PC_ASM *osm = (PC_ASM*)pc->data;
840: 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");
842: if (n_local) *n_local = osm->n_local_true;
843: if (first_local) {
844: MPI_Scan(&osm->n_local_true,first_local,1,MPIU_INT,MPI_SUM,PetscObjectComm((PetscObject)pc));
845: *first_local -= osm->n_local_true;
846: }
847: if (ksp) {
848: /* Assume that local solves are now different; not necessarily
849: true though! This flag is used only for PCView_ASM() */
850: *ksp = osm->ksp;
851: osm->same_local_solves = PETSC_FALSE;
852: }
853: return(0);
854: }
856: static PetscErrorCode PCASMGetSubMatType_ASM(PC pc,MatType *sub_mat_type)
857: {
858: PC_ASM *osm = (PC_ASM*)pc->data;
863: *sub_mat_type = osm->sub_mat_type;
864: return(0);
865: }
867: static PetscErrorCode PCASMSetSubMatType_ASM(PC pc,MatType sub_mat_type)
868: {
869: PetscErrorCode ierr;
870: PC_ASM *osm = (PC_ASM*)pc->data;
874: PetscFree(osm->sub_mat_type);
875: PetscStrallocpy(sub_mat_type,(char**)&osm->sub_mat_type);
876: return(0);
877: }
879: /*@C
880: PCASMSetLocalSubdomains - Sets the local subdomains (for this processor only) for the additive Schwarz preconditioner.
882: Collective on PC
884: Input Parameters:
885: + pc - the preconditioner context
886: . n - the number of subdomains for this processor (default value = 1)
887: . is - the index set that defines the subdomains for this processor
888: (or NULL for PETSc to determine subdomains)
889: - is_local - the index sets that define the local part of the subdomains for this processor, not used unless PCASMType is PC_ASM_RESTRICT
890: (or NULL to not provide these)
892: Notes:
893: The IS numbering is in the parallel, global numbering of the vector for both is and is_local
895: By default the ASM preconditioner uses 1 block per processor.
897: Use PCASMSetTotalSubdomains() to set the subdomains for all processors.
899: If is_local is provided and PCASMType is PC_ASM_RESTRICT then the solution only over the is_local region is interpolated
900: back to form the global solution (this is the standard restricted additive Schwarz method)
902: If the is_local is provided and PCASMType is PC_ASM_INTERPOLATE or PC_ASM_NONE then an error is generated since there is
903: no code to handle that case.
905: Level: advanced
907: .keywords: PC, ASM, set, local, subdomains, additive Schwarz
909: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
910: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), PCASMType, PCASMSetType()
911: @*/
912: PetscErrorCode PCASMSetLocalSubdomains(PC pc,PetscInt n,IS is[],IS is_local[])
913: {
918: PetscTryMethod(pc,"PCASMSetLocalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,n,is,is_local));
919: return(0);
920: }
922: /*@C
923: PCASMSetTotalSubdomains - Sets the subdomains for all processors for the
924: additive Schwarz preconditioner. Either all or no processors in the
925: PC communicator must call this routine, with the same index sets.
927: Collective on PC
929: Input Parameters:
930: + pc - the preconditioner context
931: . N - the number of subdomains for all processors
932: . is - the index sets that define the subdomains for all processors
933: (or NULL to ask PETSc to determine the subdomains)
934: - is_local - the index sets that define the local part of the subdomains for this processor
935: (or NULL to not provide this information)
937: Options Database Key:
938: To set the total number of subdomain blocks rather than specify the
939: index sets, use the option
940: . -pc_asm_blocks <blks> - Sets total blocks
942: Notes:
943: Currently you cannot use this to set the actual subdomains with the argument is or is_local.
945: By default the ASM preconditioner uses 1 block per processor.
947: These index sets cannot be destroyed until after completion of the
948: linear solves for which the ASM preconditioner is being used.
950: Use PCASMSetLocalSubdomains() to set local subdomains.
952: The IS numbering is in the parallel, global numbering of the vector for both is and is_local
954: Level: advanced
956: .keywords: PC, ASM, set, total, global, subdomains, additive Schwarz
958: .seealso: PCASMSetLocalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
959: PCASMCreateSubdomains2D()
960: @*/
961: PetscErrorCode PCASMSetTotalSubdomains(PC pc,PetscInt N,IS is[],IS is_local[])
962: {
967: PetscTryMethod(pc,"PCASMSetTotalSubdomains_C",(PC,PetscInt,IS[],IS[]),(pc,N,is,is_local));
968: return(0);
969: }
971: /*@
972: PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
973: additive Schwarz preconditioner. Either all or no processors in the
974: PC communicator must call this routine.
976: Logically Collective on PC
978: Input Parameters:
979: + pc - the preconditioner context
980: - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1)
982: Options Database Key:
983: . -pc_asm_overlap <ovl> - Sets overlap
985: Notes:
986: By default the ASM preconditioner uses 1 block per processor. To use
987: multiple blocks per perocessor, see PCASMSetTotalSubdomains() and
988: PCASMSetLocalSubdomains() (and the option -pc_asm_blocks <blks>).
990: The overlap defaults to 1, so if one desires that no additional
991: overlap be computed beyond what may have been set with a call to
992: PCASMSetTotalSubdomains() or PCASMSetLocalSubdomains(), then ovl
993: must be set to be 0. In particular, if one does not explicitly set
994: the subdomains an application code, then all overlap would be computed
995: internally by PETSc, and using an overlap of 0 would result in an ASM
996: variant that is equivalent to the block Jacobi preconditioner.
998: The default algorithm used by PETSc to increase overlap is fast, but not scalable,
999: use the option -mat_increase_overlap_scalable when the problem and number of processes is large.
1001: Note that one can define initial index sets with any overlap via
1002: PCASMSetLocalSubdomains(); the routine
1003: PCASMSetOverlap() merely allows PETSc to extend that overlap further
1004: if desired.
1006: Level: intermediate
1008: .keywords: PC, ASM, set, overlap
1010: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1011: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains(), MatIncreaseOverlap()
1012: @*/
1013: PetscErrorCode PCASMSetOverlap(PC pc,PetscInt ovl)
1014: {
1020: PetscTryMethod(pc,"PCASMSetOverlap_C",(PC,PetscInt),(pc,ovl));
1021: return(0);
1022: }
1024: /*@
1025: PCASMSetType - Sets the type of restriction and interpolation used
1026: for local problems in the additive Schwarz method.
1028: Logically Collective on PC
1030: Input Parameters:
1031: + pc - the preconditioner context
1032: - type - variant of ASM, one of
1033: .vb
1034: PC_ASM_BASIC - full interpolation and restriction
1035: PC_ASM_RESTRICT - full restriction, local processor interpolation
1036: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1037: PC_ASM_NONE - local processor restriction and interpolation
1038: .ve
1040: Options Database Key:
1041: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
1043: Notes: if the is_local arguments are passed to PCASMSetLocalSubdomains() then they are used when PC_ASM_RESTRICT has been selected
1044: to limit the local processor interpolation
1046: Level: intermediate
1048: .keywords: PC, ASM, set, type
1050: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1051: PCASMCreateSubdomains2D(), PCASMType, PCASMSetLocalType(), PCASMGetLocalType()
1052: @*/
1053: PetscErrorCode PCASMSetType(PC pc,PCASMType type)
1054: {
1060: PetscTryMethod(pc,"PCASMSetType_C",(PC,PCASMType),(pc,type));
1061: return(0);
1062: }
1064: /*@
1065: PCASMGetType - Gets the type of restriction and interpolation used
1066: for local problems in the additive Schwarz method.
1068: Logically Collective on PC
1070: Input Parameter:
1071: . pc - the preconditioner context
1073: Output Parameter:
1074: . type - variant of ASM, one of
1076: .vb
1077: PC_ASM_BASIC - full interpolation and restriction
1078: PC_ASM_RESTRICT - full restriction, local processor interpolation
1079: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
1080: PC_ASM_NONE - local processor restriction and interpolation
1081: .ve
1083: Options Database Key:
1084: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
1086: Level: intermediate
1088: .keywords: PC, ASM, set, type
1090: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1091: PCASMCreateSubdomains2D(), PCASMType, PCASMSetType(), PCASMSetLocalType(), PCASMGetLocalType()
1092: @*/
1093: PetscErrorCode PCASMGetType(PC pc,PCASMType *type)
1094: {
1099: PetscUseMethod(pc,"PCASMGetType_C",(PC,PCASMType*),(pc,type));
1100: return(0);
1101: }
1103: /*@
1104: PCASMSetLocalType - Sets the type of composition used for local problems in the additive Schwarz method.
1106: Logically Collective on PC
1108: Input Parameters:
1109: + pc - the preconditioner context
1110: - type - type of composition, one of
1111: .vb
1112: PC_COMPOSITE_ADDITIVE - local additive combination
1113: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1114: .ve
1116: Options Database Key:
1117: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1119: Level: intermediate
1121: .seealso: PCASMSetType(), PCASMGetType(), PCASMGetLocalType(), PCASM, PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1122: @*/
1123: PetscErrorCode PCASMSetLocalType(PC pc, PCCompositeType type)
1124: {
1130: PetscTryMethod(pc, "PCASMSetLocalType_C", (PC, PCCompositeType), (pc, type));
1131: return(0);
1132: }
1134: /*@
1135: PCASMGetLocalType - Gets the type of composition used for local problems in the additive Schwarz method.
1137: Logically Collective on PC
1139: Input Parameter:
1140: . pc - the preconditioner context
1142: Output Parameter:
1143: . type - type of composition, one of
1144: .vb
1145: PC_COMPOSITE_ADDITIVE - local additive combination
1146: PC_COMPOSITE_MULTIPLICATIVE - local multiplicative combination
1147: .ve
1149: Options Database Key:
1150: . -pc_asm_local_type [additive,multiplicative] - Sets local solver composition type
1152: Level: intermediate
1154: .seealso: PCASMSetType(), PCASMGetType(), PCASMSetLocalType(), PCASMCreate(), PCASMType, PCASMSetType(), PCASMGetType(), PCCompositeType
1155: @*/
1156: PetscErrorCode PCASMGetLocalType(PC pc, PCCompositeType *type)
1157: {
1163: PetscUseMethod(pc, "PCASMGetLocalType_C", (PC, PCCompositeType *), (pc, type));
1164: return(0);
1165: }
1167: /*@
1168: PCASMSetSortIndices - Determines whether subdomain indices are sorted.
1170: Logically Collective on PC
1172: Input Parameters:
1173: + pc - the preconditioner context
1174: - doSort - sort the subdomain indices
1176: Level: intermediate
1178: .keywords: PC, ASM, set, type
1180: .seealso: PCASMSetLocalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
1181: PCASMCreateSubdomains2D()
1182: @*/
1183: PetscErrorCode PCASMSetSortIndices(PC pc,PetscBool doSort)
1184: {
1190: PetscTryMethod(pc,"PCASMSetSortIndices_C",(PC,PetscBool),(pc,doSort));
1191: return(0);
1192: }
1194: /*@C
1195: PCASMGetSubKSP - Gets the local KSP contexts for all blocks on
1196: this processor.
1198: Collective on PC iff first_local is requested
1200: Input Parameter:
1201: . pc - the preconditioner context
1203: Output Parameters:
1204: + n_local - the number of blocks on this processor or NULL
1205: . first_local - the global number of the first block on this processor or NULL,
1206: all processors must request or all must pass NULL
1207: - ksp - the array of KSP contexts
1209: Note:
1210: After PCASMGetSubKSP() the array of KSPes is not to be freed.
1212: You must call KSPSetUp() before calling PCASMGetSubKSP().
1214: Fortran note:
1215: 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.
1217: Level: advanced
1219: .keywords: PC, ASM, additive Schwarz, get, sub, KSP, context
1221: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap(),
1222: PCASMCreateSubdomains2D(),
1223: @*/
1224: PetscErrorCode PCASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
1225: {
1230: PetscUseMethod(pc,"PCASMGetSubKSP_C",(PC,PetscInt*,PetscInt*,KSP **),(pc,n_local,first_local,ksp));
1231: return(0);
1232: }
1234: /* -------------------------------------------------------------------------------------*/
1235: /*MC
1236: PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
1237: its own KSP object.
1239: Options Database Keys:
1240: + -pc_asm_blocks <blks> - Sets total blocks
1241: . -pc_asm_overlap <ovl> - Sets overlap
1242: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type, default is restrict
1243: - -pc_asm_local_type [additive, multiplicative] - Sets ASM type, default is additive
1245: IMPORTANT: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you
1246: will get a different convergence rate due to the default option of -pc_asm_type restrict. Use
1247: -pc_asm_type basic to use the standard ASM.
1249: Notes: Each processor can have one or more blocks, but a block cannot be shared by more
1250: than one processor. Use PCGASM for subdomains shared by multiple processes. Defaults to one block per processor.
1252: To set options on the solvers for each block append -sub_ to all the KSP, and PC
1253: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
1255: To set the options on the solvers separate for each block call PCASMGetSubKSP()
1256: and set the options directly on the resulting KSP object (you can access its PC
1257: with KSPGetPC())
1259: Level: beginner
1261: Concepts: additive Schwarz method
1263: References:
1264: + 1. - M Dryja, OB Widlund, An additive variant of the Schwarz alternating method for the case of many subregions
1265: Courant Institute, New York University Technical report
1266: - 1. - Barry Smith, Petter Bjorstad, and William Gropp, Domain Decompositions: Parallel Multilevel Methods for Elliptic Partial Differential Equations,
1267: Cambridge University Press.
1269: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
1270: PCBJACOBI, PCASMGetSubKSP(), PCASMSetLocalSubdomains(), PCASMType, PCASMGetType(), PCASMSetLocalType(), PCASMGetLocalType()
1271: PCASMSetTotalSubdomains(), PCSetModifySubmatrices(), PCASMSetOverlap(), PCASMSetType(), PCCompositeType
1273: M*/
1275: PETSC_EXTERN PetscErrorCode PCCreate_ASM(PC pc)
1276: {
1278: PC_ASM *osm;
1281: PetscNewLog(pc,&osm);
1283: osm->n = PETSC_DECIDE;
1284: osm->n_local = 0;
1285: osm->n_local_true = PETSC_DECIDE;
1286: osm->overlap = 1;
1287: osm->ksp = 0;
1288: osm->restriction = 0;
1289: osm->localization = 0;
1290: osm->prolongation = 0;
1291: osm->x = 0;
1292: osm->y = 0;
1293: osm->y_local = 0;
1294: osm->is = 0;
1295: osm->is_local = 0;
1296: osm->mat = 0;
1297: osm->pmat = 0;
1298: osm->type = PC_ASM_RESTRICT;
1299: osm->loctype = PC_COMPOSITE_ADDITIVE;
1300: osm->same_local_solves = PETSC_TRUE;
1301: osm->sort_indices = PETSC_TRUE;
1302: osm->dm_subdomains = PETSC_FALSE;
1303: osm->sub_mat_type = NULL;
1305: pc->data = (void*)osm;
1306: pc->ops->apply = PCApply_ASM;
1307: pc->ops->applytranspose = PCApplyTranspose_ASM;
1308: pc->ops->setup = PCSetUp_ASM;
1309: pc->ops->reset = PCReset_ASM;
1310: pc->ops->destroy = PCDestroy_ASM;
1311: pc->ops->setfromoptions = PCSetFromOptions_ASM;
1312: pc->ops->setuponblocks = PCSetUpOnBlocks_ASM;
1313: pc->ops->view = PCView_ASM;
1314: pc->ops->applyrichardson = 0;
1316: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",PCASMSetLocalSubdomains_ASM);
1317: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",PCASMSetTotalSubdomains_ASM);
1318: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetOverlap_C",PCASMSetOverlap_ASM);
1319: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetType_C",PCASMSetType_ASM);
1320: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetType_C",PCASMGetType_ASM);
1321: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetLocalType_C",PCASMSetLocalType_ASM);
1322: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetLocalType_C",PCASMGetLocalType_ASM);
1323: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSortIndices_C",PCASMSetSortIndices_ASM);
1324: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubKSP_C",PCASMGetSubKSP_ASM);
1325: PetscObjectComposeFunction((PetscObject)pc,"PCASMGetSubMatType_C",PCASMGetSubMatType_ASM);
1326: PetscObjectComposeFunction((PetscObject)pc,"PCASMSetSubMatType_C",PCASMSetSubMatType_ASM);
1327: return(0);
1328: }
1330: /*@C
1331: PCASMCreateSubdomains - Creates the index sets for the overlapping Schwarz
1332: preconditioner for a any problem on a general grid.
1334: Collective
1336: Input Parameters:
1337: + A - The global matrix operator
1338: - n - the number of local blocks
1340: Output Parameters:
1341: . outis - the array of index sets defining the subdomains
1343: Level: advanced
1345: Note: this generates nonoverlapping subdomains; the PCASM will generate the overlap
1346: from these if you use PCASMSetLocalSubdomains()
1348: In the Fortran version you must provide the array outis[] already allocated of length n.
1350: .keywords: PC, ASM, additive Schwarz, create, subdomains, unstructured grid
1352: .seealso: PCASMSetLocalSubdomains(), PCASMDestroySubdomains()
1353: @*/
1354: PetscErrorCode PCASMCreateSubdomains(Mat A, PetscInt n, IS* outis[])
1355: {
1356: MatPartitioning mpart;
1357: const char *prefix;
1358: void (*f)(void);
1359: PetscInt i,j,rstart,rend,bs;
1360: PetscBool isbaij = PETSC_FALSE,foundpart = PETSC_FALSE;
1361: Mat Ad = NULL, adj;
1362: IS ispart,isnumb,*is;
1363: PetscErrorCode ierr;
1368: if (n < 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"number of local blocks must be > 0, n = %D",n);
1370: /* Get prefix, row distribution, and block size */
1371: MatGetOptionsPrefix(A,&prefix);
1372: MatGetOwnershipRange(A,&rstart,&rend);
1373: MatGetBlockSize(A,&bs);
1374: 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);
1376: /* Get diagonal block from matrix if possible */
1377: MatShellGetOperation(A,MATOP_GET_DIAGONAL_BLOCK,&f);
1378: if (f) {
1379: MatGetDiagonalBlock(A,&Ad);
1380: }
1381: if (Ad) {
1382: PetscObjectTypeCompare((PetscObject)Ad,MATSEQBAIJ,&isbaij);
1383: if (!isbaij) {PetscObjectTypeCompare((PetscObject)Ad,MATSEQSBAIJ,&isbaij);}
1384: }
1385: if (Ad && n > 1) {
1386: PetscBool match,done;
1387: /* Try to setup a good matrix partitioning if available */
1388: MatPartitioningCreate(PETSC_COMM_SELF,&mpart);
1389: PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);
1390: MatPartitioningSetFromOptions(mpart);
1391: PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGCURRENT,&match);
1392: if (!match) {
1393: PetscObjectTypeCompare((PetscObject)mpart,MATPARTITIONINGSQUARE,&match);
1394: }
1395: if (!match) { /* assume a "good" partitioner is available */
1396: PetscInt na;
1397: const PetscInt *ia,*ja;
1398: MatGetRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1399: if (done) {
1400: /* Build adjacency matrix by hand. Unfortunately a call to
1401: MatConvert(Ad,MATMPIADJ,MAT_INITIAL_MATRIX,&adj) will
1402: remove the block-aij structure and we cannot expect
1403: MatPartitioning to split vertices as we need */
1404: PetscInt i,j,len,nnz,cnt,*iia=0,*jja=0;
1405: const PetscInt *row;
1406: nnz = 0;
1407: for (i=0; i<na; i++) { /* count number of nonzeros */
1408: len = ia[i+1] - ia[i];
1409: row = ja + ia[i];
1410: for (j=0; j<len; j++) {
1411: if (row[j] == i) { /* don't count diagonal */
1412: len--; break;
1413: }
1414: }
1415: nnz += len;
1416: }
1417: PetscMalloc1(na+1,&iia);
1418: PetscMalloc1(nnz,&jja);
1419: nnz = 0;
1420: iia[0] = 0;
1421: for (i=0; i<na; i++) { /* fill adjacency */
1422: cnt = 0;
1423: len = ia[i+1] - ia[i];
1424: row = ja + ia[i];
1425: for (j=0; j<len; j++) {
1426: if (row[j] != i) { /* if not diagonal */
1427: jja[nnz+cnt++] = row[j];
1428: }
1429: }
1430: nnz += cnt;
1431: iia[i+1] = nnz;
1432: }
1433: /* Partitioning of the adjacency matrix */
1434: MatCreateMPIAdj(PETSC_COMM_SELF,na,na,iia,jja,NULL,&adj);
1435: MatPartitioningSetAdjacency(mpart,adj);
1436: MatPartitioningSetNParts(mpart,n);
1437: MatPartitioningApply(mpart,&ispart);
1438: ISPartitioningToNumbering(ispart,&isnumb);
1439: MatDestroy(&adj);
1440: foundpart = PETSC_TRUE;
1441: }
1442: MatRestoreRowIJ(Ad,0,PETSC_TRUE,isbaij,&na,&ia,&ja,&done);
1443: }
1444: MatPartitioningDestroy(&mpart);
1445: }
1447: PetscMalloc1(n,&is);
1448: *outis = is;
1450: if (!foundpart) {
1452: /* Partitioning by contiguous chunks of rows */
1454: PetscInt mbs = (rend-rstart)/bs;
1455: PetscInt start = rstart;
1456: for (i=0; i<n; i++) {
1457: PetscInt count = (mbs/n + ((mbs % n) > i)) * bs;
1458: ISCreateStride(PETSC_COMM_SELF,count,start,1,&is[i]);
1459: start += count;
1460: }
1462: } else {
1464: /* Partitioning by adjacency of diagonal block */
1466: const PetscInt *numbering;
1467: PetscInt *count,nidx,*indices,*newidx,start=0;
1468: /* Get node count in each partition */
1469: PetscMalloc1(n,&count);
1470: ISPartitioningCount(ispart,n,count);
1471: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1472: for (i=0; i<n; i++) count[i] *= bs;
1473: }
1474: /* Build indices from node numbering */
1475: ISGetLocalSize(isnumb,&nidx);
1476: PetscMalloc1(nidx,&indices);
1477: for (i=0; i<nidx; i++) indices[i] = i; /* needs to be initialized */
1478: ISGetIndices(isnumb,&numbering);
1479: PetscSortIntWithPermutation(nidx,numbering,indices);
1480: ISRestoreIndices(isnumb,&numbering);
1481: if (isbaij && bs > 1) { /* adjust for the block-aij case */
1482: PetscMalloc1(nidx*bs,&newidx);
1483: for (i=0; i<nidx; i++) {
1484: for (j=0; j<bs; j++) newidx[i*bs+j] = indices[i]*bs + j;
1485: }
1486: PetscFree(indices);
1487: nidx *= bs;
1488: indices = newidx;
1489: }
1490: /* Shift to get global indices */
1491: for (i=0; i<nidx; i++) indices[i] += rstart;
1493: /* Build the index sets for each block */
1494: for (i=0; i<n; i++) {
1495: ISCreateGeneral(PETSC_COMM_SELF,count[i],&indices[start],PETSC_COPY_VALUES,&is[i]);
1496: ISSort(is[i]);
1497: start += count[i];
1498: }
1500: PetscFree(count);
1501: PetscFree(indices);
1502: ISDestroy(&isnumb);
1503: ISDestroy(&ispart);
1505: }
1506: return(0);
1507: }
1509: /*@C
1510: PCASMDestroySubdomains - Destroys the index sets created with
1511: PCASMCreateSubdomains(). Should be called after setting subdomains
1512: with PCASMSetLocalSubdomains().
1514: Collective
1516: Input Parameters:
1517: + n - the number of index sets
1518: . is - the array of index sets
1519: - is_local - the array of local index sets, can be NULL
1521: Level: advanced
1523: .keywords: PC, ASM, additive Schwarz, create, subdomains, unstructured grid
1525: .seealso: PCASMCreateSubdomains(), PCASMSetLocalSubdomains()
1526: @*/
1527: PetscErrorCode PCASMDestroySubdomains(PetscInt n, IS is[], IS is_local[])
1528: {
1529: PetscInt i;
1533: if (n <= 0) return(0);
1534: if (is) {
1536: for (i=0; i<n; i++) { ISDestroy(&is[i]); }
1537: PetscFree(is);
1538: }
1539: if (is_local) {
1541: for (i=0; i<n; i++) { ISDestroy(&is_local[i]); }
1542: PetscFree(is_local);
1543: }
1544: return(0);
1545: }
1547: /*@
1548: PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
1549: preconditioner for a two-dimensional problem on a regular grid.
1551: Not Collective
1553: Input Parameters:
1554: + m, n - the number of mesh points in the x and y directions
1555: . M, N - the number of subdomains in the x and y directions
1556: . dof - degrees of freedom per node
1557: - overlap - overlap in mesh lines
1559: Output Parameters:
1560: + Nsub - the number of subdomains created
1561: . is - array of index sets defining overlapping (if overlap > 0) subdomains
1562: - is_local - array of index sets defining non-overlapping subdomains
1564: Note:
1565: Presently PCAMSCreateSubdomains2d() is valid only for sequential
1566: preconditioners. More general related routines are
1567: PCASMSetTotalSubdomains() and PCASMSetLocalSubdomains().
1569: Level: advanced
1571: .keywords: PC, ASM, additive Schwarz, create, subdomains, 2D, regular grid
1573: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
1574: PCASMSetOverlap()
1575: @*/
1576: PetscErrorCode PCASMCreateSubdomains2D(PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt dof,PetscInt overlap,PetscInt *Nsub,IS **is,IS **is_local)
1577: {
1578: PetscInt i,j,height,width,ystart,xstart,yleft,yright,xleft,xright,loc_outer;
1580: PetscInt nidx,*idx,loc,ii,jj,count;
1583: if (dof != 1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP," ");
1585: *Nsub = N*M;
1586: PetscMalloc1(*Nsub,is);
1587: PetscMalloc1(*Nsub,is_local);
1588: ystart = 0;
1589: loc_outer = 0;
1590: for (i=0; i<N; i++) {
1591: height = n/N + ((n % N) > i); /* height of subdomain */
1592: if (height < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many N subdomains for mesh dimension n");
1593: yleft = ystart - overlap; if (yleft < 0) yleft = 0;
1594: yright = ystart + height + overlap; if (yright > n) yright = n;
1595: xstart = 0;
1596: for (j=0; j<M; j++) {
1597: width = m/M + ((m % M) > j); /* width of subdomain */
1598: if (width < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Too many M subdomains for mesh dimension m");
1599: xleft = xstart - overlap; if (xleft < 0) xleft = 0;
1600: xright = xstart + width + overlap; if (xright > m) xright = m;
1601: nidx = (xright - xleft)*(yright - yleft);
1602: PetscMalloc1(nidx,&idx);
1603: loc = 0;
1604: for (ii=yleft; ii<yright; ii++) {
1605: count = m*ii + xleft;
1606: for (jj=xleft; jj<xright; jj++) idx[loc++] = count++;
1607: }
1608: ISCreateGeneral(PETSC_COMM_SELF,nidx,idx,PETSC_COPY_VALUES,(*is)+loc_outer);
1609: if (overlap == 0) {
1610: PetscObjectReference((PetscObject)(*is)[loc_outer]);
1612: (*is_local)[loc_outer] = (*is)[loc_outer];
1613: } else {
1614: for (loc=0,ii=ystart; ii<ystart+height; ii++) {
1615: for (jj=xstart; jj<xstart+width; jj++) {
1616: idx[loc++] = m*ii + jj;
1617: }
1618: }
1619: ISCreateGeneral(PETSC_COMM_SELF,loc,idx,PETSC_COPY_VALUES,*is_local+loc_outer);
1620: }
1621: PetscFree(idx);
1622: xstart += width;
1623: loc_outer++;
1624: }
1625: ystart += height;
1626: }
1627: for (i=0; i<*Nsub; i++) { ISSort((*is)[i]); }
1628: return(0);
1629: }
1631: /*@C
1632: PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
1633: only) for the additive Schwarz preconditioner.
1635: Not Collective
1637: Input Parameter:
1638: . pc - the preconditioner context
1640: Output Parameters:
1641: + n - the number of subdomains for this processor (default value = 1)
1642: . is - the index sets that define the subdomains for this processor
1643: - is_local - the index sets that define the local part of the subdomains for this processor (can be NULL)
1646: Notes:
1647: The IS numbering is in the parallel, global numbering of the vector.
1649: Level: advanced
1651: .keywords: PC, ASM, set, local, subdomains, additive Schwarz
1653: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1654: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubmatrices()
1655: @*/
1656: PetscErrorCode PCASMGetLocalSubdomains(PC pc,PetscInt *n,IS *is[],IS *is_local[])
1657: {
1658: PC_ASM *osm = (PC_ASM*)pc->data;
1660: PetscBool match;
1666: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1667: if (!match) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"PC is not a PCASM");
1668: if (n) *n = osm->n_local_true;
1669: if (is) *is = osm->is;
1670: if (is_local) *is_local = osm->is_local;
1671: return(0);
1672: }
1674: /*@C
1675: PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1676: only) for the additive Schwarz preconditioner.
1678: Not Collective
1680: Input Parameter:
1681: . pc - the preconditioner context
1683: Output Parameters:
1684: + n - the number of matrices for this processor (default value = 1)
1685: - mat - the matrices
1687: Level: advanced
1689: Notes: Call after PCSetUp() (or KSPSetUp()) but before PCApply() (or KSPApply()) and before PCSetUpOnBlocks())
1691: Usually one would use PCSetModifySubmatrices() to change the submatrices in building the preconditioner.
1693: .keywords: PC, ASM, set, local, subdomains, additive Schwarz, block Jacobi
1695: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1696: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains(), PCSetModifySubmatrices()
1697: @*/
1698: PetscErrorCode PCASMGetLocalSubmatrices(PC pc,PetscInt *n,Mat *mat[])
1699: {
1700: PC_ASM *osm;
1702: PetscBool match;
1708: if (!pc->setupcalled) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUP() or PCSetUp().");
1709: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1710: if (!match) {
1711: if (n) *n = 0;
1712: if (mat) *mat = NULL;
1713: } else {
1714: osm = (PC_ASM*)pc->data;
1715: if (n) *n = osm->n_local_true;
1716: if (mat) *mat = osm->pmat;
1717: }
1718: return(0);
1719: }
1721: /*@
1722: PCASMSetDMSubdomains - Indicates whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1724: Logically Collective
1726: Input Parameter:
1727: + pc - the preconditioner
1728: - flg - boolean indicating whether to use subdomains defined by the DM
1730: Options Database Key:
1731: . -pc_asm_dm_subdomains
1733: Level: intermediate
1735: Notes:
1736: PCASMSetTotalSubdomains() and PCASMSetOverlap() take precedence over PCASMSetDMSubdomains(),
1737: so setting either of the first two effectively turns the latter off.
1739: .keywords: PC, ASM, DM, set, subdomains, additive Schwarz
1741: .seealso: PCASMGetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1742: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1743: @*/
1744: PetscErrorCode PCASMSetDMSubdomains(PC pc,PetscBool flg)
1745: {
1746: PC_ASM *osm = (PC_ASM*)pc->data;
1748: PetscBool match;
1753: if (pc->setupcalled) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Not for a setup PC.");
1754: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1755: if (match) {
1756: osm->dm_subdomains = flg;
1757: }
1758: return(0);
1759: }
1761: /*@
1762: PCASMGetDMSubdomains - Returns flag indicating whether to use DMCreateDomainDecomposition() to define the subdomains, whenever possible.
1763: Not Collective
1765: Input Parameter:
1766: . pc - the preconditioner
1768: Output Parameter:
1769: . flg - boolean indicating whether to use subdomains defined by the DM
1771: Level: intermediate
1773: .keywords: PC, ASM, DM, set, subdomains, additive Schwarz
1775: .seealso: PCASMSetDMSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap()
1776: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1777: @*/
1778: PetscErrorCode PCASMGetDMSubdomains(PC pc,PetscBool* flg)
1779: {
1780: PC_ASM *osm = (PC_ASM*)pc->data;
1782: PetscBool match;
1787: PetscObjectTypeCompare((PetscObject)pc,PCASM,&match);
1788: if (match) {
1789: if (flg) *flg = osm->dm_subdomains;
1790: }
1791: return(0);
1792: }
1794: /*@
1795: PCASMGetSubMatType - Gets the matrix type used for ASM subsolves, as a string.
1797: Not Collective
1799: Input Parameter:
1800: . pc - the PC
1802: Output Parameter:
1803: . -pc_asm_sub_mat_type - name of matrix type
1805: Level: advanced
1807: .keywords: PC, PCASM, MatType, set
1809: .seealso: PCASMSetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1810: @*/
1811: PetscErrorCode PCASMGetSubMatType(PC pc,MatType *sub_mat_type){
1814: PetscTryMethod(pc,"PCASMGetSubMatType_C",(PC,MatType*),(pc,sub_mat_type));
1815: return(0);
1816: }
1818: /*@
1819: PCASMSetSubMatType - Set the type of matrix used for ASM subsolves
1821: Collective on Mat
1823: Input Parameters:
1824: + pc - the PC object
1825: - sub_mat_type - matrix type
1827: Options Database Key:
1828: . -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.
1830: Notes:
1831: See "${PETSC_DIR}/include/petscmat.h" for available types
1833: Level: advanced
1835: .keywords: PC, PCASM, MatType, set
1837: .seealso: PCASMGetSubMatType(), PCASM, PCSetType(), VecSetType(), MatType, Mat
1838: @*/
1839: PetscErrorCode PCASMSetSubMatType(PC pc,MatType sub_mat_type)
1840: {
1843: PetscTryMethod(pc,"PCASMSetSubMatType_C",(PC,MatType),(pc,sub_mat_type));
1844: return(0);
1845: }