xref: /petsc/src/ksp/pc/impls/gamg/gamg.c (revision 7a28f3e5b83784e012fb5372ca74592ce67c00ce)
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, Gmat, &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       *iss;
612 
613         sz   = nASMBlocksArr[level];
614         iss   = ASMLocalIDsArr[level];
615         ierr = PCSetType(subpc, PCASM);CHKERRQ(ierr);
616         ierr = PCASMSetOverlap(subpc, 0);CHKERRQ(ierr);
617         ierr = PCASMSetType(subpc,PC_ASM_BASIC);CHKERRQ(ierr);
618         if (!sz) {
619           IS       is;
620           ierr = ISCreateGeneral(PETSC_COMM_SELF, 0, NULL, PETSC_COPY_VALUES, &is);CHKERRQ(ierr);
621           ierr = PCASMSetLocalSubdomains(subpc, 1, NULL, &is);CHKERRQ(ierr);
622           ierr = ISDestroy(&is);CHKERRQ(ierr);
623         } else {
624           PetscInt kk;
625           ierr = PCASMSetLocalSubdomains(subpc, sz, NULL, iss);CHKERRQ(ierr);
626           for (kk=0; kk<sz; kk++) {
627             ierr = ISDestroy(&iss[kk]);CHKERRQ(ierr);
628           }
629           ierr = PetscFree(iss);CHKERRQ(ierr);
630         }
631         ASMLocalIDsArr[level] = NULL;
632         nASMBlocksArr[level]  = 0;
633       } else {
634         ierr = PCSetType(subpc, PCSOR);CHKERRQ(ierr);
635       }
636     }
637     {
638       /* coarse grid */
639       KSP smoother,*k2; PC subpc,pc2; PetscInt ii,first;
640       Mat Lmat = Aarr[(level=pc_gamg->Nlevels-1)]; lidx = 0;
641       ierr = PCMGGetSmoother(pc, lidx, &smoother);CHKERRQ(ierr);
642       ierr = KSPSetOperators(smoother, Lmat, Lmat);CHKERRQ(ierr);
643       ierr = KSPSetNormType(smoother, KSP_NORM_NONE);CHKERRQ(ierr);
644       ierr = KSPGetPC(smoother, &subpc);CHKERRQ(ierr);
645       ierr = PCSetType(subpc, PCBJACOBI);CHKERRQ(ierr);
646       ierr = PCSetUp(subpc);CHKERRQ(ierr);
647       ierr = PCBJacobiGetSubKSP(subpc,&ii,&first,&k2);CHKERRQ(ierr);
648       if (ii != 1) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_PLIB,"ii %D is not one",ii);
649       ierr = KSPGetPC(k2[0],&pc2);CHKERRQ(ierr);
650       ierr = PCSetType(pc2, PCLU);CHKERRQ(ierr);
651       ierr = PCFactorSetShiftType(pc2,MAT_SHIFT_INBLOCKS);CHKERRQ(ierr);
652       ierr = KSPSetTolerances(k2[0],PETSC_DEFAULT,PETSC_DEFAULT,PETSC_DEFAULT,1);CHKERRQ(ierr);
653       ierr = KSPSetType(k2[0], KSPPREONLY);CHKERRQ(ierr);
654       /* This flag gets reset by PCBJacobiGetSubKSP(), but our BJacobi really does the same algorithm everywhere (and in
655        * fact, all but one process will have zero dofs), so we reset the flag to avoid having PCView_BJacobi attempt to
656        * view every subdomain as though they were different. */
657       ((PC_BJacobi*)subpc->data)->same_local_solves = PETSC_TRUE;
658     }
659 
660     /* should be called in PCSetFromOptions_GAMG(), but cannot be called prior to PCMGSetLevels() */
661     ierr = PetscObjectOptionsBegin((PetscObject)pc);CHKERRQ(ierr);
662     ierr = PCSetFromOptions_MG(PetscOptionsObject,pc);CHKERRQ(ierr);
663     ierr = PetscOptionsEnd();CHKERRQ(ierr);
664     if (!mg->galerkin) SETERRQ(comm,PETSC_ERR_USER,"PCGAMG must use Galerkin for coarse operators.");
665     if (mg->galerkin == 1) mg->galerkin = 2;
666 
667     /* clean up */
668     for (level=1; level<pc_gamg->Nlevels; level++) {
669       ierr = MatDestroy(&Parr[level]);CHKERRQ(ierr);
670       ierr = MatDestroy(&Aarr[level]);CHKERRQ(ierr);
671     }
672     ierr = PCSetUp_MG(pc);CHKERRQ(ierr);
673   } else {
674     KSP smoother;
675     ierr = PetscInfo(pc,"One level solver used (system is seen as DD). Using default solver.\n");CHKERRQ(ierr);
676     ierr = PCMGGetSmoother(pc, 0, &smoother);CHKERRQ(ierr);
677     ierr = KSPSetOperators(smoother, Aarr[0], Aarr[0]);CHKERRQ(ierr);
678     ierr = KSPSetType(smoother, KSPPREONLY);CHKERRQ(ierr);
679     ierr = PCSetUp_MG(pc);CHKERRQ(ierr);
680   }
681   PetscFunctionReturn(0);
682 }
683 
684 /* ------------------------------------------------------------------------- */
685 /*
686  PCDestroy_GAMG - Destroys the private context for the GAMG preconditioner
687    that was created with PCCreate_GAMG().
688 
689    Input Parameter:
690 .  pc - the preconditioner context
691 
692    Application Interface Routine: PCDestroy()
693 */
694 #undef __FUNCT__
695 #define __FUNCT__ "PCDestroy_GAMG"
696 PetscErrorCode PCDestroy_GAMG(PC pc)
697 {
698   PetscErrorCode ierr;
699   PC_MG          *mg     = (PC_MG*)pc->data;
700   PC_GAMG        *pc_gamg= (PC_GAMG*)mg->innerctx;
701 
702   PetscFunctionBegin;
703   ierr = PCReset_GAMG(pc);CHKERRQ(ierr);
704   if (pc_gamg->ops->destroy) {
705     ierr = (*pc_gamg->ops->destroy)(pc);CHKERRQ(ierr);
706   }
707   ierr = PetscRandomDestroy(&pc_gamg->random);CHKERRQ(ierr);
708   ierr = PetscFree(pc_gamg->ops);CHKERRQ(ierr);
709   ierr = PetscFree(pc_gamg->gamg_type_name);CHKERRQ(ierr);
710   ierr = PetscFree(pc_gamg);CHKERRQ(ierr);
711   ierr = PCDestroy_MG(pc);CHKERRQ(ierr);
712   PetscFunctionReturn(0);
713 }
714 
715 #undef __FUNCT__
716 #define __FUNCT__ "PCGAMGSetProcEqLim"
717 /*@
718    PCGAMGSetProcEqLim - Set number of equations to aim for per process on the coarse grids via processor reduction.
719 
720    Logically Collective on PC
721 
722    Input Parameters:
723 +  pc - the preconditioner context
724 -  n - the number of equations
725 
726 
727    Options Database Key:
728 .  -pc_gamg_process_eq_limit <limit>
729 
730    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
731           that has degrees of freedom
732 
733    Level: intermediate
734 
735    Concepts: Unstructured multigrid preconditioner
736 
737 .seealso: PCGAMGSetCoarseEqLim()
738 @*/
739 PetscErrorCode  PCGAMGSetProcEqLim(PC pc, PetscInt n)
740 {
741   PetscErrorCode ierr;
742 
743   PetscFunctionBegin;
744   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
745   ierr = PetscTryMethod(pc,"PCGAMGSetProcEqLim_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr);
746   PetscFunctionReturn(0);
747 }
748 
749 #undef __FUNCT__
750 #define __FUNCT__ "PCGAMGSetProcEqLim_GAMG"
751 static PetscErrorCode PCGAMGSetProcEqLim_GAMG(PC pc, PetscInt n)
752 {
753   PC_MG   *mg      = (PC_MG*)pc->data;
754   PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
755 
756   PetscFunctionBegin;
757   if (n>0) pc_gamg->min_eq_proc = n;
758   PetscFunctionReturn(0);
759 }
760 
761 #undef __FUNCT__
762 #define __FUNCT__ "PCGAMGSetCoarseEqLim"
763 /*@
764    PCGAMGSetCoarseEqLim - Set maximum number of equations on coarsest grid.
765 
766  Collective on PC
767 
768    Input Parameters:
769 +  pc - the preconditioner context
770 -  n - maximum number of equations to aim for
771 
772    Options Database Key:
773 .  -pc_gamg_coarse_eq_limit <limit>
774 
775    Level: intermediate
776 
777    Concepts: Unstructured multigrid preconditioner
778 
779 .seealso: PCGAMGSetProcEqLim()
780 @*/
781 PetscErrorCode PCGAMGSetCoarseEqLim(PC pc, PetscInt n)
782 {
783   PetscErrorCode ierr;
784 
785   PetscFunctionBegin;
786   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
787   ierr = PetscTryMethod(pc,"PCGAMGSetCoarseEqLim_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr);
788   PetscFunctionReturn(0);
789 }
790 
791 #undef __FUNCT__
792 #define __FUNCT__ "PCGAMGSetCoarseEqLim_GAMG"
793 static PetscErrorCode PCGAMGSetCoarseEqLim_GAMG(PC pc, PetscInt n)
794 {
795   PC_MG   *mg      = (PC_MG*)pc->data;
796   PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
797 
798   PetscFunctionBegin;
799   if (n>0) pc_gamg->coarse_eq_limit = n;
800   PetscFunctionReturn(0);
801 }
802 
803 #undef __FUNCT__
804 #define __FUNCT__ "PCGAMGSetRepartition"
805 /*@
806    PCGAMGSetRepartition - Repartition the degrees of freedom across the processors on the coarser grids
807 
808    Collective on PC
809 
810    Input Parameters:
811 +  pc - the preconditioner context
812 -  n - PETSC_TRUE or PETSC_FALSE
813 
814    Options Database Key:
815 .  -pc_gamg_repartition <true,false>
816 
817    Notes: this will generally improve the loading balancing of the work on each level
818 
819    Level: intermediate
820 
821    Concepts: Unstructured multigrid preconditioner
822 
823 .seealso: ()
824 @*/
825 PetscErrorCode PCGAMGSetRepartition(PC pc, PetscBool n)
826 {
827   PetscErrorCode ierr;
828 
829   PetscFunctionBegin;
830   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
831   ierr = PetscTryMethod(pc,"PCGAMGSetRepartition_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr);
832   PetscFunctionReturn(0);
833 }
834 
835 #undef __FUNCT__
836 #define __FUNCT__ "PCGAMGSetRepartition_GAMG"
837 static PetscErrorCode PCGAMGSetRepartition_GAMG(PC pc, PetscBool n)
838 {
839   PC_MG   *mg      = (PC_MG*)pc->data;
840   PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
841 
842   PetscFunctionBegin;
843   pc_gamg->repart = n;
844   PetscFunctionReturn(0);
845 }
846 
847 #undef __FUNCT__
848 #define __FUNCT__ "PCGAMGSetReuseInterpolation"
849 /*@
850    PCGAMGSetReuseInterpolation - Reuse prolongation when rebuilding algebraic multigrid preconditioner
851 
852    Collective on PC
853 
854    Input Parameters:
855 +  pc - the preconditioner context
856 -  n - PETSC_TRUE or PETSC_FALSE
857 
858    Options Database Key:
859 .  -pc_gamg_reuse_interpolation <true,false>
860 
861    Level: intermediate
862 
863    Notes: this may negatively affect the convergence rate of the method on new matrices if the matrix entries change a great deal, but allows
864           rebuilding the preconditioner quicker.
865 
866    Concepts: Unstructured multigrid preconditioner
867 
868 .seealso: ()
869 @*/
870 PetscErrorCode PCGAMGSetReuseInterpolation(PC pc, PetscBool n)
871 {
872   PetscErrorCode ierr;
873 
874   PetscFunctionBegin;
875   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
876   ierr = PetscTryMethod(pc,"PCGAMGSetReuseInterpolation_C",(PC,PetscBool),(pc,n));CHKERRQ(ierr);
877   PetscFunctionReturn(0);
878 }
879 
880 #undef __FUNCT__
881 #define __FUNCT__ "PCGAMGSetReuseInterpolation_GAMG"
882 static PetscErrorCode PCGAMGSetReuseInterpolation_GAMG(PC pc, PetscBool n)
883 {
884   PC_MG   *mg      = (PC_MG*)pc->data;
885   PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
886 
887   PetscFunctionBegin;
888   pc_gamg->reuse_prol = n;
889   PetscFunctionReturn(0);
890 }
891 
892 #undef __FUNCT__
893 #define __FUNCT__ "PCGAMGASMSetUseAggs"
894 /*@
895    PCGAMGASMSetUseAggs - Have the PCGAMG smoother on each level use the aggregates defined by the coarsening process as the subdomains for the additive Schwarz preconditioner.
896 
897    Collective on PC
898 
899    Input Parameters:
900 +  pc - the preconditioner context
901 -  flg - PETSC_TRUE to use aggregates, PETSC_FALSE to not
902 
903    Options Database Key:
904 .  -pc_gamg_asm_use_agg
905 
906    Level: intermediate
907 
908    Concepts: Unstructured multigrid preconditioner
909 
910 .seealso: ()
911 @*/
912 PetscErrorCode PCGAMGASMSetUseAggs(PC pc, PetscBool flg)
913 {
914   PetscErrorCode ierr;
915 
916   PetscFunctionBegin;
917   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
918   ierr = PetscTryMethod(pc,"PCGAMGASMSetUseAggs_C",(PC,PetscBool),(pc,flg));CHKERRQ(ierr);
919   PetscFunctionReturn(0);
920 }
921 
922 #undef __FUNCT__
923 #define __FUNCT__ "PCGAMGASMSetUseAggs_GAMG"
924 static PetscErrorCode PCGAMGASMSetUseAggs_GAMG(PC pc, PetscBool flg)
925 {
926   PC_MG   *mg      = (PC_MG*)pc->data;
927   PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
928 
929   PetscFunctionBegin;
930   pc_gamg->use_aggs_in_asm = flg;
931   PetscFunctionReturn(0);
932 }
933 
934 #undef __FUNCT__
935 #define __FUNCT__ "PCGAMGSetNlevels"
936 /*@
937    PCGAMGSetNlevels -  Sets the maximum number of levels PCGAMG will use
938 
939    Not collective on PC
940 
941    Input Parameters:
942 +  pc - the preconditioner
943 -  n - the maximum number of levels to use
944 
945    Options Database Key:
946 .  -pc_mg_levels
947 
948    Level: intermediate
949 
950    Concepts: Unstructured multigrid preconditioner
951 
952 .seealso: ()
953 @*/
954 PetscErrorCode PCGAMGSetNlevels(PC pc, PetscInt n)
955 {
956   PetscErrorCode ierr;
957 
958   PetscFunctionBegin;
959   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
960   ierr = PetscTryMethod(pc,"PCGAMGSetNlevels_C",(PC,PetscInt),(pc,n));CHKERRQ(ierr);
961   PetscFunctionReturn(0);
962 }
963 
964 #undef __FUNCT__
965 #define __FUNCT__ "PCGAMGSetNlevels_GAMG"
966 static PetscErrorCode PCGAMGSetNlevels_GAMG(PC pc, PetscInt n)
967 {
968   PC_MG   *mg      = (PC_MG*)pc->data;
969   PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
970 
971   PetscFunctionBegin;
972   pc_gamg->Nlevels = n;
973   PetscFunctionReturn(0);
974 }
975 
976 #undef __FUNCT__
977 #define __FUNCT__ "PCGAMGSetThreshold"
978 /*@
979    PCGAMGSetThreshold - Relative threshold to use for dropping edges in aggregation graph
980 
981    Not collective on PC
982 
983    Input Parameters:
984 +  pc - the preconditioner context
985 -  threshold - the threshold value, 0.0 means keep all nonzero entries in the graph; negative means keep even zero entries in the graph
986 
987    Options Database Key:
988 .  -pc_gamg_threshold <threshold>
989 
990    Notes: Before aggregating the graph GAMG will remove small values from the graph thus reducing the coupling in the graph and a different
991     (perhaps better) coarser set of points.
992 
993    Level: intermediate
994 
995    Concepts: Unstructured multigrid preconditioner
996 
997 .seealso: ()
998 @*/
999 PetscErrorCode PCGAMGSetThreshold(PC pc, PetscReal n)
1000 {
1001   PetscErrorCode ierr;
1002 
1003   PetscFunctionBegin;
1004   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1005   ierr = PetscTryMethod(pc,"PCGAMGSetThreshold_C",(PC,PetscReal),(pc,n));CHKERRQ(ierr);
1006   PetscFunctionReturn(0);
1007 }
1008 
1009 #undef __FUNCT__
1010 #define __FUNCT__ "PCGAMGSetThreshold_GAMG"
1011 static PetscErrorCode PCGAMGSetThreshold_GAMG(PC pc, PetscReal n)
1012 {
1013   PC_MG   *mg      = (PC_MG*)pc->data;
1014   PC_GAMG *pc_gamg = (PC_GAMG*)mg->innerctx;
1015 
1016   PetscFunctionBegin;
1017   pc_gamg->threshold = n;
1018   PetscFunctionReturn(0);
1019 }
1020 
1021 #undef __FUNCT__
1022 #define __FUNCT__ "PCGAMGSetType"
1023 /*@
1024    PCGAMGSetType - Set solution method
1025 
1026    Collective on PC
1027 
1028    Input Parameters:
1029 +  pc - the preconditioner context
1030 -  type - PCGAMGAGG, PCGAMGGEO, or PCGAMGCLASSICAL
1031 
1032    Options Database Key:
1033 .  -pc_gamg_type <agg,geo,classical> - type of algebraic multigrid to apply
1034 
1035    Level: intermediate
1036 
1037    Concepts: Unstructured multigrid preconditioner
1038 
1039 .seealso: PCGAMGGetType(), PCGAMG, PCGAMGType
1040 @*/
1041 PetscErrorCode PCGAMGSetType(PC pc, PCGAMGType type)
1042 {
1043   PetscErrorCode ierr;
1044 
1045   PetscFunctionBegin;
1046   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1047   ierr = PetscTryMethod(pc,"PCGAMGSetType_C",(PC,PCGAMGType),(pc,type));CHKERRQ(ierr);
1048   PetscFunctionReturn(0);
1049 }
1050 
1051 #undef __FUNCT__
1052 #define __FUNCT__ "PCGAMGGetType"
1053 /*@
1054    PCGAMGGetType - Get solution method
1055 
1056    Collective on PC
1057 
1058    Input Parameter:
1059 .  pc - the preconditioner context
1060 
1061    Output Parameter:
1062 .  type - the type of algorithm used
1063 
1064    Level: intermediate
1065 
1066    Concepts: Unstructured multigrid preconditioner
1067 
1068 .seealso: PCGAMGSetType(), PCGAMGType
1069 @*/
1070 PetscErrorCode PCGAMGGetType(PC pc, PCGAMGType *type)
1071 {
1072   PetscErrorCode ierr;
1073 
1074   PetscFunctionBegin;
1075   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1076   ierr = PetscUseMethod(pc,"PCGAMGGetType_C",(PC,PCGAMGType*),(pc,type));CHKERRQ(ierr);
1077   PetscFunctionReturn(0);
1078 }
1079 
1080 #undef __FUNCT__
1081 #define __FUNCT__ "PCGAMGGetType_GAMG"
1082 static PetscErrorCode PCGAMGGetType_GAMG(PC pc, PCGAMGType *type)
1083 {
1084   PC_MG          *mg      = (PC_MG*)pc->data;
1085   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
1086 
1087   PetscFunctionBegin;
1088   *type = pc_gamg->type;
1089   PetscFunctionReturn(0);
1090 }
1091 
1092 #undef __FUNCT__
1093 #define __FUNCT__ "PCGAMGSetType_GAMG"
1094 static PetscErrorCode PCGAMGSetType_GAMG(PC pc, PCGAMGType type)
1095 {
1096   PetscErrorCode ierr,(*r)(PC);
1097   PC_MG          *mg      = (PC_MG*)pc->data;
1098   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
1099 
1100   PetscFunctionBegin;
1101   pc_gamg->type = type;
1102   ierr = PetscFunctionListFind(GAMGList,type,&r);CHKERRQ(ierr);
1103   if (!r) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown GAMG type %s given",type);
1104   if (pc_gamg->ops->destroy) {
1105     ierr = (*pc_gamg->ops->destroy)(pc);CHKERRQ(ierr);
1106     ierr = PetscMemzero(pc_gamg->ops,sizeof(struct _PCGAMGOps));CHKERRQ(ierr);
1107     pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG;
1108     /* cleaning up common data in pc_gamg - this should disapear someday */
1109     pc_gamg->data_cell_cols = 0;
1110     pc_gamg->data_cell_rows = 0;
1111     pc_gamg->orig_data_cell_cols = 0;
1112     pc_gamg->orig_data_cell_rows = 0;
1113     ierr = PetscFree(pc_gamg->data);CHKERRQ(ierr);
1114     pc_gamg->data_sz = 0;
1115   }
1116   ierr = PetscFree(pc_gamg->gamg_type_name);CHKERRQ(ierr);
1117   ierr = PetscStrallocpy(type,&pc_gamg->gamg_type_name);CHKERRQ(ierr);
1118   ierr = (*r)(pc);CHKERRQ(ierr);
1119   PetscFunctionReturn(0);
1120 }
1121 
1122 #undef __FUNCT__
1123 #define __FUNCT__ "PCView_GAMG"
1124 static PetscErrorCode PCView_GAMG(PC pc,PetscViewer viewer)
1125 {
1126   PetscErrorCode ierr;
1127   PC_MG          *mg      = (PC_MG*)pc->data;
1128   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
1129 
1130   PetscFunctionBegin;
1131   ierr = PetscViewerASCIIPrintf(viewer,"    GAMG specific options\n");CHKERRQ(ierr);
1132   ierr = PetscViewerASCIIPrintf(viewer,"      Threshold for dropping small values from graph %g\n",(double)pc_gamg->threshold);CHKERRQ(ierr);
1133   if (pc_gamg->use_aggs_in_asm) {
1134     ierr = PetscViewerASCIIPrintf(viewer,"      Using aggregates from coarsening process to define subdomains for PCASM\n");CHKERRQ(ierr);
1135   }
1136   if (pc_gamg->ops->view) {
1137     ierr = (*pc_gamg->ops->view)(pc,viewer);CHKERRQ(ierr);
1138   }
1139   PetscFunctionReturn(0);
1140 }
1141 
1142 #undef __FUNCT__
1143 #define __FUNCT__ "PCSetFromOptions_GAMG"
1144 PetscErrorCode PCSetFromOptions_GAMG(PetscOptionItems *PetscOptionsObject,PC pc)
1145 {
1146   PetscErrorCode ierr;
1147   PC_MG          *mg      = (PC_MG*)pc->data;
1148   PC_GAMG        *pc_gamg = (PC_GAMG*)mg->innerctx;
1149   PetscBool      flag;
1150   MPI_Comm       comm;
1151   char           prefix[256];
1152   const char     *pcpre;
1153 
1154   PetscFunctionBegin;
1155   ierr = PetscObjectGetComm((PetscObject)pc,&comm);CHKERRQ(ierr);
1156   ierr = PetscOptionsHead(PetscOptionsObject,"GAMG options");CHKERRQ(ierr);
1157   {
1158     char tname[256];
1159     ierr = PetscOptionsFList("-pc_gamg_type","Type of AMG method","PCGAMGSetType",GAMGList, pc_gamg->gamg_type_name, tname, sizeof(tname), &flag);CHKERRQ(ierr);
1160     if (flag) {
1161       ierr = PCGAMGSetType(pc,tname);CHKERRQ(ierr);
1162     }
1163     ierr = PetscOptionsBool("-pc_gamg_repartition","Repartion coarse grids","PCGAMGSetRepartition",pc_gamg->repart,&pc_gamg->repart,NULL);CHKERRQ(ierr);
1164     ierr = PetscOptionsBool("-pc_gamg_reuse_interpolation","Reuse prolongation operator","PCGAMGReuseInterpolation",pc_gamg->reuse_prol,&pc_gamg->reuse_prol,NULL);CHKERRQ(ierr);
1165     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);
1166     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);
1167     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);
1168     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);
1169     ierr = PetscOptionsInt("-pc_mg_levels","Set number of MG levels","PCGAMGSetNlevels",pc_gamg->Nlevels,&pc_gamg->Nlevels,NULL);CHKERRQ(ierr);
1170 
1171     /* set options for subtype */
1172     if (pc_gamg->ops->setfromoptions) {ierr = (*pc_gamg->ops->setfromoptions)(PetscOptionsObject,pc);CHKERRQ(ierr);}
1173   }
1174   ierr = PCGetOptionsPrefix(pc, &pcpre);CHKERRQ(ierr);
1175   ierr = PetscSNPrintf(prefix,sizeof(prefix),"%spc_gamg_",pcpre ? pcpre : "");CHKERRQ(ierr);
1176   ierr = PetscObjectSetOptionsPrefix((PetscObject)pc_gamg->random,prefix);CHKERRQ(ierr);
1177   ierr = PetscRandomSetFromOptions(pc_gamg->random);CHKERRQ(ierr);
1178   ierr = PetscOptionsTail();CHKERRQ(ierr);
1179   PetscFunctionReturn(0);
1180 }
1181 
1182 /* -------------------------------------------------------------------------- */
1183 /*MC
1184      PCGAMG - Geometric algebraic multigrid (AMG) preconditioner
1185 
1186    Options Database Keys:
1187 +   -pc_gamg_type <type> - one of agg, geo, or classical
1188 .   -pc_gamg_repartition  <true,default=false> - repartition the degrees of freedom accross the coarse grids as they are determined
1189 .   -pc_gamg_reuse_interpolation <true,default=false> - when rebuilding the algebraic multigrid preconditioner reuse the previously computed interpolations
1190 .   -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
1191 .   -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>
1192                                         equations on each process that has degrees of freedom
1193 .   -pc_gamg_coarse_eq_limit <limit, default=50> - Set maximum number of equations on coarsest grid to aim for.
1194 -   -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
1195 
1196    Options Database Keys for default Aggregation:
1197 +  -pc_gamg_agg_nsmooths <nsmooth, default=1> - number of smoothing steps to use with smooth aggregation
1198 .  -pc_gamg_sym_graph <true,default=false> - symmetrize the graph before computing the aggregation
1199 -  -pc_gamg_square_graph <n,default=1> - number of levels to square the graph before aggregating it
1200 
1201    Multigrid options(inherited):
1202 +  -pc_mg_cycles <v>: v or w (PCMGSetCycleType())
1203 .  -pc_mg_smoothup <1>: Number of post-smoothing steps (PCMGSetNumberSmoothUp)
1204 .  -pc_mg_smoothdown <1>: Number of pre-smoothing steps (PCMGSetNumberSmoothDown)
1205 .  -pc_mg_type <multiplicative>: (one of) additive multiplicative full kascade
1206 -  -pc_mg_levels <levels> - Number of levels of multigrid to use.
1207 
1208 
1209   Notes: In order to obtain good performance for PCGAMG for vector valued problems you must
1210 $       Call MatSetBlockSize() to indicate the number of degrees of freedom per grid point
1211 $       Call MatSetNearNullSpace() (or PCSetCoordinates() if solving the equations of elasticity) to indicate the near null space of the operator
1212 $       See the Users Manual Chapter 4 for more details
1213 
1214   Level: intermediate
1215 
1216   Concepts: algebraic multigrid
1217 
1218 .seealso:  PCCreate(), PCSetType(), MatSetBlockSize(), PCMGType, PCSetCoordinates(), MatSetNearNullSpace(), PCGAMGSetType(), PCGAMGAGG, PCGAMGGEO, PCGAMGCLASSICAL, PCGAMGSetProcEqLim(),
1219            PCGAMGSetCoarseEqLim(), PCGAMGSetRepartition(), PCGAMGRegister(), PCGAMGSetReuseInterpolation(), PCGAMGASMSetUseAggs(), PCGAMGSetNlevels(), PCGAMGSetThreshold(), PCGAMGGetType(), PCGAMGSetReuseInterpolation()
1220 M*/
1221 
1222 #undef __FUNCT__
1223 #define __FUNCT__ "PCCreate_GAMG"
1224 PETSC_EXTERN PetscErrorCode PCCreate_GAMG(PC pc)
1225 {
1226   PetscErrorCode ierr;
1227   PC_GAMG        *pc_gamg;
1228   PC_MG          *mg;
1229 
1230   PetscFunctionBegin;
1231    /* register AMG type */
1232   ierr = PCGAMGInitializePackage();CHKERRQ(ierr);
1233 
1234   /* PCGAMG is an inherited class of PCMG. Initialize pc as PCMG */
1235   ierr = PCSetType(pc, PCMG);CHKERRQ(ierr);
1236   ierr = PetscObjectChangeTypeName((PetscObject)pc, PCGAMG);CHKERRQ(ierr);
1237 
1238   /* create a supporting struct and attach it to pc */
1239   ierr         = PetscNewLog(pc,&pc_gamg);CHKERRQ(ierr);
1240   mg           = (PC_MG*)pc->data;
1241   mg->galerkin = 2;             /* Use Galerkin, but it is computed externally from PCMG by GAMG code */
1242   mg->innerctx = pc_gamg;
1243 
1244   ierr = PetscNewLog(pc,&pc_gamg->ops);CHKERRQ(ierr);
1245 
1246   pc_gamg->setup_count = 0;
1247   /* these should be in subctx but repartitioning needs simple arrays */
1248   pc_gamg->data_sz = 0;
1249   pc_gamg->data    = 0;
1250 
1251   /* overwrite the pointers of PCMG by the functions of base class PCGAMG */
1252   pc->ops->setfromoptions = PCSetFromOptions_GAMG;
1253   pc->ops->setup          = PCSetUp_GAMG;
1254   pc->ops->reset          = PCReset_GAMG;
1255   pc->ops->destroy        = PCDestroy_GAMG;
1256   mg->view                = PCView_GAMG;
1257 
1258   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetProcEqLim_C",PCGAMGSetProcEqLim_GAMG);CHKERRQ(ierr);
1259   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetCoarseEqLim_C",PCGAMGSetCoarseEqLim_GAMG);CHKERRQ(ierr);
1260   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetRepartition_C",PCGAMGSetRepartition_GAMG);CHKERRQ(ierr);
1261   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetReuseInterpolation_C",PCGAMGSetReuseInterpolation_GAMG);CHKERRQ(ierr);
1262   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGASMSetUseAggs_C",PCGAMGASMSetUseAggs_GAMG);CHKERRQ(ierr);
1263   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetThreshold_C",PCGAMGSetThreshold_GAMG);CHKERRQ(ierr);
1264   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetType_C",PCGAMGSetType_GAMG);CHKERRQ(ierr);
1265   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGGetType_C",PCGAMGGetType_GAMG);CHKERRQ(ierr);
1266   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGAMGSetNlevels_C",PCGAMGSetNlevels_GAMG);CHKERRQ(ierr);
1267   pc_gamg->repart           = PETSC_FALSE;
1268   pc_gamg->reuse_prol       = PETSC_FALSE;
1269   pc_gamg->use_aggs_in_asm  = PETSC_FALSE;
1270   pc_gamg->min_eq_proc      = 50;
1271   pc_gamg->coarse_eq_limit  = 50;
1272   pc_gamg->threshold        = 0.;
1273   pc_gamg->Nlevels          = GAMG_MAXLEVELS;
1274   pc_gamg->current_level    = 0; /* don't need to init really */
1275   pc_gamg->ops->createlevel = PCGAMGCreateLevel_GAMG;
1276 
1277   ierr = PetscRandomCreate(PetscObjectComm((PetscObject)pc),&pc_gamg->random);CHKERRQ(ierr);
1278 
1279   /* PCSetUp_GAMG assumes that the type has been set, so set it to the default now */
1280   ierr = PCGAMGSetType(pc,PCGAMGAGG);CHKERRQ(ierr);
1281   PetscFunctionReturn(0);
1282 }
1283 
1284 #undef __FUNCT__
1285 #define __FUNCT__ "PCGAMGInitializePackage"
1286 /*@C
1287  PCGAMGInitializePackage - This function initializes everything in the PCGAMG package. It is called
1288     from PetscDLLibraryRegister() when using dynamic libraries, and on the first call to PCCreate_GAMG()
1289     when using static libraries.
1290 
1291  Level: developer
1292 
1293  .keywords: PC, PCGAMG, initialize, package
1294  .seealso: PetscInitialize()
1295 @*/
1296 PetscErrorCode PCGAMGInitializePackage(void)
1297 {
1298   PetscErrorCode ierr;
1299 
1300   PetscFunctionBegin;
1301   if (PCGAMGPackageInitialized) PetscFunctionReturn(0);
1302   PCGAMGPackageInitialized = PETSC_TRUE;
1303   ierr = PetscFunctionListAdd(&GAMGList,PCGAMGGEO,PCCreateGAMG_GEO);CHKERRQ(ierr);
1304   ierr = PetscFunctionListAdd(&GAMGList,PCGAMGAGG,PCCreateGAMG_AGG);CHKERRQ(ierr);
1305   ierr = PetscFunctionListAdd(&GAMGList,PCGAMGCLASSICAL,PCCreateGAMG_Classical);CHKERRQ(ierr);
1306   ierr = PetscRegisterFinalize(PCGAMGFinalizePackage);CHKERRQ(ierr);
1307 
1308   /* general events */
1309   ierr = PetscLogEventRegister("PCGAMGGraph_AGG", 0, &PC_GAMGGraph_AGG);CHKERRQ(ierr);
1310   ierr = PetscLogEventRegister("PCGAMGGraph_GEO", PC_CLASSID, &PC_GAMGGraph_GEO);CHKERRQ(ierr);
1311   ierr = PetscLogEventRegister("PCGAMGCoarse_AGG", PC_CLASSID, &PC_GAMGCoarsen_AGG);CHKERRQ(ierr);
1312   ierr = PetscLogEventRegister("PCGAMGCoarse_GEO", PC_CLASSID, &PC_GAMGCoarsen_GEO);CHKERRQ(ierr);
1313   ierr = PetscLogEventRegister("PCGAMGProl_AGG", PC_CLASSID, &PC_GAMGProlongator_AGG);CHKERRQ(ierr);
1314   ierr = PetscLogEventRegister("PCGAMGProl_GEO", PC_CLASSID, &PC_GAMGProlongator_GEO);CHKERRQ(ierr);
1315   ierr = PetscLogEventRegister("PCGAMGPOpt_AGG", PC_CLASSID, &PC_GAMGOptProlongator_AGG);CHKERRQ(ierr);
1316 
1317 #if defined PETSC_GAMG_USE_LOG
1318   ierr = PetscLogEventRegister("GAMG: createProl", PC_CLASSID, &petsc_gamg_setup_events[SET1]);CHKERRQ(ierr);
1319   ierr = PetscLogEventRegister("  Graph", PC_CLASSID, &petsc_gamg_setup_events[GRAPH]);CHKERRQ(ierr);
1320   /* PetscLogEventRegister("    G.Mat", PC_CLASSID, &petsc_gamg_setup_events[GRAPH_MAT]); */
1321   /* PetscLogEventRegister("    G.Filter", PC_CLASSID, &petsc_gamg_setup_events[GRAPH_FILTER]); */
1322   /* PetscLogEventRegister("    G.Square", PC_CLASSID, &petsc_gamg_setup_events[GRAPH_SQR]); */
1323   ierr = PetscLogEventRegister("  MIS/Agg", PC_CLASSID, &petsc_gamg_setup_events[SET4]);CHKERRQ(ierr);
1324   ierr = PetscLogEventRegister("  geo: growSupp", PC_CLASSID, &petsc_gamg_setup_events[SET5]);CHKERRQ(ierr);
1325   ierr = PetscLogEventRegister("  geo: triangle", PC_CLASSID, &petsc_gamg_setup_events[SET6]);CHKERRQ(ierr);
1326   ierr = PetscLogEventRegister("    search-set", PC_CLASSID, &petsc_gamg_setup_events[FIND_V]);CHKERRQ(ierr);
1327   ierr = PetscLogEventRegister("  SA: col data", PC_CLASSID, &petsc_gamg_setup_events[SET7]);CHKERRQ(ierr);
1328   ierr = PetscLogEventRegister("  SA: frmProl0", PC_CLASSID, &petsc_gamg_setup_events[SET8]);CHKERRQ(ierr);
1329   ierr = PetscLogEventRegister("  SA: smooth", PC_CLASSID, &petsc_gamg_setup_events[SET9]);CHKERRQ(ierr);
1330   ierr = PetscLogEventRegister("GAMG: partLevel", PC_CLASSID, &petsc_gamg_setup_events[SET2]);CHKERRQ(ierr);
1331   ierr = PetscLogEventRegister("  repartition", PC_CLASSID, &petsc_gamg_setup_events[SET12]);CHKERRQ(ierr);
1332   ierr = PetscLogEventRegister("  Invert-Sort", PC_CLASSID, &petsc_gamg_setup_events[SET13]);CHKERRQ(ierr);
1333   ierr = PetscLogEventRegister("  Move A", PC_CLASSID, &petsc_gamg_setup_events[SET14]);CHKERRQ(ierr);
1334   ierr = PetscLogEventRegister("  Move P", PC_CLASSID, &petsc_gamg_setup_events[SET15]);CHKERRQ(ierr);
1335 
1336   /* PetscLogEventRegister(" PL move data", PC_CLASSID, &petsc_gamg_setup_events[SET13]); */
1337   /* PetscLogEventRegister("GAMG: fix", PC_CLASSID, &petsc_gamg_setup_events[SET10]); */
1338   /* PetscLogEventRegister("GAMG: set levels", PC_CLASSID, &petsc_gamg_setup_events[SET11]); */
1339   /* create timer stages */
1340 #if defined GAMG_STAGES
1341   {
1342     char     str[32];
1343     PetscInt lidx;
1344     sprintf(str,"MG Level %d (finest)",0);
1345     ierr = PetscLogStageRegister(str, &gamg_stages[0]);CHKERRQ(ierr);
1346     for (lidx=1; lidx<9; lidx++) {
1347       sprintf(str,"MG Level %d",lidx);
1348       ierr = PetscLogStageRegister(str, &gamg_stages[lidx]);CHKERRQ(ierr);
1349     }
1350   }
1351 #endif
1352 #endif
1353   PetscFunctionReturn(0);
1354 }
1355 
1356 #undef __FUNCT__
1357 #define __FUNCT__ "PCGAMGFinalizePackage"
1358 /*@C
1359  PCGAMGFinalizePackage - This function frees everything from the PCGAMG package. It is
1360     called from PetscFinalize() automatically.
1361 
1362  Level: developer
1363 
1364  .keywords: Petsc, destroy, package
1365  .seealso: PetscFinalize()
1366 @*/
1367 PetscErrorCode PCGAMGFinalizePackage(void)
1368 {
1369   PetscErrorCode ierr;
1370 
1371   PetscFunctionBegin;
1372   PCGAMGPackageInitialized = PETSC_FALSE;
1373   ierr = PetscFunctionListDestroy(&GAMGList);CHKERRQ(ierr);
1374   PetscFunctionReturn(0);
1375 }
1376 
1377 #undef __FUNCT__
1378 #define __FUNCT__ "PCGAMGRegister"
1379 /*@C
1380  PCGAMGRegister - Register a PCGAMG implementation.
1381 
1382  Input Parameters:
1383  + type - string that will be used as the name of the GAMG type.
1384  - create - function for creating the gamg context.
1385 
1386   Level: advanced
1387 
1388  .seealso: PCGAMGType, PCGAMG, PCGAMGSetType()
1389 @*/
1390 PetscErrorCode PCGAMGRegister(PCGAMGType type, PetscErrorCode (*create)(PC))
1391 {
1392   PetscErrorCode ierr;
1393 
1394   PetscFunctionBegin;
1395   ierr = PCGAMGInitializePackage();CHKERRQ(ierr);
1396   ierr = PetscFunctionListAdd(&GAMGList,type,create);CHKERRQ(ierr);
1397   PetscFunctionReturn(0);
1398 }
1399 
1400