Actual source code: smg.c
2: /*
3: Additive Multigrid V Cycle routine
4: */
5: #include <petsc/private/pcmgimpl.h>
7: PetscErrorCode PCMGACycle_Private(PC pc, PC_MG_Levels **mglevels, PetscBool transpose, PetscBool matapp)
8: {
9: PetscInt i, l = mglevels[0]->levels;
11: PetscFunctionBegin;
12: /* compute RHS on each level */
13: for (i = l - 1; i > 0; i--) {
14: if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
15: if (!transpose) {
16: if (matapp) PetscCall(MatMatRestrict(mglevels[i]->restrct, mglevels[i]->B, &mglevels[i - 1]->B));
17: else PetscCall(MatRestrict(mglevels[i]->restrct, mglevels[i]->b, mglevels[i - 1]->b));
18: } else {
19: if (matapp) PetscCall(MatMatRestrict(mglevels[i]->interpolate, mglevels[i]->B, &mglevels[i - 1]->B));
20: else PetscCall(MatRestrict(mglevels[i]->interpolate, mglevels[i]->b, mglevels[i - 1]->b));
21: }
22: if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
23: }
24: /* solve separately on each level */
25: for (i = 0; i < l; i++) {
26: if (matapp) {
27: if (!mglevels[i]->X) {
28: PetscCall(MatDuplicate(mglevels[i]->B, MAT_DO_NOT_COPY_VALUES, &mglevels[i]->X));
29: } else {
30: PetscCall(MatZeroEntries(mglevels[i]->X));
31: }
32: } else {
33: PetscCall(VecZeroEntries(mglevels[i]->x));
34: }
35: if (mglevels[i]->eventsmoothsolve) PetscCall(PetscLogEventBegin(mglevels[i]->eventsmoothsolve, 0, 0, 0, 0));
36: if (!transpose) {
37: if (matapp) {
38: PetscCall(KSPMatSolve(mglevels[i]->smoothd, mglevels[i]->B, mglevels[i]->X));
39: PetscCall(KSPCheckSolve(mglevels[i]->smoothd, pc, NULL));
40: } else {
41: PetscCall(KSPSolve(mglevels[i]->smoothd, mglevels[i]->b, mglevels[i]->x));
42: PetscCall(KSPCheckSolve(mglevels[i]->smoothd, pc, mglevels[i]->x));
43: }
44: } else {
45: PetscCheck(!matapp, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not supported");
46: PetscCall(KSPSolveTranspose(mglevels[i]->smoothu, mglevels[i]->b, mglevels[i]->x));
47: PetscCall(KSPCheckSolve(mglevels[i]->smoothu, pc, mglevels[i]->x));
48: }
49: if (mglevels[i]->eventsmoothsolve) PetscCall(PetscLogEventEnd(mglevels[i]->eventsmoothsolve, 0, 0, 0, 0));
50: }
51: for (i = 1; i < l; i++) {
52: if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
53: if (!transpose) {
54: if (matapp) PetscCall(MatMatInterpolateAdd(mglevels[i]->interpolate, mglevels[i - 1]->X, mglevels[i]->X, &mglevels[i]->X));
55: else PetscCall(MatInterpolateAdd(mglevels[i]->interpolate, mglevels[i - 1]->x, mglevels[i]->x, mglevels[i]->x));
56: } else {
57: if (matapp) PetscCall(MatMatInterpolateAdd(mglevels[i]->restrct, mglevels[i - 1]->X, mglevels[i]->X, &mglevels[i]->X));
58: else PetscCall(MatInterpolateAdd(mglevels[i]->restrct, mglevels[i - 1]->x, mglevels[i]->x, mglevels[i]->x));
59: }
60: if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
61: }
62: PetscFunctionReturn(PETSC_SUCCESS);
63: }