Actual source code: hypre.c

petsc-3.6.4 2016-04-12
Report Typos and Errors
  2: /*
  3:    Provides an interface to the LLNL package hypre
  4: */

  6: /* Must use hypre 2.0.0 or more recent. */

  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>

 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";

 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;

 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);

 31:   MPI_Comm comm_hypre;
 32:   char     *hypre_type;

 34:   /* options for Pilut and BoomerAMG*/
 35:   PetscInt maxiter;
 36:   double   tol;

 38:   /* options for Pilut */
 39:   PetscInt factorrowsize;

 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;

 51:   /* options for BoomerAMG */
 52:   PetscBool printstatistics;

 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;

 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;

 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;

103: PetscErrorCode PCHYPREGetSolver(PC pc,HYPRE_Solver *hsolver)
104: {
105:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

108:   *hsolver = jac->hsolver;
109:   return(0);
110: }

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;

123:   if (!jac->hypre_type) {
124:     PCHYPRESetType(pc,"boomeramg");
125:   }

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:   }

137:   if (!jac->ij) { /* create the matrix the first time through */
138:     MatHYPRE_IJMatrixCreate(pc->pmat,&jac->ij);
139:   }
140:   if (!jac->b) { /* create the vectors the first time through */
141:     Vec x,b;
142:     MatCreateVecs(pc->pmat,&x,&b);
143:     VecHYPRE_IJVectorCreate(x,&jac->x);
144:     VecHYPRE_IJVectorCreate(b,&jac->b);
145:     VecDestroy(&x);
146:     VecDestroy(&b);
147:   }

149:   /* special case for BoomerAMG */
150:   if (jac->setup == HYPRE_BoomerAMGSetup) {
151:     MatGetBlockSize(pc->pmat,&bs);
152:     if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
153:   }

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:   MatHYPRE_IJMatrixCopy(pc->pmat,jac->ij);
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",(*jac->setup)(jac->hsolver,hmat,bv,xv););
175:   return(0);
176: }

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: }

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;

204:   PetscCitationsRegister(hypreCitation,&cite);
205:   if (!jac->applyrichardson) {VecSet(x,0.0);}
206:   VecGetArrayRead(b,&bv);
207:   VecGetArray(x,&xv);
208:   HYPREReplacePointer(jac->b,(PetscScalar*)bv,sbv);
209:   HYPREReplacePointer(jac->x,xv,sxv);

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",h(*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;);

218:   HYPREReplacePointer(jac->b,(PetscScalar*)sbv,bv);
219:   HYPREReplacePointer(jac->x,sxv,xv);
220:   VecRestoreArray(x,&xv);
221:   VecRestoreArrayRead(b,&bv);
222:   return(0);
223: }

227: static PetscErrorCode PCReset_HYPRE(PC pc)
228: {
229:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;

233:   if (jac->ij) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij)); jac->ij = NULL;
234:   if (jac->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->b)); jac->b = NULL;
235:   if (jac->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->x)); jac->x = NULL;
236:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0])); jac->coords[0] = NULL;
237:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1])); jac->coords[1] = NULL;
238:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2])); jac->coords[2] = NULL;
239:   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0])); jac->constants[0] = NULL;
240:   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1])); jac->constants[1] = NULL;
241:   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2])); jac->constants[2] = NULL;
242:   if (jac->G) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->G)); jac->G = NULL;
243:   if (jac->C) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->C)); jac->C = NULL;
244:   if (jac->alpha_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->alpha_Poisson)); jac->alpha_Poisson = NULL;
245:   if (jac->beta_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->beta_Poisson)); jac->beta_Poisson = NULL;
246:   return(0);
247: }

251: static PetscErrorCode PCDestroy_HYPRE(PC pc)
252: {
253:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;

257:   PCReset_HYPRE(pc);
258:   if (jac->destroy) PetscStackCall("HYPRE_DestroyXXX",(*jac->destroy)(jac->hsolver););
259:   PetscFree(jac->hypre_type);
260:   if (jac->comm_hypre != MPI_COMM_NULL) { MPI_Comm_free(&(jac->comm_hypre));}
261:   PetscFree(pc->data);

263:   PetscObjectChangeTypeName((PetscObject)pc,0);
264:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);
265:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);
266:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetCoordinates_C",NULL);
267:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",NULL);
268:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",NULL);
269:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetConstantEdgeVectors_C",NULL);
270:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetAlphaPoissonMatrix_C",NULL);
271:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetBetaPoissonMatrix_C",NULL);
272:   return(0);
273: }

275: /* --------------------------------------------------------------------------------------------*/
278: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PetscOptions *PetscOptionsObject,PC pc)
279: {
280:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
282:   PetscBool      flag;

285:   PetscOptionsHead(PetscOptionsObject,"HYPRE Pilut Options");
286:   PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
287:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter));
288:   PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
289:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol));
290:   PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
291:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize));
292:   PetscOptionsTail();
293:   return(0);
294: }

298: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
299: {
300:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
302:   PetscBool      iascii;

305:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
306:   if (iascii) {
307:     PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut preconditioning\n");
308:     if (jac->maxiter != PETSC_DEFAULT) {
309:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: maximum number of iterations %d\n",jac->maxiter);
310:     } else {
311:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default maximum number of iterations \n");
312:     }
313:     if (jac->tol != PETSC_DEFAULT) {
314:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: drop tolerance %g\n",(double)jac->tol);
315:     } else {
316:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default drop tolerance \n");
317:     }
318:     if (jac->factorrowsize != PETSC_DEFAULT) {
319:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: factor row size %d\n",jac->factorrowsize);
320:     } else {
321:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default factor row size \n");
322:     }
323:   }
324:   return(0);
325: }

327: /* --------------------------------------------------------------------------------------------*/

331: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
332: {
333:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
334:   PetscErrorCode     ierr;
335:   HYPRE_ParCSRMatrix hmat;
336:   PetscScalar        *xv;
337:   const PetscScalar  *bv;
338:   HYPRE_ParVector    jbv,jxv;
339:   PetscScalar        *sbv,*sxv;
340:   PetscInt           hierr;

343:   PetscCitationsRegister(hypreCitation,&cite);
344:   VecSet(x,0.0);
345:   VecGetArrayRead(b,&bv);
346:   VecGetArray(x,&xv);
347:   HYPREReplacePointer(jac->b,(PetscScalar*)bv,sbv);
348:   HYPREReplacePointer(jac->x,xv,sxv);

350:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
351:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
352:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));

354:   hHYPRE_BoomerAMGSolveT(jac->hsolver,hmat,jbv,jxv);
355:   /* error code of 1 in BoomerAMG merely means convergence not achieved */
356:   if (hierr && (hierr != 1)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
357:   if (hierr) hypre__global_error = 0;

359:   HYPREReplacePointer(jac->b,sbv,bv);
360:   HYPREReplacePointer(jac->x,sxv,xv);
361:   VecRestoreArray(x,&xv);
362:   VecRestoreArrayRead(b,&bv);
363:   return(0);
364: }

366: /* static array length */
367: #define ALEN(a) (sizeof(a)/sizeof((a)[0]))

369: static const char *HYPREBoomerAMGCycleType[]   = {"","V","W"};
370: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
371: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
372: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
373: static const char *HYPREBoomerAMGRelaxType[]   = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi",
374:                                                   "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi",
375:                                                   "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination",
376:                                                   "" /* 10 */, "" /* 11 */, "" /* 12 */, "" /* 13 */, "" /* 14 */,
377:                                                   "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"};
378: static const char *HYPREBoomerAMGInterpType[]  = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
379:                                                   "ext+i-cc", "standard", "standard-wts", "", "", "FF", "FF1"};
382: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PetscOptions *PetscOptionsObject,PC pc)
383: {
384:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
386:   PetscInt       n,indx,level;
387:   PetscBool      flg, tmp_truth;
388:   double         tmpdbl, twodbl[2];

391:   PetscOptionsHead(PetscOptionsObject,"HYPRE BoomerAMG Options");
392:   PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);
393:   if (flg) {
394:     jac->cycletype = indx+1;
395:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
396:   }
397:   PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
398:   if (flg) {
399:     if (jac->maxlevels < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
400:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
401:   }
402:   PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
403:   if (flg) {
404:     if (jac->maxiter < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
405:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
406:   }
407:   PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);
408:   if (flg) {
409:     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);
410:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
411:   }

