Actual source code: fmg.c

  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:   if (!transpose) {
 11:     /* restrict the RHS through all levels to coarsest. */
 12:     for (i = l - 1; i > 0; i--) {
 13:       if (mglevels[i]->eventinterprestrict) PetscLogEventBegin(mglevels[i]->eventinterprestrict, 0, 0, 0, 0);
 14:       if (matapp) MatMatRestrict(mglevels[i]->restrct, mglevels[i]->B, &mglevels[i - 1]->B);
 15:       else MatRestrict(mglevels[i]->restrct, mglevels[i]->b, mglevels[i - 1]->b);
 16:       if (mglevels[i]->eventinterprestrict) PetscLogEventEnd(mglevels[i]->eventinterprestrict, 0, 0, 0, 0);
 17:     }

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

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

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

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