Actual source code: mgfunc.c
1: #include <petsc/private/pcmgimpl.h>
2: #include <petsc/private/matimpl.h>
4: /*@
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: [](ch_ksp), `PCMG`, `PCMGSetResidual()`, `PCMGSetMatResidual()`
20: @*/
21: PetscErrorCode PCMGResidualDefault(Mat mat, Vec b, Vec x, Vec r)
22: {
23: PetscFunctionBegin;
24: PetscCall(MatResidual(mat, b, x, r));
25: PetscFunctionReturn(PETSC_SUCCESS);
26: }
28: /*@
29: PCMGResidualTransposeDefault - Default routine to calculate the residual of the transposed linear system
31: Collective
33: Input Parameters:
34: + mat - the matrix
35: . b - the right-hand side
36: - x - the approximate solution
38: Output Parameter:
39: . r - location to store the residual
41: Level: developer
43: .seealso: [](ch_ksp), `PCMG`, `PCMGSetResidualTranspose()`, `PCMGMatResidualTransposeDefault()`
44: @*/
45: PetscErrorCode PCMGResidualTransposeDefault(Mat mat, Vec b, Vec x, Vec r)
46: {
47: PetscFunctionBegin;
48: PetscCall(MatMultTranspose(mat, x, r));
49: PetscCall(VecAYPX(r, -1.0, b));
50: PetscFunctionReturn(PETSC_SUCCESS);
51: }
53: /*@
54: PCMGMatResidualDefault - Default routine to calculate the residual.
56: Collective
58: Input Parameters:
59: + mat - the matrix
60: . b - the right-hand side
61: - x - the approximate solution
63: Output Parameter:
64: . r - location to store the residual
66: Level: developer
68: Note:
69: When `r` is dense, the `MatProduct` computing it stays attached to `r`, so that later calls with the same `mat` and nonzero pattern only run
70: its numeric phase; it keeps references to `mat` and `x` until `MatProductClear()` is called on `r` or `r` is destroyed.
72: .seealso: [](ch_ksp), `PCMG`, `PCMGSetMatResidual()`, `PCMGResidualDefault()`
73: @*/
74: PetscErrorCode PCMGMatResidualDefault(Mat mat, Mat b, Mat x, Mat r)
75: {
76: PetscFunctionBegin;
77: PetscCall(MatProductComputeWithMat_Private(mat, x, MATPRODUCT_AB, r));
78: PetscCall(MatAYPX(r, -1.0, b, UNKNOWN_NONZERO_PATTERN));
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: /*@
83: PCMGMatResidualTransposeDefault - Default routine to calculate the residual of the transposed linear system
85: Collective
87: Input Parameters:
88: + mat - the matrix
89: . b - the right-hand side
90: - x - the approximate solution
92: Output Parameter:
93: . r - location to store the residual
95: Level: developer
97: Note:
98: When `r` is dense, the `MatProduct` computing it stays attached to `r`, so that later calls with the same `mat` and nonzero pattern only run
99: its numeric phase; it keeps references to `mat` and `x` until `MatProductClear()` is called on `r` or `r` is destroyed.
101: .seealso: [](ch_ksp), `PCMG`, `PCMGSetMatResidualTranspose()`
102: @*/
103: PetscErrorCode PCMGMatResidualTransposeDefault(Mat mat, Mat b, Mat x, Mat r)
104: {
105: PetscFunctionBegin;
106: PetscCall(MatProductComputeWithMat_Private(mat, x, MATPRODUCT_AtB, r));
107: PetscCall(MatAYPX(r, -1.0, b, UNKNOWN_NONZERO_PATTERN));
108: PetscFunctionReturn(PETSC_SUCCESS);
109: }
110: /*@
111: PCMGGetCoarseSolve - Gets the solver context to be used on the coarse grid.
113: Not Collective
115: Input Parameter:
116: . pc - the multigrid context
118: Output Parameter:
119: . ksp - the coarse grid solver context
121: Level: advanced
123: .seealso: [](ch_ksp), `PCMG`, `PCMGGetSmootherUp()`, `PCMGGetSmootherDown()`, `PCMGGetSmoother()`
124: @*/
125: PetscErrorCode PCMGGetCoarseSolve(PC pc, KSP *ksp)
126: {
127: PC_MG *mg = (PC_MG *)pc->data;
128: PC_MG_Levels **mglevels = mg->levels;
130: PetscFunctionBegin;
132: *ksp = mglevels[0]->smoothd;
133: PetscFunctionReturn(PETSC_SUCCESS);
134: }
136: /*@
137: PCMGSetResidual - Sets the function to be used to calculate the residual on the lth level.
139: Logically Collective
141: Input Parameters:
142: + pc - the multigrid context
143: . l - the level (0 is coarsest) to supply
144: . residual - function used to form residual, if none is provided the previously provide one is used, if no
145: previous one were provided then a default is used
146: - mat - matrix associated with residual
148: Level: advanced
150: .seealso: [](ch_ksp), `PCMG`, `PCMGResidualDefault()`
151: @*/
152: PetscErrorCode PCMGSetResidual(PC pc, PetscInt l, PetscErrorCode (*residual)(Mat, Vec, Vec, Vec), Mat mat)
153: {
154: PC_MG *mg = (PC_MG *)pc->data;
155: PC_MG_Levels **mglevels = mg->levels;
157: PetscFunctionBegin;
159: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
160: if (residual) mglevels[l]->residual = residual;
161: if (!mglevels[l]->residual) mglevels[l]->residual = PCMGResidualDefault;
162: mglevels[l]->matresidual = PCMGMatResidualDefault;
163: PetscCall(PetscObjectReference((PetscObject)mat));
164: PetscCall(MatDestroy(&mglevels[l]->A));
165: mglevels[l]->A = mat;
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: /*@
170: PCMGSetResidualTranspose - Sets the function to be used to calculate the residual of the transposed linear system
171: on the lth level.
173: Logically Collective
175: Input Parameters:
176: + pc - the multigrid context
177: . l - the level (0 is coarsest) to supply
178: . residualt - function used to form transpose of residual, if none is provided the previously provide one is used, if no
179: previous one were provided then a default is used
180: - mat - matrix associated with residual
182: Level: advanced
184: .seealso: [](ch_ksp), `PCMG`, `PCMGResidualTransposeDefault()`
185: @*/
186: PetscErrorCode PCMGSetResidualTranspose(PC pc, PetscInt l, PetscErrorCode (*residualt)(Mat, Vec, Vec, Vec), Mat mat)
187: {
188: PC_MG *mg = (PC_MG *)pc->data;
189: PC_MG_Levels **mglevels = mg->levels;
191: PetscFunctionBegin;
193: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
194: if (residualt) mglevels[l]->residualtranspose = residualt;
195: if (!mglevels[l]->residualtranspose) mglevels[l]->residualtranspose = PCMGResidualTransposeDefault;
196: mglevels[l]->matresidualtranspose = PCMGMatResidualTransposeDefault;
197: PetscCall(PetscObjectReference((PetscObject)mat));
198: PetscCall(MatDestroy(&mglevels[l]->A));
199: mglevels[l]->A = mat;
200: PetscFunctionReturn(PETSC_SUCCESS);
201: }
203: /*@
204: PCMGSetInterpolation - Sets the function to be used to calculate the
205: interpolation from l-1 to the lth level
207: Logically Collective
209: Input Parameters:
210: + pc - the multigrid context
211: . mat - the interpolation operator
212: - l - the level (0 is coarsest) to supply [do not supply 0]
214: Level: advanced
216: Notes:
217: Usually this is the same matrix used also to set the restriction
218: for the same level.
220: One can pass in the interpolation matrix or its transpose; PETSc figures
221: out from the matrix size which one it is.
223: .seealso: [](ch_ksp), `PCMG`, `PCMGSetRestriction()`
224: @*/
225: PetscErrorCode PCMGSetInterpolation(PC pc, PetscInt l, Mat mat)
226: {
227: PC_MG *mg = (PC_MG *)pc->data;
228: PC_MG_Levels **mglevels = mg->levels;
230: PetscFunctionBegin;
232: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
233: PetscCheck(l, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Do not set interpolation routine for coarsest level");
234: if (mat != mglevels[l]->interpolate) PetscCall(PCMGResetBlocks_Private(pc)); /* the kept block products were set up with the old matrix */
235: PetscCall(PetscObjectReference((PetscObject)mat));
236: PetscCall(MatDestroy(&mglevels[l]->interpolate));
238: mglevels[l]->interpolate = mat;
239: PetscFunctionReturn(PETSC_SUCCESS);
240: }
242: /*@
243: PCMGSetOperators - Sets operator and matrix from which to construct a preconditioner for lth level
245: Logically Collective
247: Input Parameters:
248: + pc - the multigrid context
249: . Amat - the operator
250: . Pmat - the preconditioning operator
251: - l - the level (0 is the coarsest) to supply
253: Level: advanced
255: .seealso: [](ch_ksp), `PCMG`, `PCMGSetGalerkin()`, `PCMGSetRestriction()`, `PCMGSetInterpolation()`
256: @*/
257: PetscErrorCode PCMGSetOperators(PC pc, PetscInt l, Mat Amat, Mat Pmat)
258: {
259: PC_MG *mg = (PC_MG *)pc->data;
260: PC_MG_Levels **mglevels = mg->levels;
262: PetscFunctionBegin;
266: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
267: PetscCall(KSPSetOperators(mglevels[l]->smoothd, Amat, Pmat));
268: PetscFunctionReturn(PETSC_SUCCESS);
269: }
271: /*@
272: PCMGGetInterpolation - Gets the function to be used to calculate the
273: interpolation from l-1 to the lth level
275: Logically Collective
277: Input Parameters:
278: + pc - the multigrid context
279: - l - the level (0 is coarsest) to supply [Do not supply 0]
281: Output Parameter:
282: . mat - the interpolation matrix, can be `NULL`
284: Level: advanced
286: .seealso: [](ch_ksp), `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetRScale()`
287: @*/
288: PetscErrorCode PCMGGetInterpolation(PC pc, PetscInt l, Mat *mat)
289: {
290: PC_MG *mg = (PC_MG *)pc->data;
291: PC_MG_Levels **mglevels = mg->levels;
293: PetscFunctionBegin;
295: if (mat) PetscAssertPointer(mat, 3);
296: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
297: PetscCheck(l > 0 && l < mg->nlevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Level %" PetscInt_FMT " must be in range {1,...,%" PetscInt_FMT "}", l, mg->nlevels - 1);
298: if (!mglevels[l]->interpolate && mglevels[l]->restrct) PetscCall(PCMGSetInterpolation(pc, l, mglevels[l]->restrct));
299: if (mat) *mat = mglevels[l]->interpolate;
300: PetscFunctionReturn(PETSC_SUCCESS);
301: }
303: /*@
304: PCMGSetRestriction - Sets the function to be used to restrict dual vectors
305: from level l to l-1.
307: Logically Collective
309: Input Parameters:
310: + pc - the multigrid context
311: . l - the level (0 is coarsest) to supply [Do not supply 0]
312: - mat - the restriction matrix
314: Level: advanced
316: Notes:
317: Usually this is the same matrix used also to set the interpolation
318: for the same level.
320: One can pass in the interpolation matrix or its transpose; PETSc figures
321: out from the matrix size which one it is.
323: If you do not set this, the transpose of the `Mat` set with `PCMGSetInterpolation()`
324: is used.
326: .seealso: [](ch_ksp), `PCMG`, `PCMGSetInterpolation()`
327: @*/
328: PetscErrorCode PCMGSetRestriction(PC pc, PetscInt l, Mat mat)
329: {
330: PC_MG *mg = (PC_MG *)pc->data;
331: PC_MG_Levels **mglevels = mg->levels;
333: PetscFunctionBegin;
336: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
337: PetscCheck(l, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Do not set restriction routine for coarsest level");
338: if (mat != mglevels[l]->restrct) PetscCall(PCMGResetBlocks_Private(pc)); /* the kept block products were set up with the old matrix */
339: PetscCall(PetscObjectReference((PetscObject)mat));
340: PetscCall(MatDestroy(&mglevels[l]->restrct));
342: mglevels[l]->restrct = mat;
343: PetscFunctionReturn(PETSC_SUCCESS);
344: }
346: /*@
347: PCMGGetRestriction - Gets the function to be used to restrict dual (i.e. residual) vectors
348: from level l to l-1.
350: Logically Collective
352: Input Parameters:
353: + pc - the multigrid context
354: - l - the level (0 is coarsest) to supply [Do not supply 0]
356: Output Parameter:
357: . mat - the restriction matrix
359: Level: advanced
361: .seealso: [](ch_ksp), `PCMG`, `PCMGGetInterpolation()`, `PCMGSetRestriction()`, `PCMGGetRScale()`, `PCMGGetInjection()`
362: @*/
363: PetscErrorCode PCMGGetRestriction(PC pc, PetscInt l, Mat *mat)
364: {
365: PC_MG *mg = (PC_MG *)pc->data;
366: PC_MG_Levels **mglevels = mg->levels;
368: PetscFunctionBegin;
370: if (mat) PetscAssertPointer(mat, 3);
371: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
372: PetscCheck(l > 0 && l < mg->nlevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Level %" PetscInt_FMT " must be in range {1,...,%" PetscInt_FMT "}", l, mg->nlevels - 1);
373: if (!mglevels[l]->restrct && mglevels[l]->interpolate) PetscCall(PCMGSetRestriction(pc, l, mglevels[l]->interpolate));
374: if (mat) *mat = mglevels[l]->restrct;
375: PetscFunctionReturn(PETSC_SUCCESS);
376: }
378: /*@
379: PCMGSetRScale - Sets the pointwise scaling for the restriction operator from level l to l-1.
381: Logically Collective
383: Input Parameters:
384: + pc - the multigrid context
385: . l - the level (0 is coarsest) to supply [Do not supply 0]
386: - rscale - the scaling
388: Level: advanced
390: Note:
391: 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.
392: It is preferable to use `PCMGSetInjection()` to control moving primal vectors.
394: .seealso: [](ch_ksp), `PCMG`, `PCMGSetInterpolation()`, `PCMGSetRestriction()`, `PCMGGetRScale()`, `PCMGSetInjection()`
395: @*/
396: PetscErrorCode PCMGSetRScale(PC pc, PetscInt l, Vec rscale)
397: {
398: PC_MG *mg = (PC_MG *)pc->data;
399: PC_MG_Levels **mglevels = mg->levels;
401: PetscFunctionBegin;
403: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
404: PetscCheck(l > 0 && l < mg->nlevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Level %" PetscInt_FMT " must be in range {1,...,%" PetscInt_FMT "}", l, mg->nlevels - 1);
405: PetscCall(PetscObjectReference((PetscObject)rscale));
406: PetscCall(VecDestroy(&mglevels[l]->rscale));
408: mglevels[l]->rscale = rscale;
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: /*@
413: PCMGGetRScale - Gets the pointwise scaling for the restriction operator from level l to l-1.
415: Collective
417: Input Parameters:
418: + pc - the multigrid context
419: . rscale - the scaling
420: - l - the level (0 is coarsest) to supply [Do not supply 0]
422: Level: advanced
424: Note:
425: 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.
426: It is preferable to use `PCMGGetInjection()` to control moving primal vectors.
428: The scaling is computed using vectors compatible with the matrix returned by `PCMGGetRestriction()`.
429: If this matrix is square, it is treated as an interpolation matrix and its transpose is used. For a square restriction matrix,
430: set the scaling explicitly with `PCMGSetRScale()`.
432: .seealso: [](ch_ksp), `PCMG`, `PCMGSetInterpolation()`, `PCMGGetRestriction()`, `PCMGGetInjection()`, `PCMGSetRScale()`
433: @*/
434: PetscErrorCode PCMGGetRScale(PC pc, PetscInt l, Vec *rscale)
435: {
436: PC_MG *mg = (PC_MG *)pc->data;
437: PC_MG_Levels **mglevels = mg->levels;
439: PetscFunctionBegin;
441: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
442: PetscCheck(l > 0 && l < mg->nlevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Level %" PetscInt_FMT " must be in range {1,...,%" PetscInt_FMT "}", l, mg->nlevels - 1);
443: if (!mglevels[l]->rscale) {
444: Mat R;
445: Vec X, Y, coarse, fine;
446: PetscInt M, N;
448: PetscCall(PCMGGetRestriction(pc, l, &R));
449: PetscCall(MatCreateVecs(R, &X, &Y));
450: PetscCall(MatGetSize(R, &M, &N));
451: if (M < N) {
452: fine = X;
453: coarse = Y;
454: } else {
455: fine = Y;
456: coarse = X;
457: }
458: PetscCall(VecSet(fine, 1.));
459: PetscCall(MatRestrict(R, fine, coarse));
460: PetscCall(VecDestroy(&fine));
461: PetscCall(VecReciprocal(coarse));
462: mglevels[l]->rscale = coarse;
463: }
464: *rscale = mglevels[l]->rscale;
465: PetscFunctionReturn(PETSC_SUCCESS);
466: }
468: /*@
469: PCMGSetInjection - Sets the function to be used to inject primal (i.e. solution) vectors
470: from level l to l-1.
472: Logically Collective
474: Input Parameters:
475: + pc - the multigrid context
476: . l - the level (0 is coarsest) to supply [Do not supply 0]
477: - mat - the injection matrix
479: Level: advanced
481: .seealso: [](ch_ksp), `PCMG`, `PCMGSetRestriction()`
482: @*/
483: PetscErrorCode PCMGSetInjection(PC pc, PetscInt l, Mat mat)
484: {
485: PC_MG *mg = (PC_MG *)pc->data;
486: PC_MG_Levels **mglevels = mg->levels;
488: PetscFunctionBegin;
491: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
492: PetscCheck(l, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Do not set restriction routine for coarsest level");
493: PetscCall(PetscObjectReference((PetscObject)mat));
494: PetscCall(MatDestroy(&mglevels[l]->inject));
496: mglevels[l]->inject = mat;
497: PetscFunctionReturn(PETSC_SUCCESS);
498: }
500: /*@
501: PCMGGetInjection - Gets the function to be used to inject primal vectors (i.e. solutions)
502: from level l to l-1.
504: Logically Collective
506: Input Parameters:
507: + pc - the multigrid context
508: - l - the level (0 is coarsest) to supply [Do not supply 0]
510: Output Parameter:
511: . mat - the restriction matrix (may be `NULL` if no injection is available).
513: Level: advanced
515: .seealso: [](ch_ksp), `PCMG`, `PCMGSetInjection()`, `PCMGetGetRestriction()`
516: @*/
517: PetscErrorCode PCMGGetInjection(PC pc, PetscInt l, Mat *mat)
518: {
519: PC_MG *mg = (PC_MG *)pc->data;
520: PC_MG_Levels **mglevels = mg->levels;
522: PetscFunctionBegin;
524: if (mat) PetscAssertPointer(mat, 3);
525: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
526: PetscCheck(l > 0 && l < mg->nlevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Level %" PetscInt_FMT " must be in range {1,...,%" PetscInt_FMT "}", l, mg->nlevels - 1);
527: if (mat) *mat = mglevels[l]->inject;
528: PetscFunctionReturn(PETSC_SUCCESS);
529: }
531: /*@
532: PCMGGetSmoother - Gets the `KSP` context to be used as smoother for
533: both pre- and post-smoothing. Call both `PCMGGetSmootherUp()` and
534: `PCMGGetSmootherDown()` to use different functions for pre- and
535: post-smoothing.
537: Not Collective, ksp returned is parallel if pc is
539: Input Parameters:
540: + pc - the multigrid context
541: - l - the level (0 is coarsest) to supply
543: Output Parameter:
544: . ksp - the smoother
546: Note:
547: Once you have called this routine, you can call `KSPSetOperators()` on the resulting ksp to provide the operators for the smoother for this level.
548: You can also modify smoother options by calling the various KSPSetXXX() options on this ksp. In addition you can call `KSPGetPC`(ksp,&pc)
549: and modify PC options for the smoother; for example `PCSetType`(pc,`PCSOR`); to use SOR smoothing.
551: Level: advanced
553: .seealso: [](ch_ksp), `PCMG`, `PCMGGetSmootherUp()`, `PCMGGetSmootherDown()`, `PCMGGetCoarseSolve()`
554: @*/
555: PetscErrorCode PCMGGetSmoother(PC pc, PetscInt l, KSP *ksp)
556: {
557: PC_MG *mg = (PC_MG *)pc->data;
558: PC_MG_Levels **mglevels = mg->levels;
560: PetscFunctionBegin;
562: *ksp = mglevels[l]->smoothd;
563: PetscFunctionReturn(PETSC_SUCCESS);
564: }
566: /*@
567: PCMGGetSmootherUp - Gets the KSP context to be used as smoother after
568: coarse grid correction (post-smoother).
570: Not Collective, ksp returned is parallel if pc is
572: Input Parameters:
573: + pc - the multigrid context
574: - l - the level (0 is coarsest) to supply
576: Output Parameter:
577: . ksp - the smoother
579: Level: advanced
581: Note:
582: Calling this will result in a different pre and post smoother so you may need to set options on the pre smoother also
584: .seealso: [](ch_ksp), `PCMG`, `PCMGGetSmootherDown()`
585: @*/
586: PetscErrorCode PCMGGetSmootherUp(PC pc, PetscInt l, KSP *ksp)
587: {
588: PC_MG *mg = (PC_MG *)pc->data;
589: PC_MG_Levels **mglevels = mg->levels;
590: const char *prefix;
591: MPI_Comm comm;
593: PetscFunctionBegin;
595: /*
596: This is called only if user wants a different pre-smoother from post.
597: Thus we check if a different one has already been allocated,
598: if not we allocate it.
599: */
600: PetscCheck(l, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "There is no such thing as a up smoother on the coarse grid");
601: if (mglevels[l]->smoothu == mglevels[l]->smoothd) {
602: KSPType ksptype;
603: PCType pctype;
604: PC ipc;
605: PetscReal rtol, abstol, dtol;
606: PetscInt maxits;
607: KSPNormType normtype;
608: PetscCall(PetscObjectGetComm((PetscObject)mglevels[l]->smoothd, &comm));
609: PetscCall(KSPGetOptionsPrefix(mglevels[l]->smoothd, &prefix));
610: PetscCall(KSPGetTolerances(mglevels[l]->smoothd, &rtol, &abstol, &dtol, &maxits));
611: PetscCall(KSPGetType(mglevels[l]->smoothd, &ksptype));
612: PetscCall(KSPGetNormType(mglevels[l]->smoothd, &normtype));
613: PetscCall(KSPGetPC(mglevels[l]->smoothd, &ipc));
614: PetscCall(PCGetType(ipc, &pctype));
616: PetscCall(KSPCreate(comm, &mglevels[l]->smoothu));
617: PetscCall(KSPSetNestLevel(mglevels[l]->smoothu, pc->kspnestlevel));
618: PetscCall(KSPSetErrorIfNotConverged(mglevels[l]->smoothu, pc->erroriffailure));
619: PetscCall(PetscObjectIncrementTabLevel((PetscObject)mglevels[l]->smoothu, (PetscObject)pc, mglevels[0]->levels - l));
620: PetscCall(KSPSetOptionsPrefix(mglevels[l]->smoothu, prefix));
621: PetscCall(KSPSetTolerances(mglevels[l]->smoothu, rtol, abstol, dtol, maxits));
622: PetscCall(KSPSetType(mglevels[l]->smoothu, ksptype));
623: PetscCall(KSPSetNormType(mglevels[l]->smoothu, normtype));
624: PetscCall(KSPSetConvergenceTest(mglevels[l]->smoothu, KSPConvergedSkip, NULL, NULL));
625: PetscCall(KSPGetPC(mglevels[l]->smoothu, &ipc));
626: PetscCall(PCSetType(ipc, pctype));
627: PetscCall(PetscObjectComposedDataSetInt((PetscObject)mglevels[l]->smoothu, PetscMGLevelId, mglevels[l]->level));
628: }
629: if (ksp) *ksp = mglevels[l]->smoothu;
630: PetscFunctionReturn(PETSC_SUCCESS);
631: }
633: /*@
634: PCMGGetSmootherDown - Gets the `KSP` context to be used as smoother before
635: coarse grid correction (pre-smoother).
637: Not Collective, ksp returned is parallel if pc is
639: Input Parameters:
640: + pc - the multigrid context
641: - l - the level (0 is coarsest) to supply
643: Output Parameter:
644: . ksp - the smoother
646: Level: advanced
648: Note:
649: Calling this will result in a different pre and post smoother so you may need to
650: set options on the post smoother also
652: .seealso: [](ch_ksp), `PCMG`, `PCMGGetSmootherUp()`, `PCMGGetSmoother()`
653: @*/
654: PetscErrorCode PCMGGetSmootherDown(PC pc, PetscInt l, KSP *ksp)
655: {
656: PC_MG *mg = (PC_MG *)pc->data;
657: PC_MG_Levels **mglevels = mg->levels;
659: PetscFunctionBegin;
661: /* make sure smoother up and down are different */
662: if (l) PetscCall(PCMGGetSmootherUp(pc, l, NULL));
663: *ksp = mglevels[l]->smoothd;
664: PetscFunctionReturn(PETSC_SUCCESS);
665: }
667: /*@
668: PCMGSetCycleTypeOnLevel - Sets the type of cycle (aka cycle index) to run on the specified level.
670: Logically Collective
672: Input Parameters:
673: + pc - the multigrid context
674: . l - the level (0 is coarsest)
675: - c - either `PC_MG_CYCLE_V` or `PC_MG_CYCLE_W`
677: Level: advanced
679: .seealso: [](ch_ksp), `PCMG`, `PCMGCycleType`, `PCMGSetCycleType()`
680: @*/
681: PetscErrorCode PCMGSetCycleTypeOnLevel(PC pc, PetscInt l, PCMGCycleType c)
682: {
683: PC_MG *mg = (PC_MG *)pc->data;
684: PC_MG_Levels **mglevels = mg->levels;
686: PetscFunctionBegin;
688: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
691: mglevels[l]->cycles = c;
692: PetscFunctionReturn(PETSC_SUCCESS);
693: }
695: /*@
696: PCMGSetRhs - Sets the vector to be used to store the right-hand side on a particular level.
698: Logically Collective
700: Input Parameters:
701: + pc - the multigrid context
702: . l - the level (0 is coarsest) this is to be used for
703: - c - the Vec
705: Level: advanced
707: Note:
708: 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.
710: .seealso: [](ch_ksp), `PCMG`, `PCMGSetX()`, `PCMGSetR()`
711: @*/
712: PetscErrorCode PCMGSetRhs(PC pc, PetscInt l, Vec c)
713: {
714: PC_MG *mg = (PC_MG *)pc->data;
715: PC_MG_Levels **mglevels = mg->levels;
717: PetscFunctionBegin;
719: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
720: PetscCheck(l != mglevels[0]->levels - 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Do not set rhs for finest level");
721: PetscCall(PetscObjectReference((PetscObject)c));
722: PetscCall(VecDestroy(&mglevels[l]->b));
724: mglevels[l]->b = c;
725: PetscFunctionReturn(PETSC_SUCCESS);
726: }
728: /*@
729: PCMGSetX - Sets the vector to be used to store the solution on a particular level.
731: Logically Collective
733: Input Parameters:
734: + pc - the multigrid context
735: . l - the level (0 is coarsest) this is to be used for (do not supply the finest level)
736: - c - the Vec
738: Level: advanced
740: Note:
741: 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.
743: .seealso: [](ch_ksp), `PCMG`, `PCMGSetRhs()`, `PCMGSetR()`
744: @*/
745: PetscErrorCode PCMGSetX(PC pc, PetscInt l, Vec c)
746: {
747: PC_MG *mg = (PC_MG *)pc->data;
748: PC_MG_Levels **mglevels = mg->levels;
750: PetscFunctionBegin;
752: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
753: PetscCheck(l != mglevels[0]->levels - 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Do not set x for finest level");
754: PetscCall(PetscObjectReference((PetscObject)c));
755: PetscCall(VecDestroy(&mglevels[l]->x));
757: mglevels[l]->x = c;
758: PetscFunctionReturn(PETSC_SUCCESS);
759: }
761: /*@
762: PCMGSetR - Sets the vector to be used to store the residual on a particular level.
764: Logically Collective
766: Input Parameters:
767: + pc - the multigrid context
768: . l - the level (0 is coarsest) this is to be used for
769: - c - the Vec
771: Level: advanced
773: Note:
774: 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.
776: .seealso: [](ch_ksp), `PCMG`, `PCMGSetRhs()`, `PCMGSetX()`
777: @*/
778: PetscErrorCode PCMGSetR(PC pc, PetscInt l, Vec c)
779: {
780: PC_MG *mg = (PC_MG *)pc->data;
781: PC_MG_Levels **mglevels = mg->levels;
783: PetscFunctionBegin;
785: PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
786: PetscCheck(l, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Need not set residual vector for coarse grid");
787: PetscCall(PetscObjectReference((PetscObject)c));
788: PetscCall(VecDestroy(&mglevels[l]->r));
790: mglevels[l]->r = c;
791: PetscFunctionReturn(PETSC_SUCCESS);
792: }