Actual source code: hypre.c

petsc-master 2015-07-30
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:   PetscInt  nodal_coarsening;
 77:   PetscInt  vec_interp_variant;
 78:   HYPRE_IJVector  *hmnull;
 79:   HYPRE_ParVector *phmnull;  /* near null space passed to hypre */
 80:   PetscInt        n_hmnull;
 81:   Vec             hmnull_constant;

 83:   /* options for AS (Auxiliary Space preconditioners) */
 84:   PetscInt  as_print;
 85:   PetscInt  as_max_iter;
 86:   PetscReal as_tol;
 87:   PetscInt  as_relax_type;
 88:   PetscInt  as_relax_times;
 89:   PetscReal as_relax_weight;
 90:   PetscReal as_omega;
 91:   PetscInt  as_amg_alpha_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for vector Poisson (AMS) or Curl problem (ADS) */
 92:   PetscReal as_amg_alpha_theta;   /* AMG strength for vector Poisson (AMS) or Curl problem (ADS) */
 93:   PetscInt  as_amg_beta_opts[5];  /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for scalar Poisson (AMS) or vector Poisson (ADS) */
 94:   PetscReal as_amg_beta_theta;    /* AMG strength for scalar Poisson (AMS) or vector Poisson (ADS)  */
 95:   PetscInt  ams_cycle_type;
 96:   PetscInt  ads_cycle_type;

 98:   /* additional data */
 99:   HYPRE_IJVector coords[3];
100:   HYPRE_IJVector constants[3];
101:   HYPRE_IJMatrix G;
102:   HYPRE_IJMatrix C;
103:   HYPRE_IJMatrix alpha_Poisson;
104:   HYPRE_IJMatrix beta_Poisson;
105:   PetscBool      ams_beta_is_zero;
106: } PC_HYPRE;

110: PetscErrorCode PCHYPREGetSolver(PC pc,HYPRE_Solver *hsolver)
111: {
112:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

115:   *hsolver = jac->hsolver;
116:   return(0);
117: }

121: static PetscErrorCode PCSetUp_HYPRE(PC pc)
122: {
123:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
124:   PetscErrorCode     ierr;
125:   HYPRE_ParCSRMatrix hmat;
126:   HYPRE_ParVector    bv,xv;
127:   PetscInt           bs;

130:   if (!jac->hypre_type) {
131:     PCHYPRESetType(pc,"boomeramg");
132:   }

134:   if (pc->setupcalled) {
135:     /* always destroy the old matrix and create a new memory;
136:        hope this does not churn the memory too much. The problem
137:        is I do not know if it is possible to put the matrix back to
138:        its initial state so that we can directly copy the values
139:        the second time through. */
140:     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij));
141:     jac->ij = 0;
142:   }

144:   if (!jac->ij) { /* create the matrix the first time through */
145:     MatHYPRE_IJMatrixCreate(pc->pmat,&jac->ij);
146:   }
147:   if (!jac->b) { /* create the vectors the first time through */
148:     Vec x,b;
149:     MatCreateVecs(pc->pmat,&x,&b);
150:     VecHYPRE_IJVectorCreate(x,&jac->x);
151:     VecHYPRE_IJVectorCreate(b,&jac->b);
152:     VecDestroy(&x);
153:     VecDestroy(&b);
154:   }

156:   /* special case for BoomerAMG */
157:   if (jac->setup == HYPRE_BoomerAMGSetup) {
158:     MatNullSpace    mnull;
159:     PetscBool       has_const;
160:     PetscInt        nvec,i;
161:     const Vec       *vecs;

163:     MatGetBlockSize(pc->pmat,&bs);
164:     if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
165:     MatGetNearNullSpace(pc->mat, &mnull);
166:     if (mnull) {
167:       MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);
168:       PetscMalloc1(nvec+1,&jac->hmnull);
169:       PetscMalloc1(nvec+1,&jac->phmnull);
170:       for (i=0; i<nvec; i++) {
171:         VecHYPRE_IJVectorCreate(vecs[i],&jac->hmnull[i]);
172:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[i],(void**)&jac->phmnull[i]));
173:       }
174:       if (has_const) {
175:         MatCreateVecs(pc->pmat,&jac->hmnull_constant,NULL);
176:         VecSet(jac->hmnull_constant,1);
177:         VecNormalize(jac->hmnull_constant,NULL);
178:         VecHYPRE_IJVectorCreate(jac->hmnull_constant,&jac->hmnull[nvec]);
179:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[nvec],(void**)&jac->phmnull[nvec]));
180:         nvec++;
181:       }
182:       PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVectors,(jac->hsolver,nvec,jac->phmnull));
183:       jac->n_hmnull = nvec;
184:     }
185:   }

187:   /* special case for AMS */
188:   if (jac->setup == HYPRE_AMSSetup) {
189:     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()");
190:     if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs discrete gradient operator via PCHYPRESetDiscreteGradient");
191:   }
192:   /* special case for ADS */
193:   if (jac->setup == HYPRE_ADSSetup) {
194:     if (!jac->coords[0]) {
195:       SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs coordinate vectors via PCSetCoordinates()");
196:     } else if (!jac->coords[1] || !jac->coords[2]) {
197:       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");
198:     }
199:     if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs discrete gradient operator via PCHYPRESetDiscreteGradient");
200:     if (!jac->C) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs discrete curl operator via PCHYPRESetDiscreteGradient");
201:   }
202:   MatHYPRE_IJMatrixCopy(pc->pmat,jac->ij);
203:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
204:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&bv));
205:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&xv));
206:   PetscStackCall("HYPRE_SetupXXX",(*jac->setup)(jac->hsolver,hmat,bv,xv););
207:   return(0);
208: }

210: /*
211:     Replaces the address where the HYPRE vector points to its data with the address of
212:   PETSc's data. Saves the old address so it can be reset when we are finished with it.
213:   Allows use to get the data into a HYPRE vector without the cost of memcopies
214: */
215: #define HYPREReplacePointer(b,newvalue,savedvalue) { \
216:     hypre_ParVector *par_vector   = (hypre_ParVector*)hypre_IJVectorObject(((hypre_IJVector*)b)); \
217:     hypre_Vector    *local_vector = hypre_ParVectorLocalVector(par_vector); \
218:     savedvalue         = local_vector->data; \
219:     local_vector->data = newvalue;          \
220: }

224: static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
225: {
226:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
227:   PetscErrorCode     ierr;
228:   HYPRE_ParCSRMatrix hmat;
229:   PetscScalar        *xv;
230:   const PetscScalar  *bv,*sbv;
231:   HYPRE_ParVector    jbv,jxv;
232:   PetscScalar        *sxv;
233:   PetscInt           hierr;

236:   PetscCitationsRegister(hypreCitation,&cite);
237:   if (!jac->applyrichardson) {VecSet(x,0.0);}
238:   VecGetArrayRead(b,&bv);
239:   VecGetArray(x,&xv);
240:   HYPREReplacePointer(jac->b,(PetscScalar*)bv,sbv);
241:   HYPREReplacePointer(jac->x,xv,sxv);

243:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
244:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
245:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));
246:   PetscStackCall("Hypre solve",h(*jac->solve)(jac->hsolver,hmat,jbv,jxv);
247:                                if (hierr && hierr != HYPRE_ERROR_CONV) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
248:                                if (hierr) hypre__global_error = 0;);

