Actual source code: mgfunc.c
petsc-3.7.5 2017-01-01
2: #include <petsc/private/pcmgimpl.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, can be NULL
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 PCMGSetRestriction()");
178: PCMGSetInterpolation(pc,l,mglevels[l]->restrct);
179: }
180: if (mat) *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: if (mat) *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: MatCreateVecs(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: Notes:
381: 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.
382: You can also modify smoother options by calling the various KSPSetXXX() options on this ksp. In addition you can call KSPGetPC(ksp,&pc)
383: and modify PC options for the smoother; for example PCSetType(pc,PCSOR); to use SOR smoothing.
385: Level: advanced
387: .keywords: MG, get, multigrid, level, smoother, pre-smoother, post-smoother
389: .seealso: PCMGGetSmootherUp(), PCMGGetSmootherDown()
390: @*/
391: PetscErrorCode PCMGGetSmoother(PC pc,PetscInt l,KSP *ksp)
392: {
393: PC_MG *mg = (PC_MG*)pc->data;
394: PC_MG_Levels **mglevels = mg->levels;
398: *ksp = mglevels[l]->smoothd;
399: return(0);
400: }
404: /*@
405: PCMGGetSmootherUp - Gets the KSP context to be used as smoother after
406: coarse grid correction (post-smoother).
408: Not Collective, KSP returned is parallel if PC is
410: Input Parameters:
411: + pc - the multigrid context
412: - l - the level (0 is coarsest) to supply
414: Ouput Parameters:
415: . ksp - the smoother
417: Level: advanced
419: Notes: calling this will result in a different pre and post smoother so you may need to
420: set options on the pre smoother also
422: .keywords: MG, multigrid, get, smoother, up, post-smoother, level
424: .seealso: PCMGGetSmootherUp(), PCMGGetSmootherDown()
425: @*/
426: PetscErrorCode PCMGGetSmootherUp(PC pc,PetscInt l,KSP *ksp)
427: {
428: PC_MG *mg = (PC_MG*)pc->data;
429: PC_MG_Levels **mglevels = mg->levels;
431: const char *prefix;
432: MPI_Comm comm;
436: /*
437: This is called only if user wants a different pre-smoother from post.
438: Thus we check if a different one has already been allocated,
439: if not we allocate it.
440: */
441: if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"There is no such thing as a up smoother on the coarse grid");
442: if (mglevels[l]->smoothu == mglevels[l]->smoothd) {
443: KSPType ksptype;
444: PCType pctype;
445: PC ipc;
446: PetscReal rtol,abstol,dtol;
447: PetscInt maxits;
448: KSPNormType normtype;
449: PetscObjectGetComm((PetscObject)mglevels[l]->smoothd,&comm);
450: KSPGetOptionsPrefix(mglevels[l]->smoothd,&prefix);
451: KSPGetTolerances(mglevels[l]->smoothd,&rtol,&abstol,&dtol,&maxits);
452: KSPGetType(mglevels[l]->smoothd,&ksptype);
453: KSPGetNormType(mglevels[l]->smoothd,&normtype);
454: KSPGetPC(mglevels[l]->smoothd,&ipc);
455: PCGetType(ipc,&pctype);
457: KSPCreate(comm,&mglevels[l]->smoothu);
458: KSPSetErrorIfNotConverged(mglevels[l]->smoothu,pc->erroriffailure);
459: PetscObjectIncrementTabLevel((PetscObject)mglevels[l]->smoothu,(PetscObject)pc,mglevels[0]->levels-l);
460: KSPSetOptionsPrefix(mglevels[l]->smoothu,prefix);
461: KSPSetTolerances(mglevels[l]->smoothu,rtol,abstol,dtol,maxits);
462: KSPSetType(mglevels[l]->smoothu,ksptype);
463: KSPSetNormType(mglevels[l]->smoothu,normtype);
464: KSPSetConvergenceTest(mglevels[l]->smoothu,KSPConvergedSkip,NULL,NULL);
465: KSPGetPC(mglevels[l]->smoothu,&ipc);
466: PCSetType(ipc,pctype);
467: PetscLogObjectParent((PetscObject)pc,(PetscObject)mglevels[l]->smoothu);
468: PetscObjectComposedDataSetInt((PetscObject) mglevels[l]->smoothu, PetscMGLevelId, mglevels[l]->level);
469: }
470: if (ksp) *ksp = mglevels[l]->smoothu;
471: return(0);
472: }
476: /*@
477: PCMGGetSmootherDown - Gets the KSP context to be used as smoother before
478: coarse grid correction (pre-smoother).
480: Not Collective, KSP returned is parallel if PC is
482: Input Parameters:
483: + pc - the multigrid context
484: - l - the level (0 is coarsest) to supply
486: Ouput Parameters:
487: . ksp - the smoother
489: Level: advanced
491: Notes: calling this will result in a different pre and post smoother so you may need to
492: set options on the post smoother also
494: .keywords: MG, multigrid, get, smoother, down, pre-smoother, level
496: .seealso: PCMGGetSmootherUp(), PCMGGetSmoother()
497: @*/
498: PetscErrorCode PCMGGetSmootherDown(PC pc,PetscInt l,KSP *ksp)
499: {
501: PC_MG *mg = (PC_MG*)pc->data;
502: PC_MG_Levels **mglevels = mg->levels;
506: /* make sure smoother up and down are different */
507: if (l) {
508: PCMGGetSmootherUp(pc,l,NULL);
509: }
510: *ksp = mglevels[l]->smoothd;
511: return(0);
512: }
516: /*@
517: PCMGSetCyclesOnLevel - Sets the number of cycles to run on this level.
519: Logically Collective on PC
521: Input Parameters:
522: + pc - the multigrid context
523: . l - the level (0 is coarsest) this is to be used for
524: - n - the number of cycles
526: Level: advanced
528: .keywords: MG, multigrid, set, cycles, V-cycle, W-cycle, level
530: .seealso: PCMGSetCycles()
531: @*/
532: PetscErrorCode PCMGSetCyclesOnLevel(PC pc,PetscInt l,PetscInt c)
533: {
534: PC_MG *mg = (PC_MG*)pc->data;
535: PC_MG_Levels **mglevels = mg->levels;
539: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
542: mglevels[l]->cycles = c;
543: return(0);
544: }
548: /*@
549: PCMGSetRhs - Sets the vector space to be used to store the right-hand side
550: on a particular level.
552: Logically Collective on PC and Vec
554: Input Parameters:
555: + pc - the multigrid context
556: . l - the level (0 is coarsest) this is to be used for
557: - c - the space
559: Level: advanced
561: Notes: If this is not provided PETSc will automatically generate one.
563: You do not need to keep a reference to this vector if you do
564: not need it PCDestroy() will properly free it.
566: .keywords: MG, multigrid, set, right-hand-side, rhs, level
568: .seealso: PCMGSetX(), PCMGSetR()
569: @*/
570: PetscErrorCode PCMGSetRhs(PC pc,PetscInt l,Vec c)
571: {
573: PC_MG *mg = (PC_MG*)pc->data;
574: PC_MG_Levels **mglevels = mg->levels;
578: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
579: if (l == mglevels[0]->levels-1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Do not set rhs for finest level");
580: PetscObjectReference((PetscObject)c);
581: VecDestroy(&mglevels[l]->b);
583: mglevels[l]->b = c;
584: return(0);
585: }
589: /*@
590: PCMGSetX - Sets the vector space to be used to store the solution on a
591: particular level.
593: Logically Collective on PC and Vec
595: Input Parameters:
596: + pc - the multigrid context
597: . l - the level (0 is coarsest) this is to be used for (do not supply the finest level)
598: - c - the space
600: Level: advanced
602: Notes: If this is not provided PETSc will automatically generate one.
604: You do not need to keep a reference to this vector if you do
605: not need it PCDestroy() will properly free it.
607: .keywords: MG, multigrid, set, solution, level
609: .seealso: PCMGSetRhs(), PCMGSetR()
610: @*/
611: PetscErrorCode PCMGSetX(PC pc,PetscInt l,Vec c)
612: {
614: PC_MG *mg = (PC_MG*)pc->data;
615: PC_MG_Levels **mglevels = mg->levels;
619: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
620: if (l == mglevels[0]->levels-1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Do not set x for finest level");
621: PetscObjectReference((PetscObject)c);
622: VecDestroy(&mglevels[l]->x);
624: mglevels[l]->x = c;
625: return(0);
626: }
630: /*@
631: PCMGSetR - Sets the vector space to be used to store the residual on a
632: particular level.
634: Logically Collective on PC and Vec
636: Input Parameters:
637: + pc - the multigrid context
638: . l - the level (0 is coarsest) this is to be used for
639: - c - the space
641: Level: advanced
643: Notes: If this is not provided PETSc will automatically generate one.
645: You do not need to keep a reference to this vector if you do
646: not need it PCDestroy() will properly free it.
648: .keywords: MG, multigrid, set, residual, level
649: @*/
650: PetscErrorCode PCMGSetR(PC pc,PetscInt l,Vec c)
651: {
653: PC_MG *mg = (PC_MG*)pc->data;
654: PC_MG_Levels **mglevels = mg->levels;
658: if (!mglevels) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONGSTATE,"Must set MG levels before calling");
659: if (!l) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Need not set residual vector for coarse grid");
660: PetscObjectReference((PetscObject)c);
661: VecDestroy(&mglevels[l]->r);
663: mglevels[l]->r = c;
664: return(0);
665: }