Actual source code: hypre.c

petsc-master 2016-12-09
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: /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */
 10:  #include <petsc/private/matimpl.h>
 11:  #include <../src/vec/vec/impls/hypre/vhyp.h>
 12:  #include <../src/mat/impls/hypre/mhypre.h>
 13:  #include <../src/dm/impls/da/hypre/mhyp.h>
 14: #include <_hypre_parcsr_ls.h>

 16: static PetscBool cite = PETSC_FALSE;
 17: 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";

 19: /*
 20:    Private context (data structure) for the  preconditioner.
 21: */
 22: typedef struct {
 23:   HYPRE_Solver   hsolver;
 24:   Mat            hpmat; /* MatHYPRE */

 26:   HYPRE_Int (*destroy)(HYPRE_Solver);
 27:   HYPRE_Int (*solve)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
 28:   HYPRE_Int (*setup)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);

 30:   MPI_Comm comm_hypre;
 31:   char     *hypre_type;

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

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

 40:   /* options for ParaSails */
 41:   PetscInt nlevels;
 42:   double   threshhold;
 43:   double   filter;
 44:   PetscInt sym;
 45:   double   loadbal;
 46:   PetscInt logging;
 47:   PetscInt ruse;
 48:   PetscInt symt;

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

 53:   /* options for BoomerAMG */
 54:   PetscInt  cycletype;
 55:   PetscInt  maxlevels;
 56:   double    strongthreshold;
 57:   double    maxrowsum;
 58:   PetscInt  gridsweeps[3];
 59:   PetscInt  coarsentype;
 60:   PetscInt  measuretype;
 61:   PetscInt  smoothtype;
 62:   PetscInt  smoothnumlevels;
 63:   PetscInt  eu_level;   /* Number of levels for ILU(k) in Euclid */
 64:   double    eu_droptolerance; /* Drop tolerance for ILU(k) in Euclid */
 65:   PetscInt  eu_bj;      /* Defines use of Block Jacobi ILU in Euclid */
 66:   PetscInt  relaxtype[3];
 67:   double    relaxweight;
 68:   double    outerrelaxweight;
 69:   PetscInt  relaxorder;
 70:   double    truncfactor;
 71:   PetscBool applyrichardson;
 72:   PetscInt  pmax;
 73:   PetscInt  interptype;
 74:   PetscInt  agg_nl;
 75:   PetscInt  agg_num_paths;
 76:   PetscInt  nodal_coarsen;
 77:   PetscBool nodal_relax;
 78:   PetscInt  nodal_relax_levels;

 80:   PetscInt  nodal_coarsening;
 81:   PetscInt  vec_interp_variant;
 82:   HYPRE_IJVector  *hmnull;
 83:   HYPRE_ParVector *phmnull;  /* near null space passed to hypre */
 84:   PetscInt        n_hmnull;
 85:   Vec             hmnull_constant;
 86:   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 */

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

103:   /* additional data */
104:   Mat G;             /* MatHYPRE */
105:   Mat C;             /* MatHYPRE */
106:   Mat alpha_Poisson; /* MatHYPRE */
107:   Mat beta_Poisson;  /* MatHYPRE */

109:   /* extra information for AMS */
110:   PetscInt       dim; /* geometrical dimension */
111:   HYPRE_IJVector coords[3];
112:   HYPRE_IJVector constants[3];
113:   PetscBool      ams_beta_is_zero;
114:   PetscBool      ams_beta_is_zero_part;
115:   PetscInt       ams_proj_freq;
116: } PC_HYPRE;

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

125:   *hsolver = jac->hsolver;
126:   return(0);
127: }

131: static PetscErrorCode PCSetUp_HYPRE(PC pc)
132: {
133:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
134:   Mat_HYPRE          *hjac;
135:   HYPRE_ParCSRMatrix hmat;
136:   HYPRE_ParVector    bv,xv;
137:   PetscBool          ishypre;
138:   PetscErrorCode     ierr;

141:   if (!jac->hypre_type) {
142:     PCHYPRESetType(pc,"boomeramg");
143:   }

145:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRE,&ishypre);
146:   if (!ishypre) {
147:     MatReuse reuse;
148:     if (pc->setupcalled) reuse = MAT_REUSE_MATRIX;
149:     else reuse = MAT_INITIAL_MATRIX;
150:     MatConvert(pc->pmat,MATHYPRE,reuse,&jac->hpmat);
151:   } else {
152:     PetscObjectReference((PetscObject)pc->pmat);
153:     MatDestroy(&jac->hpmat);
154:     jac->hpmat = pc->pmat;
155:   }
156:   hjac = (Mat_HYPRE*)(jac->hpmat->data);

158:   /* special case for BoomerAMG */
159:   if (jac->setup == HYPRE_BoomerAMGSetup) {
160:     MatNullSpace    mnull;
161:     PetscBool       has_const;
162:     PetscInt        bs,nvec,i;
163:     const Vec       *vecs;
164:     PetscScalar     *petscvecarray;
165: 
166:     MatGetBlockSize(pc->pmat,&bs);
167:     if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
168:     MatGetNearNullSpace(pc->mat, &mnull);
169:     if (mnull) {
170:       MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);
171:       PetscMalloc1(nvec+1,&jac->hmnull);
172:       PetscMalloc1(nvec+1,&jac->hmnull_hypre_data_array);
173:       PetscMalloc1(nvec+1,&jac->phmnull);
174:       for (i=0; i<nvec; i++) {
175:         VecHYPRE_IJVectorCreate(vecs[i],&jac->hmnull[i]);
176:         VecGetArrayRead(vecs[i],(const PetscScalar **)&petscvecarray);
177:         VecHYPRE_ParVectorReplacePointer(jac->hmnull[i],petscvecarray,jac->hmnull_hypre_data_array[i]);
178:         VecRestoreArrayRead(vecs[i],(const PetscScalar **)&petscvecarray);
179:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[i],(void**)&jac->phmnull[i]));
180:       }
181:       if (has_const) {
182:         MatCreateVecs(pc->pmat,&jac->hmnull_constant,NULL);
183:         VecSet(jac->hmnull_constant,1);
184:         VecNormalize(jac->hmnull_constant,NULL);
185:         VecHYPRE_IJVectorCreate(jac->hmnull_constant,&jac->hmnull[nvec]);
186:         VecGetArrayRead(jac->hmnull_constant,(const PetscScalar **)&petscvecarray);
187:         VecHYPRE_ParVectorReplacePointer(jac->hmnull[nvec],petscvecarray,jac->hmnull_hypre_data_array[nvec]);
188:         VecRestoreArrayRead(jac->hmnull_constant,(const PetscScalar **)&petscvecarray);
189:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[nvec],(void**)&jac->phmnull[nvec]));
190:         nvec++;
191:       }
192:       PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVectors,(jac->hsolver,nvec,jac->phmnull));
193:       jac->n_hmnull = nvec;
194:     }
195:   }