250:   HYPREReplacePointer(jac->b,(PetscScalar*)sbv,bv);
251:   HYPREReplacePointer(jac->x,sxv,xv);
252:   VecRestoreArray(x,&xv);
253:   VecRestoreArrayRead(b,&bv);
254:   return(0);
255: }

259: static PetscErrorCode PCDestroy_HYPRE(PC pc)
260: {
261:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;

265:   if (jac->ij) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij));
266:   if (jac->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->b));
267:   if (jac->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->x));
268:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0]));
269:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1]));
270:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2]));
271:   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0]));
272:   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1]));
273:   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2]));
274:   if (jac->G) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->G));
275:   if (jac->C) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->C));
276:   if (jac->alpha_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->alpha_Poisson));
277:   if (jac->beta_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->beta_Poisson));
278:   if (jac->n_hmnull) {
279:     PetscInt i;

281:     for (i=0; i<jac->n_hmnull; i++) {
282:       PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->hmnull[i]));
283:     }
284:     PetscFree(jac->hmnull);
285:     PetscFree(jac->phmnull);
286:     VecDestroy(&jac->hmnull_constant);
287:   }
288:   if (jac->destroy) PetscStackCall("HYPRE_DestroyXXX",(*jac->destroy)(jac->hsolver););
289:   PetscFree(jac->hypre_type);
290:   if (jac->comm_hypre != MPI_COMM_NULL) { MPI_Comm_free(&(jac->comm_hypre));}
291:   PetscFree(pc->data);

293:   PetscObjectChangeTypeName((PetscObject)pc,0);
294:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);
295:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);
296:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetCoordinates_C",NULL);
297:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",NULL);
298:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",NULL);
299:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetConstantEdgeVectors_C",NULL);
300:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetAlphaPoissonMatrix_C",NULL);
301:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetBetaPoissonMatrix_C",NULL);
302:   return(0);
303: }

305: /* --------------------------------------------------------------------------------------------*/
308: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PetscOptions *PetscOptionsObject,PC pc)
309: {
310:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
312:   PetscBool      flag;

315:   PetscOptionsHead(PetscOptionsObject,"HYPRE Pilut Options");
316:   PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
317:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter));
318:   PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
319:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol));
320:   PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
321:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize));
322:   PetscOptionsTail();
323:   return(0);
324: }

328: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
329: {
330:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
332:   PetscBool      iascii;

335:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
336:   if (iascii) {
337:     PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut preconditioning\n");
338:     if (jac->maxiter != PETSC_DEFAULT) {
339:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: maximum number of iterations %d\n",jac->maxiter);
340:     } else {
341:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default maximum number of iterations \n");
342:     }
343:     if (jac->tol != PETSC_DEFAULT) {
344:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: drop tolerance %g\n",(double)jac->tol);
345:     } else {
346:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default drop tolerance \n");
347:     }
348:     if (jac->factorrowsize != PETSC_DEFAULT) {
349:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: factor row size %d\n",jac->factorrowsize);
350:     } else {
351:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default factor row size \n");
352:     }
353:   }
354:   return(0);
355: }

357: /* --------------------------------------------------------------------------------------------*/

361: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
362: {
363:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
364:   PetscErrorCode     ierr;
365:   HYPRE_ParCSRMatrix hmat;
366:   PetscScalar        *xv;
367:   const PetscScalar  *bv;
368:   HYPRE_ParVector    jbv,jxv;
369:   PetscScalar        *sbv,*sxv;
370:   PetscInt           hierr;

373:   PetscCitationsRegister(hypreCitation,&cite);
374:   VecSet(x,0.0);
375:   VecGetArrayRead(b,&bv);
376:   VecGetArray(x,&xv);
377:   HYPREReplacePointer(jac->b,(PetscScalar*)bv,sbv);
378:   HYPREReplacePointer(jac->x,xv,sxv);

380:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
381:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
382:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));

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

389:   HYPREReplacePointer(jac->b,sbv,bv);
390:   HYPREReplacePointer(jac->x,sxv,xv);
391:   VecRestoreArray(x,&xv);
392:   VecRestoreArrayRead(b,&bv);
393:   return(0);
394: }

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

399: static const char *HYPREBoomerAMGCycleType[]   = {"","V","W"};
400: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
401: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
402: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
403: static const char *HYPREBoomerAMGRelaxType[]   = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi",
404:                                                   "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi",
405:                                                   "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination",
406:                                                   "" /* 10 */, "" /* 11 */, "" /* 12 */, "" /* 13 */, "" /* 14 */,
407:                                                   "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"};
408: static const char *HYPREBoomerAMGInterpType[]  = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
409:                                                   "ext+i-cc", "standard", "standard-wts", "", "", "FF", "FF1"};
412: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PetscOptions *PetscOptionsObject,PC pc)
413: {
414:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
416:   PetscInt       n,indx,level;
417:   PetscBool      flg, tmp_truth;
418:   double         tmpdbl, twodbl[2];

421:   PetscOptionsHead(PetscOptionsObject,"HYPRE BoomerAMG Options");
422:   PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);
423:   if (flg) {
424:     jac->cycletype = indx+1;
425:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
426:   }
427:   PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
428:   if (flg) {
429:     if (jac->maxlevels < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
430:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
431:   }
432:   PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
433:   if (flg) {
434:     if (jac->maxiter < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
435:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
436:   }
437:   PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);
438:   if (flg) {
439:     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);
440:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
441:   }

443:   PetscOptionsScalar("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);
444:   if (flg) {
445:     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);
446:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
447:   }

449:   PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);
450:   if (flg) {
451:     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);
452:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
453:   }

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

459:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
460:   }


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

467:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
468:   }


471:   PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
472:   if (flg) {
473:     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);
474:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
475:   }
476:   PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
477:   if (flg) {
478:     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);
479:     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);
480:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
481:   }

483:   /* Grid sweeps */
484:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);
485:   if (flg) {
486:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx));
487:     /* modify the jac structure so we can view the updated options with PC_View */
488:     jac->gridsweeps[0] = indx;
489:     jac->gridsweeps[1] = indx;
490:     /*defaults coarse to 1 */
491:     jac->gridsweeps[2] = 1;
492:   }

494:   PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen","Use a nodal based coarsening 1-6","HYPRE_BoomerAMGSetNodal",jac->nodal_coarsening, &jac->nodal_coarsening,&flg);
495:   if (flg) {
496:     PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,jac->nodal_coarsening));
497:   }

499:   PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_variant","Variant of algorithm 1-4","HYPRE_BoomerAMGSetInterpVecVariant",jac->vec_interp_variant, &jac->vec_interp_variant,&flg);
500:   if (flg) {
501:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecVariant,(jac->hsolver,jac->vec_interp_variant));
502:   }

504:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);
505:   if (flg) {
506:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1));
507:     jac->gridsweeps[0] = indx;
508:   }
509:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);
510:   if (flg) {
511:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2));
512:     jac->gridsweeps[1] = indx;
513:   }
514:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);
515:   if (flg) {
516:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3));
517:     jac->gridsweeps[2] = indx;
518:   }

520:   /* Relax type */
521:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
522:   if (flg) {
523:     jac->relaxtype[0] = jac->relaxtype[1]  = indx;
524:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx));
525:     /* by default, coarse type set to 9 */
526:     jac->relaxtype[2] = 9;

528:   }
529:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
530:   if (flg) {
531:     jac->relaxtype[0] = indx;
532:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1));
533:   }
534:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
535:   if (flg) {
536:     jac->relaxtype[1] = indx;
537:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2));
538:   }
539:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);
540:   if (flg) {
541:     jac->relaxtype[2] = indx;
542:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3));
543:   }

545:   /* Relaxation Weight */
546:   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);
547:   if (flg) {
548:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl));
549:     jac->relaxweight = tmpdbl;
550:   }

552:   n         = 2;
553:   twodbl[0] = twodbl[1] = 1.0;
554:   PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
555:   if (flg) {
556:     if (n == 2) {
557:       indx =  (int)PetscAbsReal(twodbl[1]);
558:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx));
559:     } 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);
560:   }

562:   /* Outer relaxation Weight */
563:   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);
564:   if (flg) {
565:     PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl));
566:     jac->outerrelaxweight = tmpdbl;
567:   }

569:   n         = 2;
570:   twodbl[0] = twodbl[1] = 1.0;
571:   PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
572:   if (flg) {
573:     if (n == 2) {
574:       indx =  (int)PetscAbsReal(twodbl[1]);
575:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx));
576:     } 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);
577:   }

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

582:   if (flg && tmp_truth) {
583:     jac->relaxorder = 0;
584:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
585:   }
586:   PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);
587:   if (flg) {
588:     jac->measuretype = indx;
589:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
590:   }
591:   /* update list length 3/07 */
592:   PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);
593:   if (flg) {
594:     jac->coarsentype = indx;
595:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
596:   }

598:   /* new 3/07 */
599:   PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);
600:   if (flg) {
601:     jac->interptype = indx;
602:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
603:   }

605:   PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
606:   if (flg) {
607:     level = 3;
608:     PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);

610:     jac->printstatistics = PETSC_TRUE;
611:     PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level));
612:   }

614:   PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);
615:   if (flg) {
616:     level = 3;
617:     PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);

619:     jac->printstatistics = PETSC_TRUE;
620:     PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level));
621:   }

623:   PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);
624:   if (flg && tmp_truth) {
625:     PetscInt tmp_int;
626:     PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);
627:     if (flg) jac->nodal_relax_levels = tmp_int;
628:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6));
629:     PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1));
630:     PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0));
631:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels));
632:   }

634:   PetscOptionsTail();
635:   return(0);
636: }

640: 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)
641: {
642:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
644:   PetscInt       oits;

647:   PetscCitationsRegister(hypreCitation,&cite);
648:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter));
649:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol));
650:   jac->applyrichardson = PETSC_TRUE;
651:   PCApply_HYPRE(pc,b,y);
652:   jac->applyrichardson = PETSC_FALSE;
653:   PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
654:   *outits = oits;
655:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
656:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
657:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
658:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
659:   return(0);
660: }


665: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
666: {
667:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
669:   PetscBool      iascii;

672:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
673:   if (iascii) {
674:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG preconditioning\n");
675:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);
676:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of levels %d\n",jac->maxlevels);
677:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of iterations PER hypre call %d\n",jac->maxiter);
678:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Convergence tolerance PER hypre call %g\n",(double)jac->tol);
679:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Threshold for strong coupling %g\n",(double)jac->strongthreshold);
680:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation truncation factor %g\n",(double)jac->truncfactor);
681:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation: max elements per row %d\n",jac->pmax);
682:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Number of levels of aggressive coarsening %d\n",jac->agg_nl);
683:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Number of paths for aggressive coarsening %d\n",jac->agg_num_paths);

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

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

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

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

698:     if (jac->relaxorder) {
699:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Using CF-relaxation\n");
700:     } else {
701:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Not using CF-relaxation\n");
702:     }
703:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Measure type        %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
704:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Coarsen type        %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
705:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation type  %s\n",HYPREBoomerAMGInterpType[jac->interptype]);
706:     if (jac->nodal_coarsening) {
707:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Using nodal coarsening (with HYPRE_BOOMERAMGSetNodal() %D\n",jac->nodal_coarsening);
708:     }
709:     if (jac->vec_interp_variant) {
710:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: HYPRE_BoomerAMGSetInterpVecVariant() %D\n",jac->vec_interp_variant);
711:     }
712:     if (jac->nodal_relax) {
713:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Using nodal relaxation via Schwarz smoothing on levels %d\n",jac->nodal_relax_levels);
714:     }
715:   }
716:   return(0);
717: }

719: /* --------------------------------------------------------------------------------------------*/
722: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PetscOptions *PetscOptionsObject,PC pc)
723: {
724:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
726:   PetscInt       indx;
727:   PetscBool      flag;
728:   const char     *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};

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

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

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

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

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

748:   PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);
749:   if (flag) {
750:     jac->symt = indx;
751:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
752:   }

754:   PetscOptionsTail();
755:   return(0);
756: }

