Actual source code: mg.c

petsc-dev 2014-07-22
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:   mglevels[levels-1]->b = b;
 71:   mglevels[levels-1]->x = x;

 73:   mg->rtol   = rtol;
 74:   mg->abstol = abstol;
 75:   mg->dtol   = dtol;
 76:   if (rtol) {
 77:     /* compute initial residual norm for relative convergence test */
 78:     PetscReal rnorm;
 79:     if (zeroguess) {
 80:       VecNorm(b,NORM_2,&rnorm);
 81:     } else {
 82:       (*mglevels[levels-1]->residual)(mglevels[levels-1]->A,b,x,w);
 83:       VecNorm(w,NORM_2,&rnorm);
 84:     }
 85:     mg->ttol = PetscMax(rtol*rnorm,abstol);
 86:   } else if (abstol) mg->ttol = abstol;
 87:   else mg->ttol = 0.0;

 89:   /* since smoother is applied to full system, not just residual we need to make sure that smoothers don't
 90:      stop prematurely due to small residual */
 91:   for (i=1; i<levels; i++) {
 92:     KSPSetTolerances(mglevels[i]->smoothu,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
 93:     if (mglevels[i]->smoothu != mglevels[i]->smoothd) {
 94:       KSPSetTolerances(mglevels[i]->smoothd,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
 95:     }
 96:   }

 98:   *reason = (PCRichardsonConvergedReason)0;
 99:   for (i=0; i<its; i++) {
100:     PCMGMCycle_Private(pc,mglevels+levels-1,reason);
101:     if (*reason) break;
102:   }
103:   if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
104:   *outits = i;
105:   return(0);
106: }

110: PetscErrorCode PCReset_MG(PC pc)
111: {
112:   PC_MG          *mg        = (PC_MG*)pc->data;
113:   PC_MG_Levels   **mglevels = mg->levels;
115:   PetscInt       i,n;

118:   if (mglevels) {
119:     n = mglevels[0]->levels;
120:     for (i=0; i<n-1; i++) {
121:       VecDestroy(&mglevels[i+1]->r);
122:       VecDestroy(&mglevels[i]->b);
123:       VecDestroy(&mglevels[i]->x);
124:       MatDestroy(&mglevels[i+1]->restrct);
125:       MatDestroy(&mglevels[i+1]->interpolate);
126:       VecDestroy(&mglevels[i+1]->rscale);
127:     }

129:     for (i=0; i<n; i++) {
130:       MatDestroy(&mglevels[i]->A);
131:       if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
132:         KSPReset(mglevels[i]->smoothd);
133:       }
134:       KSPReset(mglevels[i]->smoothu);
135:     }
136:   }
137:   return(0);
138: }

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

146:    Logically Collective on PC

148:    Input Parameters:
149: +  pc - the preconditioner context
150: .  levels - the number of levels
151: -  comms - optional communicators for each level; this is to allow solving the coarser problems
152:            on smaller sets of processors. Use NULL_OBJECT for default in Fortran

154:    Level: intermediate

156:    Notes:
157:      If the number of levels is one then the multigrid uses the -mg_levels prefix
158:   for setting the level options rather than the -mg_coarse prefix.

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

162: .seealso: PCMGSetType(), PCMGGetLevels()
163: @*/
164: PetscErrorCode  PCMGSetLevels(PC pc,PetscInt levels,MPI_Comm *comms)
165: {
167:   PC_MG          *mg        = (PC_MG*)pc->data;
168:   MPI_Comm       comm;
169:   PC_MG_Levels   **mglevels = mg->levels;
170:   PetscInt       i;
171:   PetscMPIInt    size;
172:   const char     *prefix;
173:   PC             ipc;
174:   PetscInt       n;

179:   PetscObjectGetComm((PetscObject)pc,&comm);
180:   if (mg->nlevels == levels) return(0);
181:   if (mglevels) {
182:     /* changing the number of levels so free up the previous stuff */
183:     PCReset_MG(pc);
184:     n    = mglevels[0]->levels;
185:     for (i=0; i<n; i++) {
186:       if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
187:         KSPDestroy(&mglevels[i]->smoothd);
188:       }
189:       KSPDestroy(&mglevels[i]->smoothu);
190:       PetscFree(mglevels[i]);
191:     }
192:     PetscFree(mg->levels);
193:   }

