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