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: {
 10:   PetscInt       i,l = mglevels[0]->levels;

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