195:   mg->nlevels = levels;

197:   PetscMalloc1(levels,&mglevels);
198:   PetscLogObjectMemory((PetscObject)pc,levels*(sizeof(PC_MG*)));

200:   PCGetOptionsPrefix(pc,&prefix);

202:   mg->stageApply = 0;
203:   for (i=0; i<levels; i++) {
204:     PetscNewLog(pc,&mglevels[i]);

206:     mglevels[i]->level               = i;
207:     mglevels[i]->levels              = levels;
208:     mglevels[i]->cycles              = PC_MG_CYCLE_V;
209:     mg->default_smoothu              = 2;
210:     mg->default_smoothd              = 2;
211:     mglevels[i]->eventsmoothsetup    = 0;
212:     mglevels[i]->eventsmoothsolve    = 0;
213:     mglevels[i]->eventresidual       = 0;
214:     mglevels[i]->eventinterprestrict = 0;

216:     if (comms) comm = comms[i];
217:     KSPCreate(comm,&mglevels[i]->smoothd);
218:     KSPSetType(mglevels[i]->smoothd,KSPCHEBYSHEV);
219:     KSPSetConvergenceTest(mglevels[i]->smoothd,KSPConvergedSkip,NULL,NULL);
220:     KSPSetNormType(mglevels[i]->smoothd,KSP_NORM_NONE);
221:     KSPGetPC(mglevels[i]->smoothd,&ipc);
222:     PCSetType(ipc,PCSOR);
223:     PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd,(PetscObject)pc,levels-i);
224:     KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, i ? mg->default_smoothd : 1);
225:     KSPSetOptionsPrefix(mglevels[i]->smoothd,prefix);

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

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

253:     mglevels[i]->smoothu = mglevels[i]->smoothd;
254:     mg->rtol             = 0.0;
255:     mg->abstol           = 0.0;
256:     mg->dtol             = 0.0;
257:     mg->ttol             = 0.0;
258:     mg->cyclesperpcapply = 1;
259:   }
260:   mg->am                   = PC_MG_MULTIPLICATIVE;
261:   mg->levels               = mglevels;
262:   pc->ops->applyrichardson = PCApplyRichardson_MG;
263:   return(0);
264: }


269: PetscErrorCode PCDestroy_MG(PC pc)
270: {
272:   PC_MG          *mg        = (PC_MG*)pc->data;
273:   PC_MG_Levels   **mglevels = mg->levels;
274:   PetscInt       i,n;

277:   PCReset_MG(pc);
278:   if (mglevels) {
279:     n = mglevels[0]->levels;
280:     for (i=0; i<n; i++) {
281:       if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
282:         KSPDestroy(&mglevels[i]->smoothd);
283:       }
284:       KSPDestroy(&mglevels[i]->smoothu);
285:       PetscFree(mglevels[i]);
286:     }
287:     PetscFree(mg->levels);
288:   }
289:   PetscFree(pc->data);
290:   return(0);
291: }



295: extern PetscErrorCode PCMGACycle_Private(PC,PC_MG_Levels**);
296: extern PetscErrorCode PCMGFCycle_Private(PC,PC_MG_Levels**);
297: extern PetscErrorCode PCMGKCycle_Private(PC,PC_MG_Levels**);