760: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
761: {
762:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
764:   PetscBool      iascii;
765:   const char     *symt = 0;;

768:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
769:   if (iascii) {
770:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails preconditioning\n");
771:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: nlevels %d\n",jac->nlevels);
772:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: threshold %g\n",(double)jac->threshhold);
773:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: filter %g\n",(double)jac->filter);
774:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: load balance %g\n",(double)jac->loadbal);
775:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: reuse nonzero structure %s\n",PetscBools[jac->ruse]);
776:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: print info to screen %s\n",PetscBools[jac->logging]);
777:     if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
778:     else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
779:     else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
780:     else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
781:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: %s\n",symt);
782:   }
783:   return(0);
784: }
785: /* --------------------------------------------------------------------------------------------*/
788: static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PetscOptions *PetscOptionsObject,PC pc)
789: {
790:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
792:   PetscInt       n;
793:   PetscBool      flag,flag2,flag3,flag4;

796:   PetscOptionsHead(PetscOptionsObject,"HYPRE AMS Options");
797:   PetscOptionsInt("-pc_hypre_ams_print_level","Debugging output level for AMS","None",jac->as_print,&jac->as_print,&flag);
798:   if (flag) PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
799:   PetscOptionsInt("-pc_hypre_ams_max_iter","Maximum number of AMS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
800:   if (flag) PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
801:   PetscOptionsInt("-pc_hypre_ams_cycle_type","Cycle type for AMS multigrid","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag);
802:   if (flag) PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
803:   PetscOptionsReal("-pc_hypre_ams_tol","Error tolerance for AMS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
804:   if (flag) PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
805:   PetscOptionsInt("-pc_hypre_ams_relax_type","Relaxation type for AMS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
806:   PetscOptionsInt("-pc_hypre_ams_relax_times","Number of relaxation steps for AMS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
807:   PetscOptionsReal("-pc_hypre_ams_relax_weight","Relaxation weight for AMS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
808:   PetscOptionsReal("-pc_hypre_ams_omega","SSOR coefficient for AMS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
809:   if (flag || flag2 || flag3 || flag4) {
810:     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
811:                                                                       jac->as_relax_times,
812:                                                                       jac->as_relax_weight,
813:                                                                       jac->as_omega));
814:   }
815:   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);
816:   n = 5;
817:   PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options","AMG options for vector Poisson","None",jac->as_amg_alpha_opts,&n,&flag2);
818:   if (flag || flag2) {
819:     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
820:                                                                      jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
821:                                                                      jac->as_amg_alpha_opts[2],       /* AMG relax_type */
822:                                                                      jac->as_amg_alpha_theta,
823:                                                                      jac->as_amg_alpha_opts[3],       /* AMG interp_type */
824:                                                                      jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
825:   }
826:   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);
827:   n = 5;
828:   PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options","AMG options for scalar Poisson solver","None",jac->as_amg_beta_opts,&n,&flag2);
829:   if (flag || flag2) {
830:     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
831:                                                                     jac->as_amg_beta_opts[1],       /* AMG agg_levels */
832:                                                                     jac->as_amg_beta_opts[2],       /* AMG relax_type */
833:                                                                     jac->as_amg_beta_theta,
834:                                                                     jac->as_amg_beta_opts[3],       /* AMG interp_type */
835:                                                                     jac->as_amg_beta_opts[4]));     /* AMG Pmax */
836:   }
837:   PetscOptionsTail();
838:   return(0);
839: }

843: static PetscErrorCode PCView_HYPRE_AMS(PC pc,PetscViewer viewer)
844: {
845:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
847:   PetscBool      iascii;

850:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
851:   if (iascii) {
852:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS preconditioning\n");
853:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace iterations per application %d\n",jac->as_max_iter);
854:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace cycle type %d\n",jac->ams_cycle_type);
855:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace iteration tolerance %g\n",jac->as_tol);
856:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother type %d\n",jac->as_relax_type);
857:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: number of smoothing steps %d\n",jac->as_relax_times);
858:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother weight %g\n",jac->as_relax_weight);
859:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother omega %g\n",jac->as_omega);
860:     if (jac->alpha_Poisson) {
861:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: vector Poisson solver (passed in by user)\n");
862:     } else {
863:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: vector Poisson solver (computed) \n");
864:     }
865:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);
866:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
867:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);
868:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);
869:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
870:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);
871:     if (!jac->ams_beta_is_zero) {
872:       if (jac->beta_Poisson) {
873:         PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: scalar Poisson solver (passed in by user)\n");
874:       } else {
875:         PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: scalar Poisson solver (computed) \n");
876:       }
877:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);
878:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
879:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);
880:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);
881:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
882:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);
883:     }
884:   }
885:   return(0);
886: }

890: static PetscErrorCode PCSetFromOptions_HYPRE_ADS(PetscOptions *PetscOptionsObject,PC pc)
891: {
892:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
894:   PetscInt       n;
895:   PetscBool      flag,flag2,flag3,flag4;

898:   PetscOptionsHead(PetscOptionsObject,"HYPRE ADS Options");
899:   PetscOptionsInt("-pc_hypre_ads_print_level","Debugging output level for ADS","None",jac->as_print,&jac->as_print,&flag);
900:   if (flag) PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
901:   PetscOptionsInt("-pc_hypre_ads_max_iter","Maximum number of ADS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
902:   if (flag) PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
903:   PetscOptionsInt("-pc_hypre_ads_cycle_type","Cycle type for ADS multigrid","None",jac->ads_cycle_type,&jac->ads_cycle_type,&flag);
904:   if (flag) PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ads_cycle_type));
905:   PetscOptionsReal("-pc_hypre_ads_tol","Error tolerance for ADS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
906:   if (flag) PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
907:   PetscOptionsInt("-pc_hypre_ads_relax_type","Relaxation type for ADS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
908:   PetscOptionsInt("-pc_hypre_ads_relax_times","Number of relaxation steps for ADS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
909:   PetscOptionsReal("-pc_hypre_ads_relax_weight","Relaxation weight for ADS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
910:   PetscOptionsReal("-pc_hypre_ads_omega","SSOR coefficient for ADS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
911:   if (flag || flag2 || flag3 || flag4) {
912:     PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
913:                                                                       jac->as_relax_times,
914:                                                                       jac->as_relax_weight,
915:                                                                       jac->as_omega));
916:   }
917:   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);
918:   n = 5;
919:   PetscOptionsIntArray("-pc_hypre_ads_ams_options","AMG options for AMS solver inside ADS","None",jac->as_amg_alpha_opts,&n,&flag2);
920:   PetscOptionsInt("-pc_hypre_ads_ams_cycle_type","Cycle type for AMS solver inside ADS","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag3);
921:   if (flag || flag2 || flag3) {
922:     PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type,             /* AMS cycle type */
923:                                                                 jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
924:                                                                 jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
925:                                                                 jac->as_amg_alpha_opts[2],       /* AMG relax_type */
926:                                                                 jac->as_amg_alpha_theta,
927:                                                                 jac->as_amg_alpha_opts[3],       /* AMG interp_type */
928:                                                                 jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
929:   }
930:   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);
931:   n = 5;
932:   PetscOptionsIntArray("-pc_hypre_ads_amg_options","AMG options for vector AMG solver inside ADS","None",jac->as_amg_beta_opts,&n,&flag2);
933:   if (flag || flag2) {
934:     PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
935:                                                                 jac->as_amg_beta_opts[1],       /* AMG agg_levels */
936:                                                                 jac->as_amg_beta_opts[2],       /* AMG relax_type */
937:                                                                 jac->as_amg_beta_theta,
938:                                                                 jac->as_amg_beta_opts[3],       /* AMG interp_type */
939:                                                                 jac->as_amg_beta_opts[4]));     /* AMG Pmax */
940:   }
941:   PetscOptionsTail();
942:   return(0);
943: }

947: static PetscErrorCode PCView_HYPRE_ADS(PC pc,PetscViewer viewer)
948: {
949:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
951:   PetscBool      iascii;

954:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
955:   if (iascii) {
956:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS preconditioning\n");
957:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: subspace iterations per application %d\n",jac->as_max_iter);
958:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: subspace cycle type %d\n",jac->ads_cycle_type);
959:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: subspace iteration tolerance %g\n",jac->as_tol);
960:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: smoother type %d\n",jac->as_relax_type);
961:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: number of smoothing steps %d\n",jac->as_relax_times);
962:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: smoother weight %g\n",jac->as_relax_weight);
963:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: smoother omega %g\n",jac->as_omega);
964:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: AMS solver\n");
965:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     subspace cycle type %d\n",jac->ams_cycle_type);
966:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);
967:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
968:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);
969:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);
970:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
971:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);
972:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: vector Poisson solver\n");
973:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);
974:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
975:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);
976:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);
977:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
978:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);
979:   }
980:   return(0);
981: }

