Actual source code: mg.c

petsc-3.6.1 2015-07-22
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:     KSPSetErrorIfNotConverged(mglevels[i]->smoothd,pc->erroriffailure);
226:     KSPSetType(mglevels[i]->smoothd,KSPCHEBYSHEV);
227:     KSPSetConvergenceTest(mglevels[i]->smoothd,KSPConvergedSkip,NULL,NULL);
228:     KSPSetNormType(mglevels[i]->smoothd,KSP_NORM_NONE);
229:     KSPGetPC(mglevels[i]->smoothd,&ipc);
230:     PCSetType(ipc,PCSOR);
231:     PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd,(PetscObject)pc,levels-i);
232:     KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, i ? mg->default_smoothd : 1);
233:     KSPSetOptionsPrefix(mglevels[i]->smoothd,prefix);

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

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

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


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

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



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

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

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

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

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


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

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

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

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

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

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

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

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

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

535: #include <petsc/private/dmimpl.h>
536: #include <petsc/private/kspimpl.h>

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

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


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

582:   if (!opsset) {
583:     PCGetUseAmat(pc,&use_amat);
584:     if(use_amat){
585:       PetscInfo(pc,"Using outer 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->mat,pc->pmat);
587:     }
588:     else {
589:       PetscInfo(pc,"Using matrix (pmat) operators to define finest grid operator \n  because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
590:       KSPSetOperators(mglevels[n-1]->smoothd,pc->pmat,pc->pmat);
591:     }
592:   }

594:   for (i=n-1; i>0; i--) {
595:     if (!(mglevels[i]->interpolate || mglevels[i]->restrct)) {
596:       missinginterpolate = PETSC_TRUE;
597:       continue;
598:     }
599:   }
600:   /*
601:    Skipping if user has provided all interpolation/restriction needed (since DM might not be able to produce them (when coming from SNES/TS)
602:    Skipping for galerkin==2 (externally managed hierarchy such as ML and GAMG). Cleaner logic here would be great. Wrap ML/GAMG as DMs?
603:   */
604:   if (missinginterpolate && pc->dm && mg->galerkin != 2 && !pc->setupcalled) {
605:     /* construct the interpolation from the DMs */
606:     Mat p;
607:     Vec rscale;
608:     PetscMalloc1(n,&dms);
609:     dms[n-1] = pc->dm;
610:     /* Separately create them so we do not get DMKSP interference between levels */
611:     for (i=n-2; i>-1; i--) {DMCoarsen(dms[i+1],MPI_COMM_NULL,&dms[i]);}
612:     for (i=n-2; i>-1; i--) {
613:       DMKSP kdm;
614:       KSPSetDM(mglevels[i]->smoothd,dms[i]);
615:       if (mg->galerkin) {KSPSetDMActive(mglevels[i]->smoothd,PETSC_FALSE);}
616:       DMGetDMKSPWrite(dms[i],&kdm);
617:       /* Ugly hack so that the next KSPSetUp() will use the RHS that we set. A better fix is to change dmActive to take
618:        * a bitwise OR of computing the matrix, RHS, and initial iterate. */
619:       kdm->ops->computerhs = NULL;
620:       kdm->rhsctx          = NULL;
621:       if (!mglevels[i+1]->interpolate) {
622:         DMCreateInterpolation(dms[i],dms[i+1],&p,&rscale);
623:         PCMGSetInterpolation(pc,i+1,p);
624:         if (rscale) {PCMGSetRScale(pc,i+1,rscale);}
625:         VecDestroy(&rscale);
626:         MatDestroy(&p);
627:       }
628:     }

630:     for (i=n-2; i>-1; i--) {DMDestroy(&dms[i]);}
631:     PetscFree(dms);
632:   }

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

640:   if (mg->galerkin == 1) {
641:     Mat B;
642:     /* currently only handle case where mat and pmat are the same on coarser levels */
643:     KSPGetOperators(mglevels[n-1]->smoothd,&dA,&dB);
644:     if (!pc->setupcalled) {
645:       for (i=n-2; i>-1; i--) {
646:         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");
647:         if (!mglevels[i+1]->interpolate) {
648:           PCMGSetInterpolation(pc,i+1,mglevels[i+1]->restrct);
649:         }
650:         if (!mglevels[i+1]->restrct) {
651:           PCMGSetRestriction(pc,i+1,mglevels[i+1]->interpolate);
652:         }
653:         if (mglevels[i+1]->interpolate == mglevels[i+1]->restrct) {
654:           MatPtAP(dB,mglevels[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);
655:         } else {
656:           MatMatMatMult(mglevels[i+1]->restrct,dB,mglevels[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);
657:         }
658:         KSPSetOperators(mglevels[i]->smoothd,B,B);
659:         if (i != n-2) {PetscObjectDereference((PetscObject)dB);}
660:         dB = B;
661:       }
662:       if (n > 1) {PetscObjectDereference((PetscObject)dB);}
663:     } else {
664:       for (i=n-2; i>-1; i--) {
665:         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");
666:         if (!mglevels[i+1]->interpolate) {
667:           PCMGSetInterpolation(pc,i+1,mglevels[i+1]->restrct);
668:         }
669:         if (!mglevels[i+1]->restrct) {
670:           PCMGSetRestriction(pc,i+1,mglevels[i+1]->interpolate);
671:         }
672:         KSPGetOperators(mglevels[i]->smoothd,NULL,&B);
673:         if (mglevels[i+1]->interpolate == mglevels[i+1]->restrct) {
674:           MatPtAP(dB,mglevels[i+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);
675:         } else {
676:           MatMatMatMult(mglevels[i+1]->restrct,dB,mglevels[i+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);
677:         }
678:         KSPSetOperators(mglevels[i]->smoothd,B,B);
679:         dB   = B;
680:       }
681:     }
682:   } else if (!mg->galerkin && pc->dm && pc->dm->x) {
683:     /* need to restrict Jacobian location to coarser meshes for evaluation */
684:     for (i=n-2; i>-1; i--) {
685:       Mat R;
686:       Vec rscale;
687:       if (!mglevels[i]->smoothd->dm->x) {
688:         Vec *vecs;
689:         KSPCreateVecs(mglevels[i]->smoothd,1,&vecs,0,NULL);

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

693:         PetscFree(vecs);
694:       }
695:       PCMGGetRestriction(pc,i+1,&R);
696:       PCMGGetRScale(pc,i+1,&rscale);
697:       MatRestrict(R,mglevels[i+1]->smoothd->dm->x,mglevels[i]->smoothd->dm->x);
698:       VecPointwiseMult(mglevels[i]->smoothd->dm->x,mglevels[i]->smoothd->dm->x,rscale);
699:     }
700:   }
701:   if (!mg->galerkin && pc->dm) {
702:     for (i=n-2; i>=0; i--) {
703:       DM  dmfine,dmcoarse;
704:       Mat Restrict,Inject;
705:       Vec rscale;
706:       KSPGetDM(mglevels[i+1]->smoothd,&dmfine);
707:       KSPGetDM(mglevels[i]->smoothd,&dmcoarse);
708:       PCMGGetRestriction(pc,i+1,&Restrict);
709:       PCMGGetRScale(pc,i+1,&rscale);
710:       Inject = NULL;      /* Callback should create it if it needs Injection */
711:       DMRestrict(dmfine,Restrict,rscale,Inject,dmcoarse);
712:     }
713:   }

715:   if (!pc->setupcalled) {
716:     for (i=0; i<n; i++) {
717:       KSPSetFromOptions(mglevels[i]->smoothd);
718:     }
719:     for (i=1; i<n; i++) {
720:       if (mglevels[i]->smoothu && (mglevels[i]->smoothu != mglevels[i]->smoothd)) {
721:         KSPSetFromOptions(mglevels[i]->smoothu);
722:       }
723:     }
724:     for (i=1; i<n; i++) {
725:       PCMGGetInterpolation(pc,i,&mglevels[i]->interpolate);
726:       PCMGGetRestriction(pc,i,&mglevels[i]->restrct);
727:     }
728:     for (i=0; i<n-1; i++) {
729:       if (!mglevels[i]->b) {
730:         Vec *vec;
731:         KSPCreateVecs(mglevels[i]->smoothd,1,&vec,0,NULL);
732:         PCMGSetRhs(pc,i,*vec);
733:         VecDestroy(vec);
734:         PetscFree(vec);
735:       }
736:       if (!mglevels[i]->r && i) {
737:         VecDuplicate(mglevels[i]->b,&tvec);
738:         PCMGSetR(pc,i,tvec);
739:         VecDestroy(&tvec);
740:       }
741:       if (!mglevels[i]->x) {
742:         VecDuplicate(mglevels[i]->b,&tvec);
743:         PCMGSetX(pc,i,tvec);
744:         VecDestroy(&tvec);
745:       }
746:     }
747:     if (n != 1 && !mglevels[n-1]->r) {
748:       /* PCMGSetR() on the finest level if user did not supply it */
749:       Vec *vec;
750:       KSPCreateVecs(mglevels[n-1]->smoothd,1,&vec,0,NULL);
751:       PCMGSetR(pc,n-1,*vec);
752:       VecDestroy(vec);
753:       PetscFree(vec);
754:     }
755:   }

757:   if (pc->dm) {
758:     /* need to tell all the coarser levels to rebuild the matrix using the DM for that level */
759:     for (i=0; i<n-1; i++) {
760:       if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
761:     }
762:   }

764:   for (i=1; i<n; i++) {
765:     if (mglevels[i]->smoothu == mglevels[i]->smoothd || mg->am == PC_MG_FULL || mg->am == PC_MG_KASKADE || mg->cyclesperpcapply > 1){
766:       /* if doing only down then initial guess is zero */
767:       KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);
768:     }
769:     if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
770:     KSPSetUp(mglevels[i]->smoothd);
771:     if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
772:     if (!mglevels[i]->residual) {
773:       Mat mat;
774:       KSPGetOperators(mglevels[i]->smoothd,NULL,&mat);
775:       PCMGSetResidual(pc,i,PCMGResidualDefault,mat);
776:     }
777:   }
778:   for (i=1; i<n; i++) {
779:     if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) {
780:       Mat          downmat,downpmat;

782:       /* check if operators have been set for up, if not use down operators to set them */
783:       KSPGetOperatorsSet(mglevels[i]->smoothu,&opsset,NULL);
784:       if (!opsset) {
785:         KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat);
786:         KSPSetOperators(mglevels[i]->smoothu,downmat,downpmat);
787:       }

789:       KSPSetInitialGuessNonzero(mglevels[i]->smoothu,PETSC_TRUE);
790:       if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
791:       KSPSetUp(mglevels[i]->smoothu);
792:       if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
793:     }
794:   }

796:   /*
797:       If coarse solver is not direct method then DO NOT USE preonly
798:   */
799:   PetscObjectTypeCompare((PetscObject)mglevels[0]->smoothd,KSPPREONLY,&preonly);
800:   if (preonly) {
801:     PetscObjectTypeCompare((PetscObject)cpc,PCLU,&lu);
802:     PetscObjectTypeCompare((PetscObject)cpc,PCREDUNDANT,&redundant);
803:     PetscObjectTypeCompare((PetscObject)cpc,PCCHOLESKY,&cholesky);
804:     PetscObjectTypeCompare((PetscObject)cpc,PCSVD,&svd);
805:     if (!lu && !redundant && !cholesky && !svd) {
806:       KSPSetType(mglevels[0]->smoothd,KSPGMRES);
807:     }
808:   }

810:   if (mglevels[0]->eventsmoothsetup) {PetscLogEventBegin(mglevels[0]->eventsmoothsetup,0,0,0,0);}
811:   KSPSetUp(mglevels[0]->smoothd);
812:   if (mglevels[0]->eventsmoothsetup) {PetscLogEventEnd(mglevels[0]->eventsmoothsetup,0,0,0,0);}

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

819:    Only support one or the other at the same time.
820:   */
821: #if defined(PETSC_USE_SOCKET_VIEWER)
822:   PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_mg_dump_matlab",&dump,NULL);
823:   if (dump) viewer = PETSC_VIEWER_SOCKET_(PetscObjectComm((PetscObject)pc));
824:   dump = PETSC_FALSE;
825: #endif
826:   PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_mg_dump_binary",&dump,NULL);
827:   if (dump) viewer = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)pc));