299: /*
300:    PCApply_MG - Runs either an additive, multiplicative, Kaskadic
301:              or full cycle of multigrid.

303:   Note:
304:   A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle().
305: */
308: static PetscErrorCode PCApply_MG(PC pc,Vec b,Vec x)
309: {
310:   PC_MG          *mg        = (PC_MG*)pc->data;
311:   PC_MG_Levels   **mglevels = mg->levels;
313:   PetscInt       levels = mglevels[0]->levels,i;

316:   if (mg->stageApply) {PetscLogStagePush(mg->stageApply);}
317:   /* When the DM is supplying the matrix then it will not exist until here */
318:   for (i=0; i<levels; i++) {
319:     if (!mglevels[i]->A) {
320:       KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);
321:       PetscObjectReference((PetscObject)mglevels[i]->A);
322:     }
323:   }

325:   mglevels[levels-1]->b = b;
326:   mglevels[levels-1]->x = x;
327:   if (mg->am == PC_MG_MULTIPLICATIVE) {
328:     VecSet(x,0.0);
329:     for (i=0; i<mg->cyclesperpcapply; i++) {
330:       PCMGMCycle_Private(pc,mglevels+levels-1,NULL);
331:     }
332:   } else if (mg->am == PC_MG_ADDITIVE) {
333:     PCMGACycle_Private(pc,mglevels);
334:   } else if (mg->am == PC_MG_KASKADE) {
335:     PCMGKCycle_Private(pc,mglevels);
336:   } else {
337:     PCMGFCycle_Private(pc,mglevels);
338:   }
339:   if (mg->stageApply) {PetscLogStagePop();}
340:   return(0);
341: }


346: PetscErrorCode PCSetFromOptions_MG(PC pc)
347: {
349:   PetscInt       m,levels = 1,cycles;
350:   PetscBool      flg,set;
351:   PC_MG          *mg        = (PC_MG*)pc->data;
352:   PC_MG_Levels   **mglevels = mg->levels;
353:   PCMGType       mgtype;
354:   PCMGCycleType  mgctype;

357:   PetscOptionsHead("Multigrid options");
358:   if (!mg->levels) {
359:     PetscOptionsInt("-pc_mg_levels","Number of Levels","PCMGSetLevels",levels,&levels,&flg);
360:     if (!flg && pc->dm) {
361:       DMGetRefineLevel(pc->dm,&levels);
362:       levels++;
363:       mg->usedmfornumberoflevels = PETSC_TRUE;
364:     }
365:     PCMGSetLevels(pc,levels,NULL);
366:   }
367:   mglevels = mg->levels;

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

418: #if defined(PETSC_USE_LOG)
419:     {
420:       const char    *sname = "MG Apply";
421:       PetscStageLog stageLog;
422:       PetscInt      st;

425:       PetscLogGetStageLog(&stageLog);
426:       for (st = 0; st < stageLog->numStages; ++st) {
427:         PetscBool same;

429:         PetscStrcmp(stageLog->stageInfo[st].name, sname, &same);
430:         if (same) mg->stageApply = st;
431:       }
432:       if (!mg->stageApply) {
433:         PetscLogStageRegister(sname, &mg->stageApply);
434:       }
435:     }
436: #endif
437:   }
438:   PetscOptionsTail();
439:   return(0);
440: }

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

