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:   /* compute RHS on each level */
 12:   for (i=l-1; i>0; i--) {
 13:     if (mglevels[i]->eventinterprestrict) PetscLogEventBegin(mglevels[i]->eventinterprestrict,0,0,0,0);
 14:     if (!transpose) {
 15:       if (matapp) MatMatRestrict(mglevels[i]->restrct,mglevels[i]->B,&mglevels[i-1]->B);
 16:       else MatRestrict(mglevels[i]->restrct,mglevels[i]->b,mglevels[i-1]->b);
 17:     } else {
 18:       if (matapp) MatMatRestrict(mglevels[i]->interpolate,mglevels[i]->B,&mglevels[i-1]->B);
 19:       else MatRestrict(mglevels[i]->interpolate,mglevels[i]->b,mglevels[i-1]->b);
 20:     }
 21:     if (mglevels[i]->eventinterprestrict) 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:         MatDuplicate(mglevels[i]->B,MAT_DO_NOT_COPY_VALUES,&mglevels[i]->X);
 28:       } else {
 29:         MatZeroEntries(mglevels[i]->X);
 30:       }
 31:     } else {
 32:       VecZeroEntries(mglevels[i]->x);
 33:     }
 34:     if (mglevels[i]->eventsmoothsolve) PetscLogEventBegin(mglevels[i]->eventsmoothsolve,0,0,0,0);
 35:     if (!transpose) {
 36:       if (matapp) {
 37:         KSPMatSolve(mglevels[i]->smoothd,mglevels[i]->B,mglevels[i]->X);
 38:         KSPCheckSolve(mglevels[i]->smoothd,pc,NULL);
 39:       } else {
 40:         KSPSolve(mglevels[i]->smoothd,mglevels[i]->b,mglevels[i]->x);
 41:         KSPCheckSolve(mglevels[i]->smoothd,pc,mglevels[i]->x);
 42:       }
 43:     } else {
 45:       KSPSolveTranspose(mglevels[i]->smoothu,mglevels[i]->b,mglevels[i]->x);
 46:       KSPCheckSolve(mglevels[i]->smoothu,pc,mglevels[i]->x);
 47:     }
 48:     if (mglevels[i]->eventsmoothsolve) PetscLogEventEnd(mglevels[i]->eventsmoothsolve,0,0,0,0);
 49:   }
 50:   for (i=1; i<l; i++) {
 51:     if (mglevels[i]->eventinterprestrict) PetscLogEventBegin(mglevels[i]->eventinterprestrict,0,0,0,0);
 52:     if (!transpose) {
 53:       if (matapp) MatMatInterpolateAdd(mglevels[i]->interpolate,mglevels[i-1]->X,mglevels[i]->X,&mglevels[i]->X);
 54:       else MatInterpolateAdd(mglevels[i]->interpolate,mglevels[i-1]->x,mglevels[i]->x,mglevels[i]->x);
 55:     } else {
 56:       if (matapp) MatMatInterpolateAdd(mglevels[i]->restrct,mglevels[i-1]->X,mglevels[i]->X,&mglevels[i]->X);
 57:       else MatInterpolateAdd(mglevels[i]->restrct,mglevels[i-1]->x,mglevels[i]->x,mglevels[i]->x);
 58:     }
 59:     if (mglevels[i]->eventinterprestrict) PetscLogEventEnd(mglevels[i]->eventinterprestrict,0,0,0,0);
 60:   }
 61:   return 0;
 62: }