1: /*
  2:      Full multigrid using either additive or multiplicative V or W cycle
  3: */
  4: #include <petsc/private/pcmgimpl.h>

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

 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:     }

 20:     /* work our way up through the levels */
 21:     if (matapp) {
 22:       if (!mglevels[0]->X) PetscCall(MatDuplicate(mglevels[0]->B, MAT_DO_NOT_COPY_VALUES, &mglevels[0]->X));
 23:       else PetscCall(MatZeroEntries(mglevels[0]->X));
 24:     } else PetscCall(VecZeroEntries(mglevels[0]->x));
 25:     for (i = 0; i < l - 1; i++) {
 26:       PetscCall(PCMGMCycle_Private(pc, &mglevels[i], transpose, matapp, NULL));
 27:       if (mglevels[i + 1]->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels[i + 1]->eventinterprestrict, 0, 0, 0, 0));
 28:       if (matapp) PetscCall(MatMatInterpolate(mglevels[i + 1]->interpolate, mglevels[i]->X, &mglevels[i + 1]->X));
 29:       else PetscCall(MatInterpolate(mglevels[i + 1]->interpolate, mglevels[i]->x, mglevels[i + 1]->x));
 30:       if (mglevels[i + 1]->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels[i + 1]->eventinterprestrict, 0, 0, 0, 0));
 31:     }
 32:     PetscCall(PCMGMCycle_Private(pc, &mglevels[l - 1], transpose, matapp, NULL));
 33:   } else {
 34:     PetscCall(PCMGMCycle_Private(pc, &mglevels[l - 1], transpose, matapp, NULL));
 35:     for (i = l - 2; i >= 0; i--) {
 36:       if (mglevels[i + 1]->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels[i + 1]->eventinterprestrict, 0, 0, 0, 0));
 37:       if (matapp) PetscCall(MatMatRestrict(mglevels[i + 1]->interpolate, mglevels[i + 1]->X, &mglevels[i]->X));
 38:       else PetscCall(MatRestrict(mglevels[i + 1]->interpolate, mglevels[i + 1]->x, mglevels[i]->x));
 39:       if (mglevels[i + 1]->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels[i + 1]->eventinterprestrict, 0, 0, 0, 0));
 40:       PetscCall(PCMGMCycle_Private(pc, &mglevels[i], transpose, matapp, NULL));
 41:     }
 42:     for (i = 1; i < l; i++) {
 43:       if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
 44:       if (matapp) PetscCall(MatMatInterpolate(mglevels[i]->restrct, mglevels[i - 1]->B, &mglevels[i]->B));
 45:       else PetscCall(MatInterpolate(mglevels[i]->restrct, mglevels[i - 1]->b, mglevels[i]->b));
 46:       if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
 47:     }
 48:   }
 49:   PetscFunctionReturn(PETSC_SUCCESS);
 50: }

 52: PetscErrorCode PCMGKCycle_Private(PC pc, PC_MG_Levels **mglevels, PetscBool transpose, PetscBool matapp)
 53: {
 54:   PetscInt i, l = mglevels[0]->levels;

 56:   PetscFunctionBegin;
 57:   /* restrict the RHS through all levels to coarsest. */
 58:   for (i = l - 1; i > 0; i--) {
 59:     if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
 60:     if (matapp) PetscCall(MatMatRestrict(mglevels[i]->restrct, mglevels[i]->B, &mglevels[i - 1]->B));
 61:     else PetscCall(MatRestrict(mglevels[i]->restrct, mglevels[i]->b, mglevels[i - 1]->b));
 62:     if (mglevels[i]->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels[i]->eventinterprestrict, 0, 0, 0, 0));
 63:   }

 65:   /* work our way up through the levels */
 66:   if (matapp) {
 67:     if (!mglevels[0]->X) PetscCall(MatDuplicate(mglevels[0]->B, MAT_DO_NOT_COPY_VALUES, &mglevels[0]->X));
 68:     else {
 69:       PetscCall(MatZeroEntries(mglevels[0]->X));
 70:     }
 71:   } else {
 72:     PetscCall(VecZeroEntries(mglevels[0]->x));
 73:   }
 74:   for (i = 0; i < l - 1; i++) {
 75:     if (mglevels[i]->eventsmoothsolve) PetscCall(PetscLogEventBegin(mglevels[i]->eventsmoothsolve, 0, 0, 0, 0));
 76:     if (matapp) {
 77:       PetscCall(KSPMatSolve(mglevels[i]->smoothd, mglevels[i]->B, mglevels[i]->X));
 78:       PetscCall(KSPCheckSolve(mglevels[i]->smoothd, pc, NULL));
 79:     } else {
 80:       PetscCall(KSPSolve(mglevels[i]->smoothd, mglevels[i]->b, mglevels[i]->x));
 81:       PetscCall(KSPCheckSolve(mglevels[i]->smoothd, pc, mglevels[i]->x));
 82:     }
 83:     if (mglevels[i]->eventsmoothsolve) PetscCall(PetscLogEventEnd(mglevels[i]->eventsmoothsolve, 0, 0, 0, 0));
 84:     if (mglevels[i + 1]->eventinterprestrict) PetscCall(PetscLogEventBegin(mglevels[i + 1]->eventinterprestrict, 0, 0, 0, 0));
 85:     if (matapp) PetscCall(MatMatInterpolate(mglevels[i + 1]->interpolate, mglevels[i]->X, &mglevels[i + 1]->X));
 86:     else PetscCall(MatInterpolate(mglevels[i + 1]->interpolate, mglevels[i]->x, mglevels[i + 1]->x));
 87:     if (mglevels[i + 1]->eventinterprestrict) PetscCall(PetscLogEventEnd(mglevels[i + 1]->eventinterprestrict, 0, 0, 0, 0));
 88:   }
 89:   if (mglevels[l - 1]->eventsmoothsolve) PetscCall(PetscLogEventBegin(mglevels[l - 1]->eventsmoothsolve, 0, 0, 0, 0));
 90:   if (matapp) {
 91:     PetscCall(KSPMatSolve(mglevels[l - 1]->smoothd, mglevels[l - 1]->B, mglevels[l - 1]->X));
 92:     PetscCall(KSPCheckSolve(mglevels[l - 1]->smoothd, pc, NULL));
 93:   } else {
 94:     PetscCall(KSPSolve(mglevels[l - 1]->smoothd, mglevels[l - 1]->b, mglevels[l - 1]->x));
 95:     PetscCall(KSPCheckSolve(mglevels[l - 1]->smoothd, pc, mglevels[l - 1]->x));
 96:   }
 97:   if (mglevels[l - 1]->eventsmoothsolve) PetscCall(PetscLogEventEnd(mglevels[l - 1]->eventsmoothsolve, 0, 0, 0, 0));
 98:   PetscFunctionReturn(PETSC_SUCCESS);
 99: }