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