Actual source code: fmg.c
1: /*
2: Full multigrid using either additive or multiplicative V or W cycle
3: */
4: #include <petsc/private/pcmgimpl.h>
6: /*
7: Transpose of the kaskade cycle, and of the full cycle when full is PETSC_TRUE.
9: The kaskade cycle applies M_i = E_i P_i M_{i-1} R_i + N_i on level i, where N_i is the level solve and
10: E_i = I - N_i A_i, hence M_i^T = R_i^T M_{i-1}^T P_i^T E_i^T + N_i^T with M_0^T = N_0^T. N_i^T is applied by
11: KSPSolveTranspose() from a zero initial guess and E_i^T by the transposed residual that follows it.
13: The full cycle is the same recursion with the level multiplicative cycle in place of the single smoother solve,
14: so the transposed multiplicative cycle plays the smoother's role there.
15: */
16: static PetscErrorCode PCMGKCycleTranspose_Private(PC pc, PC_MG_Levels **mglevels, PetscBool full)
17: {
18: PetscInt l = mglevels[0]->levels;
20: PetscFunctionBegin;
21: /* apply the transposed level solves from the finest level down, restricting the transposed residual with P^T */
22: for (PetscInt i = l - 1; i >= 0; i--) {
23: PetscCall(VecZeroEntries(mglevels[i]->x));
24: if (full) PetscCall(PCMGMCycle_Private(pc, &mglevels[i], PETSC_TRUE, PETSC_FALSE, NULL));
25: else {
26: if (mglevels[i]->eventsmoothsolve) PetscCall(PetscLogEventBegin(mglevels[i]->eventsmoothsolve, 0, 0, 0, 0));
27: PetscCall(KSPSolveTranspose(mglevels[i]->smoothd, mglevels[i]->b, mglevels[i]->x));
28: PetscCall(KSPCheckSolve(mglevels[i]->smoothd, pc, mglevels[i]->x));
29: if (mglevels[i]->eventsmoothsolve) PetscCall(PetscLogEventEnd(mglevels[i]->eventsmoothsolve, 0, 0, 0, 0));
30: }
31: if (i) {
32: if (mglevels[i]->eventresidual) PetscCall(PetscLogEventBegin(mglevels[i]->eventresidual, 0, 0, 0, 0));
33: PetscCall((*mglevels[i]->residualtranspose)(mglevels[i]->A, mglevels[i]->b, mglevels[i]->x, mglevels[i]->r));
34: if (mglevels[i]->eventresidual) PetscCall(PetscLogEventEnd(mglevels[i]->eventresidual, 0, 0, 0, 0));
35: if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
36: PetscCall(MatRestrict(mglevels[i]->interpolate, mglevels[i]->r, mglevels[i - 1]->b));
37: if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
38: }
39: }
41: /* accumulate the coarse contributions back up through the levels with R^T */
42: for (PetscInt i = 1; i < l; i++) {
43: if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
44: PetscCall(MatInterpolateAdd(mglevels[i]->restrct, mglevels[i - 1]->x, mglevels[i]->x, mglevels[i]->x));
45: if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
46: }
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: PetscErrorCode PCMGFCycle_Private(PC pc, PC_MG_Levels **mglevels, PetscBool transpose, PetscBool matapp)
51: {
52: PetscInt l = mglevels[0]->levels;
54: PetscFunctionBegin;
55: if (!transpose) {
56: /* restrict the RHS through all levels to coarsest. */
57: for (PetscInt i = l - 1; i > 0; i--) {
58: if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
59: if (matapp) PetscCall(MatMatRestrict(mglevels[i]->restrct, mglevels[i]->B, &mglevels[i - 1]->B));
60: else PetscCall(MatRestrict(mglevels[i]->restrct, mglevels[i]->b, mglevels[i - 1]->b));
61: if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
62: }
64: /* work our way up through the levels */
65: if (matapp) {
66: if (!mglevels[0]->X) PetscCall(MatDuplicate(mglevels[0]->B, MAT_DO_NOT_COPY_VALUES, &mglevels[0]->X));
67: else PetscCall(MatZeroEntries(mglevels[0]->X));
68: } else PetscCall(VecZeroEntries(mglevels[0]->x));
69: for (PetscInt i = 0; i < l - 1; i++) {
70: PetscCall(PCMGMCycle_Private(pc, &mglevels[i], transpose, matapp, NULL));
71: if (mglevels[i + 1]->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels[i + 1]->eventinterprestrict, 0, 0, 0, 0));
72: if (matapp) PetscCall(MatMatInterpolate(mglevels[i + 1]->interpolate, mglevels[i]->X, &mglevels[i + 1]->X));
73: else PetscCall(MatInterpolate(mglevels[i + 1]->interpolate, mglevels[i]->x, mglevels[i + 1]->x));
74: if (mglevels[i + 1]->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels[i + 1]->eventinterprestrict, 0, 0, 0, 0));
75: }
76: PetscCall(PCMGMCycle_Private(pc, &mglevels[l - 1], transpose, matapp, NULL));
77: } else {
78: PetscCheck(!matapp, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not supported");
79: PetscCall(PCMGKCycleTranspose_Private(pc, mglevels, PETSC_TRUE));
80: }
81: PetscFunctionReturn(PETSC_SUCCESS);
82: }
84: PetscErrorCode PCMGKCycle_Private(PC pc, PC_MG_Levels **mglevels, PetscBool transpose, PetscBool matapp)
85: {
86: PetscInt l = mglevels[0]->levels;
88: PetscFunctionBegin;
89: if (transpose) {
90: PetscCheck(!matapp, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not supported");
91: PetscCall(PCMGKCycleTranspose_Private(pc, mglevels, PETSC_FALSE));
92: PetscFunctionReturn(PETSC_SUCCESS);
93: }
94: /* restrict the RHS through all levels to coarsest. */
95: for (PetscInt i = l - 1; i > 0; i--) {
96: if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
97: if (matapp) PetscCall(MatMatRestrict(mglevels[i]->restrct, mglevels[i]->B, &mglevels[i - 1]->B));
98: else PetscCall(MatRestrict(mglevels[i]->restrct, mglevels[i]->b, mglevels[i - 1]->b));
99: if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
100: }
102: /* work our way up through the levels */
103: if (matapp) {
104: if (!mglevels[0]->X) PetscCall(MatDuplicate(mglevels[0]->B, MAT_DO_NOT_COPY_VALUES, &mglevels[0]->X));
105: else {
106: PetscCall(MatZeroEntries(mglevels[0]->X));
107: }
108: } else {
109: PetscCall(VecZeroEntries(mglevels[0]->x));
110: }
111: for (PetscInt i = 0; i < l - 1; i++) {
112: if (mglevels[i]->eventsmoothsolve) PetscCall(PetscLogEventBegin(mglevels[i]->eventsmoothsolve, 0, 0, 0, 0));
113: if (matapp) {
114: PetscCall(KSPMatSolve(mglevels[i]->smoothd, mglevels[i]->B, mglevels[i]->X));
115: PetscCall(KSPCheckMatSolve(mglevels[i]->smoothd, pc, mglevels[i]->X));
116: } else {
117: PetscCall(KSPSolve(mglevels[i]->smoothd, mglevels[i]->b, mglevels[i]->x));
118: PetscCall(KSPCheckSolve(mglevels[i]->smoothd, pc, mglevels[i]->x));
119: }
120: if (mglevels[i]->eventsmoothsolve) PetscCall(PetscLogEventEnd(mglevels[i]->eventsmoothsolve, 0, 0, 0, 0));
121: if (mglevels[i + 1]->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels[i + 1]->eventinterprestrict, 0, 0, 0, 0));
122: if (matapp) PetscCall(MatMatInterpolate(mglevels[i + 1]->interpolate, mglevels[i]->X, &mglevels[i + 1]->X));
123: else PetscCall(MatInterpolate(mglevels[i + 1]->interpolate, mglevels[i]->x, mglevels[i + 1]->x));
124: if (mglevels[i + 1]->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels[i + 1]->eventinterprestrict, 0, 0, 0, 0));
125: }
126: if (mglevels[l - 1]->eventsmoothsolve) PetscCall(PetscLogEventBegin(mglevels[l - 1]->eventsmoothsolve, 0, 0, 0, 0));
127: if (matapp) {
128: PetscCall(KSPMatSolve(mglevels[l - 1]->smoothd, mglevels[l - 1]->B, mglevels[l - 1]->X));
129: PetscCall(KSPCheckMatSolve(mglevels[l - 1]->smoothd, pc, mglevels[l - 1]->X));
130: } else {
131: PetscCall(KSPSolve(mglevels[l - 1]->smoothd, mglevels[l - 1]->b, mglevels[l - 1]->x));
132: PetscCall(KSPCheckSolve(mglevels[l - 1]->smoothd, pc, mglevels[l - 1]->x));
133: }
134: if (mglevels[l - 1]->eventsmoothsolve) PetscCall(PetscLogEventEnd(mglevels[l - 1]->eventsmoothsolve, 0, 0, 0, 0));
135: PetscFunctionReturn(PETSC_SUCCESS);
136: }