985: static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE(PC pc, Mat G)
986: {
987:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
988:   HYPRE_ParCSRMatrix parcsr_G;
989:   PetscErrorCode     ierr;

992:   /* throw away any discrete gradient if already set */
993:   if (jac->G) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->G));
994:   MatHYPRE_IJMatrixCreate(G,&jac->G);
995:   MatHYPRE_IJMatrixCopy(G,jac->G);
996:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->G,(void**)(&parcsr_G)));
997:   PetscStackCall("Hypre set gradient",(*jac->setdgrad)(jac->hsolver,parcsr_G););
998:   return(0);
999: }

1003: /*@
1004:  PCHYPRESetDiscreteGradient - Set discrete gradient matrix

1006:    Collective on PC

1008:    Input Parameters:
1009: +  pc - the preconditioning context
1010: -  G - the discrete gradient

1012:    Level: intermediate

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

1017: .seealso:
1018: @*/
1019: PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G)
1020: {

1027:   PetscTryMethod(pc,"PCHYPRESetDiscreteGradient_C",(PC,Mat),(pc,G));
1028:   return(0);
1029: }

1033: static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C)
1034: {
1035:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1036:   HYPRE_ParCSRMatrix parcsr_C;
1037:   PetscErrorCode     ierr;

1040:   /* throw away any discrete curl if already set */
1041:   if (jac->C) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->C));
1042:   MatHYPRE_IJMatrixCreate(C,&jac->C);
1043:   MatHYPRE_IJMatrixCopy(C,jac->C);
1044:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->C,(void**)(&parcsr_C)));
1045:   PetscStackCall("Hypre set curl",(*jac->setdcurl)(jac->hsolver,parcsr_C););
1046:   return(0);
1047: }

1051: /*@
1052:  PCHYPRESetDiscreteCurl - Set discrete curl matrix

1054:    Collective on PC

1056:    Input Parameters:
1057: +  pc - the preconditioning context
1058: -  C - the discrete curl

1060:    Level: intermediate

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

1065: .seealso:
1066: @*/
1067: PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C)
1068: {

1075:   PetscTryMethod(pc,"PCHYPRESetDiscreteCurl_C",(PC,Mat),(pc,C));
1076:   return(0);
1077: }

1081: static PetscErrorCode PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS(PC pc, Mat A)
1082: {
1083:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1084:   HYPRE_ParCSRMatrix parcsr_alpha_Poisson;
1085:   PetscErrorCode     ierr;

1088:   /* throw away any matrix if already set */
1089:   if (jac->alpha_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->alpha_Poisson));
1090:   MatHYPRE_IJMatrixCreate(A,&jac->alpha_Poisson);
1091:   MatHYPRE_IJMatrixCopy(A,jac->alpha_Poisson);
1092:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->alpha_Poisson,(void**)(&parcsr_alpha_Poisson)));
1093:   PetscStackCallStandard(HYPRE_AMSSetAlphaPoissonMatrix,(jac->hsolver,parcsr_alpha_Poisson));
1094:   return(0);
1095: }

1099: /*@
1100:  PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix

1102:    Collective on PC

1104:    Input Parameters:
1105: +  pc - the preconditioning context
1106: -  A - the matrix

1108:    Level: intermediate

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

1112: .seealso:
1113: @*/
1114: PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A)
1115: {

1122:   PetscTryMethod(pc,"PCHYPRESetAlphaPoissonMatrix_C",(PC,Mat),(pc,A));
1123:   return(0);
1124: }

1128: static PetscErrorCode PCHYPRESetBetaPoissonMatrix_HYPRE_AMS(PC pc, Mat A)
1129: {
1130:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1131:   HYPRE_ParCSRMatrix parcsr_beta_Poisson;
1132:   PetscErrorCode     ierr;

1135:   if (!A) {
1136:     PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,NULL));
1137:     jac->ams_beta_is_zero = PETSC_TRUE;
1138:     return(0);
1139:   }
1140:   jac->ams_beta_is_zero = PETSC_FALSE;
1141:   /* throw away any matrix if already set */
1142:   if (jac->beta_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->beta_Poisson));
1143:   MatHYPRE_IJMatrixCreate(A,&jac->beta_Poisson);
1144:   MatHYPRE_IJMatrixCopy(A,jac->beta_Poisson);
1145:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->beta_Poisson,(void**)(&parcsr_beta_Poisson)));
1146:   PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,parcsr_beta_Poisson));
1147:   return(0);
1148: }

1152: /*@
1153:  PCHYPRESetBetaPoissonMatrix - Set Poisson matrix

1155:    Collective on PC

1157:    Input Parameters:
1158: +  pc - the preconditioning context
1159: -  A - the matrix

1161:    Level: intermediate

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

1166: .seealso:
1167: @*/
1168: PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A)
1169: {

1174:   if (A) {
1177:   }
1178:   PetscTryMethod(pc,"PCHYPRESetBetaPoissonMatrix_C",(PC,Mat),(pc,A));
1179:   return(0);
1180: }

1184: static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE_AMS(PC pc,Vec ozz, Vec zoz, Vec zzo)
1185: {
1186:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1187:   HYPRE_ParVector    par_ozz,par_zoz,par_zzo;
1188:   PetscInt           dim;
1189:   PetscErrorCode     ierr;

1192:   /* throw away any vector if already set */
1193:   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0]));
1194:   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1]));
1195:   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2]));
1196:   jac->constants[0] = NULL;
1197:   jac->constants[1] = NULL;
1198:   jac->constants[2] = NULL;
1199:   VecHYPRE_IJVectorCreate(ozz,&jac->constants[0]);
1200:   VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);
1201:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[0],(void**)(&par_ozz)));
1202:   VecHYPRE_IJVectorCreate(zoz,&jac->constants[1]);
1203:   VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);
1204:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[1],(void**)(&par_zoz)));
1205:   dim = 2;
1206:   if (zzo) {
1207:     VecHYPRE_IJVectorCreate(zzo,&jac->constants[2]);
1208:     VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);
1209:     PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[2],(void**)(&par_zzo)));
1210:     dim++;
1211:   }
1212:   PetscStackCallStandard(HYPRE_AMSSetEdgeConstantVectors,(jac->hsolver,par_ozz,par_zoz,par_zzo));
1213:   PetscStackCallStandard(HYPRE_AMSSetDimension,(jac->hsolver,dim));
1214:   return(0);
1215: }

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

1222:    Collective on PC

1224:    Input Parameters:
1225: +  pc - the preconditioning context
1226: -  ozz - vector representing (1,0,0) (or (1,0) in 2D)
1227: -  zoz - vector representing (0,1,0) (or (0,1) in 2D)
1228: -  zzo - vector representing (0,0,1) (use NULL in 2D)

1230:    Level: intermediate

1232:    Notes:

1234: .seealso:
1235: @*/
1236: PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo)
1237: {

1248:   PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));
1249:   return(0);
1250: }

