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