Actual source code: mg.c

  1: #define PETSCKSP_DLL

  3: /*
  4:     Defines the multigrid preconditioner interface.
  5: */
 6:  #include src/ksp/pc/impls/mg/mgimpl.h


 11: PetscErrorCode PCMGMCycle_Private(PC_MG **mglevels,PetscTruth *converged)
 12: {
 13:   PC_MG          *mg = *mglevels,*mgc;
 15:   PetscInt       cycles = (PetscInt) mg->cycles;

 18:   if (converged) *converged = PETSC_FALSE;

 21:   KSPSolve(mg->smoothd,mg->b,mg->x);  /* pre-smooth */
 23:   if (mg->level) {  /* not the coarsest grid */
 24:     (*mg->residual)(mg->A,mg->b,mg->x,mg->r);

 26:     /* if on finest level and have convergence criteria set */
 27:     if (mg->level == mg->levels-1 && mg->ttol) {
 28:       PetscReal rnorm;
 29:       VecNorm(mg->r,NORM_2,&rnorm);
 30:       if (rnorm <= mg->ttol) {
 31:         *converged = PETSC_TRUE;
 32:         if (rnorm < mg->abstol) {
 33:           PetscInfo2(0,"Linear solver has converged. Residual norm %G is less than absolute tolerance %G\n",rnorm,mg->abstol);
 34:         } else {
 35:           PetscInfo2(0,"Linear solver has converged. Residual norm %G is less than relative tolerance times initial residual norm %G\n",rnorm,mg->ttol);
 36:         }
 37:         return(0);
 38:       }
 39:     }

 41:     mgc = *(mglevels - 1);
 42:     MatRestrict(mg->restrct,mg->r,mgc->b);
 43:     VecSet(mgc->x,0.0);
 44:     while (cycles--) {
 45:       PCMGMCycle_Private(mglevels-1,converged);
 46:     }
 47:     MatInterpolateAdd(mg->interpolate,mgc->x,mg->x,mg->x);
 49:     KSPSolve(mg->smoothu,mg->b,mg->x);    /* post smooth */
 51:   }
 52:   return(0);
 53: }

 55: /*
 56:        PCMGCreate_Private - Creates a PC_MG structure for use with the
 57:                multigrid code. Level 0 is the coarsest. (But the 
 58:                finest level is stored first in the array).

 60: */
 63: static PetscErrorCode PCMGCreate_Private(MPI_Comm comm,PetscInt levels,PC pc,MPI_Comm *comms,PC_MG ***result)
 64: {
 65:   PC_MG          **mg;
 67:   PetscInt       i;
 68:   PetscMPIInt    size;
 69:   const char     *prefix;
 70:   PC             ipc;

 73:   PetscMalloc(levels*sizeof(PC_MG*),&mg);
 74:   PetscLogObjectMemory(pc,levels*(sizeof(PC_MG*)+sizeof(PC_MG)));

 76:   PCGetOptionsPrefix(pc,&prefix);

 78:   for (i=0; i<levels; i++) {
 79:     PetscNew(PC_MG,&mg[i]);
 80:     mg[i]->level           = i;
 81:     mg[i]->levels          = levels;
 82:     mg[i]->cycles          = PC_MG_CYCLE_V;
 83:     mg[i]->galerkin        = PETSC_FALSE;
 84:     mg[i]->galerkinused    = PETSC_FALSE;
 85:     mg[i]->default_smoothu = 1;
 86:     mg[i]->default_smoothd = 1;

 88:     if (comms) comm = comms[i];
 89:     KSPCreate(comm,&mg[i]->smoothd);
 90:     KSPSetTolerances(mg[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, mg[i]->default_smoothd);
 91:     KSPSetOptionsPrefix(mg[i]->smoothd,prefix);

 93:     /* do special stuff for coarse grid */
 94:     if (!i && levels > 1) {
 95:       KSPAppendOptionsPrefix(mg[0]->smoothd,"mg_coarse_");

 97:       /* coarse solve is (redundant) LU by default */
 98:       KSPSetType(mg[0]->smoothd,KSPPREONLY);
 99:       KSPGetPC(mg[0]->smoothd,&ipc);
100:       MPI_Comm_size(comm,&size);
101:       if (size > 1) {
102:         PCSetType(ipc,PCREDUNDANT);
103:         PCRedundantGetPC(ipc,&ipc);
104:       }
105:       PCSetType(ipc,PCLU);

107:     } else {
108:       char tprefix[128];
109:       sprintf(tprefix,"mg_levels_%d_",(int)i);
110:       KSPAppendOptionsPrefix(mg[i]->smoothd,tprefix);
111:     }
112:     PetscLogObjectParent(pc,mg[i]->smoothd);
113:     mg[i]->smoothu           = mg[i]->smoothd;
114:     mg[i]->rtol              = 0.0;
115:     mg[i]->abstol            = 0.0;
116:     mg[i]->dtol              = 0.0;
117:     mg[i]->ttol              = 0.0;
118:     mg[i]->eventsetup        = 0;
119:     mg[i]->eventsolve        = 0;
120:     mg[i]->cyclesperpcapply  = 1;
121:   }
122:   *result = mg;
123:   return(0);
124: }

128: static PetscErrorCode PCDestroy_MG(PC pc)
129: {
130:   PC_MG          **mg = (PC_MG**)pc->data;
132:   PetscInt       i,n = mg[0]->levels;

135:   for (i=0; i<n-1; i++) {
136:     if (mg[i+1]->r) {VecDestroy(mg[i+1]->r);}
137:     if (mg[i]->b) {VecDestroy(mg[i]->b);}
138:     if (mg[i]->x) {VecDestroy(mg[i]->x);}
139:     if (mg[i+1]->restrct) {MatDestroy(mg[i+1]->restrct);}
140:     if (mg[i+1]->interpolate) {MatDestroy(mg[i+1]->interpolate);}
141:   }

143:   for (i=0; i<n; i++) {
144:     if (mg[i]->smoothd != mg[i]->smoothu) {
145:       KSPDestroy(mg[i]->smoothd);
146:     }
147:     KSPDestroy(mg[i]->smoothu);
148:     PetscFree(mg[i]);
149:   }
150:   PetscFree(mg);
151:   return(0);
152: }



156: EXTERN PetscErrorCode PCMGACycle_Private(PC_MG**);
157: EXTERN PetscErrorCode PCMGFCycle_Private(PC_MG**);
158: EXTERN PetscErrorCode PCMGKCycle_Private(PC_MG**);

160: /*
161:    PCApply_MG - Runs either an additive, multiplicative, Kaskadic
162:              or full cycle of multigrid. 

164:   Note: 
165:   A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle(). 
166: */
169: static PetscErrorCode PCApply_MG(PC pc,Vec b,Vec x)
170: {
171:   PC_MG          **mg = (PC_MG**)pc->data;
173:   PetscInt       levels = mg[0]->levels,i;

176:   mg[levels-1]->b = b;
177:   mg[levels-1]->x = x;
178:   if (!mg[levels-1]->r && mg[0]->am != PC_MG_ADDITIVE && levels > 1) {
179:     Vec tvec;
180:     VecDuplicate(mg[levels-1]->b,&tvec);
181:     PCMGSetR(pc,levels-1,tvec);
182:     VecDestroy(tvec);
183:   }
184:   if (mg[0]->am == PC_MG_MULTIPLICATIVE) {
185:     VecSet(x,0.0);
186:     for (i=0; i<mg[0]->cyclesperpcapply; i++) {
187:       PCMGMCycle_Private(mg+levels-1,PETSC_NULL);
188:     }
189:   }
190:   else if (mg[0]->am == PC_MG_ADDITIVE) {
191:     PCMGACycle_Private(mg);
192:   }
193:   else if (mg[0]->am == PC_MG_KASKADE) {
194:     PCMGKCycle_Private(mg);
195:   }
196:   else {
197:     PCMGFCycle_Private(mg);
198:   }
199:   return(0);
200: }

204: static PetscErrorCode PCApplyRichardson_MG(PC pc,Vec b,Vec x,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its)
205: {
206:   PC_MG          **mg = (PC_MG**)pc->data;
208:   PetscInt       levels = mg[0]->levels;
209:   PetscTruth     converged = PETSC_FALSE;

212:   mg[levels-1]->b    = b;
213:   mg[levels-1]->x    = x;

215:   mg[levels-1]->rtol = rtol;
216:   mg[levels-1]->abstol = abstol;
217:   mg[levels-1]->dtol = dtol;
218:   if (rtol) {
219:     /* compute initial residual norm for relative convergence test */
220:     PetscReal rnorm;
221:     (*mg[levels-1]->residual)(mg[levels-1]->A,b,x,w);
222:     VecNorm(w,NORM_2,&rnorm);
223:     mg[levels-1]->ttol = PetscMax(rtol*rnorm,abstol);
224:   } else if (abstol) {
225:     mg[levels-1]->ttol = abstol;
226:   } else {
227:     mg[levels-1]->ttol = 0.0;
228:   }

230:   while (its-- && !converged) {
231:     PCMGMCycle_Private(mg+levels-1,&converged);
232:   }
233:   return(0);
234: }

238: PetscErrorCode PCSetFromOptions_MG(PC pc)
239: {
241:   PetscInt       m,levels = 1,cycles;
242:   PetscTruth     flg;
243:   PC_MG          **mg = (PC_MG**)pc->data;
244:   PCMGType       mgtype = PC_MG_ADDITIVE;
245:   PCMGCycleType  mgctype;

248:   PetscOptionsHead("Multigrid options");
249:     if (!pc->data) {
250:       PetscOptionsInt("-pc_mg_levels","Number of Levels","PCMGSetLevels",levels,&levels,&flg);
251:       PCMGSetLevels(pc,levels,PETSC_NULL);
252:       mg = (PC_MG**)pc->data;
253:     }
254:     mgctype = (PCMGCycleType) mg[0]->cycles;
255:     PetscOptionsEnum("-pc_mg_cycle_type","V cycle or for W-cycle","PCMGSetCycleType",PCMGCycleTypes,(PetscEnum)mgctype,(PetscEnum*)&mgctype,&flg);
256:     if (flg) {
257:       PCMGSetCycleType(pc,mgctype);
258:     };
259:     PetscOptionsName("-pc_mg_galerkin","Use Galerkin process to compute coarser operators","PCMGSetGalerkin",&flg);
260:     if (flg) {
261:       PCMGSetGalerkin(pc);
262:     }
263:     PetscOptionsInt("-pc_mg_smoothup","Number of post-smoothing steps","PCMGSetNumberSmoothUp",1,&m,&flg);
264:     if (flg) {
265:       PCMGSetNumberSmoothUp(pc,m);
266:     }
267:     PetscOptionsInt("-pc_mg_smoothdown","Number of pre-smoothing steps","PCMGSetNumberSmoothDown",1,&m,&flg);
268:     if (flg) {
269:       PCMGSetNumberSmoothDown(pc,m);
270:     }
271:     PetscOptionsEnum("-pc_mg_type","Multigrid type","PCMGSetType",PCMGTypes,(PetscEnum)mgtype,(PetscEnum*)&mgtype,&flg);
272:     if (flg) {
273:       PCMGSetType(pc,mgtype);
274:     }
275:     if (mg[0]->am == PC_MG_MULTIPLICATIVE) {
276:       PetscOptionsInt("-pc_mg_multiplicative_cycles","Number of cycles for each preconditioner step","PCMGSetLevels",mg[0]->cyclesperpcapply,&cycles,&flg);
277:       if (flg) {
278:         PCMGMultiplicativeSetCycles(pc,cycles);
279:       }
280:     }
281:     PetscOptionsName("-pc_mg_log","Log times for each multigrid level","None",&flg);
282:     if (flg) {
283:       PetscInt i;
284:       char     eventname[128];
285:       if (!mg) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
286:       levels = mg[0]->levels;
287:       for (i=0; i<levels; i++) {
288:         sprintf(eventname,"MSetup Level %d",(int)i);
290:         sprintf(eventname,"MGSolve Level %d to 0",(int)i);
292:       }
293:     }
294:   PetscOptionsTail();
295:   return(0);
296: }

298: const char *PCMGTypes[] = {"MULTIPLICATIVE","ADDITIVE","FULL","KASKADE","PCMGType","PC_MG",0};
299: const char *PCMGCycleTypes[] = {"invalid","v","w","PCMGCycleType","PC_MG_CYCLE",0};

303: static PetscErrorCode PCView_MG(PC pc,PetscViewer viewer)
304: {
305:   PC_MG          **mg = (PC_MG**)pc->data;
307:   PetscInt       levels = mg[0]->levels,i;
308:   PetscTruth     iascii;

311:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
312:   if (iascii) {
313:     PetscViewerASCIIPrintf(viewer,"  MG: type is %s, levels=%D cycles=%s, pre-smooths=%D, post-smooths=%D\n",
314:                                   PCMGTypes[mg[0]->am],levels,(mg[0]->cycles == PC_MG_CYCLE_V) ? "v" : "w",
315:                                   mg[0]->default_smoothd,mg[0]->default_smoothu);
316:     if (mg[0]->galerkin) {
317:       PetscViewerASCIIPrintf(viewer,"    Using Galerkin computed coarse grid matrices\n");
318:     }
319:     for (i=0; i<levels; i++) {
320:       if (!i) {
321:         PetscViewerASCIIPrintf(viewer,"Coarse gride solver -- level %D -------------------------------\n",i);
322:       } else {
323:         PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);
324:       }
325:       PetscViewerASCIIPushTab(viewer);
326:       KSPView(mg[i]->smoothd,viewer);
327:       PetscViewerASCIIPopTab(viewer);
328:       if (i && mg[i]->smoothd == mg[i]->smoothu) {
329:         PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");
330:       } else if (i){
331:         PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);
332:         PetscViewerASCIIPushTab(viewer);
333:         KSPView(mg[i]->smoothu,viewer);
334:         PetscViewerASCIIPopTab(viewer);
335:       }
336:     }
337:   } else {
338:     SETERRQ1(PETSC_ERR_SUP,"Viewer type %s not supported for PCMG",((PetscObject)viewer)->type_name);
339:   }
340:   return(0);
341: }

