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: PetscErrorCode PCMGFCycle_Private(PC pc, PC_MG_Levels **mglevels, PetscBool transpose, PetscBool matapp)
7: {
8: PetscInt i, l = mglevels[0]->levels;
10: if (!transpose) {
11: /* restrict the RHS through all levels to coarsest. */
12: for (i = l - 1; i > 0; i--) {
13: if (mglevels[i]->eventinterprestrict) PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0);
14: if (matapp) MatMatRestrict(mglevels[i]->restrct, mglevels[i]->B, &mglevels[i - 1]->B);
15: else MatRestrict(mglevels[i]->restrct, mglevels[i]->b, mglevels[i - 1]->b);
16: if (mglevels[i]->eventinterprestrict) PetscLogEventEnd(mglevels[i]->eventinterprestrict, 0, 0, 0, 0);
17: }
19: /* work our way up through the levels */
20: if (matapp) {
21: if (!mglevels[0]->X) {
22: MatDuplicate(mglevels[0]->B, MAT_DO_NOT_COPY_VALUES, &mglevels[0]->X);
23: } else {
24: MatZeroEntries(mglevels[0]->X);
25: }
26: } else {
27: VecZeroEntries(mglevels[0]->x);
28: }
29: for (i = 0; i < l - 1; i++) {
30: PCMGMCycle_Private(pc, &mglevels[i], transpose, matapp, NULL);
31: if (mglevels[i + 1]->eventinterprestrict) PetscLogEventBegin(mglevels[i + 1]->eventinterprestrict, 0, 0, 0, 0);
32: if (matapp) MatMatInterpolate(mglevels[i + 1]->interpolate, mglevels[i]->X, &mglevels[i + 1]->X);
33: else MatInterpolate(mglevels[i + 1]->interpolate, mglevels[i]->x, mglevels[i + 1]->x);
34: if (mglevels[i + 1]->eventinterprestrict) PetscLogEventEnd(mglevels[i + 1]->eventinterprestrict, 0, 0, 0, 0);
35: }
36: PCMGMCycle_Private(pc, &mglevels[l - 1], transpose, matapp, NULL);
37: } else {
38: PCMGMCycle_Private(pc, &mglevels[l - 1], transpose, matapp, NULL);
39: for (i = l - 2; i >= 0; i--) {
40: if (mglevels[i + 1]->eventinterprestrict) PetscLogEventBegin(mglevels[i + 1]->eventinterprestrict, 0, 0, 0, 0);
41: if (matapp) MatMatRestrict(mglevels[i + 1]->interpolate, mglevels[i + 1]->X, &mglevels[i]->X);
42: else MatRestrict(mglevels[i + 1]->interpolate, mglevels[i + 1]->x, mglevels[i]->x);
43: if (mglevels[i + 1]->eventinterprestrict) PetscLogEventEnd(mglevels[i + 1]->eventinterprestrict, 0, 0, 0, 0);
44: PCMGMCycle_Private(pc, &mglevels[i], transpose, matapp, NULL);
45: }
46: for (i = 1; i < l; i++) {
47: if (mglevels[i]->eventinterprestrict) PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0);
48: if (matapp) MatMatInterpolate(mglevels[i]->restrct, mglevels[i - 1]->B, &mglevels[i]->B);
49: else MatInterpolate(mglevels[i]->restrct, mglevels[i - 1]->b, mglevels[i]->b);
50: if (mglevels[i]->eventinterprestrict) PetscLogEventEnd(mglevels[i]->eventinterprestrict, 0, 0, 0, 0);
51: }
52: }
53: return 0;
54: }
56: PetscErrorCode PCMGKCycle_Private(PC pc, PC_MG_Levels **mglevels, PetscBool transpose, PetscBool matapp)
57: {
58: PetscInt i, l = mglevels[0]->levels;
60: /* restrict the RHS through all levels to coarsest. */
61: for (i = l - 1; i > 0; i--) {
62: if (mglevels[i]->eventinterprestrict) PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0);
63: if (matapp) MatMatRestrict(mglevels[i]->restrct, mglevels[i]->B, &mglevels[i - 1]->B);
64: else MatRestrict(mglevels[i]->restrct, mglevels[i]->b, mglevels[i - 1]->b);
65: if (mglevels[i]->eventinterprestrict) PetscLogEventEnd(mglevels[i]->eventinterprestrict, 0, 0, 0, 0);
66: }
68: /* work our way up through the levels */
69: if (matapp) {
70: if (!mglevels[0]->X) {
71: MatDuplicate(mglevels[0]->B, MAT_DO_NOT_COPY_VALUES, &mglevels[0]->X);
72: } else {
73: MatZeroEntries(mglevels[0]->X);
74: }
75: } else {
76: VecZeroEntries(mglevels[0]->x);
77: }
78: for (i = 0; i < l - 1; i++) {
79: if (mglevels[i]->eventsmoothsolve) PetscLogEventBegin(mglevels[i]->eventsmoothsolve, 0, 0, 0, 0);
80: if (matapp) {
81: KSPMatSolve(mglevels[i]->smoothd, mglevels[i]->B, mglevels[i]->X);
82: KSPCheckSolve(mglevels[i]->smoothd, pc, NULL);
83: } else {
84: KSPSolve(mglevels[i]->smoothd, mglevels[i]->b, mglevels[i]->x);
85: KSPCheckSolve(mglevels[i]->smoothd, pc, mglevels[i]->x);
86: }
87: if (mglevels[i]->eventsmoothsolve) PetscLogEventEnd(mglevels[i]->eventsmoothsolve, 0, 0, 0, 0);
88: if (mglevels[i + 1]->eventinterprestrict) PetscLogEventBegin(mglevels[i + 1]->eventinterprestrict, 0, 0, 0, 0);
89: if (matapp) MatMatInterpolate(mglevels[i + 1]->interpolate, mglevels[i]->X, &mglevels[i + 1]->X);
90: else MatInterpolate(mglevels[i + 1]->interpolate, mglevels[i]->x, mglevels[i + 1]->x);
91: if (mglevels[i + 1]->eventinterprestrict) PetscLogEventEnd(mglevels[i + 1]->eventinterprestrict, 0, 0, 0, 0);
92: }
93: if (mglevels[l - 1]->eventsmoothsolve) PetscLogEventBegin(mglevels[l - 1]->eventsmoothsolve, 0, 0, 0, 0);
94: if (matapp) {
95: KSPMatSolve(mglevels[l - 1]->smoothd, mglevels[l - 1]->B, mglevels[l - 1]->X);
96: KSPCheckSolve(mglevels[l - 1]->smoothd, pc, NULL);
97: } else {
98: KSPSolve(mglevels[l - 1]->smoothd, mglevels[l - 1]->b, mglevels[l - 1]->x);
99: KSPCheckSolve(mglevels[l - 1]->smoothd, pc, mglevels[l - 1]->x);
100: }
101: if (mglevels[l - 1]->eventsmoothsolve) PetscLogEventEnd(mglevels[l - 1]->eventsmoothsolve, 0, 0, 0, 0);
102: return 0;
103: }