Actual source code: mg.c

petsc-master 2015-03-04
Report Typos and Errors
  2: /*
  3:     Defines the multigrid preconditioner interface.
  4: */
  5: #include <petsc-private/pcmgimpl.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(PetscOptions *PetscOptionsObject,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(PetscOptionsObject,"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",mg->default_smoothu,&m,&flg);
387:   if (flg) {
388:     PCMGSetNumberSmoothUp(pc,m);
389:   }
390:   PetscOptionsInt("-pc_mg_smoothdown","Number of pre-smoothing steps","PCMGSetNumberSmoothDown",mg->default_smoothd,&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:     /* Separately create them so we do not get DMKSP interference between levels */
598:     for (i=n-2; i>-1; i--) {DMCoarsen(dms[i+1],MPI_COMM_NULL,&dms[i]);}
599:     for (i=n-2; i>-1; i--) {
600:       DMKSP kdm;
601:       KSPSetDM(mglevels[i]->smoothd,dms[i]);
602:       if (mg->galerkin) {KSPSetDMActive(mglevels[i]->smoothd,PETSC_FALSE);}
603:       DMGetDMKSPWrite(dms[i],&kdm);
604:       /* Ugly hack so that the next KSPSetUp() will use the RHS that we set. A better fix is to change dmActive to take
605:        * a bitwise OR of computing the matrix, RHS, and initial iterate. */
606:       kdm->ops->computerhs = NULL;
607:       kdm->rhsctx          = NULL;
608:       if (!mglevels[i+1]->interpolate) {
609:         DMCreateInterpolation(dms[i],dms[i+1],&p,&rscale);
610:         PCMGSetInterpolation(pc,i+1,p);
611:         if (rscale) {PCMGSetRScale(pc,i+1,rscale);}
612:         VecDestroy(&rscale);
613:         MatDestroy(&p);
614:       }
615:     }

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

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

627:   if (mg->galerkin == 1) {
628:     Mat B;
629:     /* currently only handle case where mat and pmat are the same on coarser levels */
630:     KSPGetOperators(mglevels[n-1]->smoothd,&dA,&dB);
631:     if (!pc->setupcalled) {
632:       for (i=n-2; i>-1; i--) {
633:         if (!mglevels[i+1]->restrct && !mglevels[i+1]->interpolate) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must provide interpolation or restriction for each MG level except level 0");
634:         if (!mglevels[i+1]->interpolate) {
635:           PCMGSetInterpolation(pc,i+1,mglevels[i+1]->restrct);
636:         }
637:         if (!mglevels[i+1]->restrct) {
638:           PCMGSetRestriction(pc,i+1,mglevels[i+1]->interpolate);
639:         }
640:         if (mglevels[i+1]->interpolate == mglevels[i+1]->restrct) {
641:           MatPtAP(dB,mglevels[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);
642:         } else {
643:           MatMatMatMult(mglevels[i+1]->restrct,dB,mglevels[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);
644:         }
645:         KSPSetOperators(mglevels[i]->smoothd,B,B);
646:         if (i != n-2) {PetscObjectDereference((PetscObject)dB);}
647:         dB = B;
648:       }
649:       if (n > 1) {PetscObjectDereference((PetscObject)dB);}
650:     } else {
651:       for (i=n-2; i>-1; i--) {
652:         if (!mglevels[i+1]->restrct && !mglevels[i+1]->interpolate) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must provide interpolation or restriction for each MG level except level 0");
653:         if (!mglevels[i+1]->interpolate) {
654:           PCMGSetInterpolation(pc,i+1,mglevels[i+1]->restrct);
655:         }
656:         if (!mglevels[i+1]->restrct) {
657:           PCMGSetRestriction(pc,i+1,mglevels[i+1]->interpolate);
658:         }
659:         KSPGetOperators(mglevels[i]->smoothd,NULL,&B);
660:         if (mglevels[i+1]->interpolate == mglevels[i+1]->restrct) {
661:           MatPtAP(dB,mglevels[i+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);
662:         } else {
663:           MatMatMatMult(mglevels[i+1]->restrct,dB,mglevels[i+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);
664:         }
665:         KSPSetOperators(mglevels[i]->smoothd,B,B);
666:         dB   = B;
667:       }
668:     }
669:   } else if (!mg->galerkin && pc->dm && pc->dm->x) {
670:     /* need to restrict Jacobian location to coarser meshes for evaluation */
671:     for (i=n-2; i>-1; i--) {
672:       Mat R;
673:       Vec rscale;
674:       if (!mglevels[i]->smoothd->dm->x) {
675:         Vec *vecs;
676:         KSPCreateVecs(mglevels[i]->smoothd,1,&vecs,0,NULL);

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

680:         PetscFree(vecs);
681:       }
682:       PCMGGetRestriction(pc,i+1,&R);
683:       PCMGGetRScale(pc,i+1,&rscale);
684:       MatRestrict(R,mglevels[i+1]->smoothd->dm->x,mglevels[i]->smoothd->dm->x);
685:       VecPointwiseMult(mglevels[i]->smoothd->dm->x,mglevels[i]->smoothd->dm->x,rscale);
686:     }
687:   }
688:   if (!mg->galerkin && pc->dm) {
689:     for (i=n-2; i>=0; i--) {
690:       DM  dmfine,dmcoarse;
691:       Mat Restrict,Inject;
692:       Vec rscale;
693:       KSPGetDM(mglevels[i+1]->smoothd,&dmfine);
694:       KSPGetDM(mglevels[i]->smoothd,&dmcoarse);
695:       PCMGGetRestriction(pc,i+1,&Restrict);
696:       PCMGGetRScale(pc,i+1,&rscale);
697:       Inject = NULL;      /* Callback should create it if it needs Injection */
698:       DMRestrict(dmfine,Restrict,rscale,Inject,dmcoarse);
699:     }
700:   }

702:   if (!pc->setupcalled) {
703:     for (i=0; i<n; i++) {
704:       KSPSetFromOptions(mglevels[i]->smoothd);
705:     }
706:     for (i=1; i<n; i++) {
707:       if (mglevels[i]->smoothu && (mglevels[i]->smoothu != mglevels[i]->smoothd)) {
708:         KSPSetFromOptions(mglevels[i]->smoothu);
709:       }
710:     }
711:     for (i=1; i<n; i++) {
712:       PCMGGetInterpolation(pc,i,&mglevels[i]->interpolate);
713:       PCMGGetRestriction(pc,i,&mglevels[i]->restrct);
714:     }
715:     for (i=0; i<n-1; i++) {
716:       if (!mglevels[i]->b) {
717:         Vec *vec;
718:         KSPCreateVecs(mglevels[i]->smoothd,1,&vec,0,NULL);
719:         PCMGSetRhs(pc,i,*vec);
720:         VecDestroy(vec);
721:         PetscFree(vec);
722:       }
723:       if (!mglevels[i]->r && i) {
724:         VecDuplicate(mglevels[i]->b,&tvec);
725:         PCMGSetR(pc,i,tvec);
726:         VecDestroy(&tvec);
727:       }
728:       if (!mglevels[i]->x) {
729:         VecDuplicate(mglevels[i]->b,&tvec);
730:         PCMGSetX(pc,i,tvec);
731:         VecDestroy(&tvec);
732:       }
733:     }
734:     if (n != 1 && !mglevels[n-1]->r) {
735:       /* PCMGSetR() on the finest level if user did not supply it */
736:       Vec *vec;
737:       KSPCreateVecs(mglevels[n-1]->smoothd,1,&vec,0,NULL);
738:       PCMGSetR(pc,n-1,*vec);
739:       VecDestroy(vec);
740:       PetscFree(vec);
741:     }
742:   }

744:   if (pc->dm) {
745:     /* need to tell all the coarser levels to rebuild the matrix using the DM for that level */
746:     for (i=0; i<n-1; i++) {
747:       if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
748:     }
749:   }

751:   for (i=1; i<n; i++) {
752:     if (mglevels[i]->smoothu == mglevels[i]->smoothd || mg->am == PC_MG_FULL || mg->am == PC_MG_KASKADE || mg->cyclesperpcapply > 1){
753:       /* if doing only down then initial guess is zero */
754:       KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);
755:     }
756:     if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
757:     KSPSetUp(mglevels[i]->smoothd);
758:     if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
759:     if (!mglevels[i]->residual) {
760:       Mat mat;
761:       KSPGetOperators(mglevels[i]->smoothd,NULL,&mat);
762:       PCMGSetResidual(pc,i,PCMGResidualDefault,mat);
763:     }
764:   }
765:   for (i=1; i<n; i++) {
766:     if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) {
767:       Mat          downmat,downpmat;

769:       /* check if operators have been set for up, if not use down operators to set them */
770:       KSPGetOperatorsSet(mglevels[i]->smoothu,&opsset,NULL);
771:       if (!opsset) {
772:         KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat);
773:         KSPSetOperators(mglevels[i]->smoothu,downmat,downpmat);
774:       }

776:       KSPSetInitialGuessNonzero(mglevels[i]->smoothu,PETSC_TRUE);
777:       if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
778:       KSPSetUp(mglevels[i]->smoothu);
779:       if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
780:     }
781:   }

783:   /*
784:       If coarse solver is not direct method then DO NOT USE preonly
785:   */
786:   PetscObjectTypeCompare((PetscObject)mglevels[0]->smoothd,KSPPREONLY,&preonly);
787:   if (preonly) {
788:     PetscObjectTypeCompare((PetscObject)cpc,PCLU,&lu);
789:     PetscObjectTypeCompare((PetscObject)cpc,PCREDUNDANT,&redundant);
790:     PetscObjectTypeCompare((PetscObject)cpc,PCCHOLESKY,&cholesky);
791:     PetscObjectTypeCompare((PetscObject)cpc,PCSVD,&svd);
792:     if (!lu && !redundant && !cholesky && !svd) {
793:       KSPSetType(mglevels[0]->smoothd,KSPGMRES);
794:     }
795:   }

797:   if (mglevels[0]->eventsmoothsetup) {PetscLogEventBegin(mglevels[0]->eventsmoothsetup,0,0,0,0);}
798:   KSPSetUp(mglevels[0]->smoothd);
799:   if (mglevels[0]->eventsmoothsetup) {PetscLogEventEnd(mglevels[0]->eventsmoothsetup,0,0,0,0);}

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

806:    Only support one or the other at the same time.
807:   */
808: #if defined(PETSC_USE_SOCKET_VIEWER)
809:   PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_mg_dump_matlab",&dump,NULL);
810:   if (dump) viewer = PETSC_VIEWER_SOCKET_(PetscObjectComm((PetscObject)pc));
811:   dump = PETSC_FALSE;
812: #endif
813:   PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_mg_dump_binary",&dump,NULL);
814:   if (dump) viewer = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)pc));