1254: static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1255: {
1256:   PC_HYPRE        *jac = (PC_HYPRE*)pc->data;
1257:   Vec             tv;
1258:   HYPRE_ParVector par_coords[3];
1259:   PetscInt        i;
1260:   PetscErrorCode  ierr;

1263:   /* throw away any coordinate vector if already set */
1264:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0]));
1265:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1]));
1266:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2]));
1267:   /* set problem's dimension */
1268:   if (jac->setdim) {
1269:     PetscStackCall("Hypre set dim",(*jac->setdim)(jac->hsolver,dim););
1270:   }
1271:   /* compute IJ vector for coordinates */
1272:   VecCreate(PetscObjectComm((PetscObject)pc),&tv);
1273:   VecSetType(tv,VECSTANDARD);
1274:   VecSetSizes(tv,nloc,PETSC_DECIDE);
1275:   for (i=0;i<dim;i++) {
1276:     PetscScalar *array;
1277:     PetscInt    j;

1279:     VecHYPRE_IJVectorCreate(tv,&jac->coords[i]);
1280:     VecGetArray(tv,&array);
1281:     for (j=0;j<nloc;j++) {
1282:       array[j] = coords[j*dim+i];
1283:     }
1284:     PetscStackCallStandard(HYPRE_IJVectorSetValues,(jac->coords[i],nloc,NULL,array));
1285:     PetscStackCallStandard(HYPRE_IJVectorAssemble,(jac->coords[i]));
1286:     VecRestoreArray(tv,&array);
1287:   }
1288:   VecDestroy(&tv);
1289:   /* pass parCSR vectors to AMS solver */
1290:   par_coords[0] = NULL;
1291:   par_coords[1] = NULL;
1292:   par_coords[2] = NULL;
1293:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0],(void**)(&par_coords[0])));
1294:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1],(void**)(&par_coords[1])));
1295:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2],(void**)(&par_coords[2])));
1296:   PetscStackCall("Hypre set coords",(*jac->setcoord)(jac->hsolver,par_coords[0],par_coords[1],par_coords[2]););
1297:   return(0);
1298: }

1300: /* ---------------------------------------------------------------------------------*/

1304: static PetscErrorCode  PCHYPREGetType_HYPRE(PC pc,const char *name[])
1305: {
1306:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

1309:   *name = jac->hypre_type;
1310:   return(0);
1311: }

1315: static PetscErrorCode  PCHYPRESetType_HYPRE(PC pc,const char name[])
1316: {
1317:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1319:   PetscBool      flag;

1322:   if (jac->hypre_type) {
1323:     PetscStrcmp(jac->hypre_type,name,&flag);
1324:     if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set");
1325:     return(0);
1326:   } else {
1327:     PetscStrallocpy(name, &jac->hypre_type);
1328:   }

1330:   jac->maxiter         = PETSC_DEFAULT;
1331:   jac->tol             = PETSC_DEFAULT;
1332:   jac->printstatistics = PetscLogPrintInfo;

1334:   PetscStrcmp("pilut",jac->hypre_type,&flag);
1335:   if (flag) {
1336:     PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver));
1337:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
1338:     pc->ops->view           = PCView_HYPRE_Pilut;
1339:     jac->destroy            = HYPRE_ParCSRPilutDestroy;
1340:     jac->setup              = HYPRE_ParCSRPilutSetup;
1341:     jac->solve              = HYPRE_ParCSRPilutSolve;
1342:     jac->factorrowsize      = PETSC_DEFAULT;
1343:     return(0);
1344:   }
1345:   PetscStrcmp("parasails",jac->hypre_type,&flag);
1346:   if (flag) {
1347:     PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver));
1348:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
1349:     pc->ops->view           = PCView_HYPRE_ParaSails;
1350:     jac->destroy            = HYPRE_ParaSailsDestroy;
1351:     jac->setup              = HYPRE_ParaSailsSetup;
1352:     jac->solve              = HYPRE_ParaSailsSolve;
1353:     /* initialize */
1354:     jac->nlevels    = 1;
1355:     jac->threshhold = .1;
1356:     jac->filter     = .1;
1357:     jac->loadbal    = 0;
1358:     if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE;
1359:     else jac->logging = (int) PETSC_FALSE;

1361:     jac->ruse = (int) PETSC_FALSE;
1362:     jac->symt = 0;
1363:     PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
1364:     PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
1365:     PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
1366:     PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
1367:     PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
1368:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
1369:     return(0);
1370:   }
1371:   PetscStrcmp("boomeramg",jac->hypre_type,&flag);
1372:   if (flag) {
1373:     HYPRE_BoomerAMGCreate(&jac->hsolver);
1374:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_BoomerAMG;
1375:     pc->ops->view            = PCView_HYPRE_BoomerAMG;
1376:     pc->ops->applytranspose  = PCApplyTranspose_HYPRE_BoomerAMG;
1377:     pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
1378:     jac->destroy             = HYPRE_BoomerAMGDestroy;
1379:     jac->setup               = HYPRE_BoomerAMGSetup;
1380:     jac->solve               = HYPRE_BoomerAMGSolve;
1381:     jac->applyrichardson     = PETSC_FALSE;
1382:     /* these defaults match the hypre defaults */
1383:     jac->cycletype        = 1;
1384:     jac->maxlevels        = 25;
1385:     jac->maxiter          = 1;
1386:     jac->tol              = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
1387:     jac->truncfactor      = 0.0;
1388:     jac->strongthreshold  = .25;
1389:     jac->maxrowsum        = .9;
1390:     jac->coarsentype      = 6;
1391:     jac->measuretype      = 0;
1392:     jac->gridsweeps[0]    = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
1393:     jac->relaxtype[0]     = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */
1394:     jac->relaxtype[2]     = 9; /*G.E. */
1395:     jac->relaxweight      = 1.0;
1396:     jac->outerrelaxweight = 1.0;
1397:     jac->relaxorder       = 1;
1398:     jac->interptype       = 0;
1399:     jac->agg_nl           = 0;
1400:     jac->pmax             = 0;
1401:     jac->truncfactor      = 0.0;
1402:     jac->agg_num_paths    = 1;

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

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

1566: /*
1567:     It only gets here if the HYPRE type has not been set before the call to
1568:    ...SetFromOptions() which actually is most of the time
1569: */
1572: static PetscErrorCode PCSetFromOptions_HYPRE(PetscOptions *PetscOptionsObject,PC pc)
1573: {
1575:   PetscInt       indx;
1576:   const char     *type[] = {"pilut","parasails","boomeramg","ams","ads"};
1577:   PetscBool      flg;

1580:   PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");
1581:   PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"boomeramg",&indx,&flg);
1582:   if (flg) {
1583:     PCHYPRESetType_HYPRE(pc,type[indx]);
1584:   } else {
1585:     PCHYPRESetType_HYPRE(pc,"boomeramg");
1586:   }
1587:   if (pc->ops->setfromoptions) {
1588:     pc->ops->setfromoptions(PetscOptionsObject,pc);
1589:   }
1590:   PetscOptionsTail();
1591:   return(0);
1592: }

