Actual source code: mgfunc.c
2: #include <petsc/private/pcmgimpl.h>
4: /*@C
5: PCMGResidualDefault - Default routine to calculate the residual.
7: Collective
9: Input Parameters:
10: + mat - the matrix
11: . b - the right-hand-side
12: - x - the approximate solution
14: Output Parameter:
15: . r - location to store the residual
17: Level: developer
19: .seealso: `PCMG`, `PCMGSetResidual()`, `PCMGSetMatResidual()`
20: @*/
21: PetscErrorCode PCMGResidualDefault(Mat mat, Vec b, Vec x, Vec r)
22: {
23: MatResidual(mat, b, x, r);
24: return 0;
25: }
27: /*@C
28: PCMGResidualTransposeDefault - Default routine to calculate the residual of the transposed linear system
30: Collective
32: Input Parameters:
33: + mat - the matrix
34: . b - the right-hand-side
35: - x - the approximate solution
37: Output Parameter:
38: . r - location to store the residual
40: Level: developer
42: .seealso: `PCMG`, `PCMGSetResidualTranspose()`, `PCMGMatResidualTransposeDefault()`
43: @*/
44: PetscErrorCode PCMGResidualTransposeDefault(Mat mat, Vec b, Vec x, Vec r)
45: {
46: MatMultTranspose(mat, x, r);
47: VecAYPX(r, -1.0, b);
48: return 0;
49: }
51: /*@C
52: PCMGMatResidualDefault - Default routine to calculate the residual.
54: Collective
56: Input Parameters:
57: + mat - the matrix
58: . b - the right-hand-side
59: - x - the approximate solution
61: Output Parameter:
62: . r - location to store the residual
64: Level: developer
66: .seealso: `PCMG`, `PCMGSetMatResidual()`, `PCMGResidualDefault()`
67: @*/
68: PetscErrorCode PCMGMatResidualDefault(Mat mat, Mat b, Mat x, Mat r)
69: {
70: MatMatMult(mat, x, MAT_REUSE_MATRIX, PETSC_DEFAULT, &r);
71: MatAYPX(r, -1.0, b, UNKNOWN_NONZERO_PATTERN);
72: return 0;
73: }
75: /*@C
76: PCMGMatResidualTransposeDefault - Default routine to calculate the residual of the transposed linear system
78: Collective
80: Input Parameters:
81: + mat - the matrix
82: . b - the right-hand-side
83: - x - the approximate solution
85: Output Parameter:
86: . r - location to store the residual
88: Level: developer
90: .seealso: `PCMG`, `PCMGSetMatResidualTranspose()`
91: @*/
92: PetscErrorCode PCMGMatResidualTransposeDefault(Mat mat, Mat b, Mat x, Mat r)
93: {
94: MatTransposeMatMult(mat, x, MAT_REUSE_MATRIX, PETSC_DEFAULT, &r);
95: MatAYPX(r, -1.0, b, UNKNOWN_NONZERO_PATTERN);
96: return 0;
97: }
98: /*@
99: PCMGGetCoarseSolve - Gets the solver context to be used on the coarse grid.
101: Not Collective
103: Input Parameter:
104: . pc - the multigrid context
106: Output Parameter:
107: . ksp - the coarse grid solver context
109: Level: advanced
111: .seealso: `PCMG`, `PCMGGetSmootherUp()`, `PCMGGetSmootherDown()`, `PCMGGetSmoother()`
112: @*/
113: PetscErrorCode PCMGGetCoarseSolve(PC pc, KSP *ksp)
114: {
115: PC_MG *mg = (PC_MG *)pc->data;
116: PC_MG_Levels **mglevels = mg->levels;
119: *ksp = mglevels[0]->smoothd;
120: return 0;
121: }
123: /*@C
124: PCMGSetResidual - Sets the function to be used to calculate the residual on the lth level.
126: Logically Collective
128: Input Parameters:
129: + pc - the multigrid context
130: . l - the level (0 is coarsest) to supply
131: . residual - function used to form residual, if none is provided the previously provide one is used, if no
132: previous one were provided then a default is used
133: - mat - matrix associated with residual
135: Level: advanced
137: .seealso: `PCMG`, `PCMGResidualDefault()`
138: @*/
139: PetscErrorCode PCMGSetResidual(PC pc, PetscInt l, PetscErrorCode (*residual)(Mat, Vec, Vec, Vec), Mat mat)
140: {
141: PC_MG *mg = (PC_MG *)pc->data;
142: PC_MG_Levels **mglevels = mg->levels;
146: if (residual) mglevels[l]->residual = residual;
147: if (!mglevels[l]->residual) mglevels[l]->residual = PCMGResidualDefault;
148: mglevels[l]->matresidual = PCMGMatResidualDefault;
149: if (mat) PetscObjectReference((PetscObject)mat);
150: MatDestroy(&mglevels[l]->A);
151: mglevels[l]->A = mat;
152: return 0;
153: }
155: /*@C
156: PCMGSetResidualTranspose - Sets the function to be used to calculate the residual of the transposed linear system
157: on the lth level.
159: Logically Collective
161: Input Parameters:
162: + pc - the multigrid context
163: . l - the level (0 is coarsest) to supply
164: . residualt - function used to form transpose of residual, if none is provided the previously provide one is used, if no
165: previous one were provided then a default is used
166: - mat - matrix associated with residual
168: Level: advanced
170: .seealso: `PCMG`, `PCMGResidualTransposeDefault()`
171: @*/
172: PetscErrorCode PCMGSetResidualTranspose(PC pc, PetscInt l, PetscErrorCode (*residualt)(Mat, Vec, Vec, Vec), Mat mat)
173: {
174: PC_MG *mg = (PC_MG *)pc->data;
175: PC_MG_Levels **mglevels = mg->levels;
179: if (residualt) mglevels[l]->residualtranspose = residualt;
180: if (!mglevels[l]->residualtranspose) mglevels[l]->residualtranspose = PCMGResidualTransposeDefault;
181: mglevels[l]->matresidualtranspose = PCMGMatResidualTransposeDefault;
182: if (mat) PetscObjectReference((PetscObject)mat);
183: MatDestroy(&mglevels[l]->A);
184: mglevels[l]->A = mat;
185: return 0;
186: }
188: /*@
189: PCMGSetInterpolation - Sets the function to be used to calculate the
190: interpolation from l-1 to the lth level
192: Logically Collective
194: Input Parameters:
195: + pc - the multigrid context
196: . mat - the interpolation operator
197: - l - the level (0 is coarsest) to supply [do not supply 0]
199: Level: advanced
201: Notes:
202: Usually this is the same matrix used also to set the restriction
203: for the same level.
205: One can pass in the interpolation matrix or its transpose; PETSc figures
206: out from the matrix size which one it is.
208: .seealso: `PCMG`, `PCMGSetRestriction()`
209: @*/
210: PetscErrorCode PCMGSetInterpolation(PC pc, PetscInt l, Mat mat)
211: {
212: PC_MG *mg = (PC_MG *)pc->data;
213: PC_MG_Levels **mglevels = mg->levels;
218: PetscObjectReference((PetscObject)mat);
219: MatDestroy(&mglevels[l]->interpolate);
221: mglevels[l]->interpolate = mat;
222: return 0;
223: }
225: /*@
226: PCMGSetOperators - Sets operator and preconditioning matrix for lth level
228: Logically Collective
230: Input Parameters:
231: + pc - the multigrid context
232: . Amat - the operator
233: . pmat - the preconditioning operator
234: - l - the level (0 is the coarsest) to supply
236: Level: advanced
238: .keywords: multigrid, set, interpolate, level
240: .seealso: `PCMG`, `PCMGSetGalerkin()`, `PCMGSetRestriction()`, `PCMGSetInterpolation()`
241: @*/
242: PetscErrorCode PCMGSetOperators(PC pc, PetscInt l, Mat Amat, Mat Pmat)
243: {
244: PC_MG *mg = (PC_MG *)pc->data;
245: PC_MG_Levels **mglevels = mg->levels;
251: KSPSetOperators(mglevels[l]->smoothd, Amat, Pmat);
252: return 0;
253: }
255: /*@
256: PCMGGetInterpolation - Gets the function to be used to calculate the
257: interpolation from l-1 to the lth level
259: Logically Collective
261: Input Parameters:
262: + pc - the multigrid context
263: - l - the level (0 is coarsest) to supply [Do not supply 0]
265: Output Parameter:
266: . mat - the interpolation matrix, can be NULL
268: Level: advanced
270: .seealso: `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetRScale()`
271: @*/
272: PetscErrorCode PCMGGetInterpolation(PC pc, PetscInt l, Mat *mat)
273: {
274: PC_MG *mg = (PC_MG *)pc->data;
275: PC_MG_Levels **mglevels = mg->levels;
281: if (!mglevels[l]->interpolate && mglevels[l]->restrct) PCMGSetInterpolation(pc, l, mglevels[l]->restrct);
282: if (mat) *mat = mglevels[l]->interpolate;
283: return 0;
284: }
286: /*@
287: PCMGSetRestriction - Sets the function to be used to restrict dual vectors
288: from level l to l-1.
290: Logically Collective
292: Input Parameters:
293: + pc - the multigrid context
294: . l - the level (0 is coarsest) to supply [Do not supply 0]
295: - mat - the restriction matrix
297: Level: advanced
299: Notes:
300: Usually this is the same matrix used also to set the interpolation
301: for the same level.
303: One can pass in the interpolation matrix or its transpose; PETSc figures
304: out from the matrix size which one it is.
306: If you do not set this, the transpose of the `Mat` set with `PCMGSetInterpolation()`
307: is used.
309: .seealso: `PCMG`, `PCMGSetInterpolation()`
310: @*/
311: PetscErrorCode PCMGSetRestriction(PC pc, PetscInt l, Mat mat)
312: {
313: PC_MG *mg = (PC_MG *)pc->data;
314: PC_MG_Levels **mglevels = mg->levels;
320: PetscObjectReference((PetscObject)mat);
321: MatDestroy(&mglevels[l]->restrct);
323: mglevels[l]->restrct = mat;
324: return 0;
325: }
327: /*@
328: PCMGGetRestriction - Gets the function to be used to restrict dual vectors
329: from level l to l-1.
331: Logically Collective
333: Input Parameters:
334: + pc - the multigrid context
335: - l - the level (0 is coarsest) to supply [Do not supply 0]
337: Output Parameter:
338: . mat - the restriction matrix
340: Level: advanced
342: .seealso: `PCMG`, `PCMGGetInterpolation()`, `PCMGSetRestriction()`, `PCMGGetRScale()`, `PCMGGetInjection()`
343: @*/
344: PetscErrorCode PCMGGetRestriction(PC pc, PetscInt l, Mat *mat)
345: {
346: PC_MG *mg = (PC_MG *)pc->data;
347: PC_MG_Levels **mglevels = mg->levels;
353: if (!mglevels[l]->restrct && mglevels[l]->interpolate) PCMGSetRestriction(pc, l, mglevels[l]->interpolate);
354: if (mat) *mat = mglevels[l]->restrct;
355: return 0;
356: }
358: /*@
359: PCMGSetRScale - Sets the pointwise scaling for the restriction operator from level l to l-1.
361: Logically Collective
363: Input Parameters:
364: + pc - the multigrid context
365: - l - the level (0 is coarsest) to supply [Do not supply 0]
366: . rscale - the scaling
368: Level: advanced
370: Note:
371: 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.
372: It is preferable to use `PCMGSetInjection()` to control moving primal vectors.
374: .seealso: `PCMG`, `PCMGSetInterpolation()`, `PCMGSetRestriction()`, `PCMGGetRScale()`, `PCMGSetInjection()`
375: @*/
376: PetscErrorCode PCMGSetRScale(PC pc, PetscInt l, Vec rscale)
377: {
378: PC_MG *mg = (PC_MG *)pc->data;
379: PC_MG_Levels **mglevels = mg->levels;
384: PetscObjectReference((PetscObject)rscale);
385: VecDestroy(&mglevels[l]->rscale);
387: mglevels[l]->rscale = rscale;
388: return 0;
389: }
391: /*@
392: PCMGGetRScale - Gets the pointwise scaling for the restriction operator from level l to l-1.
394: Collective
396: Input Parameters:
397: + pc - the multigrid context
398: . rscale - the scaling
399: - l - the level (0 is coarsest) to supply [Do not supply 0]
401: Level: advanced
403: Note:
404: 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.
405: It is preferable to use `PCMGGetInjection()` to control moving primal vectors.
407: .seealso: `PCMG`, `PCMGSetInterpolation()`, `PCMGGetRestriction()`, `PCMGGetInjection()`
408: @*/
409: PetscErrorCode PCMGGetRScale(PC pc, PetscInt l, Vec *rscale)
410: {
411: PC_MG *mg = (PC_MG *)pc->data;
412: PC_MG_Levels **mglevels = mg->levels;
417: if (!mglevels[l]->rscale) {
418: Mat R;
419: Vec X, Y, coarse, fine;
420: PetscInt M, N;
421: PCMGGetRestriction(pc, l, &R);
422: MatCreateVecs(R, &X, &Y);
423: MatGetSize(R, &M, &N);
424: if (M < N) {
425: fine = X;
426: coarse = Y;
427: } else if (N < M) {
428: fine = Y;
429: coarse = X;
430: } else SETERRQ(PetscObjectComm((PetscObject)R), PETSC_ERR_SUP, "Restriction matrix is square, cannot determine which Vec is coarser");
431: VecSet(fine, 1.);
432: MatRestrict(R, fine, coarse);
433: VecDestroy(&fine);
434: VecReciprocal(coarse);
435: mglevels[l]->rscale = coarse;
436: }
437: *rscale = mglevels[l]->rscale;
438: return 0;
439: }
441: /*@
442: PCMGSetInjection - Sets the function to be used to inject primal vectors
443: from level l to l-1.
445: Logically Collective
447: Input Parameters:
448: + pc - the multigrid context
449: . l - the level (0 is coarsest) to supply [Do not supply 0]
450: - mat - the injection matrix
452: Level: advanced
454: .seealso: `PCMG`, `PCMGSetRestriction()`
455: @*/
456: PetscErrorCode PCMGSetInjection(PC pc, PetscInt l, Mat mat)
457: {
458: PC_MG *mg = (PC_MG *)pc->data;
459: PC_MG_Levels **mglevels = mg->levels;
465: PetscObjectReference((PetscObject)mat);
466: MatDestroy(&mglevels[l]->inject);
468: mglevels[l]->inject = mat;
469: return 0;
470: }
472: /*@
473: PCMGGetInjection - Gets the function to be used to inject primal vectors
474: from level l to l-1.
476: Logically Collective
478: Input Parameters:
479: + pc - the multigrid context
480: - l - the level (0 is coarsest) to supply [Do not supply 0]
482: Output Parameter:
483: . mat - the restriction matrix (may be NULL if no injection is available).
485: Level: advanced
487: .seealso: `PCMG`, `PCMGSetInjection()`, `PCMGetGetRestriction()`
488: @*/
489: PetscErrorCode PCMGGetInjection(PC pc, PetscInt l, Mat *mat)
490: {
491: PC_MG *mg = (PC_MG *)pc->data;
492: PC_MG_Levels **mglevels = mg->levels;
498: if (mat) *mat = mglevels[l]->inject;
499: return 0;
500: }
502: /*@
503: PCMGGetSmoother - Gets the `KSP` context to be used as smoother for
504: both pre- and post-smoothing. Call both `PCMGGetSmootherUp()` and
505: `PCMGGetSmootherDown()` to use different functions for pre- and
506: post-smoothing.
508: Not Collective, ksp returned is parallel if pc is
510: Input Parameters:
511: + pc - the multigrid context
512: - l - the level (0 is coarsest) to supply
514: Output Parameter:
515: . ksp - the smoother
517: Note:
518: Once you have called this routine, you can call `KSPSetOperators()` on the resulting ksp to provide the operators for the smoother for this level.
519: You can also modify smoother options by calling the various KSPSetXXX() options on this ksp. In addition you can call `KSPGetPC`(ksp,&pc)
520: and modify PC options for the smoother; for example `PCSetType`(pc,`PCSOR`); to use SOR smoothing.
522: Level: advanced
524: .seealso: PCMG`, ``PCMGGetSmootherUp()`, `PCMGGetSmootherDown()`, `PCMGGetCoarseSolve()`
525: @*/
526: PetscErrorCode PCMGGetSmoother(PC pc, PetscInt l, KSP *ksp)
527: {
528: PC_MG *mg = (PC_MG *)pc->data;
529: PC_MG_Levels **mglevels = mg->levels;
532: *ksp = mglevels[l]->smoothd;
533: return 0;
534: }
536: /*@
537: PCMGGetSmootherUp - Gets the KSP context to be used as smoother after
538: coarse grid correction (post-smoother).
540: Not Collective, ksp returned is parallel if pc is
542: Input Parameters:
543: + pc - the multigrid context
544: - l - the level (0 is coarsest) to supply
546: Output Parameter:
547: . ksp - the smoother
549: Level: advanced
551: Note:
552: Calling this will result in a different pre and post smoother so you may need to set options on the pre smoother also
554: .seealso: `PCMG`, `PCMGGetSmootherUp()`, `PCMGGetSmootherDown()`
555: @*/
556: PetscErrorCode PCMGGetSmootherUp(PC pc, PetscInt l, KSP *ksp)
557: {
558: PC_MG *mg = (PC_MG *)pc->data;
559: PC_MG_Levels **mglevels = mg->levels;
560: const char *prefix;
561: MPI_Comm comm;
564: /*
565: This is called only if user wants a different pre-smoother from post.
566: Thus we check if a different one has already been allocated,
567: if not we allocate it.
568: */
570: if (mglevels[l]->smoothu == mglevels[l]->smoothd) {
571: KSPType ksptype;
572: PCType pctype;
573: PC ipc;
574: PetscReal rtol, abstol, dtol;
575: PetscInt maxits;
576: KSPNormType normtype;
577: PetscObjectGetComm((PetscObject)mglevels[l]->smoothd, &comm);
578: KSPGetOptionsPrefix(mglevels[l]->smoothd, &prefix);
579: KSPGetTolerances(mglevels[l]->smoothd, &rtol, &abstol, &dtol, &maxits);
580: KSPGetType(mglevels[l]->smoothd, &ksptype);
581: KSPGetNormType(mglevels[l]->smoothd, &normtype);
582: KSPGetPC(mglevels[l]->smoothd, &ipc);
583: PCGetType(ipc, &pctype);
585: KSPCreate(comm, &mglevels[l]->smoothu);
586: KSPSetErrorIfNotConverged(mglevels[l]->smoothu, pc->erroriffailure);
587: PetscObjectIncrementTabLevel((PetscObject)mglevels[l]->smoothu, (PetscObject)pc, mglevels[0]->levels - l);
588: KSPSetOptionsPrefix(mglevels[l]->smoothu, prefix);
589: KSPSetTolerances(mglevels[l]->smoothu, rtol, abstol, dtol, maxits);
590: KSPSetType(mglevels[l]->smoothu, ksptype);
591: KSPSetNormType(mglevels[l]->smoothu, normtype);
592: KSPSetConvergenceTest(mglevels[l]->smoothu, KSPConvergedSkip, NULL, NULL);
593: KSPGetPC(mglevels[l]->smoothu, &ipc);
594: PCSetType(ipc, pctype);
595: PetscObjectComposedDataSetInt((PetscObject)mglevels[l]->smoothu, PetscMGLevelId, mglevels[l]->level);
596: }
597: if (ksp) *ksp = mglevels[l]->smoothu;
598: return 0;
599: }
601: /*@
602: PCMGGetSmootherDown - Gets the `KSP` context to be used as smoother before
603: coarse grid correction (pre-smoother).
605: Not Collective, ksp returned is parallel if pc is
607: Input Parameters:
608: + pc - the multigrid context
609: - l - the level (0 is coarsest) to supply
611: Output Parameter:
612: . ksp - the smoother
614: Level: advanced
616: Note:
617: Calling this will result in a different pre and post smoother so you may need to
618: set options on the post smoother also
620: .seealso: `PCMG`, `PCMGGetSmootherUp()`, `PCMGGetSmoother()`
621: @*/
622: PetscErrorCode PCMGGetSmootherDown(PC pc, PetscInt l, KSP *ksp)
623: {
624: PC_MG *mg = (PC_MG *)pc->data;
625: PC_MG_Levels **mglevels = mg->levels;
628: /* make sure smoother up and down are different */
629: if (l) PCMGGetSmootherUp(pc, l, NULL);
630: *ksp = mglevels[l]->smoothd;
631: return 0;
632: }
634: /*@
635: PCMGSetCycleTypeOnLevel - Sets the type of cycle (aka cycle index) to run on the specified level.
637: Logically Collective
639: Input Parameters:
640: + pc - the multigrid context
641: . l - the level (0 is coarsest)
642: - c - either `PC_MG_CYCLE_V` or `PC_MG_CYCLE_W`
644: Level: advanced
646: .seealso: `PCMG`, PCMGCycleType`, `PCMGSetCycleType()`
647: @*/
648: PetscErrorCode PCMGSetCycleTypeOnLevel(PC pc, PetscInt l, PCMGCycleType c)
649: {
650: PC_MG *mg = (PC_MG *)pc->data;
651: PC_MG_Levels **mglevels = mg->levels;
657: mglevels[l]->cycles = c;
658: return 0;
659: }
661: /*@
662: PCMGSetRhs - Sets the vector to be used to store the right-hand side on a particular level.
664: Logically Collective
666: Input Parameters:
667: + pc - the multigrid context
668: . l - the level (0 is coarsest) this is to be used for
669: - c - the Vec
671: Level: advanced
673: Note:
674: If this is not provided PETSc will automatically generate one. You do not need to keep a reference to this vector if you do not need it. `PCDestroy()` will properly free it.
676: .seealso: `PCMG`, `PCMGSetX()`, `PCMGSetR()`
677: @*/
678: PetscErrorCode PCMGSetRhs(PC pc, PetscInt l, Vec c)
679: {
680: PC_MG *mg = (PC_MG *)pc->data;
681: PC_MG_Levels **mglevels = mg->levels;
686: PetscObjectReference((PetscObject)c);
687: VecDestroy(&mglevels[l]->b);
689: mglevels[l]->b = c;
690: return 0;
691: }
693: /*@
694: PCMGSetX - Sets the vector to be used to store the solution on a particular level.
696: Logically Collective
698: Input Parameters:
699: + pc - the multigrid context
700: . l - the level (0 is coarsest) this is to be used for (do not supply the finest level)
701: - c - the Vec
703: Level: advanced
705: Note:
706: If this is not provided PETSc will automatically generate one. You do not need to keep a reference to this vector if you do not need it. `PCDestroy()` will properly free it.
708: .seealso: `PCMG`, `PCMGSetRhs()`, `PCMGSetR()`
709: @*/
710: PetscErrorCode PCMGSetX(PC pc, PetscInt l, Vec c)
711: {
712: PC_MG *mg = (PC_MG *)pc->data;
713: PC_MG_Levels **mglevels = mg->levels;
718: PetscObjectReference((PetscObject)c);
719: VecDestroy(&mglevels[l]->x);
721: mglevels[l]->x = c;
722: return 0;
723: }
725: /*@
726: PCMGSetR - Sets the vector to be used to store the residual on a particular level.
728: Logically Collective
730: Input Parameters:
731: + pc - the multigrid context
732: . l - the level (0 is coarsest) this is to be used for
733: - c - the Vec
735: Level: advanced
737: Note:
738: If this is not provided PETSc will automatically generate one. You do not need to keep a reference to this vector if you do not need it. `PCDestroy()` will properly free it.
740: .seealso: `PCMG`, `PCMGSetRhs()`, `PCMGSetX()`
741: @*/
742: PetscErrorCode PCMGSetR(PC pc, PetscInt l, Vec c)
743: {
744: PC_MG *mg = (PC_MG *)pc->data;
745: PC_MG_Levels **mglevels = mg->levels;
750: PetscObjectReference((PetscObject)c);
751: VecDestroy(&mglevels[l]->r);
753: mglevels[l]->r = c;
754: return 0;
755: }