Actual source code: mg.c

petsc-3.5.2 2014-09-08
Report Typos and Errors
  2: /*
  3:     Defines the multigrid preconditioner interface.
  4: */
  5: #include <../src/ksp/pc/impls/mg/mgimpl.h>                    /*I "petscksp.h" I*/
  6: #include <petscdm.h>

 10: PetscErrorCode PCMGMCycle_Private(PC pc,PC_MG_Levels **mglevelsin,PCRichardsonConvergedReason *reason)
 11: {
 12:   PC_MG          *mg = (PC_MG*)pc->data;
 13:   PC_MG_Levels   *mgc,*mglevels = *mglevelsin;
 15:   PetscInt       cycles = (mglevels->level == 1) ? 1 : (PetscInt) mglevels->cycles;

 18:   if (mglevels->eventsmoothsolve) {PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);}
 19:   KSPSolve(mglevels->smoothd,mglevels->b,mglevels->x);  /* pre-smooth */
 20:   if (mglevels->eventsmoothsolve) {PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);}
 21:   if (mglevels->level) {  /* not the coarsest grid */
 22:     if (mglevels->eventresidual) {PetscLogEventBegin(mglevels->eventresidual,0,0,0,0);}
 23:     (*mglevels->residual)(mglevels->A,mglevels->b,mglevels->x,mglevels->r);
 24:     if (mglevels->eventresidual) {PetscLogEventEnd(mglevels->eventresidual,0,0,0,0);}

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

 42:     mgc = *(mglevelsin - 1);
 43:     if (mglevels->eventinterprestrict) {PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);}
 44:     MatRestrict(mglevels->restrct,mglevels->r,mgc->b);
 45:     if (mglevels->eventinterprestrict) {PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);}
 46:     VecSet(mgc->x,0.0);
 47:     while (cycles--) {
 48:       PCMGMCycle_Private(pc,mglevelsin-1,reason);
 49:     }
 50:     if (mglevels->eventinterprestrict) {PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);}
 51:     MatInterpolateAdd(mglevels->interpolate,mgc->x,mglevels->x,mglevels->x);
 52:     if (mglevels->eventinterprestrict) {PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);}
 53:     if (mglevels->eventsmoothsolve) {PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);}
 54:     KSPSolve(mglevels->smoothu,mglevels->b,mglevels->x);    /* post smooth */
 55:     if (mglevels->eventsmoothsolve) {PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);}
 56:   }
 57:   return(0);
 58: }

 62: static PetscErrorCode PCApplyRichardson_MG(PC pc,Vec b,Vec x,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool zeroguess,PetscInt *outits,PCRichardsonConvergedReason *reason)
 63: {
 64:   PC_MG          *mg        = (PC_MG*)pc->data;
 65:   PC_MG_Levels   **mglevels = mg->levels;
 67:   PetscInt       levels = mglevels[0]->levels,i;

 70:   /* When the DM is supplying the matrix then it will not exist until here */
 71:   for (i=0; i<levels; i++) {
 72:     if (!mglevels[i]->A) {
 73:       KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);
 74:       PetscObjectReference((PetscObject)mglevels[i]->A);
 75:     }
 76:   }
 77:   mglevels[levels-1]->b = b;
 78:   mglevels[levels-1]->x = x;

 80:   mg->rtol   = rtol;
 81:   mg->abstol = abstol;
 82:   mg->dtol   = dtol;
 83:   if (rtol) {
 84:     /* compute initial residual norm for relative convergence test */
 85:     PetscReal rnorm;
 86:     if (zeroguess) {
 87:       VecNorm(b,NORM_2,&rnorm);
 88:     } else {
 89:       (*mglevels[levels-1]->residual)(mglevels[levels-1]->A,b,x,w);
 90:       VecNorm(w,NORM_2,&rnorm);
 91:     }
 92:     mg->ttol = PetscMax(rtol*rnorm,abstol);
 93:   } else if (abstol) mg->ttol = abstol;
 94:   else mg->ttol = 0.0;

 96:   /* since smoother is applied to full system, not just residual we need to make sure that smoothers don't
 97:      stop prematurely due to small residual */
 98:   for (i=1; i<levels; i++) {
 99:     KSPSetTolerances(mglevels[i]->smoothu,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
100:     if (mglevels[i]->smoothu != mglevels[i]->smoothd) {
101:       KSPSetTolerances(mglevels[i]->smoothd,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
102:     }
103:   }

105:   *reason = (PCRichardsonConvergedReason)0;
106:   for (i=0; i<its; i++) {
107:     PCMGMCycle_Private(pc,mglevels+levels-1,reason);
108:     if (*reason) break;
109:   }
110:   if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
111:   *outits = i;
112:   return(0);
113: }

117: PetscErrorCode PCReset_MG(PC pc)
118: {
119:   PC_MG          *mg        = (PC_MG*)pc->data;
120:   PC_MG_Levels   **mglevels = mg->levels;
122:   PetscInt       i,n;

125:   if (mglevels) {
126:     n = mglevels[0]->levels;
127:     for (i=0; i<n-1; i++) {
128:       VecDestroy(&mglevels[i+1]->r);
129:       VecDestroy(&mglevels[i]->b);
130:       VecDestroy(&mglevels[i]->x);
131:       MatDestroy(&mglevels[i+1]->restrct);
132:       MatDestroy(&mglevels[i+1]->interpolate);
133:       VecDestroy(&mglevels[i+1]->rscale);
134:     }

136:     for (i=0; i<n; i++) {
137:       MatDestroy(&mglevels[i]->A);
138:       if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
139:         KSPReset(mglevels[i]->smoothd);
140:       }
141:       KSPReset(mglevels[i]->smoothu);
142:     }
143:   }
144:   return(0);
145: }

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

