Actual source code: mgfunc.c

  1: #include <petsc/private/pcmgimpl.h>

  3: /*@C
  4:   PCMGResidualDefault - Default routine to calculate the residual.

  6:   Collective

  8:   Input Parameters:
  9: + mat - the matrix
 10: . b   - the right-hand side
 11: - x   - the approximate solution

 13:   Output Parameter:
 14: . r - location to store the residual

 16:   Level: developer

 18: .seealso: [](ch_ksp), `PCMG`, `PCMGSetResidual()`, `PCMGSetMatResidual()`
 19: @*/
 20: PetscErrorCode PCMGResidualDefault(Mat mat, Vec b, Vec x, Vec r)
 21: {
 22:   PetscFunctionBegin;
 23:   PetscCall(MatResidual(mat, b, x, r));
 24:   PetscFunctionReturn(PETSC_SUCCESS);
 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: [](ch_ksp), `PCMG`, `PCMGSetResidualTranspose()`, `PCMGMatResidualTransposeDefault()`
 43: @*/
 44: PetscErrorCode PCMGResidualTransposeDefault(Mat mat, Vec b, Vec x, Vec r)
 45: {
 46:   PetscFunctionBegin;
 47:   PetscCall(MatMultTranspose(mat, x, r));
 48:   PetscCall(VecAYPX(r, -1.0, b));
 49:   PetscFunctionReturn(PETSC_SUCCESS);
 50: }

 52: /*@C
 53:   PCMGMatResidualDefault - Default routine to calculate the residual.

 55:   Collective

 57:   Input Parameters:
 58: + mat - the matrix
 59: . b   - the right-hand side
 60: - x   - the approximate solution

 62:   Output Parameter:
 63: . r - location to store the residual

 65:   Level: developer

 67: .seealso: [](ch_ksp), `PCMG`, `PCMGSetMatResidual()`, `PCMGResidualDefault()`
 68: @*/
 69: PetscErrorCode PCMGMatResidualDefault(Mat mat, Mat b, Mat x, Mat r)
 70: {
 71:   PetscFunctionBegin;
 72:   PetscCall(MatMatMult(mat, x, MAT_REUSE_MATRIX, PETSC_DEFAULT, &r));
 73:   PetscCall(MatAYPX(r, -1.0, b, UNKNOWN_NONZERO_PATTERN));
 74:   PetscFunctionReturn(PETSC_SUCCESS);
 75: }

 77: /*@C
 78:   PCMGMatResidualTransposeDefault - Default routine to calculate the residual of the transposed linear system

 80:   Collective

 82:   Input Parameters:
 83: + mat - the matrix
 84: . b   - the right-hand side
 85: - x   - the approximate solution

 87:   Output Parameter:
 88: . r - location to store the residual

 90:   Level: developer

 92: .seealso: [](ch_ksp), `PCMG`, `PCMGSetMatResidualTranspose()`
 93: @*/
 94: PetscErrorCode PCMGMatResidualTransposeDefault(Mat mat, Mat b, Mat x, Mat r)
 95: {
 96:   PetscFunctionBegin;
 97:   PetscCall(MatTransposeMatMult(mat, x, MAT_REUSE_MATRIX, PETSC_DEFAULT, &r));
 98:   PetscCall(MatAYPX(r, -1.0, b, UNKNOWN_NONZERO_PATTERN));
 99:   PetscFunctionReturn(PETSC_SUCCESS);
100: }
101: /*@
102:   PCMGGetCoarseSolve - Gets the solver context to be used on the coarse grid.

104:   Not Collective

106:   Input Parameter:
107: . pc - the multigrid context

109:   Output Parameter:
110: . ksp - the coarse grid solver context

112:   Level: advanced

114: .seealso: [](ch_ksp), `PCMG`, `PCMGGetSmootherUp()`, `PCMGGetSmootherDown()`, `PCMGGetSmoother()`
115: @*/
116: PetscErrorCode PCMGGetCoarseSolve(PC pc, KSP *ksp)
117: {
118:   PC_MG         *mg       = (PC_MG *)pc->data;
119:   PC_MG_Levels **mglevels = mg->levels;

121:   PetscFunctionBegin;
123:   *ksp = mglevels[0]->smoothd;
124:   PetscFunctionReturn(PETSC_SUCCESS);
125: }

127: /*@C
128:   PCMGSetResidual - Sets the function to be used to calculate the residual on the lth level.

130:   Logically Collective

132:   Input Parameters:
133: + pc       - the multigrid context
134: . l        - the level (0 is coarsest) to supply
135: . residual - function used to form residual, if none is provided the previously provide one is used, if no
136:               previous one were provided then a default is used
137: - mat      - matrix associated with residual

139:   Level: advanced

141: .seealso: [](ch_ksp), `PCMG`, `PCMGResidualDefault()`
142: @*/
143: PetscErrorCode PCMGSetResidual(PC pc, PetscInt l, PetscErrorCode (*residual)(Mat, Vec, Vec, Vec), Mat mat)
144: {
145:   PC_MG         *mg       = (PC_MG *)pc->data;
146:   PC_MG_Levels **mglevels = mg->levels;

148:   PetscFunctionBegin;
150:   PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
151:   if (residual) mglevels[l]->residual = residual;
152:   if (!mglevels[l]->residual) mglevels[l]->residual = PCMGResidualDefault;
153:   mglevels[l]->matresidual = PCMGMatResidualDefault;
154:   if (mat) PetscCall(PetscObjectReference((PetscObject)mat));
155:   PetscCall(MatDestroy(&mglevels[l]->A));
156:   mglevels[l]->A = mat;
157:   PetscFunctionReturn(PETSC_SUCCESS);
158: }

160: /*@C
161:   PCMGSetResidualTranspose - Sets the function to be used to calculate the residual of the transposed linear system
162:   on the lth level.

164:   Logically Collective

166:   Input Parameters:
167: + pc        - the multigrid context
168: . l         - the level (0 is coarsest) to supply
169: . residualt - function used to form transpose of residual, if none is provided the previously provide one is used, if no
170:                previous one were provided then a default is used
171: - mat       - matrix associated with residual

173:   Level: advanced

175: .seealso: [](ch_ksp), `PCMG`, `PCMGResidualTransposeDefault()`
176: @*/
177: PetscErrorCode PCMGSetResidualTranspose(PC pc, PetscInt l, PetscErrorCode (*residualt)(Mat, Vec, Vec, Vec), Mat mat)
178: {
179:   PC_MG         *mg       = (PC_MG *)pc->data;
180:   PC_MG_Levels **mglevels = mg->levels;

182:   PetscFunctionBegin;
184:   PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
185:   if (residualt) mglevels[l]->residualtranspose = residualt;
186:   if (!mglevels[l]->residualtranspose) mglevels[l]->residualtranspose = PCMGResidualTransposeDefault;
187:   mglevels[l]->matresidualtranspose = PCMGMatResidualTransposeDefault;
188:   if (mat) PetscCall(PetscObjectReference((PetscObject)mat));
189:   PetscCall(MatDestroy(&mglevels[l]->A));
190:   mglevels[l]->A = mat;
191:   PetscFunctionReturn(PETSC_SUCCESS);
192: }

194: /*@
195:   PCMGSetInterpolation - Sets the function to be used to calculate the
196:   interpolation from l-1 to the lth level

198:   Logically Collective

200:   Input Parameters:
201: + pc  - the multigrid context
202: . mat - the interpolation operator
203: - l   - the level (0 is coarsest) to supply [do not supply 0]

205:   Level: advanced

207:   Notes:
208:   Usually this is the same matrix used also to set the restriction
209:   for the same level.

211:   One can pass in the interpolation matrix or its transpose; PETSc figures
212:   out from the matrix size which one it is.

214: .seealso: [](ch_ksp), `PCMG`, `PCMGSetRestriction()`
215: @*/
216: PetscErrorCode PCMGSetInterpolation(PC pc, PetscInt l, Mat mat)
217: {
218:   PC_MG         *mg       = (PC_MG *)pc->data;
219:   PC_MG_Levels **mglevels = mg->levels;

221:   PetscFunctionBegin;
223:   PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
224:   PetscCheck(l, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Do not set interpolation routine for coarsest level");
225:   PetscCall(PetscObjectReference((PetscObject)mat));
226:   PetscCall(MatDestroy(&mglevels[l]->interpolate));

228:   mglevels[l]->interpolate = mat;
229:   PetscFunctionReturn(PETSC_SUCCESS);
230: }

232: /*@
233:   PCMGSetOperators - Sets operator and preconditioning matrix for lth level

235:   Logically Collective

237:   Input Parameters:
238: + pc   - the multigrid context
239: . Amat - the operator
240: . Pmat - the preconditioning operator
241: - l    - the level (0 is the coarsest) to supply

243:   Level: advanced

245: .seealso: [](ch_ksp), `PCMG`, `PCMGSetGalerkin()`, `PCMGSetRestriction()`, `PCMGSetInterpolation()`
246: @*/
247: PetscErrorCode PCMGSetOperators(PC pc, PetscInt l, Mat Amat, Mat Pmat)
248: {
249:   PC_MG         *mg       = (PC_MG *)pc->data;
250:   PC_MG_Levels **mglevels = mg->levels;

252:   PetscFunctionBegin;
256:   PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
257:   PetscCall(KSPSetOperators(mglevels[l]->smoothd, Amat, Pmat));
258:   PetscFunctionReturn(PETSC_SUCCESS);
259: }

261: /*@
262:   PCMGGetInterpolation - Gets the function to be used to calculate the
263:   interpolation from l-1 to the lth level

265:   Logically Collective

267:   Input Parameters:
268: + pc - the multigrid context
269: - l  - the level (0 is coarsest) to supply [Do not supply 0]

271:   Output Parameter:
272: . mat - the interpolation matrix, can be `NULL`

274:   Level: advanced

276: .seealso: [](ch_ksp), `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetRScale()`
277: @*/
278: PetscErrorCode PCMGGetInterpolation(PC pc, PetscInt l, Mat *mat)
279: {
280:   PC_MG         *mg       = (PC_MG *)pc->data;
281:   PC_MG_Levels **mglevels = mg->levels;

283:   PetscFunctionBegin;
285:   if (mat) PetscAssertPointer(mat, 3);
286:   PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
287:   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);
288:   if (!mglevels[l]->interpolate && mglevels[l]->restrct) PetscCall(PCMGSetInterpolation(pc, l, mglevels[l]->restrct));
289:   if (mat) *mat = mglevels[l]->interpolate;
290:   PetscFunctionReturn(PETSC_SUCCESS);
291: }

293: /*@
294:   PCMGSetRestriction - Sets the function to be used to restrict dual vectors
295:   from level l to l-1.

297:   Logically Collective

299:   Input Parameters:
300: + pc  - the multigrid context
301: . l   - the level (0 is coarsest) to supply [Do not supply 0]
302: - mat - the restriction matrix

304:   Level: advanced

306:   Notes:
307:   Usually this is the same matrix used also to set the interpolation
308:   for the same level.

310:   One can pass in the interpolation matrix or its transpose; PETSc figures
311:   out from the matrix size which one it is.

313:   If you do not set this, the transpose of the `Mat` set with `PCMGSetInterpolation()`
314:   is used.

316: .seealso: [](ch_ksp), `PCMG`, `PCMGSetInterpolation()`
317: @*/
318: PetscErrorCode PCMGSetRestriction(PC pc, PetscInt l, Mat mat)
319: {
320:   PC_MG         *mg       = (PC_MG *)pc->data;
321:   PC_MG_Levels **mglevels = mg->levels;

323:   PetscFunctionBegin;
326:   PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
327:   PetscCheck(l, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Do not set restriction routine for coarsest level");
328:   PetscCall(PetscObjectReference((PetscObject)mat));
329:   PetscCall(MatDestroy(&mglevels[l]->restrct));

331:   mglevels[l]->restrct = mat;
332:   PetscFunctionReturn(PETSC_SUCCESS);
333: }

335: /*@
336:   PCMGGetRestriction - Gets the function to be used to restrict dual vectors
337:   from level l to l-1.

339:   Logically Collective

341:   Input Parameters:
342: + pc - the multigrid context
343: - l  - the level (0 is coarsest) to supply [Do not supply 0]

345:   Output Parameter:
346: . mat - the restriction matrix

348:   Level: advanced

350: .seealso: [](ch_ksp), `PCMG`, `PCMGGetInterpolation()`, `PCMGSetRestriction()`, `PCMGGetRScale()`, `PCMGGetInjection()`
351: @*/
352: PetscErrorCode PCMGGetRestriction(PC pc, PetscInt l, Mat *mat)
353: {
354:   PC_MG         *mg       = (PC_MG *)pc->data;
355:   PC_MG_Levels **mglevels = mg->levels;

357:   PetscFunctionBegin;
359:   if (mat) PetscAssertPointer(mat, 3);
360:   PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
361:   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);
362:   if (!mglevels[l]->restrct && mglevels[l]->interpolate) PetscCall(PCMGSetRestriction(pc, l, mglevels[l]->interpolate));
363:   if (mat) *mat = mglevels[l]->restrct;
364:   PetscFunctionReturn(PETSC_SUCCESS);
365: }

367: /*@
368:   PCMGSetRScale - Sets the pointwise scaling for the restriction operator from level l to l-1.

370:   Logically Collective

372:   Input Parameters:
373: + pc     - the multigrid context
374: . l      - the level (0 is coarsest) to supply [Do not supply 0]
375: - rscale - the scaling

377:   Level: advanced

379:   Note:
380:   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.
381:   It is preferable to use `PCMGSetInjection()` to control moving primal vectors.

383: .seealso: [](ch_ksp), `PCMG`, `PCMGSetInterpolation()`, `PCMGSetRestriction()`, `PCMGGetRScale()`, `PCMGSetInjection()`
384: @*/
385: PetscErrorCode PCMGSetRScale(PC pc, PetscInt l, Vec rscale)
386: {
387:   PC_MG         *mg       = (PC_MG *)pc->data;
388:   PC_MG_Levels **mglevels = mg->levels;

390:   PetscFunctionBegin;
392:   PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
393:   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);
394:   PetscCall(PetscObjectReference((PetscObject)rscale));
395:   PetscCall(VecDestroy(&mglevels[l]->rscale));

397:   mglevels[l]->rscale = rscale;
398:   PetscFunctionReturn(PETSC_SUCCESS);
399: }

