Actual source code: hypre.c

petsc-master 2015-04-19
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);

 27:   MPI_Comm comm_hypre;
 28:   char     *hypre_type;

 30:   /* options for Pilut and BoomerAMG*/
 31:   PetscInt maxiter;
 32:   double   tol;

 34:   /* options for Pilut */
 35:   PetscInt factorrowsize;

 37:   /* options for ParaSails */
 38:   PetscInt nlevels;
 39:   double   threshhold;
 40:   double   filter;
 41:   PetscInt sym;
 42:   double   loadbal;
 43:   PetscInt logging;
 44:   PetscInt ruse;
 45:   PetscInt symt;

 47:   /* options for BoomerAMG */
 48:   PetscBool printstatistics;

 50:   /* options for BoomerAMG */
 51:   PetscInt  cycletype;
 52:   PetscInt  maxlevels;
 53:   double    strongthreshold;
 54:   double    maxrowsum;
 55:   PetscInt  gridsweeps[3];
 56:   PetscInt  coarsentype;
 57:   PetscInt  measuretype;
 58:   PetscInt  relaxtype[3];
 59:   double    relaxweight;
 60:   double    outerrelaxweight;
 61:   PetscInt  relaxorder;
 62:   double    truncfactor;
 63:   PetscBool applyrichardson;
 64:   PetscInt  pmax;
 65:   PetscInt  interptype;
 66:   PetscInt  agg_nl;
 67:   PetscInt  agg_num_paths;
 68:   PetscInt  nodal_coarsen;
 69:   PetscBool nodal_relax;
 70:   PetscInt  nodal_relax_levels;

 72:   /* options for AMS */
 73:   PetscInt  ams_print;
 74:   PetscInt  ams_max_iter;
 75:   PetscInt  ams_cycle_type;
 76:   PetscReal ams_tol;
 77:   PetscInt  ams_relax_type;
 78:   PetscInt  ams_relax_times;
 79:   PetscReal ams_relax_weight;
 80:   PetscReal ams_omega;
 81:   PetscInt  ams_amg_alpha_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for vector Poisson */
 82:   PetscReal ams_amg_alpha_theta;   /* AMG strength for vector Poisson */
 83:   PetscInt  ams_amg_beta_opts[5];  /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for scalar Poisson */
 84:   PetscReal ams_amg_beta_theta;    /* AMG strength for scalar Poisson */

 86:   /* additional data */
 87:   HYPRE_IJVector coords[3];
 88:   HYPRE_IJVector constants[3];
 89:   HYPRE_IJMatrix G;
 90:   HYPRE_IJMatrix alpha_Poisson;
 91:   HYPRE_IJMatrix beta_Poisson;
 92:   PetscBool      ams_beta_is_zero;
 93: } PC_HYPRE;

 97: PetscErrorCode PCHYPREGetSolver(PC pc,HYPRE_Solver *hsolver)
 98: {
 99:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

102:   *hsolver = jac->hsolver;
103:   return(0);
104: }

108: static PetscErrorCode PCSetUp_HYPRE(PC pc)
109: {
110:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
111:   PetscErrorCode     ierr;
112:   HYPRE_ParCSRMatrix hmat;
113:   HYPRE_ParVector    bv,xv;
114:   PetscInt           bs;

117:   if (!jac->hypre_type) {
118:     PCHYPRESetType(pc,"boomeramg");
119:   }

121:   if (pc->setupcalled) {
122:     /* always destroy the old matrix and create a new memory;
123:        hope this does not churn the memory too much. The problem
124:        is I do not know if it is possible to put the matrix back to
125:        its initial state so that we can directly copy the values
126:        the second time through. */
127:     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij));
128:     jac->ij = 0;
129:   }

131:   if (!jac->ij) { /* create the matrix the first time through */
132:     MatHYPRE_IJMatrixCreate(pc->pmat,&jac->ij);
133:   }
134:   if (!jac->b) { /* create the vectors the first time through */
135:     Vec x,b;
136:     MatCreateVecs(pc->pmat,&x,&b);
137:     VecHYPRE_IJVectorCreate(x,&jac->x);
138:     VecHYPRE_IJVectorCreate(b,&jac->b);
139:     VecDestroy(&x);
140:     VecDestroy(&b);
141:   }

143:   /* special case for BoomerAMG */
144:   if (jac->setup == HYPRE_BoomerAMGSetup) {
145:     MatGetBlockSize(pc->pmat,&bs);
146:     if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
147:   }
148:   /* special case for AMS */
149:   if (jac->setup == HYPRE_AMSSetup) {
150:     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()");
151:     if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs discrete gradient operator via PCHYPRESetDiscreteGradient");
152:   }
153:   MatHYPRE_IJMatrixCopy(pc->pmat,jac->ij);
154:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
155:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&bv));
156:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&xv));
157:   PetscStackCall("HYPRE_SetupXXX",(*jac->setup)(jac->hsolver,hmat,bv,xv););
158:   return(0);
159: }

161: /*
162:     Replaces the address where the HYPRE vector points to its data with the address of
163:   PETSc's data. Saves the old address so it can be reset when we are finished with it.
164:   Allows use to get the data into a HYPRE vector without the cost of memcopies
165: */
166: #define HYPREReplacePointer(b,newvalue,savedvalue) { \
167:     hypre_ParVector *par_vector   = (hypre_ParVector*)hypre_IJVectorObject(((hypre_IJVector*)b)); \
168:     hypre_Vector    *local_vector = hypre_ParVectorLocalVector(par_vector); \
169:     savedvalue         = local_vector->data; \
170:     local_vector->data = newvalue;          \
171: }

175: static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
176: {
177:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
178:   PetscErrorCode     ierr;
179:   HYPRE_ParCSRMatrix hmat;
180:   PetscScalar        *xv;
181:   const PetscScalar  *bv,*sbv;
182:   HYPRE_ParVector    jbv,jxv;
183:   PetscScalar        *sxv;
184:   PetscInt           hierr;

187:   PetscCitationsRegister(hypreCitation,&cite);
188:   if (!jac->applyrichardson) {VecSet(x,0.0);}
189:   VecGetArrayRead(b,&bv);
190:   VecGetArray(x,&xv);
191:   HYPREReplacePointer(jac->b,(PetscScalar*)bv,sbv);
192:   HYPREReplacePointer(jac->x,xv,sxv);

194:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
195:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
196:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));
197:   PetscStackCall("Hypre solve",h(*jac->solve)(jac->hsolver,hmat,jbv,jxv);
198:                                if (hierr && hierr != HYPRE_ERROR_CONV) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
199:                                if (hierr) hypre__global_error = 0;);

201:   HYPREReplacePointer(jac->b,(PetscScalar*)sbv,bv);
202:   HYPREReplacePointer(jac->x,sxv,xv);
203:   VecRestoreArray(x,&xv);
204:   VecRestoreArrayRead(b,&bv);
205:   return(0);
206: }