197:   /* special case for AMS */
198:   if (jac->setup == HYPRE_AMSSetup) {
199:     Mat_HYPRE          *hm;
200:     HYPRE_ParCSRMatrix parcsr;
201:     if (!jac->coords[0] && !jac->constants[0]) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs either the coordinate vectors via PCSetCoordinates() or the edge constant vectors via PCHYPRESetEdgeConstantVectors()");
202:     if (jac->dim) {
203:       PetscStackCallStandard(HYPRE_AMSSetDimension,(jac->hsolver,jac->dim));
204:     }
205:     if (jac->constants[0]) {
206:       HYPRE_ParVector ozz,zoz,zzo = NULL;
207:       PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[0],(void**)(&ozz)));
208:       PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[1],(void**)(&zoz)));
209:       if (jac->constants[2]) {
210:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[2],(void**)(&zzo)));
211:       }
212:       PetscStackCallStandard(HYPRE_AMSSetEdgeConstantVectors,(jac->hsolver,ozz,zoz,zzo));
213:     }
214:     if (jac->coords[0]) {
215:       HYPRE_ParVector coords[3];
216:       coords[0] = NULL;
217:       coords[1] = NULL;
218:       coords[2] = NULL;
219:       if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0],(void**)(&coords[0])));
220:       if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1],(void**)(&coords[1])));
221:       if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2],(void**)(&coords[2])));
222:       PetscStackCallStandard(HYPRE_AMSSetCoordinateVectors,(jac->hsolver,coords[0],coords[1],coords[2]));
223:     }
224:     if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient");
225:     hm = (Mat_HYPRE*)(jac->G->data);
226:     PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
227:     PetscStackCallStandard(HYPRE_AMSSetDiscreteGradient,(jac->hsolver,parcsr));
228:     if (jac->alpha_Poisson) {
229:       hm = (Mat_HYPRE*)(jac->alpha_Poisson->data);
230:       PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
231:       PetscStackCallStandard(HYPRE_AMSSetAlphaPoissonMatrix,(jac->hsolver,parcsr));
232:     }
233:     if (jac->ams_beta_is_zero) {
234:       PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,NULL));
235:     } else if (jac->beta_Poisson) {
236:       hm = (Mat_HYPRE*)(jac->beta_Poisson->data);
237:       PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
238:       PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,parcsr));
239:     }
240:   }
241:   /* special case for ADS */
242:   if (jac->setup == HYPRE_ADSSetup) {
243:     Mat_HYPRE          *hm;
244:     HYPRE_ParCSRMatrix parcsr;
245:     if (!jac->coords[0]) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs the coordinate vectors via PCSetCoordinates()");
246:     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");
247:     if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient");
248:     if (!jac->C) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs the discrete curl operator via PCHYPRESetDiscreteGradient");
249:     if (jac->coords[0]) {
250:       HYPRE_ParVector coords[3];
251:       coords[0] = NULL;
252:       coords[1] = NULL;
253:       coords[2] = NULL;
254:       if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0],(void**)(&coords[0])));
255:       if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1],(void**)(&coords[1])));
256:       if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2],(void**)(&coords[2])));
257:       PetscStackCallStandard(HYPRE_ADSSetCoordinateVectors,(jac->hsolver,coords[0],coords[1],coords[2]));
258:     }
259:     hm = (Mat_HYPRE*)(jac->G->data);
260:     PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
261:     PetscStackCallStandard(HYPRE_ADSSetDiscreteGradient,(jac->hsolver,parcsr));
262:     hm = (Mat_HYPRE*)(jac->C->data);
263:     PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
264:     PetscStackCallStandard(HYPRE_ADSSetDiscreteCurl,(jac->hsolver,parcsr));
265:   }
266:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
267:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b,(void**)&bv));
268:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x,(void**)&xv));
269:   PetscStackCall("HYPRE_SetupXXX",(*jac->setup)(jac->hsolver,hmat,bv,xv););
270:   return(0);
271: }

275: static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
276: {
277:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
278:   Mat_HYPRE          *hjac = (Mat_HYPRE*)(jac->hpmat->data);
279:   PetscErrorCode     ierr;
280:   HYPRE_ParCSRMatrix hmat;
281:   PetscScalar        *xv;
282:   const PetscScalar  *bv,*sbv;
283:   HYPRE_ParVector    jbv,jxv;
284:   PetscScalar        *sxv;
285:   PetscInt           hierr;

288:   PetscCitationsRegister(hypreCitation,&cite);
289:   if (!jac->applyrichardson) {VecSet(x,0.0);}
290:   VecGetArrayRead(b,&bv);
291:   VecGetArray(x,&xv);
292:   VecHYPRE_ParVectorReplacePointer(hjac->b,(PetscScalar*)bv,sbv);
293:   VecHYPRE_ParVectorReplacePointer(hjac->x,xv,sxv);

295:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
296:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b,(void**)&jbv));
297:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x,(void**)&jxv));
298:   PetscStackCall("Hypre solve",h(*jac->solve)(jac->hsolver,hmat,jbv,jxv);
299:                                if (hierr && hierr != HYPRE_ERROR_CONV) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
300:                                if (hierr) hypre__global_error = 0;);

302:   if (jac->setup == HYPRE_AMSSetup && jac->ams_beta_is_zero_part) {
303:     PetscStackCallStandard(HYPRE_AMSProjectOutGradients,(jac->hsolver,jxv));
304:   }
305:   VecHYPRE_ParVectorReplacePointer(hjac->b,(PetscScalar*)sbv,bv);
306:   VecHYPRE_ParVectorReplacePointer(hjac->x,sxv,xv);
307:   VecRestoreArray(x,&xv);
308:   VecRestoreArrayRead(b,&bv);
309:   return(0);
310: }

314: static PetscErrorCode PCReset_HYPRE(PC pc)
315: {
316:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;

320:   MatDestroy(&jac->hpmat);
321:   MatDestroy(&jac->G);
322:   MatDestroy(&jac->C);
323:   MatDestroy(&jac->alpha_Poisson);
324:   MatDestroy(&jac->beta_Poisson);
325:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0])); jac->coords[0] = NULL;
326:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1])); jac->coords[1] = NULL;
327:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2])); jac->coords[2] = NULL;
328:   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0])); jac->constants[0] = NULL;
329:   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1])); jac->constants[1] = NULL;
330:   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2])); jac->constants[2] = NULL;
331:   if (jac->n_hmnull && jac->hmnull) {
332:     PetscInt                 i;
333:     PETSC_UNUSED PetscScalar *petscvecarray;

335:     for (i=0; i<jac->n_hmnull; i++) {
336:       VecHYPRE_ParVectorReplacePointer(jac->hmnull[i],jac->hmnull_hypre_data_array[i],petscvecarray);
337:       PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->hmnull[i]));
338:     }
339:     PetscFree(jac->hmnull);
340:     PetscFree(jac->hmnull_hypre_data_array);
341:     PetscFree(jac->phmnull);
342:     VecDestroy(&jac->hmnull_constant);
343:   }
344:   jac->ams_beta_is_zero = PETSC_FALSE;
345:   jac->dim = 0;
346:   return(0);
347: }

351: static PetscErrorCode PCDestroy_HYPRE(PC pc)
352: {
353:   PC_HYPRE                 *jac = (PC_HYPRE*)pc->data;
354:   PetscErrorCode           ierr;

357:   PCReset_HYPRE(pc);
358:   if (jac->destroy) PetscStackCall("HYPRE_DestroyXXX",(*jac->destroy)(jac->hsolver););
359:   PetscFree(jac->hypre_type);
360:   if (jac->comm_hypre != MPI_COMM_NULL) { MPI_Comm_free(&(jac->comm_hypre));}
361:   PetscFree(pc->data);

363:   PetscObjectChangeTypeName((PetscObject)pc,0);
364:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);
365:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);
366:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetCoordinates_C",NULL);
367:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",NULL);
368:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",NULL);
369:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetConstantEdgeVectors_C",NULL);
370:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",NULL);
371:   return(0);
372: }

374: /* --------------------------------------------------------------------------------------------*/
377: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PetscOptionItems *PetscOptionsObject,PC pc)
378: {
379:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
381:   PetscBool      flag;

384:   PetscOptionsHead(PetscOptionsObject,"HYPRE Pilut Options");
385:   PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
386:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter));
387:   PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
388:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol));
389:   PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
390:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize));
391:   PetscOptionsTail();
392:   return(0);
393: }

397: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
398: {
399:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
401:   PetscBool      iascii;

404:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
405:   if (iascii) {
406:     PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut preconditioning\n");
407:     if (jac->maxiter != PETSC_DEFAULT) {
408:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: maximum number of iterations %d\n",jac->maxiter);
409:     } else {
410:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default maximum number of iterations \n");
411:     }
412:     if (jac->tol != PETSC_DEFAULT) {
413:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: drop tolerance %g\n",(double)jac->tol);
414:     } else {
415:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default drop tolerance \n");
416:     }
417:     if (jac->factorrowsize != PETSC_DEFAULT) {
418:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: factor row size %d\n",jac->factorrowsize);
419:     } else {
420:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default factor row size \n");
421:     }
422:   }
423:   return(0);
424: }

