Actual source code: asm.c
1: #define PETSCKSP_DLL
3: /*
4: This file defines an additive Schwarz preconditioner for any Mat implementation.
6: Note that each processor may have any number of subdomains. But in order to
7: deal easily with the VecScatter(), we treat each processor as if it has the
8: same number of subdomains.
10: n - total number of true subdomains on all processors
11: n_local_true - actual number of subdomains on this processor
12: n_local = maximum over all processors of n_local_true
13: */
14: #include private/pcimpl.h
16: typedef struct {
17: PetscInt n,n_local,n_local_true;
18: PetscTruth is_flg; /* flg set to 1 if the IS created in pcsetup */
19: PetscInt overlap; /* overlap requested by user */
20: KSP *ksp; /* 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 type_set; /* if user set this value (so won't change it for symmetric problems) */
27: PetscTruth same_local_solves; /* flag indicating whether all local solvers are same */
28: PetscTruth inplace; /* indicates that the sub-matrices are deleted after
29: PCSetUpOnBlocks() is done. Similar to inplace
30: factorization in the case of LU and ILU */
31: } PC_ASM;
35: static PetscErrorCode PCView_ASM(PC pc,PetscViewer viewer)
36: {
37: PC_ASM *jac = (PC_ASM*)pc->data;
39: PetscMPIInt rank;
40: PetscInt i;
41: PetscTruth iascii,isstring;
42: PetscViewer sviewer;
46: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
47: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_STRING,&isstring);
48: if (iascii) {
49: if (jac->n > 0) {
50: PetscViewerASCIIPrintf(viewer," Additive Schwarz: total subdomain blocks = %D, amount of overlap = %D\n",jac->n,jac->overlap);
51: } else {
52: PetscViewerASCIIPrintf(viewer," Additive Schwarz: total subdomain blocks not yet set, amount of overlap = %D\n",jac->overlap);
53: }
54: PetscViewerASCIIPrintf(viewer," Additive Schwarz: restriction/interpolation type - %s\n",PCASMTypes[jac->type]);
55: MPI_Comm_rank(pc->comm,&rank);
56: if (jac->same_local_solves) {
57: PetscViewerASCIIPrintf(viewer," Local solve is same for all blocks, in the following KSP and PC objects:\n");
58: PetscViewerGetSingleton(viewer,&sviewer);
59: if (!rank && jac->ksp) {
60: PetscViewerASCIIPushTab(viewer);
61: KSPView(jac->ksp[0],sviewer);
62: PetscViewerASCIIPopTab(viewer);
63: }
64: PetscViewerRestoreSingleton(viewer,&sviewer);
65: } else {
66: PetscViewerASCIIPrintf(viewer," Local solve info for each block is in the following KSP and PC objects:\n");
67: PetscViewerASCIISynchronizedPrintf(viewer,"[%d] number of local blocks = %D\n",rank,jac->n_local_true);
68: PetscViewerASCIIPushTab(viewer);
69: for (i=0; i<jac->n_local; i++) {
70: PetscViewerGetSingleton(viewer,&sviewer);
71: if (i < jac->n_local_true) {
72: PetscViewerASCIISynchronizedPrintf(sviewer,"[%d] local block number %D\n",rank,i);
73: KSPView(jac->ksp[i],sviewer);
74: PetscViewerASCIISynchronizedPrintf(viewer,"- - - - - - - - - - - - - - - - - -\n");
75: }
76: PetscViewerRestoreSingleton(viewer,&sviewer);
77: }
78: PetscViewerASCIIPopTab(viewer);
79: PetscViewerFlush(viewer);
80: }
81: } else if (isstring) {
82: PetscViewerStringSPrintf(viewer," blks=%D, overlap=%D, type=%D",jac->n,jac->overlap,jac->type);
83: PetscViewerGetSingleton(viewer,&sviewer);
84: if (jac->ksp) {KSPView(jac->ksp[0],sviewer);}
85: PetscViewerGetSingleton(viewer,&sviewer);
86: } else {
87: SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCASM",((PetscObject)viewer)->type_name);
88: }
89: return(0);
90: }
94: static PetscErrorCode PCSetUp_ASM(PC pc)
95: {
96: PC_ASM *osm = (PC_ASM*)pc->data;
98: PetscInt i,m,n_local = osm->n_local,n_local_true = osm->n_local_true;
99: PetscInt start,start_val,end_val,sz,bs;
100: PetscMPIInt size;
101: MatReuse scall = MAT_REUSE_MATRIX;
102: IS isl;
103: KSP ksp;
104: PC subpc;
105: const char *prefix,*pprefix;
106: Vec vec;
109: MatGetVecs(pc->pmat,&vec,0);
110: if (!pc->setupcalled) {
111: if (osm->n == PETSC_DECIDE && osm->n_local_true == PETSC_DECIDE) {
112: /* no subdomains given, use one per processor */
113: osm->n_local_true = osm->n_local = 1;
114: MPI_Comm_size(pc->comm,&size);
115: osm->n = size;
116: } else if (osm->n == PETSC_DECIDE) { /* determine global number of subdomains */
117: PetscInt inwork[2],outwork[2];
118: inwork[0] = inwork[1] = osm->n_local_true;
119: MPI_Allreduce(inwork,outwork,1,MPIU_2INT,PetscMaxSum_Op,pc->comm);
120: osm->n_local = outwork[0];
121: osm->n = outwork[1];
122: }
123: n_local = osm->n_local;
124: n_local_true = osm->n_local_true;
125: if (!osm->is){ /* build the index sets */
126: PetscMalloc((n_local_true+1)*sizeof(IS **),&osm->is);
127: MatGetOwnershipRange(pc->pmat,&start_val,&end_val);
128: MatGetBlockSize(pc->pmat,&bs);
129: sz = end_val - start_val;
130: start = start_val;
131: if (end_val/bs*bs != end_val || start_val/bs*bs != start_val) {
132: SETERRQ(PETSC_ERR_ARG_WRONG,"Bad distribution for matrix block size");
133: }
134: for (i=0; i<n_local_true; i++){
135: size = ((sz/bs)/n_local_true + (((sz/bs) % n_local_true) > i))*bs;
136: ISCreateStride(PETSC_COMM_SELF,size,start,1,&isl);
137: start += size;
138: osm->is[i] = isl;
139: }
140: osm->is_flg = PETSC_TRUE;
141: }
143: PetscMalloc((n_local_true+1)*sizeof(KSP **),&osm->ksp);
144: PetscMalloc(n_local*sizeof(VecScatter **),&osm->scat);
145: PetscMalloc(2*n_local*sizeof(Vec **),&osm->x);
146: osm->y = osm->x + n_local;
148: /* Extend the "overlapping" regions by a number of steps */
149: MatIncreaseOverlap(pc->pmat,n_local_true,osm->is,osm->overlap);
150: for (i=0; i<n_local_true; i++) {
151: ISSort(osm->is[i]);
152: }
154: /* create the local work vectors and scatter contexts */
155: for (i=0; i<n_local_true; i++) {
156: ISGetLocalSize(osm->is[i],&m);
157: VecCreateSeq(PETSC_COMM_SELF,m,&osm->x[i]);
158: VecDuplicate(osm->x[i],&osm->y[i]);
159: ISCreateStride(PETSC_COMM_SELF,m,0,1,&isl);
160: VecScatterCreate(vec,osm->is[i],osm->x[i],isl,&osm->scat[i]);
161: ISDestroy(isl);
162: }
163: for (i=n_local_true; i<n_local; i++) {
164: VecCreateSeq(PETSC_COMM_SELF,0,&osm->x[i]);
165: VecDuplicate(osm->x[i],&osm->y[i]);
166: ISCreateStride(PETSC_COMM_SELF,0,0,1,&isl);
167: VecScatterCreate(vec,isl,osm->x[i],isl,&osm->scat[i]);
168: ISDestroy(isl);
169: }
171: /*
172: Create the local solvers.
173: */
174: for (i=0; i<n_local_true; i++) {
175: KSPCreate(PETSC_COMM_SELF,&ksp);
176: PetscLogObjectParent(pc,ksp);
177: KSPSetType(ksp,KSPPREONLY);
178: KSPGetPC(ksp,&subpc);
179: PCGetOptionsPrefix(pc,&prefix);
180: KSPSetOptionsPrefix(ksp,prefix);
181: KSPAppendOptionsPrefix(ksp,"sub_");
182: osm->ksp[i] = ksp;
183: }
184: scall = MAT_INITIAL_MATRIX;
185: } else {
186: /*
187: Destroy the blocks from the previous iteration
188: */
189: if (pc->flag == DIFFERENT_NONZERO_PATTERN) {
190: MatDestroyMatrices(osm->n_local_true,&osm->pmat);
191: scall = MAT_INITIAL_MATRIX;
192: }
193: }
195: /* extract out the submatrices */
196: MatGetSubMatrices(pc->pmat,osm->n_local_true,osm->is,osm->is,scall,&osm->pmat);
198: /* Return control to the user so that the submatrices can be modified (e.g., to apply
199: different boundary conditions for the submatrices than for the global problem) */
200: PCModifySubMatrices(pc,osm->n_local,osm->is,osm->is,osm->pmat,pc->modifysubmatricesP);
202: /* loop over subdomains putting them into local ksp */
203: PetscObjectGetOptionsPrefix((PetscObject)pc->pmat,&pprefix);
204: for (i=0; i<n_local_true; i++) {
205: PetscObjectSetOptionsPrefix((PetscObject)osm->pmat[i],pprefix);
206: PetscLogObjectParent(pc,osm->pmat[i]);
207: KSPSetOperators(osm->ksp[i],osm->pmat[i],osm->pmat[i],pc->flag);
208: KSPSetFromOptions(osm->ksp[i]);
209: }
210: VecDestroy(vec);
211: return(0);
212: }
216: static PetscErrorCode PCSetUpOnBlocks_ASM(PC pc)
217: {
218: PC_ASM *osm = (PC_ASM*)pc->data;
220: PetscInt i;
223: for (i=0; i<osm->n_local_true; i++) {
224: KSPSetUp(osm->ksp[i]);
225: }
226: /*
227: If inplace flag is set, then destroy the matrix after the setup
228: on blocks is done.
229: */
230: if (osm->inplace && osm->n_local_true > 0) {
231: MatDestroyMatrices(osm->n_local_true,&osm->pmat);
232: }
233: return(0);
234: }
238: static PetscErrorCode PCApply_ASM(PC pc,Vec x,Vec y)
239: {
240: PC_ASM *osm = (PC_ASM*)pc->data;
242: PetscInt i,n_local = osm->n_local,n_local_true = osm->n_local_true;
243: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
246: /*
247: Support for limiting the restriction or interpolation to only local
248: subdomain values (leaving the other values 0).
249: */
250: if (!(osm->type & PC_ASM_RESTRICT)) {
251: forward = SCATTER_FORWARD_LOCAL;
252: /* have to zero the work RHS since scatter may leave some slots empty */
253: for (i=0; i<n_local_true; i++) {
254: VecSet(osm->x[i],0.0);
255: }
256: }
257: if (!(osm->type & PC_ASM_INTERPOLATE)) {
258: reverse = SCATTER_REVERSE_LOCAL;
259: }
261: for (i=0; i<n_local; i++) {
262: VecScatterBegin(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
263: }
264: VecSet(y,0.0);
265: /* do the local solves */
266: for (i=0; i<n_local_true; i++) {
267: VecScatterEnd(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
268: KSPSolve(osm->ksp[i],osm->x[i],osm->y[i]);
269: VecScatterBegin(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
270: }
271: /* handle the rest of the scatters that do not have local solves */
272: for (i=n_local_true; i<n_local; i++) {
273: VecScatterEnd(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
274: VecScatterBegin(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
275: }
276: for (i=0; i<n_local; i++) {
277: VecScatterEnd(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
278: }
279: return(0);
280: }
284: static PetscErrorCode PCApplyTranspose_ASM(PC pc,Vec x,Vec y)
285: {
286: PC_ASM *osm = (PC_ASM*)pc->data;
288: PetscInt i,n_local = osm->n_local,n_local_true = osm->n_local_true;
289: ScatterMode forward = SCATTER_FORWARD,reverse = SCATTER_REVERSE;
292: /*
293: Support for limiting the restriction or interpolation to only local
294: subdomain values (leaving the other values 0).
296: Note: these are reversed from the PCApply_ASM() because we are applying the
297: transpose of the three terms
298: */
299: if (!(osm->type & PC_ASM_INTERPOLATE)) {
300: forward = SCATTER_FORWARD_LOCAL;
301: /* have to zero the work RHS since scatter may leave some slots empty */
302: for (i=0; i<n_local_true; i++) {
303: VecSet(osm->x[i],0.0);
304: }
305: }
306: if (!(osm->type & PC_ASM_RESTRICT)) {
307: reverse = SCATTER_REVERSE_LOCAL;
308: }
310: for (i=0; i<n_local; i++) {
311: VecScatterBegin(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
312: }
313: VecSet(y,0.0);
314: /* do the local solves */
315: for (i=0; i<n_local_true; i++) {
316: VecScatterEnd(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
317: KSPSolveTranspose(osm->ksp[i],osm->x[i],osm->y[i]);
318: VecScatterBegin(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
319: }
320: /* handle the rest of the scatters that do not have local solves */
321: for (i=n_local_true; i<n_local; i++) {
322: VecScatterEnd(x,osm->x[i],INSERT_VALUES,forward,osm->scat[i]);
323: VecScatterBegin(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
324: }
325: for (i=0; i<n_local; i++) {
326: VecScatterEnd(osm->y[i],y,ADD_VALUES,reverse,osm->scat[i]);
327: }
328: return(0);
329: }
333: static PetscErrorCode PCDestroy_ASM(PC pc)
334: {
335: PC_ASM *osm = (PC_ASM*)pc->data;
337: PetscInt i;
340: for (i=0; i<osm->n_local; i++) {
341: VecScatterDestroy(osm->scat[i]);
342: VecDestroy(osm->x[i]);
343: VecDestroy(osm->y[i]);
344: }
345: if (osm->n_local_true > 0 && !osm->inplace && osm->pmat) {
346: MatDestroyMatrices(osm->n_local_true,&osm->pmat);
347: }
348: if (osm->ksp) {
349: for (i=0; i<osm->n_local_true; i++) {
350: KSPDestroy(osm->ksp[i]);
351: }
352: }
353: if (osm->is_flg) {
354: for (i=0; i<osm->n_local_true; i++) {ISDestroy(osm->is[i]);}
355: PetscFree(osm->is);
356: }
357: PetscFree(osm->ksp);
358: PetscFree(osm->scat);
359: PetscFree(osm->x);
360: PetscFree(osm);
361: return(0);
362: }
366: static PetscErrorCode PCSetFromOptions_ASM(PC pc)
367: {
368: PC_ASM *osm = (PC_ASM*)pc->data;
370: PetscInt blocks,ovl;
371: PetscTruth flg,set,sym;
374: /* set the type to symmetric if matrix is symmetric */
375: if (pc->pmat && !osm->type_set) {
376: MatIsSymmetricKnown(pc->pmat,&set,&sym);
377: if (set && sym) {
378: osm->type = PC_ASM_BASIC;
379: }
380: }
381: PetscOptionsHead("Additive Schwarz options");
382: PetscOptionsInt("-pc_asm_blocks","Number of subdomains","PCASMSetTotalSubdomains",osm->n,&blocks,&flg);
383: if (flg) {PCASMSetTotalSubdomains(pc,blocks,PETSC_NULL); }
384: PetscOptionsInt("-pc_asm_overlap","Number of grid points overlap","PCASMSetOverlap",osm->overlap,&ovl,&flg);
385: if (flg) {PCASMSetOverlap(pc,ovl); }
386: PetscOptionsName("-pc_asm_in_place","Perform matrix factorization inplace","PCASMSetUseInPlace",&flg);
387: if (flg) {PCASMSetUseInPlace(pc); }
388: PetscOptionsEnum("-pc_asm_type","Type of restriction/extension","PCASMSetType",PCASMTypes,(PetscEnum)osm->type,(PetscEnum*)&osm->type,&flg);
389: PetscOptionsTail();
390: return(0);
391: }
393: /*------------------------------------------------------------------------------------*/
398: PetscErrorCode PCASMSetLocalSubdomains_ASM(PC pc,PetscInt n,IS is[])
399: {
400: PC_ASM *osm = (PC_ASM*)pc->data;
403: if (n < 1) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Each process must have 0 or more blocks");
405: if (pc->setupcalled && n != osm->n_local_true) {
406: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"PCASMSetLocalSubdomains() should be called before calling PCSetup().");
407: }
408: if (!pc->setupcalled){
409: osm->n_local_true = n;
410: osm->is = is;
411: }
412: return(0);
413: }
419: PetscErrorCode PCASMSetTotalSubdomains_ASM(PC pc,PetscInt N,IS *is)
420: {
421: PC_ASM *osm = (PC_ASM*)pc->data;
423: PetscMPIInt rank,size;
424: PetscInt n;
428: if (is) SETERRQ(PETSC_ERR_SUP,"Use PCASMSetLocalSubdomains() to set specific index sets\n\they cannot be set globally yet.");
430: /*
431: Split the subdomains equally amoung all processors
432: */
433: MPI_Comm_rank(pc->comm,&rank);
434: MPI_Comm_size(pc->comm,&size);
435: n = N/size + ((N % size) > rank);
436: if (pc->setupcalled && n != osm->n_local_true) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"PCASMSetTotalSubdomains() should be called before PCSetup().");
437: if (!pc->setupcalled){
438: osm->n_local_true = n;
439: if (!n) SETERRQ2(PETSC_ERR_ARG_OUTOFRANGE,"Each process must have at least one block: total processors %d total blocks %d",size,(int)n);
440: osm->is = 0;
441: }
442: return(0);
443: }
449: PetscErrorCode PCASMSetOverlap_ASM(PC pc,PetscInt ovl)
450: {
451: PC_ASM *osm;
454: if (ovl < 0) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Negative overlap value requested");
456: osm = (PC_ASM*)pc->data;
457: osm->overlap = ovl;
458: return(0);
459: }
465: PetscErrorCode PCASMSetType_ASM(PC pc,PCASMType type)
466: {
467: PC_ASM *osm;
470: osm = (PC_ASM*)pc->data;
471: osm->type = type;
472: osm->type_set = PETSC_TRUE;
473: return(0);
474: }
480: PetscErrorCode PCASMGetSubKSP_ASM(PC pc,PetscInt *n_local,PetscInt *first_local,KSP **ksp)
481: {
482: PC_ASM *jac = (PC_ASM*)pc->data;
486: if (jac->n_local_true < 0) {
487: SETERRQ(PETSC_ERR_ORDER,"Need to call PCSetUP() on PC (or KSPSetUp() on the outer KSP object) before calling here");
488: }
490: if (n_local) *n_local = jac->n_local_true;
491: if (first_local) {
492: MPI_Scan(&jac->n_local_true,first_local,1,MPIU_INT,MPI_SUM,pc->comm);
493: *first_local -= jac->n_local_true;
494: }
495: *ksp = jac->ksp;
496: jac->same_local_solves = PETSC_FALSE; /* Assume that local solves are now different;
497: not necessarily true though! This flag is
498: used only for PCView_ASM() */
499: return(0);
500: }
506: PetscErrorCode PCASMSetUseInPlace_ASM(PC pc)
507: {
508: PC_ASM *dir;
511: dir = (PC_ASM*)pc->data;
512: dir->inplace = PETSC_TRUE;
513: return(0);
514: }
517: /*----------------------------------------------------------------------------*/
520: /*@
521: PCASMSetUseInPlace - Tells the system to destroy the matrix after setup is done.
523: Collective on PC
525: Input Parameters:
526: . pc - the preconditioner context
528: Options Database Key:
529: . -pc_asm_in_place - Activates in-place factorization
531: Note:
532: PCASMSetUseInplace() can only be used with the KSP method KSPPREONLY, and
533: when the original matrix is not required during the Solve process.
534: This destroys the matrix, early thus, saving on memory usage.
536: Level: intermediate
538: .keywords: PC, set, factorization, direct, inplace, in-place, ASM
540: .seealso: PCFactorSetUseInPlace()
541: @*/
542: PetscErrorCode PCASMSetUseInPlace(PC pc)
543: {
544: PetscErrorCode ierr,(*f)(PC);
548: PetscObjectQueryFunction((PetscObject)pc,"PCASMSetUseInPlace_C",(void (**)(void))&f);
549: if (f) {
550: (*f)(pc);
551: }
552: return(0);
553: }
554: /*----------------------------------------------------------------------------*/
558: /*@C
559: PCASMSetLocalSubdomains - Sets the local subdomains (for this processor
560: only) for the additive Schwarz preconditioner.
562: Collective on PC
564: Input Parameters:
565: + pc - the preconditioner context
566: . n - the number of subdomains for this processor (default value = 1)
567: - is - the index sets that define the subdomains for this processor
568: (or PETSC_NULL for PETSc to determine subdomains)
570: Notes:
571: The IS numbering is in the parallel, global numbering of the vector.
573: By default the ASM preconditioner uses 1 block per processor.
575: These index sets cannot be destroyed until after completion of the
576: linear solves for which the ASM preconditioner is being used.
578: Use PCASMSetTotalSubdomains() to set the subdomains for all processors.
580: Level: advanced
582: .keywords: PC, ASM, set, local, subdomains, additive Schwarz
584: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
585: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains()
586: @*/
587: PetscErrorCode PCASMSetLocalSubdomains(PC pc,PetscInt n,IS is[])
588: {
589: PetscErrorCode ierr,(*f)(PC,PetscInt,IS[]);
593: PetscObjectQueryFunction((PetscObject)pc,"PCASMSetLocalSubdomains_C",(void (**)(void))&f);
594: if (f) {
595: (*f)(pc,n,is);
596: }
597: return(0);
598: }
602: /*@C
603: PCASMSetTotalSubdomains - Sets the subdomains for all processor for the
604: additive Schwarz preconditioner. Either all or no processors in the
605: PC communicator must call this routine, with the same index sets.
607: Collective on PC
609: Input Parameters:
610: + pc - the preconditioner context
611: . n - the number of subdomains for all processors
612: - is - the index sets that define the subdomains for all processor
613: (or PETSC_NULL for PETSc to determine subdomains)
615: Options Database Key:
616: To set the total number of subdomain blocks rather than specify the
617: index sets, use the option
618: . -pc_asm_blocks <blks> - Sets total blocks
620: Notes:
621: Currently you cannot use this to set the actual subdomains with the argument is.
623: By default the ASM preconditioner uses 1 block per processor.
625: These index sets cannot be destroyed until after completion of the
626: linear solves for which the ASM preconditioner is being used.
628: Use PCASMSetLocalSubdomains() to set local subdomains.
630: Level: advanced
632: .keywords: PC, ASM, set, total, global, subdomains, additive Schwarz
634: .seealso: PCASMSetLocalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
635: PCASMCreateSubdomains2D()
636: @*/
637: PetscErrorCode PCASMSetTotalSubdomains(PC pc,PetscInt N,IS *is)
638: {
639: PetscErrorCode ierr,(*f)(PC,PetscInt,IS *);
643: PetscObjectQueryFunction((PetscObject)pc,"PCASMSetTotalSubdomains_C",(void (**)(void))&f);
644: if (f) {
645: (*f)(pc,N,is);
646: }
647: return(0);
648: }
652: /*@
653: PCASMSetOverlap - Sets the overlap between a pair of subdomains for the
654: additive Schwarz preconditioner. Either all or no processors in the
655: PC communicator must call this routine.
657: Collective on PC
659: Input Parameters:
660: + pc - the preconditioner context
661: - ovl - the amount of overlap between subdomains (ovl >= 0, default value = 1)
663: Options Database Key:
664: . -pc_asm_overlap <ovl> - Sets overlap
666: Notes:
667: By default the ASM preconditioner uses 1 block per processor. To use
668: multiple blocks per perocessor, see PCASMSetTotalSubdomains() and
669: PCASMSetLocalSubdomains() (and the option -pc_asm_blocks <blks>).
671: The overlap defaults to 1, so if one desires that no additional
672: overlap be computed beyond what may have been set with a call to
673: PCASMSetTotalSubdomains() or PCASMSetLocalSubdomains(), then ovl
674: must be set to be 0. In particular, if one does not explicitly set
675: the subdomains an application code, then all overlap would be computed
676: internally by PETSc, and using an overlap of 0 would result in an ASM
677: variant that is equivalent to the block Jacobi preconditioner.
679: Note that one can define initial index sets with any overlap via
680: PCASMSetTotalSubdomains() or PCASMSetLocalSubdomains(); the routine
681: PCASMSetOverlap() merely allows PETSc to extend that overlap further
682: if desired.
684: Level: intermediate
686: .keywords: PC, ASM, set, overlap
688: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
689: PCASMCreateSubdomains2D(), PCASMGetLocalSubdomains()
690: @*/
691: PetscErrorCode PCASMSetOverlap(PC pc,PetscInt ovl)
692: {
693: PetscErrorCode ierr,(*f)(PC,PetscInt);
697: PetscObjectQueryFunction((PetscObject)pc,"PCASMSetOverlap_C",(void (**)(void))&f);
698: if (f) {
699: (*f)(pc,ovl);
700: }
701: return(0);
702: }
706: /*@
707: PCASMSetType - Sets the type of restriction and interpolation used
708: for local problems in the additive Schwarz method.
710: Collective on PC
712: Input Parameters:
713: + pc - the preconditioner context
714: - type - variant of ASM, one of
715: .vb
716: PC_ASM_BASIC - full interpolation and restriction
717: PC_ASM_RESTRICT - full restriction, local processor interpolation
718: PC_ASM_INTERPOLATE - full interpolation, local processor restriction
719: PC_ASM_NONE - local processor restriction and interpolation
720: .ve
722: Options Database Key:
723: . -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
725: Level: intermediate
727: .keywords: PC, ASM, set, type
729: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMGetSubKSP(),
730: PCASMCreateSubdomains2D()
731: @*/
732: PetscErrorCode PCASMSetType(PC pc,PCASMType type)
733: {
734: PetscErrorCode ierr,(*f)(PC,PCASMType);
738: PetscObjectQueryFunction((PetscObject)pc,"PCASMSetType_C",(void (**)(void))&f);
739: if (f) {
740: (*f)(pc,type);
741: }
742: return(0);
743: }
747: /*@C
748: PCASMGetSubKSP - Gets the local KSP contexts for all blocks on
749: this processor.
750:
751: Collective on PC iff first_local is requested
753: Input Parameter:
754: . pc - the preconditioner context
756: Output Parameters:
757: + n_local - the number of blocks on this processor or PETSC_NULL
758: . first_local - the global number of the first block on this processor or PETSC_NULL,
759: all processors must request or all must pass PETSC_NULL
760: - ksp - the array of KSP contexts
762: Note:
763: After PCASMGetSubKSP() the array of KSPes is not to be freed
765: Currently for some matrix implementations only 1 block per processor
766: is supported.
767:
768: You must call KSPSetUp() before calling PCASMGetSubKSP().
770: Level: advanced
772: .keywords: PC, ASM, additive Schwarz, get, sub, KSP, context
774: .seealso: PCASMSetTotalSubdomains(), PCASMSetTotalSubdomains(), PCASMSetOverlap(),
775: PCASMCreateSubdomains2D(),
776: @*/
777: PetscErrorCode PCASMGetSubKSP(PC pc,PetscInt *n_local,PetscInt *first_local,KSP *ksp[])
778: {
779: PetscErrorCode ierr,(*f)(PC,PetscInt*,PetscInt*,KSP **);
783: PetscObjectQueryFunction((PetscObject)pc,"PCASMGetSubKSP_C",(void (**)(void))&f);
784: if (f) {
785: (*f)(pc,n_local,first_local,ksp);
786: } else {
787: SETERRQ(PETSC_ERR_ARG_WRONG,"Cannot get subksp for this type of PC");
788: }
790: return(0);
791: }
793: /* -------------------------------------------------------------------------------------*/
794: /*MC
795: PCASM - Use the (restricted) additive Schwarz method, each block is (approximately) solved with
796: its own KSP object.
798: Options Database Keys:
799: + -pc_asm_truelocal - Activates PCASMSetUseTrueLocal()
800: . -pc_asm_in_place - Activates in-place factorization
801: . -pc_asm_blocks <blks> - Sets total blocks
802: . -pc_asm_overlap <ovl> - Sets overlap
803: - -pc_asm_type [basic,restrict,interpolate,none] - Sets ASM type
805: IMPORTANT: If you run with, for example, 3 blocks on 1 processor or 3 blocks on 3 processors you
806: will get a different convergence rate due to the default option of -pc_asm_type restrict. Use
807: -pc_asm_type basic to use the standard ASM.
809: Notes: Each processor can have one or more blocks, but a block cannot be shared by more
810: than one processor. Defaults to one block per processor.
812: To set options on the solvers for each block append -sub_ to all the KSP, and PC
813: options database keys. For example, -sub_pc_type ilu -sub_pc_factor_levels 1 -sub_ksp_type preonly
814:
815: To set the options on the solvers separate for each block call PCASMGetSubKSP()
816: and set the options directly on the resulting KSP object (you can access its PC
817: with KSPGetPC())
820: Level: beginner
822: Concepts: additive Schwarz method
824: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
825: PCBJACOBI, PCASMSetUseTrueLocal(), PCASMGetSubKSP(), PCASMSetLocalSubdomains(),
826: PCASMSetTotalSubdomains(), PCSetModifySubmatrices(), PCASMSetOverlap(), PCASMSetType(),
827: PCASMSetUseInPlace()
828: M*/
833: PetscErrorCode PCCreate_ASM(PC pc)
834: {
836: PC_ASM *osm;
839: PetscNew(PC_ASM,&osm);
840: PetscLogObjectMemory(pc,sizeof(PC_ASM));
841: osm->n = PETSC_DECIDE;
842: osm->n_local = 0;
843: osm->n_local_true = PETSC_DECIDE;
844: osm->overlap = 1;
845: osm->is_flg = PETSC_FALSE;
846: osm->ksp = 0;
847: osm->scat = 0;
848: osm->is = 0;
849: osm->mat = 0;
850: osm->pmat = 0;
851: osm->type = PC_ASM_RESTRICT;
852: osm->same_local_solves = PETSC_TRUE;
853: osm->inplace = PETSC_FALSE;
854: pc->data = (void*)osm;
856: pc->ops->apply = PCApply_ASM;
857: pc->ops->applytranspose = PCApplyTranspose_ASM;
858: pc->ops->setup = PCSetUp_ASM;
859: pc->ops->destroy = PCDestroy_ASM;
860: pc->ops->setfromoptions = PCSetFromOptions_ASM;
861: pc->ops->setuponblocks = PCSetUpOnBlocks_ASM;
862: pc->ops->view = PCView_ASM;
863: pc->ops->applyrichardson = 0;
865: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMSetLocalSubdomains_C","PCASMSetLocalSubdomains_ASM",
866: PCASMSetLocalSubdomains_ASM);
867: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMSetTotalSubdomains_C","PCASMSetTotalSubdomains_ASM",
868: PCASMSetTotalSubdomains_ASM);
869: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMSetOverlap_C","PCASMSetOverlap_ASM",
870: PCASMSetOverlap_ASM);
871: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMSetType_C","PCASMSetType_ASM",
872: PCASMSetType_ASM);
873: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMGetSubKSP_C","PCASMGetSubKSP_ASM",
874: PCASMGetSubKSP_ASM);
875: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCASMSetUseInPlace_C","PCASMSetUseInPlace_ASM",
876: PCASMSetUseInPlace_ASM);
877: return(0);
878: }
884: /*@
885: PCASMCreateSubdomains2D - Creates the index sets for the overlapping Schwarz
886: preconditioner for a two-dimensional problem on a regular grid.
888: Not Collective
890: Input Parameters:
891: + m, n - the number of mesh points in the x and y directions
892: . M, N - the number of subdomains in the x and y directions
893: . dof - degrees of freedom per node
894: - overlap - overlap in mesh lines
896: Output Parameters:
897: + Nsub - the number of subdomains created
898: - is - the array of index sets defining the subdomains
900: Note:
901: Presently PCAMSCreateSubdomains2d() is valid only for sequential
902: preconditioners. More general related routines are
903: PCASMSetTotalSubdomains() and PCASMSetLocalSubdomains().
905: Level: advanced
907: .keywords: PC, ASM, additive Schwarz, create, subdomains, 2D, regular grid
909: .seealso: PCASMSetTotalSubdomains(), PCASMSetLocalSubdomains(), PCASMGetSubKSP(),
910: PCASMSetOverlap()
911: @*/
912: PetscErrorCode PCASMCreateSubdomains2D(PetscInt m,PetscInt n,PetscInt M,PetscInt N,PetscInt dof,PetscInt overlap,PetscInt *Nsub,IS **is)
913: {
914: PetscInt i,j,height,width,ystart,xstart,yleft,yright,xleft,xright,loc_outter;
916: PetscInt nidx,*idx,loc,ii,jj,count;
919: if (dof != 1) SETERRQ(PETSC_ERR_SUP," ");
921: *Nsub = N*M;
922: PetscMalloc((*Nsub)*sizeof(IS **),is);
923: ystart = 0;
924: loc_outter = 0;
925: for (i=0; i<N; i++) {
926: height = n/N + ((n % N) > i); /* height of subdomain */
927: if (height < 2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Too many N subdomains for mesh dimension n");
928: yleft = ystart - overlap; if (yleft < 0) yleft = 0;
929: yright = ystart + height + overlap; if (yright > n) yright = n;
930: xstart = 0;
931: for (j=0; j<M; j++) {
932: width = m/M + ((m % M) > j); /* width of subdomain */
933: if (width < 2) SETERRQ(PETSC_ERR_ARG_OUTOFRANGE,"Too many M subdomains for mesh dimension m");
934: xleft = xstart - overlap; if (xleft < 0) xleft = 0;
935: xright = xstart + width + overlap; if (xright > m) xright = m;
936: nidx = (xright - xleft)*(yright - yleft);
937: PetscMalloc(nidx*sizeof(PetscInt),&idx);
938: loc = 0;
939: for (ii=yleft; ii<yright; ii++) {
940: count = m*ii + xleft;
941: for (jj=xleft; jj<xright; jj++) {
942: idx[loc++] = count++;
943: }
944: }
945: ISCreateGeneral(PETSC_COMM_SELF,nidx,idx,(*is)+loc_outter++);
946: PetscFree(idx);
947: xstart += width;
948: }
949: ystart += height;
950: }
951: for (i=0; i<*Nsub; i++) { ISSort((*is)[i]); }
952: return(0);
953: }
957: /*@C
958: PCASMGetLocalSubdomains - Gets the local subdomains (for this processor
959: only) for the additive Schwarz preconditioner.
961: Collective on PC
963: Input Parameter:
964: . pc - the preconditioner context
966: Output Parameters:
967: + n - the number of subdomains for this processor (default value = 1)
968: - is - the index sets that define the subdomains for this processor
969:
971: Notes:
972: The IS numbering is in the parallel, global numbering of the vector.
974: Level: advanced
976: .keywords: PC, ASM, set, local, subdomains, additive Schwarz
978: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
979: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubmatrices()
980: @*/
981: PetscErrorCode PCASMGetLocalSubdomains(PC pc,PetscInt *n,IS *is[])
982: {
983: PC_ASM *osm;
988: if (!pc->setupcalled) {
989: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUP() or PCSetUp().");
990: }
992: osm = (PC_ASM*)pc->data;
993: if (n) *n = osm->n_local_true;
994: if (is) *is = osm->is;
995: return(0);
996: }
1000: /*@C
1001: PCASMGetLocalSubmatrices - Gets the local submatrices (for this processor
1002: only) for the additive Schwarz preconditioner.
1004: Collective on PC
1006: Input Parameter:
1007: . pc - the preconditioner context
1009: Output Parameters:
1010: + n - the number of matrices for this processor (default value = 1)
1011: - mat - the matrices
1012:
1014: Level: advanced
1016: .keywords: PC, ASM, set, local, subdomains, additive Schwarz, block Jacobi
1018: .seealso: PCASMSetTotalSubdomains(), PCASMSetOverlap(), PCASMGetSubKSP(),
1019: PCASMCreateSubdomains2D(), PCASMSetLocalSubdomains(), PCASMGetLocalSubdomains()
1020: @*/
1021: PetscErrorCode PCASMGetLocalSubmatrices(PC pc,PetscInt *n,Mat *mat[])
1022: {
1023: PC_ASM *osm;
1028: if (!pc->setupcalled) {
1029: SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must call after KSPSetUP() or PCSetUp().");
1030: }
1032: osm = (PC_ASM*)pc->data;
1033: if (n) *n = osm->n_local_true;
1034: if (mat) *mat = osm->pmat;
1035: return(0);
1036: }