Actual source code: hypre.c

petsc-3.11.2 2019-05-18
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:   PetscBool nodal_relax;
 77:   PetscInt  nodal_relax_levels;

 79:   PetscInt  nodal_coarsening;
 80:   PetscInt  nodal_coarsening_diag;
 81:   PetscInt  vec_interp_variant;
 82:   PetscInt  vec_interp_qmax;
 83:   PetscBool vec_interp_smooth;
 84:   PetscInt  interp_refine;

 86:   HYPRE_IJVector  *hmnull;
 87:   HYPRE_ParVector *phmnull;  /* near null space passed to hypre */
 88:   PetscInt        n_hmnull;
 89:   Vec             hmnull_constant;
 90:   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 */

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

107:   /* additional data */
108:   Mat G;             /* MatHYPRE */
109:   Mat C;             /* MatHYPRE */
110:   Mat alpha_Poisson; /* MatHYPRE */
111:   Mat beta_Poisson;  /* MatHYPRE */

113:   /* extra information for AMS */
114:   PetscInt       dim; /* geometrical dimension */
115:   HYPRE_IJVector coords[3];
116:   HYPRE_IJVector constants[3];
117:   Mat            RT_PiFull, RT_Pi[3];
118:   Mat            ND_PiFull, ND_Pi[3];
119:   PetscBool      ams_beta_is_zero;
120:   PetscBool      ams_beta_is_zero_part;
121:   PetscInt       ams_proj_freq;
122: } PC_HYPRE;

124: PetscErrorCode PCHYPREGetSolver(PC pc,HYPRE_Solver *hsolver)
125: {
126:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

129:   *hsolver = jac->hsolver;
130:   return(0);
131: }

133: /* Resets (frees) Hypre's representation of the near null space */
134: static PetscErrorCode PCHYPREResetNearNullSpace_Private(PC pc)
135: {
136:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
137:   PetscInt       i;
139:   PETSC_UNUSED PetscScalar *petscvecarray;

142:   for (i=0; i<jac->n_hmnull; i++) {
143:     VecHYPRE_ParVectorReplacePointer(jac->hmnull[i],jac->hmnull_hypre_data_array[i],petscvecarray);
144:     PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->hmnull[i]));
145:   }
146:   PetscFree(jac->hmnull);
147:   PetscFree(jac->hmnull_hypre_data_array);
148:   PetscFree(jac->phmnull);
149:   VecDestroy(&jac->hmnull_constant);
150:   jac->n_hmnull = 0;
151:   return(0);
152: }

154: static PetscErrorCode PCSetUp_HYPRE(PC pc)
155: {
156:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
157:   Mat_HYPRE          *hjac;
158:   HYPRE_ParCSRMatrix hmat;
159:   HYPRE_ParVector    bv,xv;
160:   PetscBool          ishypre;
161:   PetscErrorCode     ierr;

164:   if (!jac->hypre_type) {
165:     PCHYPRESetType(pc,"boomeramg");
166:   }

168:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRE,&ishypre);
169:   if (!ishypre) {
170:     MatDestroy(&jac->hpmat);
171:     MatConvert(pc->pmat,MATHYPRE,MAT_INITIAL_MATRIX,&jac->hpmat);
172:   } else {
173:     PetscObjectReference((PetscObject)pc->pmat);
174:     MatDestroy(&jac->hpmat);
175:     jac->hpmat = pc->pmat;
176:   }
177:   hjac = (Mat_HYPRE*)(jac->hpmat->data);

179:   /* special case for BoomerAMG */
180:   if (jac->setup == HYPRE_BoomerAMGSetup) {
181:     MatNullSpace    mnull;
182:     PetscBool       has_const;
183:     PetscInt        bs,nvec,i;
184:     const Vec       *vecs;
185:     PetscScalar     *petscvecarray;

187:     MatGetBlockSize(pc->pmat,&bs);
188:     if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
189:     MatGetNearNullSpace(pc->mat, &mnull);
190:     if (mnull) {
191:       PCHYPREResetNearNullSpace_Private(pc);
192:       MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);
193:       PetscMalloc1(nvec+1,&jac->hmnull);
194:       PetscMalloc1(nvec+1,&jac->hmnull_hypre_data_array);
195:       PetscMalloc1(nvec+1,&jac->phmnull);
196:       for (i=0; i<nvec; i++) {
197:         VecHYPRE_IJVectorCreate(vecs[i],&jac->hmnull[i]);
198:         VecGetArrayRead(vecs[i],(const PetscScalar **)&petscvecarray);
199:         VecHYPRE_ParVectorReplacePointer(jac->hmnull[i],petscvecarray,jac->hmnull_hypre_data_array[i]);
200:         VecRestoreArrayRead(vecs[i],(const PetscScalar **)&petscvecarray);
201:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[i],(void**)&jac->phmnull[i]));
202:       }
203:       if (has_const) {
204:         MatCreateVecs(pc->pmat,&jac->hmnull_constant,NULL);
205:         VecSet(jac->hmnull_constant,1);
206:         VecNormalize(jac->hmnull_constant,NULL);
207:         VecHYPRE_IJVectorCreate(jac->hmnull_constant,&jac->hmnull[nvec]);
208:         VecGetArrayRead(jac->hmnull_constant,(const PetscScalar **)&petscvecarray);
209:         VecHYPRE_ParVectorReplacePointer(jac->hmnull[nvec],petscvecarray,jac->hmnull_hypre_data_array[nvec]);
210:         VecRestoreArrayRead(jac->hmnull_constant,(const PetscScalar **)&petscvecarray);
211:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[nvec],(void**)&jac->phmnull[nvec]));
212:         nvec++;
213:       }
214:       PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVectors,(jac->hsolver,nvec,jac->phmnull));
215:       jac->n_hmnull = nvec;
216:     }
217:   }

219:   /* special case for AMS */
220:   if (jac->setup == HYPRE_AMSSetup) {
221:     Mat_HYPRE          *hm;
222:     HYPRE_ParCSRMatrix parcsr;
223:     if (!jac->coords[0] && !jac->constants[0] && !(jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]))) {
224:       SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs either the coordinate vectors via PCSetCoordinates() or the edge constant vectors via PCHYPRESetEdgeConstantVectors() or the interpolation matrix via PCHYPRESetInterpolations");
225:     }
226:     if (jac->dim) {
227:       PetscStackCallStandard(HYPRE_AMSSetDimension,(jac->hsolver,jac->dim));
228:     }
229:     if (jac->constants[0]) {
230:       HYPRE_ParVector ozz,zoz,zzo = NULL;
231:       PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[0],(void**)(&ozz)));
232:       PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[1],(void**)(&zoz)));
233:       if (jac->constants[2]) {
234:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[2],(void**)(&zzo)));
235:       }
236:       PetscStackCallStandard(HYPRE_AMSSetEdgeConstantVectors,(jac->hsolver,ozz,zoz,zzo));
237:     }
238:     if (jac->coords[0]) {
239:       HYPRE_ParVector coords[3];
240:       coords[0] = NULL;
241:       coords[1] = NULL;
242:       coords[2] = NULL;
243:       if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0],(void**)(&coords[0])));
244:       if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1],(void**)(&coords[1])));
245:       if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2],(void**)(&coords[2])));
246:       PetscStackCallStandard(HYPRE_AMSSetCoordinateVectors,(jac->hsolver,coords[0],coords[1],coords[2]));
247:     }
248:     if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient");
249:     hm = (Mat_HYPRE*)(jac->G->data);
250:     PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
251:     PetscStackCallStandard(HYPRE_AMSSetDiscreteGradient,(jac->hsolver,parcsr));
252:     if (jac->alpha_Poisson) {
253:       hm = (Mat_HYPRE*)(jac->alpha_Poisson->data);
254:       PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
255:       PetscStackCallStandard(HYPRE_AMSSetAlphaPoissonMatrix,(jac->hsolver,parcsr));
256:     }
257:     if (jac->ams_beta_is_zero) {
258:       PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,NULL));
259:     } else if (jac->beta_Poisson) {
260:       hm = (Mat_HYPRE*)(jac->beta_Poisson->data);
261:       PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
262:       PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,parcsr));
263:     }
264:     if (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1])) {
265:       PetscInt           i;
266:       HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3];
267:       if (jac->ND_PiFull) {
268:         hm = (Mat_HYPRE*)(jac->ND_PiFull->data);
269:         PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&nd_parcsrfull)));
270:       } else {
271:         nd_parcsrfull = NULL;
272:       }
273:       for (i=0;i<3;++i) {
274:         if (jac->ND_Pi[i]) {
275:           hm = (Mat_HYPRE*)(jac->ND_Pi[i]->data);
276:           PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&nd_parcsr[i])));
277:         } else {
278:           nd_parcsr[i] = NULL;
279:         }
280:       }
281:       PetscStackCallStandard(HYPRE_AMSSetInterpolations,(jac->hsolver,nd_parcsrfull,nd_parcsr[0],nd_parcsr[1],nd_parcsr[2]));
282:     }
283:   }
284:   /* special case for ADS */
285:   if (jac->setup == HYPRE_ADSSetup) {
286:     Mat_HYPRE          *hm;
287:     HYPRE_ParCSRMatrix parcsr;
288:     if (!jac->coords[0] && !((jac->RT_PiFull || (jac->RT_Pi[0] && jac->RT_Pi[1])) && (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1])))) {
289:       SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs either the coordinate vectors via PCSetCoordinates() or the interpolation matrices via PCHYPRESetInterpolations");
290:     }
291:     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");
292:     if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient");
293:     if (!jac->C) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs the discrete curl operator via PCHYPRESetDiscreteGradient");
294:     if (jac->coords[0]) {
295:       HYPRE_ParVector coords[3];
296:       coords[0] = NULL;
297:       coords[1] = NULL;
298:       coords[2] = NULL;
299:       if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0],(void**)(&coords[0])));
300:       if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1],(void**)(&coords[1])));
301:       if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2],(void**)(&coords[2])));
302:       PetscStackCallStandard(HYPRE_ADSSetCoordinateVectors,(jac->hsolver,coords[0],coords[1],coords[2]));
303:     }
304:     hm = (Mat_HYPRE*)(jac->G->data);
305:     PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
306:     PetscStackCallStandard(HYPRE_ADSSetDiscreteGradient,(jac->hsolver,parcsr));
307:     hm = (Mat_HYPRE*)(jac->C->data);
308:     PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
309:     PetscStackCallStandard(HYPRE_ADSSetDiscreteCurl,(jac->hsolver,parcsr));
310:     if ((jac->RT_PiFull || (jac->RT_Pi[0] && jac->RT_Pi[1])) && (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]))) {
311:       PetscInt           i;
312:       HYPRE_ParCSRMatrix rt_parcsrfull, rt_parcsr[3];
313:       HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3];
314:       if (jac->RT_PiFull) {
315:         hm = (Mat_HYPRE*)(jac->RT_PiFull->data);
316:         PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&rt_parcsrfull)));
317:       } else {
318:         rt_parcsrfull = NULL;
319:       }
320:       for (i=0;i<3;++i) {
321:         if (jac->RT_Pi[i]) {
322:           hm = (Mat_HYPRE*)(jac->RT_Pi[i]->data);
323:           PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&rt_parcsr[i])));
324:         } else {
325:           rt_parcsr[i] = NULL;
326:         }
327:       }
328:       if (jac->ND_PiFull) {
329:         hm = (Mat_HYPRE*)(jac->ND_PiFull->data);
330:         PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&nd_parcsrfull)));
331:       } else {
332:         nd_parcsrfull = NULL;
333:       }
334:       for (i=0;i<3;++i) {
335:         if (jac->ND_Pi[i]) {
336:           hm = (Mat_HYPRE*)(jac->ND_Pi[i]->data);
337:           PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&nd_parcsr[i])));
338:         } else {
339:           nd_parcsr[i] = NULL;
340:         }
341:       }
342:       PetscStackCallStandard(HYPRE_ADSSetInterpolations,(jac->hsolver,rt_parcsrfull,rt_parcsr[0],rt_parcsr[1],rt_parcsr[2],nd_parcsrfull,nd_parcsr[0],nd_parcsr[1],nd_parcsr[2]));
343:     }
344:   }
345:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
346:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b,(void**)&bv));
347:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x,(void**)&xv));
348:   PetscStackCallStandard(jac->setup,(jac->hsolver,hmat,bv,xv));
349:   return(0);
350: }