210: static PetscErrorCode PCDestroy_HYPRE(PC pc)
211: {
212:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;

216:   if (jac->ij) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij));
217:   if (jac->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->b));
218:   if (jac->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->x));
219:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0]));
220:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1]));
221:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2]));
222:   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0]));
223:   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1]));
224:   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2]));
225:   if (jac->G) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->G));
226:   if (jac->alpha_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->alpha_Poisson));
227:   if (jac->beta_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->beta_Poisson));
228:   if (jac->destroy) PetscStackCall("HYPRE_DestroyXXX",(*jac->destroy)(jac->hsolver););
229:   PetscFree(jac->hypre_type);
230:   if (jac->comm_hypre != MPI_COMM_NULL) { MPI_Comm_free(&(jac->comm_hypre));}
231:   PetscFree(pc->data);

233:   PetscObjectChangeTypeName((PetscObject)pc,0);
234:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);
235:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);
236:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetCoordinates_C",NULL);
237:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",NULL);
238:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetConstantEdgeVectors_C",NULL);
239:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetAlphaPoissonMatrix_C",NULL);
240:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetBetaPoissonMatrix_C",NULL);
241:   return(0);
242: }

244: /* --------------------------------------------------------------------------------------------*/
247: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PetscOptions *PetscOptionsObject,PC pc)
248: {
249:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
251:   PetscBool      flag;

254:   PetscOptionsHead(PetscOptionsObject,"HYPRE Pilut Options");
255:   PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
256:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter));
257:   PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
258:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol));
259:   PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
260:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize));
261:   PetscOptionsTail();
262:   return(0);
263: }

267: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
268: {
269:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
271:   PetscBool      iascii;

274:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
275:   if (iascii) {
276:     PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut preconditioning\n");
277:     if (jac->maxiter != PETSC_DEFAULT) {
278:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: maximum number of iterations %d\n",jac->maxiter);
279:     } else {
280:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default maximum number of iterations \n");
281:     }
282:     if (jac->tol != PETSC_DEFAULT) {
283:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: drop tolerance %g\n",(double)jac->tol);
284:     } else {
285:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default drop tolerance \n");
286:     }
287:     if (jac->factorrowsize != PETSC_DEFAULT) {
288:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: factor row size %d\n",jac->factorrowsize);
289:     } else {
290:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default factor row size \n");
291:     }
292:   }
293:   return(0);
294: }

296: /* --------------------------------------------------------------------------------------------*/

300: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
301: {
302:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
303:   PetscErrorCode     ierr;
304:   HYPRE_ParCSRMatrix hmat;
305:   PetscScalar        *xv;
306:   const PetscScalar  *bv;
307:   HYPRE_ParVector    jbv,jxv;
308:   PetscScalar        *sbv,*sxv;
309:   PetscInt           hierr;

312:   PetscCitationsRegister(hypreCitation,&cite);
313:   VecSet(x,0.0);
314:   VecGetArrayRead(b,&bv);
315:   VecGetArray(x,&xv);
316:   HYPREReplacePointer(jac->b,(PetscScalar*)bv,sbv);
317:   HYPREReplacePointer(jac->x,xv,sxv);

319:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
320:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
321:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));

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

328:   HYPREReplacePointer(jac->b,sbv,bv);
329:   HYPREReplacePointer(jac->x,sxv,xv);
330:   VecRestoreArray(x,&xv);
331:   VecRestoreArrayRead(b,&bv);
332:   return(0);
333: }

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

338: static const char *HYPREBoomerAMGCycleType[]   = {"","V","W"};
339: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
340: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
341: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
342: static const char *HYPREBoomerAMGRelaxType[]   = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi",
343:                                                   "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi",
344:                                                   "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination",
345:                                                   "" /* 10 */, "" /* 11 */, "" /* 12 */, "" /* 13 */, "" /* 14 */,
346:                                                   "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"};
347: static const char *HYPREBoomerAMGInterpType[]  = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
348:                                                   "ext+i-cc", "standard", "standard-wts", "", "", "FF", "FF1"};
351: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PetscOptions *PetscOptionsObject,PC pc)
352: {
353:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
355:   PetscInt       n,indx,level;
356:   PetscBool      flg, tmp_truth;
357:   double         tmpdbl, twodbl[2];

360:   PetscOptionsHead(PetscOptionsObject,"HYPRE BoomerAMG Options");
361:   PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);
362:   if (flg) {
363:     jac->cycletype = indx+1;
364:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
365:   }
366:   PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
367:   if (flg) {
368:     if (jac->maxlevels < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
369:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
370:   }
371:   PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
372:   if (flg) {
373:     if (jac->maxiter < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
374:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
375:   }
376:   PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);
377:   if (flg) {
378:     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);
379:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
380:   }

382:   PetscOptionsScalar("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);
383:   if (flg) {
384:     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);
385:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
386:   }

388:   PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);
389:   if (flg) {
390:     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);
391:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
392:   }

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

398:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
399:   }


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

406:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
407:   }


410:   PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
411:   if (flg) {
412:     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);
413:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
414:   }
415:   PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
416:   if (flg) {
417:     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);
418:     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);
419:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
420:   }

422:   /* Grid sweeps */
423:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);
424:   if (flg) {
425:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx));
426:     /* modify the jac structure so we can view the updated options with PC_View */
427:     jac->gridsweeps[0] = indx;
428:     jac->gridsweeps[1] = indx;
429:     /*defaults coarse to 1 */
430:     jac->gridsweeps[2] = 1;
431:   }

433:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);
434:   if (flg) {
435:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1));
436:     jac->gridsweeps[0] = indx;
437:   }
438:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);
439:   if (flg) {
440:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2));
441:     jac->gridsweeps[1] = indx;
442:   }
443:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);
444:   if (flg) {
445:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3));
446:     jac->gridsweeps[2] = indx;
447:   }

449:   /* Relax type */
450:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
451:   if (flg) {
452:     jac->relaxtype[0] = jac->relaxtype[1]  = indx;
453:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx));
454:     /* by default, coarse type set to 9 */
455:     jac->relaxtype[2] = 9;

457:   }
458:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
459:   if (flg) {
460:     jac->relaxtype[0] = indx;
461:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1));
462:   }
463:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
464:   if (flg) {
465:     jac->relaxtype[1] = indx;
466:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2));
467:   }
468:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);
469:   if (flg) {
470:     jac->relaxtype[2] = indx;
471:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3));
472:   }

474:   /* Relaxation Weight */
475:   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);
476:   if (flg) {
477:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl));
478:     jac->relaxweight = tmpdbl;
479:   }

481:   n         = 2;
482:   twodbl[0] = twodbl[1] = 1.0;
483:   PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
484:   if (flg) {
485:     if (n == 2) {
486:       indx =  (int)PetscAbsReal(twodbl[1]);
487:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx));
488:     } 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);
489:   }

491:   /* Outer relaxation Weight */
492:   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);
493:   if (flg) {
494:     PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl));
495:     jac->outerrelaxweight = tmpdbl;
496:   }