829:   if (viewer) {
830:     for (i=1; i<n; i++) {
831:       MatView(mglevels[i]->restrct,viewer);
832:     }
833:     for (i=0; i<n; i++) {
834:       KSPGetPC(mglevels[i]->smoothd,&pc);
835:       MatView(pc->mat,viewer);
836:     }
837:   }
838:   return(0);
839: }

841: /* -------------------------------------------------------------------------------------*/

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

848:    Not Collective

850:    Input Parameter:
851: .  pc - the preconditioner context

853:    Output parameter:
854: .  levels - the number of levels

856:    Level: advanced

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

860: .seealso: PCMGSetLevels()
861: @*/
862: PetscErrorCode  PCMGGetLevels(PC pc,PetscInt *levels)
863: {
864:   PC_MG *mg = (PC_MG*)pc->data;

869:   *levels = mg->nlevels;
870:   return(0);
871: }

875: /*@
876:    PCMGSetType - Determines the form of multigrid to use:
877:    multiplicative, additive, full, or the Kaskade algorithm.

879:    Logically Collective on PC

881:    Input Parameters:
882: +  pc - the preconditioner context
883: -  form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,
884:    PC_MG_FULL, PC_MG_KASKADE

886:    Options Database Key:
887: .  -pc_mg_type <form> - Sets <form>, one of multiplicative,
888:    additive, full, kaskade

890:    Level: advanced

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

894: .seealso: PCMGSetLevels()
895: @*/
896: PetscErrorCode  PCMGSetType(PC pc,PCMGType form)
897: {
898:   PC_MG *mg = (PC_MG*)pc->data;

903:   mg->am = form;
904:   if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG;
905:   else pc->ops->applyrichardson = NULL;
906:   return(0);
907: }