153:    Logically Collective on PC

155:    Input Parameters:
156: +  pc - the preconditioner context
157: .  levels - the number of levels
158: -  comms - optional communicators for each level; this is to allow solving the coarser problems
159:            on smaller sets of processors. Use NULL_OBJECT for default in Fortran

161:    Level: intermediate

163:    Notes:
164:      If the number of levels is one then the multigrid uses the -mg_levels prefix
165:   for setting the level options rather than the -mg_coarse prefix.

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

169: .seealso: PCMGSetType(), PCMGGetLevels()
170: @*/
171: PetscErrorCode  PCMGSetLevels(PC pc,PetscInt levels,MPI_Comm *comms)
172: {
174:   PC_MG          *mg        = (PC_MG*)pc->data;
175:   MPI_Comm       comm;
176:   PC_MG_Levels   **mglevels = mg->levels;
177:   PetscInt       i;
178:   PetscMPIInt    size;
179:   const char     *prefix;
180:   PC             ipc;
181:   PetscInt       n;

186:   PetscObjectGetComm((PetscObject)pc,&comm);
187:   if (mg->nlevels == levels) return(0);
188:   if (mglevels) {
189:     /* changing the number of levels so free up the previous stuff */
190:     PCReset_MG(pc);
191:     n    = mglevels[0]->levels;
192:     for (i=0; i<n; i++) {
193:       if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
194:         KSPDestroy(&mglevels[i]->smoothd);
195:       }
196:       KSPDestroy(&mglevels[i]->smoothu);
197:       PetscFree(mglevels[i]);
198:     }
199:     PetscFree(mg->levels);
200:   }

202:   mg->nlevels = levels;

204:   PetscMalloc1(levels,&mglevels);
205:   PetscLogObjectMemory((PetscObject)pc,levels*(sizeof(PC_MG*)));

207:   PCGetOptionsPrefix(pc,&prefix);

209:   mg->stageApply = 0;
210:   for (i=0; i<levels; i++) {
211:     PetscNewLog(pc,&mglevels[i]);

213:     mglevels[i]->level               = i;
214:     mglevels[i]->levels              = levels;
215:     mglevels[i]->cycles              = PC_MG_CYCLE_V;
216:     mg->default_smoothu              = 2;
217:     mg->default_smoothd              = 2;
218:     mglevels[i]->eventsmoothsetup    = 0;
219:     mglevels[i]->eventsmoothsolve    = 0;
220:     mglevels[i]->eventresidual       = 0;
221:     mglevels[i]->eventinterprestrict = 0;

223:     if (comms) comm = comms[i];
224:     KSPCreate(comm,&mglevels[i]->smoothd);
225:     KSPSetType(mglevels[i]->smoothd,KSPCHEBYSHEV);
226:     KSPSetConvergenceTest(mglevels[i]->smoothd,KSPConvergedSkip,NULL,NULL);
227:     KSPSetNormType(mglevels[i]->smoothd,KSP_NORM_NONE);
228:     KSPGetPC(mglevels[i]->smoothd,&ipc);
229:     PCSetType(ipc,PCSOR);
230:     PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd,(PetscObject)pc,levels-i);
231:     KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, i ? mg->default_smoothd : 1);
232:     KSPSetOptionsPrefix(mglevels[i]->smoothd,prefix);

234:     /* do special stuff for coarse grid */
235:     if (!i && levels > 1) {
236:       KSPAppendOptionsPrefix(mglevels[0]->smoothd,"mg_coarse_");

238:       /* coarse solve is (redundant) LU by default; set shifttype NONZERO to avoid annoying zero-pivot in LU preconditioner */
239:       KSPSetType(mglevels[0]->smoothd,KSPPREONLY);
240:       KSPGetPC(mglevels[0]->smoothd,&ipc);
241:       MPI_Comm_size(comm,&size);
242:       if (size > 1) {
243:         KSP innerksp;
244:         PC  innerpc;
245:         PCSetType(ipc,PCREDUNDANT);
246:         PCRedundantGetKSP(ipc,&innerksp);
247:         KSPGetPC(innerksp,&innerpc);
248:         PCFactorSetShiftType(innerpc,MAT_SHIFT_INBLOCKS);
249:       } else {
250:         PCSetType(ipc,PCLU);
251:         PCFactorSetShiftType(ipc,MAT_SHIFT_INBLOCKS);
252:       }
253:     } else {
254:       char tprefix[128];
255:       sprintf(tprefix,"mg_levels_%d_",(int)i);
256:       KSPAppendOptionsPrefix(mglevels[i]->smoothd,tprefix);
257:     }
258:     PetscLogObjectParent((PetscObject)pc,(PetscObject)mglevels[i]->smoothd);

260:     mglevels[i]->smoothu = mglevels[i]->smoothd;
261:     mg->rtol             = 0.0;
262:     mg->abstol           = 0.0;
263:     mg->dtol             = 0.0;
264:     mg->ttol             = 0.0;
265:     mg->cyclesperpcapply = 1;
266:   }
267:   mg->am                   = PC_MG_MULTIPLICATIVE;
268:   mg->levels               = mglevels;
269:   pc->ops->applyrichardson = PCApplyRichardson_MG;
270:   return(0);
271: }