413:   PetscOptionsScalar("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);
414:   if (flg) {
415:     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);
416:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
417:   }

419:   PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);
420:   if (flg) {
421:     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);
422:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
423:   }

425:   PetscOptionsInt("-pc_hypre_boomeramg_agg_nl","Number of levels of aggressive coarsening","None",jac->agg_nl,&jac->agg_nl,&flg);
426:   if (flg) {
427:     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);

429:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
430:   }


433:   PetscOptionsInt("-pc_hypre_boomeramg_agg_num_paths","Number of paths for aggressive coarsening","None",jac->agg_num_paths,&jac->agg_num_paths,&flg);
434:   if (flg) {
435:     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);

437:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
438:   }


441:   PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
442:   if (flg) {
443:     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);
444:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
445:   }
446:   PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
447:   if (flg) {
448:     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);
449:     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);
450:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
451:   }

453:   /* Grid sweeps */
454:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);
455:   if (flg) {
456:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx));
457:     /* modify the jac structure so we can view the updated options with PC_View */
458:     jac->gridsweeps[0] = indx;
459:     jac->gridsweeps[1] = indx;
460:     /*defaults coarse to 1 */
461:     jac->gridsweeps[2] = 1;
462:   }

464:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);
465:   if (flg) {
466:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1));
467:     jac->gridsweeps[0] = indx;
468:   }
469:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);
470:   if (flg) {
471:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2));
472:     jac->gridsweeps[1] = indx;
473:   }
474:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);
475:   if (flg) {
476:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3));
477:     jac->gridsweeps[2] = indx;
478:   }

480:   /* Relax type */
481:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
482:   if (flg) {
483:     jac->relaxtype[0] = jac->relaxtype[1]  = indx;
484:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx));
485:     /* by default, coarse type set to 9 */
486:     jac->relaxtype[2] = 9;

488:   }
489:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
490:   if (flg) {
491:     jac->relaxtype[0] = indx;
492:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1));
493:   }
494:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
495:   if (flg) {
496:     jac->relaxtype[1] = indx;
497:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2));
498:   }
499:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);
500:   if (flg) {
501:     jac->relaxtype[2] = indx;
502:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3));
503:   }

505:   /* Relaxation Weight */
506:   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);
507:   if (flg) {
508:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl));
509:     jac->relaxweight = tmpdbl;
510:   }

512:   n         = 2;
513:   twodbl[0] = twodbl[1] = 1.0;
514:   PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
515:   if (flg) {
516:     if (n == 2) {
517:       indx =  (int)PetscAbsReal(twodbl[1]);
518:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx));
519:     } 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);
520:   }

522:   /* Outer relaxation Weight */
523:   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);
524:   if (flg) {
525:     PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl));
526:     jac->outerrelaxweight = tmpdbl;
527:   }

529:   n         = 2;
530:   twodbl[0] = twodbl[1] = 1.0;
531:   PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
532:   if (flg) {
533:     if (n == 2) {
534:       indx =  (int)PetscAbsReal(twodbl[1]);
535:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx));
536:     } 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);
537:   }

539:   /* the Relax Order */
540:   PetscOptionsBool("-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", PETSC_FALSE, &tmp_truth, &flg);

542:   if (flg && tmp_truth) {
543:     jac->relaxorder = 0;
544:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
545:   }
546:   PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);
547:   if (flg) {
548:     jac->measuretype = indx;
549:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
550:   }
551:   /* update list length 3/07 */
552:   PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);
553:   if (flg) {
554:     jac->coarsentype = indx;
555:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
556:   }

558:   /* new 3/07 */
559:   PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);
560:   if (flg) {
561:     jac->interptype = indx;
562:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
563:   }

565:   PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
566:   if (flg) {
567:     level = 3;
568:     PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);

570:     jac->printstatistics = PETSC_TRUE;
571:     PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level));
572:   }

574:   PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);
575:   if (flg) {
576:     level = 3;
577:     PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);

579:     jac->printstatistics = PETSC_TRUE;
580:     PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level));
581:   }

583:   PetscOptionsBool("-pc_hypre_boomeramg_nodal_coarsen", "HYPRE_BoomerAMGSetNodal()", "None", jac->nodal_coarsen ? PETSC_TRUE : PETSC_FALSE, &tmp_truth, &flg);
584:   if (flg && tmp_truth) {
585:     jac->nodal_coarsen = 1;
586:     PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,1));
587:   }

589:   PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);
590:   if (flg && tmp_truth) {
591:     PetscInt tmp_int;
592:     PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);
593:     if (flg) jac->nodal_relax_levels = tmp_int;
594:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6));
595:     PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1));
596:     PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0));
597:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels));
598:   }

600:   PetscOptionsTail();
601:   return(0);
602: }

606: 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)
607: {
608:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
610:   PetscInt       oits;

613:   PetscCitationsRegister(hypreCitation,&cite);
614:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter));
615:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol));
616:   jac->applyrichardson = PETSC_TRUE;
617:   PCApply_HYPRE(pc,b,y);
618:   jac->applyrichardson = PETSC_FALSE;
619:   PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
620:   *outits = oits;
621:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
622:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
623:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
624:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
625:   return(0);
626: }


631: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
632: {
633:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
635:   PetscBool      iascii;

638:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
639:   if (iascii) {
640:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG preconditioning\n");
641:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);
642:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of levels %d\n",jac->maxlevels);
643:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of iterations PER hypre call %d\n",jac->maxiter);
644:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Convergence tolerance PER hypre call %g\n",(double)jac->tol);
645:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Threshold for strong coupling %g\n",(double)jac->strongthreshold);
646:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation truncation factor %g\n",(double)jac->truncfactor);
647:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation: max elements per row %d\n",jac->pmax);
648:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Number of levels of aggressive coarsening %d\n",jac->agg_nl);
649:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Number of paths for aggressive coarsening %d\n",jac->agg_num_paths);

651:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum row sums %g\n",(double)jac->maxrowsum);

653:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps down         %d\n",jac->gridsweeps[0]);
654:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps up           %d\n",jac->gridsweeps[1]);
655:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps on coarse    %d\n",jac->gridsweeps[2]);

657:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax down          %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[0]]);
658:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax up            %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[1]]);
659:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax on coarse     %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[2]]);

661:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax weight  (all)      %g\n",(double)jac->relaxweight);
662:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Outer relax weight (all) %g\n",(double)jac->outerrelaxweight);

664:     if (jac->relaxorder) {
665:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Using CF-relaxation\n");
666:     } else {
667:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Not using CF-relaxation\n");
668:     }
669:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Measure type        %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
670:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Coarsen type        %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
671:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation type  %s\n",HYPREBoomerAMGInterpType[jac->interptype]);
672:     if (jac->nodal_coarsen) {
673:       PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal coarsening (with HYPRE_BOOMERAMGSetNodal())\n");
674:     }
675:     if (jac->nodal_relax) {
676:       PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal relaxation via Schwarz smoothing on levels %d\n",jac->nodal_relax_levels);
677:     }
678:   }
679:   return(0);
680: }

682: /* --------------------------------------------------------------------------------------------*/
685: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PetscOptions *PetscOptionsObject,PC pc)
686: {
687:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
689:   PetscInt       indx;
690:   PetscBool      flag;
691:   const char     *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};

694:   PetscOptionsHead(PetscOptionsObject,"HYPRE ParaSails Options");
695:   PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);
696:   PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshhold,&jac->threshhold,&flag);
697:   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));

699:   PetscOptionsReal("-pc_hypre_parasails_filter","filter","None",jac->filter,&jac->filter,&flag);
700:   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));

