1: /*
  2:      Additive Multigrid V Cycle routine
  3: */
  4: #include <petsc/private/pcmgimpl.h>

  6: PetscErrorCode PCMGACycle_Private(PC pc, PC_MG_Levels **mglevels, PetscBool transpose, PetscBool matapp)
  7: {
  8:   PetscInt i, l = mglevels[0]->levels;

 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) PetscCall(MatDuplicate(mglevels[i]->B, MAT_DO_NOT_COPY_VALUES, &mglevels[i]->X));
 27:       else PetscCall(MatZeroEntries(mglevels[i]->X));
 28:     } else PetscCall(VecZeroEntries(mglevels[i]->x));
 29:     if (mglevels[i]->eventsmoothsolve) PetscCall(PetscLogEventBegin(mglevels[i]->eventsmoothsolve, 0, 0, 0, 0));
 30:     if (!transpose) {
 31:       if (matapp) {
 32:         PetscCall(KSPMatSolve(mglevels[i]->smoothd, mglevels[i]->B, mglevels[i]->X));
 33:         PetscCall(KSPCheckSolve(mglevels[i]->smoothd, pc, NULL));
 34:       } else {
 35:         PetscCall(KSPSolve(mglevels[i]->smoothd, mglevels[i]->b, mglevels[i]->x));
 36:         PetscCall(KSPCheckSolve(mglevels[i]->smoothd, pc, mglevels[i]->x));
 37:       }
 38:     } else {
 39:       PetscCheck(!matapp, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Not supported");
 40:       PetscCall(KSPSolveTranspose(mglevels[i]->smoothu, mglevels[i]->b, mglevels[i]->x));
 41:       PetscCall(KSPCheckSolve(mglevels[i]->smoothu, pc, mglevels[i]->x));
 42:     }
 43:     if (mglevels[i]->eventsmoothsolve) PetscCall(PetscLogEventEnd(mglevels[i]->eventsmoothsolve, 0, 0, 0, 0));
 44:   }
 45:   for (i = 1; i < l; i++) {
 46:     if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
 47:     if (!transpose) {
 48:       if (matapp) PetscCall(MatMatInterpolateAdd(mglevels[i]->interpolate, mglevels[i - 1]->X, mglevels[i]->X, &mglevels[i]->X));
 49:       else PetscCall(MatInterpolateAdd(mglevels[i]->interpolate, mglevels[i - 1]->x, mglevels[i]->x, mglevels[i]->x));
 50:     } else {
 51:       if (matapp) PetscCall(MatMatInterpolateAdd(mglevels[i]->restrct, mglevels[i - 1]->X, mglevels[i]->X, &mglevels[i]->X));
 52:       else PetscCall(MatInterpolateAdd(mglevels[i]->restrct, mglevels[i - 1]->x, mglevels[i]->x, mglevels[i]->x));
 53:     }
 54:     if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
 55:   }
 56:   PetscFunctionReturn(PETSC_SUCCESS);
 57: }