Actual source code: hypre.c

petsc-master 2016-08-27
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>
 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  smoothtype;
 63:   PetscInt  smoothnumlevels;
 64:   PetscInt  eu_level;   /* Number of levels for ILU(k) in Euclid */
 65:   double    eu_droptolerance; /* Drop tolerance for ILU(k) in Euclid */
 66:   PetscInt  eu_bj;      /* Defines use of Block Jacobi ILU in Euclid */
 67:   PetscInt  relaxtype[3];
 68:   double    relaxweight;
 69:   double    outerrelaxweight;
 70:   PetscInt  relaxorder;
 71:   double    truncfactor;
 72:   PetscBool applyrichardson;
 73:   PetscInt  pmax;
 74:   PetscInt  interptype;
 75:   PetscInt  agg_nl;
 76:   PetscInt  agg_num_paths;
 77:   PetscInt  nodal_coarsen;
 78:   PetscBool nodal_relax;
 79:   PetscInt  nodal_relax_levels;

 81:   PetscInt  nodal_coarsening;
 82:   PetscInt  vec_interp_variant;
 83:   HYPRE_IJVector  *hmnull;
 84:   HYPRE_ParVector *phmnull;  /* near null space passed to hypre */
 85:   PetscInt        n_hmnull;
 86:   Vec             hmnull_constant;
 87:   PetscScalar     **hmnull_hypre_data_array;   /* this is the space in hmnull that was allocated by hypre, it is restored to hypre just before freeing the phmnull vectors */

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

104:   /* additional data */
105:   HYPRE_IJVector coords[3];
106:   HYPRE_IJVector constants[3];
107:   HYPRE_IJMatrix G;
108:   HYPRE_IJMatrix C;
109:   HYPRE_IJMatrix alpha_Poisson;
110:   HYPRE_IJMatrix beta_Poisson;
111:   PetscBool      ams_beta_is_zero;
112:   PetscBool      ams_beta_is_zero_part;
113:   PetscInt       ams_proj_freq;
114: } PC_HYPRE;

118: PetscErrorCode PCHYPREGetSolver(PC pc,HYPRE_Solver *hsolver)
119: {
120:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

123:   *hsolver = jac->hsolver;
124:   return(0);
125: }

127: /*
128:     Replaces the address where the HYPRE vector points to its data with the address of
129:   PETSc's data. Saves the old address so it can be reset when we are finished with it.
130:   Allows use to get the data into a HYPRE vector without the cost of memcopies
131: */
132: #define HYPREReplacePointer(b,newvalue,savedvalue) { \
133:     hypre_ParVector *par_vector   = (hypre_ParVector*)hypre_IJVectorObject(((hypre_IJVector*)b)); \
134:     hypre_Vector    *local_vector = hypre_ParVectorLocalVector(par_vector); \
135:     savedvalue         = local_vector->data; \
136:     local_vector->data = newvalue;          \
137: }

141: static PetscErrorCode PCSetUp_HYPRE(PC pc)
142: {
143:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
144:   PetscErrorCode     ierr;
145:   HYPRE_ParCSRMatrix hmat;
146:   HYPRE_ParVector    bv,xv;
147:   PetscInt           bs;

150:   if (!jac->hypre_type) {
151:     PCHYPRESetType(pc,"boomeramg");
152:   }

154:   if (pc->setupcalled) {
155:     /* always destroy the old matrix and create a new memory;
156:        hope this does not churn the memory too much. The problem
157:        is I do not know if it is possible to put the matrix back to
158:        its initial state so that we can directly copy the values
159:        the second time through. */
160:     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij));
161:     jac->ij = 0;
162:   }

164:   if (!jac->ij) { /* create the matrix the first time through */
165:     MatHYPRE_IJMatrixCreate(pc->pmat,&jac->ij);
166:   }
167:   if (!jac->b) { /* create the vectors the first time through */
168:     Vec x,b;
169:     MatCreateVecs(pc->pmat,&x,&b);
170:     VecHYPRE_IJVectorCreate(x,&jac->x);
171:     VecHYPRE_IJVectorCreate(b,&jac->b);
172:     VecDestroy(&x);
173:     VecDestroy(&b);
174:   }

176:   /* special case for BoomerAMG */
177:   if (jac->setup == HYPRE_BoomerAMGSetup) {
178:     MatNullSpace    mnull;
179:     PetscBool       has_const;
180:     PetscInt        nvec,i;
181:     const Vec       *vecs;
182:     PetscScalar     *petscvecarray;
183: 
184:     MatGetBlockSize(pc->pmat,&bs);
185:     if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
186:     MatGetNearNullSpace(pc->mat, &mnull);
187:     if (mnull) {
188:       MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);
189:       PetscMalloc1(nvec+1,&jac->hmnull);
190:       PetscMalloc1(nvec+1,&jac->hmnull_hypre_data_array);
191:       PetscMalloc1(nvec+1,&jac->phmnull);
192:       for (i=0; i<nvec; i++) {
193:         VecHYPRE_IJVectorCreate(vecs[i],&jac->hmnull[i]);
194:         VecGetArrayRead(vecs[i],(const PetscScalar **)&petscvecarray);
195:         HYPREReplacePointer(jac->hmnull[i],petscvecarray,jac->hmnull_hypre_data_array[i]);
196:         VecRestoreArrayRead(vecs[i],(const PetscScalar **)&petscvecarray);
197:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[i],(void**)&jac->phmnull[i]));
198:       }
199:       if (has_const) {
200:         MatCreateVecs(pc->pmat,&jac->hmnull_constant,NULL);
201:         VecSet(jac->hmnull_constant,1);
202:         VecNormalize(jac->hmnull_constant,NULL);
203:         VecHYPRE_IJVectorCreate(jac->hmnull_constant,&jac->hmnull[nvec]);
204:         VecGetArrayRead(jac->hmnull_constant,(const PetscScalar **)&petscvecarray);
205:         HYPREReplacePointer(jac->hmnull[nvec],petscvecarray,jac->hmnull_hypre_data_array[nvec]);
206:         VecRestoreArrayRead(jac->hmnull_constant,(const PetscScalar **)&petscvecarray);
207:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[nvec],(void**)&jac->phmnull[nvec]));
208:         nvec++;
209:       }
210:       PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVectors,(jac->hsolver,nvec,jac->phmnull));
211:       jac->n_hmnull = nvec;
212:     }
213:   }

215:   /* special case for AMS */
216:   if (jac->setup == HYPRE_AMSSetup) {
217:     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()");
218:     if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs discrete gradient operator via PCHYPRESetDiscreteGradient");
219:   }
220:   /* special case for ADS */
221:   if (jac->setup == HYPRE_ADSSetup) {
222:     if (!jac->coords[0]) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs coordinate vectors via PCSetCoordinates()");
223:     else if (!jac->coords[1] || !jac->coords[2]) 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");
224:     if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs discrete gradient operator via PCHYPRESetDiscreteGradient");
225:     if (!jac->C) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs discrete curl operator via PCHYPRESetDiscreteGradient");
226:   }
227:   MatHYPRE_IJMatrixCopy(pc->pmat,jac->ij);
228:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
229:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&bv));
230:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&xv));
231:   PetscStackCall("HYPRE_SetupXXX",(*jac->setup)(jac->hsolver,hmat,bv,xv););
232:   return(0);
233: }

237: static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
238: {
239:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
240:   PetscErrorCode     ierr;
241:   HYPRE_ParCSRMatrix hmat;
242:   PetscScalar        *xv;
243:   const PetscScalar  *bv,*sbv;
244:   HYPRE_ParVector    jbv,jxv;
245:   PetscScalar        *sxv;
246:   PetscInt           hierr;

249:   PetscCitationsRegister(hypreCitation,&cite);
250:   if (!jac->applyrichardson) {VecSet(x,0.0);}
251:   VecGetArrayRead(b,&bv);
252:   VecGetArray(x,&xv);
253:   HYPREReplacePointer(jac->b,(PetscScalar*)bv,sbv);
254:   HYPREReplacePointer(jac->x,xv,sxv);

256:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
257:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
258:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));
259:   PetscStackCall("Hypre solve",h(*jac->solve)(jac->hsolver,hmat,jbv,jxv);
260:                                if (hierr && hierr != HYPRE_ERROR_CONV) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
261:                                if (hierr) hypre__global_error = 0;);

263:   if (jac->setup == HYPRE_AMSSetup && jac->ams_beta_is_zero_part) {
264:     PetscStackCall("HYPRE_AMSProjectOutGradients",HYPRE_AMSProjectOutGradients(jac->hsolver,jxv););
265:   }
266:   HYPREReplacePointer(jac->b,(PetscScalar*)sbv,bv);
267:   HYPREReplacePointer(jac->x,sxv,xv);
268:   VecRestoreArray(x,&xv);
269:   VecRestoreArrayRead(b,&bv);
270:   return(0);
271: }