702:   PetscOptionsReal("-pc_hypre_parasails_loadbal","Load balance","None",jac->loadbal,&jac->loadbal,&flag);
703:   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));

705:   PetscOptionsBool("-pc_hypre_parasails_logging","Print info to screen","None",(PetscBool)jac->logging,(PetscBool*)&jac->logging,&flag);
706:   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));

708:   PetscOptionsBool("-pc_hypre_parasails_reuse","Reuse nonzero pattern in preconditioner","None",(PetscBool)jac->ruse,(PetscBool*)&jac->ruse,&flag);
709:   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));

711:   PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);
712:   if (flag) {
713:     jac->symt = indx;
714:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
715:   }

717:   PetscOptionsTail();
718:   return(0);
719: }

723: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
724: {
725:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
727:   PetscBool      iascii;
728:   const char     *symt = 0;;

731:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
732:   if (iascii) {
733:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails preconditioning\n");
734:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: nlevels %d\n",jac->nlevels);
735:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: threshold %g\n",(double)jac->threshhold);
736:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: filter %g\n",(double)jac->filter);
737:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: load balance %g\n",(double)jac->loadbal);
738:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: reuse nonzero structure %s\n",PetscBools[jac->ruse]);
739:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: print info to screen %s\n",PetscBools[jac->logging]);
740:     if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
741:     else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
742:     else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
743:     else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
744:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: %s\n",symt);
745:   }
746:   return(0);
747: }
748: /* --------------------------------------------------------------------------------------------*/
751: static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PetscOptions *PetscOptionsObject,PC pc)
752: {
753:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
755:   PetscInt       n;
756:   PetscBool      flag,flag2,flag3,flag4;

759:   PetscOptionsHead(PetscOptionsObject,"HYPRE AMS Options");
760:   PetscOptionsInt("-pc_hypre_ams_print_level","Debugging output level for AMS","None",jac->as_print,&jac->as_print,&flag);
761:   if (flag) PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
762:   PetscOptionsInt("-pc_hypre_ams_max_iter","Maximum number of AMS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
763:   if (flag) PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
764:   PetscOptionsInt("-pc_hypre_ams_cycle_type","Cycle type for AMS multigrid","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag);
765:   if (flag) PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
766:   PetscOptionsReal("-pc_hypre_ams_tol","Error tolerance for AMS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
767:   if (flag) PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
768:   PetscOptionsInt("-pc_hypre_ams_relax_type","Relaxation type for AMS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
769:   PetscOptionsInt("-pc_hypre_ams_relax_times","Number of relaxation steps for AMS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
770:   PetscOptionsReal("-pc_hypre_ams_relax_weight","Relaxation weight for AMS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
771:   PetscOptionsReal("-pc_hypre_ams_omega","SSOR coefficient for AMS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
772:   if (flag || flag2 || flag3 || flag4) {
773:     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
774:                                                                       jac->as_relax_times,
775:                                                                       jac->as_relax_weight,
776:                                                                       jac->as_omega));
777:   }
778:   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);
779:   n = 5;
780:   PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options","AMG options for vector Poisson","None",jac->as_amg_alpha_opts,&n,&flag2);
781:   if (flag || flag2) {
782:     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
783:                                                                      jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
784:                                                                      jac->as_amg_alpha_opts[2],       /* AMG relax_type */
785:                                                                      jac->as_amg_alpha_theta,
786:                                                                      jac->as_amg_alpha_opts[3],       /* AMG interp_type */
787:                                                                      jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
788:   }
789:   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);
790:   n = 5;
791:   PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options","AMG options for scalar Poisson solver","None",jac->as_amg_beta_opts,&n,&flag2);
792:   if (flag || flag2) {
793:     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
794:                                                                     jac->as_amg_beta_opts[1],       /* AMG agg_levels */
795:                                                                     jac->as_amg_beta_opts[2],       /* AMG relax_type */
796:                                                                     jac->as_amg_beta_theta,
797:                                                                     jac->as_amg_beta_opts[3],       /* AMG interp_type */
798:                                                                     jac->as_amg_beta_opts[4]));     /* AMG Pmax */
799:   }
800:   PetscOptionsTail();
801:   return(0);
802: }

806: static PetscErrorCode PCView_HYPRE_AMS(PC pc,PetscViewer viewer)
807: {
808:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
810:   PetscBool      iascii;

813:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
814:   if (iascii) {
815:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS preconditioning\n");
816:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace iterations per application %d\n",jac->as_max_iter);
817:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace cycle type %d\n",jac->ams_cycle_type);
818:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace iteration tolerance %g\n",jac->as_tol);
819:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother type %d\n",jac->as_relax_type);
820:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: number of smoothing steps %d\n",jac->as_relax_times);
821:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother weight %g\n",jac->as_relax_weight);
822:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother omega %g\n",jac->as_omega);
823:     if (jac->alpha_Poisson) {
824:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: vector Poisson solver (passed in by user)\n");
825:     } else {
826:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: vector Poisson solver (computed) \n");
827:     }
828:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);
829:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
830:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);
831:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);
832:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
833:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);
834:     if (!jac->ams_beta_is_zero) {
835:       if (jac->beta_Poisson) {
836:         PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: scalar Poisson solver (passed in by user)\n");
837:       } else {
838:         PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: scalar Poisson solver (computed) \n");
839:       }
840:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);
841:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
842:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);
843:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);
844:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
845:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);
846:     }
847:   }
848:   return(0);
849: }

853: static PetscErrorCode PCSetFromOptions_HYPRE_ADS(PetscOptions *PetscOptionsObject,PC pc)
854: {
855:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
857:   PetscInt       n;
858:   PetscBool      flag,flag2,flag3,flag4;

861:   PetscOptionsHead(PetscOptionsObject,"HYPRE ADS Options");
862:   PetscOptionsInt("-pc_hypre_ads_print_level","Debugging output level for ADS","None",jac->as_print,&jac->as_print,&flag);
863:   if (flag) PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
864:   PetscOptionsInt("-pc_hypre_ads_max_iter","Maximum number of ADS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
865:   if (flag) PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
866:   PetscOptionsInt("-pc_hypre_ads_cycle_type","Cycle type for ADS multigrid","None",jac->ads_cycle_type,&jac->ads_cycle_type,&flag);
867:   if (flag) PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ads_cycle_type));
868:   PetscOptionsReal("-pc_hypre_ads_tol","Error tolerance for ADS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
869:   if (flag) PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
870:   PetscOptionsInt("-pc_hypre_ads_relax_type","Relaxation type for ADS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
871:   PetscOptionsInt("-pc_hypre_ads_relax_times","Number of relaxation steps for ADS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
872:   PetscOptionsReal("-pc_hypre_ads_relax_weight","Relaxation weight for ADS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
873:   PetscOptionsReal("-pc_hypre_ads_omega","SSOR coefficient for ADS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
874:   if (flag || flag2 || flag3 || flag4) {
875:     PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
876:                                                                       jac->as_relax_times,
877:                                                                       jac->as_relax_weight,
878:                                                                       jac->as_omega));
879:   }
880:   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);
881:   n = 5;
882:   PetscOptionsIntArray("-pc_hypre_ads_ams_options","AMG options for AMS solver inside ADS","None",jac->as_amg_alpha_opts,&n,&flag2);
883:   PetscOptionsInt("-pc_hypre_ads_ams_cycle_type","Cycle type for AMS solver inside ADS","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag3);
884:   if (flag || flag2 || flag3) {
885:     PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type,             /* AMS cycle type */
886:                                                                 jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
887:                                                                 jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
888:                                                                 jac->as_amg_alpha_opts[2],       /* AMG relax_type */
889:                                                                 jac->as_amg_alpha_theta,
890:                                                                 jac->as_amg_alpha_opts[3],       /* AMG interp_type */
891:                                                                 jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
892:   }
893:   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);
894:   n = 5;
895:   PetscOptionsIntArray("-pc_hypre_ads_amg_options","AMG options for vector AMG solver inside ADS","None",jac->as_amg_beta_opts,&n,&flag2);
896:   if (flag || flag2) {
897:     PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
898:                                                                 jac->as_amg_beta_opts[1],       /* AMG agg_levels */
899:                                                                 jac->as_amg_beta_opts[2],       /* AMG relax_type */
900:                                                                 jac->as_amg_beta_theta,
901:                                                                 jac->as_amg_beta_opts[3],       /* AMG interp_type */
902:                                                                 jac->as_amg_beta_opts[4]));     /* AMG Pmax */
903:   }
904:   PetscOptionsTail();
905:   return(0);
906: }