276: PetscErrorCode PCDestroy_MG(PC pc)
277: {
279:   PC_MG          *mg        = (PC_MG*)pc->data;
280:   PC_MG_Levels   **mglevels = mg->levels;
281:   PetscInt       i,n;

284:   PCReset_MG(pc);
285:   if (mglevels) {
286:     n = mglevels[0]->levels;
287:     for (i=0; i<n; i++) {
288:       if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
289:         KSPDestroy(&mglevels[i]->smoothd);
290:       }
291:       KSPDestroy(&mglevels[i]->smoothu);
292:       PetscFree(mglevels[i]);
293:     }
294:     PetscFree(mg->levels);
295:   }
296:   PetscFree(pc->data);
297:   return(0);
298: }



302: extern PetscErrorCode PCMGACycle_Private(PC,PC_MG_Levels**);
303: extern PetscErrorCode PCMGFCycle_Private(PC,PC_MG_Levels**);
304: extern PetscErrorCode PCMGKCycle_Private(PC,PC_MG_Levels**);

306: /*
307:    PCApply_MG - Runs either an additive, multiplicative, Kaskadic
308:              or full cycle of multigrid.

310:   Note:
311:   A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle().
312: */
315: static PetscErrorCode PCApply_MG(PC pc,Vec b,Vec x)
316: {
317:   PC_MG          *mg        = (PC_MG*)pc->data;
318:   PC_MG_Levels   **mglevels = mg->levels;
320:   PetscInt       levels = mglevels[0]->levels,i;

323:   if (mg->stageApply) {PetscLogStagePush(mg->stageApply);}
324:   /* When the DM is supplying the matrix then it will not exist until here */
325:   for (i=0; i<levels; i++) {
326:     if (!mglevels[i]->A) {
327:       KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);
328:       PetscObjectReference((PetscObject)mglevels[i]->A);
329:     }
330:   }

332:   mglevels[levels-1]->b = b;
333:   mglevels[levels-1]->x = x;
334:   if (mg->am == PC_MG_MULTIPLICATIVE) {
335:     VecSet(x,0.0);
336:     for (i=0; i<mg->cyclesperpcapply; i++) {
337:       PCMGMCycle_Private(pc,mglevels+levels-1,NULL);
338:     }
339:   } else if (mg->am == PC_MG_ADDITIVE) {
340:     PCMGACycle_Private(pc,mglevels);
341:   } else if (mg->am == PC_MG_KASKADE) {
342:     PCMGKCycle_Private(pc,mglevels);
343:   } else {
344:     PCMGFCycle_Private(pc,mglevels);
345:   }
346:   if (mg->stageApply) {PetscLogStagePop();}
347:   return(0);
348: }


353: PetscErrorCode PCSetFromOptions_MG(PC pc)
354: {
356:   PetscInt       m,levels = 1,cycles;
357:   PetscBool      flg,set;
358:   PC_MG          *mg        = (PC_MG*)pc->data;
359:   PC_MG_Levels   **mglevels = mg->levels;
360:   PCMGType       mgtype;
361:   PCMGCycleType  mgctype;

364:   PetscOptionsHead("Multigrid options");
365:   if (!mg->levels) {
366:     PetscOptionsInt("-pc_mg_levels","Number of Levels","PCMGSetLevels",levels,&levels,&flg);
367:     if (!flg && pc->dm) {
368:       DMGetRefineLevel(pc->dm,&levels);
369:       levels++;
370:       mg->usedmfornumberoflevels = PETSC_TRUE;
371:     }
372:     PCMGSetLevels(pc,levels,NULL);
373:   }
374:   mglevels = mg->levels;

376:   mgctype = (PCMGCycleType) mglevels[0]->cycles;
377:   PetscOptionsEnum("-pc_mg_cycle_type","V cycle or for W-cycle","PCMGSetCycleType",PCMGCycleTypes,(PetscEnum)mgctype,(PetscEnum*)&mgctype,&flg);
378:   if (flg) {
379:     PCMGSetCycleType(pc,mgctype);
380:   }
381:   flg  = PETSC_FALSE;
382:   PetscOptionsBool("-pc_mg_galerkin","Use Galerkin process to compute coarser operators","PCMGSetGalerkin",flg,&flg,&set);
383:   if (set) {
384:     PCMGSetGalerkin(pc,flg);
385:   }
386:   PetscOptionsInt("-pc_mg_smoothup","Number of post-smoothing steps","PCMGSetNumberSmoothUp",1,&m,&flg);
387:   if (flg) {
388:     PCMGSetNumberSmoothUp(pc,m);
389:   }
390:   PetscOptionsInt("-pc_mg_smoothdown","Number of pre-smoothing steps","PCMGSetNumberSmoothDown",1,&m,&flg);
391:   if (flg) {
392:     PCMGSetNumberSmoothDown(pc,m);
393:   }
394:   mgtype = mg->am;
395:   PetscOptionsEnum("-pc_mg_type","Multigrid type","PCMGSetType",PCMGTypes,(PetscEnum)mgtype,(PetscEnum*)&mgtype,&flg);
396:   if (flg) {
397:     PCMGSetType(pc,mgtype);
398:   }
399:   if (mg->am == PC_MG_MULTIPLICATIVE) {
400:     PetscOptionsInt("-pc_mg_multiplicative_cycles","Number of cycles for each preconditioner step","PCMGSetLevels",mg->cyclesperpcapply,&cycles,&flg);
401:     if (flg) {
402:       PCMGMultiplicativeSetCycles(pc,cycles);
403:     }
404:   }
405:   flg  = PETSC_FALSE;
406:   PetscOptionsBool("-pc_mg_log","Log times for each multigrid level","None",flg,&flg,NULL);
407:   if (flg) {
408:     PetscInt i;
409:     char     eventname[128];
410:     if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
411:     levels = mglevels[0]->levels;
412:     for (i=0; i<levels; i++) {
413:       sprintf(eventname,"MGSetup Level %d",(int)i);
414:       PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsetup);
415:       sprintf(eventname,"MGSmooth Level %d",(int)i);
416:       PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsolve);
417:       if (i) {
418:         sprintf(eventname,"MGResid Level %d",(int)i);
419:         PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventresidual);
420:         sprintf(eventname,"MGInterp Level %d",(int)i);
421:         PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventinterprestrict);
422:       }
423:     }