1596: /*@C
1597:      PCHYPRESetType - Sets which hypre preconditioner you wish to use

1599:    Input Parameters:
1600: +     pc - the preconditioner context
1601: -     name - either  pilut, parasails, boomeramg, ams, ads

1603:    Options Database Keys:
1604:    -pc_hypre_type - One of pilut, parasails, boomeramg, ams, ads

1606:    Level: intermediate

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

1611: @*/
1612: PetscErrorCode  PCHYPRESetType(PC pc,const char name[])
1613: {

1619:   PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));
1620:   return(0);
1621: }

1625: /*@C
1626:      PCHYPREGetType - Gets which hypre preconditioner you are using

1628:    Input Parameter:
1629: .     pc - the preconditioner context

1631:    Output Parameter:
1632: .     name - either  pilut, parasails, boomeramg, ams, ads

1634:    Level: intermediate

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

1639: @*/
1640: PetscErrorCode  PCHYPREGetType(PC pc,const char *name[])
1641: {

1647:   PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));
1648:   return(0);
1649: }

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

1654:    Options Database Keys:
1655: +   -pc_hypre_type - One of pilut, parasails, boomeramg, ams, ads
1656: -   Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX
1657:           preconditioner

1659:    Level: intermediate

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

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

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

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

1683:           MatSetNearNullSpace() - if you provide a near null space to your matrix it is ignored by hypre UNLESS you also use
1684:           the two options:
1685: +   -pc_hypre_boomeramg_nodal_coarsen <n> - where n is from 1 to 6 (see HYPRE_BOOMERAMGSetNodal())
1686: -   -pc_hypre_boomeramg_vec_interp_variant <v> where v is from 1 to 4 (see HYPRE_BoomerAMGSetInterpVecVariant())

1688:           Depending on the linear system you may see the same or different convergence depending on the values you use.

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

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

1695: M*/

1699: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
1700: {
1701:   PC_HYPRE       *jac;

1705:   PetscNewLog(pc,&jac);

1707:   pc->data                = jac;
1708:   pc->ops->destroy        = PCDestroy_HYPRE;
1709:   pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
1710:   pc->ops->setup          = PCSetUp_HYPRE;
1711:   pc->ops->apply          = PCApply_HYPRE;
1712:   jac->comm_hypre         = MPI_COMM_NULL;
1713:   jac->hypre_type         = NULL;
1714:   jac->coords[0]          = NULL;
1715:   jac->coords[1]          = NULL;
1716:   jac->coords[2]          = NULL;
1717:   jac->constants[0]       = NULL;
1718:   jac->constants[1]       = NULL;
1719:   jac->constants[2]       = NULL;
1720:   jac->G                  = NULL;
1721:   jac->C                  = NULL;
1722:   jac->alpha_Poisson      = NULL;
1723:   jac->beta_Poisson       = NULL;
1724:   jac->setdim             = NULL;
1725:   jac->hmnull             = NULL;
1726:   jac->n_hmnull           = 0;
1727:   /* duplicate communicator for hypre */
1728:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1729:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);
1730:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);
1731:   return(0);
1732: }

1734: /* ---------------------------------------------------------------------------------------------------------------------------------*/

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

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

1743:   /* keep copy of PFMG options used so may view them */
1744:   PetscInt its;
1745:   double   tol;
1746:   PetscInt relax_type;
1747:   PetscInt rap_type;
1748:   PetscInt num_pre_relax,num_post_relax;
1749:   PetscInt max_levels;
1750: } PC_PFMG;

1754: PetscErrorCode PCDestroy_PFMG(PC pc)
1755: {
1757:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

1760:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1761:   MPI_Comm_free(&ex->hcomm);
1762:   PetscFree(pc->data);
1763:   return(0);
1764: }

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

1771: PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
1772: {
1774:   PetscBool      iascii;
1775:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

1778:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1779:   if (iascii) {
1780:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG preconditioning\n");
1781:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: max iterations %d\n",ex->its);
1782:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: tolerance %g\n",ex->tol);
1783:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
1784:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: RAP type %s\n",PFMGRAPType[ex->rap_type]);
1785:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
1786:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: max levels %d\n",ex->max_levels);
1787:   }
1788:   return(0);
1789: }


1794: PetscErrorCode PCSetFromOptions_PFMG(PetscOptions *PetscOptionsObject,PC pc)
1795: {
1797:   PC_PFMG        *ex = (PC_PFMG*) pc->data;
1798:   PetscBool      flg = PETSC_FALSE;

1801:   PetscOptionsHead(PetscOptionsObject,"PFMG options");
1802:   PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);
1803:   if (flg) {
1804:     int level=3;
1805:     PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level));
1806:   }
1807:   PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);
1808:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its));
1809:   PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
1810:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax));
1811:   PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
1812:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax));

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

1817:   PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);
1818:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol));
1819:   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);
1820:   PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type));
1821:   PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);
1822:   PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type));
1823:   PetscOptionsTail();
1824:   return(0);
1825: }

1829: PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
1830: {
1831:   PetscErrorCode    ierr;
1832:   PC_PFMG           *ex = (PC_PFMG*) pc->data;
1833:   PetscScalar       *yy;
1834:   const PetscScalar *xx;
1835:   PetscInt          ilower[3],iupper[3];
1836:   Mat_HYPREStruct   *mx = (Mat_HYPREStruct*)(pc->pmat->data);

1839:   PetscCitationsRegister(hypreCitation,&cite);
1840:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
1841:   iupper[0] += ilower[0] - 1;
1842:   iupper[1] += ilower[1] - 1;
1843:   iupper[2] += ilower[2] - 1;

1845:   /* copy x values over to hypre */
1846:   PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
1847:   VecGetArrayRead(x,&xx);
1848:   PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx));
1849:   VecRestoreArrayRead(x,&xx);
1850:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
1851:   PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx));

1853:   /* copy solution values back to PETSc */
1854:   VecGetArray(y,&yy);
1855:   PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy));
1856:   VecRestoreArray(y,&yy);
1857:   return(0);
1858: }

1862: 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)
1863: {
1864:   PC_PFMG        *jac = (PC_PFMG*)pc->data;
1866:   PetscInt       oits;

1869:   PetscCitationsRegister(hypreCitation,&cite);
1870:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its));
1871:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol));

1873:   PCApply_PFMG(pc,b,y);
1874:   PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
1875:   *outits = oits;
1876:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1877:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
1878:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol));
1879:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its));
1880:   return(0);
1881: }


1886: PetscErrorCode PCSetUp_PFMG(PC pc)
1887: {
1888:   PetscErrorCode  ierr;
1889:   PC_PFMG         *ex = (PC_PFMG*) pc->data;
1890:   Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
1891:   PetscBool       flg;

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

1897:   /* create the hypre solver object and set its information */
1898:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1899:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1900:   PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
1901:   PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver));
1902:   return(0);
1903: }


1906: /*MC
1907:      PCPFMG - the hypre PFMG multigrid solver

1909:    Level: advanced

1911:    Options Database:
1912: + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner
1913: . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
1914: . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid
1915: . -pc_pfmg_tol <tol> tolerance of PFMG
1916: . -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
1917: - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin

1919:    Notes:  This is for CELL-centered descretizations

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

1924: .seealso:  PCMG, MATHYPRESTRUCT
1925: M*/

