xref: /petsc/src/ksp/pc/impls/hypre/hypre.c (revision 21df291b450e4ebf27ac0fc2b4dbb8de3f025202)
1 
2 /*
3    Provides an interface to the LLNL package hypre
4 */
5 
6 /* Must use hypre 2.0.0 or more recent. */
7 
8 #include <petsc/private/pcimpl.h>          /*I "petscpc.h" I*/
9 #include <../src/dm/impls/da/hypre/mhyp.h>
10 #include <_hypre_parcsr_ls.h>
11 
12 static PetscBool cite = PETSC_FALSE;
13 static const char hypreCitation[] = "@manual{hypre-web-page,\n  title  = {{\\sl hypre}: High Performance Preconditioners},\n  organization = {Lawrence Livermore National Laboratory},\n  note  = {\\url{http://www.llnl.gov/CASC/hypre/}}\n}\n";
14 
15 /*
16    Private context (data structure) for the  preconditioner.
17 */
18 typedef struct {
19   HYPRE_Solver   hsolver;
20   HYPRE_IJMatrix ij;
21   HYPRE_IJVector b,x;
22 
23   HYPRE_Int (*destroy)(HYPRE_Solver);
24   HYPRE_Int (*solve)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
25   HYPRE_Int (*setup)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
26   HYPRE_Int (*setdgrad)(HYPRE_Solver,HYPRE_ParCSRMatrix);
27   HYPRE_Int (*setdcurl)(HYPRE_Solver,HYPRE_ParCSRMatrix);
28   HYPRE_Int (*setcoord)(HYPRE_Solver,HYPRE_ParVector,HYPRE_ParVector,HYPRE_ParVector);
29   HYPRE_Int (*setdim)(HYPRE_Solver,HYPRE_Int);
30 
31   MPI_Comm comm_hypre;
32   char     *hypre_type;
33 
34   /* options for Pilut and BoomerAMG*/
35   PetscInt maxiter;
36   double   tol;
37 
38   /* options for Pilut */
39   PetscInt factorrowsize;
40 
41   /* options for ParaSails */
42   PetscInt nlevels;
43   double   threshhold;
44   double   filter;
45   PetscInt sym;
46   double   loadbal;
47   PetscInt logging;
48   PetscInt ruse;
49   PetscInt symt;
50 
51   /* options for BoomerAMG */
52   PetscBool printstatistics;
53 
54   /* options for BoomerAMG */
55   PetscInt  cycletype;
56   PetscInt  maxlevels;
57   double    strongthreshold;
58   double    maxrowsum;
59   PetscInt  gridsweeps[3];
60   PetscInt  coarsentype;
61   PetscInt  measuretype;
62   PetscInt  relaxtype[3];
63   double    relaxweight;
64   double    outerrelaxweight;
65   PetscInt  relaxorder;
66   double    truncfactor;
67   PetscBool applyrichardson;
68   PetscInt  pmax;
69   PetscInt  interptype;
70   PetscInt  agg_nl;
71   PetscInt  agg_num_paths;
72   PetscInt  nodal_coarsen;
73   PetscBool nodal_relax;
74   PetscInt  nodal_relax_levels;
75 
76   /* options for AS (Auxiliary Space preconditioners) */
77   PetscInt  as_print;
78   PetscInt  as_max_iter;
79   PetscReal as_tol;
80   PetscInt  as_relax_type;
81   PetscInt  as_relax_times;
82   PetscReal as_relax_weight;
83   PetscReal as_omega;
84   PetscInt  as_amg_alpha_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for vector Poisson (AMS) or Curl problem (ADS) */
85   PetscReal as_amg_alpha_theta;   /* AMG strength for vector Poisson (AMS) or Curl problem (ADS) */
86   PetscInt  as_amg_beta_opts[5];  /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for scalar Poisson (AMS) or vector Poisson (ADS) */
87   PetscReal as_amg_beta_theta;    /* AMG strength for scalar Poisson (AMS) or vector Poisson (ADS)  */
88   PetscInt  ams_cycle_type;
89   PetscInt  ads_cycle_type;
90 
91   /* additional data */
92   HYPRE_IJVector coords[3];
93   HYPRE_IJVector constants[3];
94   HYPRE_IJMatrix G;
95   HYPRE_IJMatrix C;
96   HYPRE_IJMatrix alpha_Poisson;
97   HYPRE_IJMatrix beta_Poisson;
98   PetscBool      ams_beta_is_zero;
99 } PC_HYPRE;
100 
101 #undef __FUNCT__
102 #define __FUNCT__ "PCHYPREGetSolver"
103 PetscErrorCode PCHYPREGetSolver(PC pc,HYPRE_Solver *hsolver)
104 {
105   PC_HYPRE *jac = (PC_HYPRE*)pc->data;
106 
107   PetscFunctionBegin;
108   *hsolver = jac->hsolver;
109   PetscFunctionReturn(0);
110 }
111 
112 #undef __FUNCT__
113 #define __FUNCT__ "PCSetUp_HYPRE"
114 static PetscErrorCode PCSetUp_HYPRE(PC pc)
115 {
116   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
117   PetscErrorCode     ierr;
118   HYPRE_ParCSRMatrix hmat;
119   HYPRE_ParVector    bv,xv;
120   PetscInt           bs;
121 
122   PetscFunctionBegin;
123   if (!jac->hypre_type) {
124     ierr = PCHYPRESetType(pc,"boomeramg");CHKERRQ(ierr);
125   }
126 
127   if (pc->setupcalled) {
128     /* always destroy the old matrix and create a new memory;
129        hope this does not churn the memory too much. The problem
130        is I do not know if it is possible to put the matrix back to
131        its initial state so that we can directly copy the values
132        the second time through. */
133     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij));
134     jac->ij = 0;
135   }
136 
137   if (!jac->ij) { /* create the matrix the first time through */
138     ierr = MatHYPRE_IJMatrixCreate(pc->pmat,&jac->ij);CHKERRQ(ierr);
139   }
140   if (!jac->b) { /* create the vectors the first time through */
141     Vec x,b;
142     ierr = MatCreateVecs(pc->pmat,&x,&b);CHKERRQ(ierr);
143     ierr = VecHYPRE_IJVectorCreate(x,&jac->x);CHKERRQ(ierr);
144     ierr = VecHYPRE_IJVectorCreate(b,&jac->b);CHKERRQ(ierr);
145     ierr = VecDestroy(&x);CHKERRQ(ierr);
146     ierr = VecDestroy(&b);CHKERRQ(ierr);
147   }
148 
149   /* special case for BoomerAMG */
150   if (jac->setup == HYPRE_BoomerAMGSetup) {
151     ierr = MatGetBlockSize(pc->pmat,&bs);CHKERRQ(ierr);
152     if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
153   }
154 
155   /* special case for AMS */
156   if (jac->setup == HYPRE_AMSSetup) {
157     if (!jac->coords[0] && !jac->constants[0]) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs either coordinate vectors via PCSetCoordinates() or edge constant vectors via PCHYPRESetEdgeConstantVectors()");
158     if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs discrete gradient operator via PCHYPRESetDiscreteGradient");
159   }
160   /* special case for ADS */
161   if (jac->setup == HYPRE_ADSSetup) {
162     if (!jac->coords[0]) {
163       SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs coordinate vectors via PCSetCoordinates()");
164     } else if (!jac->coords[1] || !jac->coords[2]) {
165       SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner has been designed for three dimensional problems! For two dimensional problems, use HYPRE AMS instead");
166     }
167     if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs discrete gradient operator via PCHYPRESetDiscreteGradient");
168     if (!jac->C) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs discrete curl operator via PCHYPRESetDiscreteGradient");
169   }
170   ierr = MatHYPRE_IJMatrixCopy(pc->pmat,jac->ij);CHKERRQ(ierr);
171   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
172   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&bv));
173   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&xv));
174   PetscStackCall("HYPRE_SetupXXX",ierr = (*jac->setup)(jac->hsolver,hmat,bv,xv);CHKERRQ(ierr););
175   PetscFunctionReturn(0);
176 }
177 
178 /*
179     Replaces the address where the HYPRE vector points to its data with the address of
180   PETSc's data. Saves the old address so it can be reset when we are finished with it.
181   Allows use to get the data into a HYPRE vector without the cost of memcopies
182 */
183 #define HYPREReplacePointer(b,newvalue,savedvalue) { \
184     hypre_ParVector *par_vector   = (hypre_ParVector*)hypre_IJVectorObject(((hypre_IJVector*)b)); \
185     hypre_Vector    *local_vector = hypre_ParVectorLocalVector(par_vector); \
186     savedvalue         = local_vector->data; \
187     local_vector->data = newvalue;          \
188 }
189 
190 #undef __FUNCT__
191 #define __FUNCT__ "PCApply_HYPRE"
192 static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
193 {
194   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
195   PetscErrorCode     ierr;
196   HYPRE_ParCSRMatrix hmat;
197   PetscScalar        *xv;
198   const PetscScalar  *bv,*sbv;
199   HYPRE_ParVector    jbv,jxv;
200   PetscScalar        *sxv;
201   PetscInt           hierr;
202 
203   PetscFunctionBegin;
204   ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr);
205   if (!jac->applyrichardson) {ierr = VecSet(x,0.0);CHKERRQ(ierr);}
206   ierr = VecGetArrayRead(b,&bv);CHKERRQ(ierr);
207   ierr = VecGetArray(x,&xv);CHKERRQ(ierr);
208   HYPREReplacePointer(jac->b,(PetscScalar*)bv,sbv);
209   HYPREReplacePointer(jac->x,xv,sxv);
210 
211   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
212   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
213   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));
214   PetscStackCall("Hypre solve",hierr = (*jac->solve)(jac->hsolver,hmat,jbv,jxv);
215                                if (hierr && hierr != HYPRE_ERROR_CONV) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
216                                if (hierr) hypre__global_error = 0;);
217 
218   if (jac->setup == HYPRE_AMSSetup && jac->ams_beta_is_zero) {
219     PetscStackCall("HYPRE_AMSProjectOutGradients",ierr = HYPRE_AMSProjectOutGradients(jac->hsolver,jxv);CHKERRQ(ierr););
220   }
221   HYPREReplacePointer(jac->b,(PetscScalar*)sbv,bv);
222   HYPREReplacePointer(jac->x,sxv,xv);
223   ierr = VecRestoreArray(x,&xv);CHKERRQ(ierr);
224   ierr = VecRestoreArrayRead(b,&bv);CHKERRQ(ierr);
225   PetscFunctionReturn(0);
226 }
227 
228 #undef __FUNCT__
229 #define __FUNCT__ "PCDestroy_HYPRE"
230 static PetscErrorCode PCDestroy_HYPRE(PC pc)
231 {
232   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
233   PetscErrorCode ierr;
234 
235   PetscFunctionBegin;
236   if (jac->ij) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij));
237   if (jac->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->b));
238   if (jac->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->x));
239   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0]));
240   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1]));
241   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2]));
242   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0]));
243   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1]));
244   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2]));
245   if (jac->G) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->G));
246   if (jac->C) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->C));
247   if (jac->alpha_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->alpha_Poisson));
248   if (jac->beta_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->beta_Poisson));
249   if (jac->destroy) PetscStackCall("HYPRE_DestroyXXX",ierr = (*jac->destroy)(jac->hsolver);CHKERRQ(ierr););
250   ierr = PetscFree(jac->hypre_type);CHKERRQ(ierr);
251   if (jac->comm_hypre != MPI_COMM_NULL) { ierr = MPI_Comm_free(&(jac->comm_hypre));CHKERRQ(ierr);}
252   ierr = PetscFree(pc->data);CHKERRQ(ierr);
253 
254   ierr = PetscObjectChangeTypeName((PetscObject)pc,0);CHKERRQ(ierr);
255   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);CHKERRQ(ierr);
256   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);CHKERRQ(ierr);
257   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetCoordinates_C",NULL);CHKERRQ(ierr);
258   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",NULL);CHKERRQ(ierr);
259   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",NULL);CHKERRQ(ierr);
260   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetConstantEdgeVectors_C",NULL);CHKERRQ(ierr);
261   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetAlphaPoissonMatrix_C",NULL);CHKERRQ(ierr);
262   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetBetaPoissonMatrix_C",NULL);CHKERRQ(ierr);
263   PetscFunctionReturn(0);
264 }
265 
266 /* --------------------------------------------------------------------------------------------*/
267 #undef __FUNCT__
268 #define __FUNCT__ "PCSetFromOptions_HYPRE_Pilut"
269 static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PetscOptions *PetscOptionsObject,PC pc)
270 {
271   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
272   PetscErrorCode ierr;
273   PetscBool      flag;
274 
275   PetscFunctionBegin;
276   ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE Pilut Options");CHKERRQ(ierr);
277   ierr = PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);CHKERRQ(ierr);
278   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter));
279   ierr = PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);CHKERRQ(ierr);
280   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol));
281   ierr = PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);CHKERRQ(ierr);
282   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize));
283   ierr = PetscOptionsTail();CHKERRQ(ierr);
284   PetscFunctionReturn(0);
285 }
286 
287 #undef __FUNCT__
288 #define __FUNCT__ "PCView_HYPRE_Pilut"
289 static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
290 {
291   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
292   PetscErrorCode ierr;
293   PetscBool      iascii;
294 
295   PetscFunctionBegin;
296   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
297   if (iascii) {
298     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut preconditioning\n");CHKERRQ(ierr);
299     if (jac->maxiter != PETSC_DEFAULT) {
300       ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: maximum number of iterations %d\n",jac->maxiter);CHKERRQ(ierr);
301     } else {
302       ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default maximum number of iterations \n");CHKERRQ(ierr);
303     }
304     if (jac->tol != PETSC_DEFAULT) {
305       ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: drop tolerance %g\n",(double)jac->tol);CHKERRQ(ierr);
306     } else {
307       ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default drop tolerance \n");CHKERRQ(ierr);
308     }
309     if (jac->factorrowsize != PETSC_DEFAULT) {
310       ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: factor row size %d\n",jac->factorrowsize);CHKERRQ(ierr);
311     } else {
312       ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default factor row size \n");CHKERRQ(ierr);
313     }
314   }
315   PetscFunctionReturn(0);
316 }
317 
318 /* --------------------------------------------------------------------------------------------*/
319 
320 #undef __FUNCT__
321 #define __FUNCT__ "PCApplyTranspose_HYPRE_BoomerAMG"
322 static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
323 {
324   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
325   PetscErrorCode     ierr;
326   HYPRE_ParCSRMatrix hmat;
327   PetscScalar        *xv;
328   const PetscScalar  *bv;
329   HYPRE_ParVector    jbv,jxv;
330   PetscScalar        *sbv,*sxv;
331   PetscInt           hierr;
332 
333   PetscFunctionBegin;
334   ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr);
335   ierr = VecSet(x,0.0);CHKERRQ(ierr);
336   ierr = VecGetArrayRead(b,&bv);CHKERRQ(ierr);
337   ierr = VecGetArray(x,&xv);CHKERRQ(ierr);
338   HYPREReplacePointer(jac->b,(PetscScalar*)bv,sbv);
339   HYPREReplacePointer(jac->x,xv,sxv);
340 
341   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
342   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
343   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));
344 
345   hierr = HYPRE_BoomerAMGSolveT(jac->hsolver,hmat,jbv,jxv);
346   /* error code of 1 in BoomerAMG merely means convergence not achieved */
347   if (hierr && (hierr != 1)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
348   if (hierr) hypre__global_error = 0;
349 
350   HYPREReplacePointer(jac->b,sbv,bv);
351   HYPREReplacePointer(jac->x,sxv,xv);
352   ierr = VecRestoreArray(x,&xv);CHKERRQ(ierr);
353   ierr = VecRestoreArrayRead(b,&bv);CHKERRQ(ierr);
354   PetscFunctionReturn(0);
355 }
356 
357 /* static array length */
358 #define ALEN(a) (sizeof(a)/sizeof((a)[0]))
359 
360 static const char *HYPREBoomerAMGCycleType[]   = {"","V","W"};
361 static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
362 static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
363 /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
364 static const char *HYPREBoomerAMGRelaxType[]   = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi",
365                                                   "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi",
366                                                   "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination",
367                                                   "" /* 10 */, "" /* 11 */, "" /* 12 */, "" /* 13 */, "" /* 14 */,
368                                                   "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"};
369 static const char *HYPREBoomerAMGInterpType[]  = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
370                                                   "ext+i-cc", "standard", "standard-wts", "", "", "FF", "FF1"};
371 #undef __FUNCT__
372 #define __FUNCT__ "PCSetFromOptions_HYPRE_BoomerAMG"
373 static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PetscOptions *PetscOptionsObject,PC pc)
374 {
375   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
376   PetscErrorCode ierr;
377   PetscInt       n,indx,level;
378   PetscBool      flg, tmp_truth;
379   double         tmpdbl, twodbl[2];
380 
381   PetscFunctionBegin;
382   ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE BoomerAMG Options");CHKERRQ(ierr);
383   ierr = PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);CHKERRQ(ierr);
384   if (flg) {
385     jac->cycletype = indx+1;
386     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
387   }
388   ierr = PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);CHKERRQ(ierr);
389   if (flg) {
390     if (jac->maxlevels < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
391     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
392   }
393   ierr = PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);CHKERRQ(ierr);
394   if (flg) {
395     if (jac->maxiter < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
396     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
397   }
398   ierr = PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);CHKERRQ(ierr);
399   if (flg) {
400     if (jac->tol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance %g must be greater than or equal to zero",(double)jac->tol);
401     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
402   }
403 
404   ierr = PetscOptionsScalar("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);CHKERRQ(ierr);
405   if (flg) {
406     if (jac->truncfactor < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Truncation factor %g must be great than or equal zero",(double)jac->truncfactor);
407     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
408   }
409 
410   ierr = PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);CHKERRQ(ierr);
411   if (flg) {
412     if (jac->pmax < 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"P_max %g must be greater than or equal to zero",(double)jac->pmax);
413     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
414   }
415 
416   ierr = PetscOptionsInt("-pc_hypre_boomeramg_agg_nl","Number of levels of aggressive coarsening","None",jac->agg_nl,&jac->agg_nl,&flg);CHKERRQ(ierr);
417   if (flg) {
418     if (jac->agg_nl < 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %g must be greater than or equal to zero",(double)jac->agg_nl);
419 
420     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
421   }
422 
423 
424   ierr = PetscOptionsInt("-pc_hypre_boomeramg_agg_num_paths","Number of paths for aggressive coarsening","None",jac->agg_num_paths,&jac->agg_num_paths,&flg);CHKERRQ(ierr);
425   if (flg) {
426     if (jac->agg_num_paths < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of paths %g must be greater than or equal to 1",(double)jac->agg_num_paths);
427 
428     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
429   }
430 
431 
432   ierr = PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);CHKERRQ(ierr);
433   if (flg) {
434     if (jac->strongthreshold < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Strong threshold %g must be great than or equal zero",(double)jac->strongthreshold);
435     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
436   }
437   ierr = PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);CHKERRQ(ierr);
438   if (flg) {
439     if (jac->maxrowsum < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %g must be greater than zero",(double)jac->maxrowsum);
440     if (jac->maxrowsum > 1.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %g must be less than or equal one",(double)jac->maxrowsum);
441     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
442   }
443 
444   /* Grid sweeps */
445   ierr = PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);CHKERRQ(ierr);
446   if (flg) {
447     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx));
448     /* modify the jac structure so we can view the updated options with PC_View */
449     jac->gridsweeps[0] = indx;
450     jac->gridsweeps[1] = indx;
451     /*defaults coarse to 1 */
452     jac->gridsweeps[2] = 1;
453   }
454 
455   ierr = PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);CHKERRQ(ierr);
456   if (flg) {
457     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1));
458     jac->gridsweeps[0] = indx;
459   }
460   ierr = PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);CHKERRQ(ierr);
461   if (flg) {
462     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2));
463     jac->gridsweeps[1] = indx;
464   }
465   ierr = PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);CHKERRQ(ierr);
466   if (flg) {
467     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3));
468     jac->gridsweeps[2] = indx;
469   }
470 
471   /* Relax type */
472   ierr = PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);CHKERRQ(ierr);
473   if (flg) {
474     jac->relaxtype[0] = jac->relaxtype[1]  = indx;
475     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx));
476     /* by default, coarse type set to 9 */
477     jac->relaxtype[2] = 9;
478 
479   }
480   ierr = PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);CHKERRQ(ierr);
481   if (flg) {
482     jac->relaxtype[0] = indx;
483     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1));
484   }
485   ierr = PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);CHKERRQ(ierr);
486   if (flg) {
487     jac->relaxtype[1] = indx;
488     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2));
489   }
490   ierr = PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);CHKERRQ(ierr);
491   if (flg) {
492     jac->relaxtype[2] = indx;
493     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3));
494   }
495 
496   /* Relaxation Weight */
497   ierr = PetscOptionsReal("-pc_hypre_boomeramg_relax_weight_all","Relaxation weight for all levels (0 = hypre estimates, -k = determined with k CG steps)","None",jac->relaxweight, &tmpdbl,&flg);CHKERRQ(ierr);
498   if (flg) {
499     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl));
500     jac->relaxweight = tmpdbl;
501   }
502 
503   n         = 2;
504   twodbl[0] = twodbl[1] = 1.0;
505   ierr      = PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);CHKERRQ(ierr);
506   if (flg) {
507     if (n == 2) {
508       indx =  (int)PetscAbsReal(twodbl[1]);
509       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx));
510     } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relax weight level: you must provide 2 values separated by a comma (and no space), you provided %d",n);
511   }
512 
513   /* Outer relaxation Weight */
514   ierr = PetscOptionsReal("-pc_hypre_boomeramg_outer_relax_weight_all","Outer relaxation weight for all levels (-k = determined with k CG steps)","None",jac->outerrelaxweight, &tmpdbl,&flg);CHKERRQ(ierr);
515   if (flg) {
516     PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl));
517     jac->outerrelaxweight = tmpdbl;
518   }
519 
520   n         = 2;
521   twodbl[0] = twodbl[1] = 1.0;
522   ierr      = PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);CHKERRQ(ierr);
523   if (flg) {
524     if (n == 2) {
525       indx =  (int)PetscAbsReal(twodbl[1]);
526       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx));
527     } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relax weight outer level: You must provide 2 values separated by a comma (and no space), you provided %d",n);
528   }
529 
530   /* the Relax Order */
531   ierr = PetscOptionsBool("-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", PETSC_FALSE, &tmp_truth, &flg);CHKERRQ(ierr);
532 
533   if (flg && tmp_truth) {
534     jac->relaxorder = 0;
535     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
536   }
537   ierr = PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);CHKERRQ(ierr);
538   if (flg) {
539     jac->measuretype = indx;
540     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
541   }
542   /* update list length 3/07 */
543   ierr = PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);CHKERRQ(ierr);
544   if (flg) {
545     jac->coarsentype = indx;
546     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
547   }
548 
549   /* new 3/07 */
550   ierr = PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);CHKERRQ(ierr);
551   if (flg) {
552     jac->interptype = indx;
553     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
554   }
555 
556   ierr = PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);CHKERRQ(ierr);
557   if (flg) {
558     level = 3;
559     ierr = PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);CHKERRQ(ierr);
560 
561     jac->printstatistics = PETSC_TRUE;
562     PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level));
563   }
564 
565   ierr = PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);CHKERRQ(ierr);
566   if (flg) {
567     level = 3;
568     ierr = PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);CHKERRQ(ierr);
569 
570     jac->printstatistics = PETSC_TRUE;
571     PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level));
572   }
573 
574   ierr = PetscOptionsBool("-pc_hypre_boomeramg_nodal_coarsen", "HYPRE_BoomerAMGSetNodal()", "None", jac->nodal_coarsen ? PETSC_TRUE : PETSC_FALSE, &tmp_truth, &flg);CHKERRQ(ierr);
575   if (flg && tmp_truth) {
576     jac->nodal_coarsen = 1;
577     PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,1));
578   }
579 
580   ierr = PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);CHKERRQ(ierr);
581   if (flg && tmp_truth) {
582     PetscInt tmp_int;
583     ierr = PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);CHKERRQ(ierr);
584     if (flg) jac->nodal_relax_levels = tmp_int;
585     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6));
586     PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1));
587     PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0));
588     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels));
589   }
590 
591   ierr = PetscOptionsTail();CHKERRQ(ierr);
592   PetscFunctionReturn(0);
593 }
594 
595 #undef __FUNCT__
596 #define __FUNCT__ "PCApplyRichardson_HYPRE_BoomerAMG"
597 static PetscErrorCode PCApplyRichardson_HYPRE_BoomerAMG(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
598 {
599   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
600   PetscErrorCode ierr;
601   PetscInt       oits;
602 
603   PetscFunctionBegin;
604   ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr);
605   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter));
606   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol));
607   jac->applyrichardson = PETSC_TRUE;
608   ierr                 = PCApply_HYPRE(pc,b,y);CHKERRQ(ierr);
609   jac->applyrichardson = PETSC_FALSE;
610   PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
611   *outits = oits;
612   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
613   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
614   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
615   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
616   PetscFunctionReturn(0);
617 }
618 
619 
620 #undef __FUNCT__
621 #define __FUNCT__ "PCView_HYPRE_BoomerAMG"
622 static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
623 {
624   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
625   PetscErrorCode ierr;
626   PetscBool      iascii;
627 
628   PetscFunctionBegin;
629   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
630   if (iascii) {
631     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG preconditioning\n");CHKERRQ(ierr);
632     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);CHKERRQ(ierr);
633     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of levels %d\n",jac->maxlevels);CHKERRQ(ierr);
634     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of iterations PER hypre call %d\n",jac->maxiter);CHKERRQ(ierr);
635     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Convergence tolerance PER hypre call %g\n",(double)jac->tol);CHKERRQ(ierr);
636     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Threshold for strong coupling %g\n",(double)jac->strongthreshold);CHKERRQ(ierr);
637     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation truncation factor %g\n",(double)jac->truncfactor);CHKERRQ(ierr);
638     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation: max elements per row %d\n",jac->pmax);CHKERRQ(ierr);
639     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Number of levels of aggressive coarsening %d\n",jac->agg_nl);CHKERRQ(ierr);
640     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Number of paths for aggressive coarsening %d\n",jac->agg_num_paths);CHKERRQ(ierr);
641 
642     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum row sums %g\n",(double)jac->maxrowsum);CHKERRQ(ierr);
643 
644     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps down         %d\n",jac->gridsweeps[0]);CHKERRQ(ierr);
645     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps up           %d\n",jac->gridsweeps[1]);CHKERRQ(ierr);
646     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps on coarse    %d\n",jac->gridsweeps[2]);CHKERRQ(ierr);
647 
648     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax down          %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[0]]);CHKERRQ(ierr);
649     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax up            %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[1]]);CHKERRQ(ierr);
650     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax on coarse     %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[2]]);CHKERRQ(ierr);
651 
652     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax weight  (all)      %g\n",(double)jac->relaxweight);CHKERRQ(ierr);
653     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Outer relax weight (all) %g\n",(double)jac->outerrelaxweight);CHKERRQ(ierr);
654 
655     if (jac->relaxorder) {
656       ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Using CF-relaxation\n");CHKERRQ(ierr);
657     } else {
658       ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Not using CF-relaxation\n");CHKERRQ(ierr);
659     }
660     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Measure type        %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);CHKERRQ(ierr);
661     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Coarsen type        %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);CHKERRQ(ierr);
662     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation type  %s\n",HYPREBoomerAMGInterpType[jac->interptype]);CHKERRQ(ierr);
663     if (jac->nodal_coarsen) {
664       ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal coarsening (with HYPRE_BOOMERAMGSetNodal())\n");CHKERRQ(ierr);
665     }
666     if (jac->nodal_relax) {
667       ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal relaxation via Schwarz smoothing on levels %d\n",jac->nodal_relax_levels);CHKERRQ(ierr);
668     }
669   }
670   PetscFunctionReturn(0);
671 }
672 
673 /* --------------------------------------------------------------------------------------------*/
674 #undef __FUNCT__
675 #define __FUNCT__ "PCSetFromOptions_HYPRE_ParaSails"
676 static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PetscOptions *PetscOptionsObject,PC pc)
677 {
678   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
679   PetscErrorCode ierr;
680   PetscInt       indx;
681   PetscBool      flag;
682   const char     *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};
683 
684   PetscFunctionBegin;
685   ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE ParaSails Options");CHKERRQ(ierr);
686   ierr = PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);CHKERRQ(ierr);
687   ierr = PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshhold,&jac->threshhold,&flag);CHKERRQ(ierr);
688   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
689 
690   ierr = PetscOptionsReal("-pc_hypre_parasails_filter","filter","None",jac->filter,&jac->filter,&flag);CHKERRQ(ierr);
691   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
692 
693   ierr = PetscOptionsReal("-pc_hypre_parasails_loadbal","Load balance","None",jac->loadbal,&jac->loadbal,&flag);CHKERRQ(ierr);
694   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
695 
696   ierr = PetscOptionsBool("-pc_hypre_parasails_logging","Print info to screen","None",(PetscBool)jac->logging,(PetscBool*)&jac->logging,&flag);CHKERRQ(ierr);
697   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
698 
699   ierr = PetscOptionsBool("-pc_hypre_parasails_reuse","Reuse nonzero pattern in preconditioner","None",(PetscBool)jac->ruse,(PetscBool*)&jac->ruse,&flag);CHKERRQ(ierr);
700   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
701 
702   ierr = PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);CHKERRQ(ierr);
703   if (flag) {
704     jac->symt = indx;
705     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
706   }
707 
708   ierr = PetscOptionsTail();CHKERRQ(ierr);
709   PetscFunctionReturn(0);
710 }
711 
712 #undef __FUNCT__
713 #define __FUNCT__ "PCView_HYPRE_ParaSails"
714 static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
715 {
716   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
717   PetscErrorCode ierr;
718   PetscBool      iascii;
719   const char     *symt = 0;;
720 
721   PetscFunctionBegin;
722   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
723   if (iascii) {
724     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails preconditioning\n");CHKERRQ(ierr);
725     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: nlevels %d\n",jac->nlevels);CHKERRQ(ierr);
726     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: threshold %g\n",(double)jac->threshhold);CHKERRQ(ierr);
727     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: filter %g\n",(double)jac->filter);CHKERRQ(ierr);
728     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: load balance %g\n",(double)jac->loadbal);CHKERRQ(ierr);
729     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: reuse nonzero structure %s\n",PetscBools[jac->ruse]);CHKERRQ(ierr);
730     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: print info to screen %s\n",PetscBools[jac->logging]);CHKERRQ(ierr);
731     if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
732     else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
733     else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
734     else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
735     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: %s\n",symt);CHKERRQ(ierr);
736   }
737   PetscFunctionReturn(0);
738 }
739 /* --------------------------------------------------------------------------------------------*/
740 #undef __FUNCT__
741 #define __FUNCT__ "PCSetFromOptions_HYPRE_AMS"
742 static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PetscOptions *PetscOptionsObject,PC pc)
743 {
744   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
745   PetscErrorCode ierr;
746   PetscInt       n;
747   PetscBool      flag,flag2,flag3,flag4;
748 
749   PetscFunctionBegin;
750   ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE AMS Options");CHKERRQ(ierr);
751   ierr = PetscOptionsInt("-pc_hypre_ams_print_level","Debugging output level for AMS","None",jac->as_print,&jac->as_print,&flag);CHKERRQ(ierr);
752   if (flag) PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
753   ierr = PetscOptionsInt("-pc_hypre_ams_max_iter","Maximum number of AMS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);CHKERRQ(ierr);
754   if (flag) PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
755   ierr = PetscOptionsInt("-pc_hypre_ams_cycle_type","Cycle type for AMS multigrid","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag);CHKERRQ(ierr);
756   if (flag) PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
757   ierr = PetscOptionsReal("-pc_hypre_ams_tol","Error tolerance for AMS multigrid","None",jac->as_tol,&jac->as_tol,&flag);CHKERRQ(ierr);
758   if (flag) PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
759   ierr = PetscOptionsInt("-pc_hypre_ams_relax_type","Relaxation type for AMS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);CHKERRQ(ierr);
760   ierr = PetscOptionsInt("-pc_hypre_ams_relax_times","Number of relaxation steps for AMS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);CHKERRQ(ierr);
761   ierr = PetscOptionsReal("-pc_hypre_ams_relax_weight","Relaxation weight for AMS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);CHKERRQ(ierr);
762   ierr = PetscOptionsReal("-pc_hypre_ams_omega","SSOR coefficient for AMS smoother","None",jac->as_omega,&jac->as_omega,&flag4);CHKERRQ(ierr);
763   if (flag || flag2 || flag3 || flag4) {
764     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
765                                                                       jac->as_relax_times,
766                                                                       jac->as_relax_weight,
767                                                                       jac->as_omega));
768   }
769   ierr = PetscOptionsReal("-pc_hypre_ams_amg_alpha_theta","Threshold for strong coupling of vector Poisson AMG solver","None",jac->as_amg_alpha_theta,&jac->as_amg_alpha_theta,&flag);CHKERRQ(ierr);
770   n = 5;
771   ierr = PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options","AMG options for vector Poisson","None",jac->as_amg_alpha_opts,&n,&flag2);CHKERRQ(ierr);
772   if (flag || flag2) {
773     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
774                                                                      jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
775                                                                      jac->as_amg_alpha_opts[2],       /* AMG relax_type */
776                                                                      jac->as_amg_alpha_theta,
777                                                                      jac->as_amg_alpha_opts[3],       /* AMG interp_type */
778                                                                      jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
779   }
780   ierr = PetscOptionsReal("-pc_hypre_ams_amg_beta_theta","Threshold for strong coupling of scalar Poisson AMG solver","None",jac->as_amg_beta_theta,&jac->as_amg_beta_theta,&flag);CHKERRQ(ierr);
781   n = 5;
782   ierr = PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options","AMG options for scalar Poisson solver","None",jac->as_amg_beta_opts,&n,&flag2);CHKERRQ(ierr);
783   if (flag || flag2) {
784     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
785                                                                     jac->as_amg_beta_opts[1],       /* AMG agg_levels */
786                                                                     jac->as_amg_beta_opts[2],       /* AMG relax_type */
787                                                                     jac->as_amg_beta_theta,
788                                                                     jac->as_amg_beta_opts[3],       /* AMG interp_type */
789                                                                     jac->as_amg_beta_opts[4]));     /* AMG Pmax */
790   }
791   ierr = PetscOptionsTail();CHKERRQ(ierr);
792   PetscFunctionReturn(0);
793 }
794 
795 #undef __FUNCT__
796 #define __FUNCT__ "PCView_HYPRE_AMS"
797 static PetscErrorCode PCView_HYPRE_AMS(PC pc,PetscViewer viewer)
798 {
799   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
800   PetscErrorCode ierr;
801   PetscBool      iascii;
802 
803   PetscFunctionBegin;
804   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
805   if (iascii) {
806     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS preconditioning\n");CHKERRQ(ierr);
807     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace iterations per application %d\n",jac->as_max_iter);CHKERRQ(ierr);
808     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace cycle type %d\n",jac->ams_cycle_type);CHKERRQ(ierr);
809     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace iteration tolerance %g\n",jac->as_tol);CHKERRQ(ierr);
810     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother type %d\n",jac->as_relax_type);CHKERRQ(ierr);
811     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: number of smoothing steps %d\n",jac->as_relax_times);CHKERRQ(ierr);
812     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother weight %g\n",jac->as_relax_weight);CHKERRQ(ierr);
813     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother omega %g\n",jac->as_omega);CHKERRQ(ierr);
814     if (jac->alpha_Poisson) {
815       ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: vector Poisson solver (passed in by user)\n");CHKERRQ(ierr);
816     } else {
817       ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: vector Poisson solver (computed) \n");CHKERRQ(ierr);
818     }
819     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);CHKERRQ(ierr);
820     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);CHKERRQ(ierr);
821     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);CHKERRQ(ierr);
822     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);CHKERRQ(ierr);
823     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);CHKERRQ(ierr);
824     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);CHKERRQ(ierr);
825     if (!jac->ams_beta_is_zero) {
826       if (jac->beta_Poisson) {
827         ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: scalar Poisson solver (passed in by user)\n");CHKERRQ(ierr);
828       } else {
829         ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: scalar Poisson solver (computed) \n");CHKERRQ(ierr);
830       }
831       ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);CHKERRQ(ierr);
832       ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);CHKERRQ(ierr);
833       ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);CHKERRQ(ierr);
834       ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);CHKERRQ(ierr);
835       ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);CHKERRQ(ierr);
836       ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);CHKERRQ(ierr);
837     }
838   }
839   PetscFunctionReturn(0);
840 }
841 
842 #undef __FUNCT__
843 #define __FUNCT__ "PCSetFromOptions_HYPRE_ADS"
844 static PetscErrorCode PCSetFromOptions_HYPRE_ADS(PetscOptions *PetscOptionsObject,PC pc)
845 {
846   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
847   PetscErrorCode ierr;
848   PetscInt       n;
849   PetscBool      flag,flag2,flag3,flag4;
850 
851   PetscFunctionBegin;
852   ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE ADS Options");CHKERRQ(ierr);
853   ierr = PetscOptionsInt("-pc_hypre_ads_print_level","Debugging output level for ADS","None",jac->as_print,&jac->as_print,&flag);CHKERRQ(ierr);
854   if (flag) PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
855   ierr = PetscOptionsInt("-pc_hypre_ads_max_iter","Maximum number of ADS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);CHKERRQ(ierr);
856   if (flag) PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
857   ierr = PetscOptionsInt("-pc_hypre_ads_cycle_type","Cycle type for ADS multigrid","None",jac->ads_cycle_type,&jac->ads_cycle_type,&flag);CHKERRQ(ierr);
858   if (flag) PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ads_cycle_type));
859   ierr = PetscOptionsReal("-pc_hypre_ads_tol","Error tolerance for ADS multigrid","None",jac->as_tol,&jac->as_tol,&flag);CHKERRQ(ierr);
860   if (flag) PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
861   ierr = PetscOptionsInt("-pc_hypre_ads_relax_type","Relaxation type for ADS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);CHKERRQ(ierr);
862   ierr = PetscOptionsInt("-pc_hypre_ads_relax_times","Number of relaxation steps for ADS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);CHKERRQ(ierr);
863   ierr = PetscOptionsReal("-pc_hypre_ads_relax_weight","Relaxation weight for ADS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);CHKERRQ(ierr);
864   ierr = PetscOptionsReal("-pc_hypre_ads_omega","SSOR coefficient for ADS smoother","None",jac->as_omega,&jac->as_omega,&flag4);CHKERRQ(ierr);
865   if (flag || flag2 || flag3 || flag4) {
866     PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
867                                                                       jac->as_relax_times,
868                                                                       jac->as_relax_weight,
869                                                                       jac->as_omega));
870   }
871   ierr = PetscOptionsReal("-pc_hypre_ads_ams_theta","Threshold for strong coupling of AMS solver inside ADS","None",jac->as_amg_alpha_theta,&jac->as_amg_alpha_theta,&flag);CHKERRQ(ierr);
872   n = 5;
873   ierr = PetscOptionsIntArray("-pc_hypre_ads_ams_options","AMG options for AMS solver inside ADS","None",jac->as_amg_alpha_opts,&n,&flag2);CHKERRQ(ierr);
874   ierr = PetscOptionsInt("-pc_hypre_ads_ams_cycle_type","Cycle type for AMS solver inside ADS","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag3);CHKERRQ(ierr);
875   if (flag || flag2 || flag3) {
876     PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type,             /* AMS cycle type */
877                                                                 jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
878                                                                 jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
879                                                                 jac->as_amg_alpha_opts[2],       /* AMG relax_type */
880                                                                 jac->as_amg_alpha_theta,
881                                                                 jac->as_amg_alpha_opts[3],       /* AMG interp_type */
882                                                                 jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
883   }
884   ierr = PetscOptionsReal("-pc_hypre_ads_amg_theta","Threshold for strong coupling of vector AMG solver inside ADS","None",jac->as_amg_beta_theta,&jac->as_amg_beta_theta,&flag);CHKERRQ(ierr);
885   n = 5;
886   ierr = PetscOptionsIntArray("-pc_hypre_ads_amg_options","AMG options for vector AMG solver inside ADS","None",jac->as_amg_beta_opts,&n,&flag2);CHKERRQ(ierr);
887   if (flag || flag2) {
888     PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
889                                                                 jac->as_amg_beta_opts[1],       /* AMG agg_levels */
890                                                                 jac->as_amg_beta_opts[2],       /* AMG relax_type */
891                                                                 jac->as_amg_beta_theta,
892                                                                 jac->as_amg_beta_opts[3],       /* AMG interp_type */
893                                                                 jac->as_amg_beta_opts[4]));     /* AMG Pmax */
894   }
895   ierr = PetscOptionsTail();CHKERRQ(ierr);
896   PetscFunctionReturn(0);
897 }
898 
899 #undef __FUNCT__
900 #define __FUNCT__ "PCView_HYPRE_ADS"
901 static PetscErrorCode PCView_HYPRE_ADS(PC pc,PetscViewer viewer)
902 {
903   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
904   PetscErrorCode ierr;
905   PetscBool      iascii;
906 
907   PetscFunctionBegin;
908   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
909   if (iascii) {
910     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS preconditioning\n");CHKERRQ(ierr);
911     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: subspace iterations per application %d\n",jac->as_max_iter);CHKERRQ(ierr);
912     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: subspace cycle type %d\n",jac->ads_cycle_type);CHKERRQ(ierr);
913     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: subspace iteration tolerance %g\n",jac->as_tol);CHKERRQ(ierr);
914     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: smoother type %d\n",jac->as_relax_type);CHKERRQ(ierr);
915     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: number of smoothing steps %d\n",jac->as_relax_times);CHKERRQ(ierr);
916     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: smoother weight %g\n",jac->as_relax_weight);CHKERRQ(ierr);
917     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: smoother omega %g\n",jac->as_omega);CHKERRQ(ierr);
918     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: AMS solver\n");CHKERRQ(ierr);
919     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     subspace cycle type %d\n",jac->ams_cycle_type);CHKERRQ(ierr);
920     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);CHKERRQ(ierr);
921     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);CHKERRQ(ierr);
922     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);CHKERRQ(ierr);
923     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);CHKERRQ(ierr);
924     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);CHKERRQ(ierr);
925     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);CHKERRQ(ierr);
926     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: vector Poisson solver\n");CHKERRQ(ierr);
927     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);CHKERRQ(ierr);
928     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);CHKERRQ(ierr);
929     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);CHKERRQ(ierr);
930     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);CHKERRQ(ierr);
931     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);CHKERRQ(ierr);
932     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);CHKERRQ(ierr);
933   }
934   PetscFunctionReturn(0);
935 }
936 
937 #undef __FUNCT__
938 #define __FUNCT__ "PCHYPRESetDiscreteGradient_HYPRE"
939 static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE(PC pc, Mat G)
940 {
941   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
942   HYPRE_ParCSRMatrix parcsr_G;
943   PetscErrorCode     ierr;
944 
945   PetscFunctionBegin;
946   /* throw away any discrete gradient if already set */
947   if (jac->G) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->G));
948   ierr = MatHYPRE_IJMatrixCreate(G,&jac->G);CHKERRQ(ierr);
949   ierr = MatHYPRE_IJMatrixCopy(G,jac->G);CHKERRQ(ierr);
950   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->G,(void**)(&parcsr_G)));
951   PetscStackCall("Hypre set gradient",ierr = (*jac->setdgrad)(jac->hsolver,parcsr_G);CHKERRQ(ierr););
952   PetscFunctionReturn(0);
953 }
954 
955 #undef __FUNCT__
956 #define __FUNCT__ "PCHYPRESetDiscreteGradient"
957 /*@
958  PCHYPRESetDiscreteGradient - Set discrete gradient matrix
959 
960    Collective on PC
961 
962    Input Parameters:
963 +  pc - the preconditioning context
964 -  G - the discrete gradient
965 
966    Level: intermediate
967 
968    Notes: G should have as many rows as the number of edges and as many columns as the number of vertices in the mesh
969           Each row of G has 2 nonzeros, with column indexes being the global indexes of edge's endpoints: matrix entries are +1 and -1 depending on edge orientation
970 
971 .seealso:
972 @*/
973 PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G)
974 {
975   PetscErrorCode ierr;
976 
977   PetscFunctionBegin;
978   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
979   PetscValidHeaderSpecific(G,MAT_CLASSID,2);
980   PetscCheckSameComm(pc,1,G,2);
981   ierr = PetscTryMethod(pc,"PCHYPRESetDiscreteGradient_C",(PC,Mat),(pc,G));CHKERRQ(ierr);
982   PetscFunctionReturn(0);
983 }
984 
985 #undef __FUNCT__
986 #define __FUNCT__ "PCHYPRESetDiscreteCurl_HYPRE"
987 static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C)
988 {
989   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
990   HYPRE_ParCSRMatrix parcsr_C;
991   PetscErrorCode     ierr;
992 
993   PetscFunctionBegin;
994   /* throw away any discrete curl if already set */
995   if (jac->C) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->C));
996   ierr = MatHYPRE_IJMatrixCreate(C,&jac->C);CHKERRQ(ierr);
997   ierr = MatHYPRE_IJMatrixCopy(C,jac->C);CHKERRQ(ierr);
998   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->C,(void**)(&parcsr_C)));
999   PetscStackCall("Hypre set curl",ierr = (*jac->setdcurl)(jac->hsolver,parcsr_C);CHKERRQ(ierr););
1000   PetscFunctionReturn(0);
1001 }
1002 
1003 #undef __FUNCT__
1004 #define __FUNCT__ "PCHYPRESetDiscreteCurl"
1005 /*@
1006  PCHYPRESetDiscreteCurl - Set discrete curl matrix
1007 
1008    Collective on PC
1009 
1010    Input Parameters:
1011 +  pc - the preconditioning context
1012 -  C - the discrete curl
1013 
1014    Level: intermediate
1015 
1016    Notes: C should have as many rows as the number of faces and as many columns as the number of edges in the mesh
1017           Each row of G has as many nonzeros as the number of edges of a face, with column indexes being the global indexes of the corresponding edge: matrix entries are +1 and -1 depending on edge orientation with respect to the face orientation
1018 
1019 .seealso:
1020 @*/
1021 PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C)
1022 {
1023   PetscErrorCode ierr;
1024 
1025   PetscFunctionBegin;
1026   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1027   PetscValidHeaderSpecific(C,MAT_CLASSID,2);
1028   PetscCheckSameComm(pc,1,C,2);
1029   ierr = PetscTryMethod(pc,"PCHYPRESetDiscreteCurl_C",(PC,Mat),(pc,C));CHKERRQ(ierr);
1030   PetscFunctionReturn(0);
1031 }
1032 
1033 #undef __FUNCT__
1034 #define __FUNCT__ "PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS"
1035 static PetscErrorCode PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS(PC pc, Mat A)
1036 {
1037   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1038   HYPRE_ParCSRMatrix parcsr_alpha_Poisson;
1039   PetscErrorCode     ierr;
1040 
1041   PetscFunctionBegin;
1042   /* throw away any matrix if already set */
1043   if (jac->alpha_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->alpha_Poisson));
1044   ierr = MatHYPRE_IJMatrixCreate(A,&jac->alpha_Poisson);CHKERRQ(ierr);
1045   ierr = MatHYPRE_IJMatrixCopy(A,jac->alpha_Poisson);CHKERRQ(ierr);
1046   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->alpha_Poisson,(void**)(&parcsr_alpha_Poisson)));
1047   PetscStackCallStandard(HYPRE_AMSSetAlphaPoissonMatrix,(jac->hsolver,parcsr_alpha_Poisson));
1048   PetscFunctionReturn(0);
1049 }
1050 
1051 #undef __FUNCT__
1052 #define __FUNCT__ "PCHYPRESetAlphaPoissonMatrix"
1053 /*@
1054  PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix
1055 
1056    Collective on PC
1057 
1058    Input Parameters:
1059 +  pc - the preconditioning context
1060 -  A - the matrix
1061 
1062    Level: intermediate
1063 
1064    Notes: A should be obtained by discretizing the vector valued Poisson problem with linear finite elements
1065 
1066 .seealso:
1067 @*/
1068 PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A)
1069 {
1070   PetscErrorCode ierr;
1071 
1072   PetscFunctionBegin;
1073   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1074   PetscValidHeaderSpecific(A,MAT_CLASSID,2);
1075   PetscCheckSameComm(pc,1,A,2);
1076   ierr = PetscTryMethod(pc,"PCHYPRESetAlphaPoissonMatrix_C",(PC,Mat),(pc,A));CHKERRQ(ierr);
1077   PetscFunctionReturn(0);
1078 }
1079 
1080 #undef __FUNCT__
1081 #define __FUNCT__ "PCHYPRESetBetaPoissonMatrix_HYPRE_AMS"
1082 static PetscErrorCode PCHYPRESetBetaPoissonMatrix_HYPRE_AMS(PC pc, Mat A)
1083 {
1084   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1085   HYPRE_ParCSRMatrix parcsr_beta_Poisson;
1086   PetscErrorCode     ierr;
1087 
1088   PetscFunctionBegin;
1089   if (!A) {
1090     PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,NULL));
1091     jac->ams_beta_is_zero = PETSC_TRUE;
1092     PetscFunctionReturn(0);
1093   }
1094   jac->ams_beta_is_zero = PETSC_FALSE;
1095   /* throw away any matrix if already set */
1096   if (jac->beta_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->beta_Poisson));
1097   ierr = MatHYPRE_IJMatrixCreate(A,&jac->beta_Poisson);CHKERRQ(ierr);
1098   ierr = MatHYPRE_IJMatrixCopy(A,jac->beta_Poisson);CHKERRQ(ierr);
1099   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->beta_Poisson,(void**)(&parcsr_beta_Poisson)));
1100   PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,parcsr_beta_Poisson));
1101   PetscFunctionReturn(0);
1102 }
1103 
1104 #undef __FUNCT__
1105 #define __FUNCT__ "PCHYPRESetBetaPoissonMatrix"
1106 /*@
1107  PCHYPRESetBetaPoissonMatrix - Set Poisson matrix
1108 
1109    Collective on PC
1110 
1111    Input Parameters:
1112 +  pc - the preconditioning context
1113 -  A - the matrix
1114 
1115    Level: intermediate
1116 
1117    Notes: A should be obtained by discretizing the Poisson problem with linear finite elements.
1118           Following HYPRE convention, the scalar Poisson solver of AMS can be turned off by passing NULL.
1119 
1120 .seealso:
1121 @*/
1122 PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A)
1123 {
1124   PetscErrorCode ierr;
1125 
1126   PetscFunctionBegin;
1127   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1128   if (A) {
1129     PetscValidHeaderSpecific(A,MAT_CLASSID,2);
1130     PetscCheckSameComm(pc,1,A,2);
1131   }
1132   ierr = PetscTryMethod(pc,"PCHYPRESetBetaPoissonMatrix_C",(PC,Mat),(pc,A));CHKERRQ(ierr);
1133   PetscFunctionReturn(0);
1134 }
1135 
1136 #undef __FUNCT__
1137 #define __FUNCT__ "PCHYPRESetEdgeConstantVectors_HYPRE_AMS"
1138 static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE_AMS(PC pc,Vec ozz, Vec zoz, Vec zzo)
1139 {
1140   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1141   HYPRE_ParVector    par_ozz,par_zoz,par_zzo;
1142   PetscInt           dim;
1143   PetscErrorCode     ierr;
1144 
1145   PetscFunctionBegin;
1146   /* throw away any vector if already set */
1147   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0]));
1148   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1]));
1149   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2]));
1150   jac->constants[0] = NULL;
1151   jac->constants[1] = NULL;
1152   jac->constants[2] = NULL;
1153   ierr = VecHYPRE_IJVectorCreate(ozz,&jac->constants[0]);CHKERRQ(ierr);
1154   ierr = VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);CHKERRQ(ierr);
1155   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[0],(void**)(&par_ozz)));
1156   ierr = VecHYPRE_IJVectorCreate(zoz,&jac->constants[1]);CHKERRQ(ierr);
1157   ierr = VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);CHKERRQ(ierr);
1158   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[1],(void**)(&par_zoz)));
1159   dim = 2;
1160   if (zzo) {
1161     ierr = VecHYPRE_IJVectorCreate(zzo,&jac->constants[2]);CHKERRQ(ierr);
1162     ierr = VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);CHKERRQ(ierr);
1163     PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[2],(void**)(&par_zzo)));
1164     dim++;
1165   }
1166   PetscStackCallStandard(HYPRE_AMSSetEdgeConstantVectors,(jac->hsolver,par_ozz,par_zoz,par_zzo));
1167   PetscStackCallStandard(HYPRE_AMSSetDimension,(jac->hsolver,dim));
1168   PetscFunctionReturn(0);
1169 }
1170 
1171 #undef __FUNCT__
1172 #define __FUNCT__ "PCHYPRESetEdgeConstantVectors"
1173 /*@
1174  PCHYPRESetEdgeConstantVectors - Set the representation of the constant vector fields in edge element basis
1175 
1176    Collective on PC
1177 
1178    Input Parameters:
1179 +  pc - the preconditioning context
1180 -  ozz - vector representing (1,0,0) (or (1,0) in 2D)
1181 -  zoz - vector representing (0,1,0) (or (0,1) in 2D)
1182 -  zzo - vector representing (0,0,1) (use NULL in 2D)
1183 
1184    Level: intermediate
1185 
1186    Notes:
1187 
1188 .seealso:
1189 @*/
1190 PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo)
1191 {
1192   PetscErrorCode ierr;
1193 
1194   PetscFunctionBegin;
1195   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1196   PetscValidHeaderSpecific(ozz,VEC_CLASSID,2);
1197   PetscValidHeaderSpecific(zoz,VEC_CLASSID,3);
1198   if (zzo) PetscValidHeaderSpecific(zzo,VEC_CLASSID,4);
1199   PetscCheckSameComm(pc,1,ozz,2);
1200   PetscCheckSameComm(pc,1,zoz,3);
1201   if (zzo) PetscCheckSameComm(pc,1,zzo,4);
1202   ierr = PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));CHKERRQ(ierr);
1203   PetscFunctionReturn(0);
1204 }
1205 
1206 #undef __FUNCT__
1207 #define __FUNCT__ "PCSetCoordinates_HYPRE"
1208 static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1209 {
1210   PC_HYPRE        *jac = (PC_HYPRE*)pc->data;
1211   Vec             tv;
1212   HYPRE_ParVector par_coords[3];
1213   PetscInt        i;
1214   PetscErrorCode  ierr;
1215 
1216   PetscFunctionBegin;
1217   /* throw away any coordinate vector if already set */
1218   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0]));
1219   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1]));
1220   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2]));
1221   /* set problem's dimension */
1222   if (jac->setdim) {
1223     PetscStackCall("Hypre set dim",ierr = (*jac->setdim)(jac->hsolver,dim);CHKERRQ(ierr););
1224   }
1225   /* compute IJ vector for coordinates */
1226   ierr = VecCreate(PetscObjectComm((PetscObject)pc),&tv);CHKERRQ(ierr);
1227   ierr = VecSetType(tv,VECSTANDARD);CHKERRQ(ierr);
1228   ierr = VecSetSizes(tv,nloc,PETSC_DECIDE);CHKERRQ(ierr);
1229   for (i=0;i<dim;i++) {
1230     PetscScalar *array;
1231     PetscInt    j;
1232 
1233     ierr = VecHYPRE_IJVectorCreate(tv,&jac->coords[i]);CHKERRQ(ierr);
1234     ierr = VecGetArray(tv,&array);CHKERRQ(ierr);
1235     for (j=0;j<nloc;j++) {
1236       array[j] = coords[j*dim+i];
1237     }
1238     PetscStackCallStandard(HYPRE_IJVectorSetValues,(jac->coords[i],nloc,NULL,array));
1239     PetscStackCallStandard(HYPRE_IJVectorAssemble,(jac->coords[i]));
1240     ierr = VecRestoreArray(tv,&array);CHKERRQ(ierr);
1241   }
1242   ierr = VecDestroy(&tv);CHKERRQ(ierr);
1243   /* pass parCSR vectors to AMS solver */
1244   par_coords[0] = NULL;
1245   par_coords[1] = NULL;
1246   par_coords[2] = NULL;
1247   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0],(void**)(&par_coords[0])));
1248   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1],(void**)(&par_coords[1])));
1249   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2],(void**)(&par_coords[2])));
1250   PetscStackCall("Hypre set coords",ierr = (*jac->setcoord)(jac->hsolver,par_coords[0],par_coords[1],par_coords[2]);CHKERRQ(ierr););
1251   PetscFunctionReturn(0);
1252 }
1253 
1254 /* ---------------------------------------------------------------------------------*/
1255 
1256 #undef __FUNCT__
1257 #define __FUNCT__ "PCHYPREGetType_HYPRE"
1258 static PetscErrorCode  PCHYPREGetType_HYPRE(PC pc,const char *name[])
1259 {
1260   PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1261 
1262   PetscFunctionBegin;
1263   *name = jac->hypre_type;
1264   PetscFunctionReturn(0);
1265 }
1266 
1267 #undef __FUNCT__
1268 #define __FUNCT__ "PCHYPRESetType_HYPRE"
1269 static PetscErrorCode  PCHYPRESetType_HYPRE(PC pc,const char name[])
1270 {
1271   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1272   PetscErrorCode ierr;
1273   PetscBool      flag;
1274 
1275   PetscFunctionBegin;
1276   if (jac->hypre_type) {
1277     ierr = PetscStrcmp(jac->hypre_type,name,&flag);CHKERRQ(ierr);
1278     if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set");
1279     PetscFunctionReturn(0);
1280   } else {
1281     ierr = PetscStrallocpy(name, &jac->hypre_type);CHKERRQ(ierr);
1282   }
1283 
1284   jac->maxiter         = PETSC_DEFAULT;
1285   jac->tol             = PETSC_DEFAULT;
1286   jac->printstatistics = PetscLogPrintInfo;
1287 
1288   ierr = PetscStrcmp("pilut",jac->hypre_type,&flag);CHKERRQ(ierr);
1289   if (flag) {
1290     PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver));
1291     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
1292     pc->ops->view           = PCView_HYPRE_Pilut;
1293     jac->destroy            = HYPRE_ParCSRPilutDestroy;
1294     jac->setup              = HYPRE_ParCSRPilutSetup;
1295     jac->solve              = HYPRE_ParCSRPilutSolve;
1296     jac->factorrowsize      = PETSC_DEFAULT;
1297     PetscFunctionReturn(0);
1298   }
1299   ierr = PetscStrcmp("parasails",jac->hypre_type,&flag);CHKERRQ(ierr);
1300   if (flag) {
1301     PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver));
1302     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
1303     pc->ops->view           = PCView_HYPRE_ParaSails;
1304     jac->destroy            = HYPRE_ParaSailsDestroy;
1305     jac->setup              = HYPRE_ParaSailsSetup;
1306     jac->solve              = HYPRE_ParaSailsSolve;
1307     /* initialize */
1308     jac->nlevels    = 1;
1309     jac->threshhold = .1;
1310     jac->filter     = .1;
1311     jac->loadbal    = 0;
1312     if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE;
1313     else jac->logging = (int) PETSC_FALSE;
1314 
1315     jac->ruse = (int) PETSC_FALSE;
1316     jac->symt = 0;
1317     PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
1318     PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
1319     PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
1320     PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
1321     PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
1322     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
1323     PetscFunctionReturn(0);
1324   }
1325   ierr = PetscStrcmp("boomeramg",jac->hypre_type,&flag);CHKERRQ(ierr);
1326   if (flag) {
1327     ierr                     = HYPRE_BoomerAMGCreate(&jac->hsolver);
1328     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_BoomerAMG;
1329     pc->ops->view            = PCView_HYPRE_BoomerAMG;
1330     pc->ops->applytranspose  = PCApplyTranspose_HYPRE_BoomerAMG;
1331     pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
1332     jac->destroy             = HYPRE_BoomerAMGDestroy;
1333     jac->setup               = HYPRE_BoomerAMGSetup;
1334     jac->solve               = HYPRE_BoomerAMGSolve;
1335     jac->applyrichardson     = PETSC_FALSE;
1336     /* these defaults match the hypre defaults */
1337     jac->cycletype        = 1;
1338     jac->maxlevels        = 25;
1339     jac->maxiter          = 1;
1340     jac->tol              = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
1341     jac->truncfactor      = 0.0;
1342     jac->strongthreshold  = .25;
1343     jac->maxrowsum        = .9;
1344     jac->coarsentype      = 6;
1345     jac->measuretype      = 0;
1346     jac->gridsweeps[0]    = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
1347     jac->relaxtype[0]     = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */
1348     jac->relaxtype[2]     = 9; /*G.E. */
1349     jac->relaxweight      = 1.0;
1350     jac->outerrelaxweight = 1.0;
1351     jac->relaxorder       = 1;
1352     jac->interptype       = 0;
1353     jac->agg_nl           = 0;
1354     jac->pmax             = 0;
1355     jac->truncfactor      = 0.0;
1356     jac->agg_num_paths    = 1;
1357 
1358     jac->nodal_coarsen      = 0;
1359     jac->nodal_relax        = PETSC_FALSE;
1360     jac->nodal_relax_levels = 1;
1361     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
1362     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
1363     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
1364     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
1365     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
1366     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
1367     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
1368     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
1369     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
1370     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
1371     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
1372     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
1373     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
1374     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
1375     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0]));  /*defaults coarse to 9*/
1376     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */
1377     PetscFunctionReturn(0);
1378   }
1379   ierr = PetscStrcmp("ams",jac->hypre_type,&flag);CHKERRQ(ierr);
1380   if (flag) {
1381     ierr                     = HYPRE_AMSCreate(&jac->hsolver);
1382     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_AMS;
1383     pc->ops->view            = PCView_HYPRE_AMS;
1384     jac->destroy             = HYPRE_AMSDestroy;
1385     jac->setup               = HYPRE_AMSSetup;
1386     jac->solve               = HYPRE_AMSSolve;
1387     jac->setdgrad            = HYPRE_AMSSetDiscreteGradient;
1388     jac->setcoord            = HYPRE_AMSSetCoordinateVectors;
1389     jac->setdim              = HYPRE_AMSSetDimension;
1390     jac->coords[0]           = NULL;
1391     jac->coords[1]           = NULL;
1392     jac->coords[2]           = NULL;
1393     jac->G                   = NULL;
1394     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */
1395     jac->as_print           = 0;
1396     jac->as_max_iter        = 1; /* used as a preconditioner */
1397     jac->as_tol             = 0.; /* used as a preconditioner */
1398     jac->ams_cycle_type     = 13;
1399     /* Smoothing options */
1400     jac->as_relax_type      = 2;
1401     jac->as_relax_times     = 1;
1402     jac->as_relax_weight    = 1.0;
1403     jac->as_omega           = 1.0;
1404     /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1405     jac->as_amg_alpha_opts[0] = 10;
1406     jac->as_amg_alpha_opts[1] = 1;
1407     jac->as_amg_alpha_opts[2] = 8;
1408     jac->as_amg_alpha_opts[3] = 6;
1409     jac->as_amg_alpha_opts[4] = 4;
1410     jac->as_amg_alpha_theta   = 0.25;
1411     /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1412     jac->ams_beta_is_zero = PETSC_FALSE;
1413     jac->as_amg_beta_opts[0] = 10;
1414     jac->as_amg_beta_opts[1] = 1;
1415     jac->as_amg_beta_opts[2] = 8;
1416     jac->as_amg_beta_opts[3] = 6;
1417     jac->as_amg_beta_opts[4] = 4;
1418     jac->as_amg_beta_theta   = 0.25;
1419     PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
1420     PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1421     PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1422     PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
1423     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1424                                                                       jac->as_relax_times,
1425                                                                       jac->as_relax_weight,
1426                                                                       jac->as_omega));
1427     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1428                                                                      jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1429                                                                      jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1430                                                                      jac->as_amg_alpha_theta,
1431                                                                      jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1432                                                                      jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1433     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1434                                                                     jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1435                                                                     jac->as_amg_beta_opts[2],       /* AMG relax_type */
1436                                                                     jac->as_amg_beta_theta,
1437                                                                     jac->as_amg_beta_opts[3],       /* AMG interp_type */
1438                                                                     jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1439     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);CHKERRQ(ierr);
1440     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);CHKERRQ(ierr);
1441     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE_AMS);CHKERRQ(ierr);
1442     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetAlphaPoissonMatrix_C",PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS);CHKERRQ(ierr);
1443     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetBetaPoissonMatrix_C",PCHYPRESetBetaPoissonMatrix_HYPRE_AMS);CHKERRQ(ierr);
1444     PetscFunctionReturn(0);
1445   }
1446   ierr = PetscStrcmp("ads",jac->hypre_type,&flag);CHKERRQ(ierr);
1447   if (flag) {
1448     ierr                     = HYPRE_ADSCreate(&jac->hsolver);
1449     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_ADS;
1450     pc->ops->view            = PCView_HYPRE_ADS;
1451     jac->destroy             = HYPRE_ADSDestroy;
1452     jac->setup               = HYPRE_ADSSetup;
1453     jac->solve               = HYPRE_ADSSolve;
1454     jac->setdgrad            = HYPRE_ADSSetDiscreteGradient;
1455     jac->setdcurl            = HYPRE_ADSSetDiscreteCurl;
1456     jac->setcoord            = HYPRE_ADSSetCoordinateVectors;
1457     jac->coords[0]           = NULL;
1458     jac->coords[1]           = NULL;
1459     jac->coords[2]           = NULL;
1460     jac->G                   = NULL;
1461     jac->C                   = NULL;
1462     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */
1463     jac->as_print           = 0;
1464     jac->as_max_iter        = 1; /* used as a preconditioner */
1465     jac->as_tol             = 0.; /* used as a preconditioner */
1466     jac->ads_cycle_type     = 13;
1467     /* Smoothing options */
1468     jac->as_relax_type      = 2;
1469     jac->as_relax_times     = 1;
1470     jac->as_relax_weight    = 1.0;
1471     jac->as_omega           = 1.0;
1472     /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */
1473     jac->ams_cycle_type       = 14;
1474     jac->as_amg_alpha_opts[0] = 10;
1475     jac->as_amg_alpha_opts[1] = 1;
1476     jac->as_amg_alpha_opts[2] = 6;
1477     jac->as_amg_alpha_opts[3] = 6;
1478     jac->as_amg_alpha_opts[4] = 4;
1479     jac->as_amg_alpha_theta   = 0.25;
1480     /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1481     jac->as_amg_beta_opts[0] = 10;
1482     jac->as_amg_beta_opts[1] = 1;
1483     jac->as_amg_beta_opts[2] = 6;
1484     jac->as_amg_beta_opts[3] = 6;
1485     jac->as_amg_beta_opts[4] = 4;
1486     jac->as_amg_beta_theta   = 0.25;
1487     PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
1488     PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1489     PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1490     PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
1491     PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1492                                                                       jac->as_relax_times,
1493                                                                       jac->as_relax_weight,
1494                                                                       jac->as_omega));
1495     PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type,             /* AMG coarsen type */
1496                                                                 jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1497                                                                 jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1498                                                                 jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1499                                                                 jac->as_amg_alpha_theta,
1500                                                                 jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1501                                                                 jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1502     PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1503                                                                 jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1504                                                                 jac->as_amg_beta_opts[2],       /* AMG relax_type */
1505                                                                 jac->as_amg_beta_theta,
1506                                                                 jac->as_amg_beta_opts[3],       /* AMG interp_type */
1507                                                                 jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1508     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);CHKERRQ(ierr);
1509     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);CHKERRQ(ierr);
1510     ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",PCHYPRESetDiscreteCurl_HYPRE);CHKERRQ(ierr);
1511     PetscFunctionReturn(0);
1512   }
1513   ierr = PetscFree(jac->hypre_type);CHKERRQ(ierr);
1514 
1515   jac->hypre_type = NULL;
1516   SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are pilut, parasails, boomeramg, ams",name);
1517   PetscFunctionReturn(0);
1518 }
1519 
1520 /*
1521     It only gets here if the HYPRE type has not been set before the call to
1522    ...SetFromOptions() which actually is most of the time
1523 */
1524 #undef __FUNCT__
1525 #define __FUNCT__ "PCSetFromOptions_HYPRE"
1526 static PetscErrorCode PCSetFromOptions_HYPRE(PetscOptions *PetscOptionsObject,PC pc)
1527 {
1528   PetscErrorCode ierr;
1529   PetscInt       indx;
1530   const char     *type[] = {"pilut","parasails","boomeramg","ams","ads"};
1531   PetscBool      flg;
1532 
1533   PetscFunctionBegin;
1534   ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");CHKERRQ(ierr);
1535   ierr = PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"boomeramg",&indx,&flg);CHKERRQ(ierr);
1536   if (flg) {
1537     ierr = PCHYPRESetType_HYPRE(pc,type[indx]);CHKERRQ(ierr);
1538   } else {
1539     ierr = PCHYPRESetType_HYPRE(pc,"boomeramg");CHKERRQ(ierr);
1540   }
1541   if (pc->ops->setfromoptions) {
1542     ierr = pc->ops->setfromoptions(PetscOptionsObject,pc);CHKERRQ(ierr);
1543   }
1544   ierr = PetscOptionsTail();CHKERRQ(ierr);
1545   PetscFunctionReturn(0);
1546 }
1547 
1548 #undef __FUNCT__
1549 #define __FUNCT__ "PCHYPRESetType"
1550 /*@C
1551      PCHYPRESetType - Sets which hypre preconditioner you wish to use
1552 
1553    Input Parameters:
1554 +     pc - the preconditioner context
1555 -     name - either  pilut, parasails, boomeramg, ams, ads
1556 
1557    Options Database Keys:
1558    -pc_hypre_type - One of pilut, parasails, boomeramg, ams, ads
1559 
1560    Level: intermediate
1561 
1562 .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
1563            PCHYPRE
1564 
1565 @*/
1566 PetscErrorCode  PCHYPRESetType(PC pc,const char name[])
1567 {
1568   PetscErrorCode ierr;
1569 
1570   PetscFunctionBegin;
1571   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1572   PetscValidCharPointer(name,2);
1573   ierr = PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));CHKERRQ(ierr);
1574   PetscFunctionReturn(0);
1575 }
1576 
1577 #undef __FUNCT__
1578 #define __FUNCT__ "PCHYPREGetType"
1579 /*@C
1580      PCHYPREGetType - Gets which hypre preconditioner you are using
1581 
1582    Input Parameter:
1583 .     pc - the preconditioner context
1584 
1585    Output Parameter:
1586 .     name - either  pilut, parasails, boomeramg, ams, ads
1587 
1588    Level: intermediate
1589 
1590 .seealso:  PCCreate(), PCHYPRESetType(), PCType (for list of available types), PC,
1591            PCHYPRE
1592 
1593 @*/
1594 PetscErrorCode  PCHYPREGetType(PC pc,const char *name[])
1595 {
1596   PetscErrorCode ierr;
1597 
1598   PetscFunctionBegin;
1599   PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1600   PetscValidPointer(name,2);
1601   ierr = PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));CHKERRQ(ierr);
1602   PetscFunctionReturn(0);
1603 }
1604 
1605 /*MC
1606      PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre
1607 
1608    Options Database Keys:
1609 +   -pc_hypre_type - One of pilut, parasails, boomeramg, ams, ads
1610 -   Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX
1611           preconditioner
1612 
1613    Level: intermediate
1614 
1615    Notes: Apart from pc_hypre_type (for which there is PCHYPRESetType()),
1616           the many hypre options can ONLY be set via the options database (e.g. the command line
1617           or with PetscOptionsSetValue(), there are no functions to set them)
1618 
1619           The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_rtol refer to the number of iterations
1620           (V-cycles) and tolerance that boomeramg does EACH time it is called. So for example, if
1621           -pc_hypre_boomeramg_max_iter is set to 2 then 2-V-cycles are being used to define the preconditioner
1622           (-pc_hypre_boomeramg_rtol should be set to 0.0 - the default - to strictly use a fixed number of
1623           iterations per hypre call). -ksp_max_it and -ksp_rtol STILL determine the total number of iterations
1624           and tolerance for the Krylov solver. For example, if -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10
1625           then AT MOST twenty V-cycles of boomeramg will be called.
1626 
1627            Note that the option -pc_hypre_boomeramg_relax_type_all defaults to symmetric relaxation
1628            (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry.
1629            Otherwise, you may want to use -pc_hypre_boomeramg_relax_type_all SOR/Jacobi.
1630           If you wish to use BoomerAMG WITHOUT a Krylov method use -ksp_type richardson NOT -ksp_type preonly
1631           and use -ksp_max_it to control the number of V-cycles.
1632           (see the PETSc FAQ.html at the PETSc website under the Documentation tab).
1633 
1634           2007-02-03 Using HYPRE-1.11.1b, the routine HYPRE_BoomerAMGSolveT and the option
1635           -pc_hypre_parasails_reuse were failing with SIGSEGV. Dalcin L.
1636 
1637           See PCPFMG for access to the hypre Struct PFMG solver
1638 
1639 .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
1640            PCHYPRESetType(), PCPFMG
1641 
1642 M*/
1643 
1644 #undef __FUNCT__
1645 #define __FUNCT__ "PCCreate_HYPRE"
1646 PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
1647 {
1648   PC_HYPRE       *jac;
1649   PetscErrorCode ierr;
1650 
1651   PetscFunctionBegin;
1652   ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr);
1653 
1654   pc->data                = jac;
1655   pc->ops->destroy        = PCDestroy_HYPRE;
1656   pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
1657   pc->ops->setup          = PCSetUp_HYPRE;
1658   pc->ops->apply          = PCApply_HYPRE;
1659   jac->comm_hypre         = MPI_COMM_NULL;
1660   jac->hypre_type         = NULL;
1661   jac->coords[0]          = NULL;
1662   jac->coords[1]          = NULL;
1663   jac->coords[2]          = NULL;
1664   jac->constants[0]       = NULL;
1665   jac->constants[1]       = NULL;
1666   jac->constants[2]       = NULL;
1667   jac->G                  = NULL;
1668   jac->C                  = NULL;
1669   jac->alpha_Poisson      = NULL;
1670   jac->beta_Poisson       = NULL;
1671   jac->setdim             = NULL;
1672   /* duplicate communicator for hypre */
1673   ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));CHKERRQ(ierr);
1674   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);CHKERRQ(ierr);
1675   ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);CHKERRQ(ierr);
1676   PetscFunctionReturn(0);
1677 }
1678 
1679 /* ---------------------------------------------------------------------------------------------------------------------------------*/
1680 
1681 /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */
1682 #include <petsc/private/matimpl.h>
1683 
1684 typedef struct {
1685   MPI_Comm           hcomm;        /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */
1686   HYPRE_StructSolver hsolver;
1687 
1688   /* keep copy of PFMG options used so may view them */
1689   PetscInt its;
1690   double   tol;
1691   PetscInt relax_type;
1692   PetscInt rap_type;
1693   PetscInt num_pre_relax,num_post_relax;
1694   PetscInt max_levels;
1695 } PC_PFMG;
1696 
1697 #undef __FUNCT__
1698 #define __FUNCT__ "PCDestroy_PFMG"
1699 PetscErrorCode PCDestroy_PFMG(PC pc)
1700 {
1701   PetscErrorCode ierr;
1702   PC_PFMG        *ex = (PC_PFMG*) pc->data;
1703 
1704   PetscFunctionBegin;
1705   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1706   ierr = MPI_Comm_free(&ex->hcomm);CHKERRQ(ierr);
1707   ierr = PetscFree(pc->data);CHKERRQ(ierr);
1708   PetscFunctionReturn(0);
1709 }
1710 
1711 static const char *PFMGRelaxType[] = {"Jacobi","Weighted-Jacobi","symmetric-Red/Black-Gauss-Seidel","Red/Black-Gauss-Seidel"};
1712 static const char *PFMGRAPType[] = {"Galerkin","non-Galerkin"};
1713 
1714 #undef __FUNCT__
1715 #define __FUNCT__ "PCView_PFMG"
1716 PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
1717 {
1718   PetscErrorCode ierr;
1719   PetscBool      iascii;
1720   PC_PFMG        *ex = (PC_PFMG*) pc->data;
1721 
1722   PetscFunctionBegin;
1723   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
1724   if (iascii) {
1725     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG preconditioning\n");CHKERRQ(ierr);
1726     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: max iterations %d\n",ex->its);CHKERRQ(ierr);
1727     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: tolerance %g\n",ex->tol);CHKERRQ(ierr);
1728     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);CHKERRQ(ierr);
1729     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: RAP type %s\n",PFMGRAPType[ex->rap_type]);CHKERRQ(ierr);
1730     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);CHKERRQ(ierr);
1731     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: max levels %d\n",ex->max_levels);CHKERRQ(ierr);
1732   }
1733   PetscFunctionReturn(0);
1734 }
1735 
1736 
1737 #undef __FUNCT__
1738 #define __FUNCT__ "PCSetFromOptions_PFMG"
1739 PetscErrorCode PCSetFromOptions_PFMG(PetscOptions *PetscOptionsObject,PC pc)
1740 {
1741   PetscErrorCode ierr;
1742   PC_PFMG        *ex = (PC_PFMG*) pc->data;
1743   PetscBool      flg = PETSC_FALSE;
1744 
1745   PetscFunctionBegin;
1746   ierr = PetscOptionsHead(PetscOptionsObject,"PFMG options");CHKERRQ(ierr);
1747   ierr = PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);CHKERRQ(ierr);
1748   if (flg) {
1749     int level=3;
1750     PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level));
1751   }
1752   ierr = PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);CHKERRQ(ierr);
1753   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its));
1754   ierr = PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);CHKERRQ(ierr);
1755   PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax));
1756   ierr = PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);CHKERRQ(ierr);
1757   PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax));
1758 
1759   ierr = PetscOptionsInt("-pc_pfmg_max_levels","Max Levels for MG hierarchy","HYPRE_StructPFMGSetMaxLevels",ex->max_levels,&ex->max_levels,NULL);CHKERRQ(ierr);
1760   PetscStackCallStandard(HYPRE_StructPFMGSetMaxLevels,(ex->hsolver,ex->max_levels));
1761 
1762   ierr = PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);CHKERRQ(ierr);
1763   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol));
1764   ierr = PetscOptionsEList("-pc_pfmg_relax_type","Relax type for the up and down cycles","HYPRE_StructPFMGSetRelaxType",PFMGRelaxType,ALEN(PFMGRelaxType),PFMGRelaxType[ex->relax_type],&ex->relax_type,NULL);CHKERRQ(ierr);
1765   PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type));
1766   ierr = PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);CHKERRQ(ierr);
1767   PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type));
1768   ierr = PetscOptionsTail();CHKERRQ(ierr);
1769   PetscFunctionReturn(0);
1770 }
1771 
1772 #undef __FUNCT__
1773 #define __FUNCT__ "PCApply_PFMG"
1774 PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
1775 {
1776   PetscErrorCode    ierr;
1777   PC_PFMG           *ex = (PC_PFMG*) pc->data;
1778   PetscScalar       *yy;
1779   const PetscScalar *xx;
1780   PetscInt          ilower[3],iupper[3];
1781   Mat_HYPREStruct   *mx = (Mat_HYPREStruct*)(pc->pmat->data);
1782 
1783   PetscFunctionBegin;
1784   ierr       = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr);
1785   ierr       = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr);
1786   iupper[0] += ilower[0] - 1;
1787   iupper[1] += ilower[1] - 1;
1788   iupper[2] += ilower[2] - 1;
1789 
1790   /* copy x values over to hypre */
1791   PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
1792   ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
1793   PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx));
1794   ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
1795   PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
1796   PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
1797 
1798   /* copy solution values back to PETSc */
1799   ierr = VecGetArray(y,&yy);CHKERRQ(ierr);
1800   PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy));
1801   ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr);
1802   PetscFunctionReturn(0);
1803 }
1804 
1805 #undef __FUNCT__
1806 #define __FUNCT__ "PCApplyRichardson_PFMG"
1807 static PetscErrorCode PCApplyRichardson_PFMG(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
1808 {
1809   PC_PFMG        *jac = (PC_PFMG*)pc->data;
1810   PetscErrorCode ierr;
1811   PetscInt       oits;
1812 
1813   PetscFunctionBegin;
1814   ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr);
1815   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its));
1816   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol));
1817 
1818   ierr = PCApply_PFMG(pc,b,y);CHKERRQ(ierr);
1819   PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
1820   *outits = oits;
1821   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1822   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
1823   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol));
1824   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its));
1825   PetscFunctionReturn(0);
1826 }
1827 
1828 
1829 #undef __FUNCT__
1830 #define __FUNCT__ "PCSetUp_PFMG"
1831 PetscErrorCode PCSetUp_PFMG(PC pc)
1832 {
1833   PetscErrorCode  ierr;
1834   PC_PFMG         *ex = (PC_PFMG*) pc->data;
1835   Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
1836   PetscBool       flg;
1837 
1838   PetscFunctionBegin;
1839   ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESTRUCT,&flg);CHKERRQ(ierr);
1840   if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESTRUCT with this preconditioner");
1841 
1842   /* create the hypre solver object and set its information */
1843   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1844   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1845   PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
1846   PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver));
1847   PetscFunctionReturn(0);
1848 }
1849 
1850 
1851 /*MC
1852      PCPFMG - the hypre PFMG multigrid solver
1853 
1854    Level: advanced
1855 
1856    Options Database:
1857 + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner
1858 . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
1859 . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid
1860 . -pc_pfmg_tol <tol> tolerance of PFMG
1861 . -pc_pfmg_relax_type -relaxation type for the up and down cycles, one of Jacobi,Weighted-Jacobi,symmetric-Red/Black-Gauss-Seidel,Red/Black-Gauss-Seidel
1862 - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin
1863 
1864    Notes:  This is for CELL-centered descretizations
1865 
1866            This must be used with the MATHYPRESTRUCT matrix type.
1867            This is less general than in hypre, it supports only one block per process defined by a PETSc DMDA.
1868 
1869 .seealso:  PCMG, MATHYPRESTRUCT
1870 M*/
1871 
1872 #undef __FUNCT__
1873 #define __FUNCT__ "PCCreate_PFMG"
1874 PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
1875 {
1876   PetscErrorCode ierr;
1877   PC_PFMG        *ex;
1878 
1879   PetscFunctionBegin;
1880   ierr     = PetscNew(&ex);CHKERRQ(ierr); \
1881   pc->data = ex;
1882 
1883   ex->its            = 1;
1884   ex->tol            = 1.e-8;
1885   ex->relax_type     = 1;
1886   ex->rap_type       = 0;
1887   ex->num_pre_relax  = 1;
1888   ex->num_post_relax = 1;
1889   ex->max_levels     = 0;
1890 
1891   pc->ops->setfromoptions  = PCSetFromOptions_PFMG;
1892   pc->ops->view            = PCView_PFMG;
1893   pc->ops->destroy         = PCDestroy_PFMG;
1894   pc->ops->apply           = PCApply_PFMG;
1895   pc->ops->applyrichardson = PCApplyRichardson_PFMG;
1896   pc->ops->setup           = PCSetUp_PFMG;
1897 
1898   ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));CHKERRQ(ierr);
1899   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1900   PetscFunctionReturn(0);
1901 }
1902 
1903 /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/
1904 
1905 /* we know we are working with a HYPRE_SStructMatrix */
1906 typedef struct {
1907   MPI_Comm            hcomm;       /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
1908   HYPRE_SStructSolver ss_solver;
1909 
1910   /* keep copy of SYSPFMG options used so may view them */
1911   PetscInt its;
1912   double   tol;
1913   PetscInt relax_type;
1914   PetscInt num_pre_relax,num_post_relax;
1915 } PC_SysPFMG;
1916 
1917 #undef __FUNCT__
1918 #define __FUNCT__ "PCDestroy_SysPFMG"
1919 PetscErrorCode PCDestroy_SysPFMG(PC pc)
1920 {
1921   PetscErrorCode ierr;
1922   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;
1923 
1924   PetscFunctionBegin;
1925   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
1926   ierr = MPI_Comm_free(&ex->hcomm);CHKERRQ(ierr);
1927   ierr = PetscFree(pc->data);CHKERRQ(ierr);
1928   PetscFunctionReturn(0);
1929 }
1930 
1931 static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi","Red/Black-Gauss-Seidel"};
1932 
1933 #undef __FUNCT__
1934 #define __FUNCT__ "PCView_SysPFMG"
1935 PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
1936 {
1937   PetscErrorCode ierr;
1938   PetscBool      iascii;
1939   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;
1940 
1941   PetscFunctionBegin;
1942   ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
1943   if (iascii) {
1944     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG preconditioning\n");CHKERRQ(ierr);
1945     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: max iterations %d\n",ex->its);CHKERRQ(ierr);
1946     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: tolerance %g\n",ex->tol);CHKERRQ(ierr);
1947     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);CHKERRQ(ierr);
1948     ierr = PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);CHKERRQ(ierr);
1949   }
1950   PetscFunctionReturn(0);
1951 }
1952 
1953 
1954 #undef __FUNCT__
1955 #define __FUNCT__ "PCSetFromOptions_SysPFMG"
1956 PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptions *PetscOptionsObject,PC pc)
1957 {
1958   PetscErrorCode ierr;
1959   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;
1960   PetscBool      flg = PETSC_FALSE;
1961 
1962   PetscFunctionBegin;
1963   ierr = PetscOptionsHead(PetscOptionsObject,"SysPFMG options");CHKERRQ(ierr);
1964   ierr = PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);CHKERRQ(ierr);
1965   if (flg) {
1966     int level=3;
1967     PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level));
1968   }
1969   ierr = PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);CHKERRQ(ierr);
1970   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its));
1971   ierr = PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);CHKERRQ(ierr);
1972   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax));
1973   ierr = PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);CHKERRQ(ierr);
1974   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax));
1975 
1976   ierr = PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);CHKERRQ(ierr);
1977   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol));
1978   ierr = PetscOptionsEList("-pc_syspfmg_relax_type","Relax type for the up and down cycles","HYPRE_SStructSysPFMGSetRelaxType",SysPFMGRelaxType,4,SysPFMGRelaxType[ex->relax_type],&ex->relax_type,NULL);CHKERRQ(ierr);
1979   PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type));
1980   ierr = PetscOptionsTail();CHKERRQ(ierr);
1981   PetscFunctionReturn(0);
1982 }
1983 
1984 #undef __FUNCT__
1985 #define __FUNCT__ "PCApply_SysPFMG"
1986 PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
1987 {
1988   PetscErrorCode    ierr;
1989   PC_SysPFMG        *ex = (PC_SysPFMG*) pc->data;
1990   PetscScalar       *yy;
1991   const PetscScalar *xx;
1992   PetscInt          ilower[3],iupper[3];
1993   Mat_HYPRESStruct  *mx     = (Mat_HYPRESStruct*)(pc->pmat->data);
1994   PetscInt          ordering= mx->dofs_order;
1995   PetscInt          nvars   = mx->nvars;
1996   PetscInt          part    = 0;
1997   PetscInt          size;
1998   PetscInt          i;
1999 
2000   PetscFunctionBegin;
2001   ierr       = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr);
2002   ierr       = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr);
2003   iupper[0] += ilower[0] - 1;
2004   iupper[1] += ilower[1] - 1;
2005   iupper[2] += ilower[2] - 1;
2006 
2007   size = 1;
2008   for (i= 0; i< 3; i++) size *= (iupper[i]-ilower[i]+1);
2009 
2010   /* copy x values over to hypre for variable ordering */
2011   if (ordering) {
2012     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2013     ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
2014     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)xx+(size*i)));
2015     ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
2016     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2017     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
2018     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2019 
2020     /* copy solution values back to PETSc */
2021     ierr = VecGetArray(y,&yy);CHKERRQ(ierr);
2022     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i)));
2023     ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr);
2024   } else {      /* nodal ordering must be mapped to variable ordering for sys_pfmg */
2025     PetscScalar *z;
2026     PetscInt    j, k;
2027 
2028     ierr = PetscMalloc1(nvars*size,&z);CHKERRQ(ierr);
2029     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2030     ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
2031 
2032     /* transform nodal to hypre's variable ordering for sys_pfmg */
2033     for (i= 0; i< size; i++) {
2034       k= i*nvars;
2035       for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
2036     }
2037     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
2038     ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
2039     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2040     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2041 
2042     /* copy solution values back to PETSc */
2043     ierr = VecGetArray(y,&yy);CHKERRQ(ierr);
2044     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
2045     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
2046     for (i= 0; i< size; i++) {
2047       k= i*nvars;
2048       for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
2049     }
2050     ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr);
2051     ierr = PetscFree(z);CHKERRQ(ierr);
2052   }
2053   PetscFunctionReturn(0);
2054 }
2055 
2056 #undef __FUNCT__
2057 #define __FUNCT__ "PCApplyRichardson_SysPFMG"
2058 static PetscErrorCode PCApplyRichardson_SysPFMG(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
2059 {
2060   PC_SysPFMG     *jac = (PC_SysPFMG*)pc->data;
2061   PetscErrorCode ierr;
2062   PetscInt       oits;
2063 
2064   PetscFunctionBegin;
2065   ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr);
2066   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its));
2067   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol));
2068   ierr = PCApply_SysPFMG(pc,b,y);CHKERRQ(ierr);
2069   PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,(HYPRE_Int *)&oits));
2070   *outits = oits;
2071   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2072   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
2073   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol));
2074   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its));
2075   PetscFunctionReturn(0);
2076 }
2077 
2078 
2079 #undef __FUNCT__
2080 #define __FUNCT__ "PCSetUp_SysPFMG"
2081 PetscErrorCode PCSetUp_SysPFMG(PC pc)
2082 {
2083   PetscErrorCode   ierr;
2084   PC_SysPFMG       *ex = (PC_SysPFMG*) pc->data;
2085   Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
2086   PetscBool        flg;
2087 
2088   PetscFunctionBegin;
2089   ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESSTRUCT,&flg);CHKERRQ(ierr);
2090   if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESSTRUCT with this preconditioner");
2091 
2092   /* create the hypre sstruct solver object and set its information */
2093   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2094   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2095   PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver));
2096   PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2097   PetscFunctionReturn(0);
2098 }
2099 
2100 
2101 /*MC
2102      PCSysPFMG - the hypre SysPFMG multigrid solver
2103 
2104    Level: advanced
2105 
2106    Options Database:
2107 + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner
2108 . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2109 . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2110 . -pc_syspfmg_tol <tol> tolerance of SysPFMG
2111 . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel
2112 
2113    Notes:  This is for CELL-centered descretizations
2114 
2115            This must be used with the MATHYPRESSTRUCT matrix type.
2116            This is less general than in hypre, it supports only one part, and one block per process defined by a PETSc DMDA.
2117            Also, only cell-centered variables.
2118 
2119 .seealso:  PCMG, MATHYPRESSTRUCT
2120 M*/
2121 
2122 #undef __FUNCT__
2123 #define __FUNCT__ "PCCreate_SysPFMG"
2124 PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
2125 {
2126   PetscErrorCode ierr;
2127   PC_SysPFMG     *ex;
2128 
2129   PetscFunctionBegin;
2130   ierr     = PetscNew(&ex);CHKERRQ(ierr); \
2131   pc->data = ex;
2132 
2133   ex->its            = 1;
2134   ex->tol            = 1.e-8;
2135   ex->relax_type     = 1;
2136   ex->num_pre_relax  = 1;
2137   ex->num_post_relax = 1;
2138 
2139   pc->ops->setfromoptions  = PCSetFromOptions_SysPFMG;
2140   pc->ops->view            = PCView_SysPFMG;
2141   pc->ops->destroy         = PCDestroy_SysPFMG;
2142   pc->ops->apply           = PCApply_SysPFMG;
2143   pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
2144   pc->ops->setup           = PCSetUp_SysPFMG;
2145 
2146   ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));CHKERRQ(ierr);
2147   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2148   PetscFunctionReturn(0);
2149 }
2150