401: /*@
402:   PCMGGetRScale - Gets the pointwise scaling for the restriction operator from level l to l-1.

404:   Collective

406:   Input Parameters:
407: + pc     - the multigrid context
408: . rscale - the scaling
409: - l      - the level (0 is coarsest) to supply [Do not supply 0]

411:   Level: advanced

413:   Note:
414:   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.
415:   It is preferable to use `PCMGGetInjection()` to control moving primal vectors.

417: .seealso: [](ch_ksp), `PCMG`, `PCMGSetInterpolation()`, `PCMGGetRestriction()`, `PCMGGetInjection()`
418: @*/
419: PetscErrorCode PCMGGetRScale(PC pc, PetscInt l, Vec *rscale)
420: {
421:   PC_MG         *mg       = (PC_MG *)pc->data;
422:   PC_MG_Levels **mglevels = mg->levels;

424:   PetscFunctionBegin;
426:   PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
427:   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);
428:   if (!mglevels[l]->rscale) {
429:     Mat      R;
430:     Vec      X, Y, coarse, fine;
431:     PetscInt M, N;

433:     PetscCall(PCMGGetRestriction(pc, l, &R));
434:     PetscCall(MatCreateVecs(R, &X, &Y));
435:     PetscCall(MatGetSize(R, &M, &N));
436:     PetscCheck(N != M, PetscObjectComm((PetscObject)R), PETSC_ERR_SUP, "Restriction matrix is square, cannot determine which Vec is coarser");
437:     if (M < N) {
438:       fine   = X;
439:       coarse = Y;
440:     } else {
441:       fine   = Y;
442:       coarse = X;
443:     }
444:     PetscCall(VecSet(fine, 1.));
445:     PetscCall(MatRestrict(R, fine, coarse));
446:     PetscCall(VecDestroy(&fine));
447:     PetscCall(VecReciprocal(coarse));
448:     mglevels[l]->rscale = coarse;
449:   }
450:   *rscale = mglevels[l]->rscale;
451:   PetscFunctionReturn(PETSC_SUCCESS);
452: }