425: #if defined(PETSC_USE_LOG)
426:     {
427:       const char    *sname = "MG Apply";
428:       PetscStageLog stageLog;
429:       PetscInt      st;

432:       PetscLogGetStageLog(&stageLog);
433:       for (st = 0; st < stageLog->numStages; ++st) {
434:         PetscBool same;

436:         PetscStrcmp(stageLog->stageInfo[st].name, sname, &same);
437:         if (same) mg->stageApply = st;
438:       }
439:       if (!mg->stageApply) {
440:         PetscLogStageRegister(sname, &mg->stageApply);
441:       }
442:     }
443: #endif
444:   }
445:   PetscOptionsTail();
446:   return(0);
447: }

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

452: #include <petscdraw.h>
455: PetscErrorCode PCView_MG(PC pc,PetscViewer viewer)
456: {
457:   PC_MG          *mg        = (PC_MG*)pc->data;
458:   PC_MG_Levels   **mglevels = mg->levels;
460:   PetscInt       levels = mglevels ? mglevels[0]->levels : 0,i;
461:   PetscBool      iascii,isbinary,isdraw;

464:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
465:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
466:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
467:   if (iascii) {
468:     const char *cyclename = levels ? (mglevels[0]->cycles == PC_MG_CYCLE_V ? "v" : "w") : "unknown";
469:     PetscViewerASCIIPrintf(viewer,"  MG: type is %s, levels=%D cycles=%s\n", PCMGTypes[mg->am],levels,cyclename);
470:     if (mg->am == PC_MG_MULTIPLICATIVE) {
471:       PetscViewerASCIIPrintf(viewer,"    Cycles per PCApply=%d\n",mg->cyclesperpcapply);
472:     }
473:     if (mg->galerkin) {
474:       PetscViewerASCIIPrintf(viewer,"    Using Galerkin computed coarse grid matrices\n");
475:     } else {
476:       PetscViewerASCIIPrintf(viewer,"    Not using Galerkin computed coarse grid matrices\n");
477:     }
478:     for (i=0; i<levels; i++) {
479:       if (!i) {
480:         PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level -------------------------------\n",i);
481:       } else {
482:         PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);
483:       }
484:       PetscViewerASCIIPushTab(viewer);
485:       KSPView(mglevels[i]->smoothd,viewer);
486:       PetscViewerASCIIPopTab(viewer);
487:       if (i && mglevels[i]->smoothd == mglevels[i]->smoothu) {
488:         PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");
489:       } else if (i) {
490:         PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);
491:         PetscViewerASCIIPushTab(viewer);
492:         KSPView(mglevels[i]->smoothu,viewer);
493:         PetscViewerASCIIPopTab(viewer);
494:       }
495:     }
496:   } else if (isbinary) {
497:     for (i=levels-1; i>=0; i--) {
498:       KSPView(mglevels[i]->smoothd,viewer);
499:       if (i && mglevels[i]->smoothd != mglevels[i]->smoothu) {
500:         KSPView(mglevels[i]->smoothu,viewer);
501:       }
502:     }
503:   } else if (isdraw) {
504:     PetscDraw draw;
505:     PetscReal x,w,y,bottom,th;
506:     PetscViewerDrawGetDraw(viewer,0,&draw);
507:     PetscDrawGetCurrentPoint(draw,&x,&y);
508:     PetscDrawStringGetSize(draw,NULL,&th);
509:     bottom = y - th;
510:     for (i=levels-1; i>=0; i--) {
511:       if (!mglevels[i]->smoothu || (mglevels[i]->smoothu == mglevels[i]->smoothd)) {
512:         PetscDrawPushCurrentPoint(draw,x,bottom);
513:         KSPView(mglevels[i]->smoothd,viewer);
514:         PetscDrawPopCurrentPoint(draw);
515:       } else {
516:         w    = 0.5*PetscMin(1.0-x,x);
517:         PetscDrawPushCurrentPoint(draw,x+w,bottom);
518:         KSPView(mglevels[i]->smoothd,viewer);
519:         PetscDrawPopCurrentPoint(draw);
520:         PetscDrawPushCurrentPoint(draw,x-w,bottom);
521:         KSPView(mglevels[i]->smoothu,viewer);
522:         PetscDrawPopCurrentPoint(draw);
523:       }
524:       PetscDrawGetBoundingBox(draw,NULL,&bottom,NULL,NULL);
525:       bottom -= th;
526:     }
527:   }
528:   return(0);
529: }

531: #include <petsc-private/dmimpl.h>
532: #include <petsc-private/kspimpl.h>

