Actual source code: smg.c


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

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

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