454: /*@
455:   PCMGSetInjection - Sets the function to be used to inject primal vectors
456:   from level l to l-1.

458:   Logically Collective

460:   Input Parameters:
461: + pc  - the multigrid context
462: . l   - the level (0 is coarsest) to supply [Do not supply 0]
463: - mat - the injection matrix

465:   Level: advanced

467: .seealso: [](ch_ksp), `PCMG`, `PCMGSetRestriction()`
468: @*/
469: PetscErrorCode PCMGSetInjection(PC pc, PetscInt l, Mat mat)
470: {
471:   PC_MG         *mg       = (PC_MG *)pc->data;
472:   PC_MG_Levels **mglevels = mg->levels;

474:   PetscFunctionBegin;
477:   PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
478:   PetscCheck(l, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Do not set restriction routine for coarsest level");
479:   PetscCall(PetscObjectReference((PetscObject)mat));
480:   PetscCall(MatDestroy(&mglevels[l]->inject));

482:   mglevels[l]->inject = mat;
483:   PetscFunctionReturn(PETSC_SUCCESS);
484: }

486: /*@
487:   PCMGGetInjection - Gets the function to be used to inject primal vectors
488:   from level l to l-1.

490:   Logically Collective

492:   Input Parameters:
493: + pc - the multigrid context
494: - l  - the level (0 is coarsest) to supply [Do not supply 0]

496:   Output Parameter:
497: . mat - the restriction matrix (may be `NULL` if no injection is available).

499:   Level: advanced

501: .seealso: [](ch_ksp), `PCMG`, `PCMGSetInjection()`, `PCMGetGetRestriction()`
502: @*/
503: PetscErrorCode PCMGGetInjection(PC pc, PetscInt l, Mat *mat)
504: {
505:   PC_MG         *mg       = (PC_MG *)pc->data;
506:   PC_MG_Levels **mglevels = mg->levels;

508:   PetscFunctionBegin;
510:   if (mat) PetscAssertPointer(mat, 3);
511:   PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
512:   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);
513:   if (mat) *mat = mglevels[l]->inject;
514:   PetscFunctionReturn(PETSC_SUCCESS);
515: }

