1 /* 2 GAMG geometric-algebric multiogrid PC - Mark Adams 2011 3 */ 4 #include <private/pcimpl.h> /*I "petscpc.h" I*/ 5 #include <../src/ksp/pc/impls/mg/mgimpl.h> /*I "petscpcmg.h" I*/ 6 #include <../src/mat/impls/aij/seq/aij.h> 7 #include <../src/mat/impls/aij/mpi/mpiaij.h> 8 9 /* Private context for the GAMG preconditioner */ 10 typedef struct gamg_TAG { 11 PetscInt m_dim; 12 PetscInt m_Nlevels; 13 PetscInt m_data_sz; 14 PetscReal *m_data; /* blocked vector of vertex data on fine grid (coordinates) */ 15 } PC_GAMG; 16 17 /* -----------------------------------------------------------------------------*/ 18 #undef __FUNCT__ 19 #define __FUNCT__ "PCReset_GAMG" 20 PetscErrorCode PCReset_GAMG(PC pc) 21 { 22 PetscErrorCode ierr; 23 PC_MG *mg = (PC_MG*)pc->data; 24 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 25 26 PetscFunctionBegin; 27 ierr = PetscFree(pc_gamg->m_data);CHKERRQ(ierr); 28 PetscFunctionReturn(0); 29 } 30 31 /* -------------------------------------------------------------------------- */ 32 /* 33 PCSetCoordinates_GAMG 34 35 Input Parameter: 36 . pc - the preconditioner context 37 */ 38 #undef __FUNCT__ 39 #define __FUNCT__ "PCSetCoordinates_GAMG" 40 PetscErrorCode PCSetCoordinates_GAMG(PC pc, const int ndm,PetscReal *coords ) 41 { 42 PC_MG *mg = (PC_MG*)pc->data; 43 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 44 PetscErrorCode ierr; 45 PetscInt bs, my0, tt; 46 Mat mat = pc->pmat; 47 PetscInt arrsz; 48 49 PetscFunctionBegin; 50 ierr = MatGetBlockSize( mat, &bs ); CHKERRQ( ierr ); 51 ierr = MatGetOwnershipRange( mat, &my0, &tt ); CHKERRQ(ierr); 52 arrsz = (tt-my0)/bs*ndm; 53 54 // put coordinates 55 if (!pc_gamg->m_data || (pc_gamg->m_data_sz != arrsz)) { 56 ierr = PetscFree( pc_gamg->m_data ); CHKERRQ(ierr); 57 ierr = PetscMalloc(arrsz*sizeof(double), &pc_gamg->m_data ); CHKERRQ(ierr); 58 } 59 60 /* copy data in */ 61 for(tt=0;tt<arrsz;tt++){ 62 pc_gamg->m_data[tt] = coords[tt]; 63 } 64 pc_gamg->m_data_sz = arrsz; 65 pc_gamg->m_dim = ndm; 66 67 PetscFunctionReturn(0); 68 } 69 70 /* -------------------------------------------------------------------------- */ 71 /* 72 createCrsOp 73 74 Input Parameter: 75 . Amat - matrix on this fine level 76 . P_out - prolongation operator to the next level 77 . Acrs - coarse matrix that is created 78 */ 79 #undef __FUNCT__ 80 #define __FUNCT__ "createCrsOp" 81 PetscErrorCode createCrsOp( Mat Amat, Mat P_inout, Mat *Acrs ) 82 { 83 PetscErrorCode ierr; 84 85 PetscFunctionBegin; 86 Mat H; 87 ierr = MatPtAP( Amat, P_inout, MAT_INITIAL_MATRIX, 2.0, &H); CHKERRQ(ierr); 88 89 /* need to repartition H and move colums of P accordingly */ 90 91 92 *Acrs = H; 93 94 PetscFunctionReturn(0); 95 } 96 97 /* -------------------------------------------------------------------------- */ 98 /* 99 PCSetUp_GAMG - Prepares for the use of the GAMG preconditioner 100 by setting data structures and options. 101 102 Input Parameter: 103 . pc - the preconditioner context 104 105 Application Interface Routine: PCSetUp() 106 107 Notes: 108 The interface routine PCSetUp() is not usually called directly by 109 the user, but instead is called by PCApply() if necessary. 110 */ 111 extern PetscErrorCode PCSetFromOptions_MG(PC); 112 extern PetscErrorCode PCReset_MG(PC); 113 extern PetscErrorCode createProlongation( Mat, Mat *, PetscReal [], PetscReal **, const PetscInt); 114 #undef __FUNCT__ 115 #define __FUNCT__ "PCSetUp_GAMG" 116 PetscErrorCode PCSetUp_GAMG( PC pc ) 117 { 118 PetscErrorCode ierr; 119 PC_MG *mg = (PC_MG*)pc->data; 120 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 121 Mat Amat = pc->mat, Pmat = pc->pmat; 122 PetscBool isSeq, isMPI; 123 PetscInt fine_level, level, level1, M, N, bs, lidx; 124 MPI_Comm wcomm = ((PetscObject)pc)->comm; 125 126 PetscFunctionBegin; 127 if (pc->setupcalled){ 128 /* no state data in GAMG to destroy (now) */ 129 ierr = PCReset_MG(pc);CHKERRQ(ierr); 130 } 131 if (!pc_gamg->m_data) SETERRQ(wcomm,PETSC_ERR_SUP,"PCSetUp_GAMG called before PCSetCoordinates"); 132 /* setup special features of PCGAMG */ 133 ierr = PetscTypeCompare((PetscObject) Amat, MATSEQAIJ, &isSeq);CHKERRQ(ierr); 134 ierr = PetscTypeCompare((PetscObject) Amat, MATMPIAIJ, &isMPI);CHKERRQ(ierr); 135 if (isMPI) { 136 } else if (isSeq) { 137 } else SETERRQ1(wcomm,PETSC_ERR_ARG_WRONG, "Matrix type '%s' cannot be used with GAMG. GAMG can only handle AIJ matrices.",((PetscObject)Amat)->type_name); 138 139 /* GAMG requires input of fine-grid matrix. It determines nlevels. */ 140 ierr = MatGetSize( Amat, &M, &N );CHKERRQ(ierr); 141 ierr = MatGetBlockSize( Amat, &bs ); CHKERRQ(ierr); 142 if(bs!=1) SETERRQ1(wcomm,PETSC_ERR_ARG_WRONG, "GAMG only supports scalar prblems bs = '%d'.",bs); 143 144 /* Get A_i and R_i */ 145 #define GAMG_MAXLEVELS 10 146 Mat Aarr[GAMG_MAXLEVELS], Rarr[GAMG_MAXLEVELS]; PetscReal *coarse_crds = 0, *crds = pc_gamg->m_data; 147 for (level=0, Aarr[0] = Pmat; level < GAMG_MAXLEVELS-1; level++ ){ 148 ierr = MatGetSize( Aarr[level], &M, &N );CHKERRQ(ierr); 149 if( M < 100 ) { /* hard wire this for now */ 150 break; 151 } 152 level1 = level + 1; 153 ierr = createProlongation( Aarr[level], &Rarr[level1], crds, &coarse_crds, pc_gamg->m_dim ); 154 CHKERRQ(ierr); 155 if(level==0) Aarr[0] = Amat; /* use Pmat for finest level setup, but use mat for solver */ 156 ierr = createCrsOp( Aarr[level], Rarr[level1], &Aarr[level1] ); CHKERRQ(ierr); 157 ierr = PetscFree( crds ); CHKERRQ( ierr ); 158 crds = coarse_crds; 159 } 160 ierr = PetscFree( coarse_crds ); CHKERRQ( ierr ); 161 162 pc_gamg->m_data = 0; /* destroyed coordinate data */ 163 pc_gamg->m_Nlevels = level + 1; 164 fine_level = level; 165 ierr = PCMGSetLevels(pc,pc_gamg->m_Nlevels,PETSC_NULL);CHKERRQ(ierr); 166 167 /* set default smoothers */ 168 PetscReal emax = 2.0, emin; 169 for (level=1; level<=fine_level; level++){ 170 KSP smoother; PC subpc; 171 ierr = PCMGGetSmoother(pc,level,&smoother);CHKERRQ(ierr); 172 ierr = KSPSetType(smoother,KSPCHEBYCHEV);CHKERRQ(ierr); 173 emin = emax/5.0; 174 ierr = KSPGetPC(smoother,&subpc);CHKERRQ(ierr); 175 ierr = PCSetType(subpc,PCJACOBI);CHKERRQ(ierr); 176 ierr = KSPChebychevSetEigenvalues(smoother, emax, emin);CHKERRQ(ierr); /* need auto !!!!*/ 177 } 178 ierr = PCSetFromOptions_MG(pc); CHKERRQ(ierr); /* should be called in PCSetFromOptions_GAMG(), but cannot be called prior to PCMGSetLevels() */ 179 { 180 PetscBool galerkin; 181 ierr = PCMGGetGalerkin( pc, &galerkin); CHKERRQ(ierr); 182 if(galerkin){ 183 SETERRQ(wcomm,PETSC_ERR_ARG_WRONG, "GAMG does galerkin manually so it must not be used in PC_MG."); 184 } 185 } 186 /* create coarse level and the interpolation between the levels */ 187 for (level=0,lidx=pc_gamg->m_Nlevels-1; level<fine_level; level++,lidx--){ 188 level1 = level + 1; 189 /* PetscInt MM,NN; */ 190 /* ierr = MatGetSize( Rarr[lidx], &MM, &NN );CHKERRQ(ierr); */ 191 /* PetscPrintf(PETSC_COMM_WORLD,"%s Set P(%d,%d) on level %d (%d)\n",__FUNCT__,MM,NN,level1,lidx); */ 192 ierr = PCMGSetInterpolation(pc,level1,Rarr[lidx]);CHKERRQ(ierr); 193 if(!PETSC_TRUE) { 194 PetscViewer viewer; 195 ierr = PetscViewerASCIIOpen(PETSC_COMM_SELF, "Rmat.m", &viewer); CHKERRQ(ierr); 196 ierr = PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB); CHKERRQ(ierr); 197 ierr = MatView(Rarr[level1],viewer);CHKERRQ(ierr); 198 ierr = PetscViewerDestroy( &viewer ); 199 } 200 KSP smoother; 201 ierr = PCMGGetSmoother(pc,level,&smoother); CHKERRQ(ierr); 202 ierr = KSPSetOperators( smoother, Aarr[lidx], Aarr[lidx], DIFFERENT_NONZERO_PATTERN ); 203 CHKERRQ(ierr); 204 /* ierr = MatGetSize( Aarr[lidx], &MM, &NN ); CHKERRQ(ierr); */ 205 /* PetscPrintf(PETSC_COMM_WORLD,"%s Set A(%d,%d) on level %d (%d)\n",__FUNCT__,MM,NN,level,lidx); */ 206 } 207 { /* fine level (no P) */ 208 KSP smoother; 209 ierr = PCMGGetSmoother(pc,fine_level,&smoother); CHKERRQ(ierr); 210 ierr = KSPSetOperators( smoother, Aarr[0], Aarr[0], DIFFERENT_NONZERO_PATTERN ); 211 CHKERRQ(ierr); 212 /* PetscInt MM,NN; */ 213 /* ierr = MatGetSize( Aarr[0], &MM, &NN );CHKERRQ(ierr); */ 214 /* PetscPrintf(PETSC_COMM_WORLD,"%s Set A(%d,%d) on level %d (%d)\n",__FUNCT__,MM,NN,fine_level,0); */ 215 } 216 217 /* setupcalled is set to 0 so that MG is setup from scratch */ 218 pc->setupcalled = 0; 219 ierr = PCSetUp_MG(pc);CHKERRQ(ierr); 220 PetscFunctionReturn(0); 221 } 222 223 /* -------------------------------------------------------------------------- */ 224 /* 225 PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner 226 that was created with PCCreate_GAMG(). 227 228 Input Parameter: 229 . pc - the preconditioner context 230 231 Application Interface Routine: PCDestroy() 232 */ 233 #undef __FUNCT__ 234 #define __FUNCT__ "PCDestroy_GAMG" 235 PetscErrorCode PCDestroy_GAMG(PC pc) 236 { 237 PetscErrorCode ierr; 238 PC_MG *mg = (PC_MG*)pc->data; 239 PC_GAMG *pc_gamg= (PC_GAMG*)mg->innerctx; 240 241 PetscFunctionBegin; 242 ierr = PCReset_GAMG(pc);CHKERRQ(ierr); 243 ierr = PetscFree(pc_gamg);CHKERRQ(ierr); 244 ierr = PCDestroy_MG(pc);CHKERRQ(ierr); 245 PetscFunctionReturn(0); 246 } 247 248 #undef __FUNCT__ 249 #define __FUNCT__ "PCSetFromOptions_GAMG" 250 PetscErrorCode PCSetFromOptions_GAMG(PC pc) 251 { 252 /* PetscErrorCode ierr; */ 253 /* PC_MG *mg = (PC_MG*)pc->data; */ 254 /* PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; */ 255 /* MPI_Comm comm = ((PetscObject)pc)->comm; */ 256 257 PetscFunctionBegin; 258 PetscFunctionReturn(0); 259 } 260 261 /* -------------------------------------------------------------------------- */ 262 /* 263 PCCreate_GAMG - Creates a GAMG preconditioner context, PC_GAMG 264 265 Input Parameter: 266 . pc - the preconditioner context 267 268 Application Interface Routine: PCCreate() 269 270 */ 271 /* MC 272 PCGAMG - Use algebraic multigrid preconditioning. This preconditioner requires you provide 273 fine grid discretization matrix and coordinates on the fine grid. 274 275 Options Database Key: 276 Multigrid options(inherited) 277 + -pc_mg_cycles <1>: 1 for V cycle, 2 for W-cycle (MGSetCycles) 278 . -pc_mg_smoothup <1>: Number of post-smoothing steps (MGSetNumberSmoothUp) 279 . -pc_mg_smoothdown <1>: Number of pre-smoothing steps (MGSetNumberSmoothDown) 280 -pc_mg_type <multiplicative>: (one of) additive multiplicative full cascade kascade 281 GAMG options: 282 283 Level: intermediate 284 Concepts: multigrid 285 286 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, 287 PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), MPSetCycles(), PCMGSetNumberSmoothDown(), 288 PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(), 289 PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(), 290 PCMGSetCyclesOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR() 291 M */ 292 293 EXTERN_C_BEGIN 294 #undef __FUNCT__ 295 #define __FUNCT__ "PCCreate_GAMG" 296 PetscErrorCode PCCreate_GAMG(PC pc) 297 { 298 PetscErrorCode ierr; 299 PC_GAMG *pc_gamg; 300 PC_MG *mg; 301 302 PetscFunctionBegin; 303 /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */ 304 ierr = PCSetType(pc,PCMG);CHKERRQ(ierr); /* calls PCCreate_MG() and MGCreate_Private() */ 305 ierr = PetscObjectChangeTypeName((PetscObject)pc,PCGAMG);CHKERRQ(ierr); 306 307 /* create a supporting struct and attach it to pc */ 308 ierr = PetscNewLog(pc,PC_GAMG,&pc_gamg);CHKERRQ(ierr); 309 mg = (PC_MG*)pc->data; 310 mg->innerctx = pc_gamg; 311 312 pc_gamg->m_Nlevels = -1; 313 314 /* overwrite the pointers of PCMG by the functions of PCGAMG */ 315 pc->ops->setfromoptions = PCSetFromOptions_GAMG; 316 pc->ops->setup = PCSetUp_GAMG; 317 pc->ops->reset = PCReset_GAMG; 318 pc->ops->destroy = PCDestroy_GAMG; 319 320 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 321 "PCSetCoordinates_C", 322 "PCSetCoordinates_GAMG", 323 PCSetCoordinates_GAMG);CHKERRQ(ierr); 324 PetscFunctionReturn(0); 325 } 326 EXTERN_C_END 327