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

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

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

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