498:   n         = 2;
499:   twodbl[0] = twodbl[1] = 1.0;
500:   PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
501:   if (flg) {
502:     if (n == 2) {
503:       indx =  (int)PetscAbsReal(twodbl[1]);
504:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx));
505:     } 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);
506:   }

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

511:   if (flg && tmp_truth) {
512:     jac->relaxorder = 0;
513:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
514:   }
515:   PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);
516:   if (flg) {
517:     jac->measuretype = indx;
518:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
519:   }
520:   /* update list length 3/07 */
521:   PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);
522:   if (flg) {
523:     jac->coarsentype = indx;
524:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
525:   }

527:   /* new 3/07 */
528:   PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);
529:   if (flg) {
530:     jac->interptype = indx;
531:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
532:   }

534:   PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
535:   if (flg) {
536:     level = 3;
537:     PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);

539:     jac->printstatistics = PETSC_TRUE;
540:     PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level));
541:   }

543:   PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);
544:   if (flg) {
545:     level = 3;
546:     PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);

548:     jac->printstatistics = PETSC_TRUE;
549:     PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level));
550:   }

552:   PetscOptionsBool("-pc_hypre_boomeramg_nodal_coarsen", "HYPRE_BoomerAMGSetNodal()", "None", jac->nodal_coarsen ? PETSC_TRUE : PETSC_FALSE, &tmp_truth, &flg);
553:   if (flg && tmp_truth) {
554:     jac->nodal_coarsen = 1;
555:     PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,1));
556:   }

558:   PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);
559:   if (flg && tmp_truth) {
560:     PetscInt tmp_int;
561:     PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);
562:     if (flg) jac->nodal_relax_levels = tmp_int;
563:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6));
564:     PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1));
565:     PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0));
566:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels));
567:   }

569:   PetscOptionsTail();
570:   return(0);
571: }

575: 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)
576: {
577:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
579:   PetscInt       oits;

582:   PetscCitationsRegister(hypreCitation,&cite);
583:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter));
584:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol));
585:   jac->applyrichardson = PETSC_TRUE;
586:   PCApply_HYPRE(pc,b,y);
587:   jac->applyrichardson = PETSC_FALSE;
588:   PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
589:   *outits = oits;
590:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
591:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
592:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
593:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
594:   return(0);
595: }


600: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
601: {
602:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
604:   PetscBool      iascii;

607:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
608:   if (iascii) {
609:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG preconditioning\n");
610:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);
611:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of levels %d\n",jac->maxlevels);
612:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of iterations PER hypre call %d\n",jac->maxiter);
613:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Convergence tolerance PER hypre call %g\n",(double)jac->tol);
614:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Threshold for strong coupling %g\n",(double)jac->strongthreshold);
615:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation truncation factor %g\n",(double)jac->truncfactor);
616:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation: max elements per row %d\n",jac->pmax);
617:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Number of levels of aggressive coarsening %d\n",jac->agg_nl);
618:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Number of paths for aggressive coarsening %d\n",jac->agg_num_paths);

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

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

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

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

633:     if (jac->relaxorder) {
634:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Using CF-relaxation\n");
635:     } else {
636:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Not using CF-relaxation\n");
637:     }
638:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Measure type        %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
639:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Coarsen type        %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
640:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation type  %s\n",HYPREBoomerAMGInterpType[jac->interptype]);
641:     if (jac->nodal_coarsen) {
642:       PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal coarsening (with HYPRE_BOOMERAMGSetNodal())\n");
643:     }
644:     if (jac->nodal_relax) {
645:       PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal relaxation via Schwarz smoothing on levels %d\n",jac->nodal_relax_levels);
646:     }
647:   }
648:   return(0);
649: }

651: /* --------------------------------------------------------------------------------------------*/
654: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PetscOptions *PetscOptionsObject,PC pc)
655: {
656:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
658:   PetscInt       indx;
659:   PetscBool      flag;
660:   const char     *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};

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

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

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

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

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

680:   PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);
681:   if (flag) {
682:     jac->symt = indx;
683:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
684:   }

686:   PetscOptionsTail();
687:   return(0);
688: }

692: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
693: {
694:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
696:   PetscBool      iascii;
697:   const char     *symt = 0;;

700:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
701:   if (iascii) {
702:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails preconditioning\n");
703:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: nlevels %d\n",jac->nlevels);
704:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: threshold %g\n",(double)jac->threshhold);
705:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: filter %g\n",(double)jac->filter);
706:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: load balance %g\n",(double)jac->loadbal);
707:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: reuse nonzero structure %s\n",PetscBools[jac->ruse]);
708:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: print info to screen %s\n",PetscBools[jac->logging]);
709:     if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
710:     else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
711:     else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
712:     else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
713:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: %s\n",symt);
714:   }
715:   return(0);
716: }
717: /* --------------------------------------------------------------------------------------------*/
720: static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PetscOptions *PetscOptionsObject,PC pc)
721: {
722:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
724:   PetscInt       n;
725:   PetscBool      flag,flag2,flag3,flag4;

728:   PetscOptionsHead(PetscOptionsObject,"HYPRE AMS Options");
729:   PetscOptionsInt("-pc_hypre_ams_print_level","Debugging output level for AMS","None",jac->ams_print,&jac->ams_print,&flag);
730:   if (flag) PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->ams_print));
731:   PetscOptionsInt("-pc_hypre_ams_max_iter","Maximum number of AMS multigrid iterations within PCApply","None",jac->ams_max_iter,&jac->ams_max_iter,&flag);
732:   if (flag) PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->ams_max_iter));
733:   PetscOptionsInt("-pc_hypre_ams_cycle_type","Cycle type for AMS multigrid","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag);
734:   if (flag) PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
735:   PetscOptionsReal("-pc_hypre_ams_tol","Error tolerance for AMS multigrid","None",jac->ams_tol,&jac->ams_tol,&flag);
736:   if (flag) PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->ams_tol));
737:   PetscOptionsInt("-pc_hypre_ams_relax_type","Relaxation type for AMS smoother","None",jac->ams_relax_type,&jac->ams_relax_type,&flag);
738:   PetscOptionsInt("-pc_hypre_ams_relax_times","Number of relaxation steps for AMS smoother","None",jac->ams_relax_times,&jac->ams_relax_times,&flag2);
739:   PetscOptionsReal("-pc_hypre_ams_relax_weight","Relaxation weight for AMS smoother","None",jac->ams_relax_weight,&jac->ams_relax_weight,&flag3);
740:   PetscOptionsReal("-pc_hypre_ams_omega","SSOR coefficient for AMS smoother","None",jac->ams_omega,&jac->ams_omega,&flag4);
741:   if (flag || flag2 || flag3 || flag4) {
742:     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->ams_relax_type,
743:                                                                       jac->ams_relax_times,
744:                                                                       jac->ams_relax_weight,
745:                                                                       jac->ams_omega));
746:   }
747:   PetscOptionsReal("-pc_hypre_ams_amg_alpha_theta","Threshold for strong coupling of vector Poisson AMG solver","None",jac->ams_amg_alpha_theta,&jac->ams_amg_alpha_theta,&flag);
748:   n = 5;
749:   PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options","AMG options for vector Poisson","None",jac->ams_amg_alpha_opts,&n,&flag2);
750:   if (flag || flag2) {
751:     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->ams_amg_alpha_opts[0],       /* AMG coarsen type */
752:                                                                      jac->ams_amg_alpha_opts[1],       /* AMG agg_levels */
753:                                                                      jac->ams_amg_alpha_opts[2],       /* AMG relax_type */
754:                                                                      jac->ams_amg_alpha_theta,
755:                                                                      jac->ams_amg_alpha_opts[3],       /* AMG interp_type */
756:                                                                      jac->ams_amg_alpha_opts[4]));     /* AMG Pmax */
757:   }
758:   PetscOptionsReal("-pc_hypre_ams_amg_beta_theta","Threshold for strong coupling of scalar Poisson AMG solver","None",jac->ams_amg_beta_theta,&jac->ams_amg_beta_theta,&flag);
759:   n = 5;
760:   PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options","AMG options for scalar Poisson solver","None",jac->ams_amg_beta_opts,&n,&flag2);
761:   if (flag || flag2) {
762:     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->ams_amg_beta_opts[0],       /* AMG coarsen type */
763:                                                                     jac->ams_amg_beta_opts[1],       /* AMG agg_levels */
764:                                                                     jac->ams_amg_beta_opts[2],       /* AMG relax_type */
765:                                                                     jac->ams_amg_beta_theta,
766:                                                                     jac->ams_amg_beta_opts[3],       /* AMG interp_type */
767:                                                                     jac->ams_amg_beta_opts[4]));     /* AMG Pmax */
768:   }
769:   PetscOptionsTail();
770:   return(0);
771: }