517: /*@
518:   PCMGGetSmoother - Gets the `KSP` context to be used as smoother for
519:   both pre- and post-smoothing.  Call both `PCMGGetSmootherUp()` and
520:   `PCMGGetSmootherDown()` to use different functions for pre- and
521:   post-smoothing.

523:   Not Collective, ksp returned is parallel if pc is

525:   Input Parameters:
526: + pc - the multigrid context
527: - l  - the level (0 is coarsest) to supply

529:   Output Parameter:
530: . ksp - the smoother

532:   Note:
533:   Once you have called this routine, you can call `KSPSetOperators()` on the resulting ksp to provide the operators for the smoother for this level.
534:   You can also modify smoother options by calling the various KSPSetXXX() options on this ksp. In addition you can call `KSPGetPC`(ksp,&pc)
535:   and modify PC options for the smoother; for example `PCSetType`(pc,`PCSOR`); to use SOR smoothing.

537:   Level: advanced

539: .seealso: [](ch_ksp), `PCMG`, ``PCMGGetSmootherUp()`, `PCMGGetSmootherDown()`, `PCMGGetCoarseSolve()`
540: @*/
541: PetscErrorCode PCMGGetSmoother(PC pc, PetscInt l, KSP *ksp)
542: {
543:   PC_MG         *mg       = (PC_MG *)pc->data;
544:   PC_MG_Levels **mglevels = mg->levels;

546:   PetscFunctionBegin;
548:   *ksp = mglevels[l]->smoothd;
549:   PetscFunctionReturn(PETSC_SUCCESS);
550: }

