Actual source code: mgfunc.c
petsc-3.3-p7 2013-05-11
2: #include <../src/ksp/pc/impls/mg/mgimpl.h> /*I "petscksp.h" I*/
3: /*I "petscpcmg.h" I*/
7: /*@C
8: PCMGDefaultResidual - Default routine to calculate the residual.
10: Collective on Mat and Vec
12: Input Parameters:
13: + mat - the matrix
14: . b - the right-hand-side
15: - x - the approximate solution
16:
17: Output Parameter:
18: . r - location to store the residual
20: Level: advanced
22: .keywords: MG, default, multigrid, residual
24: .seealso: PCMGSetResidual()
25: @*/
26: PetscErrorCode PCMGDefaultResidual(Mat mat,Vec b,Vec x,Vec r)
27: {
31: MatMult(mat,x,r);
32: VecAYPX(r,-1.0,b);
33: return(0);
34: }
36: /* ---------------------------------------------------------------------------*/
40: /*@
41: PCMGGetCoarseSolve - Gets the solver context to be used on the coarse grid.
43: Not Collective
45: Input Parameter:
46: . pc - the multigrid context
48: Output Parameter:
49: . ksp - the coarse grid solver context
51: Level: advanced
53: .keywords: MG, multigrid, get, coarse grid
54: @*/
55: PetscErrorCode PCMGGetCoarseSolve(PC pc,KSP *ksp)
56: {
57: PC_MG *mg = (PC_MG*)pc->data;
58: PC_MG_Levels **mglevels = mg->levels;
62: *ksp = mglevels[0]->smoothd;
63: return(0);
64: }
68: /*@C
69: PCMGSetResidual - Sets the function to be used to calculate the residual
70: on the lth level.
72: Logically Collective on PC and Mat
74: Input Parameters:
75: + pc - the multigrid context
76: . l - the level (0 is coarsest) to supply
77: . residual - function used to form residual, if none is provided the previously provide one is used, if no
78: previous one were provided then PCMGDefaultResidual() is used
79: - mat - matrix associated with residual
81: Level: advanced
83: .keywords: MG, set, multigrid, residual, level
85: .seealso: PCMGDefaultResidual()
86: @*/
87: PetscErrorCode PCMGSetResidual(PC pc,PetscInt l,PetscErrorCode (*residual)(Mat,Vec,Vec,Vec),Mat mat)
88: {
89: PC_MG *mg = (PC_MG*)pc->data;
90: PC_MG_Levels **mglevels = mg->levels;
95: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
96: if (residual) {
97: mglevels[l]->residual = residual;
98: } if (!mglevels[l]->residual) {
99: mglevels[l]->residual = PCMGDefaultResidual;
100: }
101: if (mat) {PetscObjectReference((PetscObject)mat);}
102: MatDestroy(&mglevels[l]->A);
103: mglevels[l]->A = mat;
104: return(0);
105: }
109: /*@
110: PCMGSetInterpolation - Sets the function to be used to calculate the
111: interpolation from l-1 to the lth level
113: Logically Collective on PC and Mat
115: Input Parameters:
116: + pc - the multigrid context
117: . mat - the interpolation operator
118: - l - the level (0 is coarsest) to supply [do not supply 0]
120: Level: advanced
122: Notes:
123: Usually this is the same matrix used also to set the restriction
124: for the same level.
126: One can pass in the interpolation matrix or its transpose; PETSc figures
127: out from the matrix size which one it is.
129: .keywords: multigrid, set, interpolate, level
131: .seealso: PCMGSetRestriction()
132: @*/
133: PetscErrorCode PCMGSetInterpolation(PC pc,PetscInt l,Mat mat)
134: {
135: PC_MG *mg = (PC_MG*)pc->data;
136: PC_MG_Levels **mglevels = mg->levels;
141: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
142: if (!l) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Do not set interpolation routine for coarsest level");
143: PetscObjectReference((PetscObject)mat);
144: MatDestroy(&mglevels[l]->interpolate);
145: mglevels[l]->interpolate = mat;
146: return(0);
147: }
151: /*@
152: PCMGGetInterpolation - Gets the function to be used to calculate the
153: interpolation from l-1 to the lth level
155: Logically Collective on PC
157: Input Parameters:
158: + pc - the multigrid context
159: - l - the level (0 is coarsest) to supply [Do not supply 0]
161: Output Parameter:
162: . mat - the interpolation matrix
164: Level: advanced
166: .keywords: MG, get, multigrid, interpolation, level
168: .seealso: PCMGGetRestriction(), PCMGSetInterpolation(), PCMGGetRScale()
169: @*/
170: PetscErrorCode PCMGGetInterpolation(PC pc,PetscInt l,Mat *mat)
171: {
172: PC_MG *mg = (PC_MG*)pc->data;
173: PC_MG_Levels **mglevels = mg->levels;
179: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
180: if (l <= 0 || mg->nlevels <= l) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Level %D must be in range {1,...,%D}",l,mg->nlevels-1);
181: if (!mglevels[l]->interpolate) {
182: if (!mglevels[l]->restrct) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must call PCMGSetInterpolation() or PCMGSetInterpolation()");
183: PCMGSetInterpolation(pc,l,mglevels[l]->restrct);
184: }
185: *mat = mglevels[l]->interpolate;
186: return(0);
187: }
191: /*@
192: PCMGSetRestriction - Sets the function to be used to restrict vector
193: from level l to l-1.
195: Logically Collective on PC and Mat
197: Input Parameters:
198: + pc - the multigrid context
199: . l - the level (0 is coarsest) to supply [Do not supply 0]
200: - mat - the restriction matrix
202: Level: advanced
204: Notes:
205: Usually this is the same matrix used also to set the interpolation
206: for the same level.
208: One can pass in the interpolation matrix or its transpose; PETSc figures
209: out from the matrix size which one it is.
211: If you do not set this, the transpose of the Mat set with PCMGSetInterpolation()
212: is used.
214: .keywords: MG, set, multigrid, restriction, level
216: .seealso: PCMGSetInterpolation()
217: @*/
218: PetscErrorCode PCMGSetRestriction(PC pc,PetscInt l,Mat mat)
219: {
221: PC_MG *mg = (PC_MG*)pc->data;
222: PC_MG_Levels **mglevels = mg->levels;
227: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
228: if (!l) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Do not set restriction routine for coarsest level");
229: PetscObjectReference((PetscObject)mat);
230: MatDestroy(&mglevels[l]->restrct);
231: mglevels[l]->restrct = mat;
232: return(0);
233: }
237: /*@
238: PCMGGetRestriction - Gets the function to be used to restrict vector
239: from level l to l-1.
241: Logically Collective on PC and Mat
243: Input Parameters:
244: + pc - the multigrid context
245: - l - the level (0 is coarsest) to supply [Do not supply 0]
247: Output Parameter:
248: . mat - the restriction matrix
250: Level: advanced
252: .keywords: MG, get, multigrid, restriction, level
254: .seealso: PCMGGetInterpolation(), PCMGSetRestriction(), PCMGGetRScale()
255: @*/
256: PetscErrorCode PCMGGetRestriction(PC pc,PetscInt l,Mat *mat)
257: {
258: PC_MG *mg = (PC_MG*)pc->data;
259: PC_MG_Levels **mglevels = mg->levels;
265: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
266: if (l <= 0 || mg->nlevels <= l) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Level %D must be in range {1,...,%D}",l,mg->nlevels-1);
267: if (!mglevels[l]->restrct) {
268: if (!mglevels[l]->interpolate) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must call PCMGSetRestriction() or PCMGSetInterpolation()");
269: PCMGSetRestriction(pc,l,mglevels[l]->interpolate);
270: }
271: *mat = mglevels[l]->restrct;
272: return(0);
273: }
277: /*@
278: PCMGSetRScale - Sets the pointwise scaling for the restriction operator from level l to l-1.
280: Logically Collective on PC and Vec
282: Input Parameters:
283: + pc - the multigrid context
284: - l - the level (0 is coarsest) to supply [Do not supply 0]
285: . rscale - the scaling
287: Level: advanced
289: Notes:
290: 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.
292: .keywords: MG, set, multigrid, restriction, level
294: .seealso: PCMGSetInterpolation(), PCMGSetRestriction(), PCMGGetRScale()
295: @*/
296: PetscErrorCode PCMGSetRScale(PC pc,PetscInt l,Vec rscale)
297: {
299: PC_MG *mg = (PC_MG*)pc->data;
300: PC_MG_Levels **mglevels = mg->levels;
304: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
305: if (l <= 0 || mg->nlevels <= l) SETERRQ2(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Level %D must be in range {1,...,%D}",l,mg->nlevels-1);
306: PetscObjectReference((PetscObject)rscale);
307: 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(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
342: if (l <= 0 || mg->nlevels <= l) SETERRQ2(((PetscObject)pc)->comm,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) {fine = X; coarse = Y;}
351: else if (N < M) {fine = Y; coarse = X;}
352: else SETERRQ(((PetscObject)R)->comm,PETSC_ERR_SUP,"Restriction matrix is square, cannot determine which Vec is coarser");
353: VecSet(fine,1.);
354: MatRestrict(R,fine,coarse);
355: VecDestroy(&fine);
356: VecReciprocal(coarse);
357: mglevels[l]->rscale = coarse;
358: }
359: *rscale = mglevels[l]->rscale;
360: return(0);
361: }
365: /*@
366: PCMGGetSmoother - Gets the KSP context to be used as smoother for
367: both pre- and post-smoothing. Call both PCMGGetSmootherUp() and
368: PCMGGetSmootherDown() to use different functions for pre- and
369: post-smoothing.
371: Not Collective, KSP returned is parallel if PC is
373: Input Parameters:
374: + pc - the multigrid context
375: - l - the level (0 is coarsest) to supply
377: Ouput Parameters:
378: . ksp - the smoother
380: Level: advanced
382: .keywords: MG, get, multigrid, level, smoother, pre-smoother, post-smoother
384: .seealso: PCMGGetSmootherUp(), PCMGGetSmootherDown()
385: @*/
386: PetscErrorCode PCMGGetSmoother(PC pc,PetscInt l,KSP *ksp)
387: {
388: PC_MG *mg = (PC_MG*)pc->data;
389: PC_MG_Levels **mglevels = mg->levels;
393: *ksp = mglevels[l]->smoothd;
394: return(0);
395: }
399: /*@
400: PCMGGetSmootherUp - Gets the KSP context to be used as smoother after
401: coarse grid correction (post-smoother).
403: Not Collective, KSP returned is parallel if PC is
405: Input Parameters:
406: + pc - the multigrid context
407: - l - the level (0 is coarsest) to supply
409: Ouput Parameters:
410: . ksp - the smoother
412: Level: advanced
414: Notes: calling this will result in a different pre and post smoother so you may need to
415: set options on the pre smoother also
417: .keywords: MG, multigrid, get, smoother, up, post-smoother, level
419: .seealso: PCMGGetSmootherUp(), PCMGGetSmootherDown()
420: @*/
421: PetscErrorCode PCMGGetSmootherUp(PC pc,PetscInt l,KSP *ksp)
422: {
423: PC_MG *mg = (PC_MG*)pc->data;
424: PC_MG_Levels **mglevels = mg->levels;
426: const char *prefix;
427: MPI_Comm comm;
431: /*
432: This is called only if user wants a different pre-smoother from post.
433: Thus we check if a different one has already been allocated,
434: if not we allocate it.
435: */
436: if (!l) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"There is no such thing as a up smoother on the coarse grid");
437: if (mglevels[l]->smoothu == mglevels[l]->smoothd) {
438: const KSPType ksptype;
439: const PCType pctype;
440: PC ipc;
441: PetscReal rtol,abstol,dtol;
442: PetscInt maxits;
443: KSPNormType normtype;
444: PetscObjectGetComm((PetscObject)mglevels[l]->smoothd,&comm);
445: KSPGetOptionsPrefix(mglevels[l]->smoothd,&prefix);
446: KSPGetTolerances(mglevels[l]->smoothd,&rtol,&abstol,&dtol,&maxits);
447: KSPGetType(mglevels[l]->smoothd,&ksptype);
448: KSPGetNormType(mglevels[l]->smoothd,&normtype);
449: KSPGetPC(mglevels[l]->smoothd,&ipc);
450: PCGetType(ipc,&pctype);
452: KSPCreate(comm,&mglevels[l]->smoothu);
453: PetscObjectIncrementTabLevel((PetscObject)mglevels[l]->smoothu,(PetscObject)pc,mglevels[0]->levels-l);
454: KSPSetOptionsPrefix(mglevels[l]->smoothu,prefix);
455: KSPSetTolerances(mglevels[l]->smoothu,rtol,abstol,dtol,maxits);
456: KSPSetType(mglevels[l]->smoothu,ksptype);
457: KSPSetNormType(mglevels[l]->smoothu,normtype);
458: KSPSetConvergenceTest(mglevels[l]->smoothu,KSPSkipConverged,PETSC_NULL,PETSC_NULL);
459: KSPGetPC(mglevels[l]->smoothu,&ipc);
460: PCSetType(ipc,pctype);
461: PetscLogObjectParent(pc,mglevels[l]->smoothu);
462: }
463: if (ksp) *ksp = mglevels[l]->smoothu;
464: return(0);
465: }
469: /*@
470: PCMGGetSmootherDown - Gets the KSP context to be used as smoother before
471: coarse grid correction (pre-smoother).
473: Not Collective, KSP returned is parallel if PC is
475: Input Parameters:
476: + pc - the multigrid context
477: - l - the level (0 is coarsest) to supply
479: Ouput Parameters:
480: . ksp - the smoother
482: Level: advanced
484: Notes: calling this will result in a different pre and post smoother so you may need to
485: set options on the post smoother also
487: .keywords: MG, multigrid, get, smoother, down, pre-smoother, level
489: .seealso: PCMGGetSmootherUp(), PCMGGetSmoother()
490: @*/
491: PetscErrorCode PCMGGetSmootherDown(PC pc,PetscInt l,KSP *ksp)
492: {
494: PC_MG *mg = (PC_MG*)pc->data;
495: PC_MG_Levels **mglevels = mg->levels;
499: /* make sure smoother up and down are different */
500: if (l) {
501: PCMGGetSmootherUp(pc,l,PETSC_NULL);
502: }
503: *ksp = mglevels[l]->smoothd;
504: return(0);
505: }
509: /*@
510: PCMGSetCyclesOnLevel - Sets the number of cycles to run on this level.
512: Logically Collective on PC
514: Input Parameters:
515: + pc - the multigrid context
516: . l - the level (0 is coarsest) this is to be used for
517: - n - the number of cycles
519: Level: advanced
521: .keywords: MG, multigrid, set, cycles, V-cycle, W-cycle, level
523: .seealso: PCMGSetCycles()
524: @*/
525: PetscErrorCode PCMGSetCyclesOnLevel(PC pc,PetscInt l,PetscInt c)
526: {
527: PC_MG *mg = (PC_MG*)pc->data;
528: PC_MG_Levels **mglevels = mg->levels;
532: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
535: mglevels[l]->cycles = c;
536: return(0);
537: }
541: /*@
542: PCMGSetRhs - Sets the vector space to be used to store the right-hand side
543: on a particular level.
545: Logically Collective on PC and Vec
547: Input Parameters:
548: + pc - the multigrid context
549: . l - the level (0 is coarsest) this is to be used for
550: - c - the space
552: Level: advanced
554: Notes: If this is not provided PETSc will automatically generate one.
556: You do not need to keep a reference to this vector if you do
557: not need it PCDestroy() will properly free it.
559: .keywords: MG, multigrid, set, right-hand-side, rhs, level
561: .seealso: PCMGSetX(), PCMGSetR()
562: @*/
563: PetscErrorCode PCMGSetRhs(PC pc,PetscInt l,Vec c)
564: {
566: PC_MG *mg = (PC_MG*)pc->data;
567: PC_MG_Levels **mglevels = mg->levels;
571: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
572: if (l == mglevels[0]->levels-1) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_INCOMP,"Do not set rhs for finest level");
573: PetscObjectReference((PetscObject)c);
574: VecDestroy(&mglevels[l]->b);
575: mglevels[l]->b = c;
576: return(0);
577: }
581: /*@
582: PCMGSetX - Sets the vector space to be used to store the solution on a
583: particular level.
585: Logically Collective on PC and Vec
587: Input Parameters:
588: + pc - the multigrid context
589: . l - the level (0 is coarsest) this is to be used for (do not supply the finest level)
590: - c - the space
592: Level: advanced
594: Notes: If this is not provided PETSc will automatically generate one.
596: You do not need to keep a reference to this vector if you do
597: not need it PCDestroy() will properly free it.
599: .keywords: MG, multigrid, set, solution, level
601: .seealso: PCMGSetRhs(), PCMGSetR()
602: @*/
603: PetscErrorCode PCMGSetX(PC pc,PetscInt l,Vec c)
604: {
606: PC_MG *mg = (PC_MG*)pc->data;
607: PC_MG_Levels **mglevels = mg->levels;
611: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
612: if (l == mglevels[0]->levels-1) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_INCOMP,"Do not set x for finest level");
613: PetscObjectReference((PetscObject)c);
614: VecDestroy(&mglevels[l]->x);
615: mglevels[l]->x = c;
616: return(0);
617: }
621: /*@
622: PCMGSetR - Sets the vector space to be used to store the residual on a
623: particular level.
625: Logically Collective on PC and Vec
627: Input Parameters:
628: + pc - the multigrid context
629: . l - the level (0 is coarsest) this is to be used for
630: - c - the space
632: Level: advanced
634: Notes: If this is not provided PETSc will automatically generate one.
636: You do not need to keep a reference to this vector if you do
637: not need it PCDestroy() will properly free it.
639: .keywords: MG, multigrid, set, residual, level
640: @*/
641: PetscErrorCode PCMGSetR(PC pc,PetscInt l,Vec c)
642: {
644: PC_MG *mg = (PC_MG*)pc->data;
645: PC_MG_Levels **mglevels = mg->levels;
649: if (!mglevels) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
650: if (!l) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Need not set residual vector for coarse grid");
651: PetscObjectReference((PetscObject)c);
652: VecDestroy(&mglevels[l]->r);
653: mglevels[l]->r = c;
654: return(0);
655: }