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: /* compute RHS on each level */
12: for (i = l - 1; i > 0; i--) {
13: if (mglevels[i]->eventinterprestrict) PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0);
14: if (!transpose) {
15: if (matapp) MatMatRestrict(mglevels[i]->restrct, mglevels[i]->B, &mglevels[i - 1]->B);
16: else MatRestrict(mglevels[i]->restrct, mglevels[i]->b, mglevels[i - 1]->b);
17: } else {
18: if (matapp) MatMatRestrict(mglevels[i]->interpolate, mglevels[i]->B, &mglevels[i - 1]->B);
19: else MatRestrict(mglevels[i]->interpolate, mglevels[i]->b, mglevels[i - 1]->b);
20: }
21: if (mglevels[i]->eventinterprestrict) PetscLogEventEnd(mglevels[i]->eventinterprestrict, 0, 0, 0, 0);
22: }
23: /* solve separately on each level */
24: for (i = 0; i < l; i++) {
25: if (matapp) {
26: if (!mglevels[i]->X) {
27: MatDuplicate(mglevels[i]->B, MAT_DO_NOT_COPY_VALUES, &mglevels[i]->X);
28: } else {
29: MatZeroEntries(mglevels[i]->X);
30: }
31: } else {
32: VecZeroEntries(mglevels[i]->x);
33: }
34: if (mglevels[i]->eventsmoothsolve) PetscLogEventBegin(mglevels[i]->eventsmoothsolve, 0, 0, 0, 0);
35: if (!transpose) {
36: if (matapp) {
37: KSPMatSolve(mglevels[i]->smoothd, mglevels[i]->B, mglevels[i]->X);
38: KSPCheckSolve(mglevels[i]->smoothd, pc, NULL);
39: } else {
40: KSPSolve(mglevels[i]->smoothd, mglevels[i]->b, mglevels[i]->x);
41: KSPCheckSolve(mglevels[i]->smoothd, pc, mglevels[i]->x);
42: }
43: } else {
45: KSPSolveTranspose(mglevels[i]->smoothu, mglevels[i]->b, mglevels[i]->x);
46: KSPCheckSolve(mglevels[i]->smoothu, pc, mglevels[i]->x);
47: }
48: if (mglevels[i]->eventsmoothsolve) PetscLogEventEnd(mglevels[i]->eventsmoothsolve, 0, 0, 0, 0);
49: }
50: for (i = 1; i < l; i++) {
51: if (mglevels[i]->eventinterprestrict) PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0);
52: if (!transpose) {
53: if (matapp) MatMatInterpolateAdd(mglevels[i]->interpolate, mglevels[i - 1]->X, mglevels[i]->X, &mglevels[i]->X);
54: else MatInterpolateAdd(mglevels[i]->interpolate, mglevels[i - 1]->x, mglevels[i]->x, mglevels[i]->x);
55: } else {
56: if (matapp) MatMatInterpolateAdd(mglevels[i]->restrct, mglevels[i - 1]->X, mglevels[i]->X, &mglevels[i]->X);
57: else MatInterpolateAdd(mglevels[i]->restrct, mglevels[i - 1]->x, mglevels[i]->x, mglevels[i]->x);
58: }
59: if (mglevels[i]->eventinterprestrict) PetscLogEventEnd(mglevels[i]->eventinterprestrict, 0, 0, 0, 0);
60: }
61: return 0;
62: }