Actual source code: smg.c

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