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: }