816:   if (viewer) {
817:     for (i=1; i<n; i++) {
818:       MatView(mglevels[i]->restrct,viewer);
819:     }
820:     for (i=0; i<n; i++) {
821:       KSPGetPC(mglevels[i]->smoothd,&pc);
822:       MatView(pc->mat,viewer);
823:     }
824:   }
825:   return(0);
826: }

828: /* -------------------------------------------------------------------------------------*/

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

835:    Not Collective

837:    Input Parameter:
838: .  pc - the preconditioner context

840:    Output parameter:
841: .  levels - the number of levels

843:    Level: advanced

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

847: .seealso: PCMGSetLevels()
848: @*/
849: PetscErrorCode  PCMGGetLevels(PC pc,PetscInt *levels)
850: {
851:   PC_MG *mg = (PC_MG*)pc->data;

856:   *levels = mg->nlevels;
857:   return(0);
858: }

862: /*@
863:    PCMGSetType - Determines the form of multigrid to use:
864:    multiplicative, additive, full, or the Kaskade algorithm.

866:    Logically Collective on PC

868:    Input Parameters:
869: +  pc - the preconditioner context
870: -  form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,
871:    PC_MG_FULL, PC_MG_KASKADE

873:    Options Database Key:
874: .  -pc_mg_type <form> - Sets <form>, one of multiplicative,
875:    additive, full, kaskade

877:    Level: advanced

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

881: .seealso: PCMGSetLevels()
882: @*/
883: PetscErrorCode  PCMGSetType(PC pc,PCMGType form)
884: {
885:   PC_MG *mg = (PC_MG*)pc->data;

890:   mg->am = form;
891:   if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG;
892:   else pc->ops->applyrichardson = NULL;
893:   return(0);
894: }