552: /*@
553:   PCMGGetSmootherUp - Gets the KSP context to be used as smoother after
554:   coarse grid correction (post-smoother).

556:   Not Collective, ksp returned is parallel if pc is

558:   Input Parameters:
559: + pc - the multigrid context
560: - l  - the level (0 is coarsest) to supply

562:   Output Parameter:
563: . ksp - the smoother

565:   Level: advanced

567:   Note:
568:   Calling this will result in a different pre and post smoother so you may need to set options on the pre smoother also

570: .seealso: [](ch_ksp), `PCMG`, `PCMGGetSmootherDown()`
571: @*/
572: PetscErrorCode PCMGGetSmootherUp(PC pc, PetscInt l, KSP *ksp)
573: {
574:   PC_MG         *mg       = (PC_MG *)pc->data;
575:   PC_MG_Levels **mglevels = mg->levels;
576:   const char    *prefix;
577:   MPI_Comm       comm;

579:   PetscFunctionBegin;
581:   /*
582:      This is called only if user wants a different pre-smoother from post.
583:      Thus we check if a different one has already been allocated,
584:      if not we allocate it.
585:   */
586:   PetscCheck(l, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "There is no such thing as a up smoother on the coarse grid");
587:   if (mglevels[l]->smoothu == mglevels[l]->smoothd) {
588:     KSPType     ksptype;
589:     PCType      pctype;
590:     PC          ipc;
591:     PetscReal   rtol, abstol, dtol;
592:     PetscInt    maxits;
593:     KSPNormType normtype;
594:     PetscCall(PetscObjectGetComm((PetscObject)mglevels[l]->smoothd, &comm));
595:     PetscCall(KSPGetOptionsPrefix(mglevels[l]->smoothd, &prefix));
596:     PetscCall(KSPGetTolerances(mglevels[l]->smoothd, &rtol, &abstol, &dtol, &maxits));
597:     PetscCall(KSPGetType(mglevels[l]->smoothd, &ksptype));
598:     PetscCall(KSPGetNormType(mglevels[l]->smoothd, &normtype));
599:     PetscCall(KSPGetPC(mglevels[l]->smoothd, &ipc));
600:     PetscCall(PCGetType(ipc, &pctype));

602:     PetscCall(KSPCreate(comm, &mglevels[l]->smoothu));
603:     PetscCall(KSPSetNestLevel(mglevels[l]->smoothu, pc->kspnestlevel));
604:     PetscCall(KSPSetErrorIfNotConverged(mglevels[l]->smoothu, pc->erroriffailure));
605:     PetscCall(PetscObjectIncrementTabLevel((PetscObject)mglevels[l]->smoothu, (PetscObject)pc, mglevels[0]->levels - l));
606:     PetscCall(KSPSetOptionsPrefix(mglevels[l]->smoothu, prefix));
607:     PetscCall(KSPSetTolerances(mglevels[l]->smoothu, rtol, abstol, dtol, maxits));
608:     PetscCall(KSPSetType(mglevels[l]->smoothu, ksptype));
609:     PetscCall(KSPSetNormType(mglevels[l]->smoothu, normtype));
610:     PetscCall(KSPSetConvergenceTest(mglevels[l]->smoothu, KSPConvergedSkip, NULL, NULL));
611:     PetscCall(KSPGetPC(mglevels[l]->smoothu, &ipc));
612:     PetscCall(PCSetType(ipc, pctype));
613:     PetscCall(PetscObjectComposedDataSetInt((PetscObject)mglevels[l]->smoothu, PetscMGLevelId, mglevels[l]->level));
614:   }
615:   if (ksp) *ksp = mglevels[l]->smoothu;
616:   PetscFunctionReturn(PETSC_SUCCESS);
617: }