775: static PetscErrorCode PCView_HYPRE_AMS(PC pc,PetscViewer viewer)
776: {
777:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
779:   PetscBool      iascii;

782:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
783:   if (iascii) {
784:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS preconditioning\n");
785:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace iterations per application %d\n",jac->ams_max_iter);
786:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace cycle type %d\n",jac->ams_cycle_type);
787:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace iteration tolerance %g\n",jac->ams_tol);
788:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother type %d\n",jac->ams_relax_type);
789:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: number of smoothing steps %d\n",jac->ams_relax_times);
790:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother weight %g\n",jac->ams_relax_weight);
791:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother omega %g\n",jac->ams_omega);
792:     if (jac->alpha_Poisson) {
793:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: vector Poisson solver (passed in by user)\n");
794:     } else {
795:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: vector Poisson solver (computed) \n");
796:     }
797:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG coarsening type %d\n",jac->ams_amg_alpha_opts[0]);
798:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG levels of aggressive coarsening %d\n",jac->ams_amg_alpha_opts[1]);
799:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG relaxation type %d\n",jac->ams_amg_alpha_opts[2]);
800:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG interpolation type %d\n",jac->ams_amg_alpha_opts[3]);
801:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->ams_amg_alpha_opts[4]);
802:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG strength threshold %g\n",jac->ams_amg_alpha_theta);
803:     if (!jac->ams_beta_is_zero) {
804:       if (jac->beta_Poisson) {
805:         PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: scalar Poisson solver (passed in by user)\n");
806:       } else {
807:         PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: scalar Poisson solver (computed) \n");
808:       }
809:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG coarsening type %d\n",jac->ams_amg_beta_opts[0]);
810:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG levels of aggressive coarsening %d\n",jac->ams_amg_beta_opts[1]);
811:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG relaxation type %d\n",jac->ams_amg_beta_opts[2]);
812:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG interpolation type %d\n",jac->ams_amg_beta_opts[3]);
813:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->ams_amg_beta_opts[4]);
814:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG strength threshold %g\n",jac->ams_amg_beta_theta);
815:     }
816:   }
817:   return(0);
818: }

822: static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE_AMS(PC pc, Mat G)
823: {
824:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
825:   HYPRE_ParCSRMatrix parcsr_G;
826:   PetscErrorCode     ierr;

829:   /* throw away any discrete gradient if already set */
830:   if (jac->G) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->G));
831:   MatHYPRE_IJMatrixCreate(G,&jac->G);
832:   MatHYPRE_IJMatrixCopy(G,jac->G);
833:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->G,(void**)(&parcsr_G)));
834:   PetscStackCallStandard(HYPRE_AMSSetDiscreteGradient,(jac->hsolver,parcsr_G));
835:   return(0);
836: }

840: /*@
841:  PCHYPRESetDiscreteGradient - Set discrete gradient matrix

843:    Collective on PC

845:    Input Parameters:
846: +  pc - the preconditioning context
847: -  G - the discrete gradient

849:    Level: intermediate

851:    Notes: G should have as many rows as the number of edges and as many columns as the number of vertices in the mesh
852:           Each row of G has 2 nonzeros, with column indexes being the global indexes of edge's endpoints: matrix values are +1 and -1 depending on edge orientation

854: .seealso:
855: @*/
856: PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G)
857: {

864:   PetscTryMethod(pc,"PCHYPRESetDiscreteGradient_C",(PC,Mat),(pc,G));
865:   return(0);
866: }

870: static PetscErrorCode PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS(PC pc, Mat A)
871: {
872:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
873:   HYPRE_ParCSRMatrix parcsr_alpha_Poisson;
874:   PetscErrorCode     ierr;

877:   /* throw away any matrix if already set */
878:   if (jac->alpha_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->alpha_Poisson));
879:   MatHYPRE_IJMatrixCreate(A,&jac->alpha_Poisson);
880:   MatHYPRE_IJMatrixCopy(A,jac->alpha_Poisson);
881:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->alpha_Poisson,(void**)(&parcsr_alpha_Poisson)));
882:   PetscStackCallStandard(HYPRE_AMSSetAlphaPoissonMatrix,(jac->hsolver,parcsr_alpha_Poisson));
883:   return(0);
884: }

888: /*@
889:  PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix

891:    Collective on PC

893:    Input Parameters:
894: +  pc - the preconditioning context
895: -  A - the matrix

897:    Level: intermediate

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

901: .seealso:
902: @*/
903: PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A)
904: {

911:   PetscTryMethod(pc,"PCHYPRESetAlphaPoissonMatrix_C",(PC,Mat),(pc,A));
912:   return(0);
913: }