910: static PetscErrorCode PCView_HYPRE_ADS(PC pc,PetscViewer viewer)
911: {
912:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
914:   PetscBool      iascii;

917:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
918:   if (iascii) {
919:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS preconditioning\n");
920:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: subspace iterations per application %d\n",jac->as_max_iter);
921:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: subspace cycle type %d\n",jac->ads_cycle_type);
922:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: subspace iteration tolerance %g\n",jac->as_tol);
923:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: smoother type %d\n",jac->as_relax_type);
924:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: number of smoothing steps %d\n",jac->as_relax_times);
925:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: smoother weight %g\n",jac->as_relax_weight);
926:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: smoother omega %g\n",jac->as_omega);
927:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: AMS solver\n");
928:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     subspace cycle type %d\n",jac->ams_cycle_type);
929:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);
930:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
931:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);
932:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);
933:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
934:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);
935:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: vector Poisson solver\n");
936:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);
937:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
938:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);
939:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);
940:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
941:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);
942:   }
943:   return(0);
944: }

948: static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE(PC pc, Mat G)
949: {
950:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
951:   HYPRE_ParCSRMatrix parcsr_G;
952:   PetscErrorCode     ierr;

955:   /* throw away any discrete gradient if already set */
956:   if (jac->G) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->G));
957:   MatHYPRE_IJMatrixCreate(G,&jac->G);
958:   MatHYPRE_IJMatrixCopy(G,jac->G);
959:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->G,(void**)(&parcsr_G)));
960:   PetscStackCall("Hypre set gradient",(*jac->setdgrad)(jac->hsolver,parcsr_G););
961:   return(0);
962: }

966: /*@
967:  PCHYPRESetDiscreteGradient - Set discrete gradient matrix

969:    Collective on PC

971:    Input Parameters:
972: +  pc - the preconditioning context
973: -  G - the discrete gradient

975:    Level: intermediate

977:    Notes: G should have as many rows as the number of edges and as many columns as the number of vertices in the mesh
978:           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

980: .seealso:
981: @*/
982: PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G)
983: {

990:   PetscTryMethod(pc,"PCHYPRESetDiscreteGradient_C",(PC,Mat),(pc,G));
991:   return(0);
992: }

996: static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C)
997: {
998:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
999:   HYPRE_ParCSRMatrix parcsr_C;
1000:   PetscErrorCode     ierr;

1003:   /* throw away any discrete curl if already set */
1004:   if (jac->C) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->C));
1005:   MatHYPRE_IJMatrixCreate(C,&jac->C);
1006:   MatHYPRE_IJMatrixCopy(C,jac->C);
1007:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->C,(void**)(&parcsr_C)));
1008:   PetscStackCall("Hypre set curl",(*jac->setdcurl)(jac->hsolver,parcsr_C););
1009:   return(0);
1010: }

1014: /*@
1015:  PCHYPRESetDiscreteCurl - Set discrete curl matrix

1017:    Collective on PC

1019:    Input Parameters:
1020: +  pc - the preconditioning context
1021: -  C - the discrete curl

1023:    Level: intermediate

1025:    Notes: C should have as many rows as the number of faces and as many columns as the number of edges in the mesh
1026:           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

1028: .seealso:
1029: @*/
1030: PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C)
1031: {

1038:   PetscTryMethod(pc,"PCHYPRESetDiscreteCurl_C",(PC,Mat),(pc,C));
1039:   return(0);
1040: }

1044: static PetscErrorCode PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS(PC pc, Mat A)
1045: {
1046:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1047:   HYPRE_ParCSRMatrix parcsr_alpha_Poisson;
1048:   PetscErrorCode     ierr;

1051:   /* throw away any matrix if already set */
1052:   if (jac->alpha_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->alpha_Poisson));
1053:   MatHYPRE_IJMatrixCreate(A,&jac->alpha_Poisson);
1054:   MatHYPRE_IJMatrixCopy(A,jac->alpha_Poisson);
1055:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->alpha_Poisson,(void**)(&parcsr_alpha_Poisson)));
1056:   PetscStackCallStandard(HYPRE_AMSSetAlphaPoissonMatrix,(jac->hsolver,parcsr_alpha_Poisson));
1057:   return(0);
1058: }

1062: /*@
1063:  PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix

1065:    Collective on PC

1067:    Input Parameters:
1068: +  pc - the preconditioning context
1069: -  A - the matrix

1071:    Level: intermediate

1073:    Notes: A should be obtained by discretizing the vector valued Poisson problem with linear finite elements

1075: .seealso:
1076: @*/
1077: PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A)
1078: {

1085:   PetscTryMethod(pc,"PCHYPRESetAlphaPoissonMatrix_C",(PC,Mat),(pc,A));
1086:   return(0);
1087: }

1091: static PetscErrorCode PCHYPRESetBetaPoissonMatrix_HYPRE_AMS(PC pc, Mat A)
1092: {
1093:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1094:   HYPRE_ParCSRMatrix parcsr_beta_Poisson;
1095:   PetscErrorCode     ierr;

1098:   if (!A) {
1099:     PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,NULL));
1100:     jac->ams_beta_is_zero = PETSC_TRUE;
1101:     return(0);
1102:   }
1103:   jac->ams_beta_is_zero = PETSC_FALSE;
1104:   /* throw away any matrix if already set */
1105:   if (jac->beta_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->beta_Poisson));
1106:   MatHYPRE_IJMatrixCreate(A,&jac->beta_Poisson);
1107:   MatHYPRE_IJMatrixCopy(A,jac->beta_Poisson);
1108:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->beta_Poisson,(void**)(&parcsr_beta_Poisson)));
1109:   PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,parcsr_beta_Poisson));
1110:   return(0);
1111: }

1115: /*@
1116:  PCHYPRESetBetaPoissonMatrix - Set Poisson matrix

1118:    Collective on PC

1120:    Input Parameters:
1121: +  pc - the preconditioning context
1122: -  A - the matrix

1124:    Level: intermediate

1126:    Notes: A should be obtained by discretizing the Poisson problem with linear finite elements.
1127:           Following HYPRE convention, the scalar Poisson solver of AMS can be turned off by passing NULL.

1129: .seealso:
1130: @*/
1131: PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A)
1132: {

1137:   if (A) {
1140:   }
1141:   PetscTryMethod(pc,"PCHYPRESetBetaPoissonMatrix_C",(PC,Mat),(pc,A));
1142:   return(0);
1143: }

1147: static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE_AMS(PC pc,Vec ozz, Vec zoz, Vec zzo)
1148: {
1149:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1150:   HYPRE_ParVector    par_ozz,par_zoz,par_zzo;
1151:   PetscInt           dim;
1152:   PetscErrorCode     ierr;

1155:   /* throw away any vector if already set */
1156:   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0]));
1157:   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1]));
1158:   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2]));
1159:   jac->constants[0] = NULL;
1160:   jac->constants[1] = NULL;
1161:   jac->constants[2] = NULL;
1162:   VecHYPRE_IJVectorCreate(ozz,&jac->constants[0]);
1163:   VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);
1164:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[0],(void**)(&par_ozz)));
1165:   VecHYPRE_IJVectorCreate(zoz,&jac->constants[1]);
1166:   VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);
1167:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[1],(void**)(&par_zoz)));
1168:   dim = 2;
1169:   if (zzo) {
1170:     VecHYPRE_IJVectorCreate(zzo,&jac->constants[2]);
1171:     VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);
1172:     PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[2],(void**)(&par_zzo)));
1173:     dim++;
1174:   }
1175:   PetscStackCallStandard(HYPRE_AMSSetEdgeConstantVectors,(jac->hsolver,par_ozz,par_zoz,par_zzo));
1176:   PetscStackCallStandard(HYPRE_AMSSetDimension,(jac->hsolver,dim));
1177:   return(0);
1178: }

