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: {
9: PetscInt i,l = mglevels[0]->levels;
12: if (!transpose) {
13: /* restrict the RHS through all levels to coarsest. */
14: for (i=l-1; i>0; i--) {
15: if (mglevels[i]->eventinterprestrict) {PetscLogEventBegin(mglevels[i]->eventinterprestrict,0,0,0,0);}
16: if (matapp) { MatMatRestrict(mglevels[i]->restrct,mglevels[i]->B,&mglevels[i-1]->B); }
17: else { MatRestrict(mglevels[i]->restrct,mglevels[i]->b,mglevels[i-1]->b); }
18: if (mglevels[i]->eventinterprestrict) {PetscLogEventEnd(mglevels[i]->eventinterprestrict,0,0,0,0);}
19: }
21: /* work our way up through the levels */
22: if (matapp) {
23: if (!mglevels[0]->X) {
24: MatDuplicate(mglevels[0]->B,MAT_DO_NOT_COPY_VALUES,&mglevels[0]->X);
25: } else {
26: MatZeroEntries(mglevels[0]->X);
27: }
28: } else {
29: VecZeroEntries(mglevels[0]->x);
30: }
31: for (i=0; i<l-1; i++) {
32: PCMGMCycle_Private(pc,&mglevels[i],transpose,matapp,NULL);
33: if (mglevels[i+1]->eventinterprestrict) {PetscLogEventBegin(mglevels[i+1]->eventinterprestrict,0,0,0,0);}
34: if (matapp) { MatMatInterpolate(mglevels[i+1]->interpolate,mglevels[i]->X,&mglevels[i+1]->X); }
35: else { MatInterpolate(mglevels[i+1]->interpolate,mglevels[i]->x,mglevels[i+1]->x); }
36: if (mglevels[i+1]->eventinterprestrict) {PetscLogEventEnd(mglevels[i+1]->eventinterprestrict,0,0,0,0);}
37: }
38: PCMGMCycle_Private(pc,&mglevels[l-1],transpose,matapp,NULL);
39: } else {
40: PCMGMCycle_Private(pc,&mglevels[l-1],transpose,matapp,NULL);
41: for (i=l-2; i>=0; i--) {
42: if (mglevels[i+1]->eventinterprestrict) {PetscLogEventBegin(mglevels[i+1]->eventinterprestrict,0,0,0,0);}
43: if (matapp) { MatMatRestrict(mglevels[i+1]->interpolate,mglevels[i+1]->X,&mglevels[i]->X); }
44: else { MatRestrict(mglevels[i+1]->interpolate,mglevels[i+1]->x,mglevels[i]->x); }
45: if (mglevels[i+1]->eventinterprestrict) {PetscLogEventEnd(mglevels[i+1]->eventinterprestrict,0,0,0,0);}
46: PCMGMCycle_Private(pc,&mglevels[i],transpose,matapp,NULL);
47: }
48: for (i=1; i<l; i++) {
49: if (mglevels[i]->eventinterprestrict) {PetscLogEventBegin(mglevels[i]->eventinterprestrict,0,0,0,0);}
50: if (matapp) { MatMatInterpolate(mglevels[i]->restrct,mglevels[i-1]->B,&mglevels[i]->B); }
51: else { MatInterpolate(mglevels[i]->restrct,mglevels[i-1]->b,mglevels[i]->b); }
52: if (mglevels[i]->eventinterprestrict) {PetscLogEventEnd(mglevels[i]->eventinterprestrict,0,0,0,0);}
53: }
54: }
55: return(0);
56: }
58: PetscErrorCode PCMGKCycle_Private(PC pc,PC_MG_Levels **mglevels,PetscBool transpose,PetscBool matapp)
59: {
61: PetscInt i,l = mglevels[0]->levels;
64: /* restrict the RHS through all levels to coarsest. */
65: for (i=l-1; i>0; i--) {
66: if (mglevels[i]->eventinterprestrict) {PetscLogEventBegin(mglevels[i]->eventinterprestrict,0,0,0,0);}
67: if (matapp) { MatMatRestrict(mglevels[i]->restrct,mglevels[i]->B,&mglevels[i-1]->B); }
68: else { MatRestrict(mglevels[i]->restrct,mglevels[i]->b,mglevels[i-1]->b); }
69: if (mglevels[i]->eventinterprestrict) {PetscLogEventEnd(mglevels[i]->eventinterprestrict,0,0,0,0);}
70: }
72: /* work our way up through the levels */
73: if (matapp) {
74: if (!mglevels[0]->X) {
75: MatDuplicate(mglevels[0]->B,MAT_DO_NOT_COPY_VALUES,&mglevels[0]->X);
76: } else {
77: MatZeroEntries(mglevels[0]->X);
78: }
79: } else {
80: VecZeroEntries(mglevels[0]->x);
81: }
82: for (i=0; i<l-1; i++) {
83: if (mglevels[i]->eventsmoothsolve) {PetscLogEventBegin(mglevels[i]->eventsmoothsolve,0,0,0,0);}
84: if (matapp) {
85: KSPMatSolve(mglevels[i]->smoothd,mglevels[i]->B,mglevels[i]->X);
86: KSPCheckSolve(mglevels[i]->smoothd,pc,NULL);
87: } else {
88: KSPSolve(mglevels[i]->smoothd,mglevels[i]->b,mglevels[i]->x);
89: KSPCheckSolve(mglevels[i]->smoothd,pc,mglevels[i]->x);
90: }
91: if (mglevels[i]->eventsmoothsolve) {PetscLogEventEnd(mglevels[i]->eventsmoothsolve,0,0,0,0);}
92: if (mglevels[i+1]->eventinterprestrict) {PetscLogEventBegin(mglevels[i+1]->eventinterprestrict,0,0,0,0);}
93: if (matapp) { MatMatInterpolate(mglevels[i+1]->interpolate,mglevels[i]->X,&mglevels[i+1]->X); }
94: else { MatInterpolate(mglevels[i+1]->interpolate,mglevels[i]->x,mglevels[i+1]->x); }
95: if (mglevels[i+1]->eventinterprestrict) {PetscLogEventEnd(mglevels[i+1]->eventinterprestrict,0,0,0,0);}
96: }
97: if (mglevels[l-1]->eventsmoothsolve) {PetscLogEventBegin(mglevels[l-1]->eventsmoothsolve,0,0,0,0);}
98: if (matapp) {
99: KSPMatSolve(mglevels[l-1]->smoothd,mglevels[l-1]->B,mglevels[l-1]->X);
100: KSPCheckSolve(mglevels[l-1]->smoothd,pc,NULL);
101: } else {
102: KSPSolve(mglevels[l-1]->smoothd,mglevels[l-1]->b,mglevels[l-1]->x);
103: KSPCheckSolve(mglevels[l-1]->smoothd,pc,mglevels[l-1]->x);
104: }
105: if (mglevels[l-1]->eventsmoothsolve) {PetscLogEventEnd(mglevels[l-1]->eventsmoothsolve,0,0,0,0);}
106: return(0);
107: }