Actual source code: mg.c
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 do 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: for (i=0; i<levels; i++) {
206: PetscNewLog(pc,PC_MG_Levels,&mglevels[i]);
207: mglevels[i]->level = i;
208: mglevels[i]->levels = levels;
209: mglevels[i]->cycles = PC_MG_CYCLE_V;
210: mg->default_smoothu = 1;
211: mg->default_smoothd = 1;
212: mglevels[i]->eventsmoothsetup = 0;
213: mglevels[i]->eventsmoothsolve = 0;
214: mglevels[i]->eventresidual = 0;
215: mglevels[i]->eventinterprestrict = 0;
217: if (comms) comm = comms[i];
218: KSPCreate(comm,&mglevels[i]->smoothd);
219: PetscObjectIncrementTabLevel((PetscObject)mglevels[i]->smoothd,(PetscObject)pc,levels-i);
220: KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT, mg->default_smoothd);
221: KSPSetOptionsPrefix(mglevels[i]->smoothd,prefix);
223: /* do special stuff for coarse grid */
224: if (!i && levels > 1) {
225: KSPAppendOptionsPrefix(mglevels[0]->smoothd,"mg_coarse_");
227: /* coarse solve is (redundant) LU by default */
228: KSPSetType(mglevels[0]->smoothd,KSPPREONLY);
229: KSPGetPC(mglevels[0]->smoothd,&ipc);
230: MPI_Comm_size(comm,&size);
231: if (size > 1) {
232: PCSetType(ipc,PCREDUNDANT);
233: } else {
234: PCSetType(ipc,PCLU);
235: }
237: } else {
238: char tprefix[128];
239: sprintf(tprefix,"mg_levels_%d_",(int)i);
240: KSPAppendOptionsPrefix(mglevels[i]->smoothd,tprefix);
241: }
242: PetscLogObjectParent(pc,mglevels[i]->smoothd);
243: mglevels[i]->smoothu = mglevels[i]->smoothd;
244: mg->rtol = 0.0;
245: mg->abstol = 0.0;
246: mg->dtol = 0.0;
247: mg->ttol = 0.0;
248: mg->cyclesperpcapply = 1;
249: }
250: mg->am = PC_MG_MULTIPLICATIVE;
251: mg->levels = mglevels;
252: pc->ops->applyrichardson = PCApplyRichardson_MG;
253: return(0);
254: }
259: PetscErrorCode PCDestroy_MG(PC pc)
260: {
262: PC_MG *mg = (PC_MG*)pc->data;
263: PC_MG_Levels **mglevels = mg->levels;
264: PetscInt i,n;
267: PCReset_MG(pc);
268: if (mglevels) {
269: n = mglevels[0]->levels;
270: for (i=0; i<n; i++) {
271: if (mglevels[i]->smoothd != mglevels[i]->smoothu) {
272: KSPDestroy(&mglevels[i]->smoothd);
273: }
274: KSPDestroy(&mglevels[i]->smoothu);
275: PetscFree(mglevels[i]);
276: }
277: PetscFree(mg->levels);
278: }
279: PetscFree(pc->data);
280: return(0);
281: }
285: extern PetscErrorCode PCMGACycle_Private(PC,PC_MG_Levels**);
286: extern PetscErrorCode PCMGFCycle_Private(PC,PC_MG_Levels**);
287: extern PetscErrorCode PCMGKCycle_Private(PC,PC_MG_Levels**);
289: /*
290: PCApply_MG - Runs either an additive, multiplicative, Kaskadic
291: or full cycle of multigrid.
293: Note:
294: A simple wrapper which calls PCMGMCycle(),PCMGACycle(), or PCMGFCycle().
295: */
298: static PetscErrorCode PCApply_MG(PC pc,Vec b,Vec x)
299: {
300: PC_MG *mg = (PC_MG*)pc->data;
301: PC_MG_Levels **mglevels = mg->levels;
303: PetscInt levels = mglevels[0]->levels,i;
307: /* When the DM is supplying the matrix then it will not exist until here */
308: for (i=0; i<levels; i++) {
309: if (!mglevels[i]->A) {
310: KSPGetOperators(mglevels[i]->smoothu,&mglevels[i]->A,PETSC_NULL,PETSC_NULL);
311: PetscObjectReference((PetscObject)mglevels[i]->A);
312: }
313: }
315: mglevels[levels-1]->b = b;
316: mglevels[levels-1]->x = x;
317: if (mg->am == PC_MG_MULTIPLICATIVE) {
318: VecSet(x,0.0);
319: for (i=0; i<mg->cyclesperpcapply; i++) {
320: PCMGMCycle_Private(pc,mglevels+levels-1,PETSC_NULL);
321: }
322: }
323: else if (mg->am == PC_MG_ADDITIVE) {
324: PCMGACycle_Private(pc,mglevels);
325: }
326: else if (mg->am == PC_MG_KASKADE) {
327: PCMGKCycle_Private(pc,mglevels);
328: }
329: else {
330: PCMGFCycle_Private(pc,mglevels);
331: }
332: return(0);
333: }
338: PetscErrorCode PCSetFromOptions_MG(PC pc)
339: {
341: PetscInt m,levels = 1,cycles;
342: PetscBool flg,set;
343: PC_MG *mg = (PC_MG*)pc->data;
344: PC_MG_Levels **mglevels = mg->levels;
345: PCMGType mgtype;
346: PCMGCycleType mgctype;
349: PetscOptionsHead("Multigrid options");
350: if (!mg->levels) {
351: PetscOptionsInt("-pc_mg_levels","Number of Levels","PCMGSetLevels",levels,&levels,&flg);
352: if (!flg && pc->dm) {
353: DMGetRefineLevel(pc->dm,&levels);
354: levels++;
355: mg->usedmfornumberoflevels = PETSC_TRUE;
356: }
357: PCMGSetLevels(pc,levels,PETSC_NULL);
358: }
359: mglevels = mg->levels;
361: mgctype = (PCMGCycleType) mglevels[0]->cycles;
362: PetscOptionsEnum("-pc_mg_cycle_type","V cycle or for W-cycle","PCMGSetCycleType",PCMGCycleTypes,(PetscEnum)mgctype,(PetscEnum*)&mgctype,&flg);
363: if (flg) {
364: PCMGSetCycleType(pc,mgctype);
365: };
366: flg = PETSC_FALSE;
367: PetscOptionsBool("-pc_mg_galerkin","Use Galerkin process to compute coarser operators","PCMGSetGalerkin",flg,&flg,&set);
368: if (set) {
369: PCMGSetGalerkin(pc,flg);
370: }
371: PetscOptionsInt("-pc_mg_smoothup","Number of post-smoothing steps","PCMGSetNumberSmoothUp",1,&m,&flg);
372: if (flg) {
373: PCMGSetNumberSmoothUp(pc,m);
374: }
375: PetscOptionsInt("-pc_mg_smoothdown","Number of pre-smoothing steps","PCMGSetNumberSmoothDown",1,&m,&flg);
376: if (flg) {
377: PCMGSetNumberSmoothDown(pc,m);
378: }
379: mgtype = mg->am;
380: PetscOptionsEnum("-pc_mg_type","Multigrid type","PCMGSetType",PCMGTypes,(PetscEnum)mgtype,(PetscEnum*)&mgtype,&flg);
381: if (flg) {
382: PCMGSetType(pc,mgtype);
383: }
384: if (mg->am == PC_MG_MULTIPLICATIVE) {
385: PetscOptionsInt("-pc_mg_multiplicative_cycles","Number of cycles for each preconditioner step","PCMGSetLevels",mg->cyclesperpcapply,&cycles,&flg);
386: if (flg) {
387: PCMGMultiplicativeSetCycles(pc,cycles);
388: }
389: }
390: flg = PETSC_FALSE;
391: PetscOptionsBool("-pc_mg_log","Log times for each multigrid level","None",flg,&flg,PETSC_NULL);
392: if (flg) {
393: PetscInt i;
394: char eventname[128];
395: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
396: levels = mglevels[0]->levels;
397: for (i=0; i<levels; i++) {
398: sprintf(eventname,"MGSetup Level %d",(int)i);
399: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsetup);
400: sprintf(eventname,"MGSmooth Level %d",(int)i);
401: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventsmoothsolve);
402: if (i) {
403: sprintf(eventname,"MGResid Level %d",(int)i);
404: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventresidual);
405: sprintf(eventname,"MGInterp Level %d",(int)i);
406: PetscLogEventRegister(eventname,((PetscObject)pc)->classid,&mglevels[i]->eventinterprestrict);
407: }
408: }
409: }
410: PetscOptionsTail();
411: return(0);
412: }
414: const char *PCMGTypes[] = {"MULTIPLICATIVE","ADDITIVE","FULL","KASKADE","PCMGType","PC_MG",0};
415: const char *PCMGCycleTypes[] = {"invalid","v","w","PCMGCycleType","PC_MG_CYCLE",0};
419: PetscErrorCode PCView_MG(PC pc,PetscViewer viewer)
420: {
421: PC_MG *mg = (PC_MG*)pc->data;
422: PC_MG_Levels **mglevels = mg->levels;
424: PetscInt levels = mglevels[0]->levels,i;
425: PetscBool iascii;
428: PetscTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
429: if (iascii) {
430: PetscViewerASCIIPrintf(viewer," MG: type is %s, levels=%D cycles=%s\n", PCMGTypes[mg->am],levels,(mglevels[0]->cycles == PC_MG_CYCLE_V) ? "v" : "w");
431: if (mg->am == PC_MG_MULTIPLICATIVE) {
432: PetscViewerASCIIPrintf(viewer," Cycles per PCApply=%d\n",mg->cyclesperpcapply);
433: }
434: if (mg->galerkin) {
435: PetscViewerASCIIPrintf(viewer," Using Galerkin computed coarse grid matrices\n");
436: } else {
437: PetscViewerASCIIPrintf(viewer," Not using Galerkin computed coarse grid matrices\n");
438: }
439: for (i=0; i<levels; i++) {
440: if (!i) {
441: PetscViewerASCIIPrintf(viewer,"Coarse grid solver -- level -------------------------------\n",i);
442: } else {
443: PetscViewerASCIIPrintf(viewer,"Down solver (pre-smoother) on level %D -------------------------------\n",i);
444: }
445: PetscViewerASCIIPushTab(viewer);
446: KSPView(mglevels[i]->smoothd,viewer);
447: PetscViewerASCIIPopTab(viewer);
448: if (i && mglevels[i]->smoothd == mglevels[i]->smoothu) {
449: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) same as down solver (pre-smoother)\n");
450: } else if (i){
451: PetscViewerASCIIPrintf(viewer,"Up solver (post-smoother) on level %D -------------------------------\n",i);
452: PetscViewerASCIIPushTab(viewer);
453: KSPView(mglevels[i]->smoothu,viewer);
454: PetscViewerASCIIPopTab(viewer);
455: }
456: }
457: } else SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_SUP,"Viewer type %s not supported for PCMG",((PetscObject)viewer)->type_name);
458: return(0);
459: }
461: #include <private/dmimpl.h>
462: #include <private/kspimpl.h>
464: /*
465: Calls setup for the KSP on each level
466: */
469: PetscErrorCode PCSetUp_MG(PC pc)
470: {
471: PC_MG *mg = (PC_MG*)pc->data;
472: PC_MG_Levels **mglevels = mg->levels;
473: PetscErrorCode ierr;
474: PetscInt i,n = mglevels[0]->levels;
475: PC cpc,mpc;
476: PetscBool preonly,lu,redundant,cholesky,svd,dump = PETSC_FALSE,opsset;
477: Mat dA,dB;
478: MatStructure uflag;
479: Vec tvec;
480: DM *dms;
481: PetscViewer viewer = 0;
484: if (mg->usedmfornumberoflevels) {
485: PetscInt levels;
486: DMGetRefineLevel(pc->dm,&levels);
487: levels++;
488: if (levels > n) { /* the problem is now being solved on a finer grid */
489: PCMGSetLevels(pc,levels,PETSC_NULL);
490: n = levels;
491: PCSetFromOptions(pc); /* it is bad to call this here, but otherwise will never be called for the new hierarchy */
492: mglevels = mg->levels;
493: }
494: }
497: /* If user did not provide fine grid operators OR operator was not updated since last global KSPSetOperators() */
498: /* so use those from global PC */
499: /* Is this what we always want? What if user wants to keep old one? */
500: KSPGetOperatorsSet(mglevels[n-1]->smoothd,PETSC_NULL,&opsset);
501: KSPGetPC(mglevels[0]->smoothd,&cpc);
502: KSPGetPC(mglevels[n-1]->smoothd,&mpc);
503: if (!opsset || ((cpc->setupcalled == 1) && (mpc->setupcalled == 2)) || ((mpc == cpc) && (mpc->setupcalled == 2))) {
504: PetscInfo(pc,"Using outer operators to define finest grid operator \n because PCMGGetSmoother(pc,nlevels-1,&ksp);KSPSetOperators(ksp,...); was not called.\n");
505: KSPSetOperators(mglevels[n-1]->smoothd,pc->mat,pc->pmat,pc->flag);
506: }
508: /* Skipping this for galerkin==2 (externally managed hierarchy such as ML and GAMG). Cleaner logic here would be great. */
509: if (pc->dm && mg->galerkin != 2 && !pc->setupcalled) {
510: /* construct the interpolation from the DMs */
511: Mat p;
512: Vec rscale;
513: PetscMalloc(n*sizeof(DM),&dms);
514: dms[n-1] = pc->dm;
515: for (i=n-2; i>-1; i--) {
516: DMCoarsen(dms[i+1],PETSC_NULL,&dms[i]);
517: KSPSetDM(mglevels[i]->smoothd,dms[i]);
518: if (mg->galerkin) {KSPSetDMActive(mglevels[i]->smoothd,PETSC_FALSE);}
519: DMSetFunction(dms[i],0);
520: DMSetInitialGuess(dms[i],0);
521: if (!mglevels[i+1]->interpolate) {
522: DMCreateInterpolation(dms[i],dms[i+1],&p,&rscale);
523: PCMGSetInterpolation(pc,i+1,p);
524: if (rscale) {PCMGSetRScale(pc,i+1,rscale);}
525: VecDestroy(&rscale);
526: MatDestroy(&p);
527: }
528: }
530: for (i=n-2; i>-1; i--) {
531: DMDestroy(&dms[i]);
532: }
533: PetscFree(dms);
534: }
536: if (pc->dm && !pc->setupcalled) {
537: /* finest smoother also gets DM but it is not active, independent of whether galerkin==2 */
538: KSPSetDM(mglevels[n-1]->smoothd,pc->dm);
539: KSPSetDMActive(mglevels[n-1]->smoothd,PETSC_FALSE);
540: }
542: if (mg->galerkin == 1) {
543: Mat B;
544: /* currently only handle case where mat and pmat are the same on coarser levels */
545: KSPGetOperators(mglevels[n-1]->smoothd,&dA,&dB,&uflag);
546: if (!pc->setupcalled) {
547: for (i=n-2; i>-1; i--) {
548: MatPtAP(dB,mglevels[i+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);
549: KSPSetOperators(mglevels[i]->smoothd,B,B,uflag);
550: if (i != n-2) {PetscObjectDereference((PetscObject)dB);}
551: dB = B;
552: }
553: if (n > 1) {PetscObjectDereference((PetscObject)dB);}
554: } else {
555: for (i=n-2; i>-1; i--) {
556: KSPGetOperators(mglevels[i]->smoothd,PETSC_NULL,&B,PETSC_NULL);
557: MatPtAP(dB,mglevels[i+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);
558: KSPSetOperators(mglevels[i]->smoothd,B,B,uflag);
559: dB = B;
560: }
561: }
562: } else if (!mg->galerkin && pc->dm && pc->dm->x) {
563: /* need to restrict Jacobian location to coarser meshes for evaluation */
564: for (i=n-2;i>-1; i--) {
565: Mat R;
566: Vec rscale;
567: if (!mglevels[i]->smoothd->dm->x) {
568: Vec *vecs;
569: KSPGetVecs(mglevels[i]->smoothd,1,&vecs,0,PETSC_NULL);
570: mglevels[i]->smoothd->dm->x = vecs[0];
571: PetscFree(vecs);
572: }
573: PCMGGetRestriction(pc,i+1,&R);
574: PCMGGetRScale(pc,i+1,&rscale);
575: MatRestrict(R,mglevels[i+1]->smoothd->dm->x,mglevels[i]->smoothd->dm->x);
576: VecPointwiseMult(mglevels[i]->smoothd->dm->x,mglevels[i]->smoothd->dm->x,rscale);
577: }
578: }
580: if (!pc->setupcalled) {
581: for (i=0; i<n; i++) {
582: KSPSetFromOptions(mglevels[i]->smoothd);
583: }
584: for (i=1; i<n; i++) {
585: if (mglevels[i]->smoothu && (mglevels[i]->smoothu != mglevels[i]->smoothd)) {
586: KSPSetFromOptions(mglevels[i]->smoothu);
587: }
588: }
589: for (i=1; i<n; i++) {
590: PCMGGetInterpolation(pc,i,&mglevels[i]->interpolate);
591: PCMGGetRestriction(pc,i,&mglevels[i]->restrct);
592: }
593: for (i=0; i<n-1; i++) {
594: if (!mglevels[i]->b) {
595: Vec *vec;
596: KSPGetVecs(mglevels[i]->smoothd,1,&vec,0,PETSC_NULL);
597: PCMGSetRhs(pc,i,*vec);
598: VecDestroy(vec);
599: PetscFree(vec);
600: }
601: if (!mglevels[i]->r && i) {
602: VecDuplicate(mglevels[i]->b,&tvec);
603: PCMGSetR(pc,i,tvec);
604: VecDestroy(&tvec);
605: }
606: if (!mglevels[i]->x) {
607: VecDuplicate(mglevels[i]->b,&tvec);
608: PCMGSetX(pc,i,tvec);
609: VecDestroy(&tvec);
610: }
611: }
612: if (n != 1 && !mglevels[n-1]->r) {
613: /* PCMGSetR() on the finest level if user did not supply it */
614: Vec *vec;
615: KSPGetVecs(mglevels[n-1]->smoothd,1,&vec,0,PETSC_NULL);
616: PCMGSetR(pc,n-1,*vec);
617: VecDestroy(vec);
618: PetscFree(vec);
619: }
620: }
623: for (i=1; i<n; i++) {
624: if (mglevels[i]->smoothu == mglevels[i]->smoothd) {
625: /* if doing only down then initial guess is zero */
626: KSPSetInitialGuessNonzero(mglevels[i]->smoothd,PETSC_TRUE);
627: }
628: if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
629: KSPSetUp(mglevels[i]->smoothd);
630: if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
631: if (!mglevels[i]->residual) {
632: Mat mat;
633: KSPGetOperators(mglevels[i]->smoothd,PETSC_NULL,&mat,PETSC_NULL);
634: PCMGSetResidual(pc,i,PCMGDefaultResidual,mat);
635: }
636: }
637: for (i=1; i<n; i++) {
638: if (mglevels[i]->smoothu && mglevels[i]->smoothu != mglevels[i]->smoothd) {
639: Mat downmat,downpmat;
640: MatStructure matflag;
641: PetscBool opsset;
643: /* check if operators have been set for up, if not use down operators to set them */
644: KSPGetOperatorsSet(mglevels[i]->smoothu,&opsset,PETSC_NULL);
645: if (!opsset) {
646: KSPGetOperators(mglevels[i]->smoothd,&downmat,&downpmat,&matflag);
647: KSPSetOperators(mglevels[i]->smoothu,downmat,downpmat,matflag);
648: }
650: KSPSetInitialGuessNonzero(mglevels[i]->smoothu,PETSC_TRUE);
651: if (mglevels[i]->eventsmoothsetup) {PetscLogEventBegin(mglevels[i]->eventsmoothsetup,0,0,0,0);}
652: KSPSetUp(mglevels[i]->smoothu);
653: if (mglevels[i]->eventsmoothsetup) {PetscLogEventEnd(mglevels[i]->eventsmoothsetup,0,0,0,0);}
654: }
655: }
657: /*
658: If coarse solver is not direct method then DO NOT USE preonly
659: */
660: PetscTypeCompare((PetscObject)mglevels[0]->smoothd,KSPPREONLY,&preonly);
661: if (preonly) {
662: PetscTypeCompare((PetscObject)cpc,PCLU,&lu);
663: PetscTypeCompare((PetscObject)cpc,PCREDUNDANT,&redundant);
664: PetscTypeCompare((PetscObject)cpc,PCCHOLESKY,&cholesky);
665: PetscTypeCompare((PetscObject)cpc,PCSVD,&svd);
666: if (!lu && !redundant && !cholesky && !svd) {
667: KSPSetType(mglevels[0]->smoothd,KSPGMRES);
668: }
669: }
671: if (!pc->setupcalled) {
672: KSPSetFromOptions(mglevels[0]->smoothd);
673: }
675: if (mglevels[0]->eventsmoothsetup) {PetscLogEventBegin(mglevels[0]->eventsmoothsetup,0,0,0,0);}
676: KSPSetUp(mglevels[0]->smoothd);
677: if (mglevels[0]->eventsmoothsetup) {PetscLogEventEnd(mglevels[0]->eventsmoothsetup,0,0,0,0);}
679: /*
680: Dump the interpolation/restriction matrices plus the
681: Jacobian/stiffness on each level. This allows MATLAB users to
682: easily check if the Galerkin condition A_c = R A_f R^T is satisfied.
684: Only support one or the other at the same time.
685: */
686: #if defined(PETSC_USE_SOCKET_VIEWER)
687: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_mg_dump_matlab",&dump,PETSC_NULL);
688: if (dump) {
689: viewer = PETSC_VIEWER_SOCKET_(((PetscObject)pc)->comm);
690: }
691: dump = PETSC_FALSE;
692: #endif
693: PetscOptionsGetBool(((PetscObject)pc)->prefix,"-pc_mg_dump_binary",&dump,PETSC_NULL);
694: if (dump) {
695: viewer = PETSC_VIEWER_BINARY_(((PetscObject)pc)->comm);
696: }
698: if (viewer) {
699: for (i=1; i<n; i++) {
700: MatView(mglevels[i]->restrct,viewer);
701: }
702: for (i=0; i<n; i++) {
703: KSPGetPC(mglevels[i]->smoothd,&pc);
704: MatView(pc->mat,viewer);
705: }
706: }
707: return(0);
708: }
710: /* -------------------------------------------------------------------------------------*/
714: /*@
715: PCMGGetLevels - Gets the number of levels to use with MG.
717: Not Collective
719: Input Parameter:
720: . pc - the preconditioner context
722: Output parameter:
723: . levels - the number of levels
725: Level: advanced
727: .keywords: MG, get, levels, multigrid
729: .seealso: PCMGSetLevels()
730: @*/
731: PetscErrorCode PCMGGetLevels(PC pc,PetscInt *levels)
732: {
733: PC_MG *mg = (PC_MG*)pc->data;
738: *levels = mg->nlevels;
739: return(0);
740: }
744: /*@
745: PCMGSetType - Determines the form of multigrid to use:
746: multiplicative, additive, full, or the Kaskade algorithm.
748: Logically Collective on PC
750: Input Parameters:
751: + pc - the preconditioner context
752: - form - multigrid form, one of PC_MG_MULTIPLICATIVE, PC_MG_ADDITIVE,
753: PC_MG_FULL, PC_MG_KASKADE
755: Options Database Key:
756: . -pc_mg_type <form> - Sets <form>, one of multiplicative,
757: additive, full, kaskade
759: Level: advanced
761: .keywords: MG, set, method, multiplicative, additive, full, Kaskade, multigrid
763: .seealso: PCMGSetLevels()
764: @*/
765: PetscErrorCode PCMGSetType(PC pc,PCMGType form)
766: {
767: PC_MG *mg = (PC_MG*)pc->data;
772: mg->am = form;
773: if (form == PC_MG_MULTIPLICATIVE) pc->ops->applyrichardson = PCApplyRichardson_MG;
774: else pc->ops->applyrichardson = 0;
775: return(0);
776: }
780: /*@
781: PCMGSetCycleType - Sets the type cycles to use. Use PCMGSetCycleTypeOnLevel() for more
782: complicated cycling.
784: Logically Collective on PC
786: Input Parameters:
787: + pc - the multigrid context
788: - PC_MG_CYCLE_V or PC_MG_CYCLE_W
790: Options Database Key:
791: $ -pc_mg_cycle_type v or w
793: Level: advanced
795: .keywords: MG, set, cycles, V-cycle, W-cycle, multigrid
797: .seealso: PCMGSetCycleTypeOnLevel()
798: @*/
799: PetscErrorCode PCMGSetCycleType(PC pc,PCMGCycleType n)
800: {
801: PC_MG *mg = (PC_MG*)pc->data;
802: PC_MG_Levels **mglevels = mg->levels;
803: PetscInt i,levels;
807: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
809: levels = mglevels[0]->levels;
811: for (i=0; i<levels; i++) {
812: mglevels[i]->cycles = n;
813: }
814: return(0);
815: }
819: /*@
820: PCMGMultiplicativeSetCycles - Sets the number of cycles to use for each preconditioner step
821: of multigrid when PCMGType of PC_MG_MULTIPLICATIVE is used
823: Logically Collective on PC
825: Input Parameters:
826: + pc - the multigrid context
827: - n - number of cycles (default is 1)
829: Options Database Key:
830: $ -pc_mg_multiplicative_cycles n
832: Level: advanced
834: Notes: This is not associated with setting a v or w cycle, that is set with PCMGSetCycleType()
836: .keywords: MG, set, cycles, V-cycle, W-cycle, multigrid
838: .seealso: PCMGSetCycleTypeOnLevel(), PCMGSetCycleType()
839: @*/
840: PetscErrorCode PCMGMultiplicativeSetCycles(PC pc,PetscInt n)
841: {
842: PC_MG *mg = (PC_MG*)pc->data;
843: PC_MG_Levels **mglevels = mg->levels;
844: PetscInt i,levels;
848: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
850: levels = mglevels[0]->levels;
852: for (i=0; i<levels; i++) {
853: mg->cyclesperpcapply = n;
854: }
855: return(0);
856: }
860: /*@
861: PCMGSetGalerkin - Causes the coarser grid matrices to be computed from the
862: finest grid via the Galerkin process: A_i-1 = r_i * A_i * r_i^t
864: Logically Collective on PC
866: Input Parameters:
867: + pc - the multigrid context
868: - use - PETSC_TRUE to use the Galerkin process to compute coarse-level operators
870: Options Database Key:
871: $ -pc_mg_galerkin
873: Level: intermediate
875: .keywords: MG, set, Galerkin
877: .seealso: PCMGGetGalerkin()
879: @*/
880: PetscErrorCode PCMGSetGalerkin(PC pc,PetscBool use)
881: {
882: PC_MG *mg = (PC_MG*)pc->data;
886: mg->galerkin = (PetscInt)use;
887: return(0);
888: }
892: /*@
893: PCMGGetGalerkin - Checks if Galerkin multigrid is being used, i.e.
894: A_i-1 = r_i * A_i * r_i^t
896: Not Collective
898: Input Parameter:
899: . pc - the multigrid context
901: Output Parameter:
902: . gelerkin - PETSC_TRUE or PETSC_FALSE
904: Options Database Key:
905: $ -pc_mg_galerkin
907: Level: intermediate
909: .keywords: MG, set, Galerkin
911: .seealso: PCMGSetGalerkin()
913: @*/
914: PetscErrorCode PCMGGetGalerkin(PC pc,PetscBool *galerkin)
915: {
916: PC_MG *mg = (PC_MG*)pc->data;
920: *galerkin = (PetscBool)mg->galerkin;
921: return(0);
922: }
926: /*@
927: PCMGSetNumberSmoothDown - Sets the number of pre-smoothing steps to
928: use on all levels. Use PCMGGetSmootherDown() to set different
929: pre-smoothing steps on different levels.
931: Logically Collective on PC
933: Input Parameters:
934: + mg - the multigrid context
935: - n - the number of smoothing steps
937: Options Database Key:
938: . -pc_mg_smoothdown <n> - Sets number of pre-smoothing steps
940: Level: advanced
942: .keywords: MG, smooth, down, pre-smoothing, steps, multigrid
944: .seealso: PCMGSetNumberSmoothUp()
945: @*/
946: PetscErrorCode PCMGSetNumberSmoothDown(PC pc,PetscInt n)
947: {
948: PC_MG *mg = (PC_MG*)pc->data;
949: PC_MG_Levels **mglevels = mg->levels;
951: PetscInt i,levels;
955: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
957: levels = mglevels[0]->levels;
959: for (i=1; i<levels; i++) {
960: /* make sure smoother up and down are different */
961: PCMGGetSmootherUp(pc,i,PETSC_NULL);
962: KSPSetTolerances(mglevels[i]->smoothd,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
963: mg->default_smoothd = n;
964: }
965: return(0);
966: }
970: /*@
971: PCMGSetNumberSmoothUp - Sets the number of post-smoothing steps to use
972: on all levels. Use PCMGGetSmootherUp() to set different numbers of
973: post-smoothing steps on different levels.
975: Logically Collective on PC
977: Input Parameters:
978: + mg - the multigrid context
979: - n - the number of smoothing steps
981: Options Database Key:
982: . -pc_mg_smoothup <n> - Sets number of post-smoothing steps
984: Level: advanced
986: Note: this does not set a value on the coarsest grid, since we assume that
987: there is no separate smooth up on the coarsest grid.
989: .keywords: MG, smooth, up, post-smoothing, steps, multigrid
991: .seealso: PCMGSetNumberSmoothDown()
992: @*/
993: PetscErrorCode PCMGSetNumberSmoothUp(PC pc,PetscInt n)
994: {
995: PC_MG *mg = (PC_MG*)pc->data;
996: PC_MG_Levels **mglevels = mg->levels;
998: PetscInt i,levels;
1002: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
1004: levels = mglevels[0]->levels;
1006: for (i=1; i<levels; i++) {
1007: /* make sure smoother up and down are different */
1008: PCMGGetSmootherUp(pc,i,PETSC_NULL);
1009: KSPSetTolerances(mglevels[i]->smoothu,PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,n);
1010: mg->default_smoothu = n;
1011: }
1012: return(0);
1013: }
1015: /* ----------------------------------------------------------------------------------------*/
1017: /*MC
1018: PCMG - Use multigrid preconditioning. This preconditioner requires you provide additional
1019: information about the coarser grid matrices and restriction/interpolation operators.
1021: Options Database Keys:
1022: + -pc_mg_levels <nlevels> - number of levels including finest
1023: . -pc_mg_cycles v or w
1024: . -pc_mg_smoothup <n> - number of smoothing steps after interpolation
1025: . -pc_mg_smoothdown <n> - number of smoothing steps before applying restriction operator
1026: . -pc_mg_type <additive,multiplicative,full,cascade> - multiplicative is the default
1027: . -pc_mg_log - log information about time spent on each level of the solver
1028: . -pc_mg_monitor - print information on the multigrid convergence
1029: . -pc_mg_galerkin - use Galerkin process to compute coarser operators
1030: - -pc_mg_dump_matlab - dumps the matrices for each level and the restriction/interpolation matrices
1031: to the Socket viewer for reading from MATLAB.
1033: 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
1035: Level: intermediate
1037: Concepts: multigrid/multilevel
1039: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, PCEXOTIC
1040: PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), PCMGSetCycleType(), PCMGSetNumberSmoothDown(),
1041: PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(),
1042: PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(),
1043: PCMGSetCycleTypeOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR()
1044: M*/
1046: EXTERN_C_BEGIN
1049: PetscErrorCode PCCreate_MG(PC pc)
1050: {
1051: PC_MG *mg;
1055: PetscNewLog(pc,PC_MG,&mg);
1056: pc->data = (void*)mg;
1057: mg->nlevels = -1;
1059: pc->ops->apply = PCApply_MG;
1060: pc->ops->setup = PCSetUp_MG;
1061: pc->ops->reset = PCReset_MG;
1062: pc->ops->destroy = PCDestroy_MG;
1063: pc->ops->setfromoptions = PCSetFromOptions_MG;
1064: pc->ops->view = PCView_MG;
1065: return(0);
1066: }
1067: EXTERN_C_END