534: /*
535:     Calls setup for the KSP on each level
536: */
539: PetscErrorCode PCSetUp_MG(PC pc)
540: {
541:   PC_MG          *mg        = (PC_MG*)pc->data;
542:   PC_MG_Levels   **mglevels = mg->levels;
544:   PetscInt       i,n = mglevels[0]->levels;
545:   PC             cpc;
546:   PetscBool      preonly,lu,redundant,cholesky,svd,dump = PETSC_FALSE,opsset,use_amat;
547:   Mat            dA,dB;
548:   Vec            tvec;
549:   DM             *dms;
550:   PetscViewer    viewer = 0;

553:   /* FIX: Move this to PCSetFromOptions_MG? */
554:   if (mg->usedmfornumberoflevels) {
555:     PetscInt levels;
556:     DMGetRefineLevel(pc->dm,&levels);
557:     levels++;
558:     if (levels > n) { /* the problem is now being solved on a finer grid */
559:       PCMGSetLevels(pc,levels,NULL);
560:       n        = levels;
561:       PCSetFromOptions(pc); /* it is bad to call this here, but otherwise will never be called for the new hierarchy */
562:       mglevels =  mg->levels;
563:     }
564:   }
565:   KSPGetPC(mglevels[0]->smoothd,&cpc);


568:   /* If user did not provide fine grid operators OR operator was not updated since last global KSPSetOperators() */
569:   /* so use those from global PC */
570:   /* Is this what we always want? What if user wants to keep old one? */
571:   KSPGetOperatorsSet(mglevels[n-1]->smoothd,NULL,&opsset);
572:   if (opsset) {
573:     Mat mmat;
574:     KSPGetOperators(mglevels[n-1]->smoothd,NULL,&mmat);
575:     if (mmat == pc->pmat) opsset = PETSC_FALSE;
576:   }

578:   if (!opsset) {
579:     PCGetUseAmat(pc,&use_amat);
580:     if(use_amat){
581:       PetscInfo(pc,"Using outer operators to define finest grid operator \n  because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
582:       KSPSetOperators(mglevels[n-1]->smoothd,pc->mat,pc->pmat);
583:     }
584:     else {
585:       PetscInfo(pc,"Using matrix (pmat) operators to define finest grid operator \n  because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
586:       KSPSetOperators(mglevels[n-1]->smoothd,pc->pmat,pc->pmat);
587:     }
588:   }

590:   /* Skipping this for galerkin==2 (externally managed hierarchy such as ML and GAMG). Cleaner logic here would be great. Wrap ML/GAMG as DMs? */
591:   if (pc->dm && mg->galerkin != 2 && !pc->setupcalled) {
592:     /* construct the interpolation from the DMs */
593:     Mat p;
594:     Vec rscale;
595:     PetscMalloc1(n,&dms);
596:     dms[n-1] = pc->dm;
597:     for (i=n-2; i>-1; i--) {
598:       DMKSP kdm;
599:       DMCoarsen(dms[i+1],MPI_COMM_NULL,&dms[i]);
600:       KSPSetDM(mglevels[i]->smoothd,dms[i]);
601:       if (mg->galerkin) {KSPSetDMActive(mglevels[i]->smoothd,PETSC_FALSE);}
602:       DMGetDMKSPWrite(dms[i],&kdm);
603:       /* Ugly hack so that the next KSPSetUp() will use the RHS that we set. A better fix is to change dmActive to take
604:        * a bitwise OR of computing the matrix, RHS, and initial iterate. */
605:       kdm->ops->computerhs = NULL;
606:       kdm->rhsctx          = NULL;
607:       if (!mglevels[i+1]->interpolate) {
608:         DMCreateInterpolation(dms[i],dms[i+1],&p,&rscale);
609:         PCMGSetInterpolation(pc,i+1,p);
610:         if (rscale) {PCMGSetRScale(pc,i+1,rscale);}
611:         VecDestroy(&rscale);
612:         MatDestroy(&p);
613:       }
614:     }

616:     for (i=n-2; i>-1; i--) {
617:       DMDestroy(&dms[i]);
618:     }
619:     PetscFree(dms);
620:   }

622:   if (pc->dm && !pc->setupcalled) {
623:     /* finest smoother also gets DM but it is not active, independent of whether galerkin==2 */
624:     KSPSetDM(mglevels[n-1]->smoothd,pc->dm);
625:     KSPSetDMActive(mglevels[n-1]->smoothd,PETSC_FALSE);
626:   }

628:   if (mg->galerkin == 1) {
629:     Mat B;
630:     /* currently only handle case where mat and pmat are the same on coarser levels */
631:     KSPGetOperators(mglevels[n-1]->smoothd,&dA,&dB);
632:     if (!pc->setupcalled) {
633:       for (i=n-2; i>-1; i--) {
634:         if (mglevels[i+1]->interpolate == mglevels[i+1]->restrct || !mglevels[i+1]->restrct) {
635:           MatPtAP(dB,mglevels[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);
636:         } else {
637:           MatMatMatMult(mglevels[i+1]->restrct,dB,mglevels[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);
638:         }
639:         KSPSetOperators(mglevels[i]->smoothd,B,B);
640:         if (i != n-2) {PetscObjectDereference((PetscObject)dB);}
641:         dB = B;
642:       }
643:       if (n > 1) {PetscObjectDereference((PetscObject)dB);}
644:     } else {
645:       for (i=n-2; i>-1; i--) {
646:         KSPGetOperators(mglevels[i]->smoothd,NULL,&B);
647:         if (mglevels[i+1]->interpolate == mglevels[i+1]->restrct || !mglevels[i+1]->restrct) {
648:           MatPtAP(dB,mglevels[i+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);
649:         } else {
650:           MatMatMatMult(mglevels[i+1]->restrct,dB,mglevels[i+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);
651:         }
652:         KSPSetOperators(mglevels[i]->smoothd,B,B);
653:         dB   = B;
654:       }
655:     }
656:   } else if (!mg->galerkin && pc->dm && pc->dm->x) {
657:     /* need to restrict Jacobian location to coarser meshes for evaluation */
658:     for (i=n-2; i>-1; i--) {
659:       Mat R;
660:       Vec rscale;
661:       if (!mglevels[i]->smoothd->dm->x) {
662:         Vec *vecs;
663:         KSPGetVecs(mglevels[i]->smoothd,1,&vecs,0,NULL);

665:         mglevels[i]->smoothd->dm->x = vecs[0];

667:         PetscFree(vecs);
668:       }
669:       PCMGGetRestriction(pc,i+1,&R);
670:       PCMGGetRScale(pc,i+1,&rscale);
671:       MatRestrict(R,mglevels[i+1]->smoothd->dm->x,mglevels[i]->smoothd->dm->x);
672:       VecPointwiseMult(mglevels[i]->smoothd->dm->x,mglevels[i]->smoothd->dm->x,rscale);
673:     }
674:   }
675:   if (!mg->galerkin && pc->dm) {
676:     for (i=n-2; i>=0; i--) {
677:       DM  dmfine,dmcoarse;
678:       Mat Restrict,Inject;
679:       Vec rscale;
680:       KSPGetDM(mglevels[i+1]->smoothd,&dmfine);
681:       KSPGetDM(mglevels[i]->smoothd,&dmcoarse);
682:       PCMGGetRestriction(pc,i+1,&Restrict);
683:       PCMGGetRScale(pc,i+1,&rscale);
684:       Inject = NULL;      /* Callback should create it if it needs Injection */
685:       DMRestrict(dmfine,Restrict,rscale,Inject,dmcoarse);
686:     }
687:   }

689:   if (!pc->setupcalled) {
690:     for (i=0; i<n; i++) {
691:       KSPSetFromOptions(mglevels[i]->smoothd);
692:     }
693:     for (i=1; i<n; i++) {
694:       if (mglevels[i]->smoothu && (mglevels[i]->smoothu != mglevels[i]->smoothd)) {
695:         KSPSetFromOptions(mglevels[i]->smoothu);
696:       }
697:     }
698:     for (i=1; i<n; i++) {
699:       PCMGGetInterpolation(pc,i,&mglevels[i]->interpolate);
700:       PCMGGetRestriction(pc,i,&mglevels[i]->restrct);
701:     }
702:     for (i=0; i<n-1; i++) {
703:       if (!mglevels[i]->b) {
704:         Vec *vec;
705:         KSPGetVecs(mglevels[i]->smoothd,1,&vec,0,NULL);
706:         PCMGSetRhs(pc,i,*vec);
707:         VecDestroy(vec);
708:         PetscFree(vec);
709:       }
710:       if (!mglevels[i]->r && i) {
711:         VecDuplicate(mglevels[i]->b,&tvec);
712:         PCMGSetR(pc,i,tvec);
713:         VecDestroy(&tvec);
714:       }
715:       if (!mglevels[i]->x) {
716:         VecDuplicate(mglevels[i]->b,&tvec);
717:         PCMGSetX(pc,i,tvec);
718:         VecDestroy(&tvec);
719:       }
720:     }
721:     if (n != 1 && !mglevels[n-1]->r) {
722:       /* PCMGSetR() on the finest level if user did not supply it */
723:       Vec *vec;
724:       KSPGetVecs(mglevels[n-1]->smoothd,1,&vec,0,NULL);
725:       PCMGSetR(pc,n-1,*vec);
726:       VecDestroy(vec);
727:       PetscFree(vec);
728:     }
729:   }

731:   if (pc->dm) {
732:     /* need to tell all the coarser levels to rebuild the matrix using the DM for that level */
733:     for (i=0; i<n-1; i++) {
734:       if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
735:     }
736:   }

738:   for (i=1; i<n; i++) {
739:     if (mglevels[i]->smoothu == mglevels[i]->smoothd) {
740:       /* if doing only down then initial guess is zero */
741:       KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);
742:     }
743:     if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
744:     KSPSetUp(mglevels[i]->smoothd);
745:     if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
746:     if (!mglevels[i]->residual) {
747:       Mat mat;
748:       KSPGetOperators(mglevels[i]->smoothd,NULL,&mat);
749:       PCMGSetResidual(pc,i,PCMGResidualDefault,mat);
750:     }
751:   }
752:   for (i=1; i<n; i++) {
753:     if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) {
754:       Mat          downmat,downpmat;

756:       /* check if operators have been set for up, if not use down operators to set them */
757:       KSPGetOperatorsSet(mglevels[i]->smoothu,&opsset,NULL);
758:       if (!opsset) {
759:         KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat);
760:         KSPSetOperators(mglevels[i]->smoothu,downmat,downpmat);
761:       }

763:       KSPSetInitialGuessNonzero(mglevels[i]->smoothu,PETSC_TRUE);
764:       if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
765:       KSPSetUp(mglevels[i]->smoothu);
766:       if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
767:     }
768:   }

770:   /*
771:       If coarse solver is not direct method then DO NOT USE preonly
772:   */
773:   PetscObjectTypeCompare((PetscObject)mglevels[0]->smoothd,KSPPREONLY,&preonly);
774:   if (preonly) {
775:     PetscObjectTypeCompare((PetscObject)cpc,PCLU,&lu);
776:     PetscObjectTypeCompare((PetscObject)cpc,PCREDUNDANT,&redundant);
777:     PetscObjectTypeCompare((PetscObject)cpc,PCCHOLESKY,&cholesky);
778:     PetscObjectTypeCompare((PetscObject)cpc,PCSVD,&svd);
779:     if (!lu && !redundant && !cholesky && !svd) {
780:       KSPSetType(mglevels[0]->smoothd,KSPGMRES);
781:     }
782:   }

784:   if (!pc->setupcalled) {
785:     KSPSetFromOptions(mglevels[0]->smoothd);
786:   }

788:   if (mglevels[0]->eventsmoothsetup) {PetscLogEventBegin(mglevels[0]->eventsmoothsetup,0,0,0,0);}
789:   KSPSetUp(mglevels[0]->smoothd);
790:   if (mglevels[0]->eventsmoothsetup) {PetscLogEventEnd(mglevels[0]->eventsmoothsetup,0,0,0,0);}

792:   /*
793:      Dump the interpolation/restriction matrices plus the
794:    Jacobian/stiffness on each level. This allows MATLAB users to
795:    easily check if the Galerkin condition A_c = R A_f R^T is satisfied.

797:    Only support one or the other at the same time.
798:   */
799: #if defined(PETSC_USE_SOCKET_VIEWER)
800:   PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_mg_dump_matlab",&dump,NULL);
801:   if (dump) viewer = PETSC_VIEWER_SOCKET_(PetscObjectComm((PetscObject)pc));
802:   dump = PETSC_FALSE;
803: #endif
804:   PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_mg_dump_binary",&dump,NULL);
805:   if (dump) viewer = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)pc));

807:   if (viewer) {
808:     for (i=1; i<n; i++) {
809:       MatView(mglevels[i]->restrct,viewer);
810:     }
811:     for (i=0; i<n; i++) {
812:       KSPGetPC(mglevels[i]->smoothd,&pc);
813:       MatView(pc->mat,viewer);
814:     }
815:   }
816:   return(0);
817: }