343: /*
344:     Calls setup for the KSP on each level
345: */
348: static PetscErrorCode PCSetUp_MG(PC pc)
349: {
350:   PC_MG                   **mg = (PC_MG**)pc->data;
351:   PetscErrorCode          ierr;
352:   PetscInt                i,n = mg[0]->levels;
353:   PC                      cpc;
354:   PetscTruth              preonly,lu,redundant,cholesky,monitor = PETSC_FALSE,dump,opsset;
355:   PetscViewerASCIIMonitor ascii;
356:   PetscViewer             viewer = PETSC_NULL;
357:   MPI_Comm                comm;
358:   Mat                     dA,dB;
359:   MatStructure            uflag;
360:   Vec                     tvec;


364:   /* If user did not provide fine grid operators OR operator was not updated since last global KSPSetOperators() */
365:   /* so use those from global PC */
366:   /* Is this what we always want? What if user wants to keep old one? */
367:   KSPGetOperatorsSet(mg[n-1]->smoothd,PETSC_NULL,&opsset);
368:   KSPGetPC(mg[0]->smoothd,&cpc);
369:   if (!opsset || cpc->setupcalled == 2) {
370:     PetscInfo(pc,"Using outer operators to define finest grid operator \n  because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
371:     KSPSetOperators(mg[n-1]->smoothd,pc->mat,pc->pmat,pc->flag);
372:   }

374:   if (mg[0]->galerkin) {
375:     Mat B;
376:     mg[0]->galerkinused = PETSC_TRUE;
377:     /* currently only handle case where mat and pmat are the same on coarser levels */
378:     KSPGetOperators(mg[n-1]->smoothd,&dA,&dB,&uflag);
379:     if (!pc->setupcalled) {
380:       for (i=n-2; i>-1; i--) {
381:         MatPtAP(dB,mg[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);
382:         KSPSetOperators(mg[i]->smoothd,B,B,uflag);
383:         if (i != n-2) {PetscObjectDereference((PetscObject)dB);}
384:         dB   = B;
385:       }
386:       PetscObjectDereference((PetscObject)dB);
387:     } else {
388:       for (i=n-2; i>-1; i--) {
389:         KSPGetOperators(mg[i]->smoothd,PETSC_NULL,&B,PETSC_NULL);
390:         MatPtAP(dB,mg[i+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);
391:         KSPSetOperators(mg[i]->smoothd,B,B,uflag);
392:         dB   = B;
393:       }
394:     }
395:   }

397:   if (!pc->setupcalled) {
398:     PetscOptionsHasName(0,"-pc_mg_monitor",&monitor);
399: 
400:     for (i=0; i<n; i++) {
401:       if (monitor) {
402:         PetscObjectGetComm((PetscObject)mg[i]->smoothd,&comm);
403:         PetscViewerASCIIMonitorCreate(comm,"stdout",n-i,&ascii);
404:         KSPMonitorSet(mg[i]->smoothd,KSPMonitorDefault,ascii,(PetscErrorCode(*)(void*))PetscViewerASCIIMonitorDestroy);
405:       }
406:       KSPSetFromOptions(mg[i]->smoothd);
407:     }
408:     for (i=1; i<n; i++) {
409:       if (mg[i]->smoothu && (mg[i]->smoothu != mg[i]->smoothd)) {
410:         if (monitor) {
411:           PetscObjectGetComm((PetscObject)mg[i]->smoothu,&comm);
412:           PetscViewerASCIIMonitorCreate(comm,"stdout",n-i,&ascii);
413:           KSPMonitorSet(mg[i]->smoothu,KSPMonitorDefault,ascii,(PetscErrorCode(*)(void*))PetscViewerASCIIMonitorDestroy);
414:         }
415:         KSPSetFromOptions(mg[i]->smoothu);
416:       }
417:     }
418:     for (i=1; i<n; i++) {
419:       if (!mg[i]->residual) {
420:         Mat mat;
421:         KSPGetOperators(mg[i]->smoothd,PETSC_NULL,&mat,PETSC_NULL);
422:         PCMGSetResidual(pc,i,PCMGDefaultResidual,mat);
423:       }
424:       if (mg[i]->restrct && !mg[i]->interpolate) {
425:         PCMGSetInterpolate(pc,i,mg[i]->restrct);
426:       }
427:       if (!mg[i]->restrct && mg[i]->interpolate) {
428:         PCMGSetRestriction(pc,i,mg[i]->interpolate);
429:       }
430: #if defined(PETSC_USE_DEBUG)
431:       if (!mg[i]->restrct || !mg[i]->interpolate) {
432:         SETERRQ1(PETSC_ERR_ARG_WRONGSTATE,"Need to set restriction or interpolation on level %d",(int)i);
433:       }
434: #endif
435:     }
436:     for (i=0; i<n-1; i++) {
437:       if (!mg[i]->b) {
438:         Vec *vec;
439:         KSPGetVecs(mg[i]->smoothd,1,&vec,0,PETSC_NULL);
440:         PCMGSetRhs(pc,i,*vec);
441:         PetscFree(vec);
442:       }
443:       if (!mg[i]->r && i) {
444:         VecDuplicate(mg[i]->b,&tvec);
445:         PCMGSetR(pc,i,tvec);
446:         VecDestroy(tvec);
447:       }
448:       if (!mg[i]->x) {
449:         VecDuplicate(mg[i]->b,&tvec);
450:         PCMGSetX(pc,i,tvec);
451:         VecDestroy(tvec);
452:       }
453:     }
454:   }


457:   for (i=1; i<n; i++) {
458:     if (mg[i]->smoothu == mg[i]->smoothd) {
459:       /* if doing only down then initial guess is zero */
460:       KSPSetInitialGuessNonzero(mg[i]->smoothd,PETSC_TRUE);
461:     }
463:     KSPSetUp(mg[i]->smoothd);
465:   }
466:   for (i=1; i<n; i++) {
467:     if (mg[i]->smoothu && mg[i]->smoothu != mg[i]->smoothd) {
468:       Mat          downmat,downpmat;
469:       MatStructure matflag;
470:       PetscTruth   opsset;

472:       /* check if operators have been set for up, if not use down operators to set them */
473:       KSPGetOperatorsSet(mg[i]->smoothu,&opsset,PETSC_NULL);
474:       if (!opsset) {
475:         KSPGetOperators(mg[i]->smoothd,&downmat,&downpmat,&matflag);
476:         KSPSetOperators(mg[i]->smoothu,downmat,downpmat,matflag);
477:       }

479:       KSPSetInitialGuessNonzero(mg[i]->smoothu,PETSC_TRUE);
481:       KSPSetUp(mg[i]->smoothu);
483:     }
484:   }

486:   /*
487:       If coarse solver is not direct method then DO NOT USE preonly 
488:   */
489:   PetscTypeCompare((PetscObject)mg[0]->smoothd,KSPPREONLY,&preonly);
490:   if (preonly) {
491:     PetscTypeCompare((PetscObject)cpc,PCLU,&lu);
492:     PetscTypeCompare((PetscObject)cpc,PCREDUNDANT,&redundant);
493:     PetscTypeCompare((PetscObject)cpc,PCCHOLESKY,&cholesky);
494:     if (!lu && !redundant && !cholesky) {
495:       KSPSetType(mg[0]->smoothd,KSPGMRES);
496:     }
497:   }

499:   if (!pc->setupcalled) {
500:     if (monitor) {
501:       PetscObjectGetComm((PetscObject)mg[0]->smoothd,&comm);
502:       PetscViewerASCIIMonitorCreate(comm,"stdout",n,&ascii);
503:       KSPMonitorSet(mg[0]->smoothd,KSPMonitorDefault,ascii,(PetscErrorCode(*)(void*))PetscViewerASCIIMonitorDestroy);
504:     }
505:     KSPSetFromOptions(mg[0]->smoothd);
506:   }

509:   KSPSetUp(mg[0]->smoothd);

512:   /*
513:      Dump the interpolation/restriction matrices plus the 
514:    Jacobian/stiffness on each level. This allows Matlab users to 
515:    easily check if the Galerkin condition A_c = R A_f R^T is satisfied.

517:    Only support one or the other at the same time.
518:   */
519: #if defined(PETSC_USE_SOCKET_VIEWER)
520:   PetscOptionsHasName(pc->prefix,"-pc_mg_dump_matlab",&dump);
521:   if (dump) {
522:     viewer = PETSC_VIEWER_SOCKET_(pc->comm);
523:   }
524: #endif
525:   PetscOptionsHasName(pc->prefix,"-pc_mg_dump_binary",&dump);
526:   if (dump) {
527:     viewer = PETSC_VIEWER_BINARY_(pc->comm);
528:   }

530:   if (viewer) {
531:     for (i=1; i<n; i++) {
532:       MatView(mg[i]->restrct,viewer);
533:     }
534:     for (i=0; i<n; i++) {
535:       KSPGetPC(mg[i]->smoothd,&pc);
536:       MatView(pc->mat,viewer);
537:     }
538:   }
539:   return(0);
540: }

542: /* -------------------------------------------------------------------------------------*/

546: /*@C
547:    PCMGSetLevels - Sets the number of levels to use with MG.
548:    Must be called before any other MG routine.

550:    Collective on PC

552:    Input Parameters:
553: +  pc - the preconditioner context
554: .  levels - the number of levels
555: -  comms - optional communicators for each level; this is to allow solving the coarser problems
556:            on smaller sets of processors. Use PETSC_NULL_OBJECT for default in Fortran

558:    Level: intermediate

560:    Notes:
561:      If the number of levels is one then the multigrid uses the -mg_levels prefix
562:   for setting the level options rather than the -mg_coarse prefix.

564: .keywords: MG, set, levels, multigrid

566: .seealso: PCMGSetType(), PCMGGetLevels()
567: @*/
568: PetscErrorCode  PCMGSetLevels(PC pc,PetscInt levels,MPI_Comm *comms)
569: {
571:   PC_MG          **mg=0;


576:   if (pc->data) {
577:     SETERRQ(PETSC_ERR_ORDER,"Number levels already set for MG\n\
578:     make sure that you call PCMGSetLevels() before KSPSetFromOptions()");
579:   }
580:   PCMGCreate_Private(pc->comm,levels,pc,comms,&mg);
581:   mg[0]->am                = PC_MG_MULTIPLICATIVE;
582:   pc->data                 = (void*)mg;
583:   pc->ops->applyrichardson = PCApplyRichardson_MG;
584:   return(0);
585: }

589: /*@
590:    PCMGGetLevels - Gets the number of levels to use with MG.

592:    Not Collective

594:    Input Parameter:
595: .  pc - the preconditioner context

597:    Output parameter:
598: .  levels - the number of levels

600:    Level: advanced

602: .keywords: MG, get, levels, multigrid

604: .seealso: PCMGSetLevels()
605: @*/
606: PetscErrorCode  PCMGGetLevels(PC pc,PetscInt *levels)
607: {
608:   PC_MG  **mg;


614:   mg      = (PC_MG**)pc->data;
615:   *levels = mg[0]->levels;
616:   return(0);
617: }

621: /*@
622:    PCMGSetType - Determines the form of multigrid to use:
623:    multiplicative, additive, full, or the Kaskade algorithm.

625:    Collective on PC

627:    Input Parameters:
628: +  pc - the preconditioner context
629: -  form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,
630:    PC_MG_FULL, PC_MG_KASKADE

632:    Options Database Key:
633: .  -pc_mg_type <form> - Sets <form>, one of multiplicative,
634:    additive, full, kaskade   

636:    Level: advanced

638: .keywords: MG, set, method, multiplicative, additive, full, Kaskade, multigrid

640: .seealso: PCMGSetLevels()
641: @*/
642: PetscErrorCode  PCMGSetType(PC pc,PCMGType form)
643: {
644:   PC_MG **mg;

648:   mg = (PC_MG**)pc->data;

650:   if (!mg) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
651:   mg[0]->am = form;
652:   if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG;
653:   else pc->ops->applyrichardson = 0;
654:   return(0);
655: }

659: /*@
660:    PCMGSetCycleType - Sets the type cycles to use.  Use PCMGSetCycleTypeOnLevel() for more 
661:    complicated cycling.

663:    Collective on PC

665:    Input Parameters:
666: +  pc - the multigrid context 
667: -  PC_MG_CYCLE_V or PC_MG_CYCLE_W

669:    Options Database Key:
670: $  -pc_mg_cycle_type v or w

672:    Level: advanced

674: .keywords: MG, set, cycles, V-cycle, W-cycle, multigrid

676: .seealso: PCMGSetCycleTypeOnLevel()
677: @*/
678: PetscErrorCode  PCMGSetCycleType(PC pc,PCMGCycleType n)
679: {
680:   PC_MG    **mg;
681:   PetscInt i,levels;

685:   mg     = (PC_MG**)pc->data;
686:   if (!mg) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
687:   levels = mg[0]->levels;

689:   for (i=0; i<levels; i++) {
690:     mg[i]->cycles  = n;
691:   }
692:   return(0);
693: }

697: /*@
698:    PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step 
699:          of multigrid when PCMGType of PC_MG_MULTIPLICATIVE is used

701:    Collective on PC

703:    Input Parameters:
704: +  pc - the multigrid context 
705: -  n - number of cycles (default is 1)

707:    Options Database Key:
708: $  -pc_mg_multiplicative_cycles n

710:    Level: advanced

712:    Notes: This is not associated with setting a v or w cycle, that is set with PCMGSetCycleType()

714: .keywords: MG, set, cycles, V-cycle, W-cycle, multigrid

716: .seealso: PCMGSetCycleTypeOnLevel(), PCMGSetCycleType()
717: @*/
718: PetscErrorCode  PCMGMultiplicativeSetCycles(PC pc,PetscInt n)
719: {
720:   PC_MG    **mg;
721:   PetscInt i,levels;

725:   mg     = (PC_MG**)pc->data;
726:   if (!mg) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
727:   levels = mg[0]->levels;

729:   for (i=0; i<levels; i++) {
730:     mg[i]->cyclesperpcapply  = n;
731:   }
732:   return(0);
733: }

737: /*@
738:    PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the
739:       finest grid via the Galerkin process: A_i-1 = r_i * A_i * r_i^t

741:    Collective on PC

743:    Input Parameters:
744: .  pc - the multigrid context 

746:    Options Database Key:
747: $  -pc_mg_galerkin

749:    Level: intermediate

751: .keywords: MG, set, Galerkin

753: .seealso: PCMGGetGalerkin()

755: @*/
756: PetscErrorCode  PCMGSetGalerkin(PC pc)
757: {
758:   PC_MG    **mg;
759:   PetscInt i,levels;

763:   mg     = (PC_MG**)pc->data;
764:   if (!mg) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
765:   levels = mg[0]->levels;

767:   for (i=0; i<levels; i++) {
768:     mg[i]->galerkin = PETSC_TRUE;
769:   }
770:   return(0);
771: }

775: /*@
776:    PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e.
777:       A_i-1 = r_i * A_i * r_i^t

779:    Not Collective

781:    Input Parameter:
782: .  pc - the multigrid context 

784:    Output Parameter:
785: .  gelerkin - PETSC_TRUE or PETSC_FALSE

787:    Options Database Key:
788: $  -pc_mg_galerkin

790:    Level: intermediate

792: .keywords: MG, set, Galerkin

794: .seealso: PCMGSetGalerkin()

796: @*/
797: PetscErrorCode  PCMGGetGalerkin(PC pc,PetscTruth *galerkin)
798: {
799:   PC_MG    **mg;

803:   mg     = (PC_MG**)pc->data;
804:   if (!mg) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
805:   *galerkin = mg[0]->galerkin;
806:   return(0);
807: }

811: /*@
812:    PCMGSetNumberSmoothDown - Sets the number of pre-smoothing steps to
813:    use on all levels. Use PCMGGetSmootherDown() to set different 
814:    pre-smoothing steps on different levels.

816:    Collective on PC

818:    Input Parameters:
819: +  mg - the multigrid context 
820: -  n - the number of smoothing steps

822:    Options Database Key:
823: .  -pc_mg_smoothdown <n> - Sets number of pre-smoothing steps

825:    Level: advanced

827: .keywords: MG, smooth, down, pre-smoothing, steps, multigrid

829: .seealso: PCMGSetNumberSmoothUp()
830: @*/
831: PetscErrorCode  PCMGSetNumberSmoothDown(PC pc,PetscInt n)
832: {
833:   PC_MG          **mg;
835:   PetscInt       i,levels;

839:   mg     = (PC_MG**)pc->data;
840:   if (!mg) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
841:   levels = mg[0]->levels;

843:   for (i=1; i<levels; i++) {
844:     /* make sure smoother up and down are different */
845:     PCMGGetSmootherUp(pc,i,PETSC_NULL);
846:     KSPSetTolerances(mg[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
847:     mg[i]->default_smoothd = n;
848:   }
849:   return(0);
850: }

854: /*@
855:    PCMGSetNumberSmoothUp - Sets the number of post-smoothing steps to use 
856:    on all levels. Use PCMGGetSmootherUp() to set different numbers of 
857:    post-smoothing steps on different levels.

859:    Collective on PC

861:    Input Parameters:
862: +  mg - the multigrid context 
863: -  n - the number of smoothing steps

865:    Options Database Key:
866: .  -pc_mg_smoothup <n> - Sets number of post-smoothing steps

868:    Level: advanced

870:    Note: this does not set a value on the coarsest grid, since we assume that
871:     there is no separate smooth up on the coarsest grid.

873: .keywords: MG, smooth, up, post-smoothing, steps, multigrid

875: .seealso: PCMGSetNumberSmoothDown()
876: @*/
877: PetscErrorCode  PCMGSetNumberSmoothUp(PC pc,PetscInt n)
878: {
879:   PC_MG          **mg;
881:   PetscInt       i,levels;

885:   mg     = (PC_MG**)pc->data;
886:   if (!mg) SETERRQ(PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
887:   levels = mg[0]->levels;

889:   for (i=1; i<levels; i++) {
890:     /* make sure smoother up and down are different */
891:     PCMGGetSmootherUp(pc,i,PETSC_NULL);
892:     KSPSetTolerances(mg[i]->smoothu,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
893:     mg[i]->default_smoothu = n;
894:   }
895:   return(0);
896: }

898: /* ----------------------------------------------------------------------------------------*/

900: /*MC
901:    PCMG - Use multigrid preconditioning. This preconditioner requires you provide additional
902:     information about the coarser grid matrices and restriction/interpolation operators.

904:    Options Database Keys:
905: +  -pc_mg_levels <nlevels> - number of levels including finest
906: .  -pc_mg_cycles v or w
907: .  -pc_mg_smoothup <n> - number of smoothing steps after interpolation
908: .  -pc_mg_smoothdown <n> - number of smoothing steps before applying restriction operator
909: .  -pc_mg_type <additive,multiplicative,full,cascade> - multiplicative is the default
910: .  -pc_mg_log - log information about time spent on each level of the solver
911: .  -pc_mg_monitor - print information on the multigrid convergence
912: .  -pc_mg_galerkin - use Galerkin process to compute coarser operators
913: -  -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices
914:                         to the Socket viewer for reading from Matlab.

916:    Notes:

918:    Level: intermediate

920:    Concepts: multigrid

922: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, 
923:            PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), PCMGSetNumberSmoothDown(),
924:            PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(),
925:            PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(),
926:            PCMGSetCycleTypeOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR()           
927: M*/

932: PetscErrorCode  PCCreate_MG(PC pc)
933: {
935:   pc->ops->apply          = PCApply_MG;
936:   pc->ops->setup          = PCSetUp_MG;
937:   pc->ops->destroy        = PCDestroy_MG;
938:   pc->ops->setfromoptions = PCSetFromOptions_MG;
939:   pc->ops->view           = PCView_MG;

941:   pc->data                = (void*)0;
942:   return(0);
943: }