xref: /petsc/src/ksp/pc/impls/gamg/gamg.c (revision cab9ed1e6be87951742340267410beca15d62530)
1 /*
2  GAMG geometric-algebric multigrid PC - Mark Adams 2011
3  */
4 #include <petsc/private/matimpl.h>
5 #include <../src/ksp/pc/impls/gamg/gamg.h>           /*I "petscpc.h" I*/
6 #include <../src/ksp/pc/impls/bjacobi/bjacobi.h> /* Hack to access same_local_solves */
7 
8 #if defined PETSC_GAMG_USE_LOG
9 PetscLogEvent petsc_gamg_setup_events[NUM_SET];
10 #endif
11 
12 #if defined PETSC_USE_LOG
13 PetscLogEvent PC_GAMGGraph_AGG;
14 PetscLogEvent PC_GAMGGraph_GEO;
15 PetscLogEvent PC_GAMGCoarsen_AGG;
16 PetscLogEvent PC_GAMGCoarsen_GEO;
17 PetscLogEvent PC_GAMGProlongator_AGG;
18 PetscLogEvent PC_GAMGProlongator_GEO;
19 PetscLogEvent PC_GAMGOptProlongator_AGG;
20 #endif
21 
22 #define GAMG_MAXLEVELS 30
23 
24 /* #define GAMG_STAGES */
25 #if (defined PETSC_GAMG_USE_LOG && defined GAMG_STAGES)
26 static PetscLogStage gamg_stages[GAMG_MAXLEVELS];
27 #endif
28 
29 static PetscFunctionList GAMGList = 0;
30 static PetscBool PCGAMGPackageInitialized;
31 
32 /* ----------------------------------------------------------------------------- */
33 #undef __FUNCT__
34 #define __FUNCT__ "PCReset_GAMG"
35 PetscErrorCode PCReset_GAMG(PC pc)
36 {
37   PetscErrorCode ierr;
38   PC_MG          *mg      = (PC_MG*)pc->data;
39   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
40 
41   PetscFunctionBegin;
42   if (pc_gamg->data) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_PLIB,"This should not happen, cleaned up in SetUp\n");
43   pc_gamg->data_sz = 0;
44   ierr = PetscFree(pc_gamg->orig_data);CHKERRQ(ierr);
45   PetscFunctionReturn(0);
46 }
47 
48 /* -------------------------------------------------------------------------- */
49 /*
50    PCGAMGCreateLevel_GAMG: create coarse op with RAP.  repartition and/or reduce number
51      of active processors.
52 
53    Input Parameter:
54    . pc - parameters + side effect: coarse data in 'pc_gamg->data' and
55           'pc_gamg->data_sz' are changed via repartitioning/reduction.
56    . Amat_fine - matrix on this fine (k) level
57    . cr_bs - coarse block size
58    In/Output Parameter:
59    . a_P_inout - prolongation operator to the next level (k-->k-1)
60    . a_nactive_proc - number of active procs
61    Output Parameter:
62    . a_Amat_crs - coarse matrix that is created (k-1)
63 */
64 
65 #undef __FUNCT__
66 #define __FUNCT__ "PCGAMGCreateLevel_GAMG"
67 static PetscErrorCode PCGAMGCreateLevel_GAMG(PC pc,Mat Amat_fine,PetscInt cr_bs,Mat *a_P_inout,Mat *a_Amat_crs,PetscMPIInt *a_nactive_proc,IS * Pcolumnperm)
68 {
69   PetscErrorCode  ierr;
70   PC_MG           *mg         = (PC_MG*)pc->data;
71   PC_GAMG         *pc_gamg    = (PC_GAMG*)mg->innerctx;
72   Mat             Cmat,Pold=*a_P_inout;
73   MPI_Comm        comm;
74   PetscMPIInt     rank,size,new_size,nactive=*a_nactive_proc;
75   PetscInt        ncrs_eq,ncrs,f_bs;
76 
77   PetscFunctionBegin;
78   ierr = PetscObjectGetComm((PetscObject)Amat_fine,&comm);CHKERRQ(ierr);
79   ierr = MPI_Comm_rank(comm, &rank);CHKERRQ(ierr);
80   ierr = MPI_Comm_size(comm, &size);CHKERRQ(ierr);
81   ierr = MatGetBlockSize(Amat_fine, &f_bs);CHKERRQ(ierr);
82   ierr = MatPtAP(Amat_fine, Pold, MAT_INITIAL_MATRIX, 2.0, &Cmat);CHKERRQ(ierr);
83 
84   /* set 'ncrs' (nodes), 'ncrs_eq' (equations)*/
85   ierr = MatGetLocalSize(Cmat, &ncrs_eq, NULL);CHKERRQ(ierr);
86   if (pc_gamg->data_cell_rows>0) {
87     ncrs = pc_gamg->data_sz/pc_gamg->data_cell_cols/pc_gamg->data_cell_rows;
88   } else {
89     PetscInt  bs;
90     ierr = MatGetBlockSize(Cmat, &bs);CHKERRQ(ierr);
91     ncrs = ncrs_eq/bs;
92   }
93 
94   /* get number of PEs to make active 'new_size', reduce, can be any integer 1-P */
95   {
96     PetscInt ncrs_eq_glob;
97     ierr     = MatGetSize(Cmat, &ncrs_eq_glob, NULL);CHKERRQ(ierr);
98     new_size = (PetscMPIInt)((float)ncrs_eq_glob/(float)pc_gamg->min_eq_proc + 0.5); /* hardwire min. number of eq/proc */
99     if (!new_size) new_size = 1; /* not likely, posible? */
100     else if (new_size >= nactive) new_size = nactive; /* no change, rare */
101   }
102 
103   if (Pcolumnperm) *Pcolumnperm = NULL;
104 
105   if (!pc_gamg->repart && new_size==nactive) *a_Amat_crs = Cmat; /* output - no repartitioning or reduction - could bail here */
106   else {
107     PetscInt       *counts,*newproc_idx,ii,jj,kk,strideNew,*tidx,ncrs_new,ncrs_eq_new,nloc_old;
108     IS             is_eq_newproc,is_eq_num,is_eq_num_prim,new_eq_indices;
109 
110     nloc_old = ncrs_eq/cr_bs;
111     if (ncrs_eq % cr_bs) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ncrs_eq %D not divisible by cr_bs %D",ncrs_eq,cr_bs);
112 #if defined PETSC_GAMG_USE_LOG
113     ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET12],0,0,0,0);CHKERRQ(ierr);
114 #endif
115     /* make 'is_eq_newproc' */
116     ierr = PetscMalloc1(size, &counts);CHKERRQ(ierr);
117     if (pc_gamg->repart) {
118       /* Repartition Cmat_{k} and move colums of P^{k}_{k-1} and coordinates of primal part accordingly */
119       Mat adj;
120 
121      ierr = PetscInfo3(pc,"Repartition: size (active): %D --> %D, neq = %D\n",*a_nactive_proc,new_size,ncrs_eq);CHKERRQ(ierr);
122 
123       /* get 'adj' */
124       if (cr_bs == 1) {
125         ierr = MatConvert(Cmat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj);CHKERRQ(ierr);
126       } else {
127         /* make a scalar matrix to partition (no Stokes here) */
128         Mat               tMat;
129         PetscInt          Istart_crs,Iend_crs,ncols,jj,Ii;
130         const PetscScalar *vals;
131         const PetscInt    *idx;
132         PetscInt          *d_nnz, *o_nnz, M, N;
133         static PetscInt   llev = 0;
134         MatType           mtype;
135 
136         ierr = PetscMalloc2(ncrs, &d_nnz,ncrs, &o_nnz);CHKERRQ(ierr);
137         ierr = MatGetOwnershipRange(Cmat, &Istart_crs, &Iend_crs);CHKERRQ(ierr);
138         ierr = MatGetSize(Cmat, &M, &N);CHKERRQ(ierr);
139         for (Ii = Istart_crs, jj = 0; Ii < Iend_crs; Ii += cr_bs, jj++) {
140           ierr      = MatGetRow(Cmat,Ii,&ncols,0,0);CHKERRQ(ierr);
141           d_nnz[jj] = ncols/cr_bs;
142           o_nnz[jj] = ncols/cr_bs;
143           ierr      = MatRestoreRow(Cmat,Ii,&ncols,0,0);CHKERRQ(ierr);
144           if (d_nnz[jj] > ncrs) d_nnz[jj] = ncrs;
145           if (o_nnz[jj] > (M/cr_bs-ncrs)) o_nnz[jj] = M/cr_bs-ncrs;
146         }
147 
148         ierr = MatGetType(Amat_fine,&mtype);CHKERRQ(ierr);
149         ierr = MatCreate(comm, &tMat);CHKERRQ(ierr);
150         ierr = MatSetSizes(tMat, ncrs, ncrs,PETSC_DETERMINE, PETSC_DETERMINE);CHKERRQ(ierr);
151         ierr = MatSetType(tMat,mtype);CHKERRQ(ierr);
152         ierr = MatSeqAIJSetPreallocation(tMat,0,d_nnz);CHKERRQ(ierr);
153         ierr = MatMPIAIJSetPreallocation(tMat,0,d_nnz,0,o_nnz);CHKERRQ(ierr);
154         ierr = PetscFree2(d_nnz,o_nnz);CHKERRQ(ierr);
155 
156         for (ii = Istart_crs; ii < Iend_crs; ii++) {
157           PetscInt dest_row = ii/cr_bs;
158           ierr = MatGetRow(Cmat,ii,&ncols,&idx,&vals);CHKERRQ(ierr);
159           for (jj = 0; jj < ncols; jj++) {
160             PetscInt    dest_col = idx[jj]/cr_bs;
161             PetscScalar v        = 1.0;
162             ierr = MatSetValues(tMat,1,&dest_row,1,&dest_col,&v,ADD_VALUES);CHKERRQ(ierr);
163           }
164           ierr = MatRestoreRow(Cmat,ii,&ncols,&idx,&vals);CHKERRQ(ierr);
165         }
166         ierr = MatAssemblyBegin(tMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
167         ierr = MatAssemblyEnd(tMat,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
168 
169         if (llev++ == -1) {
170           PetscViewer viewer; char fname[32];
171           ierr = PetscSNPrintf(fname,sizeof(fname),"part_mat_%D.mat",llev);CHKERRQ(ierr);
172           PetscViewerBinaryOpen(comm,fname,FILE_MODE_WRITE,&viewer);
173           ierr = MatView(tMat, viewer);CHKERRQ(ierr);
174           ierr = PetscViewerDestroy(&viewer);CHKERRQ(ierr);
175         }
176         ierr = MatConvert(tMat, MATMPIADJ, MAT_INITIAL_MATRIX, &adj);CHKERRQ(ierr);
177         ierr = MatDestroy(&tMat);CHKERRQ(ierr);
178       } /* create 'adj' */
179 
180       { /* partition: get newproc_idx */
181         char            prefix[256];
182         const char      *pcpre;
183         const PetscInt  *is_idx;
184         MatPartitioning mpart;
185         IS              proc_is;
186         PetscInt        targetPE;
187 
188         ierr = MatPartitioningCreate(comm, &mpart);CHKERRQ(ierr);
189         ierr = MatPartitioningSetAdjacency(mpart, adj);CHKERRQ(ierr);
190         ierr = PCGetOptionsPrefix(pc, &pcpre);CHKERRQ(ierr);
191         ierr = PetscSNPrintf(prefix,sizeof(prefix),"%spc_gamg_",pcpre ? pcpre : "");CHKERRQ(ierr);
192         ierr = PetscObjectSetOptionsPrefix((PetscObject)mpart,prefix);CHKERRQ(ierr);
193         ierr = MatPartitioningSetFromOptions(mpart);CHKERRQ(ierr);
194         ierr = MatPartitioningSetNParts(mpart, new_size);CHKERRQ(ierr);
195         ierr = MatPartitioningApply(mpart, &proc_is);CHKERRQ(ierr);
196         ierr = MatPartitioningDestroy(&mpart);CHKERRQ(ierr);
197 
198         /* collect IS info */
199         ierr     = PetscMalloc1(ncrs_eq, &newproc_idx);CHKERRQ(ierr);
200         ierr     = ISGetIndices(proc_is, &is_idx);CHKERRQ(ierr);
201         targetPE = 1; /* bring to "front" of machine */
202         /*targetPE = size/new_size;*/ /* spread partitioning across machine */
203         for (kk = jj = 0 ; kk < nloc_old ; kk++) {
204           for (ii = 0 ; ii < cr_bs ; ii++, jj++) {
205             newproc_idx[jj] = is_idx[kk] * targetPE; /* distribution */
206           }
207         }
208         ierr = ISRestoreIndices(proc_is, &is_idx);CHKERRQ(ierr);
209         ierr = ISDestroy(&proc_is);CHKERRQ(ierr);
210       }
211       ierr = MatDestroy(&adj);CHKERRQ(ierr);
212 
213       ierr = ISCreateGeneral(comm, ncrs_eq, newproc_idx, PETSC_COPY_VALUES, &is_eq_newproc);CHKERRQ(ierr);
214       ierr = PetscFree(newproc_idx);CHKERRQ(ierr);
215     } else { /* simple aggreagtion of parts -- 'is_eq_newproc' */
216       PetscInt rfactor,targetPE;
217 
218       /* find factor */
219       if (new_size == 1) rfactor = size; /* easy */
220       else {
221         PetscReal best_fact = 0.;
222         jj = -1;
223         for (kk = 1 ; kk <= size ; kk++) {
224           if (!(size%kk)) { /* a candidate */
225             PetscReal nactpe = (PetscReal)size/(PetscReal)kk, fact = nactpe/(PetscReal)new_size;
226             if (fact > 1.0) fact = 1./fact; /* keep fact < 1 */
227             if (fact > best_fact) {
228               best_fact = fact; jj = kk;
229             }
230           }
231         }
232         if (jj != -1) rfactor = jj;
233         else rfactor = 1; /* does this happen .. a prime */
234       }
235       new_size = size/rfactor;
236 
237       if (new_size==nactive) {
238         *a_Amat_crs = Cmat; /* output - no repartitioning or reduction, bail out because nested here */
239         ierr        = PetscFree(counts);CHKERRQ(ierr);
240         ierr = PetscInfo2(pc,"Aggregate processors noop: new_size=%D, neq(loc)=%D\n",new_size,ncrs_eq);CHKERRQ(ierr);
241 #if defined PETSC_GAMG_USE_LOG
242         ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET12],0,0,0,0);CHKERRQ(ierr);
243 #endif
244         PetscFunctionReturn(0);
245       }
246 
247       ierr = PetscInfo1(pc,"Number of equations (loc) %D with simple aggregation\n",ncrs_eq);CHKERRQ(ierr);
248       targetPE = rank/rfactor;
249       ierr     = ISCreateStride(comm, ncrs_eq, targetPE, 0, &is_eq_newproc);CHKERRQ(ierr);
250     } /* end simple 'is_eq_newproc' */
251 
252     /*
253      Create an index set from the is_eq_newproc index set to indicate the mapping TO
254      */
255     ierr = ISPartitioningToNumbering(is_eq_newproc, &is_eq_num);CHKERRQ(ierr);
256     is_eq_num_prim = is_eq_num;
257     /*
258       Determine how many equations/vertices are assigned to each processor
259      */
260     ierr        = ISPartitioningCount(is_eq_newproc, size, counts);CHKERRQ(ierr);
261     ncrs_eq_new = counts[rank];
262     ierr        = ISDestroy(&is_eq_newproc);CHKERRQ(ierr);
263     ncrs_new = ncrs_eq_new/cr_bs; /* eqs */
264 
265     ierr = PetscFree(counts);CHKERRQ(ierr);
266 #if defined PETSC_GAMG_USE_LOG
267     ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET12],0,0,0,0);CHKERRQ(ierr);
268 #endif
269     /* data movement scope -- this could be moved to subclasses so that we don't try to cram all auxilary data into some complex abstracted thing */
270     {
271     Vec            src_crd, dest_crd;
272     const PetscInt *idx,ndata_rows=pc_gamg->data_cell_rows,ndata_cols=pc_gamg->data_cell_cols,node_data_sz=ndata_rows*ndata_cols;
273     VecScatter     vecscat;
274     PetscScalar    *array;
275     IS isscat;
276 
277     /* move data (for primal equations only) */
278     /* Create a vector to contain the newly ordered element information */
279     ierr = VecCreate(comm, &dest_crd);CHKERRQ(ierr);
280     ierr = VecSetSizes(dest_crd, node_data_sz*ncrs_new, PETSC_DECIDE);CHKERRQ(ierr);
281     ierr = VecSetType(dest_crd,VECSTANDARD);CHKERRQ(ierr); /* this is needed! */
282     /*
283      There are 'ndata_rows*ndata_cols' data items per node, (one can think of the vectors of having
284      a block size of ...).  Note, ISs are expanded into equation space by 'cr_bs'.
285      */
286     ierr = PetscMalloc1(ncrs*node_data_sz, &tidx);CHKERRQ(ierr);
287     ierr = ISGetIndices(is_eq_num_prim, &idx);CHKERRQ(ierr);
288     for (ii=0,jj=0; ii<ncrs; ii++) {
289       PetscInt id = idx[ii*cr_bs]/cr_bs; /* get node back */
290       for (kk=0; kk<node_data_sz; kk++, jj++) tidx[jj] = id*node_data_sz + kk;
291     }
292     ierr = ISRestoreIndices(is_eq_num_prim, &idx);CHKERRQ(ierr);
293     ierr = ISCreateGeneral(comm, node_data_sz*ncrs, tidx, PETSC_COPY_VALUES, &isscat);CHKERRQ(ierr);
294     ierr = PetscFree(tidx);CHKERRQ(ierr);
295     /*
296      Create a vector to contain the original vertex information for each element
297      */
298     ierr = VecCreateSeq(PETSC_COMM_SELF, node_data_sz*ncrs, &src_crd);CHKERRQ(ierr);
299     for (jj=0; jj<ndata_cols; jj++) {
300       const PetscInt stride0=ncrs*pc_gamg->data_cell_rows;
301       for (ii=0; ii<ncrs; ii++) {
302         for (kk=0; kk<ndata_rows; kk++) {
303           PetscInt    ix = ii*ndata_rows + kk + jj*stride0, jx = ii*node_data_sz + kk*ndata_cols + jj;
304           PetscScalar tt = (PetscScalar)pc_gamg->data[ix];
305           ierr = VecSetValues(src_crd, 1, &jx, &tt, INSERT_VALUES);CHKERRQ(ierr);
306         }
307       }
308     }
309     ierr = VecAssemblyBegin(src_crd);CHKERRQ(ierr);
310     ierr = VecAssemblyEnd(src_crd);CHKERRQ(ierr);
311     /*
312       Scatter the element vertex information (still in the original vertex ordering)
313       to the correct processor
314     */
315     ierr = VecScatterCreate(src_crd, NULL, dest_crd, isscat, &vecscat);CHKERRQ(ierr);
316     ierr = ISDestroy(&isscat);CHKERRQ(ierr);
317     ierr = VecScatterBegin(vecscat,src_crd,dest_crd,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
318     ierr = VecScatterEnd(vecscat,src_crd,dest_crd,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr);
319     ierr = VecScatterDestroy(&vecscat);CHKERRQ(ierr);
320     ierr = VecDestroy(&src_crd);CHKERRQ(ierr);
321     /*
322       Put the element vertex data into a new allocation of the gdata->ele
323     */
324     ierr = PetscFree(pc_gamg->data);CHKERRQ(ierr);
325     ierr = PetscMalloc1(node_data_sz*ncrs_new, &pc_gamg->data);CHKERRQ(ierr);
326 
327     pc_gamg->data_sz = node_data_sz*ncrs_new;
328     strideNew        = ncrs_new*ndata_rows;
329 
330     ierr = VecGetArray(dest_crd, &array);CHKERRQ(ierr);
331     for (jj=0; jj<ndata_cols; jj++) {
332       for (ii=0; ii<ncrs_new; ii++) {
333         for (kk=0; kk<ndata_rows; kk++) {
334           PetscInt ix = ii*ndata_rows + kk + jj*strideNew, jx = ii*node_data_sz + kk*ndata_cols + jj;
335           pc_gamg->data[ix] = PetscRealPart(array[jx]);
336         }
337       }
338     }
339     ierr = VecRestoreArray(dest_crd, &array);CHKERRQ(ierr);
340     ierr = VecDestroy(&dest_crd);CHKERRQ(ierr);
341     }
342     /* move A and P (columns) with new layout */
343 #if defined PETSC_GAMG_USE_LOG
344     ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET13],0,0,0,0);CHKERRQ(ierr);
345 #endif
346 
347     /*
348       Invert for MatGetSubMatrix
349     */
350     ierr = ISInvertPermutation(is_eq_num, ncrs_eq_new, &new_eq_indices);CHKERRQ(ierr);
351     ierr = ISSort(new_eq_indices);CHKERRQ(ierr); /* is this needed? */
352     ierr = ISSetBlockSize(new_eq_indices, cr_bs);CHKERRQ(ierr);
353     if (is_eq_num != is_eq_num_prim) {
354       ierr = ISDestroy(&is_eq_num_prim);CHKERRQ(ierr); /* could be same as 'is_eq_num' */
355     }
356     if (Pcolumnperm) {
357       ierr = PetscObjectReference((PetscObject)new_eq_indices);CHKERRQ(ierr);
358       *Pcolumnperm = new_eq_indices;
359     }
360     ierr = ISDestroy(&is_eq_num);CHKERRQ(ierr);
361 #if defined PETSC_GAMG_USE_LOG
362     ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET13],0,0,0,0);CHKERRQ(ierr);
363     ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET14],0,0,0,0);CHKERRQ(ierr);
364 #endif
365     /* 'a_Amat_crs' output */
366     {
367       Mat mat;
368       ierr        = MatGetSubMatrix(Cmat, new_eq_indices, new_eq_indices, MAT_INITIAL_MATRIX, &mat);CHKERRQ(ierr);
369       *a_Amat_crs = mat;
370     }
371     ierr = MatDestroy(&Cmat);CHKERRQ(ierr);
372 
373 #if defined PETSC_GAMG_USE_LOG
374     ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET14],0,0,0,0);CHKERRQ(ierr);
375 #endif
376     /* prolongator */
377     {
378       IS       findices;
379       PetscInt Istart,Iend;
380       Mat      Pnew;
381 
382       ierr = MatGetOwnershipRange(Pold, &Istart, &Iend);CHKERRQ(ierr);
383 #if defined PETSC_GAMG_USE_LOG
384       ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET15],0,0,0,0);CHKERRQ(ierr);
385 #endif
386       ierr = ISCreateStride(comm,Iend-Istart,Istart,1,&findices);CHKERRQ(ierr);
387       ierr = ISSetBlockSize(findices,f_bs);CHKERRQ(ierr);
388       ierr = MatGetSubMatrix(Pold, findices, new_eq_indices, MAT_INITIAL_MATRIX, &Pnew);CHKERRQ(ierr);
389       ierr = ISDestroy(&findices);CHKERRQ(ierr);
390 
391 #if defined PETSC_GAMG_USE_LOG
392       ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET15],0,0,0,0);CHKERRQ(ierr);
393 #endif
394       ierr = MatDestroy(a_P_inout);CHKERRQ(ierr);
395 
396       /* output - repartitioned */
397       *a_P_inout = Pnew;
398     }
399     ierr = ISDestroy(&new_eq_indices);CHKERRQ(ierr);
400 
401     *a_nactive_proc = new_size; /* output */
402   }
403   PetscFunctionReturn(0);
404 }
405 
406 /* -------------------------------------------------------------------------- */
407 /*
408    PCSetUp_GAMG - Prepares for the use of the GAMG preconditioner
409                     by setting data structures and options.
410 
411    Input Parameter:
412 .  pc - the preconditioner context
413 
414 */
415 #undef __FUNCT__
416 #define __FUNCT__ "PCSetUp_GAMG"
417 PetscErrorCode PCSetUp_GAMG(PC pc)
418 {
419   PetscErrorCode ierr;
420   PC_MG          *mg      = (PC_MG*)pc->data;
421   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
422   Mat            Pmat     = pc->pmat;
423   PetscInt       fine_level,level,level1,bs,M,qq,lidx,nASMBlocksArr[GAMG_MAXLEVELS];
424   MPI_Comm       comm;
425   PetscMPIInt    rank,size,nactivepe;
426   Mat            Aarr[GAMG_MAXLEVELS],Parr[GAMG_MAXLEVELS];
427   IS             *ASMLocalIDsArr[GAMG_MAXLEVELS];
428   PetscLogDouble nnz0=0.,nnztot=0.;
429   MatInfo        info;
430 
431   PetscFunctionBegin;
432   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
433   ierr = MPI_Comm_rank(comm,&rank);CHKERRQ(ierr);
434   ierr = MPI_Comm_size(comm,&size);CHKERRQ(ierr);
435 
436   if (pc_gamg->setup_count++ > 0) {
437     if ((PetscBool)(!pc_gamg->reuse_prol)) {
438       /* reset everything */
439       ierr = PCReset_MG(pc);CHKERRQ(ierr);
440       pc->setupcalled = 0;
441     } else {
442       PC_MG_Levels **mglevels = mg->levels;
443       /* just do Galerkin grids */
444       Mat          B,dA,dB;
445 
446      if (!pc->setupcalled) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_PLIB,"PCSetUp() has not been called yet");
447       if (pc_gamg->Nlevels > 1) {
448         /* currently only handle case where mat and pmat are the same on coarser levels */
449         ierr = KSPGetOperators(mglevels[pc_gamg->Nlevels-1]->smoothd,&dA,&dB);CHKERRQ(ierr);
450         /* (re)set to get dirty flag */
451         ierr = KSPSetOperators(mglevels[pc_gamg->Nlevels-1]->smoothd,dA,dB);CHKERRQ(ierr);
452 
453         for (level=pc_gamg->Nlevels-2; level>=0; level--) {
454           /* the first time through the matrix structure has changed from repartitioning */
455           if (pc_gamg->setup_count==2) {
456             ierr = MatPtAP(dB,mglevels[level+1]->interpolate,MAT_INITIAL_MATRIX,1.0,&B);CHKERRQ(ierr);
457             ierr = MatDestroy(&mglevels[level]->A);CHKERRQ(ierr);
458 
459             mglevels[level]->A = B;
460           } else {
461             ierr = KSPGetOperators(mglevels[level]->smoothd,NULL,&B);CHKERRQ(ierr);
462             ierr = MatPtAP(dB,mglevels[level+1]->interpolate,MAT_REUSE_MATRIX,1.0,&B);CHKERRQ(ierr);
463           }
464           ierr = KSPSetOperators(mglevels[level]->smoothd,B,B);CHKERRQ(ierr);
465           dB   = B;
466         }
467       }
468 
469       ierr = PCSetUp_MG(pc);CHKERRQ(ierr);
470       PetscFunctionReturn(0);
471     }
472   }
473 
474   if (!pc_gamg->data) {
475     if (pc_gamg->orig_data) {
476       ierr = MatGetBlockSize(Pmat, &bs);CHKERRQ(ierr);
477       ierr = MatGetLocalSize(Pmat, &qq, NULL);CHKERRQ(ierr);
478 
479       pc_gamg->data_sz        = (qq/bs)*pc_gamg->orig_data_cell_rows*pc_gamg->orig_data_cell_cols;
480       pc_gamg->data_cell_rows = pc_gamg->orig_data_cell_rows;
481       pc_gamg->data_cell_cols = pc_gamg->orig_data_cell_cols;
482 
483       ierr = PetscMalloc1(pc_gamg->data_sz, &pc_gamg->data);CHKERRQ(ierr);
484       for (qq=0; qq<pc_gamg->data_sz; qq++) pc_gamg->data[qq] = pc_gamg->orig_data[qq];
485     } else {
486       if (!pc_gamg->ops->createdefaultdata) SETERRQ(comm,PETSC_ERR_PLIB,"'createdefaultdata' not set(?) need to support NULL data");
487       ierr = pc_gamg->ops->createdefaultdata(pc,Pmat);CHKERRQ(ierr);
488     }
489   }
490 
491   /* cache original data for reuse */
492   if (!pc_gamg->orig_data && (PetscBool)(!pc_gamg->reuse_prol)) {
493     ierr = PetscMalloc1(pc_gamg->data_sz, &pc_gamg->orig_data);CHKERRQ(ierr);
494     for (qq=0; qq<pc_gamg->data_sz; qq++) pc_gamg->orig_data[qq] = pc_gamg->data[qq];
495     pc_gamg->orig_data_cell_rows = pc_gamg->data_cell_rows;
496     pc_gamg->orig_data_cell_cols = pc_gamg->data_cell_cols;
497   }
498 
499   /* get basic dims */
500   ierr = MatGetBlockSize(Pmat, &bs);CHKERRQ(ierr);
501   ierr = MatGetSize(Pmat, &M, &qq);CHKERRQ(ierr);
502 
503   ierr = MatGetInfo(Pmat,MAT_GLOBAL_SUM,&info);CHKERRQ(ierr); /* global reduction */
504   nnz0   = info.nz_used;
505   nnztot = info.nz_used;
506   ierr = PetscInfo6(pc,"level %d) N=%D, n data rows=%d, n data cols=%d, nnz/row (ave)=%d, np=%d\n",0,M,pc_gamg->data_cell_rows,pc_gamg->data_cell_cols,(int)(nnz0/(PetscReal)M+0.5),size);CHKERRQ(ierr);
507 
508   /* Get A_i and R_i */
509   for (level=0, Aarr[0]=Pmat, nactivepe = size; level < (pc_gamg->Nlevels-1) && (!level || M>pc_gamg->coarse_eq_limit); level++) {
510     pc_gamg->current_level = level;
511     level1 = level + 1;
512 #if defined PETSC_GAMG_USE_LOG
513     ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET1],0,0,0,0);CHKERRQ(ierr);
514 #if (defined GAMG_STAGES)
515     ierr = PetscLogStagePush(gamg_stages[level]);CHKERRQ(ierr);
516 #endif
517 #endif
518     { /* construct prolongator */
519       Mat              Gmat;
520       PetscCoarsenData *agg_lists;
521       Mat              Prol11;
522 
523       ierr = pc_gamg->ops->graph(pc,Aarr[level], &Gmat);CHKERRQ(ierr);
524       ierr = pc_gamg->ops->coarsen(pc, &Gmat, &agg_lists);CHKERRQ(ierr);
525       ierr = pc_gamg->ops->prolongator(pc,Aarr[level],Gmat,agg_lists,&Prol11);CHKERRQ(ierr);
526 
527       /* could have failed to create new level */
528       if (Prol11) {
529         /* get new block size of coarse matrices */
530         ierr = MatGetBlockSizes(Prol11, NULL, &bs);CHKERRQ(ierr);
531 
532         if (pc_gamg->ops->optprolongator) {
533           /* smooth */
534           ierr = pc_gamg->ops->optprolongator(pc, Aarr[level], &Prol11);CHKERRQ(ierr);
535         }
536 
537         Parr[level1] = Prol11;
538       } else Parr[level1] = NULL;
539 
540       if (pc_gamg->use_aggs_in_asm) {
541         PetscInt bs;
542         ierr = MatGetBlockSizes(Prol11, &bs, NULL);CHKERRQ(ierr);
543         ierr = PetscCDGetASMBlocks(agg_lists, bs, &nASMBlocksArr[level], &ASMLocalIDsArr[level]);CHKERRQ(ierr);
544       }
545 
546       ierr = MatDestroy(&Gmat);CHKERRQ(ierr);
547       ierr = PetscCDDestroy(agg_lists);CHKERRQ(ierr);
548     } /* construct prolongator scope */
549 #if defined PETSC_GAMG_USE_LOG
550     ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET1],0,0,0,0);CHKERRQ(ierr);
551 #endif
552     if (!level) Aarr[0] = Pmat; /* use Pmat for finest level setup */
553     if (!Parr[level1]) {
554       ierr =  PetscInfo1(pc,"Stop gridding, level %D\n",level);CHKERRQ(ierr);
555 #if defined PETSC_GAMG_USE_LOG && defined GAMG_STAGES
556       ierr = PetscLogStagePop();CHKERRQ(ierr);
557 #endif
558       break;
559     }
560 #if defined PETSC_GAMG_USE_LOG
561     ierr = PetscLogEventBegin(petsc_gamg_setup_events[SET2],0,0,0,0);CHKERRQ(ierr);
562 #endif
563 
564     ierr = pc_gamg->ops->createlevel(pc, Aarr[level], bs,&Parr[level1], &Aarr[level1], &nactivepe, NULL);CHKERRQ(ierr);
565 
566 #if defined PETSC_GAMG_USE_LOG
567     ierr = PetscLogEventEnd(petsc_gamg_setup_events[SET2],0,0,0,0);CHKERRQ(ierr);
568 #endif
569     ierr = MatGetSize(Aarr[level1], &M, &qq);CHKERRQ(ierr);
570     ierr = MatGetInfo(Aarr[level1], MAT_GLOBAL_SUM, &info);CHKERRQ(ierr);
571     nnztot += info.nz_used;
572     ierr = PetscInfo5(pc,"%d) N=%D, n data cols=%d, nnz/row (ave)=%d, %d active pes\n",level1,M,pc_gamg->data_cell_cols,(int)(info.nz_used/(PetscReal)M),nactivepe);CHKERRQ(ierr);
573 
574 #if (defined PETSC_GAMG_USE_LOG && defined GAMG_STAGES)
575     ierr = PetscLogStagePop();CHKERRQ(ierr);
576 #endif
577     /* stop if one node or one proc -- could pull back for singular problems */
578     if ( (pc_gamg->data_cell_cols && M/pc_gamg->data_cell_cols < 2) || (!pc_gamg->data_cell_cols && M/bs < 2) ) {
579       ierr =  PetscInfo2(pc,"HARD stop of coarsening on level %D.  Grid too small: %D block nodes\n",level,M/bs);CHKERRQ(ierr);
580       level++;
581       break;
582     }
583   } /* levels */
584   ierr                  = PetscFree(pc_gamg->data);CHKERRQ(ierr);
585 
586   ierr = PetscInfo2(pc,"%D levels, grid complexity = %g\n",level+1,nnztot/nnz0);CHKERRQ(ierr);
587   pc_gamg->Nlevels = level + 1;
588   fine_level       = level;
589   ierr             = PCMGSetLevels(pc,pc_gamg->Nlevels,NULL);CHKERRQ(ierr);
590 
591   if (pc_gamg->Nlevels > 1) { /* don't setup MG if one level */
592     /* set default smoothers & set operators */
593     for (lidx = 1, level = pc_gamg->Nlevels-2; lidx <= fine_level; lidx++, level--) {
594       KSP smoother;
595       PC  subpc;
596 
597       ierr = PCMGGetSmoother(pc, lidx, &smoother);CHKERRQ(ierr);
598       ierr = KSPGetPC(smoother, &subpc);CHKERRQ(ierr);
599 
600       ierr = KSPSetNormType(smoother, KSP_NORM_NONE);CHKERRQ(ierr);
601       /* set ops */
602       ierr = KSPSetOperators(smoother, Aarr[level], Aarr[level]);CHKERRQ(ierr);
603       ierr = PCMGSetInterpolation(pc, lidx, Parr[level+1]);CHKERRQ(ierr);
604 
605       /* set defaults */
606       ierr = KSPSetType(smoother, KSPCHEBYSHEV);CHKERRQ(ierr);
607 
608       /* set blocks for ASM smoother that uses the 'aggregates' */
609       if (pc_gamg->use_aggs_in_asm) {
610         PetscInt sz;
611         IS       *is;
612 
613         sz   = nASMBlocksArr[level];
614         is   = ASMLocalIDsArr[level];
615         ierr = PCSetType(subpc, PCASM);CHKERRQ(ierr);
616         ierr = PCASMSetOverlap(subpc, 1);CHKERRQ(ierr);
617         ierr = PCASMSetType(subpc,PC_ASM_BASIC);CHKERRQ(ierr);
618         if (!sz) {
619           IS       is;
620           PetscInt my0,kk;
621           ierr = MatGetOwnershipRange(Aarr[level], &my0, &kk);CHKERRQ(ierr);
622           ierr = ISCreateGeneral(PETSC_COMM_SELF, 1, &my0, PETSC_COPY_VALUES, &is);CHKERRQ(ierr);
623           ierr = PCASMSetLocalSubdomains(subpc, 1, NULL, &is);CHKERRQ(ierr);
624           ierr = ISDestroy(&is);CHKERRQ(ierr);
625         } else {
626           PetscInt kk;
627           ierr = PCASMSetLocalSubdomains(subpc, sz, NULL, is);CHKERRQ(ierr);
628           for (kk=0; kk<sz; kk++) {
629             ierr = ISDestroy(&is[kk]);CHKERRQ(ierr);
630           }
631           ierr = PetscFree(is);CHKERRQ(ierr);
632         }
633         ASMLocalIDsArr[level] = NULL;
634         nASMBlocksArr[level]  = 0;
635       } else {
636         ierr = PCSetType(subpc, PCSOR);CHKERRQ(ierr);
637       }
638     }
639     {
640       /* coarse grid */
641       KSP smoother,*k2; PC subpc,pc2; PetscInt ii,first;
642       Mat Lmat = Aarr[(level=pc_gamg->Nlevels-1)]; lidx = 0;
643       ierr = PCMGGetSmoother(pc, lidx, &smoother);CHKERRQ(ierr);
644       ierr = KSPSetOperators(smoother, Lmat, Lmat);CHKERRQ(ierr);
645       ierr = KSPSetNormType(smoother, KSP_NORM_NONE);CHKERRQ(ierr);
646       ierr = KSPGetPC(smoother, &subpc);CHKERRQ(ierr);
647       ierr = PCSetType(subpc, PCBJACOBI);CHKERRQ(ierr);
648       ierr = PCSetUp(subpc);CHKERRQ(ierr);
649       ierr = PCBJacobiGetSubKSP(subpc,&ii,&first,&k2);CHKERRQ(ierr);
650       if (ii != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ii %D is not one",ii);
651       ierr = KSPGetPC(k2[0],&pc2);CHKERRQ(ierr);
652       ierr = PCSetType(pc2, PCLU);CHKERRQ(ierr);
653       ierr = PCFactorSetShiftType(pc2,MAT_SHIFT_INBLOCKS);CHKERRQ(ierr);
654       ierr = KSPSetTolerances(k2[0],PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,1);CHKERRQ(ierr);
655       ierr = KSPSetType(k2[0], KSPPREONLY);CHKERRQ(ierr);
656       /* This flag gets reset by PCBJacobiGetSubKSP(), but our BJacobi really does the same algorithm everywhere (and in
657        * fact, all but one process will have zero dofs), so we reset the flag to avoid having PCView_BJacobi attempt to
658        * view every subdomain as though they were different. */
659       ((PC_BJacobi*)subpc->data)->same_local_solves = PETSC_TRUE;
660     }
661 
662     /* should be called in PCSetFromOptions_GAMG(), but cannot be called prior to PCMGSetLevels() */
663     ierr = PetscObjectOptionsBegin((PetscObject)pc);CHKERRQ(ierr);
664     ierr = PCSetFromOptions_MG(PetscOptionsObject,pc);CHKERRQ(ierr);
665     ierr = PetscOptionsEnd();CHKERRQ(ierr);
666     if (!mg->galerkin) SETERRQ(comm,PETSC_ERR_USER,"PCGAMG must use Galerkin for coarse operators.");
667     if (mg->galerkin == 1) mg->galerkin = 2;
668 
669     /* clean up */
670     for (level=1; level<pc_gamg->Nlevels; level++) {
671       ierr = MatDestroy(&Parr[level]);CHKERRQ(ierr);
672       ierr = MatDestroy(&Aarr[level]);CHKERRQ(ierr);
673     }
674     ierr = PCSetUp_MG(pc);CHKERRQ(ierr);
675   } else {
676     KSP smoother;
677     ierr = PetscInfo(pc,"One level solver used (system is seen as DD). Using default solver.\n");CHKERRQ(ierr);
678     ierr = PCMGGetSmoother(pc, 0, &smoother);CHKERRQ(ierr);
679     ierr = KSPSetOperators(smoother, Aarr[0], Aarr[0]);CHKERRQ(ierr);
680     ierr = KSPSetType(smoother, KSPPREONLY);CHKERRQ(ierr);
681     ierr = PCSetUp_MG(pc);CHKERRQ(ierr);
682   }
683   PetscFunctionReturn(0);
684 }
685 
686 /* ------------------------------------------------------------------------- */
687 /*
688  PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner
689    that was created with PCCreate_GAMG().
690 
691    Input Parameter:
692 .  pc - the preconditioner context
693 
694    Application Interface Routine: PCDestroy()
695 */
696 #undef __FUNCT__
697 #define __FUNCT__ "PCDestroy_GAMG"
698 PetscErrorCode PCDestroy_GAMG(PC pc)
699 {
700   PetscErrorCode ierr;
701   PC_MG          *mg     = (PC_MG*)pc->data;
702   PC_GAMG        *pc_gamg= (PC_GAMG*)mg->innerctx;
703 
704   PetscFunctionBegin;
705   ierr = PCReset_GAMG(pc);CHKERRQ(ierr);
706   if (pc_gamg->ops->destroy) {
707     ierr = (*pc_gamg->ops->destroy)(pc);CHKERRQ(ierr);
708   }
709   ierr = PetscRandomDestroy(&pc_gamg->random);CHKERRQ(ierr);
710   ierr = PetscFree(pc_gamg->ops);CHKERRQ(ierr);
711   ierr = PetscFree(pc_gamg->gamg_type_name);CHKERRQ(ierr);
712   ierr = PetscFree(pc_gamg);CHKERRQ(ierr);
713   ierr = PCDestroy_MG(pc);CHKERRQ(ierr);
714   PetscFunctionReturn(0);
715 }
716 
717 #undef __FUNCT__
718 #define __FUNCT__ "PCGAMGSetProcEqLim"
719 /*@
720    PCGAMGSetProcEqLim - Set number of equations to aim for per process on the coarse grids via processor reduction.
721 
722    Logically Collective on PC
723 
724    Input Parameters:
725 +  pc - the preconditioner context
726 -  n - the number of equations
727 
728 
729    Options Database Key:
730 .  -pc_gamg_process_eq_limit <limit>
731 
732    Notes: GAMG will reduce the number of MPI processes used directly on the coarse grids so that there are around <limit> equations on each process
733           that has degrees of freedom
734 
735    Level: intermediate
736 
737    Concepts: Unstructured multigrid preconditioner
738 
739 .seealso: PCGAMGSetCoarseEqLim()
740 @*/
741 PetscErrorCode  PCGAMGSetProcEqLim(PC pc, PetscInt n)
742 {
743   PetscErrorCode ierr;
744 
745   PetscFunctionBegin;
746   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
747   ierr = PetscTryMethod(pc,"PCGAMGSetProcEqLim_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr);
748   PetscFunctionReturn(0);
749 }
750 
751 #undef __FUNCT__
752 #define __FUNCT__ "PCGAMGSetProcEqLim_GAMG"
753 static PetscErrorCode PCGAMGSetProcEqLim_GAMG(PC pc, PetscInt n)
754 {
755   PC_MG   *mg      = (PC_MG*)pc->data;
756   PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
757 
758   PetscFunctionBegin;
759   if (n>0) pc_gamg->min_eq_proc = n;
760   PetscFunctionReturn(0);
761 }
762 
763 #undef __FUNCT__
764 #define __FUNCT__ "PCGAMGSetCoarseEqLim"
765 /*@
766    PCGAMGSetCoarseEqLim - Set maximum number of equations on coarsest grid.
767 
768  Collective on PC
769 
770    Input Parameters:
771 +  pc - the preconditioner context
772 -  n - maximum number of equations to aim for
773 
774    Options Database Key:
775 .  -pc_gamg_coarse_eq_limit <limit>
776 
777    Level: intermediate
778 
779    Concepts: Unstructured multigrid preconditioner
780 
781 .seealso: PCGAMGSetProcEqLim()
782 @*/
783 PetscErrorCode PCGAMGSetCoarseEqLim(PC pc, PetscInt n)
784 {
785   PetscErrorCode ierr;
786 
787   PetscFunctionBegin;
788   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
789   ierr = PetscTryMethod(pc,"PCGAMGSetCoarseEqLim_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr);
790   PetscFunctionReturn(0);
791 }
792 
793 #undef __FUNCT__
794 #define __FUNCT__ "PCGAMGSetCoarseEqLim_GAMG"
795 static PetscErrorCode PCGAMGSetCoarseEqLim_GAMG(PC pc, PetscInt n)
796 {
797   PC_MG   *mg      = (PC_MG*)pc->data;
798   PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
799 
800   PetscFunctionBegin;
801   if (n>0) pc_gamg->coarse_eq_limit = n;
802   PetscFunctionReturn(0);
803 }
804 
805 #undef __FUNCT__
806 #define __FUNCT__ "PCGAMGSetRepartition"
807 /*@
808    PCGAMGSetRepartition - Repartition the degrees of freedom across the processors on the coarser grids
809 
810    Collective on PC
811 
812    Input Parameters:
813 +  pc - the preconditioner context
814 -  n - PETSC_TRUE or PETSC_FALSE
815 
816    Options Database Key:
817 .  -pc_gamg_repartition <true,false>
818 
819    Notes: this will generally improve the loading balancing of the work on each level
820 
821    Level: intermediate
822 
823    Concepts: Unstructured multigrid preconditioner
824 
825 .seealso: ()
826 @*/
827 PetscErrorCode PCGAMGSetRepartition(PC pc, PetscBool n)
828 {
829   PetscErrorCode ierr;
830 
831   PetscFunctionBegin;
832   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
833   ierr = PetscTryMethod(pc,"PCGAMGSetRepartition_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr);
834   PetscFunctionReturn(0);
835 }
836 
837 #undef __FUNCT__
838 #define __FUNCT__ "PCGAMGSetRepartition_GAMG"
839 static PetscErrorCode PCGAMGSetRepartition_GAMG(PC pc, PetscBool n)
840 {
841   PC_MG   *mg      = (PC_MG*)pc->data;
842   PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
843 
844   PetscFunctionBegin;
845   pc_gamg->repart = n;
846   PetscFunctionReturn(0);
847 }
848 
849 #undef __FUNCT__
850 #define __FUNCT__ "PCGAMGSetReuseInterpolation"
851 /*@
852    PCGAMGSetReuseInterpolation - Reuse prolongation when rebuilding algebraic multigrid preconditioner
853 
854    Collective on PC
855 
856    Input Parameters:
857 +  pc - the preconditioner context
858 -  n - PETSC_TRUE or PETSC_FALSE
859 
860    Options Database Key:
861 .  -pc_gamg_reuse_interpolation <true,false>
862 
863    Level: intermediate
864 
865    Notes: this may negatively affect the convergence rate of the method on new matrices if the matrix entries change a great deal, but allows
866           rebuilding the preconditioner quicker.
867 
868    Concepts: Unstructured multigrid preconditioner
869 
870 .seealso: ()
871 @*/
872 PetscErrorCode PCGAMGSetReuseInterpolation(PC pc, PetscBool n)
873 {
874   PetscErrorCode ierr;
875 
876   PetscFunctionBegin;
877   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
878   ierr = PetscTryMethod(pc,"PCGAMGSetReuseInterpolation_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr);
879   PetscFunctionReturn(0);
880 }
881 
882 #undef __FUNCT__
883 #define __FUNCT__ "PCGAMGSetReuseInterpolation_GAMG"
884 static PetscErrorCode PCGAMGSetReuseInterpolation_GAMG(PC pc, PetscBool n)
885 {
886   PC_MG   *mg      = (PC_MG*)pc->data;
887   PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
888 
889   PetscFunctionBegin;
890   pc_gamg->reuse_prol = n;
891   PetscFunctionReturn(0);
892 }
893 
894 #undef __FUNCT__
895 #define __FUNCT__ "PCGAMGASMSetUseAggs"
896 /*@
897    PCGAMGASMSetUseAggs - Have the PCGAMG smoother on each level use the aggregates defined by the coarsening process as the subdomains for the additive Schwarz preconditioner.
898 
899    Collective on PC
900 
901    Input Parameters:
902 +  pc - the preconditioner context
903 -  flg - PETSC_TRUE to use aggregates, PETSC_FALSE to not
904 
905    Options Database Key:
906 .  -pc_gamg_asm_use_agg
907 
908    Level: intermediate
909 
910    Concepts: Unstructured multigrid preconditioner
911 
912 .seealso: ()
913 @*/
914 PetscErrorCode PCGAMGASMSetUseAggs(PC pc, PetscBool flg)
915 {
916   PetscErrorCode ierr;
917 
918   PetscFunctionBegin;
919   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
920   ierr = PetscTryMethod(pc,"PCGAMGASMSetUseAggs_C",(PC,PetscBool),(pc,flg));CHKERRQ(ierr);
921   PetscFunctionReturn(0);
922 }
923 
924 #undef __FUNCT__
925 #define __FUNCT__ "PCGAMGASMSetUseAggs_GAMG"
926 static PetscErrorCode PCGAMGASMSetUseAggs_GAMG(PC pc, PetscBool flg)
927 {
928   PC_MG   *mg      = (PC_MG*)pc->data;
929   PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
930 
931   PetscFunctionBegin;
932   pc_gamg->use_aggs_in_asm = flg;
933   PetscFunctionReturn(0);
934 }
935 
936 #undef __FUNCT__
937 #define __FUNCT__ "PCGAMGSetNlevels"
938 /*@
939    PCGAMGSetNlevels -  Sets the maximum number of levels PCGAMG will use
940 
941    Not collective on PC
942 
943    Input Parameters:
944 +  pc - the preconditioner
945 -  n - the maximum number of levels to use
946 
947    Options Database Key:
948 .  -pc_mg_levels
949 
950    Level: intermediate
951 
952    Concepts: Unstructured multigrid preconditioner
953 
954 .seealso: ()
955 @*/
956 PetscErrorCode PCGAMGSetNlevels(PC pc, PetscInt n)
957 {
958   PetscErrorCode ierr;
959 
960   PetscFunctionBegin;
961   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
962   ierr = PetscTryMethod(pc,"PCGAMGSetNlevels_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr);
963   PetscFunctionReturn(0);
964 }
965 
966 #undef __FUNCT__
967 #define __FUNCT__ "PCGAMGSetNlevels_GAMG"
968 static PetscErrorCode PCGAMGSetNlevels_GAMG(PC pc, PetscInt n)
969 {
970   PC_MG   *mg      = (PC_MG*)pc->data;
971   PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
972 
973   PetscFunctionBegin;
974   pc_gamg->Nlevels = n;
975   PetscFunctionReturn(0);
976 }
977 
978 #undef __FUNCT__
979 #define __FUNCT__ "PCGAMGSetThreshold"
980 /*@
981    PCGAMGSetThreshold - Relative threshold to use for dropping edges in aggregation graph
982 
983    Not collective on PC
984 
985    Input Parameters:
986 +  pc - the preconditioner context
987 -  threshold - the threshold value, 0.0 means keep all nonzero entries in the graph; negative means keep even zero entries in the graph
988 
989    Options Database Key:
990 .  -pc_gamg_threshold <threshold>
991 
992    Notes: Before aggregating the graph GAMG will remove small values from the graph thus reducing the coupling in the graph and a different
993     (perhaps better) coarser set of points.
994 
995    Level: intermediate
996 
997    Concepts: Unstructured multigrid preconditioner
998 
999 .seealso: ()
1000 @*/
1001 PetscErrorCode PCGAMGSetThreshold(PC pc, PetscReal n)
1002 {
1003   PetscErrorCode ierr;
1004 
1005   PetscFunctionBegin;
1006   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1007   ierr = PetscTryMethod(pc,"PCGAMGSetThreshold_C",(PC,PetscReal),(pc,n));CHKERRQ(ierr);
1008   PetscFunctionReturn(0);
1009 }
1010 
1011 #undef __FUNCT__
1012 #define __FUNCT__ "PCGAMGSetThreshold_GAMG"
1013 static PetscErrorCode PCGAMGSetThreshold_GAMG(PC pc, PetscReal n)
1014 {
1015   PC_MG   *mg      = (PC_MG*)pc->data;
1016   PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1017 
1018   PetscFunctionBegin;
1019   pc_gamg->threshold = n;
1020   PetscFunctionReturn(0);
1021 }
1022 
1023 #undef __FUNCT__
1024 #define __FUNCT__ "PCGAMGSetType"
1025 /*@
1026    PCGAMGSetType - Set solution method
1027 
1028    Collective on PC
1029 
1030    Input Parameters:
1031 +  pc - the preconditioner context
1032 -  type - PCGAMGAGG, PCGAMGGEO, or PCGAMGCLASSICAL
1033 
1034    Options Database Key:
1035 .  -pc_gamg_type <agg,geo,classical> - type of algebraic multigrid to apply
1036 
1037    Level: intermediate
1038 
1039    Concepts: Unstructured multigrid preconditioner
1040 
1041 .seealso: PCGAMGGetType(), PCGAMG, PCGAMGType
1042 @*/
1043 PetscErrorCode PCGAMGSetType(PC pc, PCGAMGType type)
1044 {
1045   PetscErrorCode ierr;
1046 
1047   PetscFunctionBegin;
1048   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1049   ierr = PetscTryMethod(pc,"PCGAMGSetType_C",(PC,PCGAMGType),(pc,type));CHKERRQ(ierr);
1050   PetscFunctionReturn(0);
1051 }
1052 
1053 #undef __FUNCT__
1054 #define __FUNCT__ "PCGAMGGetType"
1055 /*@
1056    PCGAMGGetType - Get solution method
1057 
1058    Collective on PC
1059 
1060    Input Parameter:
1061 .  pc - the preconditioner context
1062 
1063    Output Parameter:
1064 .  type - the type of algorithm used
1065 
1066    Level: intermediate
1067 
1068    Concepts: Unstructured multigrid preconditioner
1069 
1070 .seealso: PCGAMGSetType(), PCGAMGType
1071 @*/
1072 PetscErrorCode PCGAMGGetType(PC pc, PCGAMGType *type)
1073 {
1074   PetscErrorCode ierr;
1075 
1076   PetscFunctionBegin;
1077   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1078   ierr = PetscUseMethod(pc,"PCGAMGGetType_C",(PC,PCGAMGType*),(pc,type));CHKERRQ(ierr);
1079   PetscFunctionReturn(0);
1080 }
1081 
1082 #undef __FUNCT__
1083 #define __FUNCT__ "PCGAMGGetType_GAMG"
1084 static PetscErrorCode PCGAMGGetType_GAMG(PC pc, PCGAMGType *type)
1085 {
1086   PC_MG          *mg      = (PC_MG*)pc->data;
1087   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
1088 
1089   PetscFunctionBegin;
1090   *type = pc_gamg->type;
1091   PetscFunctionReturn(0);
1092 }
1093 
1094 #undef __FUNCT__
1095 #define __FUNCT__ "PCGAMGSetType_GAMG"
1096 static PetscErrorCode PCGAMGSetType_GAMG(PC pc, PCGAMGType type)
1097 {
1098   PetscErrorCode ierr,(*r)(PC);
1099   PC_MG          *mg      = (PC_MG*)pc->data;
1100   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
1101 
1102   PetscFunctionBegin;
1103   pc_gamg->type = type;
1104   ierr = PetscFunctionListFind(GAMGList,type,&r);CHKERRQ(ierr);
1105   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown GAMG type %s given",type);
1106   if (pc_gamg->ops->destroy) {
1107     ierr = (*pc_gamg->ops->destroy)(pc);CHKERRQ(ierr);
1108     ierr = PetscMemzero(pc_gamg->ops,sizeof(struct _PCGAMGOps));CHKERRQ(ierr);
1109     pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG;
1110     /* cleaning up common data in pc_gamg - this should disapear someday */
1111     pc_gamg->data_cell_cols = 0;
1112     pc_gamg->data_cell_rows = 0;
1113     pc_gamg->orig_data_cell_cols = 0;
1114     pc_gamg->orig_data_cell_rows = 0;
1115     ierr = PetscFree(pc_gamg->data);CHKERRQ(ierr);
1116     pc_gamg->data_sz = 0;
1117   }
1118   ierr = PetscFree(pc_gamg->gamg_type_name);CHKERRQ(ierr);
1119   ierr = PetscStrallocpy(type,&pc_gamg->gamg_type_name);CHKERRQ(ierr);
1120   ierr = (*r)(pc);CHKERRQ(ierr);
1121   PetscFunctionReturn(0);
1122 }
1123 
1124 #undef __FUNCT__
1125 #define __FUNCT__ "PCView_GAMG"
1126 static PetscErrorCode PCView_GAMG(PC pc,PetscViewer viewer)
1127 {
1128   PetscErrorCode ierr;
1129   PC_MG          *mg      = (PC_MG*)pc->data;
1130   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
1131 
1132   PetscFunctionBegin;
1133   ierr = PetscViewerASCIIPrintf(viewer,"    GAMG specific options\n");CHKERRQ(ierr);
1134   ierr = PetscViewerASCIIPrintf(viewer,"      Threshold for dropping small values from graph %g\n",(double)pc_gamg->threshold);CHKERRQ(ierr);
1135   if (pc_gamg->use_aggs_in_asm) {
1136     ierr = PetscViewerASCIIPrintf(viewer,"      Using aggregates from coarsening process to define subdomains for PCASM\n");CHKERRQ(ierr);
1137   }
1138   if (pc_gamg->ops->view) {
1139     ierr = (*pc_gamg->ops->view)(pc,viewer);CHKERRQ(ierr);
1140   }
1141   PetscFunctionReturn(0);
1142 }
1143 
1144 #undef __FUNCT__
1145 #define __FUNCT__ "PCSetFromOptions_GAMG"
1146 PetscErrorCode PCSetFromOptions_GAMG(PetscOptionItems *PetscOptionsObject,PC pc)
1147 {
1148   PetscErrorCode ierr;
1149   PC_MG          *mg      = (PC_MG*)pc->data;
1150   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
1151   PetscBool      flag;
1152   MPI_Comm       comm;
1153   char           prefix[256];
1154   const char     *pcpre;
1155 
1156   PetscFunctionBegin;
1157   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
1158   ierr = PetscOptionsHead(PetscOptionsObject,"GAMG options");CHKERRQ(ierr);
1159   {
1160     char tname[256];
1161     ierr = PetscOptionsFList("-pc_gamg_type","Type of AMG method","PCGAMGSetType",GAMGList, pc_gamg->gamg_type_name, tname, sizeof(tname), &flag);CHKERRQ(ierr);
1162     if (flag) {
1163       ierr = PCGAMGSetType(pc,tname);CHKERRQ(ierr);
1164     }
1165     ierr = PetscOptionsBool("-pc_gamg_repartition","Repartion coarse grids","PCGAMGSetRepartition",pc_gamg->repart,&pc_gamg->repart,NULL);CHKERRQ(ierr);
1166     ierr = PetscOptionsBool("-pc_gamg_reuse_interpolation","Reuse prolongation operator","PCGAMGReuseInterpolation",pc_gamg->reuse_prol,&pc_gamg->reuse_prol,NULL);CHKERRQ(ierr);
1167     ierr = PetscOptionsBool("-pc_gamg_asm_use_agg","Use aggregation agragates for ASM smoother","PCGAMGASMSetUseAggs",pc_gamg->use_aggs_in_asm,&pc_gamg->use_aggs_in_asm,NULL);CHKERRQ(ierr);
1168     ierr = PetscOptionsInt("-pc_gamg_process_eq_limit","Limit (goal) on number of equations per process on coarse grids","PCGAMGSetProcEqLim",pc_gamg->min_eq_proc,&pc_gamg->min_eq_proc,NULL);CHKERRQ(ierr);
1169     ierr = PetscOptionsInt("-pc_gamg_coarse_eq_limit","Limit on number of equations for the coarse grid","PCGAMGSetCoarseEqLim",pc_gamg->coarse_eq_limit,&pc_gamg->coarse_eq_limit,NULL);CHKERRQ(ierr);
1170     ierr = PetscOptionsReal("-pc_gamg_threshold","Relative threshold to use for dropping edges in aggregation graph","PCGAMGSetThreshold",pc_gamg->threshold,&pc_gamg->threshold,&flag);CHKERRQ(ierr);
1171     ierr = PetscOptionsInt("-pc_mg_levels","Set number of MG levels","PCGAMGSetNlevels",pc_gamg->Nlevels,&pc_gamg->Nlevels,NULL);CHKERRQ(ierr);
1172 
1173     /* set options for subtype */
1174     if (pc_gamg->ops->setfromoptions) {ierr = (*pc_gamg->ops->setfromoptions)(PetscOptionsObject,pc);CHKERRQ(ierr);}
1175   }
1176   ierr = PCGetOptionsPrefix(pc, &pcpre);CHKERRQ(ierr);
1177   ierr = PetscSNPrintf(prefix,sizeof(prefix),"%spc_gamg_",pcpre ? pcpre : "");CHKERRQ(ierr);
1178   ierr = PetscObjectSetOptionsPrefix((PetscObject)pc_gamg->random,prefix);CHKERRQ(ierr);
1179   ierr = PetscRandomSetFromOptions(pc_gamg->random);CHKERRQ(ierr);
1180   ierr = PetscOptionsTail();CHKERRQ(ierr);
1181   PetscFunctionReturn(0);
1182 }
1183 
1184 /* -------------------------------------------------------------------------- */
1185 /*MC
1186      PCGAMG - Geometric algebraic multigrid (AMG) preconditioner
1187 
1188    Options Database Keys:
1189 +   -pc_gamg_type <type> - one of agg, geo, or classical
1190 .   -pc_gamg_repartition  <true,default=false> - repartition the degrees of freedom accross the coarse grids as they are determined
1191 .   -pc_gamg_reuse_interpolation <true,default=false> - when rebuilding the algebraic multigrid preconditioner reuse the previously computed interpolations
1192 .   -pc_gamg_asm_use_agg <true,default=false> - use the aggregates from the coasening process to defined the subdomains on each level for the PCASM smoother
1193 .   -pc_gamg_process_eq_limit <limit, default=50> - GAMG will reduce the number of MPI processes used directly on the coarse grids so that there are around <limit>
1194                                         equations on each process that has degrees of freedom
1195 .   -pc_gamg_coarse_eq_limit <limit, default=50> - Set maximum number of equations on coarsest grid to aim for.
1196 -   -pc_gamg_threshold <thresh,default=0> - Before aggregating the graph GAMG will remove small values from the graph thus reducing the coupling in the graph and a different
1197 
1198    Options Database Keys for default Aggregation:
1199 +  -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use with smooth aggregation
1200 .  -pc_gamg_sym_graph <true,default=false> - symmetrize the graph before computing the aggregation
1201 -  -pc_gamg_square_graph <n,default=1> - number of levels to square the graph before aggregating it
1202 
1203    Multigrid options(inherited):
1204 +  -pc_mg_cycles <v>: v or w (PCMGSetCycleType())
1205 .  -pc_mg_smoothup <1>: Number of post-smoothing steps (PCMGSetNumberSmoothUp)
1206 .  -pc_mg_smoothdown <1>: Number of pre-smoothing steps (PCMGSetNumberSmoothDown)
1207 .  -pc_mg_type <multiplicative>: (one of) additive multiplicative full kascade
1208 -  -pc_mg_levels <levels> - Number of levels of multigrid to use.
1209 
1210 
1211   Notes: In order to obtain good performance for PCGAMG for vector valued problems you must
1212 $       Call MatSetBlockSize() to indicate the number of degrees of freedom per grid point
1213 $       Call MatSetNearNullSpace() (or PCSetCoordinates() if solving the equations of elasticity) to indicate the near null space of the operator
1214 $       See the Users Manual Chapter 4 for more details
1215 
1216   Level: intermediate
1217 
1218   Concepts: algebraic multigrid
1219 
1220 .seealso:  PCCreate(), PCSetType(), MatSetBlockSize(), PCMGType, PCSetCoordinates(), MatSetNearNullSpace(), PCGAMGSetType(), PCGAMGAGG, PCGAMGGEO, PCGAMGCLASSICAL, PCGAMGSetProcEqLim(),
1221            PCGAMGSetCoarseEqLim(), PCGAMGSetRepartition(), PCGAMGRegister(), PCGAMGSetReuseInterpolation(), PCGAMGASMSetUseAggs(), PCGAMGSetNlevels(), PCGAMGSetThreshold(), PCGAMGGetType(), PCGAMGSetReuseInterpolation()
1222 M*/
1223 
1224 #undef __FUNCT__
1225 #define __FUNCT__ "PCCreate_GAMG"
1226 PETSC_EXTERN PetscErrorCode PCCreate_GAMG(PC pc)
1227 {
1228   PetscErrorCode ierr;
1229   PC_GAMG        *pc_gamg;
1230   PC_MG          *mg;
1231 
1232   PetscFunctionBegin;
1233    /* register AMG type */
1234   ierr = PCGAMGInitializePackage();CHKERRQ(ierr);
1235 
1236   /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */
1237   ierr = PCSetType(pc, PCMG);CHKERRQ(ierr);
1238   ierr = PetscObjectChangeTypeName((PetscObject)pc, PCGAMG);CHKERRQ(ierr);
1239 
1240   /* create a supporting struct and attach it to pc */
1241   ierr         = PetscNewLog(pc,&pc_gamg);CHKERRQ(ierr);
1242   mg           = (PC_MG*)pc->data;
1243   mg->galerkin = 2;             /* Use Galerkin, but it is computed externally from PCMG by GAMG code */
1244   mg->innerctx = pc_gamg;
1245 
1246   ierr = PetscNewLog(pc,&pc_gamg->ops);CHKERRQ(ierr);
1247 
1248   pc_gamg->setup_count = 0;
1249   /* these should be in subctx but repartitioning needs simple arrays */
1250   pc_gamg->data_sz = 0;
1251   pc_gamg->data    = 0;
1252 
1253   /* overwrite the pointers of PCMG by the functions of base class PCGAMG */
1254   pc->ops->setfromoptions = PCSetFromOptions_GAMG;
1255   pc->ops->setup          = PCSetUp_GAMG;
1256   pc->ops->reset          = PCReset_GAMG;
1257   pc->ops->destroy        = PCDestroy_GAMG;
1258   mg->view                = PCView_GAMG;
1259 
1260   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetProcEqLim_C",PCGAMGSetProcEqLim_GAMG);CHKERRQ(ierr);
1261   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetCoarseEqLim_C",PCGAMGSetCoarseEqLim_GAMG);CHKERRQ(ierr);
1262   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetRepartition_C",PCGAMGSetRepartition_GAMG);CHKERRQ(ierr);
1263   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetReuseInterpolation_C",PCGAMGSetReuseInterpolation_GAMG);CHKERRQ(ierr);
1264   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGASMSetUseAggs_C",PCGAMGASMSetUseAggs_GAMG);CHKERRQ(ierr);
1265   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetThreshold_C",PCGAMGSetThreshold_GAMG);CHKERRQ(ierr);
1266   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetType_C",PCGAMGSetType_GAMG);CHKERRQ(ierr);
1267   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGGetType_C",PCGAMGGetType_GAMG);CHKERRQ(ierr);
1268   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetNlevels_C",PCGAMGSetNlevels_GAMG);CHKERRQ(ierr);
1269   pc_gamg->repart           = PETSC_FALSE;
1270   pc_gamg->reuse_prol       = PETSC_FALSE;
1271   pc_gamg->use_aggs_in_asm  = PETSC_FALSE;
1272   pc_gamg->min_eq_proc      = 50;
1273   pc_gamg->coarse_eq_limit  = 50;
1274   pc_gamg->threshold        = 0.;
1275   pc_gamg->Nlevels          = GAMG_MAXLEVELS;
1276   pc_gamg->current_level    = 0; /* don't need to init really */
1277   pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG;
1278 
1279   ierr = PetscRandomCreate(PetscObjectComm((PetscObject)pc),&pc_gamg->random);CHKERRQ(ierr);
1280 
1281   /* PCSetUp_GAMG assumes that the type has been set, so set it to the default now */
1282   ierr = PCGAMGSetType(pc,PCGAMGAGG);CHKERRQ(ierr);
1283   PetscFunctionReturn(0);
1284 }
1285 
1286 #undef __FUNCT__
1287 #define __FUNCT__ "PCGAMGInitializePackage"
1288 /*@C
1289  PCGAMGInitializePackage - This function initializes everything in the PCGAMG package. It is called
1290     from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to PCCreate_GAMG()
1291     when using static libraries.
1292 
1293  Level: developer
1294 
1295  .keywords: PC, PCGAMG, initialize, package
1296  .seealso: PetscInitialize()
1297 @*/
1298 PetscErrorCode PCGAMGInitializePackage(void)
1299 {
1300   PetscErrorCode ierr;
1301 
1302   PetscFunctionBegin;
1303   if (PCGAMGPackageInitialized) PetscFunctionReturn(0);
1304   PCGAMGPackageInitialized = PETSC_TRUE;
1305   ierr = PetscFunctionListAdd(&GAMGList,PCGAMGGEO,PCCreateGAMG_GEO);CHKERRQ(ierr);
1306   ierr = PetscFunctionListAdd(&GAMGList,PCGAMGAGG,PCCreateGAMG_AGG);CHKERRQ(ierr);
1307   ierr = PetscFunctionListAdd(&GAMGList,PCGAMGCLASSICAL,PCCreateGAMG_Classical);CHKERRQ(ierr);
1308   ierr = PetscRegisterFinalize(PCGAMGFinalizePackage);CHKERRQ(ierr);
1309 
1310   /* general events */
1311   ierr = PetscLogEventRegister("PCGAMGGraph_AGG", 0, &PC_GAMGGraph_AGG);CHKERRQ(ierr);
1312   ierr = PetscLogEventRegister("PCGAMGGraph_GEO", PC_CLASSID, &PC_GAMGGraph_GEO);CHKERRQ(ierr);
1313   ierr = PetscLogEventRegister("PCGAMGCoarse_AGG", PC_CLASSID, &PC_GAMGCoarsen_AGG);CHKERRQ(ierr);
1314   ierr = PetscLogEventRegister("PCGAMGCoarse_GEO", PC_CLASSID, &PC_GAMGCoarsen_GEO);CHKERRQ(ierr);
1315   ierr = PetscLogEventRegister("PCGAMGProl_AGG", PC_CLASSID, &PC_GAMGProlongator_AGG);CHKERRQ(ierr);
1316   ierr = PetscLogEventRegister("PCGAMGProl_GEO", PC_CLASSID, &PC_GAMGProlongator_GEO);CHKERRQ(ierr);
1317   ierr = PetscLogEventRegister("PCGAMGPOpt_AGG", PC_CLASSID, &PC_GAMGOptProlongator_AGG);CHKERRQ(ierr);
1318 
1319 #if defined PETSC_GAMG_USE_LOG
1320   ierr = PetscLogEventRegister("GAMG: createProl", PC_CLASSID, &petsc_gamg_setup_events[SET1]);CHKERRQ(ierr);
1321   ierr = PetscLogEventRegister("  Graph", PC_CLASSID, &petsc_gamg_setup_events[GRAPH]);CHKERRQ(ierr);
1322   /* PetscLogEventRegister("    G.Mat", PC_CLASSID, &petsc_gamg_setup_events[GRAPH_MAT]); */
1323   /* PetscLogEventRegister("    G.Filter", PC_CLASSID, &petsc_gamg_setup_events[GRAPH_FILTER]); */
1324   /* PetscLogEventRegister("    G.Square", PC_CLASSID, &petsc_gamg_setup_events[GRAPH_SQR]); */
1325   ierr = PetscLogEventRegister("  MIS/Agg", PC_CLASSID, &petsc_gamg_setup_events[SET4]);CHKERRQ(ierr);
1326   ierr = PetscLogEventRegister("  geo: growSupp", PC_CLASSID, &petsc_gamg_setup_events[SET5]);CHKERRQ(ierr);
1327   ierr = PetscLogEventRegister("  geo: triangle", PC_CLASSID, &petsc_gamg_setup_events[SET6]);CHKERRQ(ierr);
1328   ierr = PetscLogEventRegister("    search-set", PC_CLASSID, &petsc_gamg_setup_events[FIND_V]);CHKERRQ(ierr);
1329   ierr = PetscLogEventRegister("  SA: col data", PC_CLASSID, &petsc_gamg_setup_events[SET7]);CHKERRQ(ierr);
1330   ierr = PetscLogEventRegister("  SA: frmProl0", PC_CLASSID, &petsc_gamg_setup_events[SET8]);CHKERRQ(ierr);
1331   ierr = PetscLogEventRegister("  SA: smooth", PC_CLASSID, &petsc_gamg_setup_events[SET9]);CHKERRQ(ierr);
1332   ierr = PetscLogEventRegister("GAMG: partLevel", PC_CLASSID, &petsc_gamg_setup_events[SET2]);CHKERRQ(ierr);
1333   ierr = PetscLogEventRegister("  repartition", PC_CLASSID, &petsc_gamg_setup_events[SET12]);CHKERRQ(ierr);
1334   ierr = PetscLogEventRegister("  Invert-Sort", PC_CLASSID, &petsc_gamg_setup_events[SET13]);CHKERRQ(ierr);
1335   ierr = PetscLogEventRegister("  Move A", PC_CLASSID, &petsc_gamg_setup_events[SET14]);CHKERRQ(ierr);
1336   ierr = PetscLogEventRegister("  Move P", PC_CLASSID, &petsc_gamg_setup_events[SET15]);CHKERRQ(ierr);
1337 
1338   /* PetscLogEventRegister(" PL move data", PC_CLASSID, &petsc_gamg_setup_events[SET13]); */
1339   /* PetscLogEventRegister("GAMG: fix", PC_CLASSID, &petsc_gamg_setup_events[SET10]); */
1340   /* PetscLogEventRegister("GAMG: set levels", PC_CLASSID, &petsc_gamg_setup_events[SET11]); */
1341   /* create timer stages */
1342 #if defined GAMG_STAGES
1343   {
1344     char     str[32];
1345     PetscInt lidx;
1346     sprintf(str,"MG Level %d (finest)",0);
1347     ierr = PetscLogStageRegister(str, &gamg_stages[0]);CHKERRQ(ierr);
1348     for (lidx=1; lidx<9; lidx++) {
1349       sprintf(str,"MG Level %d",lidx);
1350       ierr = PetscLogStageRegister(str, &gamg_stages[lidx]);CHKERRQ(ierr);
1351     }
1352   }
1353 #endif
1354 #endif
1355   PetscFunctionReturn(0);
1356 }
1357 
1358 #undef __FUNCT__
1359 #define __FUNCT__ "PCGAMGFinalizePackage"
1360 /*@C
1361  PCGAMGFinalizePackage - This function frees everything from the PCGAMG package. It is
1362     called from PetscFinalize() automatically.
1363 
1364  Level: developer
1365 
1366  .keywords: Petsc, destroy, package
1367  .seealso: PetscFinalize()
1368 @*/
1369 PetscErrorCode PCGAMGFinalizePackage(void)
1370 {
1371   PetscErrorCode ierr;
1372 
1373   PetscFunctionBegin;
1374   PCGAMGPackageInitialized = PETSC_FALSE;
1375   ierr = PetscFunctionListDestroy(&GAMGList);CHKERRQ(ierr);
1376   PetscFunctionReturn(0);
1377 }
1378 
1379 #undef __FUNCT__
1380 #define __FUNCT__ "PCGAMGRegister"
1381 /*@C
1382  PCGAMGRegister - Register a PCGAMG implementation.
1383 
1384  Input Parameters:
1385  + type - string that will be used as the name of the GAMG type.
1386  - create - function for creating the gamg context.
1387 
1388   Level: advanced
1389 
1390  .seealso: PCGAMGType, PCGAMG, PCGAMGSetType()
1391 @*/
1392 PetscErrorCode PCGAMGRegister(PCGAMGType type, PetscErrorCode (*create)(PC))
1393 {
1394   PetscErrorCode ierr;
1395 
1396   PetscFunctionBegin;
1397   ierr = PCGAMGInitializePackage();CHKERRQ(ierr);
1398   ierr = PetscFunctionListAdd(&GAMGList,type,create);CHKERRQ(ierr);
1399   PetscFunctionReturn(0);
1400 }
1401 
1402