819: /* -------------------------------------------------------------------------------------*/

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

826:    Not Collective

828:    Input Parameter:
829: .  pc - the preconditioner context

831:    Output parameter:
832: .  levels - the number of levels

834:    Level: advanced

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

838: .seealso: PCMGSetLevels()
839: @*/
840: PetscErrorCode  PCMGGetLevels(PC pc,PetscInt *levels)
841: {
842:   PC_MG *mg = (PC_MG*)pc->data;

847:   *levels = mg->nlevels;
848:   return(0);
849: }

853: /*@
854:    PCMGSetType - Determines the form of multigrid to use:
855:    multiplicative, additive, full, or the Kaskade algorithm.

857:    Logically Collective on PC

859:    Input Parameters:
860: +  pc - the preconditioner context
861: -  form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,
862:    PC_MG_FULL, PC_MG_KASKADE

864:    Options Database Key:
865: .  -pc_mg_type <form> - Sets <form>, one of multiplicative,
866:    additive, full, kaskade

868:    Level: advanced

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

872: .seealso: PCMGSetLevels()
873: @*/
874: PetscErrorCode  PCMGSetType(PC pc,PCMGType form)
875: {
876:   PC_MG *mg = (PC_MG*)pc->data;

881:   mg->am = form;
882:   if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG;
883:   else pc->ops->applyrichardson = 0;
884:   return(0);
885: }