1182: /*@
1183:  PCHYPRESetEdgeConstantVectors - Set the representation of the constant vector fields in edge element basis

1185:    Collective on PC

1187:    Input Parameters:
1188: +  pc - the preconditioning context
1189: -  ozz - vector representing (1,0,0) (or (1,0) in 2D)
1190: -  zoz - vector representing (0,1,0) (or (0,1) in 2D)
1191: -  zzo - vector representing (0,0,1) (use NULL in 2D)

1193:    Level: intermediate

1195:    Notes:

1197: .seealso:
1198: @*/
1199: PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo)
1200: {

1211:   PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));
1212:   return(0);
1213: }

1217: static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1218: {
1219:   PC_HYPRE        *jac = (PC_HYPRE*)pc->data;
1220:   Vec             tv;
1221:   HYPRE_ParVector par_coords[3];
1222:   PetscInt        i;
1223:   PetscErrorCode  ierr;

1226:   /* throw away any coordinate vector if already set */
1227:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0]));
1228:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1]));
1229:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2]));
1230:   /* set problem's dimension */
1231:   if (jac->setdim) {
1232:     PetscStackCall("Hypre set dim",(*jac->setdim)(jac->hsolver,dim););
1233:   }
1234:   /* compute IJ vector for coordinates */
1235:   VecCreate(PetscObjectComm((PetscObject)pc),&tv);
1236:   VecSetType(tv,VECSTANDARD);
1237:   VecSetSizes(tv,nloc,PETSC_DECIDE);
1238:   for (i=0;i<dim;i++) {
1239:     PetscScalar *array;
1240:     PetscInt    j;

1242:     VecHYPRE_IJVectorCreate(tv,&jac->coords[i]);
1243:     VecGetArray(tv,&array);
1244:     for (j=0;j<nloc;j++) {
1245:       array[j] = coords[j*dim+i];
1246:     }
1247:     PetscStackCallStandard(HYPRE_IJVectorSetValues,(jac->coords[i],nloc,NULL,array));
1248:     PetscStackCallStandard(HYPRE_IJVectorAssemble,(jac->coords[i]));
1249:     VecRestoreArray(tv,&array);
1250:   }
1251:   VecDestroy(&tv);
1252:   /* pass parCSR vectors to AMS solver */
1253:   par_coords[0] = NULL;
1254:   par_coords[1] = NULL;
1255:   par_coords[2] = NULL;
1256:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0],(void**)(&par_coords[0])));
1257:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1],(void**)(&par_coords[1])));
1258:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2],(void**)(&par_coords[2])));
1259:   PetscStackCall("Hypre set coords",(*jac->setcoord)(jac->hsolver,par_coords[0],par_coords[1],par_coords[2]););
1260:   return(0);
1261: }

1263: /* ---------------------------------------------------------------------------------*/

1267: static PetscErrorCode  PCHYPREGetType_HYPRE(PC pc,const char *name[])
1268: {
1269:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

1272:   *name = jac->hypre_type;
1273:   return(0);
1274: }

1278: static PetscErrorCode  PCHYPRESetType_HYPRE(PC pc,const char name[])
1279: {
1280:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1282:   PetscBool      flag;

1285:   if (jac->hypre_type) {
1286:     PetscStrcmp(jac->hypre_type,name,&flag);
1287:     if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set");
1288:     return(0);
1289:   } else {
1290:     PetscStrallocpy(name, &jac->hypre_type);
1291:   }

1293:   jac->maxiter         = PETSC_DEFAULT;
1294:   jac->tol             = PETSC_DEFAULT;
1295:   jac->printstatistics = PetscLogPrintInfo;

1297:   PetscStrcmp("pilut",jac->hypre_type,&flag);
1298:   if (flag) {
1299:     PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver));
1300:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
1301:     pc->ops->view           = PCView_HYPRE_Pilut;
1302:     jac->destroy            = HYPRE_ParCSRPilutDestroy;
1303:     jac->setup              = HYPRE_ParCSRPilutSetup;
1304:     jac->solve              = HYPRE_ParCSRPilutSolve;
1305:     jac->factorrowsize      = PETSC_DEFAULT;
1306:     return(0);
1307:   }
1308:   PetscStrcmp("parasails",jac->hypre_type,&flag);
1309:   if (flag) {
1310:     PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver));
1311:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
1312:     pc->ops->view           = PCView_HYPRE_ParaSails;
1313:     jac->destroy            = HYPRE_ParaSailsDestroy;
1314:     jac->setup              = HYPRE_ParaSailsSetup;
1315:     jac->solve              = HYPRE_ParaSailsSolve;
1316:     /* initialize */
1317:     jac->nlevels    = 1;
1318:     jac->threshhold = .1;
1319:     jac->filter     = .1;
1320:     jac->loadbal    = 0;
1321:     if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE;
1322:     else jac->logging = (int) PETSC_FALSE;