619: /*@
620:   PCMGGetSmootherDown - Gets the `KSP` context to be used as smoother before
621:   coarse grid correction (pre-smoother).

623:   Not Collective, ksp returned is parallel if pc is

625:   Input Parameters:
626: + pc - the multigrid context
627: - l  - the level (0 is coarsest) to supply

629:   Output Parameter:
630: . ksp - the smoother

632:   Level: advanced

634:   Note:
635:   Calling this will result in a different pre and post smoother so you may need to
636:   set options on the post smoother also

638: .seealso: [](ch_ksp), `PCMG`, `PCMGGetSmootherUp()`, `PCMGGetSmoother()`
639: @*/
640: PetscErrorCode PCMGGetSmootherDown(PC pc, PetscInt l, KSP *ksp)
641: {
642:   PC_MG         *mg       = (PC_MG *)pc->data;
643:   PC_MG_Levels **mglevels = mg->levels;

645:   PetscFunctionBegin;
647:   /* make sure smoother up and down are different */
648:   if (l) PetscCall(PCMGGetSmootherUp(pc, l, NULL));
649:   *ksp = mglevels[l]->smoothd;
650:   PetscFunctionReturn(PETSC_SUCCESS);
651: }

653: /*@
654:   PCMGSetCycleTypeOnLevel - Sets the type of cycle (aka cycle index) to run on the specified level.

656:   Logically Collective

658:   Input Parameters:
659: + pc - the multigrid context
660: . l  - the level (0 is coarsest)
661: - c  - either `PC_MG_CYCLE_V` or `PC_MG_CYCLE_W`

663:   Level: advanced

665: .seealso: [](ch_ksp), `PCMG`, `PCMGCycleType`, `PCMGSetCycleType()`
666: @*/
667: PetscErrorCode PCMGSetCycleTypeOnLevel(PC pc, PetscInt l, PCMGCycleType c)
668: {
669:   PC_MG         *mg       = (PC_MG *)pc->data;
670:   PC_MG_Levels **mglevels = mg->levels;

672:   PetscFunctionBegin;
674:   PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
677:   mglevels[l]->cycles = c;
678:   PetscFunctionReturn(PETSC_SUCCESS);
679: }