896: /*@
897:    PCMGGetType - Determines the form of multigrid to use:
898:    multiplicative, additive, full, or the Kaskade algorithm.

900:    Logically Collective on PC

902:    Input Parameter:
903: .  pc - the preconditioner context

905:    Output Parameter:
906: .  type - one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,PC_MG_FULL, PC_MG_KASKADE


909:    Level: advanced

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

913: .seealso: PCMGSetLevels()
914: @*/
915: PetscErrorCode  PCMGGetType(PC pc,PCMGType *type)
916: {
917:   PC_MG *mg = (PC_MG*)pc->data;

921:   *type = mg->am;
922:   return(0);
923: }

927: /*@
928:    PCMGSetCycleType - Sets the type cycles to use.  Use PCMGSetCycleTypeOnLevel() for more
929:    complicated cycling.

931:    Logically Collective on PC

933:    Input Parameters:
934: +  pc - the multigrid context
935: -  PC_MG_CYCLE_V or PC_MG_CYCLE_W

937:    Options Database Key:
938: .  -pc_mg_cycle_type <v,w>

940:    Level: advanced

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

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

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

958:   for (i=0; i<levels; i++) mglevels[i]->cycles = n;
959:   return(0);
960: }

964: /*@
965:    PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step
966:          of multigrid when PCMGType of PC_MG_MULTIPLICATIVE is used

968:    Logically Collective on PC

970:    Input Parameters:
971: +  pc - the multigrid context
972: -  n - number of cycles (default is 1)

974:    Options Database Key:
975: .  -pc_mg_multiplicative_cycles n

977:    Level: advanced

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

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

983: .seealso: PCMGSetCycleTypeOnLevel(), PCMGSetCycleType()
984: @*/
985: PetscErrorCode  PCMGMultiplicativeSetCycles(PC pc,PetscInt n)
986: {
987:   PC_MG        *mg        = (PC_MG*)pc->data;

992:   mg->cyclesperpcapply = n;
993:   return(0);
994: }

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

