1 2 /* 3 Additive Multigrid V Cycle routine 4 */ 5 #include <petsc/private/pcmgimpl.h> 6 7 PetscErrorCode PCMGACycle_Private(PC pc,PC_MG_Levels **mglevels,PetscBool transpose,PetscBool matapp) 8 { 9 PetscErrorCode ierr; 10 PetscInt i,l = mglevels[0]->levels; 11 12 PetscFunctionBegin; 13 /* compute RHS on each level */ 14 for (i=l-1; i>0; i--) { 15 if (mglevels[i]->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 16 if (!transpose) { 17 if (matapp) { ierr = MatMatRestrict(mglevels[i]->restrct,mglevels[i]->B,&mglevels[i-1]->B);CHKERRQ(ierr); } 18 else { ierr = MatRestrict(mglevels[i]->restrct,mglevels[i]->b,mglevels[i-1]->b);CHKERRQ(ierr); } 19 } else { 20 if (matapp) { ierr = MatMatRestrict(mglevels[i]->interpolate,mglevels[i]->B,&mglevels[i-1]->B);CHKERRQ(ierr); } 21 else { ierr = MatRestrict(mglevels[i]->interpolate,mglevels[i]->b,mglevels[i-1]->b);CHKERRQ(ierr); } 22 } 23 if (mglevels[i]->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 24 } 25 /* solve separately on each level */ 26 for (i=0; i<l; i++) { 27 if (matapp) { 28 if (!mglevels[i]->X) { 29 ierr = MatDuplicate(mglevels[i]->B,MAT_DO_NOT_COPY_VALUES,&mglevels[i]->X);CHKERRQ(ierr); 30 } else { 31 ierr = MatZeroEntries(mglevels[i]->X);CHKERRQ(ierr); 32 } 33 } else { 34 ierr = VecZeroEntries(mglevels[i]->x);CHKERRQ(ierr); 35 } 36 if (mglevels[i]->eventsmoothsolve) {ierr = PetscLogEventBegin(mglevels[i]->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);} 37 if (!transpose) { 38 if (matapp) { 39 ierr = KSPMatSolve(mglevels[i]->smoothd,mglevels[i]->B,mglevels[i]->X);CHKERRQ(ierr); 40 ierr = KSPCheckSolve(mglevels[i]->smoothd,pc,NULL);CHKERRQ(ierr); 41 } else { 42 ierr = KSPSolve(mglevels[i]->smoothd,mglevels[i]->b,mglevels[i]->x);CHKERRQ(ierr); 43 ierr = KSPCheckSolve(mglevels[i]->smoothd,pc,mglevels[i]->x);CHKERRQ(ierr); 44 } 45 } else { 46 if (matapp) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Not supported"); 47 ierr = KSPSolveTranspose(mglevels[i]->smoothu,mglevels[i]->b,mglevels[i]->x);CHKERRQ(ierr); 48 ierr = KSPCheckSolve(mglevels[i]->smoothu,pc,mglevels[i]->x);CHKERRQ(ierr); 49 } 50 if (mglevels[i]->eventsmoothsolve) {ierr = PetscLogEventEnd(mglevels[i]->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);} 51 } 52 for (i=1; i<l; i++) { 53 if (mglevels[i]->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 54 if (!transpose) { 55 if (matapp) { ierr = MatMatInterpolateAdd(mglevels[i]->interpolate,mglevels[i-1]->X,mglevels[i]->X,&mglevels[i]->X);CHKERRQ(ierr); } 56 else { ierr = MatInterpolateAdd(mglevels[i]->interpolate,mglevels[i-1]->x,mglevels[i]->x,mglevels[i]->x);CHKERRQ(ierr); } 57 } else { 58 if (matapp) { ierr = MatMatInterpolateAdd(mglevels[i]->restrct,mglevels[i-1]->X,mglevels[i]->X,&mglevels[i]->X);CHKERRQ(ierr); } 59 else { ierr = MatInterpolateAdd(mglevels[i]->restrct,mglevels[i-1]->x,mglevels[i]->x,mglevels[i]->x);CHKERRQ(ierr); } 60 } 61 if (mglevels[i]->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 62 } 63 PetscFunctionReturn(0); 64 } 65