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