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