917: static PetscErrorCode PCHYPRESetBetaPoissonMatrix_HYPRE_AMS(PC pc, Mat A)
918: {
919:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
920:   HYPRE_ParCSRMatrix parcsr_beta_Poisson;
921:   PetscErrorCode     ierr;

924:   if (!A) {
925:     jac->ams_beta_is_zero = PETSC_TRUE;
926:     return(0);
927:   }
928:   jac->ams_beta_is_zero = PETSC_FALSE;
929:   /* throw away any matrix if already set */
930:   if (jac->beta_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->beta_Poisson));
931:   MatHYPRE_IJMatrixCreate(A,&jac->beta_Poisson);
932:   MatHYPRE_IJMatrixCopy(A,jac->beta_Poisson);
933:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->beta_Poisson,(void**)(&parcsr_beta_Poisson)));
934:   PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,parcsr_beta_Poisson));
935:   return(0);
936: }

940: /*@
941:  PCHYPRESetBetaPoissonMatrix - Set Poisson matrix

943:    Collective on PC

945:    Input Parameters:
946: +  pc - the preconditioning context
947: -  A - the matrix

949:    Level: intermediate

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

954: .seealso:
955: @*/
956: PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A)
957: {

962:   if (A) {
965:   }
966:   PetscTryMethod(pc,"PCHYPRESetBetaPoissonMatrix_C",(PC,Mat),(pc,A));
967:   return(0);
968: }

972: static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE_AMS(PC pc,Vec ozz, Vec zoz, Vec zzo)
973: {
974:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
975:   HYPRE_ParVector    par_ozz,par_zoz,par_zzo;
976:   PetscInt           dim;
977:   PetscErrorCode     ierr;

980:   /* throw away any vector if already set */
981:   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0]));
982:   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1]));
983:   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2]));
984:   jac->constants[0] = NULL;
985:   jac->constants[1] = NULL;
986:   jac->constants[2] = NULL;
987:   VecHYPRE_IJVectorCreate(ozz,&jac->constants[0]);
988:   VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);
989:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[0],(void**)(&par_ozz)));
990:   VecHYPRE_IJVectorCreate(zoz,&jac->constants[1]);
991:   VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);
992:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[1],(void**)(&par_zoz)));
993:   dim = 2;
994:   if (zzo) {
995:     VecHYPRE_IJVectorCreate(zzo,&jac->constants[2]);
996:     VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);
997:     PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[2],(void**)(&par_zzo)));
998:     dim++;
999:   }
1000:   PetscStackCallStandard(HYPRE_AMSSetEdgeConstantVectors,(jac->hsolver,par_ozz,par_zoz,par_zzo));
1001:   PetscStackCallStandard(HYPRE_AMSSetDimension,(jac->hsolver,dim));
1002:   return(0);
1003: }

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

1010:    Collective on PC

1012:    Input Parameters:
1013: +  pc - the preconditioning context
1014: -  ozz - vector representing (1,0,0) (or (1,0) in 2D)
1015: -  zoz - vector representing (0,1,0) (or (0,1) in 2D)
1016: -  zzo - vector representing (0,0,1) (use NULL in 2D)

1018:    Level: intermediate

1020:    Notes:

1022: .seealso:
1023: @*/
1024: PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo)
1025: {

1036:   PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));
1037:   return(0);
1038: }

1042: static PetscErrorCode PCSetCoordinates_HYPRE_AMS(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1043: {
1044:   PC_HYPRE        *jac = (PC_HYPRE*)pc->data;
1045:   Vec             tv;
1046:   HYPRE_ParVector par_coords[3];
1047:   PetscInt        i;
1048:   PetscErrorCode  ierr;

1051:   /* throw away any coordinate vector if already set */
1052:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0]));
1053:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1]));
1054:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2]));
1055:   /* set problem's dimension */
1056:   PetscStackCallStandard(HYPRE_AMSSetDimension,(jac->hsolver,dim));
1057:   /* compute IJ vector for coordinates */
1058:   VecCreate(PetscObjectComm((PetscObject)pc),&tv);
1059:   VecSetType(tv,VECSTANDARD);
1060:   VecSetSizes(tv,nloc,PETSC_DECIDE);
1061:   for (i=0;i<dim;i++) {
1062:     PetscScalar *array;
1063:     PetscInt    j;

1065:     VecHYPRE_IJVectorCreate(tv,&jac->coords[i]);
1066:     VecGetArray(tv,&array);
1067:     for (j=0;j<nloc;j++) {
1068:       array[j] = coords[j*dim+i];
1069:     }
1070:     PetscStackCallStandard(HYPRE_IJVectorSetValues,(jac->coords[i],nloc,NULL,array));
1071:     PetscStackCallStandard(HYPRE_IJVectorAssemble,(jac->coords[i]));
1072:     VecRestoreArray(tv,&array);
1073:   }
1074:   VecDestroy(&tv);
1075:   /* pass parCSR vectors to AMS solver */
1076:   par_coords[0] = NULL;
1077:   par_coords[1] = NULL;
1078:   par_coords[2] = NULL;
1079:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0],(void**)(&par_coords[0])));
1080:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1],(void**)(&par_coords[1])));
1081:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2],(void**)(&par_coords[2])));
1082:   PetscStackCallStandard(HYPRE_AMSSetCoordinateVectors,(jac->hsolver,par_coords[0],par_coords[1],par_coords[2]));
1083:   return(0);
1084: }

1086: /* ---------------------------------------------------------------------------------*/

1090: static PetscErrorCode  PCHYPREGetType_HYPRE(PC pc,const char *name[])
1091: {
1092:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

1095:   *name = jac->hypre_type;
1096:   return(0);
1097: }

1101: static PetscErrorCode  PCHYPRESetType_HYPRE(PC pc,const char name[])
1102: {
1103:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1105:   PetscBool      flag;

1108:   if (jac->hypre_type) {
1109:     PetscStrcmp(jac->hypre_type,name,&flag);
1110:     if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set");
1111:     return(0);
1112:   } else {
1113:     PetscStrallocpy(name, &jac->hypre_type);
1114:   }

1116:   jac->maxiter         = PETSC_DEFAULT;
1117:   jac->tol             = PETSC_DEFAULT;
1118:   jac->printstatistics = PetscLogPrintInfo;

1120:   PetscStrcmp("pilut",jac->hypre_type,&flag);
1121:   if (flag) {
1122:     PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver));
1123:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
1124:     pc->ops->view           = PCView_HYPRE_Pilut;
1125:     jac->destroy            = HYPRE_ParCSRPilutDestroy;
1126:     jac->setup              = HYPRE_ParCSRPilutSetup;
1127:     jac->solve              = HYPRE_ParCSRPilutSolve;
1128:     jac->factorrowsize      = PETSC_DEFAULT;
1129:     return(0);
1130:   }
1131:   PetscStrcmp("parasails",jac->hypre_type,&flag);
1132:   if (flag) {
1133:     PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver));
1134:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
1135:     pc->ops->view           = PCView_HYPRE_ParaSails;
1136:     jac->destroy            = HYPRE_ParaSailsDestroy;
1137:     jac->setup              = HYPRE_ParaSailsSetup;
1138:     jac->solve              = HYPRE_ParaSailsSolve;
1139:     /* initialize */
1140:     jac->nlevels    = 1;
1141:     jac->threshhold = .1;
1142:     jac->filter     = .1;
1143:     jac->loadbal    = 0;
1144:     if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE;
1145:     else jac->logging = (int) PETSC_FALSE;