352: static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
353: {
354:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
355:   Mat_HYPRE          *hjac = (Mat_HYPRE*)(jac->hpmat->data);
356:   PetscErrorCode     ierr;
357:   HYPRE_ParCSRMatrix hmat;
358:   PetscScalar        *xv;
359:   const PetscScalar  *bv,*sbv;
360:   HYPRE_ParVector    jbv,jxv;
361:   PetscScalar        *sxv;
362:   PetscInt           hierr;

365:   PetscCitationsRegister(hypreCitation,&cite);
366:   if (!jac->applyrichardson) {VecSet(x,0.0);}
367:   VecGetArrayRead(b,&bv);
368:   VecGetArray(x,&xv);
369:   VecHYPRE_ParVectorReplacePointer(hjac->b,(PetscScalar*)bv,sbv);
370:   VecHYPRE_ParVectorReplacePointer(hjac->x,xv,sxv);

372:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
373:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b,(void**)&jbv));
374:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x,(void**)&jxv));
375:   PetscStackCall("Hypre solve",h(*jac->solve)(jac->hsolver,hmat,jbv,jxv);
376:   if (hierr && hierr != HYPRE_ERROR_CONV) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
377:   if (hierr) hypre__global_error = 0;);

379:   if (jac->setup == HYPRE_AMSSetup && jac->ams_beta_is_zero_part) {
380:     PetscStackCallStandard(HYPRE_AMSProjectOutGradients,(jac->hsolver,jxv));
381:   }
382:   VecHYPRE_ParVectorReplacePointer(hjac->b,(PetscScalar*)sbv,bv);
383:   VecHYPRE_ParVectorReplacePointer(hjac->x,sxv,xv);
384:   VecRestoreArray(x,&xv);
385:   VecRestoreArrayRead(b,&bv);
386:   return(0);
387: }

389: static PetscErrorCode PCReset_HYPRE(PC pc)
390: {
391:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;

395:   MatDestroy(&jac->hpmat);
396:   MatDestroy(&jac->G);
397:   MatDestroy(&jac->C);
398:   MatDestroy(&jac->alpha_Poisson);
399:   MatDestroy(&jac->beta_Poisson);
400:   MatDestroy(&jac->RT_PiFull);
401:   MatDestroy(&jac->RT_Pi[0]);
402:   MatDestroy(&jac->RT_Pi[1]);
403:   MatDestroy(&jac->RT_Pi[2]);
404:   MatDestroy(&jac->ND_PiFull);
405:   MatDestroy(&jac->ND_Pi[0]);
406:   MatDestroy(&jac->ND_Pi[1]);
407:   MatDestroy(&jac->ND_Pi[2]);
408:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0])); jac->coords[0] = NULL;
409:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1])); jac->coords[1] = NULL;
410:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2])); jac->coords[2] = NULL;
411:   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0])); jac->constants[0] = NULL;
412:   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1])); jac->constants[1] = NULL;
413:   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2])); jac->constants[2] = NULL;
414:   PCHYPREResetNearNullSpace_Private(pc);
415:   jac->ams_beta_is_zero = PETSC_FALSE;
416:   jac->dim = 0;
417:   return(0);
418: }

420: static PetscErrorCode PCDestroy_HYPRE(PC pc)
421: {
422:   PC_HYPRE                 *jac = (PC_HYPRE*)pc->data;
423:   PetscErrorCode           ierr;

426:   PCReset_HYPRE(pc);
427:   if (jac->destroy) PetscStackCallStandard(jac->destroy,(jac->hsolver));
428:   PetscFree(jac->hypre_type);
429:   if (jac->comm_hypre != MPI_COMM_NULL) { MPI_Comm_free(&(jac->comm_hypre));}
430:   PetscFree(pc->data);

432:   PetscObjectChangeTypeName((PetscObject)pc,0);
433:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);
434:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);
435:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetCoordinates_C",NULL);
436:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",NULL);
437:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",NULL);
438:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",NULL);
439:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetConstantEdgeVectors_C",NULL);
440:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",NULL);
441:   return(0);
442: }

444: /* --------------------------------------------------------------------------------------------*/
445: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PetscOptionItems *PetscOptionsObject,PC pc)
446: {
447:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
449:   PetscBool      flag;

452:   PetscOptionsHead(PetscOptionsObject,"HYPRE Pilut Options");
453:   PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
454:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter));
455:   PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
456:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol));
457:   PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
458:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize));
459:   PetscOptionsTail();
460:   return(0);
461: }

463: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
464: {
465:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
467:   PetscBool      iascii;

470:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
471:   if (iascii) {
472:     PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut preconditioning\n");
473:     if (jac->maxiter != PETSC_DEFAULT) {
474:       PetscViewerASCIIPrintf(viewer,"    maximum number of iterations %d\n",jac->maxiter);
475:     } else {
476:       PetscViewerASCIIPrintf(viewer,"    default maximum number of iterations \n");
477:     }
478:     if (jac->tol != PETSC_DEFAULT) {
479:       PetscViewerASCIIPrintf(viewer,"    drop tolerance %g\n",(double)jac->tol);
480:     } else {
481:       PetscViewerASCIIPrintf(viewer,"    default drop tolerance \n");
482:     }
483:     if (jac->factorrowsize != PETSC_DEFAULT) {
484:       PetscViewerASCIIPrintf(viewer,"    factor row size %d\n",jac->factorrowsize);
485:     } else {
486:       PetscViewerASCIIPrintf(viewer,"    default factor row size \n");
487:     }
488:   }
489:   return(0);
490: }

492: /* --------------------------------------------------------------------------------------------*/
493: static PetscErrorCode PCSetFromOptions_HYPRE_Euclid(PetscOptionItems *PetscOptionsObject,PC pc)
494: {
495:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
497:   PetscBool      flag;

500:   PetscOptionsHead(PetscOptionsObject,"HYPRE Euclid Options");
501:   PetscOptionsInt("-pc_hypre_euclid_level","Factorization levels","None",jac->eu_level,&jac->eu_level,&flag);
502:   if (flag) PetscStackCallStandard(HYPRE_EuclidSetLevel,(jac->hsolver,jac->eu_level));
503:   PetscOptionsTail();
504:   return(0);
505: }

507: static PetscErrorCode PCView_HYPRE_Euclid(PC pc,PetscViewer viewer)
508: {
509:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
511:   PetscBool      iascii;

514:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
515:   if (iascii) {
516:     PetscViewerASCIIPrintf(viewer,"  HYPRE Euclid preconditioning\n");
517:     if (jac->eu_level != PETSC_DEFAULT) {
518:       PetscViewerASCIIPrintf(viewer,"    factorization levels %d\n",jac->eu_level);
519:     } else {
520:       PetscViewerASCIIPrintf(viewer,"    default factorization levels \n");
521:     }
522:   }
523:   return(0);
524: }

526: /* --------------------------------------------------------------------------------------------*/

528: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
529: {
530:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
531:   Mat_HYPRE          *hjac = (Mat_HYPRE*)(jac->hpmat->data);
532:   PetscErrorCode     ierr;
533:   HYPRE_ParCSRMatrix hmat;
534:   PetscScalar        *xv;
535:   const PetscScalar  *bv;
536:   HYPRE_ParVector    jbv,jxv;
537:   PetscScalar        *sbv,*sxv;
538:   PetscInt           hierr;

541:   PetscCitationsRegister(hypreCitation,&cite);
542:   VecSet(x,0.0);
543:   VecGetArrayRead(b,&bv);
544:   VecGetArray(x,&xv);
545:   VecHYPRE_ParVectorReplacePointer(hjac->b,(PetscScalar*)bv,sbv);
546:   VecHYPRE_ParVectorReplacePointer(hjac->x,xv,sxv);

548:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
549:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b,(void**)&jbv));
550:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x,(void**)&jxv));

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

557:   VecHYPRE_ParVectorReplacePointer(hjac->b,sbv,bv);
558:   VecHYPRE_ParVectorReplacePointer(hjac->x,sxv,xv);
559:   VecRestoreArray(x,&xv);
560:   VecRestoreArrayRead(b,&bv);
561:   return(0);
562: }

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

567: static const char *HYPREBoomerAMGCycleType[]   = {"","V","W"};
568: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
569: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
570: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
571: static const char *HYPREBoomerAMGSmoothType[]   = {"Schwarz-smoothers","Pilut","ParaSails","Euclid"};
572: static const char *HYPREBoomerAMGRelaxType[]   = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi",
573:                                                   "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi",
574:                                                   "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination",
575:                                                   "" /* 10 */, "" /* 11 */, "" /* 12 */, "l1-Gauss-Seidel" /* nonsymmetric */, "backward-l1-Gauss-Seidel" /* nonsymmetric */,
576:                                                   "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"};
577: static const char *HYPREBoomerAMGInterpType[]  = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
578:                                                   "ext+i-cc", "standard", "standard-wts", "block", "block-wtd", "FF", "FF1"};
579: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PetscOptionItems *PetscOptionsObject,PC pc)
580: {
581:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
583:   PetscInt       bs,n,indx,level;
584:   PetscBool      flg, tmp_truth;
585:   double         tmpdbl, twodbl[2];

588:   PetscOptionsHead(PetscOptionsObject,"HYPRE BoomerAMG Options");
589:   PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);
590:   if (flg) {
591:     jac->cycletype = indx+1;
592:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
593:   }
594:   PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
595:   if (flg) {
596:     if (jac->maxlevels < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
597:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
598:   }
599:   PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
600:   if (flg) {
601:     if (jac->maxiter < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
602:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
603:   }
604:   PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);
605:   if (flg) {
606:     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);
607:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
608:   }
609:   bs = 1;
610:   if (pc->pmat) {
611:     MatGetBlockSize(pc->pmat,&bs);
612:   }
613:   PetscOptionsInt("-pc_hypre_boomeramg_numfunctions","Number of functions","HYPRE_BoomerAMGSetNumFunctions",bs,&bs,&flg);
614:   if (flg) {
615:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
616:   }

