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