Actual source code: mgfunc.c
petsc-3.5.4 2015-05-23
2: #include <../src/ksp/pc/impls/mg/mgimpl.h> /*I "petscksp.h" I*/
4: /* ---------------------------------------------------------------------------*/
7: /*@C
8: PCMGResidualDefault - 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
17: Output Parameter:
18: . r - location to store the residual
20: Level: developer
22: .keywords: MG, default, multigrid, residual
24: .seealso: PCMGSetResidual()
25: @*/
26: PetscErrorCode PCMGResidualDefault(Mat mat,Vec b,Vec x,Vec r)
27: {
31: MatResidual(mat,b,x,r);
32: return(0);
33: }
37: /*@
38: PCMGGetCoarseSolve - Gets the solver context to be used on the coarse grid.
40: Not Collective
42: Input Parameter:
43: . pc - the multigrid context
45: Output Parameter:
46: . ksp - the coarse grid solver context
48: Level: advanced
50: .keywords: MG, multigrid, get, coarse grid
51: @*/
52: PetscErrorCode PCMGGetCoarseSolve(PC pc,KSP *ksp)
53: {
54: PC_MG *mg = (PC_MG*)pc->data;
55: PC_MG_Levels **mglevels = mg->levels;
59: *ksp = mglevels[0]->smoothd;
60: return(0);
61: }
65: /*@C
66: PCMGSetResidual - Sets the function to be used to calculate the residual
67: on the lth level.
69: Logically Collective on PC and Mat
71: Input Parameters:
72: + pc - the multigrid context
73: . l - the level (0 is coarsest) to supply
74: . residual - function used to form residual, if none is provided the previously provide one is used, if no
75: previous one were provided then a default is used
76: - mat - matrix associated with residual
78: Level: advanced
80: .keywords: MG, set, multigrid, residual, level
82: .seealso: PCMGResidualDefault()
83: @*/
84: PetscErrorCode PCMGSetResidual(PC pc,PetscInt l,PetscErrorCode (*residual)(Mat,Vec,Vec,Vec),Mat mat)
85: {
86: PC_MG *mg = (PC_MG*)pc->data;
87: PC_MG_Levels **mglevels = mg->levels;
92: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
93: if (residual) mglevels[l]->residual = residual;
94: if (!mglevels[l]->residual) mglevels[l]->residual = PCMGResidualDefault;
95: if (mat) {PetscObjectReference((PetscObject)mat);}
96: MatDestroy(&mglevels[l]->A);
97: mglevels[l]->A = mat;
98: return(0);
99: }
103: /*@
104: PCMGSetInterpolation - Sets the function to be used to calculate the
105: interpolation from l-1 to the lth level
107: Logically Collective on PC and Mat
109: Input Parameters:
110: + pc - the multigrid context
111: . mat - the interpolation operator
112: - l - the level (0 is coarsest) to supply [do not supply 0]
114: Level: advanced
116: Notes:
117: Usually this is the same matrix used also to set the restriction
118: for the same level.
120: One can pass in the interpolation matrix or its transpose; PETSc figures
121: out from the matrix size which one it is.
123: .keywords: multigrid, set, interpolate, level
125: .seealso: PCMGSetRestriction()
126: @*/
127: PetscErrorCode PCMGSetInterpolation(PC pc,PetscInt l,Mat mat)
128: {
129: PC_MG *mg = (PC_MG*)pc->data;
130: PC_MG_Levels **mglevels = mg->levels;
135: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
136: if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Do not set interpolation routine for coarsest level");
137: PetscObjectReference((PetscObject)mat);
138: MatDestroy(&mglevels[l]->interpolate);
140: mglevels[l]->interpolate = mat;
141: return(0);
142: }
146: /*@
147: PCMGGetInterpolation - Gets the function to be used to calculate the
148: interpolation from l-1 to the lth level
150: Logically Collective on PC
152: Input Parameters:
153: + pc - the multigrid context
154: - l - the level (0 is coarsest) to supply [Do not supply 0]
156: Output Parameter:
157: . mat - the interpolation matrix
159: Level: advanced
161: .keywords: MG, get, multigrid, interpolation, level
163: .seealso: PCMGGetRestriction(), PCMGSetInterpolation(), PCMGGetRScale()
164: @*/
165: PetscErrorCode PCMGGetInterpolation(PC pc,PetscInt l,Mat *mat)
166: {
167: PC_MG *mg = (PC_MG*)pc->data;
168: PC_MG_Levels **mglevels = mg->levels;
174: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
175: 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);
176: if (!mglevels[l]->interpolate) {
177: if (!mglevels[l]->restrct) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call PCMGSetInterpolation() or PCMGSetInterpolation()");
178: PCMGSetInterpolation(pc,l,mglevels[l]->restrct);
179: }
180: *mat = mglevels[l]->interpolate;
181: return(0);
182: }
186: /*@
187: PCMGSetRestriction - Sets the function to be used to restrict vector
188: from level l to l-1.
190: Logically Collective on PC and Mat
192: Input Parameters:
193: + pc - the multigrid context
194: . l - the level (0 is coarsest) to supply [Do not supply 0]
195: - mat - the restriction matrix
197: Level: advanced
199: Notes:
200: Usually this is the same matrix used also to set the interpolation
201: for the same level.
203: One can pass in the interpolation matrix or its transpose; PETSc figures
204: out from the matrix size which one it is.
206: If you do not set this, the transpose of the Mat set with PCMGSetInterpolation()
207: is used.
209: .keywords: MG, set, multigrid, restriction, level
211: .seealso: PCMGSetInterpolation()
212: @*/
213: PetscErrorCode PCMGSetRestriction(PC pc,PetscInt l,Mat mat)
214: {
216: PC_MG *mg = (PC_MG*)pc->data;
217: PC_MG_Levels **mglevels = mg->levels;
222: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
223: if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Do not set restriction routine for coarsest level");
224: PetscObjectReference((PetscObject)mat);
225: MatDestroy(&mglevels[l]->restrct);
227: mglevels[l]->restrct = mat;
228: return(0);
229: }
233: /*@
234: PCMGGetRestriction - Gets the function to be used to restrict vector
235: from level l to l-1.
237: Logically Collective on PC and Mat
239: Input Parameters:
240: + pc - the multigrid context
241: - l - the level (0 is coarsest) to supply [Do not supply 0]
243: Output Parameter:
244: . mat - the restriction matrix
246: Level: advanced
248: .keywords: MG, get, multigrid, restriction, level
250: .seealso: PCMGGetInterpolation(), PCMGSetRestriction(), PCMGGetRScale()
251: @*/
252: PetscErrorCode PCMGGetRestriction(PC pc,PetscInt l,Mat *mat)
253: {
254: PC_MG *mg = (PC_MG*)pc->data;
255: PC_MG_Levels **mglevels = mg->levels;
261: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
262: 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);
263: if (!mglevels[l]->restrct) {
264: if (!mglevels[l]->interpolate) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call PCMGSetRestriction() or PCMGSetInterpolation()");
265: PCMGSetRestriction(pc,l,mglevels[l]->interpolate);
266: }
267: *mat = mglevels[l]->restrct;
268: return(0);
269: }
273: /*@
274: PCMGSetRScale - Sets the pointwise scaling for the restriction operator from level l to l-1.
276: Logically Collective on PC and Vec
278: Input Parameters:
279: + pc - the multigrid context
280: - l - the level (0 is coarsest) to supply [Do not supply 0]
281: . rscale - the scaling
283: Level: advanced
285: Notes:
286: 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.
288: .keywords: MG, set, multigrid, restriction, level
290: .seealso: PCMGSetInterpolation(), PCMGSetRestriction(), PCMGGetRScale()
291: @*/
292: PetscErrorCode PCMGSetRScale(PC pc,PetscInt l,Vec rscale)
293: {
295: PC_MG *mg = (PC_MG*)pc->data;
296: PC_MG_Levels **mglevels = mg->levels;
300: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
301: 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);
302: PetscObjectReference((PetscObject)rscale);
303: VecDestroy(&mglevels[l]->rscale);
305: mglevels[l]->rscale = rscale;
306: return(0);
307: }
311: /*@
312: PCMGGetRScale - Gets the pointwise scaling for the restriction operator from level l to l-1.
314: Collective on PC
316: Input Parameters:
317: + pc - the multigrid context
318: . rscale - the scaling
319: - l - the level (0 is coarsest) to supply [Do not supply 0]
321: Level: advanced
323: Notes:
324: 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.
326: .keywords: MG, set, multigrid, restriction, level
328: .seealso: PCMGSetInterpolation(), PCMGGetRestriction()
329: @*/
330: PetscErrorCode PCMGGetRScale(PC pc,PetscInt l,Vec *rscale)
331: {
333: PC_MG *mg = (PC_MG*)pc->data;
334: PC_MG_Levels **mglevels = mg->levels;
338: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
339: 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);
340: if (!mglevels[l]->rscale) {
341: Mat R;
342: Vec X,Y,coarse,fine;
343: PetscInt M,N;
344: PCMGGetRestriction(pc,l,&R);
345: MatGetVecs(R,&X,&Y);
346: MatGetSize(R,&M,&N);
347: if (M < N) {
348: fine = X;
349: coarse = Y;
350: } else if (N < M) {
351: fine = Y; coarse = X;
352: } else SETERRQ(PetscObjectComm((PetscObject)R),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(PetscObjectComm((PetscObject)pc),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: KSPType ksptype;
439: 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,KSPConvergedSkip,NULL,NULL);
459: KSPGetPC(mglevels[l]->smoothu,&ipc);
460: PCSetType(ipc,pctype);
461: PetscLogObjectParent((PetscObject)pc,(PetscObject)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,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(PetscObjectComm((PetscObject)pc),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(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
572: if (l == mglevels[0]->levels-1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Do not set rhs for finest level");
573: PetscObjectReference((PetscObject)c);
574: VecDestroy(&mglevels[l]->b);
576: mglevels[l]->b = c;
577: return(0);
578: }
582: /*@
583: PCMGSetX - Sets the vector space to be used to store the solution on a
584: particular level.
586: Logically Collective on PC and Vec
588: Input Parameters:
589: + pc - the multigrid context
590: . l - the level (0 is coarsest) this is to be used for (do not supply the finest level)
591: - c - the space
593: Level: advanced
595: Notes: If this is not provided PETSc will automatically generate one.
597: You do not need to keep a reference to this vector if you do
598: not need it PCDestroy() will properly free it.
600: .keywords: MG, multigrid, set, solution, level
602: .seealso: PCMGSetRhs(), PCMGSetR()
603: @*/
604: PetscErrorCode PCMGSetX(PC pc,PetscInt l,Vec c)
605: {
607: PC_MG *mg = (PC_MG*)pc->data;
608: PC_MG_Levels **mglevels = mg->levels;
612: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
613: if (l == mglevels[0]->levels-1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Do not set x for finest level");
614: PetscObjectReference((PetscObject)c);
615: VecDestroy(&mglevels[l]->x);
617: mglevels[l]->x = c;
618: return(0);
619: }
623: /*@
624: PCMGSetR - Sets the vector space to be used to store the residual on a
625: particular level.
627: Logically Collective on PC and Vec
629: Input Parameters:
630: + pc - the multigrid context
631: . l - the level (0 is coarsest) this is to be used for
632: - c - the space
634: Level: advanced
636: Notes: If this is not provided PETSc will automatically generate one.
638: You do not need to keep a reference to this vector if you do
639: not need it PCDestroy() will properly free it.
641: .keywords: MG, multigrid, set, residual, level
642: @*/
643: PetscErrorCode PCMGSetR(PC pc,PetscInt l,Vec c)
644: {
646: PC_MG *mg = (PC_MG*)pc->data;
647: PC_MG_Levels **mglevels = mg->levels;
651: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
652: if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Need not set residual vector for coarse grid");
653: PetscObjectReference((PetscObject)c);
654: VecDestroy(&mglevels[l]->r);
656: mglevels[l]->r = c;
657: return(0);
658: }