1002:    Logically Collective on PC

1004:    Input Parameters:
1005: +  pc - the multigrid context
1006: -  use - PETSC_TRUE to use the Galerkin process to compute coarse-level operators

1008:    Options Database Key:
1009: .  -pc_mg_galerkin

1011:    Level: intermediate

1013: .keywords: MG, set, Galerkin

1015: .seealso: PCMGGetGalerkin()

1017: @*/
1018: PetscErrorCode PCMGSetGalerkin(PC pc,PetscBool use)
1019: {
1020:   PC_MG *mg = (PC_MG*)pc->data;

1024:   mg->galerkin = use ? 1 : 0;
1025:   return(0);
1026: }

1030: /*@
1031:    PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e.
1032:       A_i-1 = r_i * A_i * p_i

1034:    Not Collective

1036:    Input Parameter:
1037: .  pc - the multigrid context

1039:    Output Parameter:
1040: .  galerkin - PETSC_TRUE or PETSC_FALSE

1042:    Options Database Key:
1043: .  -pc_mg_galerkin

1045:    Level: intermediate

1047: .keywords: MG, set, Galerkin

1049: .seealso: PCMGSetGalerkin()

1051: @*/
1052: PetscErrorCode  PCMGGetGalerkin(PC pc,PetscBool  *galerkin)
1053: {
1054:   PC_MG *mg = (PC_MG*)pc->data;

1058:   *galerkin = (PetscBool)mg->galerkin;
1059:   return(0);
1060: }