681: /*@
682:   PCMGSetRhs - Sets the vector to be used to store the right-hand side on a particular level.

684:   Logically Collective

686:   Input Parameters:
687: + pc - the multigrid context
688: . l  - the level (0 is coarsest) this is to be used for
689: - c  - the Vec

691:   Level: advanced

693:   Note:
694:   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.

696: .seealso: [](ch_ksp), `PCMG`, `PCMGSetX()`, `PCMGSetR()`
697: @*/
698: PetscErrorCode PCMGSetRhs(PC pc, PetscInt l, Vec c)
699: {
700:   PC_MG         *mg       = (PC_MG *)pc->data;
701:   PC_MG_Levels **mglevels = mg->levels;

703:   PetscFunctionBegin;
705:   PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
706:   PetscCheck(l != mglevels[0]->levels - 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Do not set rhs for finest level");
707:   PetscCall(PetscObjectReference((PetscObject)c));
708:   PetscCall(VecDestroy(&mglevels[l]->b));

710:   mglevels[l]->b = c;
711:   PetscFunctionReturn(PETSC_SUCCESS);
712: }

714: /*@
715:   PCMGSetX - Sets the vector to be used to store the solution on a particular level.

717:   Logically Collective

719:   Input Parameters:
720: + pc - the multigrid context
721: . l  - the level (0 is coarsest) this is to be used for (do not supply the finest level)
722: - c  - the Vec

724:   Level: advanced

726:   Note:
727:   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.

729: .seealso: [](ch_ksp), `PCMG`, `PCMGSetRhs()`, `PCMGSetR()`
730: @*/
731: PetscErrorCode PCMGSetX(PC pc, PetscInt l, Vec c)
732: {
733:   PC_MG         *mg       = (PC_MG *)pc->data;
734:   PC_MG_Levels **mglevels = mg->levels;

736:   PetscFunctionBegin;
738:   PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
739:   PetscCheck(l != mglevels[0]->levels - 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Do not set x for finest level");
740:   PetscCall(PetscObjectReference((PetscObject)c));
741:   PetscCall(VecDestroy(&mglevels[l]->x));

743:   mglevels[l]->x = c;
744:   PetscFunctionReturn(PETSC_SUCCESS);
745: }

747: /*@
748:   PCMGSetR - Sets the vector to be used to store the residual on a particular level.

750:   Logically Collective

752:   Input Parameters:
753: + pc - the multigrid context
754: . l  - the level (0 is coarsest) this is to be used for
755: - c  - the Vec

757:   Level: advanced

759:   Note:
760:   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.

762: .seealso: [](ch_ksp), `PCMG`, `PCMGSetRhs()`, `PCMGSetX()`
763: @*/
764: PetscErrorCode PCMGSetR(PC pc, PetscInt l, Vec c)
765: {
766:   PC_MG         *mg       = (PC_MG *)pc->data;
767:   PC_MG_Levels **mglevels = mg->levels;

769:   PetscFunctionBegin;
771:   PetscCheck(mglevels, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONGSTATE, "Must set MG levels before calling");
772:   PetscCheck(l, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Need not set residual vector for coarse grid");
773:   PetscCall(PetscObjectReference((PetscObject)c));
774:   PetscCall(VecDestroy(&mglevels[l]->r));

776:   mglevels[l]->r = c;
777:   PetscFunctionReturn(PETSC_SUCCESS);
778: }