445: #include <petscdraw.h>
448: PetscErrorCode PCView_MG(PC pc,PetscViewer viewer)
449: {
450:   PC_MG          *mg        = (PC_MG*)pc->data;
451:   PC_MG_Levels   **mglevels = mg->levels;
453:   PetscInt       levels = mglevels ? mglevels[0]->levels : 0,i;
454:   PetscBool      iascii,isbinary,isdraw;

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

524: #include <petsc-private/dmimpl.h>
525: #include <petsc-private/kspimpl.h>

527: /*
528:     Calls setup for the KSP on each level
529: */
532: PetscErrorCode PCSetUp_MG(PC pc)
533: {
534:   PC_MG          *mg        = (PC_MG*)pc->data;
535:   PC_MG_Levels   **mglevels = mg->levels;
537:   PetscInt       i,n = mglevels[0]->levels;
538:   PC             cpc;
539:   PetscBool      preonly,lu,redundant,cholesky,svd,dump = PETSC_FALSE,opsset,use_amat;
540:   Mat            dA,dB;
541:   Vec            tvec;
542:   DM             *dms;
543:   PetscViewer    viewer = 0;

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


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

571:   if (!opsset) {
572:     PCGetUseAmat(pc,&use_amat);
573:     if(use_amat){
574:       PetscInfo(pc,"Using outer operators to define finest grid operator \n  because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
575:       KSPSetOperators(mglevels[n-1]->smoothd,pc->mat,pc->pmat);
576:     }
577:     else {
578:       PetscInfo(pc,"Using matrix (pmat) operators to define finest grid operator \n  because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
579:       KSPSetOperators(mglevels[n-1]->smoothd,pc->pmat,pc->pmat);
580:     }
581:   }

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

609:     for (i=n-2; i>-1; i--) {
610:       DMDestroy(&dms[i]);
611:     }
612:     PetscFree(dms);
613:   }

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

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

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

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

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

724:   if (pc->dm) {
725:     /* need to tell all the coarser levels to rebuild the matrix using the DM for that level */
726:     for (i=0; i<n-1; i++) {
727:       if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
728:     }
729:   }

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

749:       /* check if operators have been set for up, if not use down operators to set them */
750:       KSPGetOperatorsSet(mglevels[i]->smoothu,&opsset,NULL);
751:       if (!opsset) {
752:         KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat);
753:         KSPSetOperators(mglevels[i]->smoothu,downmat,downpmat);
754:       }

756:       KSPSetInitialGuessNonzero(mglevels[i]->smoothu,PETSC_TRUE);
757:       if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
758:       KSPSetUp(mglevels[i]->smoothu);
759:       if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
760:     }
761:   }

763:   /*
764:       If coarse solver is not direct method then DO NOT USE preonly
765:   */
766:   PetscObjectTypeCompare((PetscObject)mglevels[0]->smoothd,KSPPREONLY,&preonly);
767:   if (preonly) {
768:     PetscObjectTypeCompare((PetscObject)cpc,PCLU,&lu);
769:     PetscObjectTypeCompare((PetscObject)cpc,PCREDUNDANT,&redundant);
770:     PetscObjectTypeCompare((PetscObject)cpc,PCCHOLESKY,&cholesky);
771:     PetscObjectTypeCompare((PetscObject)cpc,PCSVD,&svd);
772:     if (!lu && !redundant && !cholesky && !svd) {
773:       KSPSetType(mglevels[0]->smoothd,KSPGMRES);
774:     }
775:   }

777:   if (!pc->setupcalled) {
778:     KSPSetFromOptions(mglevels[0]->smoothd);
779:   }

781:   if (mglevels[0]->eventsmoothsetup) {PetscLogEventBegin(mglevels[0]->eventsmoothsetup,0,0,0,0);}
782:   KSPSetUp(mglevels[0]->smoothd);
783:   if (mglevels[0]->eventsmoothsetup) {PetscLogEventEnd(mglevels[0]->eventsmoothsetup,0,0,0,0);}

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

790:    Only support one or the other at the same time.
791:   */
792: #if defined(PETSC_USE_SOCKET_VIEWER)
793:   PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_mg_dump_matlab",&dump,NULL);
794:   if (dump) viewer = PETSC_VIEWER_SOCKET_(PetscObjectComm((PetscObject)pc));
795:   dump = PETSC_FALSE;
796: #endif
797:   PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_mg_dump_binary",&dump,NULL);
798:   if (dump) viewer = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)pc));

800:   if (viewer) {
801:     for (i=1; i<n; i++) {
802:       MatView(mglevels[i]->restrct,viewer);
803:     }
804:     for (i=0; i<n; i++) {
805:       KSPGetPC(mglevels[i]->smoothd,&pc);
806:       MatView(pc->mat,viewer);
807:     }
808:   }
809:   return(0);
810: }

