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