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