812: /* -------------------------------------------------------------------------------------*/

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

819:    Not Collective

821:    Input Parameter:
822: .  pc - the preconditioner context

824:    Output parameter:
825: .  levels - the number of levels

827:    Level: advanced

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

831: .seealso: PCMGSetLevels()
832: @*/
833: PetscErrorCode  PCMGGetLevels(PC pc,PetscInt *levels)
834: {
835:   PC_MG *mg = (PC_MG*)pc->data;

840:   *levels = mg->nlevels;
841:   return(0);
842: }

846: /*@
847:    PCMGSetType - Determines the form of multigrid to use:
848:    multiplicative, additive, full, or the Kaskade algorithm.

850:    Logically Collective on PC

852:    Input Parameters:
853: +  pc - the preconditioner context
854: -  form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,
855:    PC_MG_FULL, PC_MG_KASKADE

857:    Options Database Key:
858: .  -pc_mg_type <form> - Sets <form>, one of multiplicative,
859:    additive, full, kaskade

861:    Level: advanced

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

865: .seealso: PCMGSetLevels()
866: @*/
867: PetscErrorCode  PCMGSetType(PC pc,PCMGType form)
868: {
869:   PC_MG *mg = (PC_MG*)pc->data;

874:   mg->am = form;
875:   if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG;
876:   else pc->ops->applyrichardson = 0;
877:   return(0);
878: }

882: /*@
883:    PCMGSetCycleType - Sets the type cycles to use.  Use PCMGSetCycleTypeOnLevel() for more
884:    complicated cycling.

886:    Logically Collective on PC

888:    Input Parameters:
889: +  pc - the multigrid context
890: -  PC_MG_CYCLE_V or PC_MG_CYCLE_W

892:    Options Database Key:
893: $  -pc_mg_cycle_type v or w

895:    Level: advanced

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

899: .seealso: PCMGSetCycleTypeOnLevel()
900: @*/
901: PetscErrorCode  PCMGSetCycleType(PC pc,PCMGCycleType n)
902: {
903:   PC_MG        *mg        = (PC_MG*)pc->data;
904:   PC_MG_Levels **mglevels = mg->levels;
905:   PetscInt     i,levels;

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

913:   for (i=0; i<levels; i++) mglevels[i]->cycles = n;
914:   return(0);
915: }

919: /*@
920:    PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step
921:          of multigrid when PCMGType of PC_MG_MULTIPLICATIVE is used

923:    Logically Collective on PC

925:    Input Parameters:
926: +  pc - the multigrid context
927: -  n - number of cycles (default is 1)

929:    Options Database Key:
930: $  -pc_mg_multiplicative_cycles n

932:    Level: advanced

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

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

938: .seealso: PCMGSetCycleTypeOnLevel(), PCMGSetCycleType()
939: @*/
940: PetscErrorCode  PCMGMultiplicativeSetCycles(PC pc,PetscInt n)
941: {
942:   PC_MG        *mg        = (PC_MG*)pc->data;
943:   PC_MG_Levels **mglevels = mg->levels;
944:   PetscInt     i,levels;

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

952:   for (i=0; i<levels; i++) mg->cyclesperpcapply = n;
953:   return(0);
954: }

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

962:    Logically Collective on PC

964:    Input Parameters:
965: +  pc - the multigrid context
966: -  use - PETSC_TRUE to use the Galerkin process to compute coarse-level operators

968:    Options Database Key:
969: $  -pc_mg_galerkin

971:    Level: intermediate

973: .keywords: MG, set, Galerkin

975: .seealso: PCMGGetGalerkin()

977: @*/
978: PetscErrorCode PCMGSetGalerkin(PC pc,PetscBool use)
979: {
980:   PC_MG *mg = (PC_MG*)pc->data;

984:   mg->galerkin = use ? 1 : 0;
985:   return(0);
986: }

