Actual source code: mg.c
petsc-3.6.4 2016-04-12
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: PCMGType mgtype = mg->am;
178: PetscInt mgctype = (PetscInt) PC_MG_CYCLE_V;
179: PetscInt i;
180: PetscMPIInt size;
181: const char *prefix;
182: PC ipc;
183: PetscInt n;
188: PetscObjectGetComm((PetscObject)pc,&comm);
189: if (mg->nlevels == levels) return(0);
190: if (mglevels) {
191: mgctype = mglevels[0]->cycles;
192: /* changing the number of levels so free up the previous stuff */
193: PCReset_MG(pc);
194: n = mglevels[0]->levels;
195: for (i=0; i<n; i++) {
196: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
197: KSPDestroy(&mglevels[i]->smoothd);
198: }
199: KSPDestroy(&mglevels[i]->smoothu);
200: PetscFree(mglevels[i]);
201: }
202: PetscFree(mg->levels);
203: }
205: mg->nlevels = levels;
207: PetscMalloc1(levels,&mglevels);
208: PetscLogObjectMemory((PetscObject)pc,levels*(sizeof(PC_MG*)));
210: PCGetOptionsPrefix(pc,&prefix);
212: mg->stageApply = 0;
213: for (i=0; i<levels; i++) {
214: PetscNewLog(pc,&mglevels[i]);
216: mglevels[i]->level = i;
217: mglevels[i]->levels = levels;
218: mglevels[i]->cycles = mgctype;
219: mg->default_smoothu = 2;
220: mg->default_smoothd = 2;
221: mglevels[i]->eventsmoothsetup = 0;
222: mglevels[i]->eventsmoothsolve = 0;
223: mglevels[i]->eventresidual = 0;
224: mglevels[i]->eventinterprestrict = 0;
226: if (comms) comm = comms[i];
227: KSPCreate(comm,&mglevels[i]->smoothd);
228: KSPSetErrorIfNotConverged(mglevels[i]->smoothd,pc->erroriffailure);
229: KSPSetType(mglevels[i]->smoothd,KSPCHEBYSHEV);
230: KSPSetConvergenceTest(mglevels[i]->smoothd,KSPConvergedSkip,NULL,NULL);
231: KSPSetNormType(mglevels[i]->smoothd,KSP_NORM_NONE);
232: KSPGetPC(mglevels[i]->smoothd,&ipc);
233: PCSetType(ipc,PCSOR);
234: PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd,(PetscObject)pc,levels-i);
235: KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, i ? mg->default_smoothd : 1);
236: KSPSetOptionsPrefix(mglevels[i]->smoothd,prefix);
238: /* do special stuff for coarse grid */
239: if (!i && levels > 1) {
240: KSPAppendOptionsPrefix(mglevels[0]->smoothd,"mg_coarse_");
242: /* coarse solve is (redundant) LU by default; set shifttype NONZERO to avoid annoying zero-pivot in LU preconditioner */
243: KSPSetType(mglevels[0]->smoothd,KSPPREONLY);
244: KSPGetPC(mglevels[0]->smoothd,&ipc);
245: MPI_Comm_size(comm,&size);
246: if (size > 1) {
247: KSP innerksp;
248: PC innerpc;
249: PCSetType(ipc,PCREDUNDANT);
250: PCRedundantGetKSP(ipc,&innerksp);
251: KSPGetPC(innerksp,&innerpc);
252: PCFactorSetShiftType(innerpc,MAT_SHIFT_INBLOCKS);
253: } else {
254: PCSetType(ipc,PCLU);
255: PCFactorSetShiftType(ipc,MAT_SHIFT_INBLOCKS);
256: }
257: } else {
258: char tprefix[128];
259: sprintf(tprefix,"mg_levels_%d_",(int)i);
260: KSPAppendOptionsPrefix(mglevels[i]->smoothd,tprefix);
261: }
262: PetscLogObjectParent((PetscObject)pc,(PetscObject)mglevels[i]->smoothd);
264: mglevels[i]->smoothu = mglevels[i]->smoothd;
265: mg->rtol = 0.0;
266: mg->abstol = 0.0;
267: mg->dtol = 0.0;
268: mg->ttol = 0.0;
269: mg->cyclesperpcapply = 1;
270: }
271: mg->levels = mglevels;
272: PCMGSetType(pc,mgtype);
273: return(0);
274: }
279: PetscErrorCode PCDestroy_MG(PC pc)
280: {
282: PC_MG *mg = (PC_MG*)pc->data;
283: PC_MG_Levels **mglevels = mg->levels;
284: PetscInt i,n;
287: PCReset_MG(pc);
288: if (mglevels) {
289: n = mglevels[0]->levels;
290: for (i=0; i<n; i++) {
291: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
292: KSPDestroy(&mglevels[i]->smoothd);
293: }
294: KSPDestroy(&mglevels[i]->smoothu);
295: PetscFree(mglevels[i]);
296: }
297: PetscFree(mg->levels);
298: }
299: PetscFree(pc->data);
300: return(0);
301: }
305: extern PetscErrorCode PCMGACycle_Private(PC,PC_MG_Levels**);
306: extern PetscErrorCode PCMGFCycle_Private(PC,PC_MG_Levels**);
307: extern PetscErrorCode PCMGKCycle_Private(PC,PC_MG_Levels**);
309: /*
310: PCApply_MG - Runs either an additive, multiplicative, Kaskadic
311: or full cycle of multigrid.
313: Note:
314: A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle().
315: */
318: static PetscErrorCode PCApply_MG(PC pc,Vec b,Vec x)
319: {
320: PC_MG *mg = (PC_MG*)pc->data;
321: PC_MG_Levels **mglevels = mg->levels;
323: PetscInt levels = mglevels[0]->levels,i;
326: if (mg->stageApply) {PetscLogStagePush(mg->stageApply);}
327: /* When the DM is supplying the matrix then it will not exist until here */
328: for (i=0; i<levels; i++) {
329: if (!mglevels[i]->A) {
330: KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);
331: PetscObjectReference((PetscObject)mglevels[i]->A);
332: }
333: }
335: mglevels[levels-1]->b = b;
336: mglevels[levels-1]->x = x;
337: if (mg->am == PC_MG_MULTIPLICATIVE) {
338: VecSet(x,0.0);
339: for (i=0; i<mg->cyclesperpcapply; i++) {
340: PCMGMCycle_Private(pc,mglevels+levels-1,NULL);
341: }
342: } else if (mg->am == PC_MG_ADDITIVE) {
343: PCMGACycle_Private(pc,mglevels);
344: } else if (mg->am == PC_MG_KASKADE) {
345: PCMGKCycle_Private(pc,mglevels);
346: } else {
347: PCMGFCycle_Private(pc,mglevels);
348: }
349: if (mg->stageApply) {PetscLogStagePop();}
350: return(0);
351: }
356: PetscErrorCode PCSetFromOptions_MG(PetscOptions *PetscOptionsObject,PC pc)
357: {
359: PetscInt m,levels = 1,cycles;
360: PetscBool flg,set;
361: PC_MG *mg = (PC_MG*)pc->data;
362: PC_MG_Levels **mglevels = mg->levels;
363: PCMGType mgtype;
364: PCMGCycleType mgctype;
367: PetscOptionsHead(PetscOptionsObject,"Multigrid options");
368: if (!mg->levels) {
369: PetscOptionsInt("-pc_mg_levels","Number of Levels","PCMGSetLevels",levels,&levels,&flg);
370: if (!flg && pc->dm) {
371: DMGetRefineLevel(pc->dm,&levels);
372: levels++;
373: mg->usedmfornumberoflevels = PETSC_TRUE;
374: }
375: PCMGSetLevels(pc,levels,NULL);
376: }
377: mglevels = mg->levels;
379: mgctype = (PCMGCycleType) mglevels[0]->cycles;
380: PetscOptionsEnum("-pc_mg_cycle_type","V cycle or for W-cycle","PCMGSetCycleType",PCMGCycleTypes,(PetscEnum)mgctype,(PetscEnum*)&mgctype,&flg);
381: if (flg) {
382: PCMGSetCycleType(pc,mgctype);
383: }
384: flg = PETSC_FALSE;
385: PetscOptionsBool("-pc_mg_galerkin","Use Galerkin process to compute coarser operators","PCMGSetGalerkin",flg,&flg,&set);
386: if (set) {
387: PCMGSetGalerkin(pc,flg);
388: }
389: PetscOptionsInt("-pc_mg_smoothup","Number of post-smoothing steps","PCMGSetNumberSmoothUp",mg->default_smoothu,&m,&flg);
390: if (flg) {
391: PCMGSetNumberSmoothUp(pc,m);
392: }
393: PetscOptionsInt("-pc_mg_smoothdown","Number of pre-smoothing steps","PCMGSetNumberSmoothDown",mg->default_smoothd,&m,&flg);
394: if (flg) {
395: PCMGSetNumberSmoothDown(pc,m);
396: }
397: mgtype = mg->am;
398: PetscOptionsEnum("-pc_mg_type","Multigrid type","PCMGSetType",PCMGTypes,(PetscEnum)mgtype,(PetscEnum*)&mgtype,&flg);
399: if (flg) {
400: PCMGSetType(pc,mgtype);
401: }
402: if (mg->am == PC_MG_MULTIPLICATIVE) {
403: PetscOptionsInt("-pc_mg_multiplicative_cycles","Number of cycles for each preconditioner step","PCMGMultiplicativeSetCycles",mg->cyclesperpcapply,&cycles,&flg);
404: if (flg) {
405: PCMGMultiplicativeSetCycles(pc,cycles);
406: }
407: }
408: flg = PETSC_FALSE;
409: PetscOptionsBool("-pc_mg_log","Log times for each multigrid level","None",flg,&flg,NULL);
410: if (flg) {
411: PetscInt i;
412: char eventname[128];
413: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
414: levels = mglevels[0]->levels;
415: for (i=0; i<levels; i++) {
416: sprintf(eventname,"MGSetup Level %d",(int)i);
417: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsetup);
418: sprintf(eventname,"MGSmooth Level %d",(int)i);
419: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsolve);
420: if (i) {
421: sprintf(eventname,"MGResid Level %d",(int)i);
422: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventresidual);
423: sprintf(eventname,"MGInterp Level %d",(int)i);
424: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventinterprestrict);
425: }
426: }
428: #if defined(PETSC_USE_LOG)
429: {
430: const char *sname = "MG Apply";
431: PetscStageLog stageLog;
432: PetscInt st;
435: PetscLogGetStageLog(&stageLog);
436: for (st = 0; st < stageLog->numStages; ++st) {
437: PetscBool same;
439: PetscStrcmp(stageLog->stageInfo[st].name, sname, &same);
440: if (same) mg->stageApply = st;
441: }
442: if (!mg->stageApply) {
443: PetscLogStageRegister(sname, &mg->stageApply);
444: }
445: }
446: #endif
447: }
448: PetscOptionsTail();
449: return(0);
450: }
452: const char *const PCMGTypes[] = {"MULTIPLICATIVE","ADDITIVE","FULL","KASKADE","PCMGType","PC_MG",0};
453: const char *const PCMGCycleTypes[] = {"invalid","v","w","PCMGCycleType","PC_MG_CYCLE",0};
455: #include <petscdraw.h>
458: PetscErrorCode PCView_MG(PC pc,PetscViewer viewer)
459: {
460: PC_MG *mg = (PC_MG*)pc->data;
461: PC_MG_Levels **mglevels = mg->levels;
463: PetscInt levels = mglevels ? mglevels[0]->levels : 0,i;
464: PetscBool iascii,isbinary,isdraw;
467: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
468: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
469: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
470: if (iascii) {
471: const char *cyclename = levels ? (mglevels[0]->cycles == PC_MG_CYCLE_V ? "v" : "w") : "unknown";
472: PetscViewerASCIIPrintf(viewer," MG: type is %s, levels=%D cycles=%s\n", PCMGTypes[mg->am],levels,cyclename);
473: if (mg->am == PC_MG_MULTIPLICATIVE) {
474: PetscViewerASCIIPrintf(viewer," Cycles per PCApply=%d\n",mg->cyclesperpcapply);
475: }
476: if (mg->galerkin) {
477: PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices\n");
478: } else {
479: PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid matrices\n");
480: }
481: if (mg->view){
482: (*mg->view)(pc,viewer);
483: }
484: for (i=0; i<levels; i++) {
485: if (!i) {
486: PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level -------------------------------\n",i);
487: } else {
488: PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);
489: }
490: PetscViewerASCIIPushTab(viewer);
491: KSPView(mglevels[i]->smoothd,viewer);
492: PetscViewerASCIIPopTab(viewer);
493: if (i && mglevels[i]->smoothd == mglevels[i]->smoothu) {
494: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");
495: } else if (i) {
496: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);
497: PetscViewerASCIIPushTab(viewer);
498: KSPView(mglevels[i]->smoothu,viewer);
499: PetscViewerASCIIPopTab(viewer);
500: }
501: }
502: } else if (isbinary) {
503: for (i=levels-1; i>=0; i--) {
504: KSPView(mglevels[i]->smoothd,viewer);
505: if (i && mglevels[i]->smoothd != mglevels[i]->smoothu) {
506: KSPView(mglevels[i]->smoothu,viewer);
507: }
508: }
509: } else if (isdraw) {
510: PetscDraw draw;
511: PetscReal x,w,y,bottom,th;
512: PetscViewerDrawGetDraw(viewer,0,&draw);
513: PetscDrawGetCurrentPoint(draw,&x,&y);
514: PetscDrawStringGetSize(draw,NULL,&th);
515: bottom = y - th;
516: for (i=levels-1; i>=0; i--) {
517: if (!mglevels[i]->smoothu || (mglevels[i]->smoothu == mglevels[i]->smoothd)) {
518: PetscDrawPushCurrentPoint(draw,x,bottom);
519: KSPView(mglevels[i]->smoothd,viewer);
520: PetscDrawPopCurrentPoint(draw);
521: } else {
522: w = 0.5*PetscMin(1.0-x,x);
523: PetscDrawPushCurrentPoint(draw,x+w,bottom);
524: KSPView(mglevels[i]->smoothd,viewer);
525: PetscDrawPopCurrentPoint(draw);
526: PetscDrawPushCurrentPoint(draw,x-w,bottom);
527: KSPView(mglevels[i]->smoothu,viewer);
528: PetscDrawPopCurrentPoint(draw);
529: }
530: PetscDrawGetBoundingBox(draw,NULL,&bottom,NULL,NULL);
531: bottom -= th;
532: }
533: }
534: return(0);
535: }
537: #include <petsc/private/dmimpl.h>
538: #include <petsc/private/kspimpl.h>
540: /*
541: Calls setup for the KSP on each level
542: */
545: PetscErrorCode PCSetUp_MG(PC pc)
546: {
547: PC_MG *mg = (PC_MG*)pc->data;
548: PC_MG_Levels **mglevels = mg->levels;
550: PetscInt i,n = mglevels[0]->levels;
551: PC cpc;
552: PetscBool preonly,lu,redundant,cholesky,svd,dump = PETSC_FALSE,opsset,use_amat,missinginterpolate = PETSC_FALSE;
553: Mat dA,dB;
554: Vec tvec;
555: DM *dms;
556: PetscViewer viewer = 0;
559: /* FIX: Move this to PCSetFromOptions_MG? */
560: if (mg->usedmfornumberoflevels) {
561: PetscInt levels;
562: DMGetRefineLevel(pc->dm,&levels);
563: levels++;
564: if (levels > n) { /* the problem is now being solved on a finer grid */
565: PCMGSetLevels(pc,levels,NULL);
566: n = levels;
567: PCSetFromOptions(pc); /* it is bad to call this here, but otherwise will never be called for the new hierarchy */
568: mglevels = mg->levels;
569: }
570: }
571: KSPGetPC(mglevels[0]->smoothd,&cpc);
574: /* If user did not provide fine grid operators OR operator was not updated since last global KSPSetOperators() */
575: /* so use those from global PC */
576: /* Is this what we always want? What if user wants to keep old one? */
577: KSPGetOperatorsSet(mglevels[n-1]->smoothd,NULL,&opsset);
578: if (opsset) {
579: Mat mmat;
580: KSPGetOperators(mglevels[n-1]->smoothd,NULL,&mmat);
581: if (mmat == pc->pmat) opsset = PETSC_FALSE;
582: }
584: if (!opsset) {
585: PCGetUseAmat(pc,&use_amat);
586: if(use_amat){
587: PetscInfo(pc,"Using outer operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
588: KSPSetOperators(mglevels[n-1]->smoothd,pc->mat,pc->pmat);
589: }
590: else {
591: PetscInfo(pc,"Using matrix (pmat) operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
592: KSPSetOperators(mglevels[n-1]->smoothd,pc->pmat,pc->pmat);
593: }
594: }
596: for (i=n-1; i>0; i--) {
597: if (!(mglevels[i]->interpolate || mglevels[i]->restrct)) {
598: missinginterpolate = PETSC_TRUE;
599: continue;
600: }
601: }
602: /*
603: Skipping if user has provided all interpolation/restriction needed (since DM might not be able to produce them (when coming from SNES/TS)
604: Skipping for galerkin==2 (externally managed hierarchy such as ML and GAMG). Cleaner logic here would be great. Wrap ML/GAMG as DMs?
605: */
606: if (missinginterpolate && pc->dm && mg->galerkin != 2 && !pc->setupcalled) {
607: /* construct the interpolation from the DMs */
608: Mat p;
609: Vec rscale;
610: PetscMalloc1(n,&dms);
611: dms[n-1] = pc->dm;
612: /* Separately create them so we do not get DMKSP interference between levels */
613: for (i=n-2; i>-1; i--) {DMCoarsen(dms[i+1],MPI_COMM_NULL,&dms[i]);}
614: for (i=n-2; i>-1; i--) {
615: DMKSP kdm;
616: KSPSetDM(mglevels[i]->smoothd,dms[i]);
617: if (mg->galerkin) {KSPSetDMActive(mglevels[i]->smoothd,PETSC_FALSE);}
618: DMGetDMKSPWrite(dms[i],&kdm);
619: /* Ugly hack so that the next KSPSetUp() will use the RHS that we set. A better fix is to change dmActive to take
620: * a bitwise OR of computing the matrix, RHS, and initial iterate. */
621: kdm->ops->computerhs = NULL;
622: kdm->rhsctx = NULL;
623: if (!mglevels[i+1]->interpolate) {
624: DMCreateInterpolation(dms[i],dms[i+1],&p,&rscale);
625: PCMGSetInterpolation(pc,i+1,p);
626: if (rscale) {PCMGSetRScale(pc,i+1,rscale);}
627: VecDestroy(&rscale);
628: MatDestroy(&p);
629: }
630: }
632: for (i=n-2; i>-1; i--) {DMDestroy(&dms[i]);}
633: PetscFree(dms);
634: }
636: if (pc->dm && !pc->setupcalled) {
637: /* finest smoother also gets DM but it is not active, independent of whether galerkin==2 */
638: KSPSetDM(mglevels[n-1]->smoothd,pc->dm);
639: KSPSetDMActive(mglevels[n-1]->smoothd,PETSC_FALSE);
640: }
642: if (mg->galerkin == 1) {
643: Mat B;
644: /* currently only handle case where mat and pmat are the same on coarser levels */
645: KSPGetOperators(mglevels[n-1]->smoothd,&dA,&dB);
646: if (!pc->setupcalled) {
647: for (i=n-2; i>-1; i--) {
648: 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");
649: if (!mglevels[i+1]->interpolate) {
650: PCMGSetInterpolation(pc,i+1,mglevels[i+1]->restrct);
651: }
652: if (!mglevels[i+1]->restrct) {
653: PCMGSetRestriction(pc,i+1,mglevels[i+1]->interpolate);
654: }
655: if (mglevels[i+1]->interpolate == mglevels[i+1]->restrct) {
656: MatPtAP(dB,mglevels[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);
657: } else {
658: MatMatMatMult(mglevels[i+1]->restrct,dB,mglevels[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);
659: }
660: KSPSetOperators(mglevels[i]->smoothd,B,B);
661: if (i != n-2) {PetscObjectDereference((PetscObject)dB);}
662: dB = B;
663: }
664: if (n > 1) {PetscObjectDereference((PetscObject)dB);}
665: } else {
666: for (i=n-2; i>-1; i--) {
667: 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");
668: if (!mglevels[i+1]->interpolate) {
669: PCMGSetInterpolation(pc,i+1,mglevels[i+1]->restrct);
670: }
671: if (!mglevels[i+1]->restrct) {
672: PCMGSetRestriction(pc,i+1,mglevels[i+1]->interpolate);
673: }
674: KSPGetOperators(mglevels[i]->smoothd,NULL,&B);
675: if (mglevels[i+1]->interpolate == mglevels[i+1]->restrct) {
676: MatPtAP(dB,mglevels[i+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);
677: } else {
678: MatMatMatMult(mglevels[i+1]->restrct,dB,mglevels[i+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);
679: }
680: KSPSetOperators(mglevels[i]->smoothd,B,B);
681: dB = B;
682: }
683: }
684: } else if (!mg->galerkin && pc->dm && pc->dm->x) {
685: /* need to restrict Jacobian location to coarser meshes for evaluation */
686: for (i=n-2; i>-1; i--) {
687: Mat R;
688: Vec rscale;
689: if (!mglevels[i]->smoothd->dm->x) {
690: Vec *vecs;
691: KSPCreateVecs(mglevels[i]->smoothd,1,&vecs,0,NULL);
693: mglevels[i]->smoothd->dm->x = vecs[0];
695: PetscFree(vecs);
696: }
697: PCMGGetRestriction(pc,i+1,&R);
698: PCMGGetRScale(pc,i+1,&rscale);
699: MatRestrict(R,mglevels[i+1]->smoothd->dm->x,mglevels[i]->smoothd->dm->x);
700: VecPointwiseMult(mglevels[i]->smoothd->dm->x,mglevels[i]->smoothd->dm->x,rscale);
701: }
702: }
703: if (!mg->galerkin && pc->dm) {
704: for (i=n-2; i>=0; i--) {
705: DM dmfine,dmcoarse;
706: Mat Restrict,Inject;
707: Vec rscale;
708: KSPGetDM(mglevels[i+1]->smoothd,&dmfine);
709: KSPGetDM(mglevels[i]->smoothd,&dmcoarse);
710: PCMGGetRestriction(pc,i+1,&Restrict);
711: PCMGGetRScale(pc,i+1,&rscale);
712: Inject = NULL; /* Callback should create it if it needs Injection */
713: DMRestrict(dmfine,Restrict,rscale,Inject,dmcoarse);
714: }
715: }
717: if (!pc->setupcalled) {
718: for (i=0; i<n; i++) {
719: KSPSetFromOptions(mglevels[i]->smoothd);
720: }
721: for (i=1; i<n; i++) {
722: if (mglevels[i]->smoothu && (mglevels[i]->smoothu != mglevels[i]->smoothd)) {
723: KSPSetFromOptions(mglevels[i]->smoothu);
724: }
725: }
726: for (i=1; i<n; i++) {
727: PCMGGetInterpolation(pc,i,&mglevels[i]->interpolate);
728: PCMGGetRestriction(pc,i,&mglevels[i]->restrct);
729: }
730: for (i=0; i<n-1; i++) {
731: if (!mglevels[i]->b) {
732: Vec *vec;
733: KSPCreateVecs(mglevels[i]->smoothd,1,&vec,0,NULL);
734: PCMGSetRhs(pc,i,*vec);
735: VecDestroy(vec);
736: PetscFree(vec);
737: }
738: if (!mglevels[i]->r && i) {
739: VecDuplicate(mglevels[i]->b,&tvec);
740: PCMGSetR(pc,i,tvec);
741: VecDestroy(&tvec);
742: }
743: if (!mglevels[i]->x) {
744: VecDuplicate(mglevels[i]->b,&tvec);
745: PCMGSetX(pc,i,tvec);
746: VecDestroy(&tvec);
747: }
748: }
749: if (n != 1 && !mglevels[n-1]->r) {
750: /* PCMGSetR() on the finest level if user did not supply it */
751: Vec *vec;
752: KSPCreateVecs(mglevels[n-1]->smoothd,1,&vec,0,NULL);
753: PCMGSetR(pc,n-1,*vec);
754: VecDestroy(vec);
755: PetscFree(vec);
756: }
757: }
759: if (pc->dm) {
760: /* need to tell all the coarser levels to rebuild the matrix using the DM for that level */
761: for (i=0; i<n-1; i++) {
762: if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
763: }
764: }
766: for (i=1; i<n; i++) {
767: if (mglevels[i]->smoothu == mglevels[i]->smoothd || mg->am == PC_MG_FULL || mg->am == PC_MG_KASKADE || mg->cyclesperpcapply > 1){
768: /* if doing only down then initial guess is zero */
769: KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);
770: }
771: if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
772: KSPSetUp(mglevels[i]->smoothd);
773: if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
774: if (!mglevels[i]->residual) {
775: Mat mat;
776: KSPGetOperators(mglevels[i]->smoothd,NULL,&mat);
777: PCMGSetResidual(pc,i,PCMGResidualDefault,mat);
778: }
779: }
780: for (i=1; i<n; i++) {
781: if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) {
782: Mat downmat,downpmat;
784: /* check if operators have been set for up, if not use down operators to set them */
785: KSPGetOperatorsSet(mglevels[i]->smoothu,&opsset,NULL);
786: if (!opsset) {
787: KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat);
788: KSPSetOperators(mglevels[i]->smoothu,downmat,downpmat);
789: }
791: KSPSetInitialGuessNonzero(mglevels[i]->smoothu,PETSC_TRUE);
792: if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
793: KSPSetUp(mglevels[i]->smoothu);
794: if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
795: }
796: }
798: /*
799: If coarse solver is not direct method then DO NOT USE preonly
800: */
801: PetscObjectTypeCompare((PetscObject)mglevels[0]->smoothd,KSPPREONLY,&preonly);
802: if (preonly) {
803: PetscObjectTypeCompare((PetscObject)cpc,PCLU,&lu);
804: PetscObjectTypeCompare((PetscObject)cpc,PCREDUNDANT,&redundant);
805: PetscObjectTypeCompare((PetscObject)cpc,PCCHOLESKY,&cholesky);
806: PetscObjectTypeCompare((PetscObject)cpc,PCSVD,&svd);
807: if (!lu && !redundant && !cholesky && !svd) {
808: KSPSetType(mglevels[0]->smoothd,KSPGMRES);
809: }
810: }
812: if (mglevels[0]->eventsmoothsetup) {PetscLogEventBegin(mglevels[0]->eventsmoothsetup,0,0,0,0);}
813: KSPSetUp(mglevels[0]->smoothd);
814: if (mglevels[0]->eventsmoothsetup) {PetscLogEventEnd(mglevels[0]->eventsmoothsetup,0,0,0,0);}
816: /*
817: Dump the interpolation/restriction matrices plus the
818: Jacobian/stiffness on each level. This allows MATLAB users to
819: easily check if the Galerkin condition A_c = R A_f R^T is satisfied.
821: Only support one or the other at the same time.
822: */
823: #if defined(PETSC_USE_SOCKET_VIEWER)
824: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_mg_dump_matlab",&dump,NULL);
825: if (dump) viewer = PETSC_VIEWER_SOCKET_(PetscObjectComm((PetscObject)pc));
826: dump = PETSC_FALSE;
827: #endif
828: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_mg_dump_binary",&dump,NULL);
829: if (dump) viewer = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)pc));
831: if (viewer) {
832: for (i=1; i<n; i++) {
833: MatView(mglevels[i]->restrct,viewer);
834: }
835: for (i=0; i<n; i++) {
836: KSPGetPC(mglevels[i]->smoothd,&pc);
837: MatView(pc->mat,viewer);
838: }
839: }
840: return(0);
841: }
843: /* -------------------------------------------------------------------------------------*/
847: /*@
848: PCMGGetLevels - Gets the number of levels to use with MG.
850: Not Collective
852: Input Parameter:
853: . pc - the preconditioner context
855: Output parameter:
856: . levels - the number of levels
858: Level: advanced
860: .keywords: MG, get, levels, multigrid
862: .seealso: PCMGSetLevels()
863: @*/
864: PetscErrorCode PCMGGetLevels(PC pc,PetscInt *levels)
865: {
866: PC_MG *mg = (PC_MG*)pc->data;
871: *levels = mg->nlevels;
872: return(0);
873: }
877: /*@
878: PCMGSetType - Determines the form of multigrid to use:
879: multiplicative, additive, full, or the Kaskade algorithm.
881: Logically Collective on PC
883: Input Parameters:
884: + pc - the preconditioner context
885: - form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,
886: PC_MG_FULL, PC_MG_KASKADE
888: Options Database Key:
889: . -pc_mg_type <form> - Sets <form>, one of multiplicative,
890: additive, full, kaskade
892: Level: advanced
894: .keywords: MG, set, method, multiplicative, additive, full, Kaskade, multigrid
896: .seealso: PCMGSetLevels()
897: @*/
898: PetscErrorCode PCMGSetType(PC pc,PCMGType form)
899: {
900: PC_MG *mg = (PC_MG*)pc->data;
905: mg->am = form;
906: if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG;
907: else pc->ops->applyrichardson = NULL;
908: return(0);
909: }
913: /*@
914: PCMGGetType - Determines the form of multigrid to use:
915: multiplicative, additive, full, or the Kaskade algorithm.
917: Logically Collective on PC
919: Input Parameter:
920: . pc - the preconditioner context
922: Output Parameter:
923: . type - one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,PC_MG_FULL, PC_MG_KASKADE
926: Level: advanced
928: .keywords: MG, set, method, multiplicative, additive, full, Kaskade, multigrid
930: .seealso: PCMGSetLevels()
931: @*/
932: PetscErrorCode PCMGGetType(PC pc,PCMGType *type)
933: {
934: PC_MG *mg = (PC_MG*)pc->data;
938: *type = mg->am;
939: return(0);
940: }
944: /*@
945: PCMGSetCycleType - Sets the type cycles to use. Use PCMGSetCycleTypeOnLevel() for more
946: complicated cycling.
948: Logically Collective on PC
950: Input Parameters:
951: + pc - the multigrid context
952: - PC_MG_CYCLE_V or PC_MG_CYCLE_W
954: Options Database Key:
955: . -pc_mg_cycle_type <v,w>
957: Level: advanced
959: .keywords: MG, set, cycles, V-cycle, W-cycle, multigrid
961: .seealso: PCMGSetCycleTypeOnLevel()
962: @*/
963: PetscErrorCode PCMGSetCycleType(PC pc,PCMGCycleType n)
964: {
965: PC_MG *mg = (PC_MG*)pc->data;
966: PC_MG_Levels **mglevels = mg->levels;
967: PetscInt i,levels;
971: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
973: levels = mglevels[0]->levels;
975: for (i=0; i<levels; i++) mglevels[i]->cycles = n;
976: return(0);
977: }
981: /*@
982: PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step
983: of multigrid when PCMGType of PC_MG_MULTIPLICATIVE is used
985: Logically Collective on PC
987: Input Parameters:
988: + pc - the multigrid context
989: - n - number of cycles (default is 1)
991: Options Database Key:
992: . -pc_mg_multiplicative_cycles n
994: Level: advanced
996: Notes: This is not associated with setting a v or w cycle, that is set with PCMGSetCycleType()
998: .keywords: MG, set, cycles, V-cycle, W-cycle, multigrid
1000: .seealso: PCMGSetCycleTypeOnLevel(), PCMGSetCycleType()
1001: @*/
1002: PetscErrorCode PCMGMultiplicativeSetCycles(PC pc,PetscInt n)
1003: {
1004: PC_MG *mg = (PC_MG*)pc->data;
1009: mg->cyclesperpcapply = n;
1010: return(0);
1011: }
1015: PetscErrorCode PCMGSetGalerkin_MG(PC pc,PetscBool use)
1016: {
1017: PC_MG *mg = (PC_MG*)pc->data;
1020: mg->galerkin = use ? 1 : 0;
1021: return(0);
1022: }
1026: /*@
1027: PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the
1028: finest grid via the Galerkin process: A_i-1 = r_i * A_i * p_i
1030: Logically Collective on PC
1032: Input Parameters:
1033: + pc - the multigrid context
1034: - use - PETSC_TRUE to use the Galerkin process to compute coarse-level operators
1036: Options Database Key:
1037: . -pc_mg_galerkin <true,false>
1039: Level: intermediate
1041: Notes: Some codes that use PCMG such as PCGAMG use Galerkin internally while constructing the hierarchy and thus do not
1042: use the PCMG construction of the coarser grids.
1044: .keywords: MG, set, Galerkin
1046: .seealso: PCMGGetGalerkin()
1048: @*/
1049: PetscErrorCode PCMGSetGalerkin(PC pc,PetscBool use)
1050: {
1055: PetscTryMethod(pc,"PCMGSetGalerkin_C",(PC,PetscBool),(pc,use));
1056: return(0);
1057: }
1061: /*@
1062: PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e.
1063: A_i-1 = r_i * A_i * p_i
1065: Not Collective
1067: Input Parameter:
1068: . pc - the multigrid context
1070: Output Parameter:
1071: . galerkin - PETSC_TRUE or PETSC_FALSE
1073: Options Database Key:
1074: . -pc_mg_galerkin
1076: Level: intermediate
1078: .keywords: MG, set, Galerkin
1080: .seealso: PCMGSetGalerkin()
1082: @*/
1083: PetscErrorCode PCMGGetGalerkin(PC pc,PetscBool *galerkin)
1084: {
1085: PC_MG *mg = (PC_MG*)pc->data;
1089: *galerkin = (PetscBool)mg->galerkin;
1090: return(0);
1091: }
1095: /*@
1096: PCMGSetNumberSmoothDown - Sets the number of pre-smoothing steps to
1097: use on all levels. Use PCMGGetSmootherDown() to set different
1098: pre-smoothing steps on different levels.
1100: Logically Collective on PC
1102: Input Parameters:
1103: + mg - the multigrid context
1104: - n - the number of smoothing steps
1106: Options Database Key:
1107: . -pc_mg_smoothdown <n> - Sets number of pre-smoothing steps
1109: Level: advanced
1111: .keywords: MG, smooth, down, pre-smoothing, steps, multigrid
1113: .seealso: PCMGSetNumberSmoothUp()
1114: @*/
1115: PetscErrorCode PCMGSetNumberSmoothDown(PC pc,PetscInt n)
1116: {
1117: PC_MG *mg = (PC_MG*)pc->data;
1118: PC_MG_Levels **mglevels = mg->levels;
1120: PetscInt i,levels;
1124: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
1126: levels = mglevels[0]->levels;
1128: for (i=1; i<levels; i++) {
1129: /* make sure smoother up and down are different */
1130: PCMGGetSmootherUp(pc,i,NULL);
1131: KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
1133: mg->default_smoothd = n;
1134: }
1135: return(0);
1136: }
1140: /*@
1141: PCMGSetNumberSmoothUp - Sets the number of post-smoothing steps to use
1142: on all levels. Use PCMGGetSmootherUp() to set different numbers of
1143: post-smoothing steps on different levels.
1145: Logically Collective on PC
1147: Input Parameters:
1148: + mg - the multigrid context
1149: - n - the number of smoothing steps
1151: Options Database Key:
1152: . -pc_mg_smoothup <n> - Sets number of post-smoothing steps
1154: Level: advanced
1156: Note: this does not set a value on the coarsest grid, since we assume that
1157: there is no separate smooth up on the coarsest grid.
1159: .keywords: MG, smooth, up, post-smoothing, steps, multigrid
1161: .seealso: PCMGSetNumberSmoothDown()
1162: @*/
1163: PetscErrorCode PCMGSetNumberSmoothUp(PC pc,PetscInt n)
1164: {
1165: PC_MG *mg = (PC_MG*)pc->data;
1166: PC_MG_Levels **mglevels = mg->levels;
1168: PetscInt i,levels;
1172: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
1174: levels = mglevels[0]->levels;
1176: for (i=1; i<levels; i++) {
1177: /* make sure smoother up and down are different */
1178: PCMGGetSmootherUp(pc,i,NULL);
1179: KSPSetTolerances(mglevels[i]->smoothu,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
1181: mg->default_smoothu = n;
1182: }
1183: return(0);
1184: }
1186: /* ----------------------------------------------------------------------------------------*/
1188: /*MC
1189: PCMG - Use multigrid preconditioning. This preconditioner requires you provide additional
1190: information about the coarser grid matrices and restriction/interpolation operators.
1192: Options Database Keys:
1193: + -pc_mg_levels <nlevels> - number of levels including finest
1194: . -pc_mg_cycles <v,w> -
1195: . -pc_mg_smoothup <n> - number of smoothing steps after interpolation
1196: . -pc_mg_smoothdown <n> - number of smoothing steps before applying restriction operator
1197: . -pc_mg_type <additive,multiplicative,full,kaskade> - multiplicative is the default
1198: . -pc_mg_log - log information about time spent on each level of the solver
1199: . -pc_mg_galerkin - use Galerkin process to compute coarser operators, i.e. Acoarse = R A R'
1200: . -pc_mg_multiplicative_cycles - number of cycles to use as the preconditioner (defaults to 1)
1201: . -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices
1202: to the Socket viewer for reading from MATLAB.
1203: - -pc_mg_dump_binary - dumps the matrices for each level and the restriction/interpolation matrices
1204: to the binary output file called binaryoutput
1206: 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
1208: When run with a single level the smoother options are used on that level NOT the coarse grid solver options
1210: Level: intermediate
1212: Concepts: multigrid/multilevel
1214: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, PCEXOTIC, PCGAMG, PCML, PCHYPRE
1215: PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), PCMGSetNumberSmoothDown(),
1216: PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(),
1217: PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(),
1218: PCMGSetCycleTypeOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR()
1219: M*/
1223: PETSC_EXTERN PetscErrorCode PCCreate_MG(PC pc)
1224: {
1225: PC_MG *mg;
1229: PetscNewLog(pc,&mg);
1230: pc->data = (void*)mg;
1231: mg->nlevels = -1;
1232: mg->am = PC_MG_MULTIPLICATIVE;
1234: pc->useAmat = PETSC_TRUE;
1236: pc->ops->apply = PCApply_MG;
1237: pc->ops->setup = PCSetUp_MG;
1238: pc->ops->reset = PCReset_MG;
1239: pc->ops->destroy = PCDestroy_MG;
1240: pc->ops->setfromoptions = PCSetFromOptions_MG;
1241: pc->ops->view = PCView_MG;
1243: PetscObjectComposeFunction((PetscObject)pc,"PCMGSetGalerkin_C",PCMGSetGalerkin_MG);
1244: return(0);
1245: }