909: /*@
910:    PCMGGetType - Determines the form of multigrid to use:
911:    multiplicative, additive, full, or the Kaskade algorithm.

913:    Logically Collective on PC

915:    Input Parameter:
916: .  pc - the preconditioner context

918:    Output Parameter:
919: .  type - one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,PC_MG_FULL, PC_MG_KASKADE


922:    Level: advanced

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

926: .seealso: PCMGSetLevels()
927: @*/
928: PetscErrorCode  PCMGGetType(PC pc,PCMGType *type)
929: {
930:   PC_MG *mg = (PC_MG*)pc->data;

934:   *type = mg->am;
935:   return(0);
936: }

940: /*@
941:    PCMGSetCycleType - Sets the type cycles to use.  Use PCMGSetCycleTypeOnLevel() for more
942:    complicated cycling.

944:    Logically Collective on PC

946:    Input Parameters:
947: +  pc - the multigrid context
948: -  PC_MG_CYCLE_V or PC_MG_CYCLE_W

950:    Options Database Key:
951: .  -pc_mg_cycle_type <v,w>

953:    Level: advanced

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

957: .seealso: PCMGSetCycleTypeOnLevel()
958: @*/
959: PetscErrorCode  PCMGSetCycleType(PC pc,PCMGCycleType n)
960: {
961:   PC_MG        *mg        = (PC_MG*)pc->data;
962:   PC_MG_Levels **mglevels = mg->levels;
963:   PetscInt     i,levels;

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

971:   for (i=0; i<levels; i++) mglevels[i]->cycles = n;
972:   return(0);
973: }