1324:     jac->ruse = (int) PETSC_FALSE;
1325:     jac->symt = 0;
1326:     PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
1327:     PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
1328:     PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
1329:     PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
1330:     PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
1331:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
1332:     return(0);
1333:   }
1334:   PetscStrcmp("boomeramg",jac->hypre_type,&flag);
1335:   if (flag) {
1336:     HYPRE_BoomerAMGCreate(&jac->hsolver);
1337:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_BoomerAMG;
1338:     pc->ops->view            = PCView_HYPRE_BoomerAMG;
1339:     pc->ops->applytranspose  = PCApplyTranspose_HYPRE_BoomerAMG;
1340:     pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
1341:     jac->destroy             = HYPRE_BoomerAMGDestroy;
1342:     jac->setup               = HYPRE_BoomerAMGSetup;
1343:     jac->solve               = HYPRE_BoomerAMGSolve;
1344:     jac->applyrichardson     = PETSC_FALSE;
1345:     /* these defaults match the hypre defaults */
1346:     jac->cycletype        = 1;
1347:     jac->maxlevels        = 25;
1348:     jac->maxiter          = 1;
1349:     jac->tol              = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
1350:     jac->truncfactor      = 0.0;
1351:     jac->strongthreshold  = .25;
1352:     jac->maxrowsum        = .9;
1353:     jac->coarsentype      = 6;
1354:     jac->measuretype      = 0;
1355:     jac->gridsweeps[0]    = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
1356:     jac->relaxtype[0]     = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */
1357:     jac->relaxtype[2]     = 9; /*G.E. */
1358:     jac->relaxweight      = 1.0;
1359:     jac->outerrelaxweight = 1.0;
1360:     jac->relaxorder       = 1;
1361:     jac->interptype       = 0;
1362:     jac->agg_nl           = 0;
1363:     jac->pmax             = 0;
1364:     jac->truncfactor      = 0.0;
1365:     jac->agg_num_paths    = 1;

1367:     jac->nodal_coarsen      = 0;
1368:     jac->nodal_relax        = PETSC_FALSE;
1369:     jac->nodal_relax_levels = 1;
1370:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
1371:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
1372:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
1373:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
1374:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
1375:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
1376:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
1377:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
1378:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
1379:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
1380:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
1381:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
1382:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
1383:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
1384:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0]));  /*defaults coarse to 9*/
1385:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */
1386:     return(0);
1387:   }
1388:   PetscStrcmp("ams",jac->hypre_type,&flag);
1389:   if (flag) {
1390:     HYPRE_AMSCreate(&jac->hsolver);
1391:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_AMS;
1392:     pc->ops->view            = PCView_HYPRE_AMS;
1393:     jac->destroy             = HYPRE_AMSDestroy;
1394:     jac->setup               = HYPRE_AMSSetup;
1395:     jac->solve               = HYPRE_AMSSolve;
1396:     jac->setdgrad            = HYPRE_AMSSetDiscreteGradient;
1397:     jac->setcoord            = HYPRE_AMSSetCoordinateVectors;
1398:     jac->setdim              = HYPRE_AMSSetDimension;
1399:     jac->coords[0]           = NULL;
1400:     jac->coords[1]           = NULL;
1401:     jac->coords[2]           = NULL;
1402:     jac->G                   = NULL;
1403:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */
1404:     jac->as_print           = 0;
1405:     jac->as_max_iter        = 1; /* used as a preconditioner */
1406:     jac->as_tol             = 0.; /* used as a preconditioner */
1407:     jac->ams_cycle_type     = 13;
1408:     /* Smoothing options */
1409:     jac->as_relax_type      = 2;
1410:     jac->as_relax_times     = 1;
1411:     jac->as_relax_weight    = 1.0;
1412:     jac->as_omega           = 1.0;
1413:     /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1414:     jac->as_amg_alpha_opts[0] = 10;
1415:     jac->as_amg_alpha_opts[1] = 1;
1416:     jac->as_amg_alpha_opts[2] = 6;
1417:     jac->as_amg_alpha_opts[3] = 6;
1418:     jac->as_amg_alpha_opts[4] = 4;
1419:     jac->as_amg_alpha_theta   = 0.25;
1420:     /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1421:     jac->ams_beta_is_zero = PETSC_FALSE;
1422:     jac->as_amg_beta_opts[0] = 10;
1423:     jac->as_amg_beta_opts[1] = 1;
1424:     jac->as_amg_beta_opts[2] = 6;
1425:     jac->as_amg_beta_opts[3] = 6;
1426:     jac->as_amg_beta_opts[4] = 4;
1427:     jac->as_amg_beta_theta   = 0.25;
1428:     PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
1429:     PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1430:     PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1431:     PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
1432:     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1433:                                                                       jac->as_relax_times,
1434:                                                                       jac->as_relax_weight,
1435:                                                                       jac->as_omega));
1436:     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1437:                                                                      jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1438:                                                                      jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1439:                                                                      jac->as_amg_alpha_theta,
1440:                                                                      jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1441:                                                                      jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1442:     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1443:                                                                     jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1444:                                                                     jac->as_amg_beta_opts[2],       /* AMG relax_type */
1445:                                                                     jac->as_amg_beta_theta,
1446:                                                                     jac->as_amg_beta_opts[3],       /* AMG interp_type */
1447:                                                                     jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1448:     PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);
1449:     PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);
1450:     PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE_AMS);
1451:     PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetAlphaPoissonMatrix_C",PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS);
1452:     PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetBetaPoissonMatrix_C",PCHYPRESetBetaPoissonMatrix_HYPRE_AMS);
1453:     return(0);
1454:   }
1455:   PetscStrcmp("ads",jac->hypre_type,&flag);
1456:   if (flag) {
1457:     HYPRE_ADSCreate(&jac->hsolver);
1458:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_ADS;
1459:     pc->ops->view            = PCView_HYPRE_ADS;
1460:     jac->destroy             = HYPRE_ADSDestroy;
1461:     jac->setup               = HYPRE_ADSSetup;
1462:     jac->solve               = HYPRE_ADSSolve;
1463:     jac->setdgrad            = HYPRE_ADSSetDiscreteGradient;
1464:     jac->setdcurl            = HYPRE_ADSSetDiscreteCurl;
1465:     jac->setcoord            = HYPRE_ADSSetCoordinateVectors;
1466:     jac->coords[0]           = NULL;
1467:     jac->coords[1]           = NULL;
1468:     jac->coords[2]           = NULL;
1469:     jac->G                   = NULL;
1470:     jac->C                   = NULL;
1471:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */
1472:     jac->as_print           = 0;
1473:     jac->as_max_iter        = 1; /* used as a preconditioner */
1474:     jac->as_tol             = 0.; /* used as a preconditioner */
1475:     jac->ads_cycle_type     = 13;
1476:     /* Smoothing options */
1477:     jac->as_relax_type      = 2;
1478:     jac->as_relax_times     = 1;
1479:     jac->as_relax_weight    = 1.0;
1480:     jac->as_omega           = 1.0;
1481:     /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */
1482:     jac->ams_cycle_type       = 14;
1483:     jac->as_amg_alpha_opts[0] = 10;
1484:     jac->as_amg_alpha_opts[1] = 1;
1485:     jac->as_amg_alpha_opts[2] = 6;
1486:     jac->as_amg_alpha_opts[3] = 6;
1487:     jac->as_amg_alpha_opts[4] = 4;
1488:     jac->as_amg_alpha_theta   = 0.25;
1489:     /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1490:     jac->as_amg_beta_opts[0] = 10;
1491:     jac->as_amg_beta_opts[1] = 1;
1492:     jac->as_amg_beta_opts[2] = 6;
1493:     jac->as_amg_beta_opts[3] = 6;
1494:     jac->as_amg_beta_opts[4] = 4;
1495:     jac->as_amg_beta_theta   = 0.25;
1496:     PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
1497:     PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1498:     PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1499:     PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
1500:     PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1501:                                                                       jac->as_relax_times,
1502:                                                                       jac->as_relax_weight,
1503:                                                                       jac->as_omega));
1504:     PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type,             /* AMG coarsen type */
1505:                                                                 jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1506:                                                                 jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1507:                                                                 jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1508:                                                                 jac->as_amg_alpha_theta,
1509:                                                                 jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1510:                                                                 jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1511:     PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1512:                                                                 jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1513:                                                                 jac->as_amg_beta_opts[2],       /* AMG relax_type */
1514:                                                                 jac->as_amg_beta_theta,
1515:                                                                 jac->as_amg_beta_opts[3],       /* AMG interp_type */
1516:                                                                 jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1517:     PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);
1518:     PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);
1519:     PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",PCHYPRESetDiscreteCurl_HYPRE);
1520:     return(0);
1521:   }
1522:   PetscFree(jac->hypre_type);

1524:   jac->hypre_type = NULL;
1525:   SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are pilut, parasails, boomeramg, ams",name);
1526:   return(0);
1527: }

1529: /*
1530:     It only gets here if the HYPRE type has not been set before the call to
1531:    ...SetFromOptions() which actually is most of the time
1532: */
1535: static PetscErrorCode PCSetFromOptions_HYPRE(PetscOptions *PetscOptionsObject,PC pc)
1536: {
1538:   PetscInt       indx;
1539:   const char     *type[] = {"pilut","parasails","boomeramg","ams","ads"};
1540:   PetscBool      flg;

1543:   PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");
1544:   PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"boomeramg",&indx,&flg);
1545:   if (flg) {
1546:     PCHYPRESetType_HYPRE(pc,type[indx]);
1547:   } else {
1548:     PCHYPRESetType_HYPRE(pc,"boomeramg");
1549:   }
1550:   if (pc->ops->setfromoptions) {
1551:     pc->ops->setfromoptions(PetscOptionsObject,pc);
1552:   }
1553:   PetscOptionsTail();
1554:   return(0);
1555: }

1559: /*@C
1560:      PCHYPRESetType - Sets which hypre preconditioner you wish to use

1562:    Input Parameters:
1563: +     pc - the preconditioner context
1564: -     name - either  pilut, parasails, boomeramg, ams, ads

1566:    Options Database Keys:
1567:    -pc_hypre_type - One of pilut, parasails, boomeramg, ams, ads

1569:    Level: intermediate

1571: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
1572:            PCHYPRE

1574: @*/
1575: PetscErrorCode  PCHYPRESetType(PC pc,const char name[])
1576: {

1582:   PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));
1583:   return(0);
1584: }

1588: /*@C
1589:      PCHYPREGetType - Gets which hypre preconditioner you are using

1591:    Input Parameter:
1592: .     pc - the preconditioner context

1594:    Output Parameter:
1595: .     name - either  pilut, parasails, boomeramg, ams, ads

1597:    Level: intermediate

1599: .seealso:  PCCreate(), PCHYPRESetType(), PCType (for list of available types), PC,
1600:            PCHYPRE

1602: @*/
1603: PetscErrorCode  PCHYPREGetType(PC pc,const char *name[])
1604: {

1610:   PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));
1611:   return(0);
1612: }