1147:     jac->ruse = (int) PETSC_FALSE;
1148:     jac->symt = 0;
1149:     PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
1150:     PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
1151:     PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
1152:     PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
1153:     PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
1154:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
1155:     return(0);
1156:   }
1157:   PetscStrcmp("boomeramg",jac->hypre_type,&flag);
1158:   if (flag) {
1159:     HYPRE_BoomerAMGCreate(&jac->hsolver);
1160:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_BoomerAMG;
1161:     pc->ops->view            = PCView_HYPRE_BoomerAMG;
1162:     pc->ops->applytranspose  = PCApplyTranspose_HYPRE_BoomerAMG;
1163:     pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
1164:     jac->destroy             = HYPRE_BoomerAMGDestroy;
1165:     jac->setup               = HYPRE_BoomerAMGSetup;
1166:     jac->solve               = HYPRE_BoomerAMGSolve;
1167:     jac->applyrichardson     = PETSC_FALSE;
1168:     /* these defaults match the hypre defaults */
1169:     jac->cycletype        = 1;
1170:     jac->maxlevels        = 25;
1171:     jac->maxiter          = 1;
1172:     jac->tol              = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
1173:     jac->truncfactor      = 0.0;
1174:     jac->strongthreshold  = .25;
1175:     jac->maxrowsum        = .9;
1176:     jac->coarsentype      = 6;
1177:     jac->measuretype      = 0;
1178:     jac->gridsweeps[0]    = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
1179:     jac->relaxtype[0]     = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */
1180:     jac->relaxtype[2]     = 9; /*G.E. */
1181:     jac->relaxweight      = 1.0;
1182:     jac->outerrelaxweight = 1.0;
1183:     jac->relaxorder       = 1;
1184:     jac->interptype       = 0;
1185:     jac->agg_nl           = 0;
1186:     jac->pmax             = 0;
1187:     jac->truncfactor      = 0.0;
1188:     jac->agg_num_paths    = 1;

1190:     jac->nodal_coarsen      = 0;
1191:     jac->nodal_relax        = PETSC_FALSE;
1192:     jac->nodal_relax_levels = 1;
1193:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
1194:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
1195:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
1196:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
1197:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
1198:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
1199:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
1200:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
1201:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
1202:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
1203:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
1204:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
1205:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
1206:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
1207:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0]));  /*defaults coarse to 9*/
1208:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */
1209:     return(0);
1210:   }
1211:   PetscStrcmp("ams",jac->hypre_type,&flag);
1212:   if (flag) {
1213:     HYPRE_AMSCreate(&jac->hsolver);
1214:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_AMS;
1215:     pc->ops->view            = PCView_HYPRE_AMS;
1216:     jac->destroy             = HYPRE_AMSDestroy;
1217:     jac->setup               = HYPRE_AMSSetup;
1218:     jac->solve               = HYPRE_AMSSolve;
1219:     jac->coords[0]           = NULL;
1220:     jac->coords[1]           = NULL;
1221:     jac->coords[2]           = NULL;
1222:     jac->G                   = NULL;
1223:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */
1224:     jac->ams_print           = 0;
1225:     jac->ams_max_iter        = 1; /* used as a preconditioner */
1226:     jac->ams_cycle_type      = 13;
1227:     jac->ams_tol             = 0.; /* used as a preconditioner */
1228:     /* Smoothing options */
1229:     jac->ams_relax_type      = 2;
1230:     jac->ams_relax_times     = 1;
1231:     jac->ams_relax_weight    = 1.0;
1232:     jac->ams_omega           = 1.0;
1233:     /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1234:     jac->ams_amg_alpha_opts[0] = 10;
1235:     jac->ams_amg_alpha_opts[1] = 1;
1236:     jac->ams_amg_alpha_opts[2] = 8;
1237:     jac->ams_amg_alpha_opts[3] = 6;
1238:     jac->ams_amg_alpha_opts[4] = 4;
1239:     jac->ams_amg_alpha_theta   = 0.25;
1240:     /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1241:     jac->ams_beta_is_zero = PETSC_FALSE;
1242:     jac->ams_amg_beta_opts[0] = 10;
1243:     jac->ams_amg_beta_opts[1] = 1;
1244:     jac->ams_amg_beta_opts[2] = 8;
1245:     jac->ams_amg_beta_opts[3] = 6;
1246:     jac->ams_amg_beta_opts[4] = 4;
1247:     jac->ams_amg_beta_theta   = 0.25;
1248:     PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->ams_print));
1249:     PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->ams_max_iter));
1250:     PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1251:     PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->ams_tol));
1252:     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->ams_relax_type,
1253:                                                                       jac->ams_relax_times,
1254:                                                                       jac->ams_relax_weight,
1255:                                                                       jac->ams_omega));
1256:     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->ams_amg_alpha_opts[0],       /* AMG coarsen type */
1257:                                                                      jac->ams_amg_alpha_opts[1],       /* AMG agg_levels */
1258:                                                                      jac->ams_amg_alpha_opts[2],       /* AMG relax_type */
1259:                                                                      jac->ams_amg_alpha_theta,
1260:                                                                      jac->ams_amg_alpha_opts[3],       /* AMG interp_type */
1261:                                                                      jac->ams_amg_alpha_opts[4]));     /* AMG Pmax */
1262:     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->ams_amg_beta_opts[0],       /* AMG coarsen type */
1263:                                                                     jac->ams_amg_beta_opts[1],       /* AMG agg_levels */
1264:                                                                     jac->ams_amg_beta_opts[2],       /* AMG relax_type */
1265:                                                                     jac->ams_amg_beta_theta,
1266:                                                                     jac->ams_amg_beta_opts[3],       /* AMG interp_type */
1267:                                                                     jac->ams_amg_beta_opts[4]));     /* AMG Pmax */
1268:     PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE_AMS);
1269:     PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE_AMS);
1270:     PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE_AMS);
1271:     PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetAlphaPoissonMatrix_C",PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS);
1272:     PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetBetaPoissonMatrix_C",PCHYPRESetBetaPoissonMatrix_HYPRE_AMS);
1273:     return(0);
1274:   }
1275:   PetscFree(jac->hypre_type);

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

1282: /*
1283:     It only gets here if the HYPRE type has not been set before the call to
1284:    ...SetFromOptions() which actually is most of the time
1285: */
1288: static PetscErrorCode PCSetFromOptions_HYPRE(PetscOptions *PetscOptionsObject,PC pc)
1289: {
1291:   PetscInt       indx;
1292:   const char     *type[] = {"pilut","parasails","boomeramg","ams"};
1293:   PetscBool      flg;

1296:   PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");
1297:   PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"boomeramg",&indx,&flg);
1298:   if (flg) {
1299:     PCHYPRESetType_HYPRE(pc,type[indx]);
1300:   } else {
1301:     PCHYPRESetType_HYPRE(pc,"boomeramg");
1302:   }
1303:   if (pc->ops->setfromoptions) {
1304:     pc->ops->setfromoptions(PetscOptionsObject,pc);
1305:   }
1306:   PetscOptionsTail();
1307:   return(0);
1308: }