889: /*@
890:    PCMGSetCycleType - Sets the type cycles to use.  Use PCMGSetCycleTypeOnLevel() for more
891:    complicated cycling.

893:    Logically Collective on PC

895:    Input Parameters:
896: +  pc - the multigrid context
897: -  PC_MG_CYCLE_V or PC_MG_CYCLE_W

899:    Options Database Key:
900: $  -pc_mg_cycle_type v or w

902:    Level: advanced

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

906: .seealso: PCMGSetCycleTypeOnLevel()
907: @*/
908: PetscErrorCode  PCMGSetCycleType(PC pc,PCMGCycleType n)
909: {
910:   PC_MG        *mg        = (PC_MG*)pc->data;
911:   PC_MG_Levels **mglevels = mg->levels;
912:   PetscInt     i,levels;

916:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
918:   levels = mglevels[0]->levels;

920:   for (i=0; i<levels; i++) mglevels[i]->cycles = n;
921:   return(0);
922: }

926: /*@
927:    PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step
928:          of multigrid when PCMGType of PC_MG_MULTIPLICATIVE is used

930:    Logically Collective on PC

932:    Input Parameters:
933: +  pc - the multigrid context
934: -  n - number of cycles (default is 1)

936:    Options Database Key:
937: $  -pc_mg_multiplicative_cycles n

939:    Level: advanced

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

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

945: .seealso: PCMGSetCycleTypeOnLevel(), PCMGSetCycleType()
946: @*/
947: PetscErrorCode  PCMGMultiplicativeSetCycles(PC pc,PetscInt n)
948: {
949:   PC_MG        *mg        = (PC_MG*)pc->data;
950:   PC_MG_Levels **mglevels = mg->levels;
951:   PetscInt     i,levels;

955:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
957:   levels = mglevels[0]->levels;

959:   for (i=0; i<levels; i++) mg->cyclesperpcapply = n;
960:   return(0);
961: }

965: /*@
966:    PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the
967:       finest grid via the Galerkin process: A_i-1 = r_i * A_i * p_i

969:    Logically Collective on PC

971:    Input Parameters:
972: +  pc - the multigrid context
973: -  use - PETSC_TRUE to use the Galerkin process to compute coarse-level operators

975:    Options Database Key:
976: $  -pc_mg_galerkin

978:    Level: intermediate

980: .keywords: MG, set, Galerkin

982: .seealso: PCMGGetGalerkin()

984: @*/
985: PetscErrorCode PCMGSetGalerkin(PC pc,PetscBool use)
986: {
987:   PC_MG *mg = (PC_MG*)pc->data;

991:   mg->galerkin = use ? 1 : 0;
992:   return(0);
993: }