1614: /*MC
1615:      PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre

1617:    Options Database Keys:
1618: +   -pc_hypre_type - One of pilut, parasails, boomeramg, ams, ads
1619: -   Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX
1620:           preconditioner

1622:    Level: intermediate

1624:    Notes: Apart from pc_hypre_type (for which there is PCHYPRESetType()),
1625:           the many hypre options can ONLY be set via the options database (e.g. the command line
1626:           or with PetscOptionsSetValue(), there are no functions to set them)

1628:           The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_rtol refer to the number of iterations
1629:           (V-cycles) and tolerance that boomeramg does EACH time it is called. So for example, if
1630:           -pc_hypre_boomeramg_max_iter is set to 2 then 2-V-cycles are being used to define the preconditioner
1631:           (-pc_hypre_boomeramg_rtol should be set to 0.0 - the default - to strictly use a fixed number of
1632:           iterations per hypre call). -ksp_max_it and -ksp_rtol STILL determine the total number of iterations
1633:           and tolerance for the Krylov solver. For example, if -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10
1634:           then AT MOST twenty V-cycles of boomeramg will be called.

1636:            Note that the option -pc_hypre_boomeramg_relax_type_all defaults to symmetric relaxation
1637:            (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry.
1638:            Otherwise, you may want to use -pc_hypre_boomeramg_relax_type_all SOR/Jacobi.
1639:           If you wish to use BoomerAMG WITHOUT a Krylov method use -ksp_type richardson NOT -ksp_type preonly
1640:           and use -ksp_max_it to control the number of V-cycles.
1641:           (see the PETSc FAQ.html at the PETSc website under the Documentation tab).

1643:           2007-02-03 Using HYPRE-1.11.1b, the routine HYPRE_BoomerAMGSolveT and the option
1644:           -pc_hypre_parasails_reuse were failing with SIGSEGV. Dalcin L.

1646:           See PCPFMG for access to the hypre Struct PFMG solver

1648: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
1649:            PCHYPRESetType(), PCPFMG

1651: M*/

1655: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
1656: {
1657:   PC_HYPRE       *jac;

1661:   PetscNewLog(pc,&jac);

1663:   pc->data                = jac;
1664:   pc->ops->reset          = PCReset_HYPRE;
1665:   pc->ops->destroy        = PCDestroy_HYPRE;
1666:   pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
1667:   pc->ops->setup          = PCSetUp_HYPRE;
1668:   pc->ops->apply          = PCApply_HYPRE;
1669:   jac->comm_hypre         = MPI_COMM_NULL;
1670:   jac->hypre_type         = NULL;
1671:   jac->coords[0]          = NULL;
1672:   jac->coords[1]          = NULL;
1673:   jac->coords[2]          = NULL;
1674:   jac->constants[0]       = NULL;
1675:   jac->constants[1]       = NULL;
1676:   jac->constants[2]       = NULL;
1677:   jac->G                  = NULL;
1678:   jac->C                  = NULL;
1679:   jac->alpha_Poisson      = NULL;
1680:   jac->beta_Poisson       = NULL;
1681:   jac->setdim             = NULL;
1682:   /* duplicate communicator for hypre */
1683:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1684:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);
1685:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);
1686:   return(0);
1687: }

1689: /* ---------------------------------------------------------------------------------------------------------------------------------*/

1691: /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */
1692: #include <petsc/private/matimpl.h>

1694: typedef struct {
1695:   MPI_Comm           hcomm;        /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */
1696:   HYPRE_StructSolver hsolver;

1698:   /* keep copy of PFMG options used so may view them */
1699:   PetscInt its;
1700:   double   tol;
1701:   PetscInt relax_type;
1702:   PetscInt rap_type;
1703:   PetscInt num_pre_relax,num_post_relax;
1704:   PetscInt max_levels;
1705: } PC_PFMG;

1709: PetscErrorCode PCDestroy_PFMG(PC pc)
1710: {
1712:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

1715:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1716:   MPI_Comm_free(&ex->hcomm);
1717:   PetscFree(pc->data);
1718:   return(0);
1719: }

1721: static const char *PFMGRelaxType[] = {"Jacobi","Weighted-Jacobi","symmetric-Red/Black-Gauss-Seidel","Red/Black-Gauss-Seidel"};
1722: static const char *PFMGRAPType[] = {"Galerkin","non-Galerkin"};

1726: PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
1727: {
1729:   PetscBool      iascii;
1730:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

1733:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1734:   if (iascii) {
1735:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG preconditioning\n");
1736:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: max iterations %d\n",ex->its);
1737:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: tolerance %g\n",ex->tol);
1738:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
1739:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: RAP type %s\n",PFMGRAPType[ex->rap_type]);
1740:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
1741:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: max levels %d\n",ex->max_levels);
1742:   }
1743:   return(0);
1744: }


1749: PetscErrorCode PCSetFromOptions_PFMG(PetscOptions *PetscOptionsObject,PC pc)
1750: {
1752:   PC_PFMG        *ex = (PC_PFMG*) pc->data;
1753:   PetscBool      flg = PETSC_FALSE;

1756:   PetscOptionsHead(PetscOptionsObject,"PFMG options");
1757:   PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);
1758:   if (flg) {
1759:     int level=3;
1760:     PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level));
1761:   }
1762:   PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);
1763:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its));
1764:   PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
1765:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax));
1766:   PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
1767:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax));

1769:   PetscOptionsInt("-pc_pfmg_max_levels","Max Levels for MG hierarchy","HYPRE_StructPFMGSetMaxLevels",ex->max_levels,&ex->max_levels,NULL);
1770:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxLevels,(ex->hsolver,ex->max_levels));

1772:   PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);
1773:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol));
1774:   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);
1775:   PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type));
1776:   PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);
1777:   PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type));
1778:   PetscOptionsTail();
1779:   return(0);
1780: }

1784: PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
1785: {
1786:   PetscErrorCode    ierr;
1787:   PC_PFMG           *ex = (PC_PFMG*) pc->data;
1788:   PetscScalar       *yy;
1789:   const PetscScalar *xx;
1790:   PetscInt          ilower[3],iupper[3];
1791:   Mat_HYPREStruct   *mx = (Mat_HYPREStruct*)(pc->pmat->data);

1794:   PetscCitationsRegister(hypreCitation,&cite);
1795:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
1796:   iupper[0] += ilower[0] - 1;
1797:   iupper[1] += ilower[1] - 1;
1798:   iupper[2] += ilower[2] - 1;

1800:   /* copy x values over to hypre */
1801:   PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
1802:   VecGetArrayRead(x,&xx);
1803:   PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx));
1804:   VecRestoreArrayRead(x,&xx);
1805:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
1806:   PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx));

1808:   /* copy solution values back to PETSc */
1809:   VecGetArray(y,&yy);
1810:   PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy));
1811:   VecRestoreArray(y,&yy);
1812:   return(0);
1813: }

1817: 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)
1818: {
1819:   PC_PFMG        *jac = (PC_PFMG*)pc->data;
1821:   PetscInt       oits;

1824:   PetscCitationsRegister(hypreCitation,&cite);
1825:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its));
1826:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol));

1828:   PCApply_PFMG(pc,b,y);
1829:   PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
1830:   *outits = oits;
1831:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1832:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
1833:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol));
1834:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its));
1835:   return(0);
1836: }


1841: PetscErrorCode PCSetUp_PFMG(PC pc)
1842: {
1843:   PetscErrorCode  ierr;
1844:   PC_PFMG         *ex = (PC_PFMG*) pc->data;
1845:   Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
1846:   PetscBool       flg;

1849:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESTRUCT,&flg);
1850:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESTRUCT with this preconditioner");

1852:   /* create the hypre solver object and set its information */
1853:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1854:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1855:   PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
1856:   PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver));
1857:   return(0);
1858: }


