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