Actual source code: mgfunc.c
petsc-3.8.4 2018-03-24
2: #include <petsc/private/pcmgimpl.h>
4: /* ---------------------------------------------------------------------------*/
5: /*@C
6: PCMGResidualDefault - Default routine to calculate the residual.
8: Collective on Mat and Vec
10: Input Parameters:
11: + mat - the matrix
12: . b - the right-hand-side
13: - x - the approximate solution
15: Output Parameter:
16: . r - location to store the residual
18: Level: developer
20: .keywords: MG, default, multigrid, residual
22: .seealso: PCMGSetResidual()
23: @*/
24: PetscErrorCode PCMGResidualDefault(Mat mat,Vec b,Vec x,Vec r)
25: {
29: MatResidual(mat,b,x,r);
30: return(0);
31: }
33: /*@
34: PCMGGetCoarseSolve - Gets the solver context to be used on the coarse grid.
36: Not Collective
38: Input Parameter:
39: . pc - the multigrid context
41: Output Parameter:
42: . ksp - the coarse grid solver context
44: Level: advanced
46: .keywords: MG, multigrid, get, coarse grid
47: @*/
48: PetscErrorCode PCMGGetCoarseSolve(PC pc,KSP *ksp)
49: {
50: PC_MG *mg = (PC_MG*)pc->data;
51: PC_MG_Levels **mglevels = mg->levels;
55: *ksp = mglevels[0]->smoothd;
56: return(0);
57: }
59: /*@C
60: PCMGSetResidual - Sets the function to be used to calculate the residual
61: on the lth level.
63: Logically Collective on PC and Mat
65: Input Parameters:
66: + pc - the multigrid context
67: . l - the level (0 is coarsest) to supply
68: . residual - function used to form residual, if none is provided the previously provide one is used, if no
69: previous one were provided then a default is used
70: - mat - matrix associated with residual
72: Level: advanced
74: .keywords: MG, set, multigrid, residual, level
76: .seealso: PCMGResidualDefault()
77: @*/
78: PetscErrorCode PCMGSetResidual(PC pc,PetscInt l,PetscErrorCode (*residual)(Mat,Vec,Vec,Vec),Mat mat)
79: {
80: PC_MG *mg = (PC_MG*)pc->data;
81: PC_MG_Levels **mglevels = mg->levels;
86: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
87: if (residual) mglevels[l]->residual = residual;
88: if (!mglevels[l]->residual) mglevels[l]->residual = PCMGResidualDefault;
89: if (mat) {PetscObjectReference((PetscObject)mat);}
90: MatDestroy(&mglevels[l]->A);
91: mglevels[l]->A = mat;
92: return(0);
93: }
95: /*@
96: PCMGSetInterpolation - Sets the function to be used to calculate the
97: interpolation from l-1 to the lth level
99: Logically Collective on PC and Mat
101: Input Parameters:
102: + pc - the multigrid context
103: . mat - the interpolation operator
104: - l - the level (0 is coarsest) to supply [do not supply 0]
106: Level: advanced
108: Notes:
109: Usually this is the same matrix used also to set the restriction
110: for the same level.
112: One can pass in the interpolation matrix or its transpose; PETSc figures
113: out from the matrix size which one it is.
115: .keywords: multigrid, set, interpolate, level
117: .seealso: PCMGSetRestriction()
118: @*/
119: PetscErrorCode PCMGSetInterpolation(PC pc,PetscInt l,Mat mat)
120: {
121: PC_MG *mg = (PC_MG*)pc->data;
122: PC_MG_Levels **mglevels = mg->levels;
127: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
128: if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Do not set interpolation routine for coarsest level");
129: PetscObjectReference((PetscObject)mat);
130: MatDestroy(&mglevels[l]->interpolate);
132: mglevels[l]->interpolate = mat;
133: return(0);
134: }
136: /*@
137: PCMGGetInterpolation - Gets the function to be used to calculate the
138: interpolation from l-1 to the lth level
140: Logically Collective on PC
142: Input Parameters:
143: + pc - the multigrid context
144: - l - the level (0 is coarsest) to supply [Do not supply 0]
146: Output Parameter:
147: . mat - the interpolation matrix, can be NULL
149: Level: advanced
151: .keywords: MG, get, multigrid, interpolation, level
153: .seealso: PCMGGetRestriction(), PCMGSetInterpolation(), PCMGGetRScale()
154: @*/
155: PetscErrorCode PCMGGetInterpolation(PC pc,PetscInt l,Mat *mat)
156: {
157: PC_MG *mg = (PC_MG*)pc->data;
158: PC_MG_Levels **mglevels = mg->levels;
164: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
165: 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);
166: if (!mglevels[l]->interpolate) {
167: if (!mglevels[l]->restrct) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call PCMGSetInterpolation() or PCMGSetRestriction()");
168: PCMGSetInterpolation(pc,l,mglevels[l]->restrct);
169: }
170: if (mat) *mat = mglevels[l]->interpolate;
171: return(0);
172: }
174: /*@
175: PCMGSetRestriction - Sets the function to be used to restrict vector
176: from level l to l-1.
178: Logically Collective on PC and Mat
180: Input Parameters:
181: + pc - the multigrid context
182: . l - the level (0 is coarsest) to supply [Do not supply 0]
183: - mat - the restriction matrix
185: Level: advanced
187: Notes:
188: Usually this is the same matrix used also to set the interpolation
189: for the same level.
191: One can pass in the interpolation matrix or its transpose; PETSc figures
192: out from the matrix size which one it is.
194: If you do not set this, the transpose of the Mat set with PCMGSetInterpolation()
195: is used.
197: .keywords: MG, set, multigrid, restriction, level
199: .seealso: PCMGSetInterpolation()
200: @*/
201: PetscErrorCode PCMGSetRestriction(PC pc,PetscInt l,Mat mat)
202: {
204: PC_MG *mg = (PC_MG*)pc->data;
205: PC_MG_Levels **mglevels = mg->levels;
210: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
211: if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Do not set restriction routine for coarsest level");
212: PetscObjectReference((PetscObject)mat);
213: MatDestroy(&mglevels[l]->restrct);
215: mglevels[l]->restrct = mat;
216: return(0);
217: }
219: /*@
220: PCMGGetRestriction - Gets the function to be used to restrict vector
221: from level l to l-1.
223: Logically Collective on PC and Mat
225: Input Parameters:
226: + pc - the multigrid context
227: - l - the level (0 is coarsest) to supply [Do not supply 0]
229: Output Parameter:
230: . mat - the restriction matrix
232: Level: advanced
234: .keywords: MG, get, multigrid, restriction, level
236: .seealso: PCMGGetInterpolation(), PCMGSetRestriction(), PCMGGetRScale()
237: @*/
238: PetscErrorCode PCMGGetRestriction(PC pc,PetscInt l,Mat *mat)
239: {
240: PC_MG *mg = (PC_MG*)pc->data;
241: PC_MG_Levels **mglevels = mg->levels;
247: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
248: 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);
249: if (!mglevels[l]->restrct) {
250: if (!mglevels[l]->interpolate) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must call PCMGSetRestriction() or PCMGSetInterpolation()");
251: PCMGSetRestriction(pc,l,mglevels[l]->interpolate);
252: }
253: if (mat) *mat = mglevels[l]->restrct;
254: return(0);
255: }
257: /*@
258: PCMGSetRScale - Sets the pointwise scaling for the restriction operator from level l to l-1.
260: Logically Collective on PC and Vec
262: Input Parameters:
263: + pc - the multigrid context
264: - l - the level (0 is coarsest) to supply [Do not supply 0]
265: . rscale - the scaling
267: Level: advanced
269: Notes:
270: 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.
272: .keywords: MG, set, multigrid, restriction, level
274: .seealso: PCMGSetInterpolation(), PCMGSetRestriction(), PCMGGetRScale()
275: @*/
276: PetscErrorCode PCMGSetRScale(PC pc,PetscInt l,Vec rscale)
277: {
279: PC_MG *mg = (PC_MG*)pc->data;
280: PC_MG_Levels **mglevels = mg->levels;
284: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
285: 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);
286: PetscObjectReference((PetscObject)rscale);
287: VecDestroy(&mglevels[l]->rscale);
289: mglevels[l]->rscale = rscale;
290: return(0);
291: }
293: /*@
294: PCMGGetRScale - Gets the pointwise scaling for the restriction operator from level l to l-1.
296: Collective on PC
298: Input Parameters:
299: + pc - the multigrid context
300: . rscale - the scaling
301: - l - the level (0 is coarsest) to supply [Do not supply 0]
303: Level: advanced
305: Notes:
306: 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.
308: .keywords: MG, set, multigrid, restriction, level
310: .seealso: PCMGSetInterpolation(), PCMGGetRestriction()
311: @*/
312: PetscErrorCode PCMGGetRScale(PC pc,PetscInt l,Vec *rscale)
313: {
315: PC_MG *mg = (PC_MG*)pc->data;
316: PC_MG_Levels **mglevels = mg->levels;
320: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
321: 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);
322: if (!mglevels[l]->rscale) {
323: Mat R;
324: Vec X,Y,coarse,fine;
325: PetscInt M,N;
326: PCMGGetRestriction(pc,l,&R);
327: MatCreateVecs(R,&X,&Y);
328: MatGetSize(R,&M,&N);
329: if (M < N) {
330: fine = X;
331: coarse = Y;
332: } else if (N < M) {
333: fine = Y; coarse = X;
334: } else SETERRQ(PetscObjectComm((PetscObject)R),PETSC_ERR_SUP,"Restriction matrix is square, cannot determine which Vec is coarser");
335: VecSet(fine,1.);
336: MatRestrict(R,fine,coarse);
337: VecDestroy(&fine);
338: VecReciprocal(coarse);
339: mglevels[l]->rscale = coarse;
340: }
341: *rscale = mglevels[l]->rscale;
342: return(0);
343: }
345: /*@
346: PCMGGetSmoother - Gets the KSP context to be used as smoother for
347: both pre- and post-smoothing. Call both PCMGGetSmootherUp() and
348: PCMGGetSmootherDown() to use different functions for pre- and
349: post-smoothing.
351: Not Collective, KSP returned is parallel if PC is
353: Input Parameters:
354: + pc - the multigrid context
355: - l - the level (0 is coarsest) to supply
357: Ouput Parameters:
358: . ksp - the smoother
360: Notes:
361: Once you have called this routine, you can call KSPSetOperators(ksp,...) on the resulting ksp to provide the operators for the smoother for this level.
362: You can also modify smoother options by calling the various KSPSetXXX() options on this ksp. In addition you can call KSPGetPC(ksp,&pc)
363: and modify PC options for the smoother; for example PCSetType(pc,PCSOR); to use SOR smoothing.
365: Level: advanced
367: .keywords: MG, get, multigrid, level, smoother, pre-smoother, post-smoother
369: .seealso: PCMGGetSmootherUp(), PCMGGetSmootherDown()
370: @*/
371: PetscErrorCode PCMGGetSmoother(PC pc,PetscInt l,KSP *ksp)
372: {
373: PC_MG *mg = (PC_MG*)pc->data;
374: PC_MG_Levels **mglevels = mg->levels;
378: *ksp = mglevels[l]->smoothd;
379: return(0);
380: }
382: /*@
383: PCMGGetSmootherUp - Gets the KSP context to be used as smoother after
384: coarse grid correction (post-smoother).
386: Not Collective, KSP returned is parallel if PC is
388: Input Parameters:
389: + pc - the multigrid context
390: - l - the level (0 is coarsest) to supply
392: Ouput Parameters:
393: . ksp - the smoother
395: Level: advanced
397: Notes: calling this will result in a different pre and post smoother so you may need to
398: set options on the pre smoother also
400: .keywords: MG, multigrid, get, smoother, up, post-smoother, level
402: .seealso: PCMGGetSmootherUp(), PCMGGetSmootherDown()
403: @*/
404: PetscErrorCode PCMGGetSmootherUp(PC pc,PetscInt l,KSP *ksp)
405: {
406: PC_MG *mg = (PC_MG*)pc->data;
407: PC_MG_Levels **mglevels = mg->levels;
409: const char *prefix;
410: MPI_Comm comm;
414: /*
415: This is called only if user wants a different pre-smoother from post.
416: Thus we check if a different one has already been allocated,
417: if not we allocate it.
418: */
419: if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"There is no such thing as a up smoother on the coarse grid");
420: if (mglevels[l]->smoothu == mglevels[l]->smoothd) {
421: KSPType ksptype;
422: PCType pctype;
423: PC ipc;
424: PetscReal rtol,abstol,dtol;
425: PetscInt maxits;
426: KSPNormType normtype;
427: PetscObjectGetComm((PetscObject)mglevels[l]->smoothd,&comm);
428: KSPGetOptionsPrefix(mglevels[l]->smoothd,&prefix);
429: KSPGetTolerances(mglevels[l]->smoothd,&rtol,&abstol,&dtol,&maxits);
430: KSPGetType(mglevels[l]->smoothd,&ksptype);
431: KSPGetNormType(mglevels[l]->smoothd,&normtype);
432: KSPGetPC(mglevels[l]->smoothd,&ipc);
433: PCGetType(ipc,&pctype);
435: KSPCreate(comm,&mglevels[l]->smoothu);
436: KSPSetErrorIfNotConverged(mglevels[l]->smoothu,pc->erroriffailure);
437: PetscObjectIncrementTabLevel((PetscObject)mglevels[l]->smoothu,(PetscObject)pc,mglevels[0]->levels-l);
438: KSPSetOptionsPrefix(mglevels[l]->smoothu,prefix);
439: KSPSetTolerances(mglevels[l]->smoothu,rtol,abstol,dtol,maxits);
440: KSPSetType(mglevels[l]->smoothu,ksptype);
441: KSPSetNormType(mglevels[l]->smoothu,normtype);
442: KSPSetConvergenceTest(mglevels[l]->smoothu,KSPConvergedSkip,NULL,NULL);
443: KSPGetPC(mglevels[l]->smoothu,&ipc);
444: PCSetType(ipc,pctype);
445: PetscLogObjectParent((PetscObject)pc,(PetscObject)mglevels[l]->smoothu);
446: PetscObjectComposedDataSetInt((PetscObject) mglevels[l]->smoothu, PetscMGLevelId, mglevels[l]->level);
447: }
448: if (ksp) *ksp = mglevels[l]->smoothu;
449: return(0);
450: }
452: /*@
453: PCMGGetSmootherDown - Gets the KSP context to be used as smoother before
454: coarse grid correction (pre-smoother).
456: Not Collective, KSP returned is parallel if PC is
458: Input Parameters:
459: + pc - the multigrid context
460: - l - the level (0 is coarsest) to supply
462: Ouput Parameters:
463: . ksp - the smoother
465: Level: advanced
467: Notes: calling this will result in a different pre and post smoother so you may need to
468: set options on the post smoother also
470: .keywords: MG, multigrid, get, smoother, down, pre-smoother, level
472: .seealso: PCMGGetSmootherUp(), PCMGGetSmoother()
473: @*/
474: PetscErrorCode PCMGGetSmootherDown(PC pc,PetscInt l,KSP *ksp)
475: {
477: PC_MG *mg = (PC_MG*)pc->data;
478: PC_MG_Levels **mglevels = mg->levels;
482: /* make sure smoother up and down are different */
483: if (l) {
484: PCMGGetSmootherUp(pc,l,NULL);
485: }
486: *ksp = mglevels[l]->smoothd;
487: return(0);
488: }
490: /*@
491: PCMGSetCycleTypeOnLevel - Sets the type of cycle (aka cycle index) to run on the specified level.
493: Logically Collective on PC
495: Input Parameters:
496: + pc - the multigrid context
497: . l - the level (0 is coarsest)
498: - c - either PC_MG_CYCLE_V or PC_MG_CYCLE_W
500: Level: advanced
502: .keywords: MG, multigrid, set, cycles, V-cycle, W-cycle, level
504: .seealso: PCMGSetCycleType()
505: @*/
506: PetscErrorCode PCMGSetCycleTypeOnLevel(PC pc,PetscInt l,PCMGCycleType c)
507: {
508: PC_MG *mg = (PC_MG*)pc->data;
509: PC_MG_Levels **mglevels = mg->levels;
513: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
516: mglevels[l]->cycles = c;
517: return(0);
518: }
520: /*@
521: PCMGSetRhs - Sets the vector space to be used to store the right-hand side
522: on a particular level.
524: Logically Collective on PC and Vec
526: Input Parameters:
527: + pc - the multigrid context
528: . l - the level (0 is coarsest) this is to be used for
529: - c - the space
531: Level: advanced
533: Notes: If this is not provided PETSc will automatically generate one.
535: You do not need to keep a reference to this vector if you do
536: not need it PCDestroy() will properly free it.
538: .keywords: MG, multigrid, set, right-hand-side, rhs, level
540: .seealso: PCMGSetX(), PCMGSetR()
541: @*/
542: PetscErrorCode PCMGSetRhs(PC pc,PetscInt l,Vec c)
543: {
545: PC_MG *mg = (PC_MG*)pc->data;
546: PC_MG_Levels **mglevels = mg->levels;
550: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
551: if (l == mglevels[0]->levels-1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Do not set rhs for finest level");
552: PetscObjectReference((PetscObject)c);
553: VecDestroy(&mglevels[l]->b);
555: mglevels[l]->b = c;
556: return(0);
557: }
559: /*@
560: PCMGSetX - Sets the vector space to be used to store the solution on a
561: particular level.
563: Logically Collective on PC and Vec
565: Input Parameters:
566: + pc - the multigrid context
567: . l - the level (0 is coarsest) this is to be used for (do not supply the finest level)
568: - c - the space
570: Level: advanced
572: Notes: If this is not provided PETSc will automatically generate one.
574: You do not need to keep a reference to this vector if you do
575: not need it PCDestroy() will properly free it.
577: .keywords: MG, multigrid, set, solution, level
579: .seealso: PCMGSetRhs(), PCMGSetR()
580: @*/
581: PetscErrorCode PCMGSetX(PC pc,PetscInt l,Vec c)
582: {
584: PC_MG *mg = (PC_MG*)pc->data;
585: PC_MG_Levels **mglevels = mg->levels;
589: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
590: if (l == mglevels[0]->levels-1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Do not set x for finest level");
591: PetscObjectReference((PetscObject)c);
592: VecDestroy(&mglevels[l]->x);
594: mglevels[l]->x = c;
595: return(0);
596: }
598: /*@
599: PCMGSetR - Sets the vector space to be used to store the residual on a
600: particular level.
602: Logically Collective on PC and Vec
604: Input Parameters:
605: + pc - the multigrid context
606: . l - the level (0 is coarsest) this is to be used for
607: - c - the space
609: Level: advanced
611: Notes: If this is not provided PETSc will automatically generate one.
613: You do not need to keep a reference to this vector if you do
614: not need it PCDestroy() will properly free it.
616: .keywords: MG, multigrid, set, residual, level
617: @*/
618: PetscErrorCode PCMGSetR(PC pc,PetscInt l,Vec c)
619: {
621: PC_MG *mg = (PC_MG*)pc->data;
622: PC_MG_Levels **mglevels = mg->levels;
626: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
627: if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Need not set residual vector for coarse grid");
628: PetscObjectReference((PetscObject)c);
629: VecDestroy(&mglevels[l]->r);
631: mglevels[l]->r = c;
632: return(0);
633: }