1861: /*MC
1862:      PCPFMG - the hypre PFMG multigrid solver

1864:    Level: advanced

1866:    Options Database:
1867: + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner
1868: . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
1869: . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid
1870: . -pc_pfmg_tol <tol> tolerance of PFMG
1871: . -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
1872: - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin

1874:    Notes:  This is for CELL-centered descretizations

1876:            This must be used with the MATHYPRESTRUCT matrix type.
1877:            This is less general than in hypre, it supports only one block per process defined by a PETSc DMDA.

1879: .seealso:  PCMG, MATHYPRESTRUCT
1880: M*/

1884: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
1885: {
1887:   PC_PFMG        *ex;

1890:   PetscNew(&ex); \
1891:   pc->data = ex;

1893:   ex->its            = 1;
1894:   ex->tol            = 1.e-8;
1895:   ex->relax_type     = 1;
1896:   ex->rap_type       = 0;
1897:   ex->num_pre_relax  = 1;
1898:   ex->num_post_relax = 1;
1899:   ex->max_levels     = 0;

1901:   pc->ops->setfromoptions  = PCSetFromOptions_PFMG;
1902:   pc->ops->view            = PCView_PFMG;
1903:   pc->ops->destroy         = PCDestroy_PFMG;
1904:   pc->ops->apply           = PCApply_PFMG;
1905:   pc->ops->applyrichardson = PCApplyRichardson_PFMG;
1906:   pc->ops->setup           = PCSetUp_PFMG;

1908:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
1909:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1910:   return(0);
1911: }

1913: /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/

1915: /* we know we are working with a HYPRE_SStructMatrix */
1916: typedef struct {
1917:   MPI_Comm            hcomm;       /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
1918:   HYPRE_SStructSolver ss_solver;

1920:   /* keep copy of SYSPFMG options used so may view them */
1921:   PetscInt its;
1922:   double   tol;
1923:   PetscInt relax_type;
1924:   PetscInt num_pre_relax,num_post_relax;
1925: } PC_SysPFMG;

1929: PetscErrorCode PCDestroy_SysPFMG(PC pc)
1930: {
1932:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

1935:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
1936:   MPI_Comm_free(&ex->hcomm);
1937:   PetscFree(pc->data);
1938:   return(0);
1939: }

1941: static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi","Red/Black-Gauss-Seidel"};

1945: PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
1946: {
1948:   PetscBool      iascii;
1949:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

1952:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1953:   if (iascii) {
1954:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG preconditioning\n");
1955:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: max iterations %d\n",ex->its);
1956:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: tolerance %g\n",ex->tol);
1957:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
1958:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
1959:   }
1960:   return(0);
1961: }


1966: PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptions *PetscOptionsObject,PC pc)
1967: {
1969:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;
1970:   PetscBool      flg = PETSC_FALSE;

1973:   PetscOptionsHead(PetscOptionsObject,"SysPFMG options");
1974:   PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);
1975:   if (flg) {
1976:     int level=3;
1977:     PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level));
1978:   }
1979:   PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);
1980:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its));
1981:   PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
1982:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax));
1983:   PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
1984:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax));

1986:   PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);
1987:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol));
1988:   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);
1989:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type));
1990:   PetscOptionsTail();
1991:   return(0);
1992: }

1996: PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
1997: {
1998:   PetscErrorCode    ierr;
1999:   PC_SysPFMG        *ex = (PC_SysPFMG*) pc->data;
2000:   PetscScalar       *yy;
2001:   const PetscScalar *xx;
2002:   PetscInt          ilower[3],iupper[3];
2003:   Mat_HYPRESStruct  *mx     = (Mat_HYPRESStruct*)(pc->pmat->data);
2004:   PetscInt          ordering= mx->dofs_order;
2005:   PetscInt          nvars   = mx->nvars;
2006:   PetscInt          part    = 0;
2007:   PetscInt          size;
2008:   PetscInt          i;

2011:   PetscCitationsRegister(hypreCitation,&cite);
2012:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
2013:   iupper[0] += ilower[0] - 1;
2014:   iupper[1] += ilower[1] - 1;
2015:   iupper[2] += ilower[2] - 1;

2017:   size = 1;
2018:   for (i= 0; i< 3; i++) size *= (iupper[i]-ilower[i]+1);

2020:   /* copy x values over to hypre for variable ordering */
2021:   if (ordering) {
2022:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2023:     VecGetArrayRead(x,&xx);
2024:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)xx+(size*i)));
2025:     VecRestoreArrayRead(x,&xx);
2026:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2027:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
2028:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

2030:     /* copy solution values back to PETSc */
2031:     VecGetArray(y,&yy);
2032:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i)));
2033:     VecRestoreArray(y,&yy);
2034:   } else {      /* nodal ordering must be mapped to variable ordering for sys_pfmg */
2035:     PetscScalar *z;
2036:     PetscInt    j, k;

2038:     PetscMalloc1(nvars*size,&z);
2039:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2040:     VecGetArrayRead(x,&xx);

2042:     /* transform nodal to hypre's variable ordering for sys_pfmg */
2043:     for (i= 0; i< size; i++) {
2044:       k= i*nvars;
2045:       for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
2046:     }
2047:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
2048:     VecRestoreArrayRead(x,&xx);
2049:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2050:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

2052:     /* copy solution values back to PETSc */
2053:     VecGetArray(y,&yy);
2054:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
2055:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
2056:     for (i= 0; i< size; i++) {
2057:       k= i*nvars;
2058:       for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
2059:     }
2060:     VecRestoreArray(y,&yy);
2061:     PetscFree(z);
2062:   }
2063:   return(0);
2064: }

2068: 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)
2069: {
2070:   PC_SysPFMG     *jac = (PC_SysPFMG*)pc->data;
2072:   PetscInt       oits;

2075:   PetscCitationsRegister(hypreCitation,&cite);
2076:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its));
2077:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol));
2078:   PCApply_SysPFMG(pc,b,y);
2079:   PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,(HYPRE_Int *)&oits));
2080:   *outits = oits;
2081:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2082:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
2083:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol));
2084:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its));
2085:   return(0);
2086: }


2091: PetscErrorCode PCSetUp_SysPFMG(PC pc)
2092: {
2093:   PetscErrorCode   ierr;
2094:   PC_SysPFMG       *ex = (PC_SysPFMG*) pc->data;
2095:   Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
2096:   PetscBool        flg;

2099:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESSTRUCT,&flg);
2100:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESSTRUCT with this preconditioner");

2102:   /* create the hypre sstruct solver object and set its information */
2103:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2104:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2105:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver));
2106:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2107:   return(0);
2108: }


2111: /*MC
2112:      PCSysPFMG - the hypre SysPFMG multigrid solver

2114:    Level: advanced

2116:    Options Database:
2117: + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner
2118: . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2119: . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2120: . -pc_syspfmg_tol <tol> tolerance of SysPFMG
2121: . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel

2123:    Notes:  This is for CELL-centered descretizations

2125:            This must be used with the MATHYPRESSTRUCT matrix type.
2126:            This is less general than in hypre, it supports only one part, and one block per process defined by a PETSc DMDA.
2127:            Also, only cell-centered variables.

2129: .seealso:  PCMG, MATHYPRESSTRUCT
2130: M*/

2134: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
2135: {
2137:   PC_SysPFMG     *ex;

2140:   PetscNew(&ex); \
2141:   pc->data = ex;

2143:   ex->its            = 1;
2144:   ex->tol            = 1.e-8;
2145:   ex->relax_type     = 1;
2146:   ex->num_pre_relax  = 1;
2147:   ex->num_post_relax = 1;

2149:   pc->ops->setfromoptions  = PCSetFromOptions_SysPFMG;
2150:   pc->ops->view            = PCView_SysPFMG;
2151:   pc->ops->destroy         = PCDestroy_SysPFMG;
2152:   pc->ops->apply           = PCApply_SysPFMG;
2153:   pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
2154:   pc->ops->setup           = PCSetUp_SysPFMG;

2156:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
2157:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2158:   return(0);
2159: }