977: /*@
978:    PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step
979:          of multigrid when PCMGType of PC_MG_MULTIPLICATIVE is used

981:    Logically Collective on PC

983:    Input Parameters:
984: +  pc - the multigrid context
985: -  n - number of cycles (default is 1)

987:    Options Database Key:
988: .  -pc_mg_multiplicative_cycles n

990:    Level: advanced

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

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

996: .seealso: PCMGSetCycleTypeOnLevel(), PCMGSetCycleType()
997: @*/
998: PetscErrorCode  PCMGMultiplicativeSetCycles(PC pc,PetscInt n)
999: {
1000:   PC_MG        *mg        = (PC_MG*)pc->data;

1005:   mg->cyclesperpcapply = n;
1006:   return(0);
1007: }

1011: PetscErrorCode PCMGSetGalerkin_MG(PC pc,PetscBool use)
1012: {
1013:   PC_MG *mg = (PC_MG*)pc->data;

1016:   mg->galerkin = use ? 1 : 0;
1017:   return(0);
1018: }

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

1026:    Logically Collective on PC

1028:    Input Parameters:
1029: +  pc - the multigrid context
1030: -  use - PETSC_TRUE to use the Galerkin process to compute coarse-level operators

1032:    Options Database Key:
1033: .  -pc_mg_galerkin <true,false>

1035:    Level: intermediate

1037:    Notes: Some codes that use PCMG such as PCGAMG use Galerkin internally while constructing the hierarchy and thus do not
1038:      use the PCMG construction of the coarser grids.

1040: .keywords: MG, set, Galerkin

1042: .seealso: PCMGGetGalerkin()

1044: @*/
1045: PetscErrorCode PCMGSetGalerkin(PC pc,PetscBool use)
1046: {

1051:   PetscTryMethod(pc,"PCMGSetGalerkin_C",(PC,PetscBool),(pc,use));
1052:   return(0);
1053: }

