Actual source code: mg.c
petsc-3.7.5 2017-01-01
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;
16: PC subpc;
17: PCFailedReason pcreason;
20: if (mglevels->eventsmoothsolve) {PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);}
21: KSPSolve(mglevels->smoothd,mglevels->b,mglevels->x); /* pre-smooth */
22: KSPGetPC(mglevels->smoothd,&subpc);
23: PCGetSetUpFailedReason(subpc,&pcreason);
24: if (pcreason) {
25: pc->failedreason = PC_SUBPC_ERROR;
26: }
27: if (mglevels->eventsmoothsolve) {PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);}
28: if (mglevels->level) { /* not the coarsest grid */
29: if (mglevels->eventresidual) {PetscLogEventBegin(mglevels->eventresidual,0,0,0,0);}
30: (*mglevels->residual)(mglevels->A,mglevels->b,mglevels->x,mglevels->r);
31: if (mglevels->eventresidual) {PetscLogEventEnd(mglevels->eventresidual,0,0,0,0);}
33: /* if on finest level and have convergence criteria set */
34: if (mglevels->level == mglevels->levels-1 && mg->ttol && reason) {
35: PetscReal rnorm;
36: VecNorm(mglevels->r,NORM_2,&rnorm);
37: if (rnorm <= mg->ttol) {
38: if (rnorm < mg->abstol) {
39: *reason = PCRICHARDSON_CONVERGED_ATOL;
40: PetscInfo2(pc,"Linear solver has converged. Residual norm %g is less than absolute tolerance %g\n",(double)rnorm,(double)mg->abstol);
41: } else {
42: *reason = PCRICHARDSON_CONVERGED_RTOL;
43: 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);
44: }
45: return(0);
46: }
47: }
49: mgc = *(mglevelsin - 1);
50: if (mglevels->eventinterprestrict) {PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);}
51: MatRestrict(mglevels->restrct,mglevels->r,mgc->b);
52: if (mglevels->eventinterprestrict) {PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);}
53: VecSet(mgc->x,0.0);
54: while (cycles--) {
55: PCMGMCycle_Private(pc,mglevelsin-1,reason);
56: }
57: if (mglevels->eventinterprestrict) {PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);}
58: MatInterpolateAdd(mglevels->interpolate,mgc->x,mglevels->x,mglevels->x);
59: if (mglevels->eventinterprestrict) {PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);}
60: if (mglevels->eventsmoothsolve) {PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);}
61: KSPSolve(mglevels->smoothu,mglevels->b,mglevels->x); /* post smooth */
62: if (mglevels->eventsmoothsolve) {PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);}
63: }
64: return(0);
65: }
69: 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)
70: {
71: PC_MG *mg = (PC_MG*)pc->data;
72: PC_MG_Levels **mglevels = mg->levels;
74: PetscInt levels = mglevels[0]->levels,i;
77: /* When the DM is supplying the matrix then it will not exist until here */
78: for (i=0; i<levels; i++) {
79: if (!mglevels[i]->A) {
80: KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);
81: PetscObjectReference((PetscObject)mglevels[i]->A);
82: }
83: }
84: mglevels[levels-1]->b = b;
85: mglevels[levels-1]->x = x;
87: mg->rtol = rtol;
88: mg->abstol = abstol;
89: mg->dtol = dtol;
90: if (rtol) {
91: /* compute initial residual norm for relative convergence test */
92: PetscReal rnorm;
93: if (zeroguess) {
94: VecNorm(b,NORM_2,&rnorm);
95: } else {
96: (*mglevels[levels-1]->residual)(mglevels[levels-1]->A,b,x,w);
97: VecNorm(w,NORM_2,&rnorm);
98: }
99: mg->ttol = PetscMax(rtol*rnorm,abstol);
100: } else if (abstol) mg->ttol = abstol;
101: else mg->ttol = 0.0;
103: /* since smoother is applied to full system, not just residual we need to make sure that smoothers don't
104: stop prematurely due to small residual */
105: for (i=1; i<levels; i++) {
106: KSPSetTolerances(mglevels[i]->smoothu,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
107: if (mglevels[i]->smoothu != mglevels[i]->smoothd) {
108: KSPSetTolerances(mglevels[i]->smoothd,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
109: }
110: }
112: *reason = (PCRichardsonConvergedReason)0;
113: for (i=0; i<its; i++) {
114: PCMGMCycle_Private(pc,mglevels+levels-1,reason);
115: if (*reason) break;
116: }
117: if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
118: *outits = i;
119: return(0);
120: }
124: PetscErrorCode PCReset_MG(PC pc)
125: {
126: PC_MG *mg = (PC_MG*)pc->data;
127: PC_MG_Levels **mglevels = mg->levels;
129: PetscInt i,n;
132: if (mglevels) {
133: n = mglevels[0]->levels;
134: for (i=0; i<n-1; i++) {
135: VecDestroy(&mglevels[i+1]->r);
136: VecDestroy(&mglevels[i]->b);
137: VecDestroy(&mglevels[i]->x);
138: MatDestroy(&mglevels[i+1]->restrct);
139: MatDestroy(&mglevels[i+1]->interpolate);
140: VecDestroy(&mglevels[i+1]->rscale);
141: }
143: for (i=0; i<n; i++) {
144: MatDestroy(&mglevels[i]->A);
145: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
146: KSPReset(mglevels[i]->smoothd);
147: }
148: KSPReset(mglevels[i]->smoothu);
149: }
150: }
151: return(0);
152: }
156: /*@C
157: PCMGSetLevels - Sets the number of levels to use with MG.
158: Must be called before any other MG routine.
160: Logically Collective on PC
162: Input Parameters:
163: + pc - the preconditioner context
164: . levels - the number of levels
165: - comms - optional communicators for each level; this is to allow solving the coarser problems
166: on smaller sets of processors. Use NULL_OBJECT for default in Fortran
168: Level: intermediate
170: Notes:
171: If the number of levels is one then the multigrid uses the -mg_levels prefix
172: for setting the level options rather than the -mg_coarse prefix.
174: .keywords: MG, set, levels, multigrid
176: .seealso: PCMGSetType(), PCMGGetLevels()
177: @*/
178: PetscErrorCode PCMGSetLevels(PC pc,PetscInt levels,MPI_Comm *comms)
179: {
181: PC_MG *mg = (PC_MG*)pc->data;
182: MPI_Comm comm;
183: PC_MG_Levels **mglevels = mg->levels;
184: PCMGType mgtype = mg->am;
185: PetscInt mgctype = (PetscInt) PC_MG_CYCLE_V;
186: PetscInt i;
187: PetscMPIInt size;
188: const char *prefix;
189: PC ipc;
190: PetscInt n;
195: PetscObjectGetComm((PetscObject)pc,&comm);
196: if (mg->nlevels == levels) return(0);
197: if (mglevels) {
198: mgctype = mglevels[0]->cycles;
199: /* changing the number of levels so free up the previous stuff */
200: PCReset_MG(pc);
201: n = mglevels[0]->levels;
202: for (i=0; i<n; i++) {
203: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
204: KSPDestroy(&mglevels[i]->smoothd);
205: }
206: KSPDestroy(&mglevels[i]->smoothu);
207: PetscFree(mglevels[i]);
208: }
209: PetscFree(mg->levels);
210: }
212: mg->nlevels = levels;
214: PetscMalloc1(levels,&mglevels);
215: PetscLogObjectMemory((PetscObject)pc,levels*(sizeof(PC_MG*)));
217: PCGetOptionsPrefix(pc,&prefix);
219: mg->stageApply = 0;
220: for (i=0; i<levels; i++) {
221: PetscNewLog(pc,&mglevels[i]);
223: mglevels[i]->level = i;
224: mglevels[i]->levels = levels;
225: mglevels[i]->cycles = mgctype;
226: mg->default_smoothu = 2;
227: mg->default_smoothd = 2;
228: mglevels[i]->eventsmoothsetup = 0;
229: mglevels[i]->eventsmoothsolve = 0;
230: mglevels[i]->eventresidual = 0;
231: mglevels[i]->eventinterprestrict = 0;
233: if (comms) comm = comms[i];
234: KSPCreate(comm,&mglevels[i]->smoothd);
235: KSPSetErrorIfNotConverged(mglevels[i]->smoothd,pc->erroriffailure);
236: PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd,(PetscObject)pc,levels-i);
237: KSPSetOptionsPrefix(mglevels[i]->smoothd,prefix);
238: PetscObjectComposedDataSetInt((PetscObject) mglevels[i]->smoothd, PetscMGLevelId, mglevels[i]->level);
239: if (i || levels == 1) {
240: char tprefix[128];
242: KSPSetType(mglevels[i]->smoothd,KSPCHEBYSHEV);
243: KSPSetConvergenceTest(mglevels[i]->smoothd,KSPConvergedSkip,NULL,NULL);
244: KSPSetNormType(mglevels[i]->smoothd,KSP_NORM_NONE);
245: KSPGetPC(mglevels[i]->smoothd,&ipc);
246: PCSetType(ipc,PCSOR);
247: KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, mg->default_smoothd);
249: sprintf(tprefix,"mg_levels_%d_",(int)i);
250: KSPAppendOptionsPrefix(mglevels[i]->smoothd,tprefix);
251: } else {
252: KSPAppendOptionsPrefix(mglevels[0]->smoothd,"mg_coarse_");
254: /* coarse solve is (redundant) LU by default; set shifttype NONZERO to avoid annoying zero-pivot in LU preconditioner */
255: KSPSetType(mglevels[0]->smoothd,KSPPREONLY);
256: KSPGetPC(mglevels[0]->smoothd,&ipc);
257: MPI_Comm_size(comm,&size);
258: if (size > 1) {
259: PCSetType(ipc,PCREDUNDANT);
260: } else {
261: PCSetType(ipc,PCLU);
262: }
263: PCFactorSetShiftType(ipc,MAT_SHIFT_INBLOCKS);
264: }
265: PetscLogObjectParent((PetscObject)pc,(PetscObject)mglevels[i]->smoothd);
267: mglevels[i]->smoothu = mglevels[i]->smoothd;
268: mg->rtol = 0.0;
269: mg->abstol = 0.0;
270: mg->dtol = 0.0;
271: mg->ttol = 0.0;
272: mg->cyclesperpcapply = 1;
273: }
274: mg->levels = mglevels;
275: PCMGSetType(pc,mgtype);
276: return(0);
277: }
282: PetscErrorCode PCDestroy_MG(PC pc)
283: {
285: PC_MG *mg = (PC_MG*)pc->data;
286: PC_MG_Levels **mglevels = mg->levels;
287: PetscInt i,n;
290: PCReset_MG(pc);
291: if (mglevels) {
292: n = mglevels[0]->levels;
293: for (i=0; i<n; i++) {
294: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
295: KSPDestroy(&mglevels[i]->smoothd);
296: }
297: KSPDestroy(&mglevels[i]->smoothu);
298: PetscFree(mglevels[i]);
299: }
300: PetscFree(mg->levels);
301: }
302: PetscFree(pc->data);
303: return(0);
304: }
308: extern PetscErrorCode PCMGACycle_Private(PC,PC_MG_Levels**);
309: extern PetscErrorCode PCMGFCycle_Private(PC,PC_MG_Levels**);
310: extern PetscErrorCode PCMGKCycle_Private(PC,PC_MG_Levels**);
312: /*
313: PCApply_MG - Runs either an additive, multiplicative, Kaskadic
314: or full cycle of multigrid.
316: Note:
317: A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle().
318: */
321: static PetscErrorCode PCApply_MG(PC pc,Vec b,Vec x)
322: {
323: PC_MG *mg = (PC_MG*)pc->data;
324: PC_MG_Levels **mglevels = mg->levels;
326: PetscInt levels = mglevels[0]->levels,i;
329: if (mg->stageApply) {PetscLogStagePush(mg->stageApply);}
330: /* When the DM is supplying the matrix then it will not exist until here */
331: for (i=0; i<levels; i++) {
332: if (!mglevels[i]->A) {
333: KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,NULL);
334: PetscObjectReference((PetscObject)mglevels[i]->A);
335: }
336: }
338: mglevels[levels-1]->b = b;
339: mglevels[levels-1]->x = x;
340: if (mg->am == PC_MG_MULTIPLICATIVE) {
341: VecSet(x,0.0);
342: for (i=0; i<mg->cyclesperpcapply; i++) {
343: PCMGMCycle_Private(pc,mglevels+levels-1,NULL);
344: }
345: } else if (mg->am == PC_MG_ADDITIVE) {
346: PCMGACycle_Private(pc,mglevels);
347: } else if (mg->am == PC_MG_KASKADE) {
348: PCMGKCycle_Private(pc,mglevels);
349: } else {
350: PCMGFCycle_Private(pc,mglevels);
351: }
352: if (mg->stageApply) {PetscLogStagePop();}
353: return(0);
354: }
359: PetscErrorCode PCSetFromOptions_MG(PetscOptionItems *PetscOptionsObject,PC pc)
360: {
362: PetscInt m,levels = 1,cycles;
363: PetscBool flg,set;
364: PC_MG *mg = (PC_MG*)pc->data;
365: PC_MG_Levels **mglevels;
366: PCMGType mgtype;
367: PCMGCycleType mgctype;
370: PetscOptionsHead(PetscOptionsObject,"Multigrid options");
371: if (!mg->levels) {
372: PetscOptionsInt("-pc_mg_levels","Number of Levels","PCMGSetLevels",levels,&levels,&flg);
373: if (!flg && pc->dm) {
374: DMGetRefineLevel(pc->dm,&levels);
375: levels++;
376: mg->usedmfornumberoflevels = PETSC_TRUE;
377: }
378: PCMGSetLevels(pc,levels,NULL);
379: }
380: mglevels = mg->levels;
382: mgctype = (PCMGCycleType) mglevels[0]->cycles;
383: PetscOptionsEnum("-pc_mg_cycle_type","V cycle or for W-cycle","PCMGSetCycleType",PCMGCycleTypes,(PetscEnum)mgctype,(PetscEnum*)&mgctype,&flg);
384: if (flg) {
385: PCMGSetCycleType(pc,mgctype);
386: }
387: flg = PETSC_FALSE;
388: PetscOptionsBool("-pc_mg_galerkin","Use Galerkin process to compute coarser operators","PCMGSetGalerkin",flg,&flg,&set);
389: if (set) {
390: PCMGSetGalerkin(pc,flg);
391: }
392: PetscOptionsInt("-pc_mg_smoothup","Number of post-smoothing steps","PCMGSetNumberSmoothUp",mg->default_smoothu,&m,&flg);
393: if (flg) {
394: PCMGSetNumberSmoothUp(pc,m);
395: }
396: PetscOptionsInt("-pc_mg_smoothdown","Number of pre-smoothing steps","PCMGSetNumberSmoothDown",mg->default_smoothd,&m,&flg);
397: if (flg) {
398: PCMGSetNumberSmoothDown(pc,m);
399: }
400: mgtype = mg->am;
401: PetscOptionsEnum("-pc_mg_type","Multigrid type","PCMGSetType",PCMGTypes,(PetscEnum)mgtype,(PetscEnum*)&mgtype,&flg);
402: if (flg) {
403: PCMGSetType(pc,mgtype);
404: }
405: if (mg->am == PC_MG_MULTIPLICATIVE) {
406: PetscOptionsInt("-pc_mg_multiplicative_cycles","Number of cycles for each preconditioner step","PCMGMultiplicativeSetCycles",mg->cyclesperpcapply,&cycles,&flg);
407: if (flg) {
408: PCMGMultiplicativeSetCycles(pc,cycles);
409: }
410: }
411: flg = PETSC_FALSE;
412: PetscOptionsBool("-pc_mg_log","Log times for each multigrid level","None",flg,&flg,NULL);
413: if (flg) {
414: PetscInt i;
415: char eventname[128];
416: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
417: levels = mglevels[0]->levels;
418: for (i=0; i<levels; i++) {
419: sprintf(eventname,"MGSetup Level %d",(int)i);
420: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsetup);
421: sprintf(eventname,"MGSmooth Level %d",(int)i);
422: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsolve);
423: if (i) {
424: sprintf(eventname,"MGResid Level %d",(int)i);
425: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventresidual);
426: sprintf(eventname,"MGInterp Level %d",(int)i);
427: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventinterprestrict);
428: }
429: }
431: #if defined(PETSC_USE_LOG)
432: {
433: const char *sname = "MG Apply";
434: PetscStageLog stageLog;
435: PetscInt st;
438: PetscLogGetStageLog(&stageLog);
439: for (st = 0; st < stageLog->numStages; ++st) {
440: PetscBool same;
442: PetscStrcmp(stageLog->stageInfo[st].name, sname, &same);
443: if (same) mg->stageApply = st;
444: }
445: if (!mg->stageApply) {
446: PetscLogStageRegister(sname, &mg->stageApply);
447: }
448: }
449: #endif
450: }
451: PetscOptionsTail();
452: return(0);
453: }
455: const char *const PCMGTypes[] = {"MULTIPLICATIVE","ADDITIVE","FULL","KASKADE","PCMGType","PC_MG",0};
456: const char *const PCMGCycleTypes[] = {"invalid","v","w","PCMGCycleType","PC_MG_CYCLE",0};
458: #include <petscdraw.h>
461: PetscErrorCode PCView_MG(PC pc,PetscViewer viewer)
462: {
463: PC_MG *mg = (PC_MG*)pc->data;
464: PC_MG_Levels **mglevels = mg->levels;
466: PetscInt levels = mglevels ? mglevels[0]->levels : 0,i;
467: PetscBool iascii,isbinary,isdraw;
470: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
471: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERBINARY,&isbinary);
472: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERDRAW,&isdraw);
473: if (iascii) {
474: const char *cyclename = levels ? (mglevels[0]->cycles == PC_MG_CYCLE_V ? "v" : "w") : "unknown";
475: PetscViewerASCIIPrintf(viewer," MG: type is %s, levels=%D cycles=%s\n", PCMGTypes[mg->am],levels,cyclename);
476: if (mg->am == PC_MG_MULTIPLICATIVE) {
477: PetscViewerASCIIPrintf(viewer," Cycles per PCApply=%d\n",mg->cyclesperpcapply);
478: }
479: if (mg->galerkin) {
480: PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices\n");
481: } else {
482: PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid matrices\n");
483: }
484: if (mg->view){
485: (*mg->view)(pc,viewer);
486: }
487: for (i=0; i<levels; i++) {
488: if (!i) {
489: PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level -------------------------------\n",i);
490: } else {
491: PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);
492: }
493: PetscViewerASCIIPushTab(viewer);
494: KSPView(mglevels[i]->smoothd,viewer);
495: PetscViewerASCIIPopTab(viewer);
496: if (i && mglevels[i]->smoothd == mglevels[i]->smoothu) {
497: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");
498: } else if (i) {
499: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);
500: PetscViewerASCIIPushTab(viewer);
501: KSPView(mglevels[i]->smoothu,viewer);
502: PetscViewerASCIIPopTab(viewer);
503: }
504: }
505: } else if (isbinary) {
506: for (i=levels-1; i>=0; i--) {
507: KSPView(mglevels[i]->smoothd,viewer);
508: if (i && mglevels[i]->smoothd != mglevels[i]->smoothu) {
509: KSPView(mglevels[i]->smoothu,viewer);
510: }
511: }
512: } else if (isdraw) {
513: PetscDraw draw;
514: PetscReal x,w,y,bottom,th;
515: PetscViewerDrawGetDraw(viewer,0,&draw);
516: PetscDrawGetCurrentPoint(draw,&x,&y);
517: PetscDrawStringGetSize(draw,NULL,&th);
518: bottom = y - th;
519: for (i=levels-1; i>=0; i--) {
520: if (!mglevels[i]->smoothu || (mglevels[i]->smoothu == mglevels[i]->smoothd)) {
521: PetscDrawPushCurrentPoint(draw,x,bottom);
522: KSPView(mglevels[i]->smoothd,viewer);
523: PetscDrawPopCurrentPoint(draw);
524: } else {
525: w = 0.5*PetscMin(1.0-x,x);
526: PetscDrawPushCurrentPoint(draw,x+w,bottom);
527: KSPView(mglevels[i]->smoothd,viewer);
528: PetscDrawPopCurrentPoint(draw);
529: PetscDrawPushCurrentPoint(draw,x-w,bottom);
530: KSPView(mglevels[i]->smoothu,viewer);
531: PetscDrawPopCurrentPoint(draw);
532: }
533: PetscDrawGetBoundingBox(draw,NULL,&bottom,NULL,NULL);
534: bottom -= th;
535: }
536: }
537: return(0);
538: }
540: #include <petsc/private/dmimpl.h>
541: #include <petsc/private/kspimpl.h>
543: /*
544: Calls setup for the KSP on each level
545: */
548: PetscErrorCode PCSetUp_MG(PC pc)
549: {
550: PC_MG *mg = (PC_MG*)pc->data;
551: PC_MG_Levels **mglevels = mg->levels;
553: PetscInt i,n = mglevels[0]->levels;
554: PC cpc;
555: PetscBool dump = PETSC_FALSE,opsset,use_amat,missinginterpolate = PETSC_FALSE;
556: Mat dA,dB;
557: Vec tvec;
558: DM *dms;
559: PetscViewer viewer = 0;
562: /* FIX: Move this to PCSetFromOptions_MG? */
563: if (mg->usedmfornumberoflevels) {
564: PetscInt levels;
565: DMGetRefineLevel(pc->dm,&levels);
566: levels++;
567: if (levels > n) { /* the problem is now being solved on a finer grid */
568: PCMGSetLevels(pc,levels,NULL);
569: n = levels;
570: PCSetFromOptions(pc); /* it is bad to call this here, but otherwise will never be called for the new hierarchy */
571: mglevels = mg->levels;
572: }
573: }
574: KSPGetPC(mglevels[0]->smoothd,&cpc);
577: /* If user did not provide fine grid operators OR operator was not updated since last global KSPSetOperators() */
578: /* so use those from global PC */
579: /* Is this what we always want? What if user wants to keep old one? */
580: KSPGetOperatorsSet(mglevels[n-1]->smoothd,NULL,&opsset);
581: if (opsset) {
582: Mat mmat;
583: KSPGetOperators(mglevels[n-1]->smoothd,NULL,&mmat);
584: if (mmat == pc->pmat) opsset = PETSC_FALSE;
585: }
587: if (!opsset) {
588: PCGetUseAmat(pc,&use_amat);
589: if(use_amat){
590: PetscInfo(pc,"Using outer operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
591: KSPSetOperators(mglevels[n-1]->smoothd,pc->mat,pc->pmat);
592: }
593: else {
594: PetscInfo(pc,"Using matrix (pmat) operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
595: KSPSetOperators(mglevels[n-1]->smoothd,pc->pmat,pc->pmat);
596: }
597: }
599: for (i=n-1; i>0; i--) {
600: if (!(mglevels[i]->interpolate || mglevels[i]->restrct)) {
601: missinginterpolate = PETSC_TRUE;
602: continue;
603: }
604: }
605: /*
606: Skipping if user has provided all interpolation/restriction needed (since DM might not be able to produce them (when coming from SNES/TS)
607: Skipping for galerkin==2 (externally managed hierarchy such as ML and GAMG). Cleaner logic here would be great. Wrap ML/GAMG as DMs?
608: */
609: if (missinginterpolate && pc->dm && mg->galerkin != 2 && !pc->setupcalled) {
610: /* construct the interpolation from the DMs */
611: Mat p;
612: Vec rscale;
613: PetscMalloc1(n,&dms);
614: dms[n-1] = pc->dm;
615: /* Separately create them so we do not get DMKSP interference between levels */
616: for (i=n-2; i>-1; i--) {DMCoarsen(dms[i+1],MPI_COMM_NULL,&dms[i]);}
617: for (i=n-2; i>-1; i--) {
618: DMKSP kdm;
619: PetscBool dmhasrestrict;
620: KSPSetDM(mglevels[i]->smoothd,dms[i]);
621: if (mg->galerkin) {KSPSetDMActive(mglevels[i]->smoothd,PETSC_FALSE);}
622: DMGetDMKSPWrite(dms[i],&kdm);
623: /* Ugly hack so that the next KSPSetUp() will use the RHS that we set. A better fix is to change dmActive to take
624: * a bitwise OR of computing the matrix, RHS, and initial iterate. */
625: kdm->ops->computerhs = NULL;
626: kdm->rhsctx = NULL;
627: if (!mglevels[i+1]->interpolate) {
628: DMCreateInterpolation(dms[i],dms[i+1],&p,&rscale);
629: PCMGSetInterpolation(pc,i+1,p);
630: if (rscale) {PCMGSetRScale(pc,i+1,rscale);}
631: VecDestroy(&rscale);
632: MatDestroy(&p);
633: }
634: DMHasCreateRestriction(dms[i],&dmhasrestrict);
635: if (dmhasrestrict && !mglevels[i+1]->restrct){
636: DMCreateRestriction(dms[i],dms[i+1],&p);
637: PCMGSetRestriction(pc,i+1,p);
638: MatDestroy(&p);
639: }
640: }
642: for (i=n-2; i>-1; i--) {DMDestroy(&dms[i]);}
643: PetscFree(dms);
644: }
646: if (pc->dm && !pc->setupcalled) {
647: /* finest smoother also gets DM but it is not active, independent of whether galerkin==2 */
648: KSPSetDM(mglevels[n-1]->smoothd,pc->dm);
649: KSPSetDMActive(mglevels[n-1]->smoothd,PETSC_FALSE);
650: }
652: if (mg->galerkin == 1) {
653: Mat B;
654: /* currently only handle case where mat and pmat are the same on coarser levels */
655: KSPGetOperators(mglevels[n-1]->smoothd,&dA,&dB);
656: if (!pc->setupcalled) {
657: for (i=n-2; i>-1; i--) {
658: 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");
659: if (!mglevels[i+1]->interpolate) {
660: PCMGSetInterpolation(pc,i+1,mglevels[i+1]->restrct);
661: }
662: if (!mglevels[i+1]->restrct) {
663: PCMGSetRestriction(pc,i+1,mglevels[i+1]->interpolate);
664: }
665: if (mglevels[i+1]->interpolate == mglevels[i+1]->restrct) {
666: MatPtAP(dB,mglevels[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);
667: } else {
668: MatMatMatMult(mglevels[i+1]->restrct,dB,mglevels[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);
669: }
670: KSPSetOperators(mglevels[i]->smoothd,B,B);
671: if (i != n-2) {PetscObjectDereference((PetscObject)dB);}
672: dB = B;
673: }
674: if (n > 1) {PetscObjectDereference((PetscObject)dB);}
675: } else {
676: for (i=n-2; i>-1; i--) {
677: 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");
678: if (!mglevels[i+1]->interpolate) {
679: PCMGSetInterpolation(pc,i+1,mglevels[i+1]->restrct);
680: }
681: if (!mglevels[i+1]->restrct) {
682: PCMGSetRestriction(pc,i+1,mglevels[i+1]->interpolate);
683: }
684: KSPGetOperators(mglevels[i]->smoothd,NULL,&B);
685: if (mglevels[i+1]->interpolate == mglevels[i+1]->restrct) {
686: MatPtAP(dB,mglevels[i+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);
687: } else {
688: MatMatMatMult(mglevels[i+1]->restrct,dB,mglevels[i+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);
689: }
690: KSPSetOperators(mglevels[i]->smoothd,B,B);
691: dB = B;
692: }
693: }
694: } else if (!mg->galerkin && pc->dm && pc->dm->x) {
695: /* need to restrict Jacobian location to coarser meshes for evaluation */
696: for (i=n-2; i>-1; i--) {
697: Mat R;
698: Vec rscale;
699: if (!mglevels[i]->smoothd->dm->x) {
700: Vec *vecs;
701: KSPCreateVecs(mglevels[i]->smoothd,1,&vecs,0,NULL);
703: mglevels[i]->smoothd->dm->x = vecs[0];
705: PetscFree(vecs);
706: }
707: PCMGGetRestriction(pc,i+1,&R);
708: PCMGGetRScale(pc,i+1,&rscale);
709: MatRestrict(R,mglevels[i+1]->smoothd->dm->x,mglevels[i]->smoothd->dm->x);
710: VecPointwiseMult(mglevels[i]->smoothd->dm->x,mglevels[i]->smoothd->dm->x,rscale);
711: }
712: }
713: if (!mg->galerkin && pc->dm) {
714: for (i=n-2; i>=0; i--) {
715: DM dmfine,dmcoarse;
716: Mat Restrict,Inject;
717: Vec rscale;
718: KSPGetDM(mglevels[i+1]->smoothd,&dmfine);
719: KSPGetDM(mglevels[i]->smoothd,&dmcoarse);
720: PCMGGetRestriction(pc,i+1,&Restrict);
721: PCMGGetRScale(pc,i+1,&rscale);
722: Inject = NULL; /* Callback should create it if it needs Injection */
723: DMRestrict(dmfine,Restrict,rscale,Inject,dmcoarse);
724: }
725: }
727: if (!pc->setupcalled) {
728: for (i=0; i<n; i++) {
729: KSPSetFromOptions(mglevels[i]->smoothd);
730: }
731: for (i=1; i<n; i++) {
732: if (mglevels[i]->smoothu && (mglevels[i]->smoothu != mglevels[i]->smoothd)) {
733: KSPSetFromOptions(mglevels[i]->smoothu);
734: }
735: }
736: /* insure that if either interpolation or restriction is set the other other one is set */
737: for (i=1; i<n; i++) {
738: PCMGGetInterpolation(pc,i,NULL);
739: PCMGGetRestriction(pc,i,NULL);
740: }
741: for (i=0; i<n-1; i++) {
742: if (!mglevels[i]->b) {
743: Vec *vec;
744: KSPCreateVecs(mglevels[i]->smoothd,1,&vec,0,NULL);
745: PCMGSetRhs(pc,i,*vec);
746: VecDestroy(vec);
747: PetscFree(vec);
748: }
749: if (!mglevels[i]->r && i) {
750: VecDuplicate(mglevels[i]->b,&tvec);
751: PCMGSetR(pc,i,tvec);
752: VecDestroy(&tvec);
753: }
754: if (!mglevels[i]->x) {
755: VecDuplicate(mglevels[i]->b,&tvec);
756: PCMGSetX(pc,i,tvec);
757: VecDestroy(&tvec);
758: }
759: }
760: if (n != 1 && !mglevels[n-1]->r) {
761: /* PCMGSetR() on the finest level if user did not supply it */
762: Vec *vec;
763: KSPCreateVecs(mglevels[n-1]->smoothd,1,&vec,0,NULL);
764: PCMGSetR(pc,n-1,*vec);
765: VecDestroy(vec);
766: PetscFree(vec);
767: }
768: }
770: if (pc->dm) {
771: /* need to tell all the coarser levels to rebuild the matrix using the DM for that level */
772: for (i=0; i<n-1; i++) {
773: if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
774: }
775: }
777: for (i=1; i<n; i++) {
778: if (mglevels[i]->smoothu == mglevels[i]->smoothd || mg->am == PC_MG_FULL || mg->am == PC_MG_KASKADE || mg->cyclesperpcapply > 1){
779: /* if doing only down then initial guess is zero */
780: KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);
781: }
782: if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
783: KSPSetUp(mglevels[i]->smoothd);
784: if (mglevels[i]->smoothd->reason == KSP_DIVERGED_PCSETUP_FAILED) {
785: pc->failedreason = PC_SUBPC_ERROR;
786: }
787: if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
788: if (!mglevels[i]->residual) {
789: Mat mat;
790: KSPGetOperators(mglevels[i]->smoothd,NULL,&mat);
791: PCMGSetResidual(pc,i,PCMGResidualDefault,mat);
792: }
793: }
794: for (i=1; i<n; i++) {
795: if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) {
796: Mat downmat,downpmat;
798: /* check if operators have been set for up, if not use down operators to set them */
799: KSPGetOperatorsSet(mglevels[i]->smoothu,&opsset,NULL);
800: if (!opsset) {
801: KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat);
802: KSPSetOperators(mglevels[i]->smoothu,downmat,downpmat);
803: }
805: KSPSetInitialGuessNonzero(mglevels[i]->smoothu,PETSC_TRUE);
806: if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
807: KSPSetUp(mglevels[i]->smoothu);
808: if (mglevels[i]->smoothu->reason == KSP_DIVERGED_PCSETUP_FAILED) {
809: pc->failedreason = PC_SUBPC_ERROR;
810: }
811: if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
812: }
813: }
815: if (mglevels[0]->eventsmoothsetup) {PetscLogEventBegin(mglevels[0]->eventsmoothsetup,0,0,0,0);}
816: KSPSetUp(mglevels[0]->smoothd);
817: if (mglevels[0]->smoothd->reason == KSP_DIVERGED_PCSETUP_FAILED) {
818: pc->failedreason = PC_SUBPC_ERROR;
819: }
820: if (mglevels[0]->eventsmoothsetup) {PetscLogEventEnd(mglevels[0]->eventsmoothsetup,0,0,0,0);}
822: /*
823: Dump the interpolation/restriction matrices plus the
824: Jacobian/stiffness on each level. This allows MATLAB users to
825: easily check if the Galerkin condition A_c = R A_f R^T is satisfied.
827: Only support one or the other at the same time.
828: */
829: #if defined(PETSC_USE_SOCKET_VIEWER)
830: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_mg_dump_matlab",&dump,NULL);
831: if (dump) viewer = PETSC_VIEWER_SOCKET_(PetscObjectComm((PetscObject)pc));
832: dump = PETSC_FALSE;
833: #endif
834: PetscOptionsGetBool(((PetscObject)pc)->options,((PetscObject)pc)->prefix,"-pc_mg_dump_binary",&dump,NULL);
835: if (dump) viewer = PETSC_VIEWER_BINARY_(PetscObjectComm((PetscObject)pc));
837: if (viewer) {
838: for (i=1; i<n; i++) {
839: MatView(mglevels[i]->restrct,viewer);
840: }
841: for (i=0; i<n; i++) {
842: KSPGetPC(mglevels[i]->smoothd,&pc);
843: MatView(pc->mat,viewer);
844: }
845: }
846: return(0);
847: }
849: /* -------------------------------------------------------------------------------------*/
853: /*@
854: PCMGGetLevels - Gets the number of levels to use with MG.
856: Not Collective
858: Input Parameter:
859: . pc - the preconditioner context
861: Output parameter:
862: . levels - the number of levels
864: Level: advanced
866: .keywords: MG, get, levels, multigrid
868: .seealso: PCMGSetLevels()
869: @*/
870: PetscErrorCode PCMGGetLevels(PC pc,PetscInt *levels)
871: {
872: PC_MG *mg = (PC_MG*)pc->data;
877: *levels = mg->nlevels;
878: return(0);
879: }
883: /*@
884: PCMGSetType - Determines the form of multigrid to use:
885: multiplicative, additive, full, or the Kaskade algorithm.
887: Logically Collective on PC
889: Input Parameters:
890: + pc - the preconditioner context
891: - form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,
892: PC_MG_FULL, PC_MG_KASKADE
894: Options Database Key:
895: . -pc_mg_type <form> - Sets <form>, one of multiplicative,
896: additive, full, kaskade
898: Level: advanced
900: .keywords: MG, set, method, multiplicative, additive, full, Kaskade, multigrid
902: .seealso: PCMGSetLevels()
903: @*/
904: PetscErrorCode PCMGSetType(PC pc,PCMGType form)
905: {
906: PC_MG *mg = (PC_MG*)pc->data;
911: mg->am = form;
912: if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG;
913: else pc->ops->applyrichardson = NULL;
914: return(0);
915: }
919: /*@
920: PCMGGetType - Determines the form of multigrid to use:
921: multiplicative, additive, full, or the Kaskade algorithm.
923: Logically Collective on PC
925: Input Parameter:
926: . pc - the preconditioner context
928: Output Parameter:
929: . type - one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,PC_MG_FULL, PC_MG_KASKADE
932: Level: advanced
934: .keywords: MG, set, method, multiplicative, additive, full, Kaskade, multigrid
936: .seealso: PCMGSetLevels()
937: @*/
938: PetscErrorCode PCMGGetType(PC pc,PCMGType *type)
939: {
940: PC_MG *mg = (PC_MG*)pc->data;
944: *type = mg->am;
945: return(0);
946: }
950: /*@
951: PCMGSetCycleType - Sets the type cycles to use. Use PCMGSetCycleTypeOnLevel() for more
952: complicated cycling.
954: Logically Collective on PC
956: Input Parameters:
957: + pc - the multigrid context
958: - PC_MG_CYCLE_V or PC_MG_CYCLE_W
960: Options Database Key:
961: . -pc_mg_cycle_type <v,w>
963: Level: advanced
965: .keywords: MG, set, cycles, V-cycle, W-cycle, multigrid
967: .seealso: PCMGSetCycleTypeOnLevel()
968: @*/
969: PetscErrorCode PCMGSetCycleType(PC pc,PCMGCycleType n)
970: {
971: PC_MG *mg = (PC_MG*)pc->data;
972: PC_MG_Levels **mglevels = mg->levels;
973: PetscInt i,levels;
977: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
979: levels = mglevels[0]->levels;
981: for (i=0; i<levels; i++) mglevels[i]->cycles = n;
982: return(0);
983: }
987: /*@
988: PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step
989: of multigrid when PCMGType of PC_MG_MULTIPLICATIVE is used
991: Logically Collective on PC
993: Input Parameters:
994: + pc - the multigrid context
995: - n - number of cycles (default is 1)
997: Options Database Key:
998: . -pc_mg_multiplicative_cycles n
1000: Level: advanced
1002: Notes: This is not associated with setting a v or w cycle, that is set with PCMGSetCycleType()
1004: .keywords: MG, set, cycles, V-cycle, W-cycle, multigrid
1006: .seealso: PCMGSetCycleTypeOnLevel(), PCMGSetCycleType()
1007: @*/
1008: PetscErrorCode PCMGMultiplicativeSetCycles(PC pc,PetscInt n)
1009: {
1010: PC_MG *mg = (PC_MG*)pc->data;
1015: mg->cyclesperpcapply = n;
1016: return(0);
1017: }
1021: PetscErrorCode PCMGSetGalerkin_MG(PC pc,PetscBool use)
1022: {
1023: PC_MG *mg = (PC_MG*)pc->data;
1026: mg->galerkin = use ? 1 : 0;
1027: return(0);
1028: }
1032: /*@
1033: PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the
1034: finest grid via the Galerkin process: A_i-1 = r_i * A_i * p_i
1036: Logically Collective on PC
1038: Input Parameters:
1039: + pc - the multigrid context
1040: - use - PETSC_TRUE to use the Galerkin process to compute coarse-level operators
1042: Options Database Key:
1043: . -pc_mg_galerkin <true,false>
1045: Level: intermediate
1047: Notes: Some codes that use PCMG such as PCGAMG use Galerkin internally while constructing the hierarchy and thus do not
1048: use the PCMG construction of the coarser grids.
1050: .keywords: MG, set, Galerkin
1052: .seealso: PCMGGetGalerkin()
1054: @*/
1055: PetscErrorCode PCMGSetGalerkin(PC pc,PetscBool use)
1056: {
1061: PetscTryMethod(pc,"PCMGSetGalerkin_C",(PC,PetscBool),(pc,use));
1062: return(0);
1063: }
1067: /*@
1068: PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e.
1069: A_i-1 = r_i * A_i * p_i
1071: Not Collective
1073: Input Parameter:
1074: . pc - the multigrid context
1076: Output Parameter:
1077: . galerkin - PETSC_TRUE or PETSC_FALSE
1079: Options Database Key:
1080: . -pc_mg_galerkin
1082: Level: intermediate
1084: .keywords: MG, set, Galerkin
1086: .seealso: PCMGSetGalerkin()
1088: @*/
1089: PetscErrorCode PCMGGetGalerkin(PC pc,PetscBool *galerkin)
1090: {
1091: PC_MG *mg = (PC_MG*)pc->data;
1095: *galerkin = (PetscBool)mg->galerkin;
1096: return(0);
1097: }
1101: /*@
1102: PCMGSetNumberSmoothDown - Sets the number of pre-smoothing steps to
1103: use on all levels. Use PCMGGetSmootherDown() to set different
1104: pre-smoothing steps on different levels.
1106: Logically Collective on PC
1108: Input Parameters:
1109: + mg - the multigrid context
1110: - n - the number of smoothing steps
1112: Options Database Key:
1113: . -pc_mg_smoothdown <n> - Sets number of pre-smoothing steps
1115: Level: advanced
1117: .keywords: MG, smooth, down, pre-smoothing, steps, multigrid
1119: .seealso: PCMGSetNumberSmoothUp()
1120: @*/
1121: PetscErrorCode PCMGSetNumberSmoothDown(PC pc,PetscInt n)
1122: {
1123: PC_MG *mg = (PC_MG*)pc->data;
1124: PC_MG_Levels **mglevels = mg->levels;
1126: PetscInt i,levels;
1130: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
1132: levels = mglevels[0]->levels;
1134: for (i=1; i<levels; i++) {
1135: /* make sure smoother up and down are different */
1136: PCMGGetSmootherUp(pc,i,NULL);
1137: KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
1139: mg->default_smoothd = n;
1140: }
1141: return(0);
1142: }
1146: /*@
1147: PCMGSetNumberSmoothUp - Sets the number of post-smoothing steps to use
1148: on all levels. Use PCMGGetSmootherUp() to set different numbers of
1149: post-smoothing steps on different levels.
1151: Logically Collective on PC
1153: Input Parameters:
1154: + mg - the multigrid context
1155: - n - the number of smoothing steps
1157: Options Database Key:
1158: . -pc_mg_smoothup <n> - Sets number of post-smoothing steps
1160: Level: advanced
1162: Note: this does not set a value on the coarsest grid, since we assume that
1163: there is no separate smooth up on the coarsest grid.
1165: .keywords: MG, smooth, up, post-smoothing, steps, multigrid
1167: .seealso: PCMGSetNumberSmoothDown()
1168: @*/
1169: PetscErrorCode PCMGSetNumberSmoothUp(PC pc,PetscInt n)
1170: {
1171: PC_MG *mg = (PC_MG*)pc->data;
1172: PC_MG_Levels **mglevels = mg->levels;
1174: PetscInt i,levels;
1178: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
1180: levels = mglevels[0]->levels;
1182: for (i=1; i<levels; i++) {
1183: /* make sure smoother up and down are different */
1184: PCMGGetSmootherUp(pc,i,NULL);
1185: KSPSetTolerances(mglevels[i]->smoothu,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
1187: mg->default_smoothu = n;
1188: }
1189: return(0);
1190: }
1192: /* ----------------------------------------------------------------------------------------*/
1194: /*MC
1195: PCMG - Use multigrid preconditioning. This preconditioner requires you provide additional
1196: information about the coarser grid matrices and restriction/interpolation operators.
1198: Options Database Keys:
1199: + -pc_mg_levels <nlevels> - number of levels including finest
1200: . -pc_mg_cycle_type <v,w> -
1201: . -pc_mg_smoothup <n> - number of smoothing steps after interpolation
1202: . -pc_mg_smoothdown <n> - number of smoothing steps before applying restriction operator
1203: . -pc_mg_type <additive,multiplicative,full,kaskade> - multiplicative is the default
1204: . -pc_mg_log - log information about time spent on each level of the solver
1205: . -pc_mg_galerkin - use Galerkin process to compute coarser operators, i.e. Acoarse = R A R'
1206: . -pc_mg_multiplicative_cycles - number of cycles to use as the preconditioner (defaults to 1)
1207: . -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices
1208: to the Socket viewer for reading from MATLAB.
1209: - -pc_mg_dump_binary - dumps the matrices for each level and the restriction/interpolation matrices
1210: to the binary output file called binaryoutput
1212: Notes: If one uses a Krylov method such GMRES or CG as the smoother than one must use KSPFGMRES, KSPGCG, or KSPRICHARDSON as the outer Krylov method
1214: When run with a single level the smoother options are used on that level NOT the coarse grid solver options
1216: Level: intermediate
1218: Concepts: multigrid/multilevel
1220: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, PCEXOTIC, PCGAMG, PCML, PCHYPRE
1221: PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), PCMGSetNumberSmoothDown(),
1222: PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(),
1223: PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(),
1224: PCMGSetCycleTypeOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR()
1225: M*/
1229: PETSC_EXTERN PetscErrorCode PCCreate_MG(PC pc)
1230: {
1231: PC_MG *mg;
1235: PetscNewLog(pc,&mg);
1236: pc->data = (void*)mg;
1237: mg->nlevels = -1;
1238: mg->am = PC_MG_MULTIPLICATIVE;
1240: pc->useAmat = PETSC_TRUE;
1242: pc->ops->apply = PCApply_MG;
1243: pc->ops->setup = PCSetUp_MG;
1244: pc->ops->reset = PCReset_MG;
1245: pc->ops->destroy = PCDestroy_MG;
1246: pc->ops->setfromoptions = PCSetFromOptions_MG;
1247: pc->ops->view = PCView_MG;
1249: PetscObjectComposeFunction((PetscObject)pc,"PCMGSetGalerkin_C",PCMGSetGalerkin_MG);
1250: return(0);
1251: }