Actual source code: asm.c
1: /*$Id: asm.c,v 1.131 2001/08/07 03:03:38 balay Exp $*/
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 src/sles/pc/pcimpl.h
14: #include petscsles.h
16: typedef struct {
17: int n,n_local,n_local_true;
18: PetscTruth is_flg; /* flg set to 1 if the IS created in pcsetup */
19: int overlap; /* overlap requested by user */
20: SLES *sles; /* linear solvers for each block */
21: VecScatter *scat; /* mapping to subregion */
22: Vec *x,*y;
23: IS *is; /* index set that defines each subdomain */
24: Mat *mat,*pmat; /* mat is not currently used */
25: PCASMType type; /* use reduced interpolation, restriction or both */
26: PetscTruth same_local_solves; /* flag indicating whether all local solvers are same */
27: PetscTruth inplace; /* indicates that the sub-matrices are deleted after
28: PCSetUpOnBlocks() is done. Similar to inplace
29: factorization in the case of LU and ILU */
30: } PC_ASM;
32: #undef __FUNCT__
34: static int PCView_ASM(PC pc,PetscViewer viewer)
35: {
36: PC_ASM *jac = (PC_ASM*)pc->data;
37: int rank,ierr,i;
38: char *cstring = 0;
39: PetscTruth isascii,isstring;
40: PetscViewer sviewer;
44: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
45: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
46: if (isascii) {
47: PetscViewerASCIIPrintf(viewer," Additive Schwarz: total subdomain blocks = %d, amount of overlap = %dn",jac->n,jac->overlap);
48: if (jac->type == PC_ASM_NONE) cstring = "limited restriction and interpolation (PC_ASM_NONE)";
49: else if (jac->type == PC_ASM_RESTRICT) cstring = "full restriction (PC_ASM_RESTRICT)";
50: else if (jac->type == PC_ASM_INTERPOLATE) cstring = "full interpolation (PC_ASM_INTERPOLATE)";
51: else if (jac->type == PC_ASM_BASIC) cstring = "full restriction and interpolation (PC_ASM_BASIC)";
52: else cstring = "Unknown ASM type";
53: PetscViewerASCIIPrintf(viewer," Additive Schwarz: type - %sn",cstring);
54: MPI_Comm_rank(pc->comm,&rank);
55: if (jac->same_local_solves) {
56: PetscViewerASCIIPrintf(viewer," Local solve is same for all blocks, in the following KSP and PC objects:n");
57: PetscViewerGetSingleton(viewer,&sviewer);
58: if (!rank && jac->sles) {
59: PetscViewerASCIIPushTab(viewer);
60: SLESView(jac->sles[0],sviewer);
61: PetscViewerASCIIPopTab(viewer);
62: }
63: PetscViewerRestoreSingleton(viewer,&sviewer);
64: } else {
65: PetscViewerASCIIPrintf(viewer," Local solve info for each block is in the following KSP and PC objects:n");
66: PetscViewerASCIISynchronizedPrintf(viewer,"Proc %d: number of local blocks = %dn",rank,jac->n_local);
67: PetscViewerASCIIPushTab(viewer);
68: for (i=0; i<jac->n_local; i++) {
69: PetscViewerASCIISynchronizedPrintf(viewer,"Proc %d: local block number %dn",rank,i);
70: PetscViewerGetSingleton(viewer,&sviewer);
71: SLESView(jac->sles[i],sviewer);
72: PetscViewerRestoreSingleton(viewer,&sviewer);
73: if (i != jac->n_local-1) {
74: PetscViewerASCIISynchronizedPrintf(viewer,"- - - - - - - - - - - - - - - - - -n");
75: }
76: }
77: PetscViewerASCIIPopTab(viewer);
78: PetscViewerFlush(viewer);
79: }
80: } else if (isstring) {
81: PetscViewerStringSPrintf(viewer," blks=%d, overlap=%d, type=%d",jac->n,jac->overlap,jac->type);
82: PetscViewerGetSingleton(viewer,&sviewer);
83: if (jac->sles) {SLESView(jac->sles[0],sviewer);}
84: PetscViewerGetSingleton(viewer,&sviewer);
85: } else {
86: SETERRQ1(1,"Viewer type %s not supported for PCASM",((PetscObject)viewer)->type_name);
87: }
88: return(0);
89: }
91: #undef __FUNCT__
93: static int PCSetUp_ASM(PC pc)
94: {
95: PC_ASM *osm = (PC_ASM*)pc->data;
96: int i,ierr,m,n_local = osm->n_local,n_local_true = osm->n_local_true;
97: int start,start_val,end_val,size,sz,bs;
98: MatReuse scall = MAT_REUSE_MATRIX;
99: IS isl;
100: SLES sles;
101: KSP subksp;
102: PC subpc;
103: char *prefix,*pprefix;
106: if (!pc->setupcalled) {
107: if (osm->n == PETSC_DECIDE && osm->n_local_true == PETSC_DECIDE) {
108: /* no subdomains given, use one per processor */
109: osm->n_local_true = osm->n_local = 1;
110: MPI_Comm_size(pc->comm,&size);
111: osm->n = size;
112: } else if (osm->n == PETSC_DECIDE) { /* determine global number of subdomains */
113: int inwork[2],outwork[2];
114: inwork[0] = inwork[1] = osm->n_local_true;
115: MPI_Allreduce(inwork,outwork,1,MPI_2INT,PetscMaxSum_Op,pc->comm);
116: osm->n_local = outwork[0];
117: osm->n = outwork[1];
118: }
119: n_local = osm->n_local;
120: n_local_true = osm->n_local_true;
121: if (!osm->is){ /* build the index sets */
122: ierr = PetscMalloc((n_local_true+1)*sizeof(IS **),&osm->is);
123: ierr = MatGetOwnershipRange(pc->pmat,&start_val,&end_val);
124: ierr = MatGetBlockSize(pc->pmat,&bs);
125: sz = end_val - start_val;
126: start = start_val;
127: if (end_val/bs*bs != end_val || start_val/bs*bs != start_val) {
128: SETERRQ(PETSC_ERR_ARG_WRONG,"Bad distribution for matrix block size");
129: }
130: for (i=0; i<n_local_true; i++){
131: size = ((sz/bs)/n_local_true + (((sz/bs) % n_local_true) > i))*bs;
132: ierr = ISCreateStride(PETSC_COMM_SELF,size,start,1,&isl);
133: start += size;
134: osm->is[i] = isl;
135: }
136: osm->is_flg = PETSC_TRUE;
137: }
139: ierr = PetscMalloc((n_local_true+1)*sizeof(SLES **),&osm->sles);
140: ierr = PetscMalloc(n_local*sizeof(VecScatter **),&osm->scat);
141: ierr = PetscMalloc(2*n_local*sizeof(Vec **),&osm->x);
142: osm->y = osm->x + n_local;
144: /* Extend the "overlapping" regions by a number of steps */
145: MatIncreaseOverlap(pc->pmat,n_local_true,osm->is,osm->overlap);
146: for (i=0; i<n_local_true; i++) {
147: ISSort(osm->is[i]);
148: }
150: /* create the local work vectors and scatter contexts */
151: for (i=0; i<n_local_true; i++) {
152: ISGetLocalSize(osm->is[i],&m);
153: VecCreateSeq(PETSC_COMM_SELF,m,&osm->x[i]);
154: VecDuplicate(osm->x[i],&osm->y[i]);
155: ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
156: VecScatterCreate(pc->vec,osm->is[i],osm->x[i],isl,&osm->scat[i]);
157: ISDestroy(isl);
158: }
159: for (i=n_local_true; i<n_local; i++) {
160: VecCreateSeq(PETSC_COMM_SELF,0,&osm->x[i]);
161: VecDuplicate(osm->x[i],&osm->y[i]);
162: ISCreateStride(PETSC_COMM_SELF,0,0,1,&isl);
163: VecScatterCreate(pc->vec,isl,osm->x[i],isl,&osm->scat[i]);
164: ISDestroy(isl);
165: }
167: /*
168: Create the local solvers.
169: */
170: for (i=0; i<n_local_true; i++) {
171: SLESCreate(PETSC_COMM_SELF,&sles);
172: PetscLogObjectParent(pc,sles);
173: SLESGetKSP(sles,&subksp);
174: KSPSetType(subksp,KSPPREONLY);
175: SLESGetPC(sles,&subpc);
176: PCSetType(subpc,PCILU);
177: PCGetOptionsPrefix(pc,&prefix);
178: SLESSetOptionsPrefix(sles,prefix);
179: SLESAppendOptionsPrefix(sles,"sub_");
180: SLESSetFromOptions(sles);
181: osm->sles[i] = sles;
182: }
183: scall = MAT_INITIAL_MATRIX;
184: } else {
185: /*
186: Destroy the blocks from the previous iteration
187: */
188: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
189: MatDestroyMatrices(osm->n_local_true,&osm->pmat);
190: scall = MAT_INITIAL_MATRIX;
191: }
192: }
194: /* extract out the submatrices */
195: MatGetSubMatrices(pc->pmat,osm->n_local_true,osm->is,osm->is,scall,&osm->pmat);
197: /* Return control to the user so that the submatrices can be modified (e.g., to apply
198: different boundary conditions for the submatrices than for the global problem) */
199: PCModifySubMatrices(pc,osm->n_local,osm->is,osm->is,osm->pmat,pc->modifysubmatricesP);
201: /* loop over subdomains putting them into local sles */
202: PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
203: for (i=0; i<n_local_true; i++) {
204: PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);
205: PetscLogObjectParent(pc,osm->pmat[i]);
206: SLESSetOperators(osm->sles[i],osm->pmat[i],osm->pmat[i],pc->flag);
207: }
208: return(0);
209: }
211: #undef __FUNCT__
213: static int PCSetUpOnBlocks_ASM(PC pc)
214: {
215: PC_ASM *osm = (PC_ASM*)pc->data;
216: int i,ierr;
219: for (i=0; i<osm->n_local_true; i++) {
220: SLESSetUp(osm->sles[i],osm->x[i],osm->y[i]);
221: }
222: /*
223: If inplace flag is set, then destroy the matrix after the setup
224: on blocks is done.
225: */
226: if (osm->inplace && osm->n_local_true > 0) {
227: MatDestroyMatrices(osm->n_local_true,&osm->pmat);
228: }
229: return(0);
230: }
232: #undef __FUNCT__
234: static int PCApply_ASM(PC pc,Vec x,Vec y)
235: {
236: PC_ASM *osm = (PC_ASM*)pc->data;
237: int i,n_local = osm->n_local,n_local_true = osm->n_local_true,ierr,its;
238: PetscScalar zero = 0.0;
239: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
242: /*
243: Support for limiting the restriction or interpolation to only local
244: subdomain values (leaving the other values 0).
245: */
246: if (!(osm->type & PC_ASM_RESTRICT)) {
247: forward = SCATTER_FORWARD_LOCAL;
248: /* have to zero the work RHS since scatter may leave some slots empty */
249: for (i=0; i<n_local; i++) {
250: VecSet(&zero,osm->x[i]);
251: }
252: }
253: if (!(osm->type & PC_ASM_INTERPOLATE)) {
254: reverse = SCATTER_REVERSE_LOCAL;
255: }
257: for (i=0; i<n_local; i++) {
258: VecScatterBegin(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
259: }
260: VecSet(&zero,y);
261: /* do the local solves */
262: for (i=0; i<n_local_true; i++) {
263: VecScatterEnd(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
264: SLESSolve(osm->sles[i],osm->x[i],osm->y[i],&its);
265: VecScatterBegin(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
266: }
267: /* handle the rest of the scatters that do not have local solves */
268: for (i=n_local_true; i<n_local; i++) {
269: VecScatterEnd(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
270: VecScatterBegin(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
271: }
272: for (i=0; i<n_local; i++) {
273: VecScatterEnd(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
274: }
275: return(0);
276: }
278: #undef __FUNCT__
280: static int PCApplyTranspose_ASM(PC pc,Vec x,Vec y)
281: {
282: PC_ASM *osm = (PC_ASM*)pc->data;
283: int i,n_local = osm->n_local,n_local_true = osm->n_local_true,ierr,its;
284: PetscScalar zero = 0.0;
285: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
288: /*
289: Support for limiting the restriction or interpolation to only local
290: subdomain values (leaving the other values 0).
292: Note: these are reversed from the PCApply_ASM() because we are applying the
293: transpose of the three terms
294: */
295: if (!(osm->type & PC_ASM_INTERPOLATE)) {
296: forward = SCATTER_FORWARD_LOCAL;
297: /* have to zero the work RHS since scatter may leave some slots empty */
298: for (i=0; i<n_local; i++) {
299: VecSet(&zero,osm->x[i]);
300: }
301: }
302: if (!(osm->type & PC_ASM_RESTRICT)) {
303: reverse = SCATTER_REVERSE_LOCAL;
304: }
306: for (i=0; i<n_local; i++) {
307: VecScatterBegin(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
308: }
309: VecSet(&zero,y);
310: /* do the local solves */
311: for (i=0; i<n_local_true; i++) {
312: VecScatterEnd(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
313: SLESSolveTranspose(osm->sles[i],osm->x[i],osm->y[i],&its);
314: VecScatterBegin(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
315: }
316: /* handle the rest of the scatters that do not have local solves */
317: for (i=n_local_true; i<n_local; i++) {
318: VecScatterEnd(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
319: VecScatterBegin(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
320: }
321: for (i=0; i<n_local; i++) {
322: VecScatterEnd(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
323: }
324: return(0);
325: }
327: #undef __FUNCT__
329: static int PCDestroy_ASM(PC pc)
330: {
331: PC_ASM *osm = (PC_ASM*)pc->data;
332: int i,ierr;
335: for (i=0; i<osm->n_local; i++) {
336: VecScatterDestroy(osm->scat[i]);
337: VecDestroy(osm->x[i]);
338: VecDestroy(osm->y[i]);
339: }
340: if (osm->n_local_true > 0 && !osm->inplace && osm->pmat) {
341: MatDestroyMatrices(osm->n_local_true,&osm->pmat);
342: }
343: if (osm->sles) {
344: for (i=0; i<osm->n_local_true; i++) {
345: SLESDestroy(osm->sles[i]);
346: }
347: }
348: if (osm->is_flg) {
349: for (i=0; i<osm->n_local_true; i++) {ISDestroy(osm->is[i]);}
350: PetscFree(osm->is);
351: }
352: if (osm->sles) {PetscFree(osm->sles);}
353: if (osm->scat) {PetscFree(osm->scat);}
354: if (osm->x) {PetscFree(osm->x);}
355: PetscFree(osm);
356: return(0);
357: }
359: #undef __FUNCT__
361: static int PCSetFromOptions_ASM(PC pc)
362: {
363: PC_ASM *osm = (PC_ASM*)pc->data;
364: int blocks,ovl,ierr;
365: PetscTruth flg;
366: char buff[16],*type[] = {"basic","restrict","interpolate","none"};
369: PetscOptionsHead("Additive Schwarz options");
370: PetscOptionsInt("-pc_asm_blocks","Number of subdomains","PCASMSetTotalSubdomains",osm->n,&blocks,&flg);
371: if (flg) {PCASMSetTotalSubdomains(pc,blocks,PETSC_NULL); }
372: PetscOptionsInt("-pc_asm_overlap","Number of grid points overlap","PCASMSetOverlap",osm->overlap,&ovl,&flg);
373: if (flg) {PCASMSetOverlap(pc,ovl); }
374: PetscOptionsName("-pc_asm_in_place","Perform matrix factorization inplace","PCASMSetUseInPlace",&flg);
375: if (flg) {PCASMSetUseInPlace(pc); }
376: PetscOptionsEList("-pc_asm_type","Type of restriction/extension","PCASMSetType",type,4,"restrict",buff,15,&flg);
377: if (flg) {
378: PCASMType atype;
379: PetscTruth isbasic,isrestrict,isinterpolate,isnone;
381: PetscStrcmp(buff,type[0],&isbasic);
382: PetscStrcmp(buff,type[1],&isrestrict);
383: PetscStrcmp(buff,type[2],&isinterpolate);
384: PetscStrcmp(buff,type[3],&isnone);
386: if (isbasic) atype = PC_ASM_BASIC;
387: else if (isrestrict) atype = PC_ASM_RESTRICT;
388: else if (isinterpolate) atype = PC_ASM_INTERPOLATE;
389: else if (isnone) atype = PC_ASM_NONE;
390: else SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Unknown type");
391: PCASMSetType(pc,atype);
392: }
393: PetscOptionsTail();
394: return(0);
395: }
397: /*------------------------------------------------------------------------------------*/
399: EXTERN_C_BEGIN
400: #undef __FUNCT__
402: int PCASMSetLocalSubdomains_ASM(PC pc,int n,IS *is)
403: {
404: PC_ASM *osm;
408: if (pc->setupcalled) {
409: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"PCASMSetLocalSubdomains() should be called before calling PCSetup().");
410: }
412: if (n < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Each process must have 0 or more blocks");
413: osm = (PC_ASM*)pc->data;
414: osm->n_local_true = n;
415: osm->is = is;
416: return(0);
417: }
418: EXTERN_C_END
420: EXTERN_C_BEGIN
421: #undef __FUNCT__
423: int PCASMSetTotalSubdomains_ASM(PC pc,int N,IS *is)
424: {
425: PC_ASM *osm;
426: int rank,size,ierr;
429: if (pc->setupcalled) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"PCASMSetTotalSubdomains() should be called before PCSetup().");
431: if (is) SETERRQ(PETSC_ERR_SUP,"Use PCASMSetLocalSubdomains to set specific index setsn
432: they cannot be set globally yet.");
434: osm = (PC_ASM*)pc->data;
435: /*
436: Split the subdomains equally amoung all processors
437: */
438: MPI_Comm_rank(pc->comm,&rank);
439: MPI_Comm_size(pc->comm,&size);
440: osm->n_local_true = N/size + ((N % size) > rank);
441: if (osm->n_local_true < 0) {
442: SETERRQ(PETSC_ERR_SUP,"Each process must have 0 or more blocks");
443: }
444: osm->is = 0;
445: return(0);
446: }
447: EXTERN_C_END
449: EXTERN_C_BEGIN
450: #undef __FUNCT__
452: int PCASMSetOverlap_ASM(PC pc,int ovl)
453: {
454: PC_ASM *osm;
457: if (ovl < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested");
459: osm = (PC_ASM*)pc->data;
460: osm->overlap = ovl;
461: return(0);
462: }
463: EXTERN_C_END
465: EXTERN_C_BEGIN
466: #undef __FUNCT__
468: int PCASMSetType_ASM(PC pc,PCASMType type)
469: {
470: PC_ASM *osm;
473: osm = (PC_ASM*)pc->data;
474: osm->type = type;
475: return(0);
476: }
477: EXTERN_C_END
479: EXTERN_C_BEGIN
480: #undef __FUNCT__
482: int PCASMGetSubSLES_ASM(PC pc,int *n_local,int *first_local,SLES **sles)
483: {
484: PC_ASM *jac = (PC_ASM*)pc->data;
485: int ierr;
488: if (jac->n_local_true < 0) {
489: SETERRQ(1,"Need to call PCSetUP() on PC (or SLESSetUp() on the outer SLES object) before calling here");
490: }
492: if (n_local) *n_local = jac->n_local_true;
493: if (first_local) {
494: MPI_Scan(&jac->n_local_true,first_local,1,MPI_INT,MPI_SUM,pc->comm);
495: *first_local -= jac->n_local_true;
496: }
497: *sles = jac->sles;
498: jac->same_local_solves = PETSC_FALSE; /* Assume that local solves are now different;
499: not necessarily true though! This flag is
500: used only for PCView_ASM() */
501: return(0);
502: }
503: EXTERN_C_END
505: EXTERN_C_BEGIN
506: #undef __FUNCT__
508: int PCASMSetUseInPlace_ASM(PC pc)
509: {
510: PC_ASM *dir;
513: dir = (PC_ASM*)pc->data;
514: dir->inplace = PETSC_TRUE;
515: return(0);
516: }
517: EXTERN_C_END
519: /*----------------------------------------------------------------------------*/
520: #undef __FUNCT__
522: /*@
523: PCASMSetUseInPlace - Tells the system to destroy the matrix after setup is done.
525: Collective on PC
527: Input Parameters:
528: . pc - the preconditioner context
530: Options Database Key:
531: . -pc_asm_in_place - Activates in-place factorization
533: Note:
534: PCASMSetUseInplace() can only be used with the KSP method KSPPREONLY, and
535: when the original matrix is not required during the Solve process.
536: This destroys the matrix, early thus, saving on memory usage.
538: Level: intermediate
540: .keywords: PC, set, factorization, direct, inplace, in-place, ASM
542: .seealso: PCILUSetUseInPlace(), PCLUSetUseInPlace ()
543: @*/
544: int PCASMSetUseInPlace(PC pc)
545: {
546: int ierr,(*f)(PC);
550: PetscObjectQueryFunction((PetscObject)pc,"PCASMSetUseInPlace_C",(void (**)(void))&f);
551: if (f) {
552: (*f)(pc);
553: }
554: return(0);
555: }
556: /*----------------------------------------------------------------------------*/
558: #undef __FUNCT__
560: /*@C
561: PCASMSetLocalSubdomains - Sets the local subdomains (for this processor
562: only) for the additive Schwarz preconditioner.
564: Collective on PC
566: Input Parameters:
567: + pc - the preconditioner context
568: . n - the number of subdomains for this processor (default value = 1)
569: - is - the index sets that define the subdomains for this processor
570: (or PETSC_NULL for PETSc to determine subdomains)
572: Notes:
573: The IS numbering is in the parallel, global numbering of the vector.
575: By default the ASM preconditioner uses 1 block per processor.
577: These index sets cannot be destroyed until after completion of the
578: linear solves for which the ASM preconditioner is being used.
580: Use PCASMSetTotalSubdomains() to set the subdomains for all processors.
582: Level: advanced
584: .keywords: PC, ASM, set, local, subdomains, additive Schwarz
586: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubSLES(),
587: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains()
588: @*/
589: int PCASMSetLocalSubdomains(PC pc,int n,IS *is)
590: {
591: int ierr,(*f)(PC,int,IS *);
595: PetscObjectQueryFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",(void (**)(void))&f);
596: if (f) {
597: (*f)(pc,n,is);
598: }
599: return(0);
600: }
602: #undef __FUNCT__
604: /*@C
605: PCASMSetTotalSubdomains - Sets the subdomains for all processor for the
606: additive Schwarz preconditioner. Either all or no processors in the
607: PC communicator must call this routine, with the same index sets.
609: Collective on PC
611: Input Parameters:
612: + pc - the preconditioner context
613: . n - the number of subdomains for all processors
614: - is - the index sets that define the subdomains for all processor
615: (or PETSC_NULL for PETSc to determine subdomains)
617: Options Database Key:
618: To set the total number of subdomain blocks rather than specify the
619: index sets, use the option
620: . -pc_asm_blocks <blks> - Sets total blocks
622: Notes:
623: Currently you cannot use this to set the actual subdomains with the argument is.
625: By default the ASM preconditioner uses 1 block per processor.
627: These index sets cannot be destroyed until after completion of the
628: linear solves for which the ASM preconditioner is being used.
630: Use PCASMSetLocalSubdomains() to set local subdomains.
632: Level: advanced
634: .keywords: PC, ASM, set, total, global, subdomains, additive Schwarz
636: .seealso: PCASMSetLocalSubdomains(), PCASMSetOverlap(), PCASMGetSubSLES(),
637: PCASMCreateSubdomains2D()
638: @*/
639: int PCASMSetTotalSubdomains(PC pc,int N,IS *is)
640: {
641: int ierr,(*f)(PC,int,IS *);
645: PetscObjectQueryFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",(void (**)(void))&f);
646: if (f) {
647: (*f)(pc,N,is);
648: }
649: return(0);
650: }
652: #undef __FUNCT__
654: /*@
655: PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
656: additive Schwarz preconditioner. Either all or no processors in the
657: PC communicator must call this routine.
659: Collective on PC
661: Input Parameters:
662: + pc - the preconditioner context
663: - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1)
665: Options Database Key:
666: . -pc_asm_overlap <ovl> - Sets overlap
668: Notes:
669: By default the ASM preconditioner uses 1 block per processor. To use
670: multiple blocks per perocessor, see PCASMSetTotalSubdomains() and
671: PCASMSetLocalSubdomains() (and the option -pc_asm_blocks <blks>).
673: The overlap defaults to 1, so if one desires that no additional
674: overlap be computed beyond what may have been set with a call to
675: PCASMSetTotalSubdomains() or PCASMSetLocalSubdomains(), then ovl
676: must be set to be 0. In particular, if one does not explicitly set
677: the subdomains an application code, then all overlap would be computed
678: internally by PETSc, and using an overlap of 0 would result in an ASM
679: variant that is equivalent to the block Jacobi preconditioner.
681: Note that one can define initial index sets with any overlap via
682: PCASMSetTotalSubdomains() or PCASMSetLocalSubdomains(); the routine
683: PCASMSetOverlap() merely allows PETSc to extend that overlap further
684: if desired.
686: Level: intermediate
688: .keywords: PC, ASM, set, overlap
690: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubSLES(),
691: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains()
692: @*/
693: int PCASMSetOverlap(PC pc,int ovl)
694: {
695: int ierr,(*f)(PC,int);
699: PetscObjectQueryFunction((PetscObject)pc,"PCASMSetOverlap_C",(void (**)(void))&f);
700: if (f) {
701: (*f)(pc,ovl);
702: }
703: return(0);
704: }
706: #undef __FUNCT__
708: /*@
709: PCASMSetType - Sets the type of restriction and interpolation used
710: for local problems in the additive Schwarz method.
712: Collective on PC
714: Input Parameters:
715: + pc - the preconditioner context
716: - type - variant of ASM, one of
717: .vb
718: PC_ASM_BASIC - full interpolation and restriction
719: PC_ASM_RESTRICT - full restriction, local processor interpolation
720: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
721: PC_ASM_NONE - local processor restriction and interpolation
722: .ve
724: Options Database Key:
725: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
727: Level: intermediate
729: .keywords: PC, ASM, set, type
731: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubSLES(),
732: PCASMCreateSubdomains2D()
733: @*/
734: int PCASMSetType(PC pc,PCASMType type)
735: {
736: int ierr,(*f)(PC,PCASMType);
740: PetscObjectQueryFunction((PetscObject)pc,"PCASMSetType_C",(void (**)(void))&f);
741: if (f) {
742: (*f)(pc,type);
743: }
744: return(0);
745: }
747: #undef __FUNCT__
749: /*@C
750: PCASMGetSubSLES - Gets the local SLES contexts for all blocks on
751: this processor.
752:
753: Collective on PC iff first_local is requested
755: Input Parameter:
756: . pc - the preconditioner context
758: Output Parameters:
759: + n_local - the number of blocks on this processor or PETSC_NULL
760: . first_local - the global number of the first block on this processor or PETSC_NULL,
761: all processors must request or all must pass PETSC_NULL
762: - sles - the array of SLES contexts
764: Note:
765: After PCASMGetSubSLES() the array of SLESes is not to be freed
767: Currently for some matrix implementations only 1 block per processor
768: is supported.
769:
770: You must call SLESSetUp() before calling PCASMGetSubSLES().
772: Level: advanced
774: .keywords: PC, ASM, additive Schwarz, get, sub, SLES, context
776: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap(),
777: PCASMCreateSubdomains2D(),
778: @*/
779: int PCASMGetSubSLES(PC pc,int *n_local,int *first_local,SLES **sles)
780: {
781: int ierr,(*f)(PC,int*,int*,SLES **);
785: PetscObjectQueryFunction((PetscObject)pc,"PCASMGetSubSLES_C",(void (**)(void))&f);
786: if (f) {
787: (*f)(pc,n_local,first_local,sles);
788: } else {
789: SETERRQ(1,"Cannot get subsles for this type of PC");
790: }
792: return(0);
793: }
795: /* -------------------------------------------------------------------------------------*/
796: /*MC
797: PCASM - Use the additive Schwarz method, each block is (approximately) solved with
798: its own SLES object.
800: Options Database Keys:
801: + -pc_asm_truelocal - Activates PCASMSetUseTrueLocal()
802: . -pc_asm_in_place - Activates in-place factorization
803: . -pc_asm_blocks <blks> - Sets total blocks
804: . -pc_asm_overlap <ovl> - Sets overlap
805: - -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
807: Notes: Each processor can have one or more blocks, but a block cannot be shared by more
808: than one processor. Defaults to one block per processor.
810: To set options on the solvers for each block append -sub_ to all the SLES, KSP, and PC
811: options database keys. For example, -sub_pc_type ilu -sub_pc_ilu_levels 1 -sub_ksp_type preonly
812:
813: To set the options on the solvers seperate for each block call PCASMGetSubSLES()
814: and set the options directly on the resulting SLES object (you can access its KSP and PC
815: with SLESGetKSP() and SLESGetPC())
817: Level: beginner
819: Concepts: additive Schwarz method
821: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
822: PCBJACOBI, PCASMSetUseTrueLocal(), PCASMGetSubSLES(), PCASMSetLocalSubdomains(),
823: PCASMSetTotalSubdomains(), PCSetModifySubmatrices(), PCASMSetOverlap(), PCASMSetType(),
824: PCASMSetUseInPlace()
825: M*/
827: EXTERN_C_BEGIN
828: #undef __FUNCT__
830: int PCCreate_ASM(PC pc)
831: {
832: int ierr;
833: PC_ASM *osm;
836: PetscNew(PC_ASM,&osm);
837: PetscLogObjectMemory(pc,sizeof(PC_ASM));
838: PetscMemzero(osm,sizeof(PC_ASM));
839: osm->n = PETSC_DECIDE;
840: osm->n_local = 0;
841: osm->n_local_true = PETSC_DECIDE;
842: osm->overlap = 1;
843: osm->is_flg = PETSC_FALSE;
844: osm->sles = 0;
845: osm->scat = 0;
846: osm->is = 0;
847: osm->mat = 0;
848: osm->pmat = 0;
849: osm->type = PC_ASM_RESTRICT;
850: osm->same_local_solves = PETSC_TRUE;
851: osm->inplace = PETSC_FALSE;
852: pc->data = (void*)osm;
854: pc->ops->apply = PCApply_ASM;
855: pc->ops->applytranspose = PCApplyTranspose_ASM;
856: pc->ops->setup = PCSetUp_ASM;
857: pc->ops->destroy = PCDestroy_ASM;
858: pc->ops->setfromoptions = PCSetFromOptions_ASM;
859: pc->ops->setuponblocks = PCSetUpOnBlocks_ASM;
860: pc->ops->view = PCView_ASM;
861: pc->ops->applyrichardson = 0;
863: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMSetLocalSubdomains_C","PCASMSetLocalSubdomains_ASM",
864: PCASMSetLocalSubdomains_ASM);
865: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMSetTotalSubdomains_C","PCASMSetTotalSubdomains_ASM",
866: PCASMSetTotalSubdomains_ASM);
867: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMSetOverlap_C","PCASMSetOverlap_ASM",
868: PCASMSetOverlap_ASM);
869: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMSetType_C","PCASMSetType_ASM",
870: PCASMSetType_ASM);
871: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMGetSubSLES_C","PCASMGetSubSLES_ASM",
872: PCASMGetSubSLES_ASM);
873: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMSetUseInPlace_C","PCASMSetUseInPlace_ASM",
874: PCASMSetUseInPlace_ASM);
875: return(0);
876: }
877: EXTERN_C_END
880: #undef __FUNCT__
882: /*@
883: PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
884: preconditioner for a two-dimensional problem on a regular grid.
886: Not Collective
888: Input Parameters:
889: + m, n - the number of mesh points in the x and y directions
890: . M, N - the number of subdomains in the x and y directions
891: . dof - degrees of freedom per node
892: - overlap - overlap in mesh lines
894: Output Parameters:
895: + Nsub - the number of subdomains created
896: - is - the array of index sets defining the subdomains
898: Note:
899: Presently PCAMSCreateSubdomains2d() is valid only for sequential
900: preconditioners. More general related routines are
901: PCASMSetTotalSubdomains() and PCASMSetLocalSubdomains().
903: Level: advanced
905: .keywords: PC, ASM, additive Schwarz, create, subdomains, 2D, regular grid
907: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubSLES(),
908: PCASMSetOverlap()
909: @*/
910: int PCASMCreateSubdomains2D(int m,int n,int M,int N,int dof,int overlap,int *Nsub,IS **is)
911: {
912: int i,j,height,width,ystart,xstart,yleft,yright,xleft,xright,loc_outter;
913: int nidx,*idx,loc,ii,jj,ierr,count;
916: if (dof != 1) SETERRQ(PETSC_ERR_SUP," ");
918: *Nsub = N*M;
919: PetscMalloc((*Nsub)*sizeof(IS **),is);
920: ystart = 0;
921: loc_outter = 0;
922: for (i=0; i<N; i++) {
923: height = n/N + ((n % N) > i); /* height of subdomain */
924: if (height < 2) SETERRQ(1,"Too many N subdomains for mesh dimension n");
925: yleft = ystart - overlap; if (yleft < 0) yleft = 0;
926: yright = ystart + height + overlap; if (yright > n) yright = n;
927: xstart = 0;
928: for (j=0; j<M; j++) {
929: width = m/M + ((m % M) > j); /* width of subdomain */
930: if (width < 2) SETERRQ(1,"Too many M subdomains for mesh dimension m");
931: xleft = xstart - overlap; if (xleft < 0) xleft = 0;
932: xright = xstart + width + overlap; if (xright > m) xright = m;
933: nidx = (xright - xleft)*(yright - yleft);
934: PetscMalloc(nidx*sizeof(int),&idx);
935: loc = 0;
936: for (ii=yleft; ii<yright; ii++) {
937: count = m*ii + xleft;
938: for (jj=xleft; jj<xright; jj++) {
939: idx[loc++] = count++;
940: }
941: }
942: ISCreateGeneral(PETSC_COMM_SELF,nidx,idx,(*is)+loc_outter++);
943: PetscFree(idx);
944: xstart += width;
945: }
946: ystart += height;
947: }
948: for (i=0; i<*Nsub; i++) { ISSort((*is)[i]); }
949: return(0);
950: }
952: #undef __FUNCT__
954: /*@C
955: PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
956: only) for the additive Schwarz preconditioner.
958: Collective on PC
960: Input Parameter:
961: . pc - the preconditioner context
963: Output Parameters:
964: + n - the number of subdomains for this processor (default value = 1)
965: - is - the index sets that define the subdomains for this processor
966:
968: Notes:
969: The IS numbering is in the parallel, global numbering of the vector.
971: Level: advanced
973: .keywords: PC, ASM, set, local, subdomains, additive Schwarz
975: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubSLES(),
976: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubmatrices()
977: @*/
978: int PCASMGetLocalSubdomains(PC pc,int *n,IS **is)
979: {
980: PC_ASM *osm;
984: if (!pc->setupcalled) {
985: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must call after SLESSetUP() or PCSetUp().");
986: }
988: osm = (PC_ASM*)pc->data;
989: if (n) *n = osm->n_local_true;
990: if (is) *is = osm->is;
991: return(0);
992: }
994: #undef __FUNCT__
996: /*@C
997: PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
998: only) for the additive Schwarz preconditioner.
1000: Collective on PC
1002: Input Parameter:
1003: . pc - the preconditioner context
1005: Output Parameters:
1006: + n - the number of matrices for this processor (default value = 1)
1007: - mat - the matrices
1008:
1010: Level: advanced
1012: .keywords: PC, ASM, set, local, subdomains, additive Schwarz, block Jacobi
1014: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubSLES(),
1015: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1016: @*/
1017: int PCASMGetLocalSubmatrices(PC pc,int *n,Mat **mat)
1018: {
1019: PC_ASM *osm;
1023: if (!pc->setupcalled) {
1024: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must call after SLESSetUP() or PCSetUp().");
1025: }
1027: osm = (PC_ASM*)pc->data;
1028: if (n) *n = osm->n_local_true;
1029: if (mat) *mat = osm->pmat;
1030: return(0);
1031: }