1057: /*@
1058:    PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e.
1059:       A_i-1 = r_i * A_i * p_i

1061:    Not Collective

1063:    Input Parameter:
1064: .  pc - the multigrid context

1066:    Output Parameter:
1067: .  galerkin - PETSC_TRUE or PETSC_FALSE

1069:    Options Database Key:
1070: .  -pc_mg_galerkin

1072:    Level: intermediate

1074: .keywords: MG, set, Galerkin

1076: .seealso: PCMGSetGalerkin()

1078: @*/
1079: PetscErrorCode  PCMGGetGalerkin(PC pc,PetscBool  *galerkin)
1080: {
1081:   PC_MG *mg = (PC_MG*)pc->data;

1085:   *galerkin = (PetscBool)mg->galerkin;
1086:   return(0);
1087: }

1091: /*@
1092:    PCMGSetNumberSmoothDown - Sets the number of pre-smoothing steps to
1093:    use on all levels. Use PCMGGetSmootherDown() to set different
1094:    pre-smoothing steps on different levels.

1096:    Logically Collective on PC

1098:    Input Parameters:
1099: +  mg - the multigrid context
1100: -  n - the number of smoothing steps

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

1105:    Level: advanced

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

1109: .seealso: PCMGSetNumberSmoothUp()
1110: @*/
1111: PetscErrorCode  PCMGSetNumberSmoothDown(PC pc,PetscInt n)
1112: {
1113:   PC_MG          *mg        = (PC_MG*)pc->data;
1114:   PC_MG_Levels   **mglevels = mg->levels;
1116:   PetscInt       i,levels;

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

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

1129:     mg->default_smoothd = n;
1130:   }
1131:   return(0);
1132: }

1136: /*@
1137:    PCMGSetNumberSmoothUp - Sets the number of post-smoothing steps to use
1138:    on all levels. Use PCMGGetSmootherUp() to set different numbers of
1139:    post-smoothing steps on different levels.

1141:    Logically Collective on PC

1143:    Input Parameters:
1144: +  mg - the multigrid context
1145: -  n - the number of smoothing steps

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

1150:    Level: advanced

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

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

1157: .seealso: PCMGSetNumberSmoothDown()
1158: @*/
1159: PetscErrorCode  PCMGSetNumberSmoothUp(PC pc,PetscInt n)
1160: {
1161:   PC_MG          *mg        = (PC_MG*)pc->data;
1162:   PC_MG_Levels   **mglevels = mg->levels;
1164:   PetscInt       i,levels;

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

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

1177:     mg->default_smoothu = n;
1178:   }
1179:   return(0);
1180: }

1182: /* ----------------------------------------------------------------------------------------*/

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

1188:    Options Database Keys:
1189: +  -pc_mg_levels <nlevels> - number of levels including finest
1190: .  -pc_mg_cycles <v,w> -
1191: .  -pc_mg_smoothup <n> - number of smoothing steps after interpolation
1192: .  -pc_mg_smoothdown <n> - number of smoothing steps before applying restriction operator
1193: .  -pc_mg_type <additive,multiplicative,full,kaskade> - multiplicative is the default
1194: .  -pc_mg_log - log information about time spent on each level of the solver
1195: .  -pc_mg_galerkin - use Galerkin process to compute coarser operators, i.e. Acoarse = R A R'
1196: .  -pc_mg_multiplicative_cycles - number of cycles to use as the preconditioner (defaults to 1)
1197: .  -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices
1198:                         to the Socket viewer for reading from MATLAB.
1199: -  -pc_mg_dump_binary - dumps the matrices for each level and the restriction/interpolation matrices
1200:                         to the binary output file called binaryoutput

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

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

1206:    Level: intermediate

1208:    Concepts: multigrid/multilevel

1210: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, PCEXOTIC, PCGAMG, PCML, PCHYPRE
1211:            PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), PCMGSetNumberSmoothDown(),
1212:            PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(),
1213:            PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(),
1214:            PCMGSetCycleTypeOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR()
1215: M*/

1219: PETSC_EXTERN PetscErrorCode PCCreate_MG(PC pc)
1220: {
1221:   PC_MG          *mg;

1225:   PetscNewLog(pc,&mg);
1226:   pc->data    = (void*)mg;
1227:   mg->nlevels = -1;

1229:   pc->useAmat = PETSC_TRUE;

1231:   pc->ops->apply          = PCApply_MG;
1232:   pc->ops->setup          = PCSetUp_MG;
1233:   pc->ops->reset          = PCReset_MG;
1234:   pc->ops->destroy        = PCDestroy_MG;
1235:   pc->ops->setfromoptions = PCSetFromOptions_MG;
1236:   pc->ops->view           = PCView_MG;

1238:   PetscObjectComposeFunction((PetscObject)pc,"PCMGSetGalerkin_C",PCMGSetGalerkin_MG);
1239:   return(0);
1240: }