618:   PetscOptionsScalar("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);
619:   if (flg) {
620:     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);
621:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
622:   }

624:   PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);
625:   if (flg) {
626:     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);
627:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
628:   }

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

634:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
635:   }

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

641:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
642:   }


645:   PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
646:   if (flg) {
647:     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);
648:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
649:   }
650:   PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
651:   if (flg) {
652:     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);
653:     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);
654:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
655:   }

657:   /* Grid sweeps */
658:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);
659:   if (flg) {
660:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx));
661:     /* modify the jac structure so we can view the updated options with PC_View */
662:     jac->gridsweeps[0] = indx;
663:     jac->gridsweeps[1] = indx;
664:     /*defaults coarse to 1 */
665:     jac->gridsweeps[2] = 1;
666:   }
667:   PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen","Use a nodal based coarsening 1-6","HYPRE_BoomerAMGSetNodal",jac->nodal_coarsening,&jac->nodal_coarsening,&flg);
668:   if (flg) {
669:     PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,jac->nodal_coarsening));
670:   }
671:   PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen_diag","Diagonal in strength matrix for nodal based coarsening 0-2","HYPRE_BoomerAMGSetNodalDiag",jac->nodal_coarsening_diag,&jac->nodal_coarsening_diag,&flg);
672:   if (flg) {
673:     PetscStackCallStandard(HYPRE_BoomerAMGSetNodalDiag,(jac->hsolver,jac->nodal_coarsening_diag));
674:   }
675:   PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_variant","Variant of algorithm 1-3","HYPRE_BoomerAMGSetInterpVecVariant",jac->vec_interp_variant, &jac->vec_interp_variant,&flg);
676:   if (flg) {
677:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecVariant,(jac->hsolver,jac->vec_interp_variant));
678:   }
679:   PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_qmax","Max elements per row for each Q","HYPRE_BoomerAMGSetInterpVecQMax",jac->vec_interp_qmax, &jac->vec_interp_qmax,&flg);
680:   if (flg) {
681:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecQMax,(jac->hsolver,jac->vec_interp_qmax));
682:   }
683:   PetscOptionsBool("-pc_hypre_boomeramg_vec_interp_smooth","Whether to smooth the interpolation vectors","HYPRE_BoomerAMGSetSmoothInterpVectors",jac->vec_interp_smooth, &jac->vec_interp_smooth,&flg);
684:   if (flg) {
685:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothInterpVectors,(jac->hsolver,jac->vec_interp_smooth));
686:   }
687:   PetscOptionsInt("-pc_hypre_boomeramg_interp_refine","Preprocess the interpolation matrix through iterative weight refinement","HYPRE_BoomerAMGSetInterpRefine",jac->interp_refine, &jac->interp_refine,&flg);
688:   if (flg) {
689:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpRefine,(jac->hsolver,jac->interp_refine));
690:   }
691:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);
692:   if (flg) {
693:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1));
694:     jac->gridsweeps[0] = indx;
695:   }
696:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);
697:   if (flg) {
698:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2));
699:     jac->gridsweeps[1] = indx;
700:   }
701:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);
702:   if (flg) {
703:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3));
704:     jac->gridsweeps[2] = indx;
705:   }

707:   /* Smooth type */
708:   PetscOptionsEList("-pc_hypre_boomeramg_smooth_type","Enable more complex smoothers","None",HYPREBoomerAMGSmoothType,ALEN(HYPREBoomerAMGSmoothType),HYPREBoomerAMGSmoothType[0],&indx,&flg);
709:   if (flg) {
710:     jac->smoothtype = indx;
711:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,indx+6));
712:     jac->smoothnumlevels = 25;
713:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,25));
714:   }

716:   /* Number of smoothing levels */
717:   PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_levels","Number of levels on which more complex smoothers are used","None",25,&indx,&flg);
718:   if (flg && (jac->smoothtype != -1)) {
719:     jac->smoothnumlevels = indx;
720:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,indx));
721:   }

723:   /* Number of levels for ILU(k) for Euclid */
724:   PetscOptionsInt("-pc_hypre_boomeramg_eu_level","Number of levels for ILU(k) in Euclid smoother","None",0,&indx,&flg);
725:   if (flg && (jac->smoothtype == 3)) {
726:     jac->eu_level = indx;
727:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,indx));
728:   }

730:   /* Filter for ILU(k) for Euclid */
731:   double droptolerance;
732:   PetscOptionsScalar("-pc_hypre_boomeramg_eu_droptolerance","Drop tolerance for ILU(k) in Euclid smoother","None",0,&droptolerance,&flg);
733:   if (flg && (jac->smoothtype == 3)) {
734:     jac->eu_droptolerance = droptolerance;
735:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,droptolerance));
736:   }

738:   /* Use Block Jacobi ILUT for Euclid */
739:   PetscOptionsBool("-pc_hypre_boomeramg_eu_bj", "Use Block Jacobi for ILU in Euclid smoother?", "None", PETSC_FALSE, &tmp_truth, &flg);
740:   if (flg && (jac->smoothtype == 3)) {
741:     jac->eu_bj = tmp_truth;
742:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuBJ,(jac->hsolver,jac->eu_bj));
743:   }

745:   /* Relax type */
746:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
747:   if (flg) {
748:     jac->relaxtype[0] = jac->relaxtype[1]  = indx;
749:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx));
750:     /* by default, coarse type set to 9 */
751:     jac->relaxtype[2] = 9;
752:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, 9, 3));
753:   }
754:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
755:   if (flg) {
756:     jac->relaxtype[0] = indx;
757:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1));
758:   }
759:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
760:   if (flg) {
761:     jac->relaxtype[1] = indx;
762:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2));
763:   }
764:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);
765:   if (flg) {
766:     jac->relaxtype[2] = indx;
767:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3));
768:   }

770:   /* Relaxation Weight */
771:   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);
772:   if (flg) {
773:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl));
774:     jac->relaxweight = tmpdbl;
775:   }

777:   n         = 2;
778:   twodbl[0] = twodbl[1] = 1.0;
779:   PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
780:   if (flg) {
781:     if (n == 2) {
782:       indx =  (int)PetscAbsReal(twodbl[1]);
783:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx));
784:     } 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);
785:   }

787:   /* Outer relaxation Weight */
788:   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);
789:   if (flg) {
790:     PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl));
791:     jac->outerrelaxweight = tmpdbl;
792:   }

794:   n         = 2;
795:   twodbl[0] = twodbl[1] = 1.0;
796:   PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
797:   if (flg) {
798:     if (n == 2) {
799:       indx =  (int)PetscAbsReal(twodbl[1]);
800:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx));
801:     } 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);
802:   }

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

807:   if (flg && tmp_truth) {
808:     jac->relaxorder = 0;
809:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
810:   }
811:   PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);
812:   if (flg) {
813:     jac->measuretype = indx;
814:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
815:   }
816:   /* update list length 3/07 */
817:   PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);
818:   if (flg) {
819:     jac->coarsentype = indx;
820:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
821:   }

823:   /* new 3/07 */
824:   PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);
825:   if (flg) {
826:     jac->interptype = indx;
827:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
828:   }

830:   PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
831:   if (flg) {
832:     level = 3;
833:     PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);

835:     jac->printstatistics = PETSC_TRUE;
836:     PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level));
837:   }

839:   PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);
840:   if (flg) {
841:     level = 3;
842:     PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);

844:     jac->printstatistics = PETSC_TRUE;
845:     PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level));
846:   }

848:   PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);
849:   if (flg && tmp_truth) {
850:     PetscInt tmp_int;
851:     PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);
852:     if (flg) jac->nodal_relax_levels = tmp_int;
853:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6));
854:     PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1));
855:     PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0));
856:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels));
857:   }

859:   PetscOptionsTail();
860:   return(0);
861: }

863: 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)
864: {
865:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
867:   PetscInt       oits;

870:   PetscCitationsRegister(hypreCitation,&cite);
871:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter));
872:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol));
873:   jac->applyrichardson = PETSC_TRUE;
874:   PCApply_HYPRE(pc,b,y);
875:   jac->applyrichardson = PETSC_FALSE;
876:   PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
877:   *outits = oits;
878:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
879:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
880:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
881:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
882:   return(0);
883: }


