1 /* 2 Full multigrid using either additive or multiplicative V or W cycle 3 */ 4 #include <petsc/private/pcmgimpl.h> 5 6 PetscErrorCode PCMGFCycle_Private(PC pc,PC_MG_Levels **mglevels,PetscBool transpose,PetscBool matapp) 7 { 8 PetscErrorCode ierr; 9 PetscInt i,l = mglevels[0]->levels; 10 11 PetscFunctionBegin; 12 if (!transpose) { 13 /* restrict the RHS through all levels to coarsest. */ 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 (matapp) { ierr = MatMatRestrict(mglevels[i]->restrct,mglevels[i]->B,&mglevels[i-1]->B);CHKERRQ(ierr); } 17 else { ierr = MatRestrict(mglevels[i]->restrct,mglevels[i]->b,mglevels[i-1]->b);CHKERRQ(ierr); } 18 if (mglevels[i]->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 19 } 20 21 /* work our way up through the levels */ 22 if (matapp) { 23 if (!mglevels[0]->X) { 24 ierr = MatDuplicate(mglevels[0]->B,MAT_DO_NOT_COPY_VALUES,&mglevels[0]->X);CHKERRQ(ierr); 25 } else { 26 ierr = MatZeroEntries(mglevels[0]->X);CHKERRQ(ierr); 27 } 28 } else { 29 ierr = VecZeroEntries(mglevels[0]->x);CHKERRQ(ierr); 30 } 31 for (i=0; i<l-1; i++) { 32 ierr = PCMGMCycle_Private(pc,&mglevels[i],transpose,matapp,NULL);CHKERRQ(ierr); 33 if (mglevels[i+1]->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels[i+1]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 34 if (matapp) { ierr = MatMatInterpolate(mglevels[i+1]->interpolate,mglevels[i]->X,&mglevels[i+1]->X);CHKERRQ(ierr); } 35 else { ierr = MatInterpolate(mglevels[i+1]->interpolate,mglevels[i]->x,mglevels[i+1]->x);CHKERRQ(ierr); } 36 if (mglevels[i+1]->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels[i+1]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 37 } 38 ierr = PCMGMCycle_Private(pc,&mglevels[l-1],transpose,matapp,NULL);CHKERRQ(ierr); 39 } else { 40 ierr = PCMGMCycle_Private(pc,&mglevels[l-1],transpose,matapp,NULL);CHKERRQ(ierr); 41 for (i=l-2; i>=0; i--) { 42 if (mglevels[i+1]->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels[i+1]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 43 if (matapp) { ierr = MatMatRestrict(mglevels[i+1]->interpolate,mglevels[i+1]->X,&mglevels[i]->X);CHKERRQ(ierr); } 44 else { ierr = MatRestrict(mglevels[i+1]->interpolate,mglevels[i+1]->x,mglevels[i]->x);CHKERRQ(ierr); } 45 if (mglevels[i+1]->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels[i+1]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 46 ierr = PCMGMCycle_Private(pc,&mglevels[i],transpose,matapp,NULL);CHKERRQ(ierr); 47 } 48 for (i=1; i<l; i++) { 49 if (mglevels[i]->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 50 if (matapp) { ierr = MatMatInterpolate(mglevels[i]->restrct,mglevels[i-1]->B,&mglevels[i]->B);CHKERRQ(ierr); } 51 else { ierr = MatInterpolate(mglevels[i]->restrct,mglevels[i-1]->b,mglevels[i]->b);CHKERRQ(ierr); } 52 if (mglevels[i]->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 53 } 54 } 55 PetscFunctionReturn(0); 56 } 57 58 PetscErrorCode PCMGKCycle_Private(PC pc,PC_MG_Levels **mglevels,PetscBool transpose,PetscBool matapp) 59 { 60 PetscErrorCode ierr; 61 PetscInt i,l = mglevels[0]->levels; 62 63 PetscFunctionBegin; 64 /* restrict the RHS through all levels to coarsest. */ 65 for (i=l-1; i>0; i--) { 66 if (mglevels[i]->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 67 if (matapp) { ierr = MatMatRestrict(mglevels[i]->restrct,mglevels[i]->B,&mglevels[i-1]->B);CHKERRQ(ierr); } 68 else { ierr = MatRestrict(mglevels[i]->restrct,mglevels[i]->b,mglevels[i-1]->b);CHKERRQ(ierr); } 69 if (mglevels[i]->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels[i]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 70 } 71 72 /* work our way up through the levels */ 73 if (matapp) { 74 if (!mglevels[0]->X) { 75 ierr = MatDuplicate(mglevels[0]->B,MAT_DO_NOT_COPY_VALUES,&mglevels[0]->X);CHKERRQ(ierr); 76 } else { 77 ierr = MatZeroEntries(mglevels[0]->X);CHKERRQ(ierr); 78 } 79 } else { 80 ierr = VecZeroEntries(mglevels[0]->x);CHKERRQ(ierr); 81 } 82 for (i=0; i<l-1; i++) { 83 if (mglevels[i]->eventsmoothsolve) {ierr = PetscLogEventBegin(mglevels[i]->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);} 84 if (matapp) { 85 ierr = KSPMatSolve(mglevels[i]->smoothd,mglevels[i]->B,mglevels[i]->X);CHKERRQ(ierr); 86 ierr = KSPCheckSolve(mglevels[i]->smoothd,pc,NULL);CHKERRQ(ierr); 87 } else { 88 ierr = KSPSolve(mglevels[i]->smoothd,mglevels[i]->b,mglevels[i]->x);CHKERRQ(ierr); 89 ierr = KSPCheckSolve(mglevels[i]->smoothd,pc,mglevels[i]->x);CHKERRQ(ierr); 90 } 91 if (mglevels[i]->eventsmoothsolve) {ierr = PetscLogEventEnd(mglevels[i]->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);} 92 if (mglevels[i+1]->eventinterprestrict) {ierr = PetscLogEventBegin(mglevels[i+1]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 93 if (matapp) { ierr = MatMatInterpolate(mglevels[i+1]->interpolate,mglevels[i]->X,&mglevels[i+1]->X);CHKERRQ(ierr); } 94 else { ierr = MatInterpolate(mglevels[i+1]->interpolate,mglevels[i]->x,mglevels[i+1]->x);CHKERRQ(ierr); } 95 if (mglevels[i+1]->eventinterprestrict) {ierr = PetscLogEventEnd(mglevels[i+1]->eventinterprestrict,0,0,0,0);CHKERRQ(ierr);} 96 } 97 if (mglevels[l-1]->eventsmoothsolve) {ierr = PetscLogEventBegin(mglevels[l-1]->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);} 98 if (matapp) { 99 ierr = KSPMatSolve(mglevels[l-1]->smoothd,mglevels[l-1]->B,mglevels[l-1]->X);CHKERRQ(ierr); 100 ierr = KSPCheckSolve(mglevels[l-1]->smoothd,pc,NULL);CHKERRQ(ierr); 101 } else { 102 ierr = KSPSolve(mglevels[l-1]->smoothd,mglevels[l-1]->b,mglevels[l-1]->x);CHKERRQ(ierr); 103 ierr = KSPCheckSolve(mglevels[l-1]->smoothd,pc,mglevels[l-1]->x);CHKERRQ(ierr); 104 } 105 if (mglevels[l-1]->eventsmoothsolve) {ierr = PetscLogEventEnd(mglevels[l-1]->eventsmoothsolve,0,0,0,0);CHKERRQ(ierr);} 106 PetscFunctionReturn(0); 107 } 108