1064: /*@
1065:    PCMGSetNumberSmoothDown - Sets the number of pre-smoothing steps to
1066:    use on all levels. Use PCMGGetSmootherDown() to set different
1067:    pre-smoothing steps on different levels.

1069:    Logically Collective on PC

1071:    Input Parameters:
1072: +  mg - the multigrid context
1073: -  n - the number of smoothing steps

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

1078:    Level: advanced

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

1082: .seealso: PCMGSetNumberSmoothUp()
1083: @*/
1084: PetscErrorCode  PCMGSetNumberSmoothDown(PC pc,PetscInt n)
1085: {
1086:   PC_MG          *mg        = (PC_MG*)pc->data;
1087:   PC_MG_Levels   **mglevels = mg->levels;
1089:   PetscInt       i,levels;

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

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

1102:     mg->default_smoothd = n;
1103:   }
1104:   return(0);
1105: }

1109: /*@
1110:    PCMGSetNumberSmoothUp - Sets the number of post-smoothing steps to use
1111:    on all levels. Use PCMGGetSmootherUp() to set different numbers of
1112:    post-smoothing steps on different levels.

1114:    Logically Collective on PC

1116:    Input Parameters:
1117: +  mg - the multigrid context
1118: -  n - the number of smoothing steps

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

1123:    Level: advanced

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

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

1130: .seealso: PCMGSetNumberSmoothDown()
1131: @*/
1132: PetscErrorCode  PCMGSetNumberSmoothUp(PC pc,PetscInt n)
1133: {
1134:   PC_MG          *mg        = (PC_MG*)pc->data;
1135:   PC_MG_Levels   **mglevels = mg->levels;
1137:   PetscInt       i,levels;

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

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

1150:     mg->default_smoothu = n;
1151:   }
1152:   return(0);
1153: }

1155: /* ----------------------------------------------------------------------------------------*/

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

1161:    Options Database Keys:
1162: +  -pc_mg_levels <nlevels> - number of levels including finest
1163: .  -pc_mg_cycles <v,w> -
1164: .  -pc_mg_smoothup <n> - number of smoothing steps after interpolation
1165: .  -pc_mg_smoothdown <n> - number of smoothing steps before applying restriction operator
1166: .  -pc_mg_type <additive,multiplicative,full,kaskade> - multiplicative is the default
1167: .  -pc_mg_log - log information about time spent on each level of the solver
1168: .  -pc_mg_monitor - print information on the multigrid convergence
1169: .  -pc_mg_galerkin - use Galerkin process to compute coarser operators, i.e. Acoarse = R A R'
1170: .  -pc_mg_multiplicative_cycles - number of cycles to use as the preconditioner (defaults to 1)
1171: .  -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices
1172:                         to the Socket viewer for reading from MATLAB.
1173: -  -pc_mg_dump_binary - dumps the matrices for each level and the restriction/interpolation matrices
1174:                         to the binary output file called binaryoutput

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

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

1180:    Level: intermediate

1182:    Concepts: multigrid/multilevel

1184: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, PCEXOTIC, PCGAMG, PCML, PCHYPRE
1185:            PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), PCMGSetNumberSmoothDown(),
1186:            PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(),
1187:            PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(),
1188:            PCMGSetCycleTypeOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR()
1189: M*/

1193: PETSC_EXTERN PetscErrorCode PCCreate_MG(PC pc)
1194: {
1195:   PC_MG          *mg;

1199:   PetscNewLog(pc,&mg);
1200:   pc->data    = (void*)mg;
1201:   mg->nlevels = -1;

1203:   pc->useAmat = PETSC_TRUE;

1205:   pc->ops->apply          = PCApply_MG;
1206:   pc->ops->setup          = PCSetUp_MG;
1207:   pc->ops->reset          = PCReset_MG;
1208:   pc->ops->destroy        = PCDestroy_MG;
1209:   pc->ops->setfromoptions = PCSetFromOptions_MG;
1210:   pc->ops->view           = PCView_MG;
1211:   return(0);
1212: }