886: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
887: {
888:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
890:   PetscBool      iascii;

893:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
894:   if (iascii) {
895:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG preconditioning\n");
896:     PetscViewerASCIIPrintf(viewer,"    Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);
897:     PetscViewerASCIIPrintf(viewer,"    Maximum number of levels %D\n",jac->maxlevels);
898:     PetscViewerASCIIPrintf(viewer,"    Maximum number of iterations PER hypre call %D\n",jac->maxiter);
899:     PetscViewerASCIIPrintf(viewer,"    Convergence tolerance PER hypre call %g\n",(double)jac->tol);
900:     PetscViewerASCIIPrintf(viewer,"    Threshold for strong coupling %g\n",(double)jac->strongthreshold);
901:     PetscViewerASCIIPrintf(viewer,"    Interpolation truncation factor %g\n",(double)jac->truncfactor);
902:     PetscViewerASCIIPrintf(viewer,"    Interpolation: max elements per row %D\n",jac->pmax);
903:     if (jac->interp_refine) {
904:       PetscViewerASCIIPrintf(viewer,"    Interpolation: number of steps of weighted refinement %D\n",jac->interp_refine);
905:     }
906:     PetscViewerASCIIPrintf(viewer,"    Number of levels of aggressive coarsening %D\n",jac->agg_nl);
907:     PetscViewerASCIIPrintf(viewer,"    Number of paths for aggressive coarsening %D\n",jac->agg_num_paths);

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

911:     PetscViewerASCIIPrintf(viewer,"    Sweeps down         %D\n",jac->gridsweeps[0]);
912:     PetscViewerASCIIPrintf(viewer,"    Sweeps up           %D\n",jac->gridsweeps[1]);
913:     PetscViewerASCIIPrintf(viewer,"    Sweeps on coarse    %D\n",jac->gridsweeps[2]);

915:     PetscViewerASCIIPrintf(viewer,"    Relax down          %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[0]]);
916:     PetscViewerASCIIPrintf(viewer,"    Relax up            %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[1]]);
917:     PetscViewerASCIIPrintf(viewer,"    Relax on coarse     %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[2]]);

919:     PetscViewerASCIIPrintf(viewer,"    Relax weight  (all)      %g\n",(double)jac->relaxweight);
920:     PetscViewerASCIIPrintf(viewer,"    Outer relax weight (all) %g\n",(double)jac->outerrelaxweight);

922:     if (jac->relaxorder) {
923:       PetscViewerASCIIPrintf(viewer,"    Using CF-relaxation\n");
924:     } else {
925:       PetscViewerASCIIPrintf(viewer,"    Not using CF-relaxation\n");
926:     }
927:     if (jac->smoothtype!=-1) {
928:       PetscViewerASCIIPrintf(viewer,"    Smooth type          %s\n",HYPREBoomerAMGSmoothType[jac->smoothtype]);
929:       PetscViewerASCIIPrintf(viewer,"    Smooth num levels    %D\n",jac->smoothnumlevels);
930:     } else {
931:       PetscViewerASCIIPrintf(viewer,"    Not using more complex smoothers.\n");
932:     }
933:     if (jac->smoothtype==3) {
934:       PetscViewerASCIIPrintf(viewer,"    Euclid ILU(k) levels %D\n",jac->eu_level);
935:       PetscViewerASCIIPrintf(viewer,"    Euclid ILU(k) drop tolerance %g\n",(double)jac->eu_droptolerance);
936:       PetscViewerASCIIPrintf(viewer,"    Euclid ILU use Block-Jacobi? %D\n",jac->eu_bj);
937:     }
938:     PetscViewerASCIIPrintf(viewer,"    Measure type        %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
939:     PetscViewerASCIIPrintf(viewer,"    Coarsen type        %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
940:     PetscViewerASCIIPrintf(viewer,"    Interpolation type  %s\n",HYPREBoomerAMGInterpType[jac->interptype]);
941:     if (jac->nodal_coarsening) {
942:       PetscViewerASCIIPrintf(viewer,"    Using nodal coarsening (with HYPRE_BOOMERAMGSetNodal() %D\n",jac->nodal_coarsening);
943:     }
944:     if (jac->vec_interp_variant) {
945:       PetscViewerASCIIPrintf(viewer,"    HYPRE_BoomerAMGSetInterpVecVariant() %D\n",jac->vec_interp_variant);
946:       PetscViewerASCIIPrintf(viewer,"    HYPRE_BoomerAMGSetInterpVecQMax() %D\n",jac->vec_interp_qmax);
947:       PetscViewerASCIIPrintf(viewer,"    HYPRE_BoomerAMGSetSmoothInterpVectors() %d\n",jac->vec_interp_smooth);
948:     }
949:     if (jac->nodal_relax) {
950:       PetscViewerASCIIPrintf(viewer,"    Using nodal relaxation via Schwarz smoothing on levels %D\n",jac->nodal_relax_levels);
951:     }
952:   }
953:   return(0);
954: }

956: /* --------------------------------------------------------------------------------------------*/
957: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PetscOptionItems *PetscOptionsObject,PC pc)
958: {
959:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
961:   PetscInt       indx;
962:   PetscBool      flag;
963:   const char     *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};

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

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

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

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

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

983:   PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);
984:   if (flag) {
985:     jac->symt = indx;
986:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
987:   }

989:   PetscOptionsTail();
990:   return(0);
991: }

993: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
994: {
995:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
997:   PetscBool      iascii;
998:   const char     *symt = 0;;

1001:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1002:   if (iascii) {
1003:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails preconditioning\n");
1004:     PetscViewerASCIIPrintf(viewer,"    nlevels %d\n",jac->nlevels);
1005:     PetscViewerASCIIPrintf(viewer,"    threshold %g\n",(double)jac->threshhold);
1006:     PetscViewerASCIIPrintf(viewer,"    filter %g\n",(double)jac->filter);
1007:     PetscViewerASCIIPrintf(viewer,"    load balance %g\n",(double)jac->loadbal);
1008:     PetscViewerASCIIPrintf(viewer,"    reuse nonzero structure %s\n",PetscBools[jac->ruse]);
1009:     PetscViewerASCIIPrintf(viewer,"    print info to screen %s\n",PetscBools[jac->logging]);
1010:     if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
1011:     else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
1012:     else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
1013:     else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
1014:     PetscViewerASCIIPrintf(viewer,"    %s\n",symt);
1015:   }
1016:   return(0);
1017: }
1018: /* --------------------------------------------------------------------------------------------*/
1019: static PetscErrorCode PCSetFromOptions_HYPRE_AMS(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 AMS Options");
1028:   PetscOptionsInt("-pc_hypre_ams_print_level","Debugging output level for AMS","None",jac->as_print,&jac->as_print,&flag);
1029:   if (flag) PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
1030:   PetscOptionsInt("-pc_hypre_ams_max_iter","Maximum number of AMS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
1031:   if (flag) PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1032:   PetscOptionsInt("-pc_hypre_ams_cycle_type","Cycle type for AMS multigrid","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag);
1033:   if (flag) PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1034:   PetscOptionsReal("-pc_hypre_ams_tol","Error tolerance for AMS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
1035:   if (flag) PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
1036:   PetscOptionsInt("-pc_hypre_ams_relax_type","Relaxation type for AMS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
1037:   PetscOptionsInt("-pc_hypre_ams_relax_times","Number of relaxation steps for AMS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
1038:   PetscOptionsReal("-pc_hypre_ams_relax_weight","Relaxation weight for AMS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
1039:   PetscOptionsReal("-pc_hypre_ams_omega","SSOR coefficient for AMS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
1040:   if (flag || flag2 || flag3 || flag4) {
1041:     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(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_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);
1047:   n = 5;
1048:   PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options","AMG options for vector Poisson","None",jac->as_amg_alpha_opts,&n,&flag2);
1049:   if (flag || flag2) {
1050:     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1051:                                                                      jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1052:                                                                      jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1053:                                                                      jac->as_amg_alpha_theta,
1054:                                                                      jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1055:                                                                      jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1056:   }
1057:   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);
1058:   n = 5;
1059:   PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options","AMG options for scalar Poisson solver","None",jac->as_amg_beta_opts,&n,&flag2);
1060:   if (flag || flag2) {
1061:     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1062:                                                                     jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1063:                                                                     jac->as_amg_beta_opts[2],       /* AMG relax_type */
1064:                                                                     jac->as_amg_beta_theta,
1065:                                                                     jac->as_amg_beta_opts[3],       /* AMG interp_type */
1066:                                                                     jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1067:   }
1068:   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);
1069:   if (flag) { /* override HYPRE's default only if the options is used */
1070:     PetscStackCallStandard(HYPRE_AMSSetProjectionFrequency,(jac->hsolver,jac->ams_proj_freq));
1071:   }
1072:   PetscOptionsTail();
1073:   return(0);
1074: }

1076: static PetscErrorCode PCView_HYPRE_AMS(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 AMS preconditioning\n");
1086:     PetscViewerASCIIPrintf(viewer,"    subspace iterations per application %d\n",jac->as_max_iter);
1087:     PetscViewerASCIIPrintf(viewer,"    subspace cycle type %d\n",jac->ams_cycle_type);
1088:     PetscViewerASCIIPrintf(viewer,"    subspace iteration tolerance %g\n",jac->as_tol);
1089:     PetscViewerASCIIPrintf(viewer,"    smoother type %d\n",jac->as_relax_type);
1090:     PetscViewerASCIIPrintf(viewer,"    number of smoothing steps %d\n",jac->as_relax_times);
1091:     PetscViewerASCIIPrintf(viewer,"    smoother weight %g\n",jac->as_relax_weight);
1092:     PetscViewerASCIIPrintf(viewer,"    smoother omega %g\n",jac->as_omega);
1093:     if (jac->alpha_Poisson) {
1094:       PetscViewerASCIIPrintf(viewer,"    vector Poisson solver (passed in by user)\n");
1095:     } else {
1096:       PetscViewerASCIIPrintf(viewer,"    vector Poisson solver (computed) \n");
1097:     }
1098:     PetscViewerASCIIPrintf(viewer,"        boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);
1099:     PetscViewerASCIIPrintf(viewer,"        boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
1100:     PetscViewerASCIIPrintf(viewer,"        boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);
1101:     PetscViewerASCIIPrintf(viewer,"        boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);
1102:     PetscViewerASCIIPrintf(viewer,"        boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
1103:     PetscViewerASCIIPrintf(viewer,"        boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);
1104:     if (!jac->ams_beta_is_zero) {
1105:       if (jac->beta_Poisson) {
1106:         PetscViewerASCIIPrintf(viewer,"    scalar Poisson solver (passed in by user)\n");
1107:       } else {
1108:         PetscViewerASCIIPrintf(viewer,"    scalar Poisson solver (computed) \n");
1109:       }
1110:       PetscViewerASCIIPrintf(viewer,"        boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);
1111:       PetscViewerASCIIPrintf(viewer,"        boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
1112:       PetscViewerASCIIPrintf(viewer,"        boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);
1113:       PetscViewerASCIIPrintf(viewer,"        boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);
1114:       PetscViewerASCIIPrintf(viewer,"        boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
1115:       PetscViewerASCIIPrintf(viewer,"        boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);
1116:       if (jac->ams_beta_is_zero_part) {
1117:         PetscViewerASCIIPrintf(viewer,"        compatible subspace projection frequency %d (-1 HYPRE uses default)\n",jac->ams_proj_freq);
1118:       }
1119:     } else {
1120:       PetscViewerASCIIPrintf(viewer,"    scalar Poisson solver not used (zero-conductivity everywhere) \n");
1121:     }
1122:   }
1123:   return(0);
1124: }

1126: static PetscErrorCode PCSetFromOptions_HYPRE_ADS(PetscOptionItems *PetscOptionsObject,PC pc)
1127: {
1128:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1130:   PetscInt       n;
1131:   PetscBool      flag,flag2,flag3,flag4;

1134:   PetscOptionsHead(PetscOptionsObject,"HYPRE ADS Options");
1135:   PetscOptionsInt("-pc_hypre_ads_print_level","Debugging output level for ADS","None",jac->as_print,&jac->as_print,&flag);
1136:   if (flag) PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
1137:   PetscOptionsInt("-pc_hypre_ads_max_iter","Maximum number of ADS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
1138:   if (flag) PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1139:   PetscOptionsInt("-pc_hypre_ads_cycle_type","Cycle type for ADS multigrid","None",jac->ads_cycle_type,&jac->ads_cycle_type,&flag);
1140:   if (flag) PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ads_cycle_type));
1141:   PetscOptionsReal("-pc_hypre_ads_tol","Error tolerance for ADS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
1142:   if (flag) PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
1143:   PetscOptionsInt("-pc_hypre_ads_relax_type","Relaxation type for ADS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
1144:   PetscOptionsInt("-pc_hypre_ads_relax_times","Number of relaxation steps for ADS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
1145:   PetscOptionsReal("-pc_hypre_ads_relax_weight","Relaxation weight for ADS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
1146:   PetscOptionsReal("-pc_hypre_ads_omega","SSOR coefficient for ADS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
1147:   if (flag || flag2 || flag3 || flag4) {
1148:     PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1149:                                                                       jac->as_relax_times,
1150:                                                                       jac->as_relax_weight,
1151:                                                                       jac->as_omega));
1152:   }
1153:   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);
1154:   n = 5;
1155:   PetscOptionsIntArray("-pc_hypre_ads_ams_options","AMG options for AMS solver inside ADS","None",jac->as_amg_alpha_opts,&n,&flag2);
1156:   PetscOptionsInt("-pc_hypre_ads_ams_cycle_type","Cycle type for AMS solver inside ADS","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag3);
1157:   if (flag || flag2 || flag3) {
1158:     PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type,             /* AMS cycle type */
1159:                                                                 jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1160:                                                                 jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1161:                                                                 jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1162:                                                                 jac->as_amg_alpha_theta,
1163:                                                                 jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1164:                                                                 jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1165:   }
1166:   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);
1167:   n = 5;
1168:   PetscOptionsIntArray("-pc_hypre_ads_amg_options","AMG options for vector AMG solver inside ADS","None",jac->as_amg_beta_opts,&n,&flag2);
1169:   if (flag || flag2) {
1170:     PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1171:                                                                 jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1172:                                                                 jac->as_amg_beta_opts[2],       /* AMG relax_type */
1173:                                                                 jac->as_amg_beta_theta,
1174:                                                                 jac->as_amg_beta_opts[3],       /* AMG interp_type */
1175:                                                                 jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1176:   }
1177:   PetscOptionsTail();
1178:   return(0);
1179: }

1181: static PetscErrorCode PCView_HYPRE_ADS(PC pc,PetscViewer viewer)
1182: {
1183:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1185:   PetscBool      iascii;

1188:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1189:   if (iascii) {
1190:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS preconditioning\n");
1191:     PetscViewerASCIIPrintf(viewer,"    subspace iterations per application %d\n",jac->as_max_iter);
1192:     PetscViewerASCIIPrintf(viewer,"    subspace cycle type %d\n",jac->ads_cycle_type);
1193:     PetscViewerASCIIPrintf(viewer,"    subspace iteration tolerance %g\n",jac->as_tol);
1194:     PetscViewerASCIIPrintf(viewer,"    smoother type %d\n",jac->as_relax_type);
1195:     PetscViewerASCIIPrintf(viewer,"    number of smoothing steps %d\n",jac->as_relax_times);
1196:     PetscViewerASCIIPrintf(viewer,"    smoother weight %g\n",jac->as_relax_weight);
1197:     PetscViewerASCIIPrintf(viewer,"    smoother omega %g\n",jac->as_omega);
1198:     PetscViewerASCIIPrintf(viewer,"    AMS solver using boomerAMG\n");
1199:     PetscViewerASCIIPrintf(viewer,"        subspace cycle type %d\n",jac->ams_cycle_type);
1200:     PetscViewerASCIIPrintf(viewer,"        coarsening type %d\n",jac->as_amg_alpha_opts[0]);
1201:     PetscViewerASCIIPrintf(viewer,"        levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
1202:     PetscViewerASCIIPrintf(viewer,"        relaxation type %d\n",jac->as_amg_alpha_opts[2]);
1203:     PetscViewerASCIIPrintf(viewer,"        interpolation type %d\n",jac->as_amg_alpha_opts[3]);
1204:     PetscViewerASCIIPrintf(viewer,"        max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
1205:     PetscViewerASCIIPrintf(viewer,"        strength threshold %g\n",jac->as_amg_alpha_theta);
1206:     PetscViewerASCIIPrintf(viewer,"    vector Poisson solver using boomerAMG\n");
1207:     PetscViewerASCIIPrintf(viewer,"        coarsening type %d\n",jac->as_amg_beta_opts[0]);
1208:     PetscViewerASCIIPrintf(viewer,"        levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
1209:     PetscViewerASCIIPrintf(viewer,"        relaxation type %d\n",jac->as_amg_beta_opts[2]);
1210:     PetscViewerASCIIPrintf(viewer,"        interpolation type %d\n",jac->as_amg_beta_opts[3]);
1211:     PetscViewerASCIIPrintf(viewer,"        max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
1212:     PetscViewerASCIIPrintf(viewer,"        strength threshold %g\n",jac->as_amg_beta_theta);
1213:   }
1214:   return(0);
1215: }

1217: static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE(PC pc, Mat G)
1218: {
1219:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1220:   PetscBool      ishypre;

1224:   PetscObjectTypeCompare((PetscObject)G,MATHYPRE,&ishypre);
1225:   if (ishypre) {
1226:     PetscObjectReference((PetscObject)G);
1227:     MatDestroy(&jac->G);
1228:     jac->G = G;
1229:   } else {
1230:     MatDestroy(&jac->G);
1231:     MatConvert(G,MATHYPRE,MAT_INITIAL_MATRIX,&jac->G);
1232:   }
1233:   return(0);
1234: }

1236: /*@
1237:  PCHYPRESetDiscreteGradient - Set discrete gradient matrix

1239:    Collective on PC

1241:    Input Parameters:
1242: +  pc - the preconditioning context
1243: -  G - the discrete gradient

1245:    Level: intermediate

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

1251: .seealso:
1252: @*/
1253: PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G)
1254: {

1261:   PetscTryMethod(pc,"PCHYPRESetDiscreteGradient_C",(PC,Mat),(pc,G));
1262:   return(0);
1263: }

1265: static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C)
1266: {
1267:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1268:   PetscBool      ishypre;

1272:   PetscObjectTypeCompare((PetscObject)C,MATHYPRE,&ishypre);
1273:   if (ishypre) {
1274:     PetscObjectReference((PetscObject)C);
1275:     MatDestroy(&jac->C);
1276:     jac->C = C;
1277:   } else {
1278:     MatDestroy(&jac->C);
1279:     MatConvert(C,MATHYPRE,MAT_INITIAL_MATRIX,&jac->C);
1280:   }
1281:   return(0);
1282: }

1284: /*@
1285:  PCHYPRESetDiscreteCurl - Set discrete curl matrix

1287:    Collective on PC

1289:    Input Parameters:
1290: +  pc - the preconditioning context
1291: -  C - the discrete curl

1293:    Level: intermediate

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

1299: .seealso:
1300: @*/
1301: PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C)
1302: {

1309:   PetscTryMethod(pc,"PCHYPRESetDiscreteCurl_C",(PC,Mat),(pc,C));
1310:   return(0);
1311: }

1313: static PetscErrorCode PCHYPRESetInterpolations_HYPRE(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[])
1314: {
1315:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1316:   PetscBool      ishypre;
1318:   PetscInt       i;

1321:   MatDestroy(&jac->RT_PiFull);
1322:   MatDestroy(&jac->ND_PiFull);
1323:   for (i=0;i<3;++i) {
1324:     MatDestroy(&jac->RT_Pi[i]);
1325:     MatDestroy(&jac->ND_Pi[i]);
1326:   }

1328:   jac->dim = dim;
1329:   if (RT_PiFull) {
1330:     PetscObjectTypeCompare((PetscObject)RT_PiFull,MATHYPRE,&ishypre);
1331:     if (ishypre) {
1332:       PetscObjectReference((PetscObject)RT_PiFull);
1333:       jac->RT_PiFull = RT_PiFull;
1334:     } else {
1335:       MatConvert(RT_PiFull,MATHYPRE,MAT_INITIAL_MATRIX,&jac->RT_PiFull);
1336:     }
1337:   }
1338:   if (RT_Pi) {
1339:     for (i=0;i<dim;++i) {
1340:       if (RT_Pi[i]) {
1341:         PetscObjectTypeCompare((PetscObject)RT_Pi[i],MATHYPRE,&ishypre);
1342:         if (ishypre) {
1343:           PetscObjectReference((PetscObject)RT_Pi[i]);
1344:           jac->RT_Pi[i] = RT_Pi[i];
1345:         } else {
1346:           MatConvert(RT_Pi[i],MATHYPRE,MAT_INITIAL_MATRIX,&jac->RT_Pi[i]);
1347:         }
1348:       }
1349:     }
1350:   }
1351:   if (ND_PiFull) {
1352:     PetscObjectTypeCompare((PetscObject)ND_PiFull,MATHYPRE,&ishypre);
1353:     if (ishypre) {
1354:       PetscObjectReference((PetscObject)ND_PiFull);
1355:       jac->ND_PiFull = ND_PiFull;
1356:     } else {
1357:       MatConvert(ND_PiFull,MATHYPRE,MAT_INITIAL_MATRIX,&jac->ND_PiFull);
1358:     }
1359:   }
1360:   if (ND_Pi) {
1361:     for (i=0;i<dim;++i) {
1362:       if (ND_Pi[i]) {
1363:         PetscObjectTypeCompare((PetscObject)ND_Pi[i],MATHYPRE,&ishypre);
1364:         if (ishypre) {
1365:           PetscObjectReference((PetscObject)ND_Pi[i]);
1366:           jac->ND_Pi[i] = ND_Pi[i];
1367:         } else {
1368:           MatConvert(ND_Pi[i],MATHYPRE,MAT_INITIAL_MATRIX,&jac->ND_Pi[i]);
1369:         }
1370:       }
1371:     }
1372:   }

1374:   return(0);
1375: }

1377: /*@
1378:  PCHYPRESetInterpolations - Set interpolation matrices for AMS/ADS preconditioner

1380:    Collective on PC

1382:    Input Parameters:
1383: +  pc - the preconditioning context
1384: -  dim - the dimension of the problem, only used in AMS
1385: -  RT_PiFull - Raviart-Thomas interpolation matrix
1386: -  RT_Pi - x/y/z component of Raviart-Thomas interpolation matrix
1387: -  ND_PiFull - Nedelec interpolation matrix
1388: -  ND_Pi - x/y/z component of Nedelec interpolation matrix

1390:    Notes:
1391:     For AMS, only Nedelec interpolation matrices are needed, the Raviart-Thomas interpolation matrices can be set to NULL.
1392:           For ADS, both type of interpolation matrices are needed.
1393:    Level: intermediate

1395: .seealso:
1396: @*/
1397: PetscErrorCode PCHYPRESetInterpolations(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[])
1398: {
1400:   PetscInt       i;

1404:   if (RT_PiFull) {
1407:   }
1408:   if (RT_Pi) {
1410:     for (i=0;i<dim;++i) {
1411:       if (RT_Pi[i]) {
1414:       }
1415:     }
1416:   }
1417:   if (ND_PiFull) {
1420:   }
1421:   if (ND_Pi) {
1423:     for (i=0;i<dim;++i) {
1424:       if (ND_Pi[i]) {
1427:       }
1428:     }
1429:   }
1430:   PetscTryMethod(pc,"PCHYPRESetInterpolations_C",(PC,PetscInt,Mat,Mat[],Mat,Mat[]),(pc,dim,RT_PiFull,RT_Pi,ND_PiFull,ND_Pi));
1431:   return(0);
1432: }

1434: static PetscErrorCode PCHYPRESetPoissonMatrix_HYPRE(PC pc, Mat A, PetscBool isalpha)
1435: {
1436:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1437:   PetscBool      ishypre;

1441:   PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
1442:   if (ishypre) {
1443:     if (isalpha) {
1444:       PetscObjectReference((PetscObject)A);
1445:       MatDestroy(&jac->alpha_Poisson);
1446:       jac->alpha_Poisson = A;
1447:     } else {
1448:       if (A) {
1449:         PetscObjectReference((PetscObject)A);
1450:       } else {
1451:         jac->ams_beta_is_zero = PETSC_TRUE;
1452:       }
1453:       MatDestroy(&jac->beta_Poisson);
1454:       jac->beta_Poisson = A;
1455:     }
1456:   } else {
1457:     if (isalpha) {
1458:       MatDestroy(&jac->alpha_Poisson);
1459:       MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&jac->alpha_Poisson);
1460:     } else {
1461:       if (A) {
1462:         MatDestroy(&jac->beta_Poisson);
1463:         MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&jac->beta_Poisson);
1464:       } else {
1465:         MatDestroy(&jac->beta_Poisson);
1466:         jac->ams_beta_is_zero = PETSC_TRUE;
1467:       }
1468:     }
1469:   }
1470:   return(0);
1471: }

1473: /*@
1474:  PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix

1476:    Collective on PC

1478:    Input Parameters:
1479: +  pc - the preconditioning context
1480: -  A - the matrix

1482:    Level: intermediate

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

1487: .seealso:
1488: @*/
1489: PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A)
1490: {

1497:   PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_TRUE));
1498:   return(0);
1499: }

1501: /*@
1502:  PCHYPRESetBetaPoissonMatrix - Set Poisson matrix

1504:    Collective on PC

1506:    Input Parameters:
1507: +  pc - the preconditioning context
1508: -  A - the matrix

1510:    Level: intermediate

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

1516: .seealso:
1517: @*/
1518: PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A)
1519: {

1524:   if (A) {
1527:   }
1528:   PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_FALSE));
1529:   return(0);
1530: }

1532: static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE(PC pc,Vec ozz, Vec zoz, Vec zzo)
1533: {
1534:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1535:   PetscErrorCode     ierr;

1538:   /* throw away any vector if already set */
1539:   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0]));
1540:   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1]));
1541:   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2]));
1542:   jac->constants[0] = NULL;
1543:   jac->constants[1] = NULL;
1544:   jac->constants[2] = NULL;
1545:   VecHYPRE_IJVectorCreate(ozz,&jac->constants[0]);
1546:   VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);
1547:   VecHYPRE_IJVectorCreate(zoz,&jac->constants[1]);
1548:   VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);
1549:   jac->dim = 2;
1550:   if (zzo) {
1551:     VecHYPRE_IJVectorCreate(zzo,&jac->constants[2]);
1552:     VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);
1553:     jac->dim++;
1554:   }
1555:   return(0);
1556: }

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

1561:    Collective on PC

1563:    Input Parameters:
1564: +  pc - the preconditioning context
1565: -  ozz - vector representing (1,0,0) (or (1,0) in 2D)
1566: -  zoz - vector representing (0,1,0) (or (0,1) in 2D)
1567: -  zzo - vector representing (0,0,1) (use NULL in 2D)

1569:    Level: intermediate

1571:    Notes:

1573: .seealso:
1574: @*/
1575: PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo)
1576: {

1587:   PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));
1588:   return(0);
1589: }