1929: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
1930: {
1932:   PC_PFMG        *ex;

1935:   PetscNew(&ex); \
1936:   pc->data = ex;

1938:   ex->its            = 1;
1939:   ex->tol            = 1.e-8;
1940:   ex->relax_type     = 1;
1941:   ex->rap_type       = 0;
1942:   ex->num_pre_relax  = 1;
1943:   ex->num_post_relax = 1;
1944:   ex->max_levels     = 0;

1946:   pc->ops->setfromoptions  = PCSetFromOptions_PFMG;
1947:   pc->ops->view            = PCView_PFMG;
1948:   pc->ops->destroy         = PCDestroy_PFMG;
1949:   pc->ops->apply           = PCApply_PFMG;
1950:   pc->ops->applyrichardson = PCApplyRichardson_PFMG;
1951:   pc->ops->setup           = PCSetUp_PFMG;

1953:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
1954:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1955:   return(0);
1956: }

1958: /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/

1960: /* we know we are working with a HYPRE_SStructMatrix */
1961: typedef struct {
1962:   MPI_Comm            hcomm;       /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
1963:   HYPRE_SStructSolver ss_solver;

1965:   /* keep copy of SYSPFMG options used so may view them */
1966:   PetscInt its;
1967:   double   tol;
1968:   PetscInt relax_type;
1969:   PetscInt num_pre_relax,num_post_relax;
1970: } PC_SysPFMG;

1974: PetscErrorCode PCDestroy_SysPFMG(PC pc)
1975: {
1977:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

1980:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
1981:   MPI_Comm_free(&ex->hcomm);
1982:   PetscFree(pc->data);
1983:   return(0);
1984: }

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

1990: PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
1991: {
1993:   PetscBool      iascii;
1994:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

1997:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1998:   if (iascii) {
1999:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG preconditioning\n");
2000:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: max iterations %d\n",ex->its);
2001:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: tolerance %g\n",ex->tol);
2002:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
2003:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
2004:   }
2005:   return(0);
2006: }


2011: PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptions *PetscOptionsObject,PC pc)
2012: {
2014:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;
2015:   PetscBool      flg = PETSC_FALSE;

2018:   PetscOptionsHead(PetscOptionsObject,"SysPFMG options");
2019:   PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);
2020:   if (flg) {
2021:     int level=3;
2022:     PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level));
2023:   }
2024:   PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);
2025:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its));
2026:   PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
2027:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax));
2028:   PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
2029:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax));

2031:   PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);
2032:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol));
2033:   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);
2034:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type));
2035:   PetscOptionsTail();
2036:   return(0);
2037: }

2041: PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
2042: {
2043:   PetscErrorCode    ierr;
2044:   PC_SysPFMG        *ex = (PC_SysPFMG*) pc->data;
2045:   PetscScalar       *yy;
2046:   const PetscScalar *xx;
2047:   PetscInt          ilower[3],iupper[3];
2048:   Mat_HYPRESStruct  *mx     = (Mat_HYPRESStruct*)(pc->pmat->data);
2049:   PetscInt          ordering= mx->dofs_order;
2050:   PetscInt          nvars   = mx->nvars;
2051:   PetscInt          part    = 0;
2052:   PetscInt          size;
2053:   PetscInt          i;

2056:   PetscCitationsRegister(hypreCitation,&cite);
2057:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
2058:   iupper[0] += ilower[0] - 1;
2059:   iupper[1] += ilower[1] - 1;
2060:   iupper[2] += ilower[2] - 1;

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

2065:   /* copy x values over to hypre for variable ordering */
2066:   if (ordering) {
2067:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2068:     VecGetArrayRead(x,&xx);
2069:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)xx+(size*i)));
2070:     VecRestoreArrayRead(x,&xx);
2071:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2072:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
2073:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

2075:     /* copy solution values back to PETSc */
2076:     VecGetArray(y,&yy);
2077:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i)));
2078:     VecRestoreArray(y,&yy);
2079:   } else {      /* nodal ordering must be mapped to variable ordering for sys_pfmg */
2080:     PetscScalar *z;
2081:     PetscInt    j, k;

2083:     PetscMalloc1(nvars*size,&z);
2084:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2085:     VecGetArrayRead(x,&xx);

2087:     /* transform nodal to hypre's variable ordering for sys_pfmg */
2088:     for (i= 0; i< size; i++) {
2089:       k= i*nvars;
2090:       for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
2091:     }
2092:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
2093:     VecRestoreArrayRead(x,&xx);
2094:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2095:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

2097:     /* copy solution values back to PETSc */
2098:     VecGetArray(y,&yy);
2099:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
2100:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
2101:     for (i= 0; i< size; i++) {
2102:       k= i*nvars;
2103:       for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
2104:     }
2105:     VecRestoreArray(y,&yy);
2106:     PetscFree(z);
2107:   }
2108:   return(0);
2109: }

2113: 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)
2114: {
2115:   PC_SysPFMG     *jac = (PC_SysPFMG*)pc->data;
2117:   PetscInt       oits;

2120:   PetscCitationsRegister(hypreCitation,&cite);
2121:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its));
2122:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol));
2123:   PCApply_SysPFMG(pc,b,y);
2124:   PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,(HYPRE_Int *)&oits));
2125:   *outits = oits;
2126:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2127:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
2128:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol));
2129:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its));
2130:   return(0);
2131: }


2136: PetscErrorCode PCSetUp_SysPFMG(PC pc)
2137: {
2138:   PetscErrorCode   ierr;
2139:   PC_SysPFMG       *ex = (PC_SysPFMG*) pc->data;
2140:   Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
2141:   PetscBool        flg;

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

2147:   /* create the hypre sstruct solver object and set its information */
2148:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2149:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2150:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver));
2151:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2152:   return(0);
2153: }


2156: /*MC
2157:      PCSysPFMG - the hypre SysPFMG multigrid solver

2159:    Level: advanced

2161:    Options Database:
2162: + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner
2163: . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2164: . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2165: . -pc_syspfmg_tol <tol> tolerance of SysPFMG
2166: . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel

2168:    Notes:  This is for CELL-centered descretizations

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

2174: .seealso:  PCMG, MATHYPRESSTRUCT
2175: M*/

2179: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
2180: {
2182:   PC_SysPFMG     *ex;

2185:   PetscNew(&ex); \
2186:   pc->data = ex;

2188:   ex->its            = 1;
2189:   ex->tol            = 1.e-8;
2190:   ex->relax_type     = 1;
2191:   ex->num_pre_relax  = 1;
2192:   ex->num_post_relax = 1;

2194:   pc->ops->setfromoptions  = PCSetFromOptions_SysPFMG;
2195:   pc->ops->view            = PCView_SysPFMG;
2196:   pc->ops->destroy         = PCDestroy_SysPFMG;
2197:   pc->ops->apply           = PCApply_SysPFMG;
2198:   pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
2199:   pc->ops->setup           = PCSetUp_SysPFMG;

2201:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
2202:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2203:   return(0);
2204: }