275: static PetscErrorCode PCReset_HYPRE(PC pc)
276: {
277:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;

281:   if (jac->ij) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij)); jac->ij = NULL;
282:   if (jac->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->b)); jac->b = NULL;
283:   if (jac->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->x)); jac->x = NULL;
284:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0])); jac->coords[0] = NULL;
285:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1])); jac->coords[1] = NULL;
286:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2])); jac->coords[2] = NULL;
287:   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0])); jac->constants[0] = NULL;
288:   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1])); jac->constants[1] = NULL;
289:   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2])); jac->constants[2] = NULL;
290:   if (jac->G) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->G)); jac->G = NULL;
291:   if (jac->C) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->C)); jac->C = NULL;
292:   if (jac->alpha_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->alpha_Poisson)); jac->alpha_Poisson = NULL;
293:   if (jac->beta_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->beta_Poisson)); jac->beta_Poisson = NULL;
294:   if (jac->n_hmnull && jac->hmnull) {
295:     PetscInt                 i;
296:     PETSC_UNUSED PetscScalar *petscvecarray;

298:     for (i=0; i<jac->n_hmnull; i++) {
299:       HYPREReplacePointer(jac->hmnull[i],jac->hmnull_hypre_data_array[i],petscvecarray);
300:       PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->hmnull[i]));
301:     }
302:     PetscFree(jac->hmnull);
303:     PetscFree(jac->hmnull_hypre_data_array);
304:     PetscFree(jac->phmnull);
305:     VecDestroy(&jac->hmnull_constant);
306:   }
307:   return(0);
308: }

312: static PetscErrorCode PCDestroy_HYPRE(PC pc)
313: {
314:   PC_HYPRE                 *jac = (PC_HYPRE*)pc->data;
315:   PetscErrorCode           ierr;

318:   PCReset_HYPRE(pc);
319:   if (jac->destroy) PetscStackCall("HYPRE_DestroyXXX",(*jac->destroy)(jac->hsolver););
320:   PetscFree(jac->hypre_type);
321:   if (jac->comm_hypre != MPI_COMM_NULL) { MPI_Comm_free(&(jac->comm_hypre));}
322:   PetscFree(pc->data);

324:   PetscObjectChangeTypeName((PetscObject)pc,0);
325:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);
326:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);
327:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetCoordinates_C",NULL);
328:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",NULL);
329:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",NULL);
330:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetConstantEdgeVectors_C",NULL);
331:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetAlphaPoissonMatrix_C",NULL);
332:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetBetaPoissonMatrix_C",NULL);
333:   return(0);
334: }

336: /* --------------------------------------------------------------------------------------------*/
339: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PetscOptionItems *PetscOptionsObject,PC pc)
340: {
341:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
343:   PetscBool      flag;

346:   PetscOptionsHead(PetscOptionsObject,"HYPRE Pilut Options");
347:   PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
348:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter));
349:   PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
350:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol));
351:   PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
352:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize));
353:   PetscOptionsTail();
354:   return(0);
355: }

359: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
360: {
361:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
363:   PetscBool      iascii;

366:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
367:   if (iascii) {
368:     PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut preconditioning\n");
369:     if (jac->maxiter != PETSC_DEFAULT) {
370:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: maximum number of iterations %d\n",jac->maxiter);
371:     } else {
372:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default maximum number of iterations \n");
373:     }
374:     if (jac->tol != PETSC_DEFAULT) {
375:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: drop tolerance %g\n",(double)jac->tol);
376:     } else {
377:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default drop tolerance \n");
378:     }
379:     if (jac->factorrowsize != PETSC_DEFAULT) {
380:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: factor row size %d\n",jac->factorrowsize);
381:     } else {
382:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default factor row size \n");
383:     }
384:   }
385:   return(0);
386: }

388: /* --------------------------------------------------------------------------------------------*/

392: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
393: {
394:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
395:   PetscErrorCode     ierr;
396:   HYPRE_ParCSRMatrix hmat;
397:   PetscScalar        *xv;
398:   const PetscScalar  *bv;
399:   HYPRE_ParVector    jbv,jxv;
400:   PetscScalar        *sbv,*sxv;
401:   PetscInt           hierr;

404:   PetscCitationsRegister(hypreCitation,&cite);
405:   VecSet(x,0.0);
406:   VecGetArrayRead(b,&bv);
407:   VecGetArray(x,&xv);
408:   HYPREReplacePointer(jac->b,(PetscScalar*)bv,sbv);
409:   HYPREReplacePointer(jac->x,xv,sxv);

411:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
412:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
413:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));

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

420:   HYPREReplacePointer(jac->b,sbv,bv);
421:   HYPREReplacePointer(jac->x,sxv,xv);
422:   VecRestoreArray(x,&xv);
423:   VecRestoreArrayRead(b,&bv);
424:   return(0);
425: }

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

430: static const char *HYPREBoomerAMGCycleType[]   = {"","V","W"};
431: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
432: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
433: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
434: static const char *HYPREBoomerAMGSmoothType[]   = {"Schwarz-smoothers","Pilut","ParaSails","Euclid"};
435: static const char *HYPREBoomerAMGRelaxType[]   = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi",
436:                                                   "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi",
437:                                                   "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination",
438:                                                   "" /* 10 */, "" /* 11 */, "" /* 12 */, "" /* 13 */, "" /* 14 */,
439:                                                   "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"};
440: static const char *HYPREBoomerAMGInterpType[]  = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
441:                                                   "ext+i-cc", "standard", "standard-wts", "", "", "FF", "FF1"};
444: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PetscOptionItems *PetscOptionsObject,PC pc)
445: {
446:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
448:   PetscInt       n,indx,level;
449:   PetscBool      flg, tmp_truth;
450:   double         tmpdbl, twodbl[2];

453:   PetscOptionsHead(PetscOptionsObject,"HYPRE BoomerAMG Options");
454:   PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);
455:   if (flg) {
456:     jac->cycletype = indx+1;
457:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
458:   }
459:   PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
460:   if (flg) {
461:     if (jac->maxlevels < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
462:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
463:   }
464:   PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
465:   if (flg) {
466:     if (jac->maxiter < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
467:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
468:   }
469:   PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);
470:   if (flg) {
471:     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);
472:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
473:   }

475:   PetscOptionsScalar("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);
476:   if (flg) {
477:     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);
478:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
479:   }

481:   PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);
482:   if (flg) {
483:     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);
484:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
485:   }

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

491:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
492:   }


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

499:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
500:   }


503:   PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
504:   if (flg) {
505:     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);
506:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
507:   }
508:   PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
509:   if (flg) {
510:     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);
511:     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);
512:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
513:   }

515:   /* Grid sweeps */
516:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);
517:   if (flg) {
518:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx));
519:     /* modify the jac structure so we can view the updated options with PC_View */
520:     jac->gridsweeps[0] = indx;
521:     jac->gridsweeps[1] = indx;
522:     /*defaults coarse to 1 */
523:     jac->gridsweeps[2] = 1;
524:   }

526:   PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen","Use a nodal based coarsening 1-6","HYPRE_BoomerAMGSetNodal",jac->nodal_coarsening, &jac->nodal_coarsening,&flg);
527:   if (flg) {
528:     PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,jac->nodal_coarsening));
529:   }

531:   PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_variant","Variant of algorithm 1-3","HYPRE_BoomerAMGSetInterpVecVariant",jac->vec_interp_variant, &jac->vec_interp_variant,&flg);
532:   if (flg) {
533:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecVariant,(jac->hsolver,jac->vec_interp_variant));
534:   }

536:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);
537:   if (flg) {
538:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1));
539:     jac->gridsweeps[0] = indx;
540:   }
541:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);
542:   if (flg) {
543:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2));
544:     jac->gridsweeps[1] = indx;
545:   }
546:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);
547:   if (flg) {
548:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3));
549:     jac->gridsweeps[2] = indx;
550:   }

552:   /* Smooth type */
553:   PetscOptionsEList("-pc_hypre_boomeramg_smooth_type","Enable more complex smoothers","None",HYPREBoomerAMGSmoothType,ALEN(HYPREBoomerAMGSmoothType),HYPREBoomerAMGSmoothType[0],&indx,&flg);
554:   if (flg) {
555:     jac->smoothtype = indx;
556:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,indx+6));
557:     jac->smoothnumlevels = 25;
558:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,25));
559:   }