1591: static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1592: {
1593:   PC_HYPRE        *jac = (PC_HYPRE*)pc->data;
1594:   Vec             tv;
1595:   PetscInt        i;
1596:   PetscErrorCode  ierr;

1599:   /* throw away any coordinate vector if already set */
1600:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0]));
1601:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1]));
1602:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2]));
1603:   jac->dim = dim;

1605:   /* compute IJ vector for coordinates */
1606:   VecCreate(PetscObjectComm((PetscObject)pc),&tv);
1607:   VecSetType(tv,VECSTANDARD);
1608:   VecSetSizes(tv,nloc,PETSC_DECIDE);
1609:   for (i=0;i<dim;i++) {
1610:     PetscScalar *array;
1611:     PetscInt    j;

1613:     VecHYPRE_IJVectorCreate(tv,&jac->coords[i]);
1614:     VecGetArray(tv,&array);
1615:     for (j=0;j<nloc;j++) {
1616:       array[j] = coords[j*dim+i];
1617:     }
1618:     PetscStackCallStandard(HYPRE_IJVectorSetValues,(jac->coords[i],nloc,NULL,array));
1619:     PetscStackCallStandard(HYPRE_IJVectorAssemble,(jac->coords[i]));
1620:     VecRestoreArray(tv,&array);
1621:   }
1622:   VecDestroy(&tv);
1623:   return(0);
1624: }

