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