997: /*@
998:    PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e.
999:       A_i-1 = r_i * A_i * p_i

1001:    Not Collective

1003:    Input Parameter:
1004: .  pc - the multigrid context

1006:    Output Parameter:
1007: .  gelerkin - PETSC_TRUE or PETSC_FALSE

1009:    Options Database Key:
1010: $  -pc_mg_galerkin

1012:    Level: intermediate

1014: .keywords: MG, set, Galerkin

1016: .seealso: PCMGSetGalerkin()

1018: @*/
1019: PetscErrorCode  PCMGGetGalerkin(PC pc,PetscBool  *galerkin)
1020: {
1021:   PC_MG *mg = (PC_MG*)pc->data;

1025:   *galerkin = (PetscBool)mg->galerkin;
1026:   return(0);
1027: }

1031: /*@
1032:    PCMGSetNumberSmoothDown - Sets the number of pre-smoothing steps to
1033:    use on all levels. Use PCMGGetSmootherDown() to set different
1034:    pre-smoothing steps on different levels.

1036:    Logically Collective on PC

1038:    Input Parameters:
1039: +  mg - the multigrid context
1040: -  n - the number of smoothing steps

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

1045:    Level: advanced

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

1049: .seealso: PCMGSetNumberSmoothUp()
1050: @*/
1051: PetscErrorCode  PCMGSetNumberSmoothDown(PC pc,PetscInt n)
1052: {
1053:   PC_MG          *mg        = (PC_MG*)pc->data;
1054:   PC_MG_Levels   **mglevels = mg->levels;
1056:   PetscInt       i,levels;

1060:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
1062:   levels = mglevels[0]->levels;

1064:   for (i=1; i<levels; i++) {
1065:     /* make sure smoother up and down are different */
1066:     PCMGGetSmootherUp(pc,i,NULL);
1067:     KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);

1069:     mg->default_smoothd = n;
1070:   }
1071:   return(0);
1072: }

1076: /*@
1077:    PCMGSetNumberSmoothUp - Sets the number of post-smoothing steps to use
1078:    on all levels. Use PCMGGetSmootherUp() to set different numbers of
1079:    post-smoothing steps on different levels.

1081:    Logically Collective on PC

1083:    Input Parameters:
1084: +  mg - the multigrid context
1085: -  n - the number of smoothing steps

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

1090:    Level: advanced

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

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

1097: .seealso: PCMGSetNumberSmoothDown()
1098: @*/
1099: PetscErrorCode  PCMGSetNumberSmoothUp(PC pc,PetscInt n)
1100: {
1101:   PC_MG          *mg        = (PC_MG*)pc->data;
1102:   PC_MG_Levels   **mglevels = mg->levels;
1104:   PetscInt       i,levels;

1108:   if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
1110:   levels = mglevels[0]->levels;

1112:   for (i=1; i<levels; i++) {
1113:     /* make sure smoother up and down are different */
1114:     PCMGGetSmootherUp(pc,i,NULL);
1115:     KSPSetTolerances(mglevels[i]->smoothu,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);

1117:     mg->default_smoothu = n;
1118:   }
1119:   return(0);
1120: }

1122: /* ----------------------------------------------------------------------------------------*/

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

1128:    Options Database Keys:
1129: +  -pc_mg_levels <nlevels> - number of levels including finest
1130: .  -pc_mg_cycles v or w
1131: .  -pc_mg_smoothup <n> - number of smoothing steps after interpolation
1132: .  -pc_mg_smoothdown <n> - number of smoothing steps before applying restriction operator
1133: .  -pc_mg_type <additive,multiplicative,full,kaskade> - multiplicative is the default
1134: .  -pc_mg_log - log information about time spent on each level of the solver
1135: .  -pc_mg_monitor - print information on the multigrid convergence
1136: .  -pc_mg_galerkin - use Galerkin process to compute coarser operators, i.e. Acoarse = R A R'
1137: .  -pc_mg_multiplicative_cycles - number of cycles to use as the preconditioner (defaults to 1)
1138: .  -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices
1139:                         to the Socket viewer for reading from MATLAB.
1140: -  -pc_mg_dump_binary - dumps the matrices for each level and the restriction/interpolation matrices
1141:                         to the binary output file called binaryoutput

1143:    Notes: By default this uses GMRES on the fine grid smoother so this should be used with KSPFGMRES or the smoother changed to not use GMRES

1145:        When run with a single level the smoother options are used on that level NOT the coarse grid solver options

1147:    Level: intermediate

1149:    Concepts: multigrid/multilevel

1151: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, PCEXOTIC, PCGAMG, PCML, PCHYPRE
1152:            PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), PCMGSetNumberSmoothDown(),
1153:            PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(),
1154:            PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(),
1155:            PCMGSetCycleTypeOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR()
1156: M*/

1160: PETSC_EXTERN PetscErrorCode PCCreate_MG(PC pc)
1161: {
1162:   PC_MG          *mg;

1166:   PetscNewLog(pc,&mg);
1167:   pc->data    = (void*)mg;
1168:   mg->nlevels = -1;

1170:   pc->useAmat = PETSC_TRUE;

1172:   pc->ops->apply          = PCApply_MG;
1173:   pc->ops->setup          = PCSetUp_MG;
1174:   pc->ops->reset          = PCReset_MG;
1175:   pc->ops->destroy        = PCDestroy_MG;
1176:   pc->ops->setfromoptions = PCSetFromOptions_MG;
1177:   pc->ops->view           = PCView_MG;
1178:   return(0);
1179: }