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 #include <assert.h> 9 10 /* Private context for the GAMG preconditioner */ 11 typedef struct gamg_TAG { 12 PetscInt m_dim; 13 PetscInt m_Nlevels; 14 PetscInt m_data_sz; 15 PetscReal *m_data; /* blocked vector of vertex data on fine grid (coordinates) */ 16 } PC_GAMG; 17 18 /* -----------------------------------------------------------------------------*/ 19 #undef __FUNCT__ 20 #define __FUNCT__ "PCReset_GAMG" 21 PetscErrorCode PCReset_GAMG(PC pc) 22 { 23 PetscErrorCode ierr; 24 PC_MG *mg = (PC_MG*)pc->data; 25 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 26 27 PetscFunctionBegin; 28 ierr = PetscFree(pc_gamg->m_data);CHKERRQ(ierr); 29 PetscFunctionReturn(0); 30 } 31 32 /* -------------------------------------------------------------------------- */ 33 /* 34 PCSetCoordinates_GAMG 35 36 Input Parameter: 37 . pc - the preconditioner context 38 */ 39 EXTERN_C_BEGIN 40 #undef __FUNCT__ 41 #define __FUNCT__ "PCSetCoordinates_GAMG" 42 PetscErrorCode PCSetCoordinates_GAMG(PC pc, const int ndm,PetscReal *coords ) 43 { 44 PC_MG *mg = (PC_MG*)pc->data; 45 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 46 PetscErrorCode ierr; 47 PetscInt bs, my0, tt; 48 Mat mat = pc->pmat; 49 PetscInt arrsz; 50 51 PetscFunctionBegin; 52 ierr = MatGetBlockSize( mat, &bs ); CHKERRQ( ierr ); 53 ierr = MatGetOwnershipRange( mat, &my0, &tt ); CHKERRQ(ierr); 54 arrsz = (tt-my0)/bs*ndm; 55 56 // put coordinates 57 if (!pc_gamg->m_data || (pc_gamg->m_data_sz != arrsz)) { 58 ierr = PetscFree( pc_gamg->m_data ); CHKERRQ(ierr); 59 ierr = PetscMalloc(arrsz*sizeof(double), &pc_gamg->m_data ); CHKERRQ(ierr); 60 } 61 62 /* copy data in */ 63 for(tt=0;tt<arrsz;tt++){ 64 pc_gamg->m_data[tt] = coords[tt]; 65 } 66 pc_gamg->m_data_sz = arrsz; 67 pc_gamg->m_dim = ndm; 68 69 PetscFunctionReturn(0); 70 } 71 EXTERN_C_END 72 73 /* -------------------------------------------------------------------------- */ 74 /* 75 partitionLevel 76 77 Input Parameter: 78 . a_Amat_fine - matrix on this fine (k) level 79 . a_dime - 2 or 3 80 In/Output Parameter: 81 . a_P_inout - prolongation operator to the next level (k-1) 82 . a_coarse_crds - coordinates that need to be moved 83 . a_active_proc - number of active procs 84 Output Parameter: 85 . a_Amat_crs - coarse matrix that is created (k-1) 86 */ 87 88 PetscMPIInt n_active_procs; // hack!!!!! 89 90 #undef __FUNCT__ 91 #define __FUNCT__ "partitionLevel" 92 PetscErrorCode partitionLevel( Mat a_Amat_fine, 93 PetscInt a_dim, 94 Mat *a_P_inout, 95 PetscReal **a_coarse_crds, 96 PetscMPIInt *a_active_proc, 97 Mat *a_Amat_crs 98 ) 99 { 100 PetscErrorCode ierr; 101 Mat Amat, Pnew, Pold = *a_P_inout; 102 IS new_indices,isnum; 103 MPI_Comm wcomm = ((PetscObject)a_Amat_fine)->comm; 104 PetscMPIInt nactive,mype,npe; 105 PetscInt Istart,Iend,Istart0,Iend0,ncrs0,ncrs_new,bs=1; /* bs ??? */ 106 107 PetscFunctionBegin; 108 ierr = MPI_Comm_rank(wcomm,&mype);CHKERRQ(ierr); 109 ierr = MPI_Comm_size(wcomm,&npe);CHKERRQ(ierr); 110 /* RAP */ 111 ierr = MatPtAP( a_Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, &Amat ); CHKERRQ(ierr); 112 ierr = MatGetOwnershipRange( Amat, &Istart0, &Iend0 ); CHKERRQ(ierr); /* x2 size of 'a_coarse_crds' */ 113 ncrs0 = (Iend0 - Istart0)/bs; 114 115 /* Repartition Amat_{k} and move colums of P^{k}_{k-1} and coordinates accordingly */ 116 { 117 PetscInt neq,N,counts[npe]; 118 IS isnewproc; 119 PetscMPIInt factstart,fact,new_npe,targ_npe; 120 121 ierr = MatGetSize( Amat, &neq, &N );CHKERRQ(ierr); 122 #define MIN_EQ_PROC 100 123 nactive = *a_active_proc; 124 targ_npe = neq/MIN_EQ_PROC; /* hardwire min. number of eq/proc */ 125 if( targ_npe == 0 || neq < 1000 ) new_npe = 1; /* chop coarsest grid */ 126 else if (targ_npe > nactive ) new_npe = nactive; /* no change */ 127 else { 128 new_npe = -9999; 129 factstart = nactive; 130 for(fact=factstart;fact>0;fact--){ /* try to find a better number of procs */ 131 if( nactive%fact==0 && neq/(nactive/fact) > MIN_EQ_PROC ) { 132 new_npe = nactive/fact; 133 } 134 } 135 assert(new_npe != -9999); 136 } 137 *a_active_proc = new_npe; /* output for next time */ 138 139 n_active_procs = new_npe; // hack!!!!! 140 PetscPrintf(PETSC_COMM_WORLD,"\t\t%s n_active_procs=%d\n",__FUNCT__,n_active_procs); 141 { /* partition: get 'isnewproc' */ 142 MatPartitioning mpart; 143 PetscMPIInt factor = npe/new_npe; 144 /* create sub communicator */ 145 MPI_Comm cm,new_comm; 146 int membershipKey = mype % factor; 147 148 ierr = MPI_Comm_split(wcomm, membershipKey, mype, &cm); CHKERRQ(ierr); 149 ierr = PetscCommDuplicate( cm, &new_comm, PETSC_NULL ); CHKERRQ(ierr); 150 ierr = MPI_Comm_free( &cm ); CHKERRQ(ierr); 151 152 /* use 'new_comm' -- todo!!! */ 153 ierr = MatPartitioningCreate( wcomm, &mpart ); CHKERRQ(ierr); 154 ierr = MatPartitioningSetAdjacency( mpart, Amat ); CHKERRQ(ierr); 155 ierr = MatPartitioningSetFromOptions( mpart ); CHKERRQ(ierr); 156 ierr = MatPartitioningApply( mpart, &isnewproc ); CHKERRQ(ierr); 157 ierr = MatPartitioningDestroy( &mpart ); CHKERRQ(ierr); 158 159 ierr = MPI_Comm_free( &new_comm ); CHKERRQ(ierr); 160 161 /* spread partitioning across machine - probably the right thing to do but machine specific */ 162 /* how to do this??? scale IS by 'factor' */ 163 } 164 165 /* 166 Create an index set from the isnewproc index set to indicate the mapping TO 167 */ 168 ierr = ISPartitioningToNumbering( isnewproc, &isnum ); CHKERRQ(ierr); 169 /* 170 Determine how many elements are assigned to each processor 171 */ 172 ierr = ISPartitioningCount( isnewproc, npe, counts ); CHKERRQ(ierr); 173 ierr = ISDestroy( &isnewproc ); CHKERRQ(ierr); 174 ncrs_new = counts[mype]; 175 } 176 { /* Create a vector to contain the newly ordered element information */ 177 const PetscInt *idx; 178 PetscInt i,j,k; 179 IS isscat; 180 PetscScalar *array; 181 Vec src_crd, dest_crd; 182 PetscReal *coords = *a_coarse_crds; 183 VecScatter vecscat; 184 185 ierr = VecCreate( wcomm, &dest_crd ); 186 ierr = VecSetSizes( dest_crd, a_dim*ncrs_new, PETSC_DECIDE ); CHKERRQ(ierr); 187 ierr = VecSetFromOptions( dest_crd ); CHKERRQ(ierr); /*funny vector-get global options?*/ 188 /* 189 There are three data items per cell (element), the integer vertex numbers of its three 190 coordinates (we convert to double to use the scatter) (one can think 191 of the vectors of having a block size of 3, then there is one index in idx[] for each element) 192 */ 193 { 194 PetscInt tidx[ncrs0*a_dim]; 195 ierr = ISGetIndices( isnum, &idx ); CHKERRQ(ierr); 196 for (i=0,j=0; i<ncrs0; i++) { 197 for (k=0; k<a_dim; k++) tidx[j++] = idx[i]*a_dim + k; 198 } 199 ierr = ISCreateGeneral( wcomm, a_dim*ncrs0, tidx, PETSC_COPY_VALUES, &isscat ); 200 CHKERRQ(ierr); 201 ierr = ISRestoreIndices( isnum, &idx ); CHKERRQ(ierr); 202 } 203 /* 204 Create a vector to contain the original vertex information for each element 205 */ 206 ierr = VecCreateSeq( PETSC_COMM_SELF, a_dim*ncrs0, &src_crd ); CHKERRQ(ierr); 207 ierr = VecGetArray( src_crd, &array ); CHKERRQ(ierr); 208 for (i=0; i<a_dim*ncrs0; i++) array[i] = coords[i]; 209 ierr = VecRestoreArray( src_crd, &array ); CHKERRQ(ierr); 210 /* 211 Scatter the element vertex information (still in the original vertex ordering) 212 to the correct processor 213 */ 214 ierr = VecScatterCreate( src_crd, PETSC_NULL, dest_crd, isscat, &vecscat); 215 CHKERRQ(ierr); 216 ierr = ISDestroy( &isscat ); CHKERRQ(ierr); 217 ierr = VecScatterBegin(vecscat,src_crd,dest_crd,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 218 ierr = VecScatterEnd(vecscat,src_crd,dest_crd,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 219 ierr = VecScatterDestroy( &vecscat ); CHKERRQ(ierr); 220 ierr = VecDestroy( &src_crd ); CHKERRQ(ierr); 221 /* 222 Put the element vertex data into a new allocation of the gdata->ele 223 */ 224 ierr = PetscFree( *a_coarse_crds ); CHKERRQ(ierr); 225 ierr = PetscMalloc( a_dim*ncrs_new*sizeof(PetscReal), a_coarse_crds ); CHKERRQ(ierr); 226 coords = *a_coarse_crds; /* convient */ 227 ierr = VecGetArray( dest_crd, &array ); CHKERRQ(ierr); 228 for (i=0; i<a_dim*ncrs_new; i++) coords[i] = PetscRealPart(array[i]); 229 ierr = VecRestoreArray( dest_crd, &array ); CHKERRQ(ierr); 230 ierr = VecDestroy( &dest_crd ); CHKERRQ(ierr); 231 } 232 /* 233 Invert for MatGetSubMatrix 234 */ 235 ierr = ISInvertPermutation( isnum, ncrs_new, &new_indices ); CHKERRQ(ierr); 236 ierr = ISSort( new_indices ); CHKERRQ(ierr); /* is this needed? */ 237 ierr = ISDestroy( &isnum ); CHKERRQ(ierr); 238 /* A_crs output */ 239 ierr = MatGetSubMatrix( Amat, new_indices, new_indices, MAT_INITIAL_MATRIX, a_Amat_crs ); 240 CHKERRQ(ierr); 241 ierr = MatDestroy( &Amat ); CHKERRQ(ierr); 242 Amat = *a_Amat_crs; 243 /* prolongator */ 244 ierr = MatGetOwnershipRange( Pold, &Istart, &Iend ); CHKERRQ(ierr); 245 { 246 IS findices; 247 ierr = ISCreateStride(wcomm,Iend-Istart,Istart,1,&findices); CHKERRQ(ierr); 248 ierr = MatGetSubMatrix( Pold, findices, new_indices, MAT_INITIAL_MATRIX, &Pnew ); 249 CHKERRQ(ierr); 250 ierr = ISDestroy( &findices ); CHKERRQ(ierr); 251 } 252 ierr = MatDestroy( a_P_inout ); CHKERRQ(ierr); 253 *a_P_inout = Pnew; /* output */ 254 255 ierr = ISDestroy( &new_indices ); CHKERRQ(ierr); 256 257 PetscFunctionReturn(0); 258 } 259 260 /* -------------------------------------------------------------------------- */ 261 /* 262 PCSetUp_GAMG - Prepares for the use of the GAMG preconditioner 263 by setting data structures and options. 264 265 Input Parameter: 266 . pc - the preconditioner context 267 268 Application Interface Routine: PCSetUp() 269 270 Notes: 271 The interface routine PCSetUp() is not usually called directly by 272 the user, but instead is called by PCApply() if necessary. 273 */ 274 extern PetscErrorCode PCSetFromOptions_MG(PC); 275 extern PetscErrorCode PCReset_MG(PC); 276 extern PetscErrorCode createProlongation( Mat, PetscReal [], const PetscInt, 277 Mat *, PetscReal **, PetscBool *a_isOK ); 278 #undef __FUNCT__ 279 #define __FUNCT__ "PCSetUp_GAMG" 280 PetscErrorCode PCSetUp_GAMG( PC pc ) 281 { 282 PetscErrorCode ierr; 283 PC_MG *mg = (PC_MG*)pc->data; 284 PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; 285 Mat Amat = pc->mat, Pmat = pc->pmat; 286 PetscBool isSeq, isMPI; 287 PetscInt fine_level, level, level1, M, N, bs, lidx; 288 MPI_Comm wcomm = ((PetscObject)pc)->comm; 289 PetscMPIInt mype,npe,nactivepe; 290 #define GAMG_MAXLEVELS 20 291 Mat Aarr[GAMG_MAXLEVELS], Rarr[GAMG_MAXLEVELS]; PetscReal *coarse_crds = 0, *crds = pc_gamg->m_data; 292 PetscBool isOK; 293 294 PetscFunctionBegin; 295 ierr = MPI_Comm_rank(wcomm,&mype);CHKERRQ(ierr); 296 ierr = MPI_Comm_size(wcomm,&npe);CHKERRQ(ierr); 297 if (pc->setupcalled){ 298 /* no state data in GAMG to destroy (now) */ 299 ierr = PCReset_MG(pc); CHKERRQ(ierr); 300 } 301 if (!pc_gamg->m_data) SETERRQ(wcomm,PETSC_ERR_SUP,"PCSetUp_GAMG called before PCSetCoordinates"); 302 /* setup special features of PCGAMG */ 303 ierr = PetscTypeCompare((PetscObject)Amat, MATSEQAIJ, &isSeq);CHKERRQ(ierr); 304 ierr = PetscTypeCompare((PetscObject)Amat, MATMPIAIJ, &isMPI);CHKERRQ(ierr); 305 if (isMPI) { 306 } else if (isSeq) { 307 } 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); 308 309 /* GAMG requires input of fine-grid matrix. It determines nlevels. */ 310 ierr = MatGetBlockSize( Amat, &bs ); CHKERRQ(ierr); 311 if(bs!=1) SETERRQ1(wcomm,PETSC_ERR_ARG_WRONG, "GAMG only supports scalar prblems bs = '%d'.",bs); 312 313 /* Get A_i and R_i */ 314 ierr = MatGetSize( Amat, &M, &N );CHKERRQ(ierr); 315 PetscPrintf(PETSC_COMM_WORLD,"[%d]%s level %d N=%d\n",0,__FUNCT__,0,N); 316 for (level=0, Aarr[0] = Pmat, nactivepe = npe; 317 level < GAMG_MAXLEVELS-1 && (level==0 || M>1000); /* hard wired stopping logic */ 318 level++ ) { 319 level1 = level + 1; 320 ierr = createProlongation( Aarr[level], crds, pc_gamg->m_dim, 321 &Rarr[level1], &coarse_crds, &isOK ); 322 CHKERRQ(ierr); 323 ierr = PetscFree( crds ); CHKERRQ( ierr ); 324 if(level==0) Aarr[0] = Amat; /* use Pmat for finest level setup, but use mat for solver */ 325 if( isOK ) { 326 ierr = partitionLevel( Aarr[level], pc_gamg->m_dim, &Rarr[level1], &coarse_crds, &nactivepe, &Aarr[level1] ); 327 CHKERRQ(ierr); 328 ierr = MatGetSize( Aarr[level1], &M, &N );CHKERRQ(ierr); 329 PetscPrintf(PETSC_COMM_WORLD,"[%d]%s done level %d N=%d, active pe = %d\n",0,__FUNCT__,level1,N,nactivepe); 330 } 331 else{ 332 break; 333 } 334 crds = coarse_crds; 335 } 336 ierr = PetscFree( coarse_crds ); CHKERRQ( ierr ); 337 PetscPrintf(PETSC_COMM_WORLD,"\t[%d]%s %d levels\n",0,__FUNCT__,level + 1); 338 pc_gamg->m_data = 0; /* destroyed coordinate data */ 339 pc_gamg->m_Nlevels = level + 1; 340 fine_level = level; 341 ierr = PCMGSetLevels(pc,pc_gamg->m_Nlevels,PETSC_NULL);CHKERRQ(ierr); 342 343 /* set default smoothers */ 344 PetscReal emax = 2.0, emin; 345 for (level=1; level<=fine_level; level++){ 346 KSP smoother; PC subpc; 347 ierr = PCMGGetSmoother(pc,level,&smoother);CHKERRQ(ierr); 348 ierr = KSPSetType(smoother,KSPCHEBYCHEV);CHKERRQ(ierr); 349 emin = emax/10.0; /* fix!!! */ 350 ierr = KSPGetPC(smoother,&subpc);CHKERRQ(ierr); 351 ierr = PCSetType(subpc,PCJACOBI);CHKERRQ(ierr); 352 ierr = KSPChebychevSetEigenvalues(smoother, emax, emin);CHKERRQ(ierr); /* need auto !!!!*/ 353 } 354 ierr = PCSetFromOptions_MG(pc); CHKERRQ(ierr); /* should be called in PCSetFromOptions_GAMG(), but cannot be called prior to PCMGSetLevels() */ 355 { 356 PetscBool galerkin; 357 ierr = PCMGGetGalerkin( pc, &galerkin); CHKERRQ(ierr); 358 if(galerkin){ 359 SETERRQ(wcomm,PETSC_ERR_ARG_WRONG, "GAMG does galerkin manually so it must not be used in PC_MG."); 360 } 361 } 362 /* create coarse level and the interpolation between the levels */ 363 for (level=0,lidx=pc_gamg->m_Nlevels-1; level<fine_level; level++,lidx--){ 364 level1 = level + 1; 365 ierr = PCMGSetInterpolation(pc,level1,Rarr[lidx]);CHKERRQ(ierr); 366 if(!PETSC_TRUE) { 367 PetscViewer viewer; char fname[32]; 368 sprintf(fname,"Amat_%d.m",lidx); 369 ierr = PetscViewerASCIIOpen( wcomm, fname, &viewer ); CHKERRQ(ierr); 370 ierr = PetscViewerSetFormat( viewer, PETSC_VIEWER_ASCII_MATLAB); CHKERRQ(ierr); 371 ierr = MatView( Aarr[lidx], viewer ); CHKERRQ(ierr); 372 ierr = PetscViewerDestroy( &viewer ); 373 } 374 ierr = MatDestroy( &Rarr[lidx] ); CHKERRQ(ierr); 375 { 376 KSP smoother; 377 ierr = PCMGGetSmoother(pc,level,&smoother); CHKERRQ(ierr); 378 ierr = KSPSetOperators( smoother, Aarr[lidx], Aarr[lidx], DIFFERENT_NONZERO_PATTERN ); 379 CHKERRQ(ierr); 380 ierr = MatDestroy( &Aarr[lidx] ); CHKERRQ(ierr); 381 } 382 } 383 { /* fine level (no P) */ 384 KSP smoother; 385 ierr = PCMGGetSmoother(pc,fine_level,&smoother); CHKERRQ(ierr); 386 ierr = KSPSetOperators( smoother, Aarr[0], Aarr[0], DIFFERENT_NONZERO_PATTERN ); 387 CHKERRQ(ierr); 388 } 389 390 /* setupcalled is set to 0 so that MG is setup from scratch */ 391 pc->setupcalled = 0; 392 ierr = PCSetUp_MG(pc);CHKERRQ(ierr); 393 PetscFunctionReturn(0); 394 } 395 396 /* -------------------------------------------------------------------------- */ 397 /* 398 PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner 399 that was created with PCCreate_GAMG(). 400 401 Input Parameter: 402 . pc - the preconditioner context 403 404 Application Interface Routine: PCDestroy() 405 */ 406 #undef __FUNCT__ 407 #define __FUNCT__ "PCDestroy_GAMG" 408 PetscErrorCode PCDestroy_GAMG(PC pc) 409 { 410 PetscErrorCode ierr; 411 PC_MG *mg = (PC_MG*)pc->data; 412 PC_GAMG *pc_gamg= (PC_GAMG*)mg->innerctx; 413 414 PetscFunctionBegin; 415 ierr = PCReset_GAMG(pc);CHKERRQ(ierr); 416 ierr = PetscFree(pc_gamg);CHKERRQ(ierr); 417 ierr = PCDestroy_MG(pc);CHKERRQ(ierr); 418 PetscFunctionReturn(0); 419 } 420 421 #undef __FUNCT__ 422 #define __FUNCT__ "PCSetFromOptions_GAMG" 423 PetscErrorCode PCSetFromOptions_GAMG(PC pc) 424 { 425 /* PetscErrorCode ierr; */ 426 /* PC_MG *mg = (PC_MG*)pc->data; */ 427 /* PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx; */ 428 /* MPI_Comm comm = ((PetscObject)pc)->comm; */ 429 430 PetscFunctionBegin; 431 PetscFunctionReturn(0); 432 } 433 434 /* -------------------------------------------------------------------------- */ 435 /* 436 PCCreate_GAMG - Creates a GAMG preconditioner context, PC_GAMG 437 438 Input Parameter: 439 . pc - the preconditioner context 440 441 Application Interface Routine: PCCreate() 442 443 */ 444 /* MC 445 PCGAMG - Use algebraic multigrid preconditioning. This preconditioner requires you provide 446 fine grid discretization matrix and coordinates on the fine grid. 447 448 Options Database Key: 449 Multigrid options(inherited) 450 + -pc_mg_cycles <1>: 1 for V cycle, 2 for W-cycle (MGSetCycles) 451 . -pc_mg_smoothup <1>: Number of post-smoothing steps (MGSetNumberSmoothUp) 452 . -pc_mg_smoothdown <1>: Number of pre-smoothing steps (MGSetNumberSmoothDown) 453 -pc_mg_type <multiplicative>: (one of) additive multiplicative full cascade kascade 454 GAMG options: 455 456 Level: intermediate 457 Concepts: multigrid 458 459 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC, PCMGType, 460 PCMGSetLevels(), PCMGGetLevels(), PCMGSetType(), MPSetCycles(), PCMGSetNumberSmoothDown(), 461 PCMGSetNumberSmoothUp(), PCMGGetCoarseSolve(), PCMGSetResidual(), PCMGSetInterpolation(), 462 PCMGSetRestriction(), PCMGGetSmoother(), PCMGGetSmootherUp(), PCMGGetSmootherDown(), 463 PCMGSetCyclesOnLevel(), PCMGSetRhs(), PCMGSetX(), PCMGSetR() 464 M */ 465 466 EXTERN_C_BEGIN 467 #undef __FUNCT__ 468 #define __FUNCT__ "PCCreate_GAMG" 469 PetscErrorCode PCCreate_GAMG(PC pc) 470 { 471 PetscErrorCode ierr; 472 PC_GAMG *pc_gamg; 473 PC_MG *mg; 474 475 PetscFunctionBegin; 476 /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */ 477 ierr = PCSetType(pc,PCMG);CHKERRQ(ierr); /* calls PCCreate_MG() and MGCreate_Private() */ 478 ierr = PetscObjectChangeTypeName((PetscObject)pc,PCGAMG);CHKERRQ(ierr); 479 480 /* create a supporting struct and attach it to pc */ 481 ierr = PetscNewLog(pc,PC_GAMG,&pc_gamg);CHKERRQ(ierr); 482 mg = (PC_MG*)pc->data; 483 mg->innerctx = pc_gamg; 484 485 pc_gamg->m_Nlevels = -1; 486 487 /* overwrite the pointers of PCMG by the functions of PCGAMG */ 488 pc->ops->setfromoptions = PCSetFromOptions_GAMG; 489 pc->ops->setup = PCSetUp_GAMG; 490 pc->ops->reset = PCReset_GAMG; 491 pc->ops->destroy = PCDestroy_GAMG; 492 493 ierr = PetscObjectComposeFunctionDynamic( (PetscObject)pc, 494 "PCSetCoordinates_C", 495 "PCSetCoordinates_GAMG", 496 PCSetCoordinates_GAMG);CHKERRQ(ierr); 497 PetscFunctionReturn(0); 498 } 499 EXTERN_C_END 500