1312: /*@C
1313:      PCHYPRESetType - Sets which hypre preconditioner you wish to use

1315:    Input Parameters:
1316: +     pc - the preconditioner context
1317: -     name - either  pilut, parasails, boomeramg, ams

1319:    Options Database Keys:
1320:    -pc_hypre_type - One of pilut, parasails, boomeramg, ams

1322:    Level: intermediate

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

1327: @*/
1328: PetscErrorCode  PCHYPRESetType(PC pc,const char name[])
1329: {

1335:   PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));
1336:   return(0);
1337: }

1341: /*@C
1342:      PCHYPREGetType - Gets which hypre preconditioner you are using

1344:    Input Parameter:
1345: .     pc - the preconditioner context

1347:    Output Parameter:
1348: .     name - either  pilut, parasails, boomeramg, ams

1350:    Level: intermediate

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

1355: @*/
1356: PetscErrorCode  PCHYPREGetType(PC pc,const char *name[])
1357: {

1363:   PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));
1364:   return(0);
1365: }

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

1370:    Options Database Keys:
1371: +   -pc_hypre_type - One of pilut, parasails, boomeramg, ams
1372: -   Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX
1373:           preconditioner

1375:    Level: intermediate

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

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

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

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

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

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

1404: M*/

1408: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
1409: {
1410:   PC_HYPRE       *jac;

1414:   PetscNewLog(pc,&jac);

1416:   pc->data                = jac;
1417:   pc->ops->destroy        = PCDestroy_HYPRE;
1418:   pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
1419:   pc->ops->setup          = PCSetUp_HYPRE;
1420:   pc->ops->apply          = PCApply_HYPRE;
1421:   jac->comm_hypre         = MPI_COMM_NULL;
1422:   jac->hypre_type         = NULL;
1423:   jac->coords[0]          = NULL;
1424:   jac->coords[1]          = NULL;
1425:   jac->coords[2]          = NULL;
1426:   jac->constants[0]       = NULL;
1427:   jac->constants[1]       = NULL;
1428:   jac->constants[2]       = NULL;
1429:   /* duplicate communicator for hypre */
1430:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1431:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);
1432:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);
1433:   return(0);
1434: }

1436: /* ---------------------------------------------------------------------------------------------------------------------------------*/

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

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

1445:   /* keep copy of PFMG options used so may view them */
1446:   PetscInt its;
1447:   double   tol;
1448:   PetscInt relax_type;
1449:   PetscInt rap_type;
1450:   PetscInt num_pre_relax,num_post_relax;
1451:   PetscInt max_levels;
1452: } PC_PFMG;

1456: PetscErrorCode PCDestroy_PFMG(PC pc)
1457: {
1459:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

1462:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1463:   MPI_Comm_free(&ex->hcomm);
1464:   PetscFree(pc->data);
1465:   return(0);
1466: }

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

1473: PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
1474: {
1476:   PetscBool      iascii;
1477:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

1480:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1481:   if (iascii) {
1482:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG preconditioning\n");
1483:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: max iterations %d\n",ex->its);
1484:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: tolerance %g\n",ex->tol);
1485:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
1486:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: RAP type %s\n",PFMGRAPType[ex->rap_type]);
1487:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
1488:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: max levels %d\n",ex->max_levels);
1489:   }
1490:   return(0);
1491: }


1496: PetscErrorCode PCSetFromOptions_PFMG(PetscOptions *PetscOptionsObject,PC pc)
1497: {
1499:   PC_PFMG        *ex = (PC_PFMG*) pc->data;
1500:   PetscBool      flg = PETSC_FALSE;

1503:   PetscOptionsHead(PetscOptionsObject,"PFMG options");
1504:   PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);
1505:   if (flg) {
1506:     int level=3;
1507:     PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level));
1508:   }
1509:   PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);
1510:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its));
1511:   PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
1512:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax));
1513:   PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
1514:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax));

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

1519:   PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);
1520:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol));
1521:   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);
1522:   PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type));
1523:   PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);
1524:   PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type));
1525:   PetscOptionsTail();
1526:   return(0);
1527: }

1531: PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
1532: {
1533:   PetscErrorCode    ierr;
1534:   PC_PFMG           *ex = (PC_PFMG*) pc->data;
1535:   PetscScalar       *yy;
1536:   const PetscScalar *xx;
1537:   PetscInt          ilower[3],iupper[3];
1538:   Mat_HYPREStruct   *mx = (Mat_HYPREStruct*)(pc->pmat->data);

1541:   PetscCitationsRegister(hypreCitation,&cite);
1542:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
1543:   iupper[0] += ilower[0] - 1;
1544:   iupper[1] += ilower[1] - 1;
1545:   iupper[2] += ilower[2] - 1;

1547:   /* copy x values over to hypre */
1548:   PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
1549:   VecGetArrayRead(x,&xx);
1550:   PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx));
1551:   VecRestoreArrayRead(x,&xx);
1552:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
1553:   PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx));

1555:   /* copy solution values back to PETSc */
1556:   VecGetArray(y,&yy);
1557:   PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy));
1558:   VecRestoreArray(y,&yy);
1559:   return(0);
1560: }

1564: 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)
1565: {
1566:   PC_PFMG        *jac = (PC_PFMG*)pc->data;
1568:   PetscInt       oits;

1571:   PetscCitationsRegister(hypreCitation,&cite);
1572:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its));
1573:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol));

1575:   PCApply_PFMG(pc,b,y);
1576:   PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
1577:   *outits = oits;
1578:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1579:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
1580:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol));
1581:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its));
1582:   return(0);
1583: }


1588: PetscErrorCode PCSetUp_PFMG(PC pc)
1589: {
1590:   PetscErrorCode  ierr;
1591:   PC_PFMG         *ex = (PC_PFMG*) pc->data;
1592:   Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
1593:   PetscBool       flg;

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

1599:   /* create the hypre solver object and set its information */
1600:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1601:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1602:   PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
1603:   PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver));
1604:   return(0);
1605: }


1608: /*MC
1609:      PCPFMG - the hypre PFMG multigrid solver

1611:    Level: advanced

1613:    Options Database:
1614: + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner
1615: . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
1616: . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid
1617: . -pc_pfmg_tol <tol> tolerance of PFMG
1618: . -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
1619: - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin

1621:    Notes:  This is for CELL-centered descretizations

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

1626: .seealso:  PCMG, MATHYPRESTRUCT
1627: M*/

