Actual source code: mgfunc.c
petsc-3.4.5 2014-06-29
2: #include <../src/ksp/pc/impls/mg/mgimpl.h> /*I "petscksp.h" I*/
6: /*@C
7: PCMGResidual_Default - Default routine to calculate the residual.
9: Collective on Mat and Vec
11: Input Parameters:
12: + mat - the matrix
13: . b - the right-hand-side
14: - x - the approximate solution
16: Output Parameter:
17: . r - location to store the residual
19: Level: developer
21: .keywords: MG, default, multigrid, residual
23: .seealso: PCMGSetResidual()
24: @*/
25: PetscErrorCode PCMGResidual_Default(Mat mat,Vec b,Vec x,Vec r)
26: {
30: MatMult(mat,x,r);
31: VecAYPX(r,-1.0,b);
32: return(0);
33: }
35: /* ---------------------------------------------------------------------------*/
39: /*@
40: PCMGGetCoarseSolve - Gets the solver context to be used on the coarse grid.
42: Not Collective
44: Input Parameter:
45: . pc - the multigrid context
47: Output Parameter:
48: . ksp - the coarse grid solver context
50: Level: advanced
52: .keywords: MG, multigrid, get, coarse grid
53: @*/
54: PetscErrorCode PCMGGetCoarseSolve(PC pc,KSP *ksp)
55: {
56: PC_MG *mg = (PC_MG*)pc->data;
57: PC_MG_Levels **mglevels = mg->levels;
61: *ksp = mglevels[0]->smoothd;
62: return(0);
63: }
67: /*@C
68: PCMGSetResidual - Sets the function to be used to calculate the residual
69: on the lth level.
71: Logically Collective on PC and Mat
73: Input Parameters:
74: + pc - the multigrid context
75: . l - the level (0 is coarsest) to supply
76: . residual - function used to form residual, if none is provided the previously provide one is used, if no
77: previous one were provided then a default is used
78: - mat - matrix associated with residual
80: Level: advanced
82: .keywords: MG, set, multigrid, residual, level
84: .seealso: PCMGResidual_Default()
85: @*/
86: PetscErrorCode PCMGSetResidual(PC pc,PetscInt l,PetscErrorCode (*residual)(Mat,Vec,Vec,Vec),Mat mat)
87: {
88: PC_MG *mg = (PC_MG*)pc->data;
89: PC_MG_Levels **mglevels = mg->levels;
94: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
95: if (residual) mglevels[l]->residual = residual;
96: if (!mglevels[l]->residual) mglevels[l]->residual = PCMGResidual_Default;
97: if (mat) {PetscObjectReference((PetscObject)mat);}
98: MatDestroy(&mglevels[l]->A);
100: mglevels[l]->A = mat;
101: return(0);
102: }
106: /*@
107: PCMGSetInterpolation - Sets the function to be used to calculate the
108: interpolation from l-1 to the lth level
110: Logically Collective on PC and Mat
112: Input Parameters:
113: + pc - the multigrid context
114: . mat - the interpolation operator
115: - l - the level (0 is coarsest) to supply [do not supply 0]
117: Level: advanced
119: Notes:
120: Usually this is the same matrix used also to set the restriction
121: for the same level.
123: One can pass in the interpolation matrix or its transpose; PETSc figures
124: out from the matrix size which one it is.
126: .keywords: multigrid, set, interpolate, level
128: .seealso: PCMGSetRestriction()
129: @*/
130: PetscErrorCode PCMGSetInterpolation(PC pc,PetscInt l,Mat mat)
131: {
132: PC_MG *mg = (PC_MG*)pc->data;
133: PC_MG_Levels **mglevels = mg->levels;
138: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
139: if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Do not set interpolation routine for coarsest level");
140: PetscObjectReference((PetscObject)mat);
141: MatDestroy(&mglevels[l]->interpolate);
143: mglevels[l]->interpolate = mat;
144: return(0);
145: }
149: /*@
150: PCMGGetInterpolation - Gets the function to be used to calculate the
151: interpolation from l-1 to the lth level
153: Logically Collective on PC
155: Input Parameters:
156: + pc - the multigrid context
157: - l - the level (0 is coarsest) to supply [Do not supply 0]
159: Output Parameter:
160: . mat - the interpolation matrix
162: Level: advanced
164: .keywords: MG, get, multigrid, interpolation, level
166: .seealso: PCMGGetRestriction(), PCMGSetInterpolation(), PCMGGetRScale()
167: @*/
168: PetscErrorCode PCMGGetInterpolation(PC pc,PetscInt l,Mat *mat)
169: {
170: PC_MG *mg = (PC_MG*)pc->data;
171: PC_MG_Levels **mglevels = mg->levels;
177: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
178: if (l <= 0 || mg->nlevels <= l) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Level %D must be in range {1,...,%D}",l,mg->nlevels-1);
179: if (!mglevels[l]->interpolate) {
180: if (!mglevels[l]->restrct) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call PCMGSetInterpolation() or PCMGSetInterpolation()");
181: PCMGSetInterpolation(pc,l,mglevels[l]->restrct);
182: }
183: *mat = mglevels[l]->interpolate;
184: return(0);
185: }
189: /*@
190: PCMGSetRestriction - Sets the function to be used to restrict vector
191: from level l to l-1.
193: Logically Collective on PC and Mat
195: Input Parameters:
196: + pc - the multigrid context
197: . l - the level (0 is coarsest) to supply [Do not supply 0]
198: - mat - the restriction matrix
200: Level: advanced
202: Notes:
203: Usually this is the same matrix used also to set the interpolation
204: for the same level.
206: One can pass in the interpolation matrix or its transpose; PETSc figures
207: out from the matrix size which one it is.
209: If you do not set this, the transpose of the Mat set with PCMGSetInterpolation()
210: is used.
212: .keywords: MG, set, multigrid, restriction, level
214: .seealso: PCMGSetInterpolation()
215: @*/
216: PetscErrorCode PCMGSetRestriction(PC pc,PetscInt l,Mat mat)
217: {
219: PC_MG *mg = (PC_MG*)pc->data;
220: PC_MG_Levels **mglevels = mg->levels;
225: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
226: if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Do not set restriction routine for coarsest level");
227: PetscObjectReference((PetscObject)mat);
228: MatDestroy(&mglevels[l]->restrct);
230: mglevels[l]->restrct = mat;
231: return(0);
232: }
236: /*@
237: PCMGGetRestriction - Gets the function to be used to restrict vector
238: from level l to l-1.
240: Logically Collective on PC and Mat
242: Input Parameters:
243: + pc - the multigrid context
244: - l - the level (0 is coarsest) to supply [Do not supply 0]
246: Output Parameter:
247: . mat - the restriction matrix
249: Level: advanced
251: .keywords: MG, get, multigrid, restriction, level
253: .seealso: PCMGGetInterpolation(), PCMGSetRestriction(), PCMGGetRScale()
254: @*/
255: PetscErrorCode PCMGGetRestriction(PC pc,PetscInt l,Mat *mat)
256: {
257: PC_MG *mg = (PC_MG*)pc->data;
258: PC_MG_Levels **mglevels = mg->levels;
264: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
265: if (l <= 0 || mg->nlevels <= l) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Level %D must be in range {1,...,%D}",l,mg->nlevels-1);
266: if (!mglevels[l]->restrct) {
267: if (!mglevels[l]->interpolate) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call PCMGSetRestriction() or PCMGSetInterpolation()");
268: PCMGSetRestriction(pc,l,mglevels[l]->interpolate);
269: }
270: *mat = mglevels[l]->restrct;
271: return(0);
272: }
276: /*@
277: PCMGSetRScale - Sets the pointwise scaling for the restriction operator from level l to l-1.
279: Logically Collective on PC and Vec
281: Input Parameters:
282: + pc - the multigrid context
283: - l - the level (0 is coarsest) to supply [Do not supply 0]
284: . rscale - the scaling
286: Level: advanced
288: Notes:
289: When evaluating a function on a coarse level one does not want to do F(R * x) one does F(rscale * R * x) where rscale is 1 over the row sums of R.
291: .keywords: MG, set, multigrid, restriction, level
293: .seealso: PCMGSetInterpolation(), PCMGSetRestriction(), PCMGGetRScale()
294: @*/
295: PetscErrorCode PCMGSetRScale(PC pc,PetscInt l,Vec rscale)
296: {
298: PC_MG *mg = (PC_MG*)pc->data;
299: PC_MG_Levels **mglevels = mg->levels;
303: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
304: if (l <= 0 || mg->nlevels <= l) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Level %D must be in range {1,...,%D}",l,mg->nlevels-1);
305: PetscObjectReference((PetscObject)rscale);
306: VecDestroy(&mglevels[l]->rscale);
308: mglevels[l]->rscale = rscale;
309: return(0);
310: }
314: /*@
315: PCMGGetRScale - Gets the pointwise scaling for the restriction operator from level l to l-1.
317: Collective on PC
319: Input Parameters:
320: + pc - the multigrid context
321: . rscale - the scaling
322: - l - the level (0 is coarsest) to supply [Do not supply 0]
324: Level: advanced
326: Notes:
327: When evaluating a function on a coarse level one does not want to do F(R * x) one does F(rscale * R * x) where rscale is 1 over the row sums of R.
329: .keywords: MG, set, multigrid, restriction, level
331: .seealso: PCMGSetInterpolation(), PCMGGetRestriction()
332: @*/
333: PetscErrorCode PCMGGetRScale(PC pc,PetscInt l,Vec *rscale)
334: {
336: PC_MG *mg = (PC_MG*)pc->data;
337: PC_MG_Levels **mglevels = mg->levels;
341: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
342: if (l <= 0 || mg->nlevels <= l) SETERRQ2(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Level %D must be in range {1,...,%D}",l,mg->nlevels-1);
343: if (!mglevels[l]->rscale) {
344: Mat R;
345: Vec X,Y,coarse,fine;
346: PetscInt M,N;
347: PCMGGetRestriction(pc,l,&R);
348: MatGetVecs(R,&X,&Y);
349: MatGetSize(R,&M,&N);
350: if (M < N) {
351: fine = X;
352: coarse = Y;
353: } else if (N < M) {
354: fine = Y; coarse = X;
355: } else SETERRQ(PetscObjectComm((PetscObject)R),PETSC_ERR_SUP,"Restriction matrix is square, cannot determine which Vec is coarser");
356: VecSet(fine,1.);
357: MatRestrict(R,fine,coarse);
358: VecDestroy(&fine);
359: VecReciprocal(coarse);
360: mglevels[l]->rscale = coarse;
361: }
362: *rscale = mglevels[l]->rscale;
363: return(0);
364: }
368: /*@
369: PCMGGetSmoother - Gets the KSP context to be used as smoother for
370: both pre- and post-smoothing. Call both PCMGGetSmootherUp() and
371: PCMGGetSmootherDown() to use different functions for pre- and
372: post-smoothing.
374: Not Collective, KSP returned is parallel if PC is
376: Input Parameters:
377: + pc - the multigrid context
378: - l - the level (0 is coarsest) to supply
380: Ouput Parameters:
381: . ksp - the smoother
383: Level: advanced
385: .keywords: MG, get, multigrid, level, smoother, pre-smoother, post-smoother
387: .seealso: PCMGGetSmootherUp(), PCMGGetSmootherDown()
388: @*/
389: PetscErrorCode PCMGGetSmoother(PC pc,PetscInt l,KSP *ksp)
390: {
391: PC_MG *mg = (PC_MG*)pc->data;
392: PC_MG_Levels **mglevels = mg->levels;
396: *ksp = mglevels[l]->smoothd;
397: return(0);
398: }
402: /*@
403: PCMGGetSmootherUp - Gets the KSP context to be used as smoother after
404: coarse grid correction (post-smoother).
406: Not Collective, KSP returned is parallel if PC is
408: Input Parameters:
409: + pc - the multigrid context
410: - l - the level (0 is coarsest) to supply
412: Ouput Parameters:
413: . ksp - the smoother
415: Level: advanced
417: Notes: calling this will result in a different pre and post smoother so you may need to
418: set options on the pre smoother also
420: .keywords: MG, multigrid, get, smoother, up, post-smoother, level
422: .seealso: PCMGGetSmootherUp(), PCMGGetSmootherDown()
423: @*/
424: PetscErrorCode PCMGGetSmootherUp(PC pc,PetscInt l,KSP *ksp)
425: {
426: PC_MG *mg = (PC_MG*)pc->data;
427: PC_MG_Levels **mglevels = mg->levels;
429: const char *prefix;
430: MPI_Comm comm;
434: /*
435: This is called only if user wants a different pre-smoother from post.
436: Thus we check if a different one has already been allocated,
437: if not we allocate it.
438: */
439: if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"There is no such thing as a up smoother on the coarse grid");
440: if (mglevels[l]->smoothu == mglevels[l]->smoothd) {
441: KSPType ksptype;
442: PCType pctype;
443: PC ipc;
444: PetscReal rtol,abstol,dtol;
445: PetscInt maxits;
446: KSPNormType normtype;
447: PetscObjectGetComm((PetscObject)mglevels[l]->smoothd,&comm);
448: KSPGetOptionsPrefix(mglevels[l]->smoothd,&prefix);
449: KSPGetTolerances(mglevels[l]->smoothd,&rtol,&abstol,&dtol,&maxits);
450: KSPGetType(mglevels[l]->smoothd,&ksptype);
451: KSPGetNormType(mglevels[l]->smoothd,&normtype);
452: KSPGetPC(mglevels[l]->smoothd,&ipc);
453: PCGetType(ipc,&pctype);
455: KSPCreate(comm,&mglevels[l]->smoothu);
456: PetscObjectIncrementTabLevel((PetscObject)mglevels[l]->smoothu,(PetscObject)pc,mglevels[0]->levels-l);
457: KSPSetOptionsPrefix(mglevels[l]->smoothu,prefix);
458: KSPSetTolerances(mglevels[l]->smoothu,rtol,abstol,dtol,maxits);
459: KSPSetType(mglevels[l]->smoothu,ksptype);
460: KSPSetNormType(mglevels[l]->smoothu,normtype);
461: KSPSetConvergenceTest(mglevels[l]->smoothu,KSPSkipConverged,NULL,NULL);
462: KSPGetPC(mglevels[l]->smoothu,&ipc);
463: PCSetType(ipc,pctype);
464: PetscLogObjectParent(pc,mglevels[l]->smoothu);
465: }
466: if (ksp) *ksp = mglevels[l]->smoothu;
467: return(0);
468: }
472: /*@
473: PCMGGetSmootherDown - Gets the KSP context to be used as smoother before
474: coarse grid correction (pre-smoother).
476: Not Collective, KSP returned is parallel if PC is
478: Input Parameters:
479: + pc - the multigrid context
480: - l - the level (0 is coarsest) to supply
482: Ouput Parameters:
483: . ksp - the smoother
485: Level: advanced
487: Notes: calling this will result in a different pre and post smoother so you may need to
488: set options on the post smoother also
490: .keywords: MG, multigrid, get, smoother, down, pre-smoother, level
492: .seealso: PCMGGetSmootherUp(), PCMGGetSmoother()
493: @*/
494: PetscErrorCode PCMGGetSmootherDown(PC pc,PetscInt l,KSP *ksp)
495: {
497: PC_MG *mg = (PC_MG*)pc->data;
498: PC_MG_Levels **mglevels = mg->levels;
502: /* make sure smoother up and down are different */
503: if (l) {
504: PCMGGetSmootherUp(pc,l,NULL);
505: }
506: *ksp = mglevels[l]->smoothd;
507: return(0);
508: }
512: /*@
513: PCMGSetCyclesOnLevel - Sets the number of cycles to run on this level.
515: Logically Collective on PC
517: Input Parameters:
518: + pc - the multigrid context
519: . l - the level (0 is coarsest) this is to be used for
520: - n - the number of cycles
522: Level: advanced
524: .keywords: MG, multigrid, set, cycles, V-cycle, W-cycle, level
526: .seealso: PCMGSetCycles()
527: @*/
528: PetscErrorCode PCMGSetCyclesOnLevel(PC pc,PetscInt l,PetscInt c)
529: {
530: PC_MG *mg = (PC_MG*)pc->data;
531: PC_MG_Levels **mglevels = mg->levels;
535: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
538: mglevels[l]->cycles = c;
539: return(0);
540: }
544: /*@
545: PCMGSetRhs - Sets the vector space to be used to store the right-hand side
546: on a particular level.
548: Logically Collective on PC and Vec
550: Input Parameters:
551: + pc - the multigrid context
552: . l - the level (0 is coarsest) this is to be used for
553: - c - the space
555: Level: advanced
557: Notes: If this is not provided PETSc will automatically generate one.
559: You do not need to keep a reference to this vector if you do
560: not need it PCDestroy() will properly free it.
562: .keywords: MG, multigrid, set, right-hand-side, rhs, level
564: .seealso: PCMGSetX(), PCMGSetR()
565: @*/
566: PetscErrorCode PCMGSetRhs(PC pc,PetscInt l,Vec c)
567: {
569: PC_MG *mg = (PC_MG*)pc->data;
570: PC_MG_Levels **mglevels = mg->levels;
574: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
575: if (l == mglevels[0]->levels-1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Do not set rhs for finest level");
576: PetscObjectReference((PetscObject)c);
577: VecDestroy(&mglevels[l]->b);
579: mglevels[l]->b = c;
580: return(0);
581: }
585: /*@
586: PCMGSetX - Sets the vector space to be used to store the solution on a
587: particular level.
589: Logically Collective on PC and Vec
591: Input Parameters:
592: + pc - the multigrid context
593: . l - the level (0 is coarsest) this is to be used for (do not supply the finest level)
594: - c - the space
596: Level: advanced
598: Notes: If this is not provided PETSc will automatically generate one.
600: You do not need to keep a reference to this vector if you do
601: not need it PCDestroy() will properly free it.
603: .keywords: MG, multigrid, set, solution, level
605: .seealso: PCMGSetRhs(), PCMGSetR()
606: @*/
607: PetscErrorCode PCMGSetX(PC pc,PetscInt l,Vec c)
608: {
610: PC_MG *mg = (PC_MG*)pc->data;
611: PC_MG_Levels **mglevels = mg->levels;
615: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
616: if (l == mglevels[0]->levels-1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Do not set x for finest level");
617: PetscObjectReference((PetscObject)c);
618: VecDestroy(&mglevels[l]->x);
620: mglevels[l]->x = c;
621: return(0);
622: }
626: /*@
627: PCMGSetR - Sets the vector space to be used to store the residual on a
628: particular level.
630: Logically Collective on PC and Vec
632: Input Parameters:
633: + pc - the multigrid context
634: . l - the level (0 is coarsest) this is to be used for
635: - c - the space
637: Level: advanced
639: Notes: If this is not provided PETSc will automatically generate one.
641: You do not need to keep a reference to this vector if you do
642: not need it PCDestroy() will properly free it.
644: .keywords: MG, multigrid, set, residual, level
645: @*/
646: PetscErrorCode PCMGSetR(PC pc,PetscInt l,Vec c)
647: {
649: PC_MG *mg = (PC_MG*)pc->data;
650: PC_MG_Levels **mglevels = mg->levels;
654: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
655: if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Need not set residual vector for coarse grid");
656: PetscObjectReference((PetscObject)c);
657: VecDestroy(&mglevels[l]->r);
659: mglevels[l]->r = c;
660: return(0);
661: }