426: /* --------------------------------------------------------------------------------------------*/

430: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
431: {
432:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
433:   Mat_HYPRE          *hjac = (Mat_HYPRE*)(jac->hpmat->data);
434:   PetscErrorCode     ierr;
435:   HYPRE_ParCSRMatrix hmat;
436:   PetscScalar        *xv;
437:   const PetscScalar  *bv;
438:   HYPRE_ParVector    jbv,jxv;
439:   PetscScalar        *sbv,*sxv;
440:   PetscInt           hierr;

443:   PetscCitationsRegister(hypreCitation,&cite);
444:   VecSet(x,0.0);
445:   VecGetArrayRead(b,&bv);
446:   VecGetArray(x,&xv);
447:   VecHYPRE_ParVectorReplacePointer(hjac->b,(PetscScalar*)bv,sbv);
448:   VecHYPRE_ParVectorReplacePointer(hjac->x,xv,sxv);

450:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
451:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b,(void**)&jbv));
452:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x,(void**)&jxv));

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

459:   VecHYPRE_ParVectorReplacePointer(hjac->b,sbv,bv);
460:   VecHYPRE_ParVectorReplacePointer(hjac->x,sxv,xv);
461:   VecRestoreArray(x,&xv);
462:   VecRestoreArrayRead(b,&bv);
463:   return(0);
464: }

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

469: static const char *HYPREBoomerAMGCycleType[]   = {"","V","W"};
470: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
471: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
472: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
473: static const char *HYPREBoomerAMGSmoothType[]   = {"Schwarz-smoothers","Pilut","ParaSails","Euclid"};
474: static const char *HYPREBoomerAMGRelaxType[]   = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi",
475:                                                   "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi",
476:                                                   "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination",
477:                                                   "" /* 10 */, "" /* 11 */, "" /* 12 */, "" /* 13 */, "" /* 14 */,
478:                                                   "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"};
479: static const char *HYPREBoomerAMGInterpType[]  = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
480:                                                   "ext+i-cc", "standard", "standard-wts", "", "", "FF", "FF1"};
483: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PetscOptionItems *PetscOptionsObject,PC pc)
484: {
485:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
487:   PetscInt       n,indx,level;
488:   PetscBool      flg, tmp_truth;
489:   double         tmpdbl, twodbl[2];

492:   PetscOptionsHead(PetscOptionsObject,"HYPRE BoomerAMG Options");
493:   PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);
494:   if (flg) {
495:     jac->cycletype = indx+1;
496:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
497:   }
498:   PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
499:   if (flg) {
500:     if (jac->maxlevels < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
501:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
502:   }
503:   PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
504:   if (flg) {
505:     if (jac->maxiter < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
506:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
507:   }
508:   PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);
509:   if (flg) {
510:     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);
511:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
512:   }

514:   PetscOptionsScalar("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);
515:   if (flg) {
516:     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);
517:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
518:   }

520:   PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);
521:   if (flg) {
522:     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);
523:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
524:   }

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

530:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
531:   }


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

538:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
539:   }


542:   PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
543:   if (flg) {
544:     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);
545:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
546:   }
547:   PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
548:   if (flg) {
549:     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);
550:     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);
551:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
552:   }

554:   /* Grid sweeps */
555:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);
556:   if (flg) {
557:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx));
558:     /* modify the jac structure so we can view the updated options with PC_View */
559:     jac->gridsweeps[0] = indx;
560:     jac->gridsweeps[1] = indx;
561:     /*defaults coarse to 1 */
562:     jac->gridsweeps[2] = 1;
563:   }

565:   PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen","Use a nodal based coarsening 1-6","HYPRE_BoomerAMGSetNodal",jac->nodal_coarsening, &jac->nodal_coarsening,&flg);
566:   if (flg) {
567:     PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,jac->nodal_coarsening));
568:   }

570:   PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_variant","Variant of algorithm 1-3","HYPRE_BoomerAMGSetInterpVecVariant",jac->vec_interp_variant, &jac->vec_interp_variant,&flg);
571:   if (flg) {
572:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecVariant,(jac->hsolver,jac->vec_interp_variant));
573:   }

575:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);
576:   if (flg) {
577:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1));
578:     jac->gridsweeps[0] = indx;
579:   }
580:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);
581:   if (flg) {
582:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2));
583:     jac->gridsweeps[1] = indx;
584:   }
585:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);
586:   if (flg) {
587:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3));
588:     jac->gridsweeps[2] = indx;
589:   }

591:   /* Smooth type */
592:   PetscOptionsEList("-pc_hypre_boomeramg_smooth_type","Enable more complex smoothers","None",HYPREBoomerAMGSmoothType,ALEN(HYPREBoomerAMGSmoothType),HYPREBoomerAMGSmoothType[0],&indx,&flg);
593:   if (flg) {
594:     jac->smoothtype = indx;
595:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,indx+6));
596:     jac->smoothnumlevels = 25;
597:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,25));
598:   }

600:   /* Number of smoothing levels */
601:   PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_levels","Number of levels on which more complex smoothers are used","None",25,&indx,&flg);
602:   if (flg && (jac->smoothtype != -1)) {
603:     jac->smoothnumlevels = indx;
604:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,indx));
605:   }

607:   /* Number of levels for ILU(k) for Euclid */
608:   PetscOptionsInt("-pc_hypre_boomeramg_eu_level","Number of levels for ILU(k) in Euclid smoother","None",0,&indx,&flg);
609:   if (flg && (jac->smoothtype == 3)) {
610:     jac->eu_level = indx;
611:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,indx));
612:   }

614:   /* Filter for ILU(k) for Euclid */
615:   double droptolerance;
616:   PetscOptionsScalar("-pc_hypre_boomeramg_eu_droptolerance","Drop tolerance for ILU(k) in Euclid smoother","None",0,&droptolerance,&flg);
617:   if (flg && (jac->smoothtype == 3)) {
618:     jac->eu_droptolerance = droptolerance;
619:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,droptolerance));
620:   }

622:   /* Use Block Jacobi ILUT for Euclid */
623:   PetscOptionsBool("-pc_hypre_boomeramg_eu_bj", "Use Block Jacobi for ILU in Euclid smoother?", "None", PETSC_FALSE, &tmp_truth, &flg);
624:   if (flg && (jac->smoothtype == 3)) {
625:     jac->eu_bj = tmp_truth;
626:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuBJ,(jac->hsolver,jac->eu_bj));
627:   }

629:   /* Relax type */
630:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
631:   if (flg) {
632:     jac->relaxtype[0] = jac->relaxtype[1]  = indx;
633:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx));
634:     /* by default, coarse type set to 9 */
635:     jac->relaxtype[2] = 9;

637:   }
638:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
639:   if (flg) {
640:     jac->relaxtype[0] = indx;
641:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1));
642:   }
643:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
644:   if (flg) {
645:     jac->relaxtype[1] = indx;
646:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2));
647:   }
648:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);
649:   if (flg) {
650:     jac->relaxtype[2] = indx;
651:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3));
652:   }

654:   /* Relaxation Weight */
655:   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);
656:   if (flg) {
657:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl));
658:     jac->relaxweight = tmpdbl;
659:   }

661:   n         = 2;
662:   twodbl[0] = twodbl[1] = 1.0;
663:   PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
664:   if (flg) {
665:     if (n == 2) {
666:       indx =  (int)PetscAbsReal(twodbl[1]);
667:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx));
668:     } 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);
669:   }

671:   /* Outer relaxation Weight */
672:   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);
673:   if (flg) {
674:     PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl));
675:     jac->outerrelaxweight = tmpdbl;
676:   }

678:   n         = 2;
679:   twodbl[0] = twodbl[1] = 1.0;
680:   PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
681:   if (flg) {
682:     if (n == 2) {
683:       indx =  (int)PetscAbsReal(twodbl[1]);
684:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx));
685:     } 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);
686:   }

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

