1 /* 2 Additive Multigrid V Cycle routine 3 */ 4 #include <petsc/private/pcmgimpl.h> 5 6 PetscErrorCode PCMGACycle_Private(PC pc, PC_MG_Levels **mglevels, PetscBool transpose, PetscBool matapp) 7 { 8 PetscInt i, l = mglevels[0]->levels; 9 10 PetscFunctionBegin; 11 /* compute RHS on each level */ 12 for (i = l - 1; i > 0; i--) { 13 if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0)); 14 if (!transpose) { 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 } else { 18 if (matapp) PetscCall(MatMatRestrict(mglevels[i]->interpolate, mglevels[i]->B, &mglevels[i - 1]->B)); 19 else PetscCall(MatRestrict(mglevels[i]->interpolate, mglevels[i]->b, mglevels[i - 1]->b)); 20 } 21 if (mglevels[i]->eventinterprestrict) PetscCall(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 PetscCall(MatDuplicate(mglevels[i]->B, MAT_DO_NOT_COPY_VALUES, &mglevels[i]->X)); 28 } else { 29 PetscCall(MatZeroEntries(mglevels[i]->X)); 30 } 31 } else { 32 PetscCall(VecZeroEntries(mglevels[i]->x)); 33 } 34 if (mglevels[i]->eventsmoothsolve) PetscCall(PetscLogEventBegin(mglevels[i]->eventsmoothsolve, 0, 0, 0, 0)); 35 if (!transpose) { 36 if (matapp) { 37 PetscCall(KSPMatSolve(mglevels[i]->smoothd, mglevels[i]->B, mglevels[i]->X)); 38 PetscCall(KSPCheckSolve(mglevels[i]->smoothd, pc, NULL)); 39 } else { 40 PetscCall(KSPSolve(mglevels[i]->smoothd, mglevels[i]->b, mglevels[i]->x)); 41 PetscCall(KSPCheckSolve(mglevels[i]->smoothd, pc, mglevels[i]->x)); 42 } 43 } else { 44 PetscCheck(!matapp, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not supported"); 45 PetscCall(KSPSolveTranspose(mglevels[i]->smoothu, mglevels[i]->b, mglevels[i]->x)); 46 PetscCall(KSPCheckSolve(mglevels[i]->smoothu, pc, mglevels[i]->x)); 47 } 48 if (mglevels[i]->eventsmoothsolve) PetscCall(PetscLogEventEnd(mglevels[i]->eventsmoothsolve, 0, 0, 0, 0)); 49 } 50 for (i = 1; i < l; i++) { 51 if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0)); 52 if (!transpose) { 53 if (matapp) PetscCall(MatMatInterpolateAdd(mglevels[i]->interpolate, mglevels[i - 1]->X, mglevels[i]->X, &mglevels[i]->X)); 54 else PetscCall(MatInterpolateAdd(mglevels[i]->interpolate, mglevels[i - 1]->x, mglevels[i]->x, mglevels[i]->x)); 55 } else { 56 if (matapp) PetscCall(MatMatInterpolateAdd(mglevels[i]->restrct, mglevels[i - 1]->X, mglevels[i]->X, &mglevels[i]->X)); 57 else PetscCall(MatInterpolateAdd(mglevels[i]->restrct, mglevels[i - 1]->x, mglevels[i]->x, mglevels[i]->x)); 58 } 59 if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels[i]->eventinterprestrict, 0, 0, 0, 0)); 60 } 61 PetscFunctionReturn(PETSC_SUCCESS); 62 } 63