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: }