1 /*
2 Additive Multigrid V Cycle routine
3 */
4 #include <petsc/private/pcmgimpl.h>
5
PCMGACycle_Private(PC pc,PC_MG_Levels ** mglevels,PetscBool transpose,PetscBool matapp)6 PetscErrorCode PCMGACycle_Private(PC pc, PC_MG_Levels **mglevels, PetscBool transpose, PetscBool matapp)
7 {
8 PetscInt i, l = mglevels[0]->levels;
9
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 }
63