1631: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
1632: {
1634:   PC_PFMG        *ex;

1637:   PetscNew(&ex); \
1638:   pc->data = ex;

1640:   ex->its            = 1;
1641:   ex->tol            = 1.e-8;
1642:   ex->relax_type     = 1;
1643:   ex->rap_type       = 0;
1644:   ex->num_pre_relax  = 1;
1645:   ex->num_post_relax = 1;
1646:   ex->max_levels     = 0;

1648:   pc->ops->setfromoptions  = PCSetFromOptions_PFMG;
1649:   pc->ops->view            = PCView_PFMG;
1650:   pc->ops->destroy         = PCDestroy_PFMG;
1651:   pc->ops->apply           = PCApply_PFMG;
1652:   pc->ops->applyrichardson = PCApplyRichardson_PFMG;
1653:   pc->ops->setup           = PCSetUp_PFMG;

1655:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
1656:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1657:   return(0);
1658: }

1660: /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/

1662: /* we know we are working with a HYPRE_SStructMatrix */
1663: typedef struct {
1664:   MPI_Comm            hcomm;       /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
1665:   HYPRE_SStructSolver ss_solver;

1667:   /* keep copy of SYSPFMG options used so may view them */
1668:   PetscInt its;
1669:   double   tol;
1670:   PetscInt relax_type;
1671:   PetscInt num_pre_relax,num_post_relax;
1672: } PC_SysPFMG;

1676: PetscErrorCode PCDestroy_SysPFMG(PC pc)
1677: {
1679:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

1682:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
1683:   MPI_Comm_free(&ex->hcomm);
1684:   PetscFree(pc->data);
1685:   return(0);
1686: }

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

1692: PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
1693: {
1695:   PetscBool      iascii;
1696:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

1699:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1700:   if (iascii) {
1701:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG preconditioning\n");
1702:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: max iterations %d\n",ex->its);
1703:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: tolerance %g\n",ex->tol);
1704:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
1705:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
1706:   }
1707:   return(0);
1708: }


1713: PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptions *PetscOptionsObject,PC pc)
1714: {
1716:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;
1717:   PetscBool      flg = PETSC_FALSE;

1720:   PetscOptionsHead(PetscOptionsObject,"SysPFMG options");
1721:   PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);
1722:   if (flg) {
1723:     int level=3;
1724:     PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level));
1725:   }
1726:   PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);
1727:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its));
1728:   PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
1729:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax));
1730:   PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
1731:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax));

1733:   PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);
1734:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol));
1735:   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);
1736:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type));
1737:   PetscOptionsTail();
1738:   return(0);
1739: }

1743: PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
1744: {
1745:   PetscErrorCode    ierr;
1746:   PC_SysPFMG        *ex = (PC_SysPFMG*) pc->data;
1747:   PetscScalar       *yy;
1748:   const PetscScalar *xx;
1749:   PetscInt          ilower[3],iupper[3];
1750:   Mat_HYPRESStruct  *mx     = (Mat_HYPRESStruct*)(pc->pmat->data);
1751:   PetscInt          ordering= mx->dofs_order;
1752:   PetscInt          nvars   = mx->nvars;
1753:   PetscInt          part    = 0;
1754:   PetscInt          size;
1755:   PetscInt          i;

1758:   PetscCitationsRegister(hypreCitation,&cite);
1759:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
1760:   iupper[0] += ilower[0] - 1;
1761:   iupper[1] += ilower[1] - 1;
1762:   iupper[2] += ilower[2] - 1;

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

1767:   /* copy x values over to hypre for variable ordering */
1768:   if (ordering) {
1769:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
1770:     VecGetArrayRead(x,&xx);
1771:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)xx+(size*i)));
1772:     VecRestoreArrayRead(x,&xx);
1773:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
1774:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
1775:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

1777:     /* copy solution values back to PETSc */
1778:     VecGetArray(y,&yy);
1779:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i)));
1780:     VecRestoreArray(y,&yy);
1781:   } else {      /* nodal ordering must be mapped to variable ordering for sys_pfmg */
1782:     PetscScalar *z;
1783:     PetscInt    j, k;

1785:     PetscMalloc1(nvars*size,&z);
1786:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
1787:     VecGetArrayRead(x,&xx);

1789:     /* transform nodal to hypre's variable ordering for sys_pfmg */
1790:     for (i= 0; i< size; i++) {
1791:       k= i*nvars;
1792:       for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
1793:     }
1794:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
1795:     VecRestoreArrayRead(x,&xx);
1796:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
1797:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

1799:     /* copy solution values back to PETSc */
1800:     VecGetArray(y,&yy);
1801:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
1802:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
1803:     for (i= 0; i< size; i++) {
1804:       k= i*nvars;
1805:       for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
1806:     }
1807:     VecRestoreArray(y,&yy);
1808:     PetscFree(z);
1809:   }
1810:   return(0);
1811: }

1815: 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)
1816: {
1817:   PC_SysPFMG     *jac = (PC_SysPFMG*)pc->data;
1819:   PetscInt       oits;

1822:   PetscCitationsRegister(hypreCitation,&cite);
1823:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its));
1824:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol));
1825:   PCApply_SysPFMG(pc,b,y);
1826:   PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,(HYPRE_Int *)&oits));
1827:   *outits = oits;
1828:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1829:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
1830:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol));
1831:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its));
1832:   return(0);
1833: }


1838: PetscErrorCode PCSetUp_SysPFMG(PC pc)
1839: {
1840:   PetscErrorCode   ierr;
1841:   PC_SysPFMG       *ex = (PC_SysPFMG*) pc->data;
1842:   Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
1843:   PetscBool        flg;

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

1849:   /* create the hypre sstruct solver object and set its information */
1850:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
1851:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
1852:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver));
1853:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
1854:   return(0);
1855: }


1858: /*MC
1859:      PCSysPFMG - the hypre SysPFMG multigrid solver

1861:    Level: advanced

1863:    Options Database:
1864: + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner
1865: . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
1866: . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid
1867: . -pc_syspfmg_tol <tol> tolerance of SysPFMG
1868: . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel

1870:    Notes:  This is for CELL-centered descretizations

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

1876: .seealso:  PCMG, MATHYPRESSTRUCT
1877: M*/

1881: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
1882: {
1884:   PC_SysPFMG     *ex;

1887:   PetscNew(&ex); \
1888:   pc->data = ex;

1890:   ex->its            = 1;
1891:   ex->tol            = 1.e-8;
1892:   ex->relax_type     = 1;
1893:   ex->num_pre_relax  = 1;
1894:   ex->num_post_relax = 1;

1896:   pc->ops->setfromoptions  = PCSetFromOptions_SysPFMG;
1897:   pc->ops->view            = PCView_SysPFMG;
1898:   pc->ops->destroy         = PCDestroy_SysPFMG;
1899:   pc->ops->apply           = PCApply_SysPFMG;
1900:   pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
1901:   pc->ops->setup           = PCSetUp_SysPFMG;

1903:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
1904:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
1905:   return(0);
1906: }