691:   if (flg && tmp_truth) {
692:     jac->relaxorder = 0;
693:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
694:   }
695:   PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);
696:   if (flg) {
697:     jac->measuretype = indx;
698:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
699:   }
700:   /* update list length 3/07 */
701:   PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);
702:   if (flg) {
703:     jac->coarsentype = indx;
704:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
705:   }

707:   /* new 3/07 */
708:   PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);
709:   if (flg) {
710:     jac->interptype = indx;
711:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
712:   }

714:   PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
715:   if (flg) {
716:     level = 3;
717:     PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);

719:     jac->printstatistics = PETSC_TRUE;
720:     PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level));
721:   }

723:   PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);
724:   if (flg) {
725:     level = 3;
726:     PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);

728:     jac->printstatistics = PETSC_TRUE;
729:     PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level));
730:   }

732:   PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);
733:   if (flg && tmp_truth) {
734:     PetscInt tmp_int;
735:     PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);
736:     if (flg) jac->nodal_relax_levels = tmp_int;
737:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6));
738:     PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1));
739:     PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0));
740:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels));
741:   }

743:   PetscOptionsTail();
744:   return(0);
745: }

749: 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)
750: {
751:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
753:   PetscInt       oits;

756:   PetscCitationsRegister(hypreCitation,&cite);
757:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter));
758:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol));
759:   jac->applyrichardson = PETSC_TRUE;
760:   PCApply_HYPRE(pc,b,y);
761:   jac->applyrichardson = PETSC_FALSE;
762:   PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
763:   *outits = oits;
764:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
765:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
766:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
767:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
768:   return(0);
769: }


774: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
775: {
776:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
778:   PetscBool      iascii;

781:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
782:   if (iascii) {
783:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG preconditioning\n");
784:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);
785:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of levels %d\n",jac->maxlevels);
786:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of iterations PER hypre call %d\n",jac->maxiter);
787:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Convergence tolerance PER hypre call %g\n",(double)jac->tol);
788:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Threshold for strong coupling %g\n",(double)jac->strongthreshold);
789:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation truncation factor %g\n",(double)jac->truncfactor);
790:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation: max elements per row %d\n",jac->pmax);
791:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Number of levels of aggressive coarsening %d\n",jac->agg_nl);
792:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Number of paths for aggressive coarsening %d\n",jac->agg_num_paths);

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

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

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

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

807:     if (jac->relaxorder) {
808:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Using CF-relaxation\n");
809:     } else {
810:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Not using CF-relaxation\n");
811:     }
812:     if (jac->smoothtype!=-1) {
813:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Smooth type          %s\n",HYPREBoomerAMGSmoothType[jac->smoothtype]);
814:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Smooth num levels    %d\n",jac->smoothnumlevels);
815:     } else {
816:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Not using more complex smoothers.\n");
817:     }
818:     if (jac->smoothtype==3) {
819:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Euclid ILU(k) levels %d\n",jac->eu_level);
820:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Euclid ILU(k) drop tolerance %g\n",jac->eu_droptolerance);
821:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Euclid ILU use Block-Jacobi? %d\n",jac->eu_bj);
822:     }
823:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Measure type        %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
824:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Coarsen type        %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
825:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation type  %s\n",HYPREBoomerAMGInterpType[jac->interptype]);
826:     if (jac->nodal_coarsening) {
827:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Using nodal coarsening (with HYPRE_BOOMERAMGSetNodal() %D\n",jac->nodal_coarsening);
828:     }
829:     if (jac->vec_interp_variant) {
830:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: HYPRE_BoomerAMGSetInterpVecVariant() %D\n",jac->vec_interp_variant);
831:     }
832:     if (jac->nodal_relax) {
833:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Using nodal relaxation via Schwarz smoothing on levels %d\n",jac->nodal_relax_levels);
834:     }
835:   }
836:   return(0);
837: }

839: /* --------------------------------------------------------------------------------------------*/
842: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PetscOptionItems *PetscOptionsObject,PC pc)
843: {
844:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
846:   PetscInt       indx;
847:   PetscBool      flag;
848:   const char     *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};

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

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

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

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

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

868:   PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);
869:   if (flag) {
870:     jac->symt = indx;
871:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
872:   }

874:   PetscOptionsTail();
875:   return(0);
876: }

880: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
881: {
882:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
884:   PetscBool      iascii;
885:   const char     *symt = 0;;

888:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
889:   if (iascii) {
890:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails preconditioning\n");
891:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: nlevels %d\n",jac->nlevels);
892:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: threshold %g\n",(double)jac->threshhold);
893:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: filter %g\n",(double)jac->filter);
894:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: load balance %g\n",(double)jac->loadbal);
895:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: reuse nonzero structure %s\n",PetscBools[jac->ruse]);
896:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: print info to screen %s\n",PetscBools[jac->logging]);
897:     if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
898:     else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
899:     else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
900:     else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
901:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: %s\n",symt);
902:   }
903:   return(0);
904: }
905: /* --------------------------------------------------------------------------------------------*/
908: static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PetscOptionItems *PetscOptionsObject,PC pc)
909: {
910:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
912:   PetscInt       n;
913:   PetscBool      flag,flag2,flag3,flag4;

916:   PetscOptionsHead(PetscOptionsObject,"HYPRE AMS Options");
917:   PetscOptionsInt("-pc_hypre_ams_print_level","Debugging output level for AMS","None",jac->as_print,&jac->as_print,&flag);
918:   if (flag) PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
919:   PetscOptionsInt("-pc_hypre_ams_max_iter","Maximum number of AMS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
920:   if (flag) PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
921:   PetscOptionsInt("-pc_hypre_ams_cycle_type","Cycle type for AMS multigrid","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag);
922:   if (flag) PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
923:   PetscOptionsReal("-pc_hypre_ams_tol","Error tolerance for AMS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
924:   if (flag) PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
925:   PetscOptionsInt("-pc_hypre_ams_relax_type","Relaxation type for AMS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
926:   PetscOptionsInt("-pc_hypre_ams_relax_times","Number of relaxation steps for AMS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
927:   PetscOptionsReal("-pc_hypre_ams_relax_weight","Relaxation weight for AMS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
928:   PetscOptionsReal("-pc_hypre_ams_omega","SSOR coefficient for AMS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
929:   if (flag || flag2 || flag3 || flag4) {
930:     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
931:                                                                       jac->as_relax_times,
932:                                                                       jac->as_relax_weight,
933:                                                                       jac->as_omega));
934:   }
935:   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);
936:   n = 5;
937:   PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options","AMG options for vector Poisson","None",jac->as_amg_alpha_opts,&n,&flag2);
938:   if (flag || flag2) {
939:     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
940:                                                                      jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
941:                                                                      jac->as_amg_alpha_opts[2],       /* AMG relax_type */
942:                                                                      jac->as_amg_alpha_theta,
943:                                                                      jac->as_amg_alpha_opts[3],       /* AMG interp_type */
944:                                                                      jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
945:   }
946:   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);
947:   n = 5;
948:   PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options","AMG options for scalar Poisson solver","None",jac->as_amg_beta_opts,&n,&flag2);
949:   if (flag || flag2) {
950:     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
951:                                                                     jac->as_amg_beta_opts[1],       /* AMG agg_levels */
952:                                                                     jac->as_amg_beta_opts[2],       /* AMG relax_type */
953:                                                                     jac->as_amg_beta_theta,
954:                                                                     jac->as_amg_beta_opts[3],       /* AMG interp_type */
955:                                                                     jac->as_amg_beta_opts[4]));     /* AMG Pmax */
956:   }
957:   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);
958:   if (flag) { /* override HYPRE's default only if the options is used */
959:     PetscStackCallStandard(HYPRE_AMSSetProjectionFrequency,(jac->hsolver,jac->ams_proj_freq));
960:   }
961:   PetscOptionsTail();
962:   return(0);
963: }

967: static PetscErrorCode PCView_HYPRE_AMS(PC pc,PetscViewer viewer)
968: {
969:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
971:   PetscBool      iascii;

974:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
975:   if (iascii) {
976:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS preconditioning\n");
977:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace iterations per application %d\n",jac->as_max_iter);
978:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace cycle type %d\n",jac->ams_cycle_type);
979:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: subspace iteration tolerance %g\n",jac->as_tol);
980:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother type %d\n",jac->as_relax_type);
981:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: number of smoothing steps %d\n",jac->as_relax_times);
982:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother weight %g\n",jac->as_relax_weight);
983:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: smoother omega %g\n",jac->as_omega);
984:     if (jac->alpha_Poisson) {
985:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: vector Poisson solver (passed in by user)\n");
986:     } else {
987:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: vector Poisson solver (computed) \n");
988:     }
989:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);
990:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
991:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);
992:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);
993:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
994:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);
995:     if (!jac->ams_beta_is_zero) {
996:       if (jac->beta_Poisson) {
997:         PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: scalar Poisson solver (passed in by user)\n");
998:       } else {
999:         PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: scalar Poisson solver (computed) \n");
1000:       }
1001:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);
1002:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
1003:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);
1004:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);
1005:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
1006:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);
1007:       if (jac->ams_beta_is_zero_part) {
1008:         PetscViewerASCIIPrintf(viewer,"  HYPRE AMS:     compatible subspace projection frequency %d (-1 HYPRE uses default)\n",jac->ams_proj_freq);
1009:       }
1010:     } else {
1011:       PetscViewerASCIIPrintf(viewer,"  HYPRE AMS: scalar Poisson solver not used (zero-conductivity everywhere) \n");
1012:     }
1013:   }
1014:   return(0);
1015: }

1019: static PetscErrorCode PCSetFromOptions_HYPRE_ADS(PetscOptionItems *PetscOptionsObject,PC pc)
1020: {
1021:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1023:   PetscInt       n;
1024:   PetscBool      flag,flag2,flag3,flag4;

1027:   PetscOptionsHead(PetscOptionsObject,"HYPRE ADS Options");
1028:   PetscOptionsInt("-pc_hypre_ads_print_level","Debugging output level for ADS","None",jac->as_print,&jac->as_print,&flag);
1029:   if (flag) PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
1030:   PetscOptionsInt("-pc_hypre_ads_max_iter","Maximum number of ADS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
1031:   if (flag) PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1032:   PetscOptionsInt("-pc_hypre_ads_cycle_type","Cycle type for ADS multigrid","None",jac->ads_cycle_type,&jac->ads_cycle_type,&flag);
1033:   if (flag) PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ads_cycle_type));
1034:   PetscOptionsReal("-pc_hypre_ads_tol","Error tolerance for ADS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
1035:   if (flag) PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
1036:   PetscOptionsInt("-pc_hypre_ads_relax_type","Relaxation type for ADS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
1037:   PetscOptionsInt("-pc_hypre_ads_relax_times","Number of relaxation steps for ADS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
1038:   PetscOptionsReal("-pc_hypre_ads_relax_weight","Relaxation weight for ADS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
1039:   PetscOptionsReal("-pc_hypre_ads_omega","SSOR coefficient for ADS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
1040:   if (flag || flag2 || flag3 || flag4) {
1041:     PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1042:                                                                       jac->as_relax_times,
1043:                                                                       jac->as_relax_weight,
1044:                                                                       jac->as_omega));
1045:   }
1046:   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);
1047:   n = 5;
1048:   PetscOptionsIntArray("-pc_hypre_ads_ams_options","AMG options for AMS solver inside ADS","None",jac->as_amg_alpha_opts,&n,&flag2);
1049:   PetscOptionsInt("-pc_hypre_ads_ams_cycle_type","Cycle type for AMS solver inside ADS","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag3);
1050:   if (flag || flag2 || flag3) {
1051:     PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type,             /* AMS cycle type */
1052:                                                                 jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1053:                                                                 jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1054:                                                                 jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1055:                                                                 jac->as_amg_alpha_theta,
1056:                                                                 jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1057:                                                                 jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1058:   }
1059:   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);
1060:   n = 5;
1061:   PetscOptionsIntArray("-pc_hypre_ads_amg_options","AMG options for vector AMG solver inside ADS","None",jac->as_amg_beta_opts,&n,&flag2);
1062:   if (flag || flag2) {
1063:     PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1064:                                                                 jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1065:                                                                 jac->as_amg_beta_opts[2],       /* AMG relax_type */
1066:                                                                 jac->as_amg_beta_theta,
1067:                                                                 jac->as_amg_beta_opts[3],       /* AMG interp_type */
1068:                                                                 jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1069:   }
1070:   PetscOptionsTail();
1071:   return(0);
1072: }

1076: static PetscErrorCode PCView_HYPRE_ADS(PC pc,PetscViewer viewer)
1077: {
1078:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1080:   PetscBool      iascii;

1083:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1084:   if (iascii) {
1085:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS preconditioning\n");
1086:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: subspace iterations per application %d\n",jac->as_max_iter);
1087:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: subspace cycle type %d\n",jac->ads_cycle_type);
1088:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: subspace iteration tolerance %g\n",jac->as_tol);
1089:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: smoother type %d\n",jac->as_relax_type);
1090:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: number of smoothing steps %d\n",jac->as_relax_times);
1091:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: smoother weight %g\n",jac->as_relax_weight);
1092:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: smoother omega %g\n",jac->as_omega);
1093:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: AMS solver\n");
1094:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     subspace cycle type %d\n",jac->ams_cycle_type);
1095:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);
1096:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
1097:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);
1098:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);
1099:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
1100:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);
1101:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS: vector Poisson solver\n");
1102:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);
1103:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
1104:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);
1105:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);
1106:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
1107:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS:     boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);
1108:   }
1109:   return(0);
1110: }

1114: static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE(PC pc, Mat G)
1115: {
1116:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1117:   PetscBool      ishypre;

1121:   PetscObjectTypeCompare((PetscObject)G,MATHYPRE,&ishypre);
1122:   if (ishypre) {
1123:     PetscObjectReference((PetscObject)G);
1124:     MatDestroy(&jac->G);
1125:     jac->G = G;
1126:   } else {
1127:     MatReuse reuse = MAT_INITIAL_MATRIX;
1128:     if (jac->G) reuse = MAT_REUSE_MATRIX;
1129:     MatConvert(G,MATHYPRE,reuse,&jac->G);
1130:   }
1131:   return(0);
1132: }

1136: /*@
1137:  PCHYPRESetDiscreteGradient - Set discrete gradient matrix

1139:    Collective on PC

1141:    Input Parameters:
1142: +  pc - the preconditioning context
1143: -  G - the discrete gradient

1145:    Level: intermediate

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

1150: .seealso:
1151: @*/
1152: PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G)
1153: {

1160:   PetscTryMethod(pc,"PCHYPRESetDiscreteGradient_C",(PC,Mat),(pc,G));
1161:   return(0);
1162: }

1166: static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C)
1167: {
1168:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1169:   PetscBool      ishypre;

1173:   PetscObjectTypeCompare((PetscObject)C,MATHYPRE,&ishypre);
1174:   if (ishypre) {
1175:     PetscObjectReference((PetscObject)C);
1176:     MatDestroy(&jac->C);
1177:     jac->C = C;
1178:   } else {
1179:     MatReuse reuse = MAT_INITIAL_MATRIX;
1180:     if (jac->C) reuse = MAT_REUSE_MATRIX;
1181:     MatConvert(C,MATHYPRE,reuse,&jac->C);
1182:   }
1183:   return(0);
1184: }

1188: /*@
1189:  PCHYPRESetDiscreteCurl - Set discrete curl matrix

1191:    Collective on PC

1193:    Input Parameters:
1194: +  pc - the preconditioning context
1195: -  C - the discrete curl

1197:    Level: intermediate

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

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

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

1218: static PetscErrorCode PCHYPRESetPoissonMatrix_HYPRE(PC pc, Mat A, PetscBool isalpha)
1219: {
1220:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1221:   PetscBool      ishypre;

1225:   PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
1226:   if (ishypre) {
1227:     if (isalpha) {
1228:       PetscObjectReference((PetscObject)A);
1229:       MatDestroy(&jac->alpha_Poisson);
1230:       jac->alpha_Poisson = A;
1231:     } else {
1232:       if (A) {
1233:         PetscObjectReference((PetscObject)A);
1234:       } else {
1235:         jac->ams_beta_is_zero = PETSC_TRUE;
1236:       }
1237:       MatDestroy(&jac->beta_Poisson);
1238:       jac->beta_Poisson = A;
1239:     }
1240:   } else {
1241:     if (isalpha) {
1242:       MatReuse reuse = MAT_INITIAL_MATRIX;
1243:       if (jac->alpha_Poisson) reuse = MAT_REUSE_MATRIX;
1244:       MatConvert(A,MATHYPRE,reuse,&jac->alpha_Poisson);
1245:     } else {
1246:       if (A) {
1247:         MatReuse reuse = MAT_INITIAL_MATRIX;
1248:         if (jac->beta_Poisson) reuse = MAT_REUSE_MATRIX;
1249:         MatConvert(A,MATHYPRE,reuse,&jac->beta_Poisson);
1250:       } else {
1251:         MatDestroy(&jac->beta_Poisson);
1252:         jac->ams_beta_is_zero = PETSC_TRUE;
1253:       }
1254:     }
1255:   }
1256:   return(0);
1257: }

1261: /*@
1262:  PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix

1264:    Collective on PC

1266:    Input Parameters:
1267: +  pc - the preconditioning context
1268: -  A - the matrix

1270:    Level: intermediate

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

1274: .seealso:
1275: @*/
1276: PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A)
1277: {

1284:   PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_TRUE));
1285:   return(0);
1286: }

1290: /*@
1291:  PCHYPRESetBetaPoissonMatrix - Set Poisson matrix

1293:    Collective on PC

1295:    Input Parameters:
1296: +  pc - the preconditioning context
1297: -  A - the matrix

1299:    Level: intermediate

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

1304: .seealso:
1305: @*/
1306: PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A)
1307: {

1312:   if (A) {
1315:   }
1316:   PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_FALSE));
1317:   return(0);
1318: }

1322: static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE(PC pc,Vec ozz, Vec zoz, Vec zzo)
1323: {
1324:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1325:   PetscErrorCode     ierr;

1328:   /* throw away any vector if already set */
1329:   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0]));
1330:   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1]));
1331:   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2]));
1332:   jac->constants[0] = NULL;
1333:   jac->constants[1] = NULL;
1334:   jac->constants[2] = NULL;
1335:   VecHYPRE_IJVectorCreate(ozz,&jac->constants[0]);
1336:   VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);
1337:   VecHYPRE_IJVectorCreate(zoz,&jac->constants[1]);
1338:   VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);
1339:   jac->dim = 2;
1340:   if (zzo) {
1341:     VecHYPRE_IJVectorCreate(zzo,&jac->constants[2]);
1342:     VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);
1343:     jac->dim++;
1344:   }
1345:   return(0);
1346: }

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

1353:    Collective on PC

1355:    Input Parameters:
1356: +  pc - the preconditioning context
1357: -  ozz - vector representing (1,0,0) (or (1,0) in 2D)
1358: -  zoz - vector representing (0,1,0) (or (0,1) in 2D)
1359: -  zzo - vector representing (0,0,1) (use NULL in 2D)

1361:    Level: intermediate

1363:    Notes:

1365: .seealso:
1366: @*/
1367: PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo)
1368: {

1379:   PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));
1380:   return(0);
1381: }

1385: static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1386: {
1387:   PC_HYPRE        *jac = (PC_HYPRE*)pc->data;
1388:   Vec             tv;
1389:   PetscInt        i;
1390:   PetscErrorCode  ierr;

1393:   /* throw away any coordinate vector if already set */
1394:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0]));
1395:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1]));
1396:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2]));
1397:   jac->dim = dim;

1399:   /* compute IJ vector for coordinates */
1400:   VecCreate(PetscObjectComm((PetscObject)pc),&tv);
1401:   VecSetType(tv,VECSTANDARD);
1402:   VecSetSizes(tv,nloc,PETSC_DECIDE);
1403:   for (i=0;i<dim;i++) {
1404:     PetscScalar *array;
1405:     PetscInt    j;

1407:     VecHYPRE_IJVectorCreate(tv,&jac->coords[i]);
1408:     VecGetArray(tv,&array);
1409:     for (j=0;j<nloc;j++) {
1410:       array[j] = coords[j*dim+i];
1411:     }
1412:     PetscStackCallStandard(HYPRE_IJVectorSetValues,(jac->coords[i],nloc,NULL,array));
1413:     PetscStackCallStandard(HYPRE_IJVectorAssemble,(jac->coords[i]));
1414:     VecRestoreArray(tv,&array);
1415:   }
1416:   VecDestroy(&tv);
1417:   return(0);
1418: }

1420: /* ---------------------------------------------------------------------------------*/

1424: static PetscErrorCode  PCHYPREGetType_HYPRE(PC pc,const char *name[])
1425: {
1426:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

1429:   *name = jac->hypre_type;
1430:   return(0);
1431: }

1435: static PetscErrorCode  PCHYPRESetType_HYPRE(PC pc,const char name[])
1436: {
1437:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1439:   PetscBool      flag;

1442:   if (jac->hypre_type) {
1443:     PetscStrcmp(jac->hypre_type,name,&flag);
1444:     if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set");
1445:     return(0);
1446:   } else {
1447:     PetscStrallocpy(name, &jac->hypre_type);
1448:   }

1450:   jac->maxiter         = PETSC_DEFAULT;
1451:   jac->tol             = PETSC_DEFAULT;
1452:   jac->printstatistics = PetscLogPrintInfo;

1454:   PetscStrcmp("pilut",jac->hypre_type,&flag);
1455:   if (flag) {
1456:     PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver));
1457:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
1458:     pc->ops->view           = PCView_HYPRE_Pilut;
1459:     jac->destroy            = HYPRE_ParCSRPilutDestroy;
1460:     jac->setup              = HYPRE_ParCSRPilutSetup;
1461:     jac->solve              = HYPRE_ParCSRPilutSolve;
1462:     jac->factorrowsize      = PETSC_DEFAULT;
1463:     return(0);
1464:   }
1465:   PetscStrcmp("parasails",jac->hypre_type,&flag);
1466:   if (flag) {
1467:     PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver));
1468:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
1469:     pc->ops->view           = PCView_HYPRE_ParaSails;
1470:     jac->destroy            = HYPRE_ParaSailsDestroy;
1471:     jac->setup              = HYPRE_ParaSailsSetup;
1472:     jac->solve              = HYPRE_ParaSailsSolve;
1473:     /* initialize */
1474:     jac->nlevels    = 1;
1475:     jac->threshhold = .1;
1476:     jac->filter     = .1;
1477:     jac->loadbal    = 0;
1478:     if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE;
1479:     else jac->logging = (int) PETSC_FALSE;

1481:     jac->ruse = (int) PETSC_FALSE;
1482:     jac->symt = 0;
1483:     PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
1484:     PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
1485:     PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
1486:     PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
1487:     PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
1488:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
1489:     return(0);
1490:   }
1491:   PetscStrcmp("boomeramg",jac->hypre_type,&flag);
1492:   if (flag) {
1493:     HYPRE_BoomerAMGCreate(&jac->hsolver);
1494:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_BoomerAMG;
1495:     pc->ops->view            = PCView_HYPRE_BoomerAMG;
1496:     pc->ops->applytranspose  = PCApplyTranspose_HYPRE_BoomerAMG;
1497:     pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
1498:     jac->destroy             = HYPRE_BoomerAMGDestroy;
1499:     jac->setup               = HYPRE_BoomerAMGSetup;
1500:     jac->solve               = HYPRE_BoomerAMGSolve;
1501:     jac->applyrichardson     = PETSC_FALSE;
1502:     /* these defaults match the hypre defaults */
1503:     jac->cycletype        = 1;
1504:     jac->maxlevels        = 25;
1505:     jac->maxiter          = 1;
1506:     jac->tol              = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
1507:     jac->truncfactor      = 0.0;
1508:     jac->strongthreshold  = .25;
1509:     jac->maxrowsum        = .9;
1510:     jac->coarsentype      = 6;
1511:     jac->measuretype      = 0;
1512:     jac->gridsweeps[0]    = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
1513:     jac->smoothtype       = -1; /* Not set by default */
1514:     jac->smoothnumlevels  = 25;
1515:     jac->eu_level         = 0;
1516:     jac->eu_droptolerance = 0;
1517:     jac->eu_bj            = 0;
1518:     jac->relaxtype[0]     = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */
1519:     jac->relaxtype[2]     = 9; /*G.E. */
1520:     jac->relaxweight      = 1.0;
1521:     jac->outerrelaxweight = 1.0;
1522:     jac->relaxorder       = 1;
1523:     jac->interptype       = 0;
1524:     jac->agg_nl           = 0;
1525:     jac->pmax             = 0;
1526:     jac->truncfactor      = 0.0;
1527:     jac->agg_num_paths    = 1;

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

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

1676: /*
1677:     It only gets here if the HYPRE type has not been set before the call to
1678:    ...SetFromOptions() which actually is most of the time
1679: */
1682: static PetscErrorCode PCSetFromOptions_HYPRE(PetscOptionItems *PetscOptionsObject,PC pc)
1683: {
1685:   PetscInt       indx;
1686:   const char     *type[] = {"pilut","parasails","boomeramg","ams","ads"};
1687:   PetscBool      flg;

1690:   PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");
1691:   PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"boomeramg",&indx,&flg);
1692:   if (flg) {
1693:     PCHYPRESetType_HYPRE(pc,type[indx]);
1694:   } else {
1695:     PCHYPRESetType_HYPRE(pc,"boomeramg");
1696:   }
1697:   if (pc->ops->setfromoptions) {
1698:     pc->ops->setfromoptions(PetscOptionsObject,pc);
1699:   }
1700:   PetscOptionsTail();
1701:   return(0);
1702: }

1706: /*@C
1707:      PCHYPRESetType - Sets which hypre preconditioner you wish to use

1709:    Input Parameters:
1710: +     pc - the preconditioner context
1711: -     name - either  pilut, parasails, boomeramg, ams, ads

1713:    Options Database Keys:
1714:    -pc_hypre_type - One of pilut, parasails, boomeramg, ams, ads

1716:    Level: intermediate

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

1721: @*/
1722: PetscErrorCode  PCHYPRESetType(PC pc,const char name[])
1723: {

1729:   PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));
1730:   return(0);
1731: }

1735: /*@C
1736:      PCHYPREGetType - Gets which hypre preconditioner you are using

1738:    Input Parameter:
1739: .     pc - the preconditioner context

1741:    Output Parameter:
1742: .     name - either  pilut, parasails, boomeramg, ams, ads

1744:    Level: intermediate

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

1749: @*/
1750: PetscErrorCode  PCHYPREGetType(PC pc,const char *name[])
1751: {

1757:   PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));
1758:   return(0);
1759: }

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

1764:    Options Database Keys:
1765: +   -pc_hypre_type - One of pilut, parasails, boomeramg, ams, ads
1766: -   Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX
1767:           preconditioner

1769:    Level: intermediate

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

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

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

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

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

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

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

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

1805: M*/

1809: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
1810: {
1811:   PC_HYPRE       *jac;

1815:   PetscNewLog(pc,&jac);

1817:   pc->data                = jac;
1818:   pc->ops->reset          = PCReset_HYPRE;
1819:   pc->ops->destroy        = PCDestroy_HYPRE;
1820:   pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
1821:   pc->ops->setup          = PCSetUp_HYPRE;
1822:   pc->ops->apply          = PCApply_HYPRE;
1823:   jac->comm_hypre         = MPI_COMM_NULL;
1824:   jac->coords[0]          = NULL;
1825:   jac->coords[1]          = NULL;
1826:   jac->coords[2]          = NULL;
1827:   jac->constants[0]       = NULL;
1828:   jac->constants[1]       = NULL;
1829:   jac->constants[2]       = NULL;
1830:   jac->hmnull             = NULL;
1831:   jac->n_hmnull           = 0;
1832:   /* duplicate communicator for hypre */
1833:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1834:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);
1835:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);
1836:   PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);
1837:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);
1838:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",PCHYPRESetDiscreteCurl_HYPRE);
1839:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE);
1840:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",PCHYPRESetPoissonMatrix_HYPRE);
1841:   return(0);
1842: }

1844: /* ---------------------------------------------------------------------------------------------------------------------------------*/

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

1850:   /* keep copy of PFMG options used so may view them */
1851:   PetscInt its;
1852:   double   tol;
1853:   PetscInt relax_type;
1854:   PetscInt rap_type;
1855:   PetscInt num_pre_relax,num_post_relax;
1856:   PetscInt max_levels;
1857: } PC_PFMG;

1861: PetscErrorCode PCDestroy_PFMG(PC pc)
1862: {
1864:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

1867:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1868:   MPI_Comm_free(&ex->hcomm);
1869:   PetscFree(pc->data);
1870:   return(0);
1871: }

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

1878: PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
1879: {
1881:   PetscBool      iascii;
1882:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

1885:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1886:   if (iascii) {
1887:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG preconditioning\n");
1888:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: max iterations %d\n",ex->its);
1889:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: tolerance %g\n",ex->tol);
1890:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
1891:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: RAP type %s\n",PFMGRAPType[ex->rap_type]);
1892:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
1893:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: max levels %d\n",ex->max_levels);
1894:   }
1895:   return(0);
1896: }


1901: PetscErrorCode PCSetFromOptions_PFMG(PetscOptionItems *PetscOptionsObject,PC pc)
1902: {
1904:   PC_PFMG        *ex = (PC_PFMG*) pc->data;
1905:   PetscBool      flg = PETSC_FALSE;

1908:   PetscOptionsHead(PetscOptionsObject,"PFMG options");
1909:   PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);
1910:   if (flg) {
1911:     int level=3;
1912:     PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level));
1913:   }
1914:   PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);
1915:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its));
1916:   PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
1917:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax));
1918:   PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
1919:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax));

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

1924:   PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);
1925:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol));
1926:   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);
1927:   PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type));
1928:   PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);
1929:   PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type));
1930:   PetscOptionsTail();
1931:   return(0);
1932: }

1936: PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
1937: {
1938:   PetscErrorCode    ierr;
1939:   PC_PFMG           *ex = (PC_PFMG*) pc->data;
1940:   PetscScalar       *yy;
1941:   const PetscScalar *xx;
1942:   PetscInt          ilower[3],iupper[3];
1943:   Mat_HYPREStruct   *mx = (Mat_HYPREStruct*)(pc->pmat->data);

1946:   PetscCitationsRegister(hypreCitation,&cite);
1947:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
1948:   iupper[0] += ilower[0] - 1;
1949:   iupper[1] += ilower[1] - 1;
1950:   iupper[2] += ilower[2] - 1;

1952:   /* copy x values over to hypre */
1953:   PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
1954:   VecGetArrayRead(x,&xx);
1955:   PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx));
1956:   VecRestoreArrayRead(x,&xx);
1957:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
1958:   PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx));

1960:   /* copy solution values back to PETSc */
1961:   VecGetArray(y,&yy);
1962:   PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy));
1963:   VecRestoreArray(y,&yy);
1964:   return(0);
1965: }

1969: 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)
1970: {
1971:   PC_PFMG        *jac = (PC_PFMG*)pc->data;
1973:   PetscInt       oits;

1976:   PetscCitationsRegister(hypreCitation,&cite);
1977:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its));
1978:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol));

1980:   PCApply_PFMG(pc,b,y);
1981:   PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
1982:   *outits = oits;
1983:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1984:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
1985:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol));
1986:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its));
1987:   return(0);
1988: }


1993: PetscErrorCode PCSetUp_PFMG(PC pc)
1994: {
1995:   PetscErrorCode  ierr;
1996:   PC_PFMG         *ex = (PC_PFMG*) pc->data;
1997:   Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
1998:   PetscBool       flg;

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

2004:   /* create the hypre solver object and set its information */
2005:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
2006:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
2007:   PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
2008:   PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver));
2009:   return(0);
2010: }


2013: /*MC
2014:      PCPFMG - the hypre PFMG multigrid solver

2016:    Level: advanced

2018:    Options Database:
2019: + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner
2020: . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2021: . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2022: . -pc_pfmg_tol <tol> tolerance of PFMG
2023: . -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
2024: - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin

2026:    Notes:  This is for CELL-centered descretizations

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

2031: .seealso:  PCMG, MATHYPRESTRUCT
2032: M*/

2036: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
2037: {
2039:   PC_PFMG        *ex;

2042:   PetscNew(&ex); \
2043:   pc->data = ex;

2045:   ex->its            = 1;
2046:   ex->tol            = 1.e-8;
2047:   ex->relax_type     = 1;
2048:   ex->rap_type       = 0;
2049:   ex->num_pre_relax  = 1;
2050:   ex->num_post_relax = 1;
2051:   ex->max_levels     = 0;

2053:   pc->ops->setfromoptions  = PCSetFromOptions_PFMG;
2054:   pc->ops->view            = PCView_PFMG;
2055:   pc->ops->destroy         = PCDestroy_PFMG;
2056:   pc->ops->apply           = PCApply_PFMG;
2057:   pc->ops->applyrichardson = PCApplyRichardson_PFMG;
2058:   pc->ops->setup           = PCSetUp_PFMG;

2060:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
2061:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
2062:   return(0);
2063: }

2065: /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/

2067: /* we know we are working with a HYPRE_SStructMatrix */
2068: typedef struct {
2069:   MPI_Comm            hcomm;       /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
2070:   HYPRE_SStructSolver ss_solver;

2072:   /* keep copy of SYSPFMG options used so may view them */
2073:   PetscInt its;
2074:   double   tol;
2075:   PetscInt relax_type;
2076:   PetscInt num_pre_relax,num_post_relax;
2077: } PC_SysPFMG;

2081: PetscErrorCode PCDestroy_SysPFMG(PC pc)
2082: {
2084:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

2087:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2088:   MPI_Comm_free(&ex->hcomm);
2089:   PetscFree(pc->data);
2090:   return(0);
2091: }

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

2097: PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
2098: {
2100:   PetscBool      iascii;
2101:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

2104:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2105:   if (iascii) {
2106:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG preconditioning\n");
2107:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: max iterations %d\n",ex->its);
2108:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: tolerance %g\n",ex->tol);
2109:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
2110:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
2111:   }
2112:   return(0);
2113: }


2118: PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptionItems *PetscOptionsObject,PC pc)
2119: {
2121:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;
2122:   PetscBool      flg = PETSC_FALSE;

2125:   PetscOptionsHead(PetscOptionsObject,"SysPFMG options");
2126:   PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);
2127:   if (flg) {
2128:     int level=3;
2129:     PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level));
2130:   }
2131:   PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);
2132:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its));
2133:   PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
2134:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax));
2135:   PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
2136:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax));

2138:   PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);
2139:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol));
2140:   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);
2141:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type));
2142:   PetscOptionsTail();
2143:   return(0);
2144: }

2148: PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
2149: {
2150:   PetscErrorCode    ierr;
2151:   PC_SysPFMG        *ex = (PC_SysPFMG*) pc->data;
2152:   PetscScalar       *yy;
2153:   const PetscScalar *xx;
2154:   PetscInt          ilower[3],iupper[3];
2155:   Mat_HYPRESStruct  *mx     = (Mat_HYPRESStruct*)(pc->pmat->data);
2156:   PetscInt          ordering= mx->dofs_order;
2157:   PetscInt          nvars   = mx->nvars;
2158:   PetscInt          part    = 0;
2159:   PetscInt          size;
2160:   PetscInt          i;

2163:   PetscCitationsRegister(hypreCitation,&cite);
2164:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
2165:   iupper[0] += ilower[0] - 1;
2166:   iupper[1] += ilower[1] - 1;
2167:   iupper[2] += ilower[2] - 1;

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

2172:   /* copy x values over to hypre for variable ordering */
2173:   if (ordering) {
2174:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2175:     VecGetArrayRead(x,&xx);
2176:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)xx+(size*i)));
2177:     VecRestoreArrayRead(x,&xx);
2178:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2179:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
2180:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

2182:     /* copy solution values back to PETSc */
2183:     VecGetArray(y,&yy);
2184:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i)));
2185:     VecRestoreArray(y,&yy);
2186:   } else {      /* nodal ordering must be mapped to variable ordering for sys_pfmg */
2187:     PetscScalar *z;
2188:     PetscInt    j, k;

2190:     PetscMalloc1(nvars*size,&z);
2191:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2192:     VecGetArrayRead(x,&xx);

2194:     /* transform nodal to hypre's variable ordering for sys_pfmg */
2195:     for (i= 0; i< size; i++) {
2196:       k= i*nvars;
2197:       for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
2198:     }
2199:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
2200:     VecRestoreArrayRead(x,&xx);
2201:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2202:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

2204:     /* copy solution values back to PETSc */
2205:     VecGetArray(y,&yy);
2206:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
2207:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
2208:     for (i= 0; i< size; i++) {
2209:       k= i*nvars;
2210:       for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
2211:     }
2212:     VecRestoreArray(y,&yy);
2213:     PetscFree(z);
2214:   }
2215:   return(0);
2216: }

2220: 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)
2221: {
2222:   PC_SysPFMG     *jac = (PC_SysPFMG*)pc->data;
2224:   PetscInt       oits;

2227:   PetscCitationsRegister(hypreCitation,&cite);
2228:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its));
2229:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol));
2230:   PCApply_SysPFMG(pc,b,y);
2231:   PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,(HYPRE_Int *)&oits));
2232:   *outits = oits;
2233:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2234:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
2235:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol));
2236:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its));
2237:   return(0);
2238: }


2243: PetscErrorCode PCSetUp_SysPFMG(PC pc)
2244: {
2245:   PetscErrorCode   ierr;
2246:   PC_SysPFMG       *ex = (PC_SysPFMG*) pc->data;
2247:   Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
2248:   PetscBool        flg;

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

2254:   /* create the hypre sstruct solver object and set its information */
2255:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2256:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2257:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver));
2258:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2259:   return(0);
2260: }


2263: /*MC
2264:      PCSysPFMG - the hypre SysPFMG multigrid solver

2266:    Level: advanced

2268:    Options Database:
2269: + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner
2270: . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2271: . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2272: . -pc_syspfmg_tol <tol> tolerance of SysPFMG
2273: . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel

2275:    Notes:  This is for CELL-centered descretizations

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

2281: .seealso:  PCMG, MATHYPRESSTRUCT
2282: M*/

2286: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
2287: {
2289:   PC_SysPFMG     *ex;

2292:   PetscNew(&ex); \
2293:   pc->data = ex;

2295:   ex->its            = 1;
2296:   ex->tol            = 1.e-8;
2297:   ex->relax_type     = 1;
2298:   ex->num_pre_relax  = 1;
2299:   ex->num_post_relax = 1;

2301:   pc->ops->setfromoptions  = PCSetFromOptions_SysPFMG;
2302:   pc->ops->view            = PCView_SysPFMG;
2303:   pc->ops->destroy         = PCDestroy_SysPFMG;
2304:   pc->ops->apply           = PCApply_SysPFMG;
2305:   pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
2306:   pc->ops->setup           = PCSetUp_SysPFMG;

2308:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
2309:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2310:   return(0);
2311: }