990: /*@
991:    PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e.
992:       A_i-1 = r_i * A_i * p_i

994:    Not Collective

996:    Input Parameter:
997: .  pc - the multigrid context

999:    Output Parameter:
1000: .  gelerkin - PETSC_TRUE or PETSC_FALSE

1002:    Options Database Key:
1003: $  -pc_mg_galerkin

1005:    Level: intermediate

1007: .keywords: MG, set, Galerkin

1009: .seealso: PCMGSetGalerkin()

1011: @*/
1012: PetscErrorCode  PCMGGetGalerkin(PC pc,PetscBool  *galerkin)
1013: {
1014:   PC_MG *mg = (PC_MG*)pc->data;

1018:   *galerkin = (PetscBool)mg->galerkin;
1019:   return(0);
1020: }

1024: /*@
1025:    PCMGSetNumberSmoothDown - Sets the number of pre-smoothing steps to
1026:    use on all levels. Use PCMGGetSmootherDown() to set different
1027:    pre-smoothing steps on different levels.

1029:    Logically Collective on PC

1031:    Input Parameters:
1032: +  mg - the multigrid context
1033: -  n - the number of smoothing steps

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

1038:    Level: advanced

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

1042: .seealso: PCMGSetNumberSmoothUp()
1043: @*/
1044: PetscErrorCode  PCMGSetNumberSmoothDown(PC pc,PetscInt n)
1045: {
1046:   PC_MG          *mg        = (PC_MG*)pc->data;
1047:   PC_MG_Levels   **mglevels = mg->levels;
1049:   PetscInt       i,levels;

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

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

1062:     mg->default_smoothd = n;
1063:   }
1064:   return(0);
1065: }

1069: /*@
1070:    PCMGSetNumberSmoothUp - Sets the number of post-smoothing steps to use
1071:    on all levels. Use PCMGGetSmootherUp() to set different numbers of
1072:    post-smoothing steps on different levels.

1074:    Logically Collective on PC

1076:    Input Parameters:
1077: +  mg - the multigrid context
1078: -  n - the number of smoothing steps

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

1083:    Level: advanced

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

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

1090: .seealso: PCMGSetNumberSmoothDown()
1091: @*/
1092: PetscErrorCode  PCMGSetNumberSmoothUp(PC pc,PetscInt n)
1093: {
1094:   PC_MG          *mg        = (PC_MG*)pc->data;
1095:   PC_MG_Levels   **mglevels = mg->levels;
1097:   PetscInt       i,levels;

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

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

1110:     mg->default_smoothu = n;
1111:   }
1112:   return(0);
1113: }

1115: /* ----------------------------------------------------------------------------------------*/

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

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

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

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

1140:    Level: intermediate

1142:    Concepts: multigrid/multilevel

1144: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, PCEXOTIC, PCGAMG, PCML, PCHYPRE
1145:            PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), PCMGSetNumberSmoothDown(),
1146:            PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(),
1147:            PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(),
1148:            PCMGSetCycleTypeOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR()
1149: M*/

1153: PETSC_EXTERN PetscErrorCode PCCreate_MG(PC pc)
1154: {
1155:   PC_MG          *mg;

1159:   PetscNewLog(pc,&mg);
1160:   pc->data    = (void*)mg;
1161:   mg->nlevels = -1;

1163:   pc->useAmat = PETSC_TRUE;

1165:   pc->ops->apply          = PCApply_MG;
1166:   pc->ops->setup          = PCSetUp_MG;
1167:   pc->ops->reset          = PCReset_MG;
1168:   pc->ops->destroy        = PCDestroy_MG;
1169:   pc->ops->setfromoptions = PCSetFromOptions_MG;
1170:   pc->ops->view           = PCView_MG;
1171:   return(0);
1172: }