561:   /* Number of smoothing levels */
562:   PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_levels","Number of levels on which more complex smoothers are used","None",25,&indx,&flg);
563:   if (flg && (jac->smoothtype != -1)) {
564:     jac->smoothnumlevels = indx;
565:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,indx));
566:   }

568:   /* Number of levels for ILU(k) for Euclid */
569:   PetscOptionsInt("-pc_hypre_boomeramg_eu_level","Number of levels for ILU(k) in Euclid smoother","None",0,&indx,&flg);
570:   if (flg && (jac->smoothtype == 3)) {
571:     jac->eu_level = indx;
572:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,indx));
573:   }

575:   /* Filter for ILU(k) for Euclid */
576:   double droptolerance;
577:   PetscOptionsScalar("-pc_hypre_boomeramg_eu_droptolerance","Drop tolerance for ILU(k) in Euclid smoother","None",0,&droptolerance,&flg);
578:   if (flg && (jac->smoothtype == 3)) {
579:     jac->eu_droptolerance = droptolerance;
580:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,droptolerance));
581:   }

583:   /* Use Block Jacobi ILUT for Euclid */
584:   PetscOptionsBool("-pc_hypre_boomeramg_eu_bj", "Use Block Jacobi for ILU in Euclid smoother?", "None", PETSC_FALSE, &tmp_truth, &flg);
585:   if (flg && (jac->smoothtype == 3)) {
586:     jac->eu_bj = tmp_truth;
587:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuBJ,(jac->hsolver,jac->eu_bj));
588:   }

590:   /* Relax type */
591:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
592:   if (flg) {
593:     jac->relaxtype[0] = jac->relaxtype[1]  = indx;
594:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx));
595:     /* by default, coarse type set to 9 */
596:     jac->relaxtype[2] = 9;

598:   }
599:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
600:   if (flg) {
601:     jac->relaxtype[0] = indx;
602:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1));
603:   }
604:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
605:   if (flg) {
606:     jac->relaxtype[1] = indx;
607:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2));
608:   }
609:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);
610:   if (flg) {
611:     jac->relaxtype[2] = indx;
612:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3));
613:   }

615:   /* Relaxation Weight */
616:   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);
617:   if (flg) {
618:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl));
619:     jac->relaxweight = tmpdbl;
620:   }

622:   n         = 2;
623:   twodbl[0] = twodbl[1] = 1.0;
624:   PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
625:   if (flg) {
626:     if (n == 2) {
627:       indx =  (int)PetscAbsReal(twodbl[1]);
628:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx));
629:     } 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);
630:   }

632:   /* Outer relaxation Weight */
633:   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);
634:   if (flg) {
635:     PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl));
636:     jac->outerrelaxweight = tmpdbl;
637:   }

639:   n         = 2;
640:   twodbl[0] = twodbl[1] = 1.0;
641:   PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
642:   if (flg) {
643:     if (n == 2) {
644:       indx =  (int)PetscAbsReal(twodbl[1]);
645:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx));
646:     } 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);
647:   }

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

652:   if (flg && tmp_truth) {
653:     jac->relaxorder = 0;
654:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
655:   }
656:   PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);
657:   if (flg) {
658:     jac->measuretype = indx;
659:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
660:   }
661:   /* update list length 3/07 */
662:   PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);
663:   if (flg) {
664:     jac->coarsentype = indx;
665:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
666:   }

668:   /* new 3/07 */
669:   PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);
670:   if (flg) {
671:     jac->interptype = indx;
672:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
673:   }

675:   PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
676:   if (flg) {
677:     level = 3;
678:     PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);

680:     jac->printstatistics = PETSC_TRUE;
681:     PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level));
682:   }

684:   PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);
685:   if (flg) {
686:     level = 3;
687:     PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);

689:     jac->printstatistics = PETSC_TRUE;
690:     PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level));
691:   }

693:   PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);
694:   if (flg && tmp_truth) {
695:     PetscInt tmp_int;
696:     PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);
697:     if (flg) jac->nodal_relax_levels = tmp_int;
698:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6));
699:     PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1));
700:     PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0));
701:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels));
702:   }

704:   PetscOptionsTail();
705:   return(0);
706: }

710: 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)
711: {
712:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
714:   PetscInt       oits;

717:   PetscCitationsRegister(hypreCitation,&cite);
718:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter));
719:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol));
720:   jac->applyrichardson = PETSC_TRUE;
721:   PCApply_HYPRE(pc,b,y);
722:   jac->applyrichardson = PETSC_FALSE;
723:   PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
724:   *outits = oits;
725:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
726:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
727:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
728:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
729:   return(0);
730: }


735: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
736: {
737:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
739:   PetscBool      iascii;

742:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
743:   if (iascii) {
744:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG preconditioning\n");
745:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);
746:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of levels %d\n",jac->maxlevels);
747:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of iterations PER hypre call %d\n",jac->maxiter);
748:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Convergence tolerance PER hypre call %g\n",(double)jac->tol);
749:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Threshold for strong coupling %g\n",(double)jac->strongthreshold);
750:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation truncation factor %g\n",(double)jac->truncfactor);
751:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation: max elements per row %d\n",jac->pmax);
752:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Number of levels of aggressive coarsening %d\n",jac->agg_nl);
753:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Number of paths for aggressive coarsening %d\n",jac->agg_num_paths);

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

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

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

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

768:     if (jac->relaxorder) {
769:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Using CF-relaxation\n");
770:     } else {
771:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Not using CF-relaxation\n");
772:     }
773:     if (jac->smoothtype!=-1) {
774:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Smooth type          %s\n",HYPREBoomerAMGSmoothType[jac->smoothtype]);
775:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Smooth num levels    %d\n",jac->smoothnumlevels);
776:     } else {
777:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Not using more complex smoothers.\n");
778:     }
779:     if (jac->smoothtype==3) {
780:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Euclid ILU(k) levels %d\n",jac->eu_level);
781:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Euclid ILU(k) drop tolerance %g\n",jac->eu_droptolerance);
782:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Euclid ILU use Block-Jacobi? %d\n",jac->eu_bj);
783:     }
784:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Measure type        %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
785:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Coarsen type        %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
786:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation type  %s\n",HYPREBoomerAMGInterpType[jac->interptype]);
787:     if (jac->nodal_coarsening) {
788:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Using nodal coarsening (with HYPRE_BOOMERAMGSetNodal() %D\n",jac->nodal_coarsening);
789:     }
790:     if (jac->vec_interp_variant) {
791:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: HYPRE_BoomerAMGSetInterpVecVariant() %D\n",jac->vec_interp_variant);
792:     }
793:     if (jac->nodal_relax) {
794:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Using nodal relaxation via Schwarz smoothing on levels %d\n",jac->nodal_relax_levels);
795:     }
796:   }
797:   return(0);
798: }

800: /* --------------------------------------------------------------------------------------------*/
803: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PetscOptionItems *PetscOptionsObject,PC pc)
804: {
805:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
807:   PetscInt       indx;
808:   PetscBool      flag;
809:   const char     *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};

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

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

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

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

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

829:   PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);
830:   if (flag) {
831:     jac->symt = indx;
832:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
833:   }

835:   PetscOptionsTail();
836:   return(0);
837: }

841: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
842: {
843:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
845:   PetscBool      iascii;
846:   const char     *symt = 0;;

849:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
850:   if (iascii) {
851:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails preconditioning\n");
852:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: nlevels %d\n",jac->nlevels);
853:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: threshold %g\n",(double)jac->threshhold);
854:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: filter %g\n",(double)jac->filter);
855:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: load balance %g\n",(double)jac->loadbal);
856:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: reuse nonzero structure %s\n",PetscBools[jac->ruse]);
857:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: print info to screen %s\n",PetscBools[jac->logging]);
858:     if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
859:     else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
860:     else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
861:     else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
862:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: %s\n",symt);
863:   }
864:   return(0);
865: }
866: /* --------------------------------------------------------------------------------------------*/
869: static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PetscOptionItems *PetscOptionsObject,PC pc)
870: {
871:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
873:   PetscInt       n;
874:   PetscBool      flag,flag2,flag3,flag4;

877:   PetscOptionsHead(PetscOptionsObject,"HYPRE AMS Options");
878:   PetscOptionsInt("-pc_hypre_ams_print_level","Debugging output level for AMS","None",jac->as_print,&jac->as_print,&flag);
879:   if (flag) PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
880:   PetscOptionsInt("-pc_hypre_ams_max_iter","Maximum number of AMS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
881:   if (flag) PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
882:   PetscOptionsInt("-pc_hypre_ams_cycle_type","Cycle type for AMS multigrid","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag);
883:   if (flag) PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
884:   PetscOptionsReal("-pc_hypre_ams_tol","Error tolerance for AMS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
885:   if (flag) PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
886:   PetscOptionsInt("-pc_hypre_ams_relax_type","Relaxation type for AMS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
887:   PetscOptionsInt("-pc_hypre_ams_relax_times","Number of relaxation steps for AMS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
888:   PetscOptionsReal("-pc_hypre_ams_relax_weight","Relaxation weight for AMS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
889:   PetscOptionsReal("-pc_hypre_ams_omega","SSOR coefficient for AMS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
890:   if (flag || flag2 || flag3 || flag4) {
891:     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
892:                                                                       jac->as_relax_times,
893:                                                                       jac->as_relax_weight,
894:                                                                       jac->as_omega));
895:   }
896:   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);
897:   n = 5;
898:   PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options","AMG options for vector Poisson","None",jac->as_amg_alpha_opts,&n,&flag2);
899:   if (flag || flag2) {
900:     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
901:                                                                      jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
902:                                                                      jac->as_amg_alpha_opts[2],       /* AMG relax_type */
903:                                                                      jac->as_amg_alpha_theta,
904:                                                                      jac->as_amg_alpha_opts[3],       /* AMG interp_type */
905:                                                                      jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
906:   }
907:   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);
908:   n = 5;
909:   PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options","AMG options for scalar Poisson solver","None",jac->as_amg_beta_opts,&n,&flag2);
910:   if (flag || flag2) {
911:     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
912:                                                                     jac->as_amg_beta_opts[1],       /* AMG agg_levels */
913:                                                                     jac->as_amg_beta_opts[2],       /* AMG relax_type */
914:                                                                     jac->as_amg_beta_theta,
915:                                                                     jac->as_amg_beta_opts[3],       /* AMG interp_type */
916:                                                                     jac->as_amg_beta_opts[4]));     /* AMG Pmax */
917:   }
918:   PetscOptionsInt("-pc_hypre_ams_projection_frequency","Frequency at which a projection onto the compatible subspace for problems with zero conductivity regions is performed","None",jac->ams_proj_freq,&jac->ams_proj_freq,&flag);
919:   if (flag) { /* override HYPRE's default only if the options is used */
920:     PetscStackCallStandard(HYPRE_AMSSetProjectionFrequency,(jac->hsolver,jac->ams_proj_freq));
921:   }
922:   PetscOptionsTail();
923:   return(0);
924: }

928: static PetscErrorCode PCView_HYPRE_AMS(PC pc,PetscViewer viewer)
929: {
930:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
932:   PetscBool      iascii;

935:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
936:   if (iascii) {
937:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS preconditioning\n");
938:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace iterations per application %d\n",jac->as_max_iter);
939:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace cycle type %d\n",jac->ams_cycle_type);
940:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace iteration tolerance %g\n",jac->as_tol);
941:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother type %d\n",jac->as_relax_type);
942:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: number of smoothing steps %d\n",jac->as_relax_times);
943:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother weight %g\n",jac->as_relax_weight);
944:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother omega %g\n",jac->as_omega);
945:     if (jac->alpha_Poisson) {
946:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: vector Poisson solver (passed in by user)\n");
947:     } else {
948:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: vector Poisson solver (computed) \n");
949:     }
950:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);
951:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
952:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);
953:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);
954:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
955:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);
956:     if (!jac->ams_beta_is_zero) {
957:       if (jac->beta_Poisson) {
958:         PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: scalar Poisson solver (passed in by user)\n");
959:       } else {
960:         PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: scalar Poisson solver (computed) \n");
961:       }
962:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);
963:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
964:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);
965:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);
966:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
967:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);
968:       if (jac->ams_beta_is_zero_part) {
969:         PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     compatible subspace projection frequency %d (-1 HYPRE uses default)\n",jac->ams_proj_freq);
970:       }
971:     } else {
972:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: scalar Poisson solver not used (zero-conductivity everywhere) \n");
973:     }
974:   }
975:   return(0);
976: }

980: static PetscErrorCode PCSetFromOptions_HYPRE_ADS(PetscOptionItems *PetscOptionsObject,PC pc)
981: {
982:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
984:   PetscInt       n;
985:   PetscBool      flag,flag2,flag3,flag4;

988:   PetscOptionsHead(PetscOptionsObject,"HYPRE ADS Options");
989:   PetscOptionsInt("-pc_hypre_ads_print_level","Debugging output level for ADS","None",jac->as_print,&jac->as_print,&flag);
990:   if (flag) PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
991:   PetscOptionsInt("-pc_hypre_ads_max_iter","Maximum number of ADS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
992:   if (flag) PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
993:   PetscOptionsInt("-pc_hypre_ads_cycle_type","Cycle type for ADS multigrid","None",jac->ads_cycle_type,&jac->ads_cycle_type,&flag);
994:   if (flag) PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ads_cycle_type));
995:   PetscOptionsReal("-pc_hypre_ads_tol","Error tolerance for ADS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
996:   if (flag) PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
997:   PetscOptionsInt("-pc_hypre_ads_relax_type","Relaxation type for ADS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
998:   PetscOptionsInt("-pc_hypre_ads_relax_times","Number of relaxation steps for ADS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
999:   PetscOptionsReal("-pc_hypre_ads_relax_weight","Relaxation weight for ADS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
1000:   PetscOptionsReal("-pc_hypre_ads_omega","SSOR coefficient for ADS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
1001:   if (flag || flag2 || flag3 || flag4) {
1002:     PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1003:                                                                       jac->as_relax_times,
1004:                                                                       jac->as_relax_weight,
1005:                                                                       jac->as_omega));
1006:   }
1007:   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);
1008:   n = 5;
1009:   PetscOptionsIntArray("-pc_hypre_ads_ams_options","AMG options for AMS solver inside ADS","None",jac->as_amg_alpha_opts,&n,&flag2);
1010:   PetscOptionsInt("-pc_hypre_ads_ams_cycle_type","Cycle type for AMS solver inside ADS","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag3);
1011:   if (flag || flag2 || flag3) {
1012:     PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type,             /* AMS cycle type */
1013:                                                                 jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1014:                                                                 jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1015:                                                                 jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1016:                                                                 jac->as_amg_alpha_theta,
1017:                                                                 jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1018:                                                                 jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1019:   }
1020:   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);
1021:   n = 5;
1022:   PetscOptionsIntArray("-pc_hypre_ads_amg_options","AMG options for vector AMG solver inside ADS","None",jac->as_amg_beta_opts,&n,&flag2);
1023:   if (flag || flag2) {
1024:     PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1025:                                                                 jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1026:                                                                 jac->as_amg_beta_opts[2],       /* AMG relax_type */
1027:                                                                 jac->as_amg_beta_theta,
1028:                                                                 jac->as_amg_beta_opts[3],       /* AMG interp_type */
1029:                                                                 jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1030:   }
1031:   PetscOptionsTail();
1032:   return(0);
1033: }

1037: static PetscErrorCode PCView_HYPRE_ADS(PC pc,PetscViewer viewer)
1038: {
1039:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1041:   PetscBool      iascii;

1044:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1045:   if (iascii) {
1046:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS preconditioning\n");
1047:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: subspace iterations per application %d\n",jac->as_max_iter);
1048:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: subspace cycle type %d\n",jac->ads_cycle_type);
1049:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: subspace iteration tolerance %g\n",jac->as_tol);
1050:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: smoother type %d\n",jac->as_relax_type);
1051:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: number of smoothing steps %d\n",jac->as_relax_times);
1052:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: smoother weight %g\n",jac->as_relax_weight);
1053:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: smoother omega %g\n",jac->as_omega);
1054:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: AMS solver\n");
1055:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     subspace cycle type %d\n",jac->ams_cycle_type);
1056:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);
1057:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
1058:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);
1059:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);
1060:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
1061:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);
1062:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: vector Poisson solver\n");
1063:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);
1064:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
1065:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);
1066:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);
1067:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
1068:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);
1069:   }
1070:   return(0);
1071: }

1075: static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE(PC pc, Mat G)
1076: {
1077:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1078:   HYPRE_ParCSRMatrix parcsr_G;
1079:   PetscErrorCode     ierr;

1082:   /* throw away any discrete gradient if already set */
1083:   if (jac->G) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->G));
1084:   MatHYPRE_IJMatrixCreate(G,&jac->G);
1085:   MatHYPRE_IJMatrixCopy(G,jac->G);
1086:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->G,(void**)(&parcsr_G)));
1087:   PetscStackCall("Hypre set gradient",(*jac->setdgrad)(jac->hsolver,parcsr_G););
1088:   return(0);
1089: }

1093: /*@
1094:  PCHYPRESetDiscreteGradient - Set discrete gradient matrix

1096:    Collective on PC

1098:    Input Parameters:
1099: +  pc - the preconditioning context
1100: -  G - the discrete gradient

1102:    Level: intermediate

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

1107: .seealso:
1108: @*/
1109: PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G)
1110: {

1117:   PetscTryMethod(pc,"PCHYPRESetDiscreteGradient_C",(PC,Mat),(pc,G));
1118:   return(0);
1119: }

1123: static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C)
1124: {
1125:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1126:   HYPRE_ParCSRMatrix parcsr_C;
1127:   PetscErrorCode     ierr;

1130:   /* throw away any discrete curl if already set */
1131:   if (jac->C) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->C));
1132:   MatHYPRE_IJMatrixCreate(C,&jac->C);
1133:   MatHYPRE_IJMatrixCopy(C,jac->C);
1134:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->C,(void**)(&parcsr_C)));
1135:   PetscStackCall("Hypre set curl",(*jac->setdcurl)(jac->hsolver,parcsr_C););
1136:   return(0);
1137: }

1141: /*@
1142:  PCHYPRESetDiscreteCurl - Set discrete curl matrix

1144:    Collective on PC

1146:    Input Parameters:
1147: +  pc - the preconditioning context
1148: -  C - the discrete curl

1150:    Level: intermediate

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

1155: .seealso:
1156: @*/
1157: PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C)
1158: {

1165:   PetscTryMethod(pc,"PCHYPRESetDiscreteCurl_C",(PC,Mat),(pc,C));
1166:   return(0);
1167: }

1171: static PetscErrorCode PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS(PC pc, Mat A)
1172: {
1173:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1174:   HYPRE_ParCSRMatrix parcsr_alpha_Poisson;
1175:   PetscErrorCode     ierr;

1178:   /* throw away any matrix if already set */
1179:   if (jac->alpha_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->alpha_Poisson));
1180:   MatHYPRE_IJMatrixCreate(A,&jac->alpha_Poisson);
1181:   MatHYPRE_IJMatrixCopy(A,jac->alpha_Poisson);
1182:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->alpha_Poisson,(void**)(&parcsr_alpha_Poisson)));
1183:   PetscStackCallStandard(HYPRE_AMSSetAlphaPoissonMatrix,(jac->hsolver,parcsr_alpha_Poisson));
1184:   return(0);
1185: }

1189: /*@
1190:  PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix

1192:    Collective on PC

1194:    Input Parameters:
1195: +  pc - the preconditioning context
1196: -  A - the matrix

1198:    Level: intermediate

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

1202: .seealso:
1203: @*/
1204: PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A)
1205: {

1212:   PetscTryMethod(pc,"PCHYPRESetAlphaPoissonMatrix_C",(PC,Mat),(pc,A));
1213:   return(0);
1214: }

1218: static PetscErrorCode PCHYPRESetBetaPoissonMatrix_HYPRE_AMS(PC pc, Mat A)
1219: {
1220:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1221:   HYPRE_ParCSRMatrix parcsr_beta_Poisson;
1222:   PetscErrorCode     ierr;

1225:   if (!A) {
1226:     PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,NULL));
1227:     jac->ams_beta_is_zero = PETSC_TRUE;
1228:     return(0);
1229:   }
1230:   /* throw away any matrix if already set */
1231:   if (jac->beta_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->beta_Poisson));
1232:   MatHYPRE_IJMatrixCreate(A,&jac->beta_Poisson);
1233:   MatHYPRE_IJMatrixCopy(A,jac->beta_Poisson);
1234:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->beta_Poisson,(void**)(&parcsr_beta_Poisson)));
1235:   PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,parcsr_beta_Poisson));
1236:   return(0);
1237: }

1241: /*@
1242:  PCHYPRESetBetaPoissonMatrix - Set Poisson matrix

1244:    Collective on PC

1246:    Input Parameters:
1247: +  pc - the preconditioning context
1248: -  A - the matrix

1250:    Level: intermediate

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

1255: .seealso:
1256: @*/
1257: PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A)
1258: {

1263:   if (A) {
1266:   }
1267:   PetscTryMethod(pc,"PCHYPRESetBetaPoissonMatrix_C",(PC,Mat),(pc,A));
1268:   return(0);
1269: }

1273: static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE_AMS(PC pc,Vec ozz, Vec zoz, Vec zzo)
1274: {
1275:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1276:   HYPRE_ParVector    par_ozz,par_zoz,par_zzo;
1277:   PetscInt           dim;
1278:   PetscErrorCode     ierr;

1281:   /* throw away any vector if already set */
1282:   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0]));
1283:   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1]));
1284:   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2]));
1285:   jac->constants[0] = NULL;
1286:   jac->constants[1] = NULL;
1287:   jac->constants[2] = NULL;
1288:   VecHYPRE_IJVectorCreate(ozz,&jac->constants[0]);
1289:   VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);
1290:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[0],(void**)(&par_ozz)));
1291:   VecHYPRE_IJVectorCreate(zoz,&jac->constants[1]);
1292:   VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);
1293:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[1],(void**)(&par_zoz)));
1294:   dim = 2;
1295:   if (zzo) {
1296:     VecHYPRE_IJVectorCreate(zzo,&jac->constants[2]);
1297:     VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);
1298:     PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[2],(void**)(&par_zzo)));
1299:     dim++;
1300:   }
1301:   PetscStackCallStandard(HYPRE_AMSSetEdgeConstantVectors,(jac->hsolver,par_ozz,par_zoz,par_zzo));
1302:   PetscStackCallStandard(HYPRE_AMSSetDimension,(jac->hsolver,dim));
1303:   return(0);
1304: }

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

1311:    Collective on PC

1313:    Input Parameters:
1314: +  pc - the preconditioning context
1315: -  ozz - vector representing (1,0,0) (or (1,0) in 2D)
1316: -  zoz - vector representing (0,1,0) (or (0,1) in 2D)
1317: -  zzo - vector representing (0,0,1) (use NULL in 2D)

1319:    Level: intermediate

1321:    Notes:

1323: .seealso:
1324: @*/
1325: PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo)
1326: {

1337:   PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));
1338:   return(0);
1339: }

1343: static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1344: {
1345:   PC_HYPRE        *jac = (PC_HYPRE*)pc->data;
1346:   Vec             tv;
1347:   HYPRE_ParVector par_coords[3];
1348:   PetscInt        i;
1349:   PetscErrorCode  ierr;

1352:   /* throw away any coordinate vector if already set */
1353:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0]));
1354:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1]));
1355:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2]));
1356:   /* set problem's dimension */
1357:   if (jac->setdim) {
1358:     PetscStackCall("Hypre set dim",(*jac->setdim)(jac->hsolver,dim););
1359:   }
1360:   /* compute IJ vector for coordinates */
1361:   VecCreate(PetscObjectComm((PetscObject)pc),&tv);
1362:   VecSetType(tv,VECSTANDARD);
1363:   VecSetSizes(tv,nloc,PETSC_DECIDE);
1364:   for (i=0;i<dim;i++) {
1365:     PetscScalar *array;
1366:     PetscInt    j;

1368:     VecHYPRE_IJVectorCreate(tv,&jac->coords[i]);
1369:     VecGetArray(tv,&array);
1370:     for (j=0;j<nloc;j++) {
1371:       array[j] = coords[j*dim+i];
1372:     }
1373:     PetscStackCallStandard(HYPRE_IJVectorSetValues,(jac->coords[i],nloc,NULL,array));
1374:     PetscStackCallStandard(HYPRE_IJVectorAssemble,(jac->coords[i]));
1375:     VecRestoreArray(tv,&array);
1376:   }
1377:   VecDestroy(&tv);
1378:   /* pass parCSR vectors to AMS solver */
1379:   par_coords[0] = NULL;
1380:   par_coords[1] = NULL;
1381:   par_coords[2] = NULL;
1382:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0],(void**)(&par_coords[0])));
1383:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1],(void**)(&par_coords[1])));
1384:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2],(void**)(&par_coords[2])));
1385:   PetscStackCall("Hypre set coords",(*jac->setcoord)(jac->hsolver,par_coords[0],par_coords[1],par_coords[2]););
1386:   return(0);
1387: }

1389: /* ---------------------------------------------------------------------------------*/

1393: static PetscErrorCode  PCHYPREGetType_HYPRE(PC pc,const char *name[])
1394: {
1395:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

1398:   *name = jac->hypre_type;
1399:   return(0);
1400: }

1404: static PetscErrorCode  PCHYPRESetType_HYPRE(PC pc,const char name[])
1405: {
1406:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1408:   PetscBool      flag;

1411:   if (jac->hypre_type) {
1412:     PetscStrcmp(jac->hypre_type,name,&flag);
1413:     if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set");
1414:     return(0);
1415:   } else {
1416:     PetscStrallocpy(name, &jac->hypre_type);
1417:   }

1419:   jac->maxiter         = PETSC_DEFAULT;
1420:   jac->tol             = PETSC_DEFAULT;
1421:   jac->printstatistics = PetscLogPrintInfo;

1423:   PetscStrcmp("pilut",jac->hypre_type,&flag);
1424:   if (flag) {
1425:     PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver));
1426:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
1427:     pc->ops->view           = PCView_HYPRE_Pilut;
1428:     jac->destroy            = HYPRE_ParCSRPilutDestroy;
1429:     jac->setup              = HYPRE_ParCSRPilutSetup;
1430:     jac->solve              = HYPRE_ParCSRPilutSolve;
1431:     jac->factorrowsize      = PETSC_DEFAULT;
1432:     return(0);
1433:   }
1434:   PetscStrcmp("parasails",jac->hypre_type,&flag);
1435:   if (flag) {
1436:     PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver));
1437:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
1438:     pc->ops->view           = PCView_HYPRE_ParaSails;
1439:     jac->destroy            = HYPRE_ParaSailsDestroy;
1440:     jac->setup              = HYPRE_ParaSailsSetup;
1441:     jac->solve              = HYPRE_ParaSailsSolve;
1442:     /* initialize */
1443:     jac->nlevels    = 1;
1444:     jac->threshhold = .1;
1445:     jac->filter     = .1;
1446:     jac->loadbal    = 0;
1447:     if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE;
1448:     else jac->logging = (int) PETSC_FALSE;

1450:     jac->ruse = (int) PETSC_FALSE;
1451:     jac->symt = 0;
1452:     PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
1453:     PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
1454:     PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
1455:     PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
1456:     PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
1457:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
1458:     return(0);
1459:   }
1460:   PetscStrcmp("boomeramg",jac->hypre_type,&flag);
1461:   if (flag) {
1462:     HYPRE_BoomerAMGCreate(&jac->hsolver);
1463:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_BoomerAMG;
1464:     pc->ops->view            = PCView_HYPRE_BoomerAMG;
1465:     pc->ops->applytranspose  = PCApplyTranspose_HYPRE_BoomerAMG;
1466:     pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
1467:     jac->destroy             = HYPRE_BoomerAMGDestroy;
1468:     jac->setup               = HYPRE_BoomerAMGSetup;
1469:     jac->solve               = HYPRE_BoomerAMGSolve;
1470:     jac->applyrichardson     = PETSC_FALSE;
1471:     /* these defaults match the hypre defaults */
1472:     jac->cycletype        = 1;
1473:     jac->maxlevels        = 25;
1474:     jac->maxiter          = 1;
1475:     jac->tol              = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
1476:     jac->truncfactor      = 0.0;
1477:     jac->strongthreshold  = .25;
1478:     jac->maxrowsum        = .9;
1479:     jac->coarsentype      = 6;
1480:     jac->measuretype      = 0;
1481:     jac->gridsweeps[0]    = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
1482:     jac->smoothtype       = -1; /* Not set by default */
1483:     jac->smoothnumlevels  = 25;
1484:     jac->eu_level         = 0;
1485:     jac->eu_droptolerance = 0;
1486:     jac->eu_bj            = 0;
1487:     jac->relaxtype[0]     = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */
1488:     jac->relaxtype[2]     = 9; /*G.E. */
1489:     jac->relaxweight      = 1.0;
1490:     jac->outerrelaxweight = 1.0;
1491:     jac->relaxorder       = 1;
1492:     jac->interptype       = 0;
1493:     jac->agg_nl           = 0;
1494:     jac->pmax             = 0;
1495:     jac->truncfactor      = 0.0;
1496:     jac->agg_num_paths    = 1;

1498:     jac->nodal_coarsen      = 0;
1499:     jac->nodal_relax        = PETSC_FALSE;
1500:     jac->nodal_relax_levels = 1;
1501:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
1502:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
1503:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
1504:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
1505:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
1506:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
1507:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
1508:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
1509:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
1510:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
1511:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
1512:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
1513:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
1514:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
1515:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0]));  /*defaults coarse to 9*/
1516:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */
1517:     return(0);
1518:   }
1519:   PetscStrcmp("ams",jac->hypre_type,&flag);
1520:   if (flag) {
1521:     HYPRE_AMSCreate(&jac->hsolver);
1522:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_AMS;
1523:     pc->ops->view            = PCView_HYPRE_AMS;
1524:     jac->destroy             = HYPRE_AMSDestroy;
1525:     jac->setup               = HYPRE_AMSSetup;
1526:     jac->solve               = HYPRE_AMSSolve;
1527:     jac->setdgrad            = HYPRE_AMSSetDiscreteGradient;
1528:     jac->setcoord            = HYPRE_AMSSetCoordinateVectors;
1529:     jac->setdim              = HYPRE_AMSSetDimension;
1530:     jac->coords[0]           = NULL;
1531:     jac->coords[1]           = NULL;
1532:     jac->coords[2]           = NULL;
1533:     jac->G                   = NULL;
1534:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */
1535:     jac->as_print           = 0;
1536:     jac->as_max_iter        = 1; /* used as a preconditioner */
1537:     jac->as_tol             = 0.; /* used as a preconditioner */
1538:     jac->ams_cycle_type     = 13;
1539:     /* Smoothing options */
1540:     jac->as_relax_type      = 2;
1541:     jac->as_relax_times     = 1;
1542:     jac->as_relax_weight    = 1.0;
1543:     jac->as_omega           = 1.0;
1544:     /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1545:     jac->as_amg_alpha_opts[0] = 10;
1546:     jac->as_amg_alpha_opts[1] = 1;
1547:     jac->as_amg_alpha_opts[2] = 6;
1548:     jac->as_amg_alpha_opts[3] = 6;
1549:     jac->as_amg_alpha_opts[4] = 4;
1550:     jac->as_amg_alpha_theta   = 0.25;
1551:     /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1552:     jac->as_amg_beta_opts[0] = 10;
1553:     jac->as_amg_beta_opts[1] = 1;
1554:     jac->as_amg_beta_opts[2] = 6;
1555:     jac->as_amg_beta_opts[3] = 6;
1556:     jac->as_amg_beta_opts[4] = 4;
1557:     jac->as_amg_beta_theta   = 0.25;
1558:     PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
1559:     PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1560:     PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1561:     PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
1562:     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1563:                                                                       jac->as_relax_times,
1564:                                                                       jac->as_relax_weight,
1565:                                                                       jac->as_omega));
1566:     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1567:                                                                      jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1568:                                                                      jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1569:                                                                      jac->as_amg_alpha_theta,
1570:                                                                      jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1571:                                                                      jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1572:     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1573:                                                                     jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1574:                                                                     jac->as_amg_beta_opts[2],       /* AMG relax_type */
1575:                                                                     jac->as_amg_beta_theta,
1576:                                                                     jac->as_amg_beta_opts[3],       /* AMG interp_type */
1577:                                                                     jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1578:     /* Zero conductivity */
1579:     jac->ams_beta_is_zero      = PETSC_FALSE;
1580:     jac->ams_beta_is_zero_part = PETSC_FALSE;
1581:     PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);
1582:     PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);
1583:     PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE_AMS);
1584:     PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetAlphaPoissonMatrix_C",PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS);
1585:     PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetBetaPoissonMatrix_C",PCHYPRESetBetaPoissonMatrix_HYPRE_AMS);
1586:     return(0);
1587:   }
1588:   PetscStrcmp("ads",jac->hypre_type,&flag);
1589:   if (flag) {
1590:     HYPRE_ADSCreate(&jac->hsolver);
1591:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_ADS;
1592:     pc->ops->view            = PCView_HYPRE_ADS;
1593:     jac->destroy             = HYPRE_ADSDestroy;
1594:     jac->setup               = HYPRE_ADSSetup;
1595:     jac->solve               = HYPRE_ADSSolve;
1596:     jac->setdgrad            = HYPRE_ADSSetDiscreteGradient;
1597:     jac->setdcurl            = HYPRE_ADSSetDiscreteCurl;
1598:     jac->setcoord            = HYPRE_ADSSetCoordinateVectors;
1599:     jac->coords[0]           = NULL;
1600:     jac->coords[1]           = NULL;
1601:     jac->coords[2]           = NULL;
1602:     jac->G                   = NULL;
1603:     jac->C                   = NULL;
1604:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */
1605:     jac->as_print           = 0;
1606:     jac->as_max_iter        = 1; /* used as a preconditioner */
1607:     jac->as_tol             = 0.; /* used as a preconditioner */
1608:     jac->ads_cycle_type     = 13;
1609:     /* Smoothing options */
1610:     jac->as_relax_type      = 2;
1611:     jac->as_relax_times     = 1;
1612:     jac->as_relax_weight    = 1.0;
1613:     jac->as_omega           = 1.0;
1614:     /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */
1615:     jac->ams_cycle_type       = 14;
1616:     jac->as_amg_alpha_opts[0] = 10;
1617:     jac->as_amg_alpha_opts[1] = 1;
1618:     jac->as_amg_alpha_opts[2] = 6;
1619:     jac->as_amg_alpha_opts[3] = 6;
1620:     jac->as_amg_alpha_opts[4] = 4;
1621:     jac->as_amg_alpha_theta   = 0.25;
1622:     /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1623:     jac->as_amg_beta_opts[0] = 10;
1624:     jac->as_amg_beta_opts[1] = 1;
1625:     jac->as_amg_beta_opts[2] = 6;
1626:     jac->as_amg_beta_opts[3] = 6;
1627:     jac->as_amg_beta_opts[4] = 4;
1628:     jac->as_amg_beta_theta   = 0.25;
1629:     PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
1630:     PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1631:     PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1632:     PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
1633:     PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1634:                                                                       jac->as_relax_times,
1635:                                                                       jac->as_relax_weight,
1636:                                                                       jac->as_omega));
1637:     PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type,             /* AMG coarsen type */
1638:                                                                 jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1639:                                                                 jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1640:                                                                 jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1641:                                                                 jac->as_amg_alpha_theta,
1642:                                                                 jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1643:                                                                 jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1644:     PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1645:                                                                 jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1646:                                                                 jac->as_amg_beta_opts[2],       /* AMG relax_type */
1647:                                                                 jac->as_amg_beta_theta,
1648:                                                                 jac->as_amg_beta_opts[3],       /* AMG interp_type */
1649:                                                                 jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1650:     PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);
1651:     PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);
1652:     PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",PCHYPRESetDiscreteCurl_HYPRE);
1653:     return(0);
1654:   }
1655:   PetscFree(jac->hypre_type);

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

1662: /*
1663:     It only gets here if the HYPRE type has not been set before the call to
1664:    ...SetFromOptions() which actually is most of the time
1665: */
1668: static PetscErrorCode PCSetFromOptions_HYPRE(PetscOptionItems *PetscOptionsObject,PC pc)
1669: {
1671:   PetscInt       indx;
1672:   const char     *type[] = {"pilut","parasails","boomeramg","ams","ads"};
1673:   PetscBool      flg;

1676:   PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");
1677:   PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"boomeramg",&indx,&flg);
1678:   if (flg) {
1679:     PCHYPRESetType_HYPRE(pc,type[indx]);
1680:   } else {
1681:     PCHYPRESetType_HYPRE(pc,"boomeramg");
1682:   }
1683:   if (pc->ops->setfromoptions) {
1684:     pc->ops->setfromoptions(PetscOptionsObject,pc);
1685:   }
1686:   PetscOptionsTail();
1687:   return(0);
1688: }

1692: /*@C
1693:      PCHYPRESetType - Sets which hypre preconditioner you wish to use

1695:    Input Parameters:
1696: +     pc - the preconditioner context
1697: -     name - either  pilut, parasails, boomeramg, ams, ads

1699:    Options Database Keys:
1700:    -pc_hypre_type - One of pilut, parasails, boomeramg, ams, ads

1702:    Level: intermediate

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

1707: @*/
1708: PetscErrorCode  PCHYPRESetType(PC pc,const char name[])
1709: {

1715:   PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));
1716:   return(0);
1717: }

1721: /*@C
1722:      PCHYPREGetType - Gets which hypre preconditioner you are using

1724:    Input Parameter:
1725: .     pc - the preconditioner context

1727:    Output Parameter:
1728: .     name - either  pilut, parasails, boomeramg, ams, ads

1730:    Level: intermediate

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

1735: @*/
1736: PetscErrorCode  PCHYPREGetType(PC pc,const char *name[])
1737: {

1743:   PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));
1744:   return(0);
1745: }

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

1750:    Options Database Keys:
1751: +   -pc_hypre_type - One of pilut, parasails, boomeramg, ams, ads
1752: -   Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX
1753:           preconditioner

1755:    Level: intermediate

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

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

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

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

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

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

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

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

1791: M*/

1795: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
1796: {
1797:   PC_HYPRE       *jac;

1801:   PetscNewLog(pc,&jac);

1803:   pc->data                = jac;
1804:   pc->ops->reset          = PCReset_HYPRE;
1805:   pc->ops->destroy        = PCDestroy_HYPRE;
1806:   pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
1807:   pc->ops->setup          = PCSetUp_HYPRE;
1808:   pc->ops->apply          = PCApply_HYPRE;
1809:   jac->comm_hypre         = MPI_COMM_NULL;
1810:   jac->hypre_type         = NULL;
1811:   jac->coords[0]          = NULL;
1812:   jac->coords[1]          = NULL;
1813:   jac->coords[2]          = NULL;
1814:   jac->constants[0]       = NULL;
1815:   jac->constants[1]       = NULL;
1816:   jac->constants[2]       = NULL;
1817:   jac->G                  = NULL;
1818:   jac->C                  = NULL;
1819:   jac->alpha_Poisson      = NULL;
1820:   jac->beta_Poisson       = NULL;
1821:   jac->setdim             = NULL;
1822:   jac->hmnull             = NULL;
1823:   jac->n_hmnull           = 0;
1824:   /* duplicate communicator for hypre */
1825:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1826:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);
1827:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);
1828:   return(0);
1829: }

1831: /* ---------------------------------------------------------------------------------------------------------------------------------*/

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

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

1840:   /* keep copy of PFMG options used so may view them */
1841:   PetscInt its;
1842:   double   tol;
1843:   PetscInt relax_type;
1844:   PetscInt rap_type;
1845:   PetscInt num_pre_relax,num_post_relax;
1846:   PetscInt max_levels;
1847: } PC_PFMG;

1851: PetscErrorCode PCDestroy_PFMG(PC pc)
1852: {
1854:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

1857:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1858:   MPI_Comm_free(&ex->hcomm);
1859:   PetscFree(pc->data);
1860:   return(0);
1861: }

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

1868: PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
1869: {
1871:   PetscBool      iascii;
1872:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

1875:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1876:   if (iascii) {
1877:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG preconditioning\n");
1878:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: max iterations %d\n",ex->its);
1879:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: tolerance %g\n",ex->tol);
1880:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
1881:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: RAP type %s\n",PFMGRAPType[ex->rap_type]);
1882:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
1883:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: max levels %d\n",ex->max_levels);
1884:   }
1885:   return(0);
1886: }


1891: PetscErrorCode PCSetFromOptions_PFMG(PetscOptionItems *PetscOptionsObject,PC pc)
1892: {
1894:   PC_PFMG        *ex = (PC_PFMG*) pc->data;
1895:   PetscBool      flg = PETSC_FALSE;

1898:   PetscOptionsHead(PetscOptionsObject,"PFMG options");
1899:   PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);
1900:   if (flg) {
1901:     int level=3;
1902:     PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level));
1903:   }
1904:   PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);
1905:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its));
1906:   PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
1907:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax));
1908:   PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
1909:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax));

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

1914:   PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);
1915:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol));
1916:   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);
1917:   PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type));
1918:   PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);
1919:   PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type));
1920:   PetscOptionsTail();
1921:   return(0);
1922: }

1926: PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
1927: {
1928:   PetscErrorCode    ierr;
1929:   PC_PFMG           *ex = (PC_PFMG*) pc->data;
1930:   PetscScalar       *yy;
1931:   const PetscScalar *xx;
1932:   PetscInt          ilower[3],iupper[3];
1933:   Mat_HYPREStruct   *mx = (Mat_HYPREStruct*)(pc->pmat->data);

1936:   PetscCitationsRegister(hypreCitation,&cite);
1937:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
1938:   iupper[0] += ilower[0] - 1;
1939:   iupper[1] += ilower[1] - 1;
1940:   iupper[2] += ilower[2] - 1;

1942:   /* copy x values over to hypre */
1943:   PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
1944:   VecGetArrayRead(x,&xx);
1945:   PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx));
1946:   VecRestoreArrayRead(x,&xx);
1947:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
1948:   PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx));

1950:   /* copy solution values back to PETSc */
1951:   VecGetArray(y,&yy);
1952:   PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy));
1953:   VecRestoreArray(y,&yy);
1954:   return(0);
1955: }

1959: 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)
1960: {
1961:   PC_PFMG        *jac = (PC_PFMG*)pc->data;
1963:   PetscInt       oits;

1966:   PetscCitationsRegister(hypreCitation,&cite);
1967:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its));
1968:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol));

1970:   PCApply_PFMG(pc,b,y);
1971:   PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
1972:   *outits = oits;
1973:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1974:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
1975:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol));
1976:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its));
1977:   return(0);
1978: }


1983: PetscErrorCode PCSetUp_PFMG(PC pc)
1984: {
1985:   PetscErrorCode  ierr;
1986:   PC_PFMG         *ex = (PC_PFMG*) pc->data;
1987:   Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
1988:   PetscBool       flg;

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

1994:   /* create the hypre solver object and set its information */
1995:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1996:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1997:   PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
1998:   PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver));
1999:   return(0);
2000: }


2003: /*MC
2004:      PCPFMG - the hypre PFMG multigrid solver

2006:    Level: advanced

2008:    Options Database:
2009: + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner
2010: . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2011: . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2012: . -pc_pfmg_tol <tol> tolerance of PFMG
2013: . -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
2014: - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin

2016:    Notes:  This is for CELL-centered descretizations

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

2021: .seealso:  PCMG, MATHYPRESTRUCT
2022: M*/

2026: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
2027: {
2029:   PC_PFMG        *ex;

2032:   PetscNew(&ex); \
2033:   pc->data = ex;

2035:   ex->its            = 1;
2036:   ex->tol            = 1.e-8;
2037:   ex->relax_type     = 1;
2038:   ex->rap_type       = 0;
2039:   ex->num_pre_relax  = 1;
2040:   ex->num_post_relax = 1;
2041:   ex->max_levels     = 0;

2043:   pc->ops->setfromoptions  = PCSetFromOptions_PFMG;
2044:   pc->ops->view            = PCView_PFMG;
2045:   pc->ops->destroy         = PCDestroy_PFMG;
2046:   pc->ops->apply           = PCApply_PFMG;
2047:   pc->ops->applyrichardson = PCApplyRichardson_PFMG;
2048:   pc->ops->setup           = PCSetUp_PFMG;

2050:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
2051:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
2052:   return(0);
2053: }

2055: /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/

2057: /* we know we are working with a HYPRE_SStructMatrix */
2058: typedef struct {
2059:   MPI_Comm            hcomm;       /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
2060:   HYPRE_SStructSolver ss_solver;

2062:   /* keep copy of SYSPFMG options used so may view them */
2063:   PetscInt its;
2064:   double   tol;
2065:   PetscInt relax_type;
2066:   PetscInt num_pre_relax,num_post_relax;
2067: } PC_SysPFMG;

2071: PetscErrorCode PCDestroy_SysPFMG(PC pc)
2072: {
2074:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

2077:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2078:   MPI_Comm_free(&ex->hcomm);
2079:   PetscFree(pc->data);
2080:   return(0);
2081: }

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

2087: PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
2088: {
2090:   PetscBool      iascii;
2091:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

2094:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2095:   if (iascii) {
2096:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG preconditioning\n");
2097:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: max iterations %d\n",ex->its);
2098:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: tolerance %g\n",ex->tol);
2099:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
2100:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
2101:   }
2102:   return(0);
2103: }


2108: PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptionItems *PetscOptionsObject,PC pc)
2109: {
2111:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;
2112:   PetscBool      flg = PETSC_FALSE;

2115:   PetscOptionsHead(PetscOptionsObject,"SysPFMG options");
2116:   PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);
2117:   if (flg) {
2118:     int level=3;
2119:     PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level));
2120:   }
2121:   PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);
2122:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its));
2123:   PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
2124:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax));
2125:   PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
2126:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax));

2128:   PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);
2129:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol));
2130:   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);
2131:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type));
2132:   PetscOptionsTail();
2133:   return(0);
2134: }

2138: PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
2139: {
2140:   PetscErrorCode    ierr;
2141:   PC_SysPFMG        *ex = (PC_SysPFMG*) pc->data;
2142:   PetscScalar       *yy;
2143:   const PetscScalar *xx;
2144:   PetscInt          ilower[3],iupper[3];
2145:   Mat_HYPRESStruct  *mx     = (Mat_HYPRESStruct*)(pc->pmat->data);
2146:   PetscInt          ordering= mx->dofs_order;
2147:   PetscInt          nvars   = mx->nvars;
2148:   PetscInt          part    = 0;
2149:   PetscInt          size;
2150:   PetscInt          i;

2153:   PetscCitationsRegister(hypreCitation,&cite);
2154:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
2155:   iupper[0] += ilower[0] - 1;
2156:   iupper[1] += ilower[1] - 1;
2157:   iupper[2] += ilower[2] - 1;

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

2162:   /* copy x values over to hypre for variable ordering */
2163:   if (ordering) {
2164:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2165:     VecGetArrayRead(x,&xx);
2166:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)xx+(size*i)));
2167:     VecRestoreArrayRead(x,&xx);
2168:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2169:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
2170:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

2172:     /* copy solution values back to PETSc */
2173:     VecGetArray(y,&yy);
2174:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i)));
2175:     VecRestoreArray(y,&yy);
2176:   } else {      /* nodal ordering must be mapped to variable ordering for sys_pfmg */
2177:     PetscScalar *z;
2178:     PetscInt    j, k;

2180:     PetscMalloc1(nvars*size,&z);
2181:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2182:     VecGetArrayRead(x,&xx);

2184:     /* transform nodal to hypre's variable ordering for sys_pfmg */
2185:     for (i= 0; i< size; i++) {
2186:       k= i*nvars;
2187:       for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
2188:     }
2189:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
2190:     VecRestoreArrayRead(x,&xx);
2191:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2192:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

2194:     /* copy solution values back to PETSc */
2195:     VecGetArray(y,&yy);
2196:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
2197:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
2198:     for (i= 0; i< size; i++) {
2199:       k= i*nvars;
2200:       for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
2201:     }
2202:     VecRestoreArray(y,&yy);
2203:     PetscFree(z);
2204:   }
2205:   return(0);
2206: }

2210: 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)
2211: {
2212:   PC_SysPFMG     *jac = (PC_SysPFMG*)pc->data;
2214:   PetscInt       oits;

2217:   PetscCitationsRegister(hypreCitation,&cite);
2218:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its));
2219:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol));
2220:   PCApply_SysPFMG(pc,b,y);
2221:   PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,(HYPRE_Int *)&oits));
2222:   *outits = oits;
2223:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2224:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
2225:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol));
2226:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its));
2227:   return(0);
2228: }


2233: PetscErrorCode PCSetUp_SysPFMG(PC pc)
2234: {
2235:   PetscErrorCode   ierr;
2236:   PC_SysPFMG       *ex = (PC_SysPFMG*) pc->data;
2237:   Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
2238:   PetscBool        flg;

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

2244:   /* create the hypre sstruct solver object and set its information */
2245:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2246:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2247:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver));
2248:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2249:   return(0);
2250: }


2253: /*MC
2254:      PCSysPFMG - the hypre SysPFMG multigrid solver

2256:    Level: advanced

2258:    Options Database:
2259: + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner
2260: . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2261: . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2262: . -pc_syspfmg_tol <tol> tolerance of SysPFMG
2263: . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel

2265:    Notes:  This is for CELL-centered descretizations

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

2271: .seealso:  PCMG, MATHYPRESSTRUCT
2272: M*/

2276: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
2277: {
2279:   PC_SysPFMG     *ex;

2282:   PetscNew(&ex); \
2283:   pc->data = ex;

2285:   ex->its            = 1;
2286:   ex->tol            = 1.e-8;
2287:   ex->relax_type     = 1;
2288:   ex->num_pre_relax  = 1;
2289:   ex->num_post_relax = 1;

2291:   pc->ops->setfromoptions  = PCSetFromOptions_SysPFMG;
2292:   pc->ops->view            = PCView_SysPFMG;
2293:   pc->ops->destroy         = PCDestroy_SysPFMG;
2294:   pc->ops->apply           = PCApply_SysPFMG;
2295:   pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
2296:   pc->ops->setup           = PCSetUp_SysPFMG;

2298:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
2299:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2300:   return(0);
2301: }