1626: /* ---------------------------------------------------------------------------------*/

1628: static PetscErrorCode  PCHYPREGetType_HYPRE(PC pc,const char *name[])
1629: {
1630:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

1633:   *name = jac->hypre_type;
1634:   return(0);
1635: }

1637: static PetscErrorCode  PCHYPRESetType_HYPRE(PC pc,const char name[])
1638: {
1639:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1641:   PetscBool      flag;

1644:   if (jac->hypre_type) {
1645:     PetscStrcmp(jac->hypre_type,name,&flag);
1646:     if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set");
1647:     return(0);
1648:   } else {
1649:     PetscStrallocpy(name, &jac->hypre_type);
1650:   }

1652:   jac->maxiter         = PETSC_DEFAULT;
1653:   jac->tol             = PETSC_DEFAULT;
1654:   jac->printstatistics = PetscLogPrintInfo;

1656:   PetscStrcmp("pilut",jac->hypre_type,&flag);
1657:   if (flag) {
1658:     MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1659:     PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver));
1660:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
1661:     pc->ops->view           = PCView_HYPRE_Pilut;
1662:     jac->destroy            = HYPRE_ParCSRPilutDestroy;
1663:     jac->setup              = HYPRE_ParCSRPilutSetup;
1664:     jac->solve              = HYPRE_ParCSRPilutSolve;
1665:     jac->factorrowsize      = PETSC_DEFAULT;
1666:     return(0);
1667:   }
1668:   PetscStrcmp("euclid",jac->hypre_type,&flag);
1669:   if (flag) {
1670:     MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1671:     PetscStackCallStandard(HYPRE_EuclidCreate,(jac->comm_hypre,&jac->hsolver));
1672:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
1673:     pc->ops->view           = PCView_HYPRE_Euclid;
1674:     jac->destroy            = HYPRE_EuclidDestroy;
1675:     jac->setup              = HYPRE_EuclidSetup;
1676:     jac->solve              = HYPRE_EuclidSolve;
1677:     jac->factorrowsize      = PETSC_DEFAULT;
1678:     jac->eu_level           = PETSC_DEFAULT; /* default */
1679:     return(0);
1680:   }
1681:   PetscStrcmp("parasails",jac->hypre_type,&flag);
1682:   if (flag) {
1683:     MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1684:     PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver));
1685:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
1686:     pc->ops->view           = PCView_HYPRE_ParaSails;
1687:     jac->destroy            = HYPRE_ParaSailsDestroy;
1688:     jac->setup              = HYPRE_ParaSailsSetup;
1689:     jac->solve              = HYPRE_ParaSailsSolve;
1690:     /* initialize */
1691:     jac->nlevels    = 1;
1692:     jac->threshhold = .1;
1693:     jac->filter     = .1;
1694:     jac->loadbal    = 0;
1695:     if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE;
1696:     else jac->logging = (int) PETSC_FALSE;

1698:     jac->ruse = (int) PETSC_FALSE;
1699:     jac->symt = 0;
1700:     PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
1701:     PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
1702:     PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
1703:     PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
1704:     PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
1705:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
1706:     return(0);
1707:   }
1708:   PetscStrcmp("boomeramg",jac->hypre_type,&flag);
1709:   if (flag) {
1710:     HYPRE_BoomerAMGCreate(&jac->hsolver);
1711:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_BoomerAMG;
1712:     pc->ops->view            = PCView_HYPRE_BoomerAMG;
1713:     pc->ops->applytranspose  = PCApplyTranspose_HYPRE_BoomerAMG;
1714:     pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
1715:     jac->destroy             = HYPRE_BoomerAMGDestroy;
1716:     jac->setup               = HYPRE_BoomerAMGSetup;
1717:     jac->solve               = HYPRE_BoomerAMGSolve;
1718:     jac->applyrichardson     = PETSC_FALSE;
1719:     /* these defaults match the hypre defaults */
1720:     jac->cycletype        = 1;
1721:     jac->maxlevels        = 25;
1722:     jac->maxiter          = 1;
1723:     jac->tol              = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
1724:     jac->truncfactor      = 0.0;
1725:     jac->strongthreshold  = .25;
1726:     jac->maxrowsum        = .9;
1727:     jac->coarsentype      = 6;
1728:     jac->measuretype      = 0;
1729:     jac->gridsweeps[0]    = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
1730:     jac->smoothtype       = -1; /* Not set by default */
1731:     jac->smoothnumlevels  = 25;
1732:     jac->eu_level         = 0;
1733:     jac->eu_droptolerance = 0;
1734:     jac->eu_bj            = 0;
1735:     jac->relaxtype[0]     = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */
1736:     jac->relaxtype[2]     = 9; /*G.E. */
1737:     jac->relaxweight      = 1.0;
1738:     jac->outerrelaxweight = 1.0;
1739:     jac->relaxorder       = 1;
1740:     jac->interptype       = 0;
1741:     jac->agg_nl           = 0;
1742:     jac->pmax             = 0;
1743:     jac->truncfactor      = 0.0;
1744:     jac->agg_num_paths    = 1;

1746:     jac->nodal_coarsening      = 0;
1747:     jac->nodal_coarsening_diag = 0;
1748:     jac->vec_interp_variant    = 0;
1749:     jac->vec_interp_qmax       = 0;
1750:     jac->vec_interp_smooth     = PETSC_FALSE;
1751:     jac->interp_refine         = 0;
1752:     jac->nodal_relax           = PETSC_FALSE;
1753:     jac->nodal_relax_levels    = 1;
1754:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
1755:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
1756:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
1757:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
1758:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
1759:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
1760:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
1761:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
1762:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
1763:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
1764:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
1765:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
1766:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
1767:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
1768:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0]));  /*defaults coarse to 9*/
1769:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */
1770:     return(0);
1771:   }
1772:   PetscStrcmp("ams",jac->hypre_type,&flag);
1773:   if (flag) {
1774:     HYPRE_AMSCreate(&jac->hsolver);
1775:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_AMS;
1776:     pc->ops->view            = PCView_HYPRE_AMS;
1777:     jac->destroy             = HYPRE_AMSDestroy;
1778:     jac->setup               = HYPRE_AMSSetup;
1779:     jac->solve               = HYPRE_AMSSolve;
1780:     jac->coords[0]           = NULL;
1781:     jac->coords[1]           = NULL;
1782:     jac->coords[2]           = NULL;
1783:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */
1784:     jac->as_print           = 0;
1785:     jac->as_max_iter        = 1; /* used as a preconditioner */
1786:     jac->as_tol             = 0.; /* used as a preconditioner */
1787:     jac->ams_cycle_type     = 13;
1788:     /* Smoothing options */
1789:     jac->as_relax_type      = 2;
1790:     jac->as_relax_times     = 1;
1791:     jac->as_relax_weight    = 1.0;
1792:     jac->as_omega           = 1.0;
1793:     /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1794:     jac->as_amg_alpha_opts[0] = 10;
1795:     jac->as_amg_alpha_opts[1] = 1;
1796:     jac->as_amg_alpha_opts[2] = 6;
1797:     jac->as_amg_alpha_opts[3] = 6;
1798:     jac->as_amg_alpha_opts[4] = 4;
1799:     jac->as_amg_alpha_theta   = 0.25;
1800:     /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1801:     jac->as_amg_beta_opts[0] = 10;
1802:     jac->as_amg_beta_opts[1] = 1;
1803:     jac->as_amg_beta_opts[2] = 6;
1804:     jac->as_amg_beta_opts[3] = 6;
1805:     jac->as_amg_beta_opts[4] = 4;
1806:     jac->as_amg_beta_theta   = 0.25;
1807:     PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
1808:     PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1809:     PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1810:     PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
1811:     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1812:                                                                       jac->as_relax_times,
1813:                                                                       jac->as_relax_weight,
1814:                                                                       jac->as_omega));
1815:     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1816:                                                                      jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1817:                                                                      jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1818:                                                                      jac->as_amg_alpha_theta,
1819:                                                                      jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1820:                                                                      jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1821:     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1822:                                                                     jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1823:                                                                     jac->as_amg_beta_opts[2],       /* AMG relax_type */
1824:                                                                     jac->as_amg_beta_theta,
1825:                                                                     jac->as_amg_beta_opts[3],       /* AMG interp_type */
1826:                                                                     jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1827:     /* Zero conductivity */
1828:     jac->ams_beta_is_zero      = PETSC_FALSE;
1829:     jac->ams_beta_is_zero_part = PETSC_FALSE;
1830:     return(0);
1831:   }
1832:   PetscStrcmp("ads",jac->hypre_type,&flag);
1833:   if (flag) {
1834:     HYPRE_ADSCreate(&jac->hsolver);
1835:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_ADS;
1836:     pc->ops->view            = PCView_HYPRE_ADS;
1837:     jac->destroy             = HYPRE_ADSDestroy;
1838:     jac->setup               = HYPRE_ADSSetup;
1839:     jac->solve               = HYPRE_ADSSolve;
1840:     jac->coords[0]           = NULL;
1841:     jac->coords[1]           = NULL;
1842:     jac->coords[2]           = NULL;
1843:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */
1844:     jac->as_print           = 0;
1845:     jac->as_max_iter        = 1; /* used as a preconditioner */
1846:     jac->as_tol             = 0.; /* used as a preconditioner */
1847:     jac->ads_cycle_type     = 13;
1848:     /* Smoothing options */
1849:     jac->as_relax_type      = 2;
1850:     jac->as_relax_times     = 1;
1851:     jac->as_relax_weight    = 1.0;
1852:     jac->as_omega           = 1.0;
1853:     /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */
1854:     jac->ams_cycle_type       = 14;
1855:     jac->as_amg_alpha_opts[0] = 10;
1856:     jac->as_amg_alpha_opts[1] = 1;
1857:     jac->as_amg_alpha_opts[2] = 6;
1858:     jac->as_amg_alpha_opts[3] = 6;
1859:     jac->as_amg_alpha_opts[4] = 4;
1860:     jac->as_amg_alpha_theta   = 0.25;
1861:     /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1862:     jac->as_amg_beta_opts[0] = 10;
1863:     jac->as_amg_beta_opts[1] = 1;
1864:     jac->as_amg_beta_opts[2] = 6;
1865:     jac->as_amg_beta_opts[3] = 6;
1866:     jac->as_amg_beta_opts[4] = 4;
1867:     jac->as_amg_beta_theta   = 0.25;
1868:     PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
1869:     PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1870:     PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1871:     PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
1872:     PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1873:                                                                       jac->as_relax_times,
1874:                                                                       jac->as_relax_weight,
1875:                                                                       jac->as_omega));
1876:     PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type,             /* AMG coarsen type */
1877:                                                                 jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1878:                                                                 jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1879:                                                                 jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1880:                                                                 jac->as_amg_alpha_theta,
1881:                                                                 jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1882:                                                                 jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1883:     PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1884:                                                                 jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1885:                                                                 jac->as_amg_beta_opts[2],       /* AMG relax_type */
1886:                                                                 jac->as_amg_beta_theta,
1887:                                                                 jac->as_amg_beta_opts[3],       /* AMG interp_type */
1888:                                                                 jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1889:     return(0);
1890:   }
1891:   PetscFree(jac->hypre_type);

1893:   jac->hypre_type = NULL;
1894:   SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are euclid, pilut, parasails, boomeramg, ams",name);
1895:   return(0);
1896: }

1898: /*
1899:     It only gets here if the HYPRE type has not been set before the call to
1900:    ...SetFromOptions() which actually is most of the time
1901: */
1902: static PetscErrorCode PCSetFromOptions_HYPRE(PetscOptionItems *PetscOptionsObject,PC pc)
1903: {
1905:   PetscInt       indx;
1906:   const char     *type[] = {"euclid","pilut","parasails","boomeramg","ams","ads"};
1907:   PetscBool      flg;

1910:   PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");
1911:   PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"boomeramg",&indx,&flg);
1912:   if (flg) {
1913:     PCHYPRESetType_HYPRE(pc,type[indx]);
1914:   } else {
1915:     PCHYPRESetType_HYPRE(pc,"boomeramg");
1916:   }
1917:   if (pc->ops->setfromoptions) {
1918:     pc->ops->setfromoptions(PetscOptionsObject,pc);
1919:   }
1920:   PetscOptionsTail();
1921:   return(0);
1922: }

1924: /*@C
1925:      PCHYPRESetType - Sets which hypre preconditioner you wish to use

1927:    Input Parameters:
1928: +     pc - the preconditioner context
1929: -     name - either  euclid, pilut, parasails, boomeramg, ams, ads

1931:    Options Database Keys:
1932:    -pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads

1934:    Level: intermediate

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

1939: @*/
1940: PetscErrorCode  PCHYPRESetType(PC pc,const char name[])
1941: {

1947:   PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));
1948:   return(0);
1949: }

1951: /*@C
1952:      PCHYPREGetType - Gets which hypre preconditioner you are using

1954:    Input Parameter:
1955: .     pc - the preconditioner context

1957:    Output Parameter:
1958: .     name - either  euclid, pilut, parasails, boomeramg, ams, ads

1960:    Level: intermediate

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

1965: @*/
1966: PetscErrorCode  PCHYPREGetType(PC pc,const char *name[])
1967: {

1973:   PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));
1974:   return(0);
1975: }

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

1980:    Options Database Keys:
1981: +   -pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads
1982: -   Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX
1983:           preconditioner

1985:    Level: intermediate

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

1992:           The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_tol refer to the number of iterations
1993:           (V-cycles) and tolerance that boomeramg does EACH time it is called. So for example, if
1994:           -pc_hypre_boomeramg_max_iter is set to 2 then 2-V-cycles are being used to define the preconditioner
1995:           (-pc_hypre_boomeramg_tol should be set to 0.0 - the default - to strictly use a fixed number of
1996:           iterations per hypre call). -ksp_max_it and -ksp_rtol STILL determine the total number of iterations
1997:           and tolerance for the Krylov solver. For example, if -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10
1998:           then AT MOST twenty V-cycles of boomeramg will be called.

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

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

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

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

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

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

2022: M*/

2024: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
2025: {
2026:   PC_HYPRE       *jac;

2030:   PetscNewLog(pc,&jac);

2032:   pc->data                = jac;
2033:   pc->ops->reset          = PCReset_HYPRE;
2034:   pc->ops->destroy        = PCDestroy_HYPRE;
2035:   pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
2036:   pc->ops->setup          = PCSetUp_HYPRE;
2037:   pc->ops->apply          = PCApply_HYPRE;
2038:   jac->comm_hypre         = MPI_COMM_NULL;
2039:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);
2040:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);
2041:   PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);
2042:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);
2043:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",PCHYPRESetDiscreteCurl_HYPRE);
2044:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",PCHYPRESetInterpolations_HYPRE);
2045:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE);
2046:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",PCHYPRESetPoissonMatrix_HYPRE);
2047:   return(0);
2048: }

2050: /* ---------------------------------------------------------------------------------------------------------------------------------*/

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

2056:   /* keep copy of PFMG options used so may view them */
2057:   PetscInt its;
2058:   double   tol;
2059:   PetscInt relax_type;
2060:   PetscInt rap_type;
2061:   PetscInt num_pre_relax,num_post_relax;
2062:   PetscInt max_levels;
2063: } PC_PFMG;

2065: PetscErrorCode PCDestroy_PFMG(PC pc)
2066: {
2068:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

2071:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
2072:   MPI_Comm_free(&ex->hcomm);
2073:   PetscFree(pc->data);
2074:   return(0);
2075: }

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

2080: PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
2081: {
2083:   PetscBool      iascii;
2084:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

2087:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2088:   if (iascii) {
2089:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG preconditioning\n");
2090:     PetscViewerASCIIPrintf(viewer,"    max iterations %d\n",ex->its);
2091:     PetscViewerASCIIPrintf(viewer,"    tolerance %g\n",ex->tol);
2092:     PetscViewerASCIIPrintf(viewer,"    relax type %s\n",PFMGRelaxType[ex->relax_type]);
2093:     PetscViewerASCIIPrintf(viewer,"    RAP type %s\n",PFMGRAPType[ex->rap_type]);
2094:     PetscViewerASCIIPrintf(viewer,"    number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
2095:     PetscViewerASCIIPrintf(viewer,"    max levels %d\n",ex->max_levels);
2096:   }
2097:   return(0);
2098: }

2100: PetscErrorCode PCSetFromOptions_PFMG(PetscOptionItems *PetscOptionsObject,PC pc)
2101: {
2103:   PC_PFMG        *ex = (PC_PFMG*) pc->data;
2104:   PetscBool      flg = PETSC_FALSE;

2107:   PetscOptionsHead(PetscOptionsObject,"PFMG options");
2108:   PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);
2109:   if (flg) {
2110:     int level=3;
2111:     PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level));
2112:   }
2113:   PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);
2114:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its));
2115:   PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
2116:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax));
2117:   PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
2118:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax));

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

2123:   PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);
2124:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol));
2125:   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);
2126:   PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type));
2127:   PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);
2128:   PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type));
2129:   PetscOptionsTail();
2130:   return(0);
2131: }

2133: PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
2134: {
2135:   PetscErrorCode    ierr;
2136:   PC_PFMG           *ex = (PC_PFMG*) pc->data;
2137:   PetscScalar       *yy;
2138:   const PetscScalar *xx;
2139:   PetscInt          ilower[3],iupper[3];
2140:   Mat_HYPREStruct   *mx = (Mat_HYPREStruct*)(pc->pmat->data);

2143:   PetscCitationsRegister(hypreCitation,&cite);
2144:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
2145:   iupper[0] += ilower[0] - 1;
2146:   iupper[1] += ilower[1] - 1;
2147:   iupper[2] += ilower[2] - 1;

2149:   /* copy x values over to hypre */
2150:   PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
2151:   VecGetArrayRead(x,&xx);
2152:   PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx));
2153:   VecRestoreArrayRead(x,&xx);
2154:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
2155:   PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx));

2157:   /* copy solution values back to PETSc */
2158:   VecGetArray(y,&yy);
2159:   PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy));
2160:   VecRestoreArray(y,&yy);
2161:   return(0);
2162: }

2164: 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)
2165: {
2166:   PC_PFMG        *jac = (PC_PFMG*)pc->data;
2168:   PetscInt       oits;

2171:   PetscCitationsRegister(hypreCitation,&cite);
2172:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its));
2173:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol));

2175:   PCApply_PFMG(pc,b,y);
2176:   PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
2177:   *outits = oits;
2178:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2179:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
2180:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol));
2181:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its));
2182:   return(0);
2183: }


2186: PetscErrorCode PCSetUp_PFMG(PC pc)
2187: {
2188:   PetscErrorCode  ierr;
2189:   PC_PFMG         *ex = (PC_PFMG*) pc->data;
2190:   Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
2191:   PetscBool       flg;

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

2197:   /* create the hypre solver object and set its information */
2198:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
2199:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
2200:   PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
2201:   PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver));
2202:   return(0);
2203: }

2205: /*MC
2206:      PCPFMG - the hypre PFMG multigrid solver

2208:    Level: advanced

2210:    Options Database:
2211: + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner
2212: . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2213: . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2214: . -pc_pfmg_tol <tol> tolerance of PFMG
2215: . -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
2216: - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin

2218:    Notes:
2219:     This is for CELL-centered descretizations

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

2224: .seealso:  PCMG, MATHYPRESTRUCT
2225: M*/

2227: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
2228: {
2230:   PC_PFMG        *ex;

2233:   PetscNew(&ex); \
2234:   pc->data = ex;

2236:   ex->its            = 1;
2237:   ex->tol            = 1.e-8;
2238:   ex->relax_type     = 1;
2239:   ex->rap_type       = 0;
2240:   ex->num_pre_relax  = 1;
2241:   ex->num_post_relax = 1;
2242:   ex->max_levels     = 0;

2244:   pc->ops->setfromoptions  = PCSetFromOptions_PFMG;
2245:   pc->ops->view            = PCView_PFMG;
2246:   pc->ops->destroy         = PCDestroy_PFMG;
2247:   pc->ops->apply           = PCApply_PFMG;
2248:   pc->ops->applyrichardson = PCApplyRichardson_PFMG;
2249:   pc->ops->setup           = PCSetUp_PFMG;

2251:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
2252:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
2253:   return(0);
2254: }

2256: /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/

2258: /* we know we are working with a HYPRE_SStructMatrix */
2259: typedef struct {
2260:   MPI_Comm            hcomm;       /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
2261:   HYPRE_SStructSolver ss_solver;

2263:   /* keep copy of SYSPFMG options used so may view them */
2264:   PetscInt its;
2265:   double   tol;
2266:   PetscInt relax_type;
2267:   PetscInt num_pre_relax,num_post_relax;
2268: } PC_SysPFMG;

2270: PetscErrorCode PCDestroy_SysPFMG(PC pc)
2271: {
2273:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

2276:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2277:   MPI_Comm_free(&ex->hcomm);
2278:   PetscFree(pc->data);
2279:   return(0);
2280: }

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

2284: PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
2285: {
2287:   PetscBool      iascii;
2288:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

2291:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2292:   if (iascii) {
2293:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG preconditioning\n");
2294:     PetscViewerASCIIPrintf(viewer,"  max iterations %d\n",ex->its);
2295:     PetscViewerASCIIPrintf(viewer,"  tolerance %g\n",ex->tol);
2296:     PetscViewerASCIIPrintf(viewer,"  relax type %s\n",PFMGRelaxType[ex->relax_type]);
2297:     PetscViewerASCIIPrintf(viewer,"  number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
2298:   }
2299:   return(0);
2300: }

2302: PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptionItems *PetscOptionsObject,PC pc)
2303: {
2305:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;
2306:   PetscBool      flg = PETSC_FALSE;

2309:   PetscOptionsHead(PetscOptionsObject,"SysPFMG options");
2310:   PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);
2311:   if (flg) {
2312:     int level=3;
2313:     PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level));
2314:   }
2315:   PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);
2316:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its));
2317:   PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
2318:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax));
2319:   PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
2320:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax));

2322:   PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);
2323:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol));
2324:   PetscOptionsEList("-pc_syspfmg_relax_type","Relax type for the up and down cycles","HYPRE_SStructSysPFMGSetRelaxType",SysPFMGRelaxType,ALEN(SysPFMGRelaxType),SysPFMGRelaxType[ex->relax_type],&ex->relax_type,NULL);
2325:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type));
2326:   PetscOptionsTail();
2327:   return(0);
2328: }

2330: PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
2331: {
2332:   PetscErrorCode    ierr;
2333:   PC_SysPFMG        *ex = (PC_SysPFMG*) pc->data;
2334:   PetscScalar       *yy;
2335:   const PetscScalar *xx;
2336:   PetscInt          ilower[3],iupper[3];
2337:   Mat_HYPRESStruct  *mx     = (Mat_HYPRESStruct*)(pc->pmat->data);
2338:   PetscInt          ordering= mx->dofs_order;
2339:   PetscInt          nvars   = mx->nvars;
2340:   PetscInt          part    = 0;
2341:   PetscInt          size;
2342:   PetscInt          i;

2345:   PetscCitationsRegister(hypreCitation,&cite);
2346:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
2347:   iupper[0] += ilower[0] - 1;
2348:   iupper[1] += ilower[1] - 1;
2349:   iupper[2] += ilower[2] - 1;

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

2354:   /* copy x values over to hypre for variable ordering */
2355:   if (ordering) {
2356:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2357:     VecGetArrayRead(x,&xx);
2358:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)xx+(size*i)));
2359:     VecRestoreArrayRead(x,&xx);
2360:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2361:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
2362:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

2364:     /* copy solution values back to PETSc */
2365:     VecGetArray(y,&yy);
2366:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i)));
2367:     VecRestoreArray(y,&yy);
2368:   } else {      /* nodal ordering must be mapped to variable ordering for sys_pfmg */
2369:     PetscScalar *z;
2370:     PetscInt    j, k;

2372:     PetscMalloc1(nvars*size,&z);
2373:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2374:     VecGetArrayRead(x,&xx);

2376:     /* transform nodal to hypre's variable ordering for sys_pfmg */
2377:     for (i= 0; i< size; i++) {
2378:       k= i*nvars;
2379:       for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
2380:     }
2381:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
2382:     VecRestoreArrayRead(x,&xx);
2383:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2384:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

2386:     /* copy solution values back to PETSc */
2387:     VecGetArray(y,&yy);
2388:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
2389:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
2390:     for (i= 0; i< size; i++) {
2391:       k= i*nvars;
2392:       for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
2393:     }
2394:     VecRestoreArray(y,&yy);
2395:     PetscFree(z);
2396:   }
2397:   return(0);
2398: }

2400: 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)
2401: {
2402:   PC_SysPFMG     *jac = (PC_SysPFMG*)pc->data;
2404:   PetscInt       oits;

2407:   PetscCitationsRegister(hypreCitation,&cite);
2408:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its));
2409:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol));
2410:   PCApply_SysPFMG(pc,b,y);
2411:   PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,(HYPRE_Int *)&oits));
2412:   *outits = oits;
2413:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2414:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
2415:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol));
2416:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its));
2417:   return(0);
2418: }

2420: PetscErrorCode PCSetUp_SysPFMG(PC pc)
2421: {
2422:   PetscErrorCode   ierr;
2423:   PC_SysPFMG       *ex = (PC_SysPFMG*) pc->data;
2424:   Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
2425:   PetscBool        flg;

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

2431:   /* create the hypre sstruct solver object and set its information */
2432:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2433:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2434:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver));
2435:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2436:   return(0);
2437: }

2439: /*MC
2440:      PCSysPFMG - the hypre SysPFMG multigrid solver

2442:    Level: advanced

2444:    Options Database:
2445: + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner
2446: . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2447: . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2448: . -pc_syspfmg_tol <tol> tolerance of SysPFMG
2449: . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel

2451:    Notes:
2452:     This is for CELL-centered descretizations

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

2458: .seealso:  PCMG, MATHYPRESSTRUCT
2459: M*/

2461: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
2462: {
2464:   PC_SysPFMG     *ex;

2467:   PetscNew(&ex); \
2468:   pc->data = ex;

2470:   ex->its            = 1;
2471:   ex->tol            = 1.e-8;
2472:   ex->relax_type     = 1;
2473:   ex->num_pre_relax  = 1;
2474:   ex->num_post_relax = 1;

2476:   pc->ops->setfromoptions  = PCSetFromOptions_SysPFMG;
2477:   pc->ops->view            = PCView_SysPFMG;
2478:   pc->ops->destroy         = PCDestroy_SysPFMG;
2479:   pc->ops->apply           = PCApply_SysPFMG;
2480:   pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
2481:   pc->ops->setup           = PCSetUp_SysPFMG;

2483:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
2484:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2485:   return(0);
2486: }