Actual source code: mg.c
petsc-3.3-p7 2013-05-11
2: /*
3: Defines the multigrid preconditioner interface.
4: */
5: #include <../src/ksp/pc/impls/mg/mgimpl.h> /*I "petscpcmg.h" I*/
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;
19: if (mglevels->eventsmoothsolve) {PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);}
20: KSPSolve(mglevels->smoothd,mglevels->b,mglevels->x); /* pre-smooth */
21: if (mglevels->eventsmoothsolve) {PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);}
22: if (mglevels->level) { /* not the coarsest grid */
23: if (mglevels->eventresidual) {PetscLogEventBegin(mglevels->eventresidual,0,0,0,0);}
24: (*mglevels->residual)(mglevels->A,mglevels->b,mglevels->x,mglevels->r);
25: if (mglevels->eventresidual) {PetscLogEventEnd(mglevels->eventresidual,0,0,0,0);}
27: /* if on finest level and have convergence criteria set */
28: if (mglevels->level == mglevels->levels-1 && mg->ttol && reason) {
29: PetscReal rnorm;
30: VecNorm(mglevels->r,NORM_2,&rnorm);
31: if (rnorm <= mg->ttol) {
32: if (rnorm < mg->abstol) {
33: *reason = PCRICHARDSON_CONVERGED_ATOL;
34: PetscInfo2(pc,"Linear solver has converged. Residual norm %G is less than absolute tolerance %G\n",rnorm,mg->abstol);
35: } else {
36: *reason = PCRICHARDSON_CONVERGED_RTOL;
37: PetscInfo2(pc,"Linear solver has converged. Residual norm %G is less than relative tolerance times initial residual norm %G\n",rnorm,mg->ttol);
38: }
39: return(0);
40: }
41: }
43: mgc = *(mglevelsin - 1);
44: if (mglevels->eventinterprestrict) {PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);}
45: MatRestrict(mglevels->restrct,mglevels->r,mgc->b);
46: if (mglevels->eventinterprestrict) {PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);}
47: VecSet(mgc->x,0.0);
48: while (cycles--) {
49: PCMGMCycle_Private(pc,mglevelsin-1,reason);
50: }
51: if (mglevels->eventinterprestrict) {PetscLogEventBegin(mglevels->eventinterprestrict,0,0,0,0);}
52: MatInterpolateAdd(mglevels->interpolate,mgc->x,mglevels->x,mglevels->x);
53: if (mglevels->eventinterprestrict) {PetscLogEventEnd(mglevels->eventinterprestrict,0,0,0,0);}
54: if (mglevels->eventsmoothsolve) {PetscLogEventBegin(mglevels->eventsmoothsolve,0,0,0,0);}
55: KSPSolve(mglevels->smoothu,mglevels->b,mglevels->x); /* post smooth */
56: if (mglevels->eventsmoothsolve) {PetscLogEventEnd(mglevels->eventsmoothsolve,0,0,0,0);}
57: }
58: return(0);
59: }
63: 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)
64: {
65: PC_MG *mg = (PC_MG*)pc->data;
66: PC_MG_Levels **mglevels = mg->levels;
68: PetscInt levels = mglevels[0]->levels,i;
71: mglevels[levels-1]->b = b;
72: mglevels[levels-1]->x = x;
74: mg->rtol = rtol;
75: mg->abstol = abstol;
76: mg->dtol = dtol;
77: if (rtol) {
78: /* compute initial residual norm for relative convergence test */
79: PetscReal rnorm;
80: if (zeroguess) {
81: VecNorm(b,NORM_2,&rnorm);
82: } else {
83: (*mglevels[levels-1]->residual)(mglevels[levels-1]->A,b,x,w);
84: VecNorm(w,NORM_2,&rnorm);
85: }
86: mg->ttol = PetscMax(rtol*rnorm,abstol);
87: } else if (abstol) {
88: mg->ttol = abstol;
89: } else {
90: mg->ttol = 0.0;
91: }
93: /* since smoother is applied to full system, not just residual we need to make sure that smoothers don't
94: stop prematurely due to small residual */
95: for (i=1; i<levels; i++) {
96: KSPSetTolerances(mglevels[i]->smoothu,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
97: if (mglevels[i]->smoothu != mglevels[i]->smoothd) {
98: KSPSetTolerances(mglevels[i]->smoothd,0,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT);
99: }
100: }
102: *reason = (PCRichardsonConvergedReason)0;
103: for (i=0; i<its; i++) {
104: PCMGMCycle_Private(pc,mglevels+levels-1,reason);
105: if (*reason) break;
106: }
107: if (!*reason) *reason = PCRICHARDSON_CONVERGED_ITS;
108: *outits = i;
109: return(0);
110: }
114: PetscErrorCode PCReset_MG(PC pc)
115: {
116: PC_MG *mg = (PC_MG*)pc->data;
117: PC_MG_Levels **mglevels = mg->levels;
119: PetscInt i,n;
122: if (mglevels) {
123: n = mglevels[0]->levels;
124: for (i=0; i<n-1; i++) {
125: VecDestroy(&mglevels[i+1]->r);
126: VecDestroy(&mglevels[i]->b);
127: VecDestroy(&mglevels[i]->x);
128: MatDestroy(&mglevels[i+1]->restrct);
129: MatDestroy(&mglevels[i+1]->interpolate);
130: VecDestroy(&mglevels[i+1]->rscale);
131: }
133: for (i=0; i<n; i++) {
134: MatDestroy(&mglevels[i]->A);
135: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
136: KSPReset(mglevels[i]->smoothd);
137: }
138: KSPReset(mglevels[i]->smoothu);
139: }
140: }
141: return(0);
142: }
146: /*@C
147: PCMGSetLevels - Sets the number of levels to use with MG.
148: Must be called before any other MG routine.
150: Logically Collective on PC
152: Input Parameters:
153: + pc - the preconditioner context
154: . levels - the number of levels
155: - comms - optional communicators for each level; this is to allow solving the coarser problems
156: on smaller sets of processors. Use PETSC_NULL_OBJECT for default in Fortran
158: Level: intermediate
160: Notes:
161: If the number of levels is one then the multigrid uses the -mg_levels prefix
162: for setting the level options rather than the -mg_coarse prefix.
164: .keywords: MG, set, levels, multigrid
166: .seealso: PCMGSetType(), PCMGGetLevels()
167: @*/
168: PetscErrorCode PCMGSetLevels(PC pc,PetscInt levels,MPI_Comm *comms)
169: {
171: PC_MG *mg = (PC_MG*)pc->data;
172: MPI_Comm comm = ((PetscObject)pc)->comm;
173: PC_MG_Levels **mglevels = mg->levels;
174: PetscInt i;
175: PetscMPIInt size;
176: const char *prefix;
177: PC ipc;
178: PetscInt n;
183: if (mg->nlevels == levels) return(0);
184: if (mglevels) {
185: /* changing the number of levels so free up the previous stuff */
186: PCReset_MG(pc);
187: n = mglevels[0]->levels;
188: for (i=0; i<n; i++) {
189: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
190: KSPDestroy(&mglevels[i]->smoothd);
191: }
192: KSPDestroy(&mglevels[i]->smoothu);
193: PetscFree(mglevels[i]);
194: }
195: PetscFree(mg->levels);
196: }
198: mg->nlevels = levels;
200: PetscMalloc(levels*sizeof(PC_MG*),&mglevels);
201: PetscLogObjectMemory(pc,levels*(sizeof(PC_MG*)));
203: PCGetOptionsPrefix(pc,&prefix);
205: mg->stageApply = 0;
206: for (i=0; i<levels; i++) {
207: PetscNewLog(pc,PC_MG_Levels,&mglevels[i]);
208: mglevels[i]->level = i;
209: mglevels[i]->levels = levels;
210: mglevels[i]->cycles = PC_MG_CYCLE_V;
211: mg->default_smoothu = 2;
212: mg->default_smoothd = 2;
213: mglevels[i]->eventsmoothsetup = 0;
214: mglevels[i]->eventsmoothsolve = 0;
215: mglevels[i]->eventresidual = 0;
216: mglevels[i]->eventinterprestrict = 0;
218: if (comms) comm = comms[i];
219: KSPCreate(comm,&mglevels[i]->smoothd);
220: KSPSetType(mglevels[i]->smoothd,KSPCHEBYSHEV);
221: KSPSetConvergenceTest(mglevels[i]->smoothd,KSPSkipConverged,PETSC_NULL,PETSC_NULL);
222: KSPSetNormType(mglevels[i]->smoothd,KSP_NORM_NONE);
223: KSPGetPC(mglevels[i]->smoothd,&ipc);
224: PCSetType(ipc,PCSOR);
225: PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd,(PetscObject)pc,levels-i);
226: KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, i?mg->default_smoothd:1);
227: KSPSetOptionsPrefix(mglevels[i]->smoothd,prefix);
229: /* do special stuff for coarse grid */
230: if (!i && levels > 1) {
231: KSPAppendOptionsPrefix(mglevels[0]->smoothd,"mg_coarse_");
233: /* coarse solve is (redundant) LU by default */
234: KSPSetType(mglevels[0]->smoothd,KSPPREONLY);
235: KSPGetPC(mglevels[0]->smoothd,&ipc);
236: MPI_Comm_size(comm,&size);
237: if (size > 1) {
238: PCSetType(ipc,PCREDUNDANT);
239: } else {
240: PCSetType(ipc,PCLU);
241: }
243: } else {
244: char tprefix[128];
245: sprintf(tprefix,"mg_levels_%d_",(int)i);
246: KSPAppendOptionsPrefix(mglevels[i]->smoothd,tprefix);
247: }
248: PetscLogObjectParent(pc,mglevels[i]->smoothd);
249: mglevels[i]->smoothu = mglevels[i]->smoothd;
250: mg->rtol = 0.0;
251: mg->abstol = 0.0;
252: mg->dtol = 0.0;
253: mg->ttol = 0.0;
254: mg->cyclesperpcapply = 1;
255: }
256: mg->am = PC_MG_MULTIPLICATIVE;
257: mg->levels = mglevels;
258: pc->ops->applyrichardson = PCApplyRichardson_MG;
259: return(0);
260: }
265: PetscErrorCode PCDestroy_MG(PC pc)
266: {
268: PC_MG *mg = (PC_MG*)pc->data;
269: PC_MG_Levels **mglevels = mg->levels;
270: PetscInt i,n;
273: PCReset_MG(pc);
274: if (mglevels) {
275: n = mglevels[0]->levels;
276: for (i=0; i<n; i++) {
277: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
278: KSPDestroy(&mglevels[i]->smoothd);
279: }
280: KSPDestroy(&mglevels[i]->smoothu);
281: PetscFree(mglevels[i]);
282: }
283: PetscFree(mg->levels);
284: }
285: PetscFree(pc->data);
286: return(0);
287: }
291: extern PetscErrorCode PCMGACycle_Private(PC,PC_MG_Levels**);
292: extern PetscErrorCode PCMGFCycle_Private(PC,PC_MG_Levels**);
293: extern PetscErrorCode PCMGKCycle_Private(PC,PC_MG_Levels**);
295: /*
296: PCApply_MG - Runs either an additive, multiplicative, Kaskadic
297: or full cycle of multigrid.
299: Note:
300: A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle().
301: */
304: static PetscErrorCode PCApply_MG(PC pc,Vec b,Vec x)
305: {
306: PC_MG *mg = (PC_MG*)pc->data;
307: PC_MG_Levels **mglevels = mg->levels;
309: PetscInt levels = mglevels[0]->levels,i;
312: if (mg->stageApply) {PetscLogStagePush(mg->stageApply);}
313: /* When the DM is supplying the matrix then it will not exist until here */
314: for (i=0; i<levels; i++) {
315: if (!mglevels[i]->A) {
316: KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,PETSC_NULL,PETSC_NULL);
317: PetscObjectReference((PetscObject)mglevels[i]->A);
318: }
319: }
321: mglevels[levels-1]->b = b;
322: mglevels[levels-1]->x = x;
323: if (mg->am == PC_MG_MULTIPLICATIVE) {
324: VecSet(x,0.0);
325: for (i=0; i<mg->cyclesperpcapply; i++) {
326: PCMGMCycle_Private(pc,mglevels+levels-1,PETSC_NULL);
327: }
328: }
329: else if (mg->am == PC_MG_ADDITIVE) {
330: PCMGACycle_Private(pc,mglevels);
331: }
332: else if (mg->am == PC_MG_KASKADE) {
333: PCMGKCycle_Private(pc,mglevels);
334: }
335: else {
336: PCMGFCycle_Private(pc,mglevels);
337: }
338: if (mg->stageApply) {PetscLogStagePop();}
339: return(0);
340: }
345: PetscErrorCode PCSetFromOptions_MG(PC pc)
346: {
348: PetscInt m,levels = 1,cycles;
349: PetscBool flg,set;
350: PC_MG *mg = (PC_MG*)pc->data;
351: PC_MG_Levels **mglevels = mg->levels;
352: PCMGType mgtype;
353: PCMGCycleType mgctype;
356: PetscOptionsHead("Multigrid options");
357: if (!mg->levels) {
358: PetscOptionsInt("-pc_mg_levels","Number of Levels","PCMGSetLevels",levels,&levels,&flg);
359: if (!flg && pc->dm) {
360: DMGetRefineLevel(pc->dm,&levels);
361: levels++;
362: mg->usedmfornumberoflevels = PETSC_TRUE;
363: }
364: PCMGSetLevels(pc,levels,PETSC_NULL);
365: }
366: mglevels = mg->levels;
368: mgctype = (PCMGCycleType) mglevels[0]->cycles;
369: PetscOptionsEnum("-pc_mg_cycle_type","V cycle or for W-cycle","PCMGSetCycleType",PCMGCycleTypes,(PetscEnum)mgctype,(PetscEnum*)&mgctype,&flg);
370: if (flg) {
371: PCMGSetCycleType(pc,mgctype);
372: };
373: flg = PETSC_FALSE;
374: PetscOptionsBool("-pc_mg_galerkin","Use Galerkin process to compute coarser operators","PCMGSetGalerkin",flg,&flg,&set);
375: if (set) {
376: PCMGSetGalerkin(pc,flg);
377: }
378: PetscOptionsInt("-pc_mg_smoothup","Number of post-smoothing steps","PCMGSetNumberSmoothUp",1,&m,&flg);
379: if (flg) {
380: PCMGSetNumberSmoothUp(pc,m);
381: }
382: PetscOptionsInt("-pc_mg_smoothdown","Number of pre-smoothing steps","PCMGSetNumberSmoothDown",1,&m,&flg);
383: if (flg) {
384: PCMGSetNumberSmoothDown(pc,m);
385: }
386: mgtype = mg->am;
387: PetscOptionsEnum("-pc_mg_type","Multigrid type","PCMGSetType",PCMGTypes,(PetscEnum)mgtype,(PetscEnum*)&mgtype,&flg);
388: if (flg) {
389: PCMGSetType(pc,mgtype);
390: }
391: if (mg->am == PC_MG_MULTIPLICATIVE) {
392: PetscOptionsInt("-pc_mg_multiplicative_cycles","Number of cycles for each preconditioner step","PCMGSetLevels",mg->cyclesperpcapply,&cycles,&flg);
393: if (flg) {
394: PCMGMultiplicativeSetCycles(pc,cycles);
395: }
396: }
397: flg = PETSC_FALSE;
398: PetscOptionsBool("-pc_mg_log","Log times for each multigrid level","None",flg,&flg,PETSC_NULL);
399: if (flg) {
400: PetscInt i;
401: char eventname[128];
402: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
403: levels = mglevels[0]->levels;
404: for (i=0; i<levels; i++) {
405: sprintf(eventname,"MGSetup Level %d",(int)i);
406: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsetup);
407: sprintf(eventname,"MGSmooth Level %d",(int)i);
408: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsolve);
409: if (i) {
410: sprintf(eventname,"MGResid Level %d",(int)i);
411: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventresidual);
412: sprintf(eventname,"MGInterp Level %d",(int)i);
413: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventinterprestrict);
414: }
415: }
417: {
418: const char *sname = "MG Apply";
419: PetscStageLog stageLog;
420: PetscInt st;
423: PetscLogGetStageLog(&stageLog);
424: for(st = 0; st < stageLog->numStages; ++st) {
425: PetscBool same;
427: PetscStrcmp(stageLog->stageInfo[st].name, sname, &same);
428: if (same) {mg->stageApply = st;}
429: }
430: if (!mg->stageApply) {
431: PetscLogStageRegister(sname, &mg->stageApply);
432: }
433: }
434: }
435: PetscOptionsTail();
436: return(0);
437: }
439: const char *PCMGTypes[] = {"MULTIPLICATIVE","ADDITIVE","FULL","KASKADE","PCMGType","PC_MG",0};
440: const char *PCMGCycleTypes[] = {"invalid","v","w","PCMGCycleType","PC_MG_CYCLE",0};
444: PetscErrorCode PCView_MG(PC pc,PetscViewer viewer)
445: {
446: PC_MG *mg = (PC_MG*)pc->data;
447: PC_MG_Levels **mglevels = mg->levels;
449: PetscInt levels = mglevels[0]->levels,i;
450: PetscBool iascii;
453: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
454: if (iascii) {
455: PetscViewerASCIIPrintf(viewer," MG: type is %s, levels=%D cycles=%s\n", PCMGTypes[mg->am],levels,(mglevels[0]->cycles == PC_MG_CYCLE_V) ? "v" : "w");
456: if (mg->am == PC_MG_MULTIPLICATIVE) {
457: PetscViewerASCIIPrintf(viewer," Cycles per PCApply=%d\n",mg->cyclesperpcapply);
458: }
459: if (mg->galerkin) {
460: PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices\n");
461: } else {
462: PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid matrices\n");
463: }
464: for (i=0; i<levels; i++) {
465: if (!i) {
466: PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level -------------------------------\n",i);
467: } else {
468: PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);
469: }
470: PetscViewerASCIIPushTab(viewer);
471: KSPView(mglevels[i]->smoothd,viewer);
472: PetscViewerASCIIPopTab(viewer);
473: if (i && mglevels[i]->smoothd == mglevels[i]->smoothu) {
474: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");
475: } else if (i){
476: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);
477: PetscViewerASCIIPushTab(viewer);
478: KSPView(mglevels[i]->smoothu,viewer);
479: PetscViewerASCIIPopTab(viewer);
480: }
481: }
482: } else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCMG",((PetscObject)viewer)->type_name);
483: return(0);
484: }
486: #include <petsc-private/dmimpl.h>
487: #include <petsc-private/kspimpl.h>
489: /*
490: Calls setup for the KSP on each level
491: */
494: PetscErrorCode PCSetUp_MG(PC pc)
495: {
496: PC_MG *mg = (PC_MG*)pc->data;
497: PC_MG_Levels **mglevels = mg->levels;
498: PetscErrorCode ierr;
499: PetscInt i,n = mglevels[0]->levels;
500: PC cpc;
501: PetscBool preonly,lu,redundant,cholesky,svd,dump = PETSC_FALSE,opsset;
502: Mat dA,dB;
503: MatStructure uflag;
504: Vec tvec;
505: DM *dms;
506: PetscViewer viewer = 0;
509: /* FIX: Move this to PCSetFromOptions_MG? */
510: if (mg->usedmfornumberoflevels) {
511: PetscInt levels;
512: DMGetRefineLevel(pc->dm,&levels);
513: levels++;
514: if (levels > n) { /* the problem is now being solved on a finer grid */
515: PCMGSetLevels(pc,levels,PETSC_NULL);
516: n = levels;
517: PCSetFromOptions(pc); /* it is bad to call this here, but otherwise will never be called for the new hierarchy */
518: mglevels = mg->levels;
519: }
520: }
521: KSPGetPC(mglevels[0]->smoothd,&cpc);
524: /* If user did not provide fine grid operators OR operator was not updated since last global KSPSetOperators() */
525: /* so use those from global PC */
526: /* Is this what we always want? What if user wants to keep old one? */
527: KSPGetOperatorsSet(mglevels[n-1]->smoothd,PETSC_NULL,&opsset);
528: if (opsset) {
529: Mat mmat;
530: KSPGetOperators(mglevels[n-1]->smoothd,PETSC_NULL,&mmat,PETSC_NULL);
531: if (mmat == pc->pmat) opsset = PETSC_FALSE;
532: }
533: if (!opsset) {
534: PetscInfo(pc,"Using outer operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
535: KSPSetOperators(mglevels[n-1]->smoothd,pc->mat,pc->pmat,pc->flag);
536: }
538: /* Skipping this for galerkin==2 (externally managed hierarchy such as ML and GAMG). Cleaner logic here would be great. Wrap ML/GAMG as DMs? */
539: if (pc->dm && mg->galerkin != 2 && !pc->setupcalled) {
540: /* construct the interpolation from the DMs */
541: Mat p;
542: Vec rscale;
543: PetscMalloc(n*sizeof(DM),&dms);
544: dms[n-1] = pc->dm;
545: for (i=n-2; i>-1; i--) {
546: KSPDM kdm;
547: DMCoarsen(dms[i+1],MPI_COMM_NULL,&dms[i]);
548: KSPSetDM(mglevels[i]->smoothd,dms[i]);
549: if (mg->galerkin) {KSPSetDMActive(mglevels[i]->smoothd,PETSC_FALSE);}
550: DMKSPGetContextWrite(dms[i],&kdm);
551: /* Ugly hack so that the next KSPSetUp() will use the RHS that we set */
552: kdm->computerhs = PETSC_NULL;
553: kdm->rhsctx = PETSC_NULL;
554: DMSetFunction(dms[i],0);
555: DMSetInitialGuess(dms[i],0);
556: if (!mglevels[i+1]->interpolate) {
557: DMCreateInterpolation(dms[i],dms[i+1],&p,&rscale);
558: PCMGSetInterpolation(pc,i+1,p);
559: if (rscale) {PCMGSetRScale(pc,i+1,rscale);}
560: VecDestroy(&rscale);
561: MatDestroy(&p);
562: }
563: }
565: for (i=n-2; i>-1; i--) {
566: DMDestroy(&dms[i]);
567: }
568: PetscFree(dms);
569: }
571: if (pc->dm && !pc->setupcalled) {
572: /* finest smoother also gets DM but it is not active, independent of whether galerkin==2 */
573: KSPSetDM(mglevels[n-1]->smoothd,pc->dm);
574: KSPSetDMActive(mglevels[n-1]->smoothd,PETSC_FALSE);
575: }
577: if (mg->galerkin == 1) {
578: Mat B;
579: /* currently only handle case where mat and pmat are the same on coarser levels */
580: KSPGetOperators(mglevels[n-1]->smoothd,&dA,&dB,&uflag);
581: if (!pc->setupcalled) {
582: for (i=n-2; i>-1; i--) {
583: MatPtAP(dB,mglevels[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);
584: KSPSetOperators(mglevels[i]->smoothd,B,B,uflag);
585: if (i != n-2) {PetscObjectDereference((PetscObject)dB);}
586: dB = B;
587: }
588: if (n > 1) {PetscObjectDereference((PetscObject)dB);}
589: } else {
590: for (i=n-2; i>-1; i--) {
591: KSPGetOperators(mglevels[i]->smoothd,PETSC_NULL,&B,PETSC_NULL);
592: MatPtAP(dB,mglevels[i+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);
593: KSPSetOperators(mglevels[i]->smoothd,B,B,uflag);
594: dB = B;
595: }
596: }
597: } else if (!mg->galerkin && pc->dm && pc->dm->x) {
598: /* need to restrict Jacobian location to coarser meshes for evaluation */
599: for (i=n-2;i>-1; i--) {
600: Mat R;
601: Vec rscale;
602: if (!mglevels[i]->smoothd->dm->x) {
603: Vec *vecs;
604: KSPGetVecs(mglevels[i]->smoothd,1,&vecs,0,PETSC_NULL);
605: mglevels[i]->smoothd->dm->x = vecs[0];
606: PetscFree(vecs);
607: }
608: PCMGGetRestriction(pc,i+1,&R);
609: PCMGGetRScale(pc,i+1,&rscale);
610: MatRestrict(R,mglevels[i+1]->smoothd->dm->x,mglevels[i]->smoothd->dm->x);
611: VecPointwiseMult(mglevels[i]->smoothd->dm->x,mglevels[i]->smoothd->dm->x,rscale);
612: }
613: }
614: if (!mg->galerkin && pc->dm) {
615: for (i=n-2;i>=0; i--) {
616: DM dmfine,dmcoarse;
617: Mat Restrict,Inject;
618: Vec rscale;
619: KSPGetDM(mglevels[i+1]->smoothd,&dmfine);
620: KSPGetDM(mglevels[i]->smoothd,&dmcoarse);
621: PCMGGetRestriction(pc,i+1,&Restrict);
622: PCMGGetRScale(pc,i+1,&rscale);
623: Inject = PETSC_NULL; /* Callback should create it if it needs Injection */
624: DMRestrict(dmfine,Restrict,rscale,Inject,dmcoarse);
625: }
626: }
628: if (!pc->setupcalled) {
629: for (i=0; i<n; i++) {
630: KSPSetFromOptions(mglevels[i]->smoothd);
631: }
632: for (i=1; i<n; i++) {
633: if (mglevels[i]->smoothu && (mglevels[i]->smoothu != mglevels[i]->smoothd)) {
634: KSPSetFromOptions(mglevels[i]->smoothu);
635: }
636: }
637: for (i=1; i<n; i++) {
638: PCMGGetInterpolation(pc,i,&mglevels[i]->interpolate);
639: PCMGGetRestriction(pc,i,&mglevels[i]->restrct);
640: }
641: for (i=0; i<n-1; i++) {
642: if (!mglevels[i]->b) {
643: Vec *vec;
644: KSPGetVecs(mglevels[i]->smoothd,1,&vec,0,PETSC_NULL);
645: PCMGSetRhs(pc,i,*vec);
646: VecDestroy(vec);
647: PetscFree(vec);
648: }
649: if (!mglevels[i]->r && i) {
650: VecDuplicate(mglevels[i]->b,&tvec);
651: PCMGSetR(pc,i,tvec);
652: VecDestroy(&tvec);
653: }
654: if (!mglevels[i]->x) {
655: VecDuplicate(mglevels[i]->b,&tvec);
656: PCMGSetX(pc,i,tvec);
657: VecDestroy(&tvec);
658: }
659: }
660: if (n != 1 && !mglevels[n-1]->r) {
661: /* PCMGSetR() on the finest level if user did not supply it */
662: Vec *vec;
663: KSPGetVecs(mglevels[n-1]->smoothd,1,&vec,0,PETSC_NULL);
664: PCMGSetR(pc,n-1,*vec);
665: VecDestroy(vec);
666: PetscFree(vec);
667: }
668: }
670: if (pc->dm) {
671: /* need to tell all the coarser levels to rebuild the matrix using the DM for that level */
672: for (i=0; i<n-1; i++){
673: if (mglevels[i]->smoothd->setupstage != KSP_SETUP_NEW) mglevels[i]->smoothd->setupstage = KSP_SETUP_NEWMATRIX;
674: }
675: }
677: for (i=1; i<n; i++) {
678: if (mglevels[i]->smoothu == mglevels[i]->smoothd) {
679: /* if doing only down then initial guess is zero */
680: KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);
681: }
682: if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
683: KSPSetUp(mglevels[i]->smoothd);
684: if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
685: if (!mglevels[i]->residual) {
686: Mat mat;
687: KSPGetOperators(mglevels[i]->smoothd,PETSC_NULL,&mat,PETSC_NULL);
688: PCMGSetResidual(pc,i,PCMGDefaultResidual,mat);
689: }
690: }
691: for (i=1; i<n; i++) {
692: if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) {
693: Mat downmat,downpmat;
694: MatStructure matflag;
696: /* check if operators have been set for up, if not use down operators to set them */
697: KSPGetOperatorsSet(mglevels[i]->smoothu,&opsset,PETSC_NULL);
698: if (!opsset) {
699: KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat,&matflag);
700: KSPSetOperators(mglevels[i]->smoothu,downmat,downpmat,matflag);
701: }
703: KSPSetInitialGuessNonzero(mglevels[i]->smoothu,PETSC_TRUE);
704: if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
705: KSPSetUp(mglevels[i]->smoothu);
706: if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
707: }
708: }
710: /*
711: If coarse solver is not direct method then DO NOT USE preonly
712: */
713: PetscObjectTypeCompare((PetscObject)mglevels[0]->smoothd,KSPPREONLY,&preonly);
714: if (preonly) {
715: PetscObjectTypeCompare((PetscObject)cpc,PCLU,&lu);
716: PetscObjectTypeCompare((PetscObject)cpc,PCREDUNDANT,&redundant);
717: PetscObjectTypeCompare((PetscObject)cpc,PCCHOLESKY,&cholesky);
718: PetscObjectTypeCompare((PetscObject)cpc,PCSVD,&svd);
719: if (!lu && !redundant && !cholesky && !svd) {
720: KSPSetType(mglevels[0]->smoothd,KSPGMRES);
721: }
722: }
724: if (!pc->setupcalled) {
725: KSPSetFromOptions(mglevels[0]->smoothd);
726: }
728: if (mglevels[0]->eventsmoothsetup) {PetscLogEventBegin(mglevels[0]->eventsmoothsetup,0,0,0,0);}
729: KSPSetUp(mglevels[0]->smoothd);
730: if (mglevels[0]->eventsmoothsetup) {PetscLogEventEnd(mglevels[0]->eventsmoothsetup,0,0,0,0);}
732: /*
733: Dump the interpolation/restriction matrices plus the
734: Jacobian/stiffness on each level. This allows MATLAB users to
735: easily check if the Galerkin condition A_c = R A_f R^T is satisfied.
737: Only support one or the other at the same time.
738: */
739: #if defined(PETSC_USE_SOCKET_VIEWER)
740: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_mg_dump_matlab",&dump,PETSC_NULL);
741: if (dump) {
742: viewer = PETSC_VIEWER_SOCKET_(((PetscObject)pc)->comm);
743: }
744: dump = PETSC_FALSE;
745: #endif
746: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_mg_dump_binary",&dump,PETSC_NULL);
747: if (dump) {
748: viewer = PETSC_VIEWER_BINARY_(((PetscObject)pc)->comm);
749: }
751: if (viewer) {
752: for (i=1; i<n; i++) {
753: MatView(mglevels[i]->restrct,viewer);
754: }
755: for (i=0; i<n; i++) {
756: KSPGetPC(mglevels[i]->smoothd,&pc);
757: MatView(pc->mat,viewer);
758: }
759: }
760: return(0);
761: }
763: /* -------------------------------------------------------------------------------------*/
767: /*@
768: PCMGGetLevels - Gets the number of levels to use with MG.
770: Not Collective
772: Input Parameter:
773: . pc - the preconditioner context
775: Output parameter:
776: . levels - the number of levels
778: Level: advanced
780: .keywords: MG, get, levels, multigrid
782: .seealso: PCMGSetLevels()
783: @*/
784: PetscErrorCode PCMGGetLevels(PC pc,PetscInt *levels)
785: {
786: PC_MG *mg = (PC_MG*)pc->data;
791: *levels = mg->nlevels;
792: return(0);
793: }
797: /*@
798: PCMGSetType - Determines the form of multigrid to use:
799: multiplicative, additive, full, or the Kaskade algorithm.
801: Logically Collective on PC
803: Input Parameters:
804: + pc - the preconditioner context
805: - form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,
806: PC_MG_FULL, PC_MG_KASKADE
808: Options Database Key:
809: . -pc_mg_type <form> - Sets <form>, one of multiplicative,
810: additive, full, kaskade
812: Level: advanced
814: .keywords: MG, set, method, multiplicative, additive, full, Kaskade, multigrid
816: .seealso: PCMGSetLevels()
817: @*/
818: PetscErrorCode PCMGSetType(PC pc,PCMGType form)
819: {
820: PC_MG *mg = (PC_MG*)pc->data;
825: mg->am = form;
826: if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG;
827: else pc->ops->applyrichardson = 0;
828: return(0);
829: }
833: /*@
834: PCMGSetCycleType - Sets the type cycles to use. Use PCMGSetCycleTypeOnLevel() for more
835: complicated cycling.
837: Logically Collective on PC
839: Input Parameters:
840: + pc - the multigrid context
841: - PC_MG_CYCLE_V or PC_MG_CYCLE_W
843: Options Database Key:
844: $ -pc_mg_cycle_type v or w
846: Level: advanced
848: .keywords: MG, set, cycles, V-cycle, W-cycle, multigrid
850: .seealso: PCMGSetCycleTypeOnLevel()
851: @*/
852: PetscErrorCode PCMGSetCycleType(PC pc,PCMGCycleType n)
853: {
854: PC_MG *mg = (PC_MG*)pc->data;
855: PC_MG_Levels **mglevels = mg->levels;
856: PetscInt i,levels;
860: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
862: levels = mglevels[0]->levels;
864: for (i=0; i<levels; i++) {
865: mglevels[i]->cycles = n;
866: }
867: return(0);
868: }
872: /*@
873: PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step
874: of multigrid when PCMGType of PC_MG_MULTIPLICATIVE is used
876: Logically Collective on PC
878: Input Parameters:
879: + pc - the multigrid context
880: - n - number of cycles (default is 1)
882: Options Database Key:
883: $ -pc_mg_multiplicative_cycles n
885: Level: advanced
887: Notes: This is not associated with setting a v or w cycle, that is set with PCMGSetCycleType()
889: .keywords: MG, set, cycles, V-cycle, W-cycle, multigrid
891: .seealso: PCMGSetCycleTypeOnLevel(), PCMGSetCycleType()
892: @*/
893: PetscErrorCode PCMGMultiplicativeSetCycles(PC pc,PetscInt n)
894: {
895: PC_MG *mg = (PC_MG*)pc->data;
896: PC_MG_Levels **mglevels = mg->levels;
897: PetscInt i,levels;
901: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
903: levels = mglevels[0]->levels;
905: for (i=0; i<levels; i++) {
906: mg->cyclesperpcapply = n;
907: }
908: return(0);
909: }
913: /*@
914: PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the
915: finest grid via the Galerkin process: A_i-1 = r_i * A_i * r_i^t
917: Logically Collective on PC
919: Input Parameters:
920: + pc - the multigrid context
921: - use - PETSC_TRUE to use the Galerkin process to compute coarse-level operators
923: Options Database Key:
924: $ -pc_mg_galerkin
926: Level: intermediate
928: .keywords: MG, set, Galerkin
930: .seealso: PCMGGetGalerkin()
932: @*/
933: PetscErrorCode PCMGSetGalerkin(PC pc,PetscBool use)
934: {
935: PC_MG *mg = (PC_MG*)pc->data;
939: mg->galerkin = use ? 1 : 0;
940: return(0);
941: }
945: /*@
946: PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e.
947: A_i-1 = r_i * A_i * r_i^t
949: Not Collective
951: Input Parameter:
952: . pc - the multigrid context
954: Output Parameter:
955: . gelerkin - PETSC_TRUE or PETSC_FALSE
957: Options Database Key:
958: $ -pc_mg_galerkin
960: Level: intermediate
962: .keywords: MG, set, Galerkin
964: .seealso: PCMGSetGalerkin()
966: @*/
967: PetscErrorCode PCMGGetGalerkin(PC pc,PetscBool *galerkin)
968: {
969: PC_MG *mg = (PC_MG*)pc->data;
973: *galerkin = (PetscBool)mg->galerkin;
974: return(0);
975: }
979: /*@
980: PCMGSetNumberSmoothDown - Sets the number of pre-smoothing steps to
981: use on all levels. Use PCMGGetSmootherDown() to set different
982: pre-smoothing steps on different levels.
984: Logically Collective on PC
986: Input Parameters:
987: + mg - the multigrid context
988: - n - the number of smoothing steps
990: Options Database Key:
991: . -pc_mg_smoothdown <n> - Sets number of pre-smoothing steps
993: Level: advanced
995: .keywords: MG, smooth, down, pre-smoothing, steps, multigrid
997: .seealso: PCMGSetNumberSmoothUp()
998: @*/
999: PetscErrorCode PCMGSetNumberSmoothDown(PC pc,PetscInt n)
1000: {
1001: PC_MG *mg = (PC_MG*)pc->data;
1002: PC_MG_Levels **mglevels = mg->levels;
1004: PetscInt i,levels;
1008: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
1010: levels = mglevels[0]->levels;
1012: for (i=1; i<levels; i++) {
1013: /* make sure smoother up and down are different */
1014: PCMGGetSmootherUp(pc,i,PETSC_NULL);
1015: KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
1016: mg->default_smoothd = n;
1017: }
1018: return(0);
1019: }
1023: /*@
1024: PCMGSetNumberSmoothUp - Sets the number of post-smoothing steps to use
1025: on all levels. Use PCMGGetSmootherUp() to set different numbers of
1026: post-smoothing steps on different levels.
1028: Logically Collective on PC
1030: Input Parameters:
1031: + mg - the multigrid context
1032: - n - the number of smoothing steps
1034: Options Database Key:
1035: . -pc_mg_smoothup <n> - Sets number of post-smoothing steps
1037: Level: advanced
1039: Note: this does not set a value on the coarsest grid, since we assume that
1040: there is no separate smooth up on the coarsest grid.
1042: .keywords: MG, smooth, up, post-smoothing, steps, multigrid
1044: .seealso: PCMGSetNumberSmoothDown()
1045: @*/
1046: PetscErrorCode PCMGSetNumberSmoothUp(PC pc,PetscInt n)
1047: {
1048: PC_MG *mg = (PC_MG*)pc->data;
1049: PC_MG_Levels **mglevels = mg->levels;
1051: PetscInt i,levels;
1055: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
1057: levels = mglevels[0]->levels;
1059: for (i=1; i<levels; i++) {
1060: /* make sure smoother up and down are different */
1061: PCMGGetSmootherUp(pc,i,PETSC_NULL);
1062: KSPSetTolerances(mglevels[i]->smoothu,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
1063: mg->default_smoothu = n;
1064: }
1065: return(0);
1066: }
1068: /* ----------------------------------------------------------------------------------------*/
1070: /*MC
1071: PCMG - Use multigrid preconditioning. This preconditioner requires you provide additional
1072: information about the coarser grid matrices and restriction/interpolation operators.
1074: Options Database Keys:
1075: + -pc_mg_levels <nlevels> - number of levels including finest
1076: . -pc_mg_cycles v or w
1077: . -pc_mg_smoothup <n> - number of smoothing steps after interpolation
1078: . -pc_mg_smoothdown <n> - number of smoothing steps before applying restriction operator
1079: . -pc_mg_type <additive,multiplicative,full,cascade> - multiplicative is the default
1080: . -pc_mg_log - log information about time spent on each level of the solver
1081: . -pc_mg_monitor - print information on the multigrid convergence
1082: . -pc_mg_galerkin - use Galerkin process to compute coarser operators, i.e. Acoarse = R A R'
1083: . -pc_mg_multiplicative_cycles - number of cycles to use as the preconditioner (defaults to 1)
1084: . -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices
1085: to the Socket viewer for reading from MATLAB.
1086: - -pc_mg_dump_binary - dumps the matrices for each level and the restriction/interpolation matrices
1087: to the binary output file called binaryoutput
1089: 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
1091: When run with a single level the smoother options are used on that level NOT the coarse grid solver options
1093: Level: intermediate
1095: Concepts: multigrid/multilevel
1097: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, PCEXOTIC, PCGAMG, PCML, PCHYPRE
1098: PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), PCMGSetNumberSmoothDown(),
1099: PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(),
1100: PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(),
1101: PCMGSetCycleTypeOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR()
1102: M*/
1104: EXTERN_C_BEGIN
1107: PetscErrorCode PCCreate_MG(PC pc)
1108: {
1109: PC_MG *mg;
1113: PetscNewLog(pc,PC_MG,&mg);
1114: pc->data = (void*)mg;
1115: mg->nlevels = -1;
1117: pc->ops->apply = PCApply_MG;
1118: pc->ops->setup = PCSetUp_MG;
1119: pc->ops->reset = PCReset_MG;
1120: pc->ops->destroy = PCDestroy_MG;
1121: pc->ops->setfromoptions = PCSetFromOptions_MG;
1122: pc->ops->view = PCView_MG;
1123: return(0);
1124: }
1125: EXTERN_C_END