Actual source code: hypre.c

  1: /*
  2:    Provides an interface to the LLNL package hypre
  3: */

  5: #include <petscpkg_version.h>
  6: #include <petsc/private/pcimpl.h>
  7: /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */
  8: #include <petsc/private/matimpl.h>
  9: #include <petsc/private/vecimpl.h>
 10: #include <../src/vec/vec/impls/hypre/vhyp.h>
 11: #include <../src/mat/impls/hypre/mhypre.h>
 12: #include <../src/dm/impls/da/hypre/mhyp.h>
 13: #include <_hypre_parcsr_ls.h>
 14: #include <petscmathypre.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{https://computation.llnl.gov/projects/hypre-scalable-linear-solvers-multigrid-methods}}\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:   PetscReal tol;

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

 40:   /* options for ParaSails */
 41:   PetscInt  nlevels;
 42:   PetscReal threshold;
 43:   PetscReal filter;
 44:   PetscReal loadbal;
 45:   PetscInt  logging;
 46:   PetscInt  ruse;
 47:   PetscInt  symt;

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

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

 76:   /* GPU */
 77:   PetscBool keeptranspose;
 78:   PetscInt  rap2;
 79:   PetscInt  mod_rap2;

 81:   /* AIR */
 82:   PetscInt  Rtype;
 83:   PetscReal Rstrongthreshold;
 84:   PetscReal Rfilterthreshold;
 85:   PetscInt  Adroptype;
 86:   PetscReal Adroptol;

 88:   PetscInt  agg_nl;
 89:   PetscInt  agg_interptype;
 90:   PetscInt  agg_num_paths;
 91:   PetscBool nodal_relax;
 92:   PetscInt  nodal_relax_levels;

 94:   PetscInt  nodal_coarsening;
 95:   PetscInt  nodal_coarsening_diag;
 96:   PetscInt  vec_interp_variant;
 97:   PetscInt  vec_interp_qmax;
 98:   PetscBool vec_interp_smooth;
 99:   PetscInt  interp_refine;

101:   /* NearNullSpace support */
102:   VecHYPRE_IJVector *hmnull;
103:   HYPRE_ParVector   *phmnull;
104:   PetscInt          n_hmnull;
105:   Vec               hmnull_constant;

107:   /* options for AS (Auxiliary Space preconditioners) */
108:   PetscInt  as_print;
109:   PetscInt  as_max_iter;
110:   PetscReal as_tol;
111:   PetscInt  as_relax_type;
112:   PetscInt  as_relax_times;
113:   PetscReal as_relax_weight;
114:   PetscReal as_omega;
115:   PetscInt  as_amg_alpha_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for vector Poisson (AMS) or Curl problem (ADS) */
116:   PetscReal as_amg_alpha_theta;   /* AMG strength for vector Poisson (AMS) or Curl problem (ADS) */
117:   PetscInt  as_amg_beta_opts[5];  /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for scalar Poisson (AMS) or vector Poisson (ADS) */
118:   PetscReal as_amg_beta_theta;    /* AMG strength for scalar Poisson (AMS) or vector Poisson (ADS)  */
119:   PetscInt  ams_cycle_type;
120:   PetscInt  ads_cycle_type;

122:   /* additional data */
123:   Mat G;             /* MatHYPRE */
124:   Mat C;             /* MatHYPRE */
125:   Mat alpha_Poisson; /* MatHYPRE */
126:   Mat beta_Poisson;  /* MatHYPRE */

128:   /* extra information for AMS */
129:   PetscInt          dim; /* geometrical dimension */
130:   VecHYPRE_IJVector coords[3];
131:   VecHYPRE_IJVector constants[3];
132:   Mat               RT_PiFull, RT_Pi[3];
133:   Mat               ND_PiFull, ND_Pi[3];
134:   PetscBool         ams_beta_is_zero;
135:   PetscBool         ams_beta_is_zero_part;
136:   PetscInt          ams_proj_freq;
137: } PC_HYPRE;

139: PetscErrorCode PCHYPREGetSolver(PC pc,HYPRE_Solver *hsolver)
140: {
141:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

144:   *hsolver = jac->hsolver;
145:   return(0);
146: }

148: /*
149:   Matrices with AIJ format are created IN PLACE with using (I,J,data) from BoomerAMG. Since the data format in hypre_ParCSRMatrix
150:   is different from that used in PETSc, the original hypre_ParCSRMatrix can not be used any more after call this routine.
151:   It is used in PCHMG. Other users should avoid using this function.
152: */
153: static PetscErrorCode PCGetCoarseOperators_BoomerAMG(PC pc,PetscInt *nlevels,Mat *operators[])
154: {
155:   PC_HYPRE             *jac  = (PC_HYPRE*)pc->data;
156:   PetscBool            same = PETSC_FALSE;
157:   PetscErrorCode       ierr;
158:   PetscInt             num_levels,l;
159:   Mat                  *mattmp;
160:   hypre_ParCSRMatrix   **A_array;

163:   PetscStrcmp(jac->hypre_type,"boomeramg",&same);
164:   if (!same) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_NOTSAMETYPE,"Hypre type is not BoomerAMG \n");
165:   num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData*) (jac->hsolver));
166:   PetscMalloc1(num_levels,&mattmp);
167:   A_array    = hypre_ParAMGDataAArray((hypre_ParAMGData*) (jac->hsolver));
168:   for (l=1; l<num_levels; l++) {
169:     MatCreateFromParCSR(A_array[l],MATAIJ,PETSC_OWN_POINTER, &(mattmp[num_levels-1-l]));
170:     /* We want to own the data, and HYPRE can not touch this matrix any more */
171:     A_array[l] = NULL;
172:   }
173:   *nlevels = num_levels;
174:   *operators = mattmp;
175:   return(0);
176: }

178: /*
179:   Matrices with AIJ format are created IN PLACE with using (I,J,data) from BoomerAMG. Since the data format in hypre_ParCSRMatrix
180:   is different from that used in PETSc, the original hypre_ParCSRMatrix can not be used any more after call this routine.
181:   It is used in PCHMG. Other users should avoid using this function.
182: */
183: static PetscErrorCode PCGetInterpolations_BoomerAMG(PC pc,PetscInt *nlevels,Mat *interpolations[])
184: {
185:   PC_HYPRE             *jac  = (PC_HYPRE*)pc->data;
186:   PetscBool            same = PETSC_FALSE;
187:   PetscErrorCode       ierr;
188:   PetscInt             num_levels,l;
189:   Mat                  *mattmp;
190:   hypre_ParCSRMatrix   **P_array;

193:   PetscStrcmp(jac->hypre_type,"boomeramg",&same);
194:   if (!same) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_NOTSAMETYPE,"Hypre type is not BoomerAMG \n");
195:   num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData*) (jac->hsolver));
196:   PetscMalloc1(num_levels,&mattmp);
197:   P_array  = hypre_ParAMGDataPArray((hypre_ParAMGData*) (jac->hsolver));
198:   for (l=1; l<num_levels; l++) {
199:     MatCreateFromParCSR(P_array[num_levels-1-l],MATAIJ,PETSC_OWN_POINTER, &(mattmp[l-1]));
200:     /* We want to own the data, and HYPRE can not touch this matrix any more */
201:     P_array[num_levels-1-l] = NULL;
202:   }
203:   *nlevels = num_levels;
204:   *interpolations = mattmp;
205:   return(0);
206: }

208: /* Resets (frees) Hypre's representation of the near null space */
209: static PetscErrorCode PCHYPREResetNearNullSpace_Private(PC pc)
210: {
211:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
212:   PetscInt       i;

216:   for (i=0; i<jac->n_hmnull; i++) {
217:     VecHYPRE_IJVectorDestroy(&jac->hmnull[i]);
218:   }
219:   PetscFree(jac->hmnull);
220:   PetscFree(jac->phmnull);
221:   VecDestroy(&jac->hmnull_constant);
222:   jac->n_hmnull = 0;
223:   return(0);
224: }

226: static PetscErrorCode PCSetUp_HYPRE(PC pc)
227: {
228:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
229:   Mat_HYPRE          *hjac;
230:   HYPRE_ParCSRMatrix hmat;
231:   HYPRE_ParVector    bv,xv;
232:   PetscBool          ishypre;
233:   PetscErrorCode     ierr;

236:   if (!jac->hypre_type) {
237:     PCHYPRESetType(pc,"boomeramg");
238:   }

240:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRE,&ishypre);
241:   if (!ishypre) {
242:     MatDestroy(&jac->hpmat);
243:     MatConvert(pc->pmat,MATHYPRE,MAT_INITIAL_MATRIX,&jac->hpmat);
244:     PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->hpmat);
245:   } else {
246:     PetscObjectReference((PetscObject)pc->pmat);
247:     MatDestroy(&jac->hpmat);
248:     jac->hpmat = pc->pmat;
249:   }
250:   /* allow debug */
251:   MatViewFromOptions(jac->hpmat,NULL,"-pc_hypre_mat_view");
252:   hjac = (Mat_HYPRE*)(jac->hpmat->data);

254:   /* special case for BoomerAMG */
255:   if (jac->setup == HYPRE_BoomerAMGSetup) {
256:     MatNullSpace mnull;
257:     PetscBool    has_const;
258:     PetscInt     bs,nvec,i;
259:     const Vec    *vecs;

261:     MatGetBlockSize(pc->pmat,&bs);
262:     if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
263:     MatGetNearNullSpace(pc->mat, &mnull);
264:     if (mnull) {
265:       PCHYPREResetNearNullSpace_Private(pc);
266:       MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);
267:       PetscMalloc1(nvec+1,&jac->hmnull);
268:       PetscMalloc1(nvec+1,&jac->phmnull);
269:       for (i=0; i<nvec; i++) {
270:         VecHYPRE_IJVectorCreate(vecs[i]->map,&jac->hmnull[i]);
271:         VecHYPRE_IJVectorCopy(vecs[i],jac->hmnull[i]);
272:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[i]->ij,(void**)&jac->phmnull[i]));
273:       }
274:       if (has_const) {
275:         MatCreateVecs(pc->pmat,&jac->hmnull_constant,NULL);
276:         PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->hmnull_constant);
277:         VecSet(jac->hmnull_constant,1);
278:         VecNormalize(jac->hmnull_constant,NULL);
279:         VecHYPRE_IJVectorCreate(jac->hmnull_constant->map,&jac->hmnull[nvec]);
280:         VecHYPRE_IJVectorCopy(jac->hmnull_constant,jac->hmnull[nvec]);
281:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[nvec]->ij,(void**)&jac->phmnull[nvec]));
282:         nvec++;
283:       }
284:       PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVectors,(jac->hsolver,nvec,jac->phmnull));
285:       jac->n_hmnull = nvec;
286:     }
287:   }

289:   /* special case for AMS */
290:   if (jac->setup == HYPRE_AMSSetup) {
291:     Mat_HYPRE          *hm;
292:     HYPRE_ParCSRMatrix parcsr;
293:     if (!jac->coords[0] && !jac->constants[0] && !(jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]))) {
294:       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");
295:     }
296:     if (jac->dim) {
297:       PetscStackCallStandard(HYPRE_AMSSetDimension,(jac->hsolver,jac->dim));
298:     }
299:     if (jac->constants[0]) {
300:       HYPRE_ParVector ozz,zoz,zzo = NULL;
301:       PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[0]->ij,(void**)(&ozz)));
302:       PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[1]->ij,(void**)(&zoz)));
303:       if (jac->constants[2]) {
304:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[2]->ij,(void**)(&zzo)));
305:       }
306:       PetscStackCallStandard(HYPRE_AMSSetEdgeConstantVectors,(jac->hsolver,ozz,zoz,zzo));
307:     }
308:     if (jac->coords[0]) {
309:       HYPRE_ParVector coords[3];
310:       coords[0] = NULL;
311:       coords[1] = NULL;
312:       coords[2] = NULL;
313:       if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0]->ij,(void**)(&coords[0])));
314:       if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1]->ij,(void**)(&coords[1])));
315:       if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2]->ij,(void**)(&coords[2])));
316:       PetscStackCallStandard(HYPRE_AMSSetCoordinateVectors,(jac->hsolver,coords[0],coords[1],coords[2]));
317:     }
318:     if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient");
319:     hm = (Mat_HYPRE*)(jac->G->data);
320:     PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
321:     PetscStackCallStandard(HYPRE_AMSSetDiscreteGradient,(jac->hsolver,parcsr));
322:     if (jac->alpha_Poisson) {
323:       hm = (Mat_HYPRE*)(jac->alpha_Poisson->data);
324:       PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
325:       PetscStackCallStandard(HYPRE_AMSSetAlphaPoissonMatrix,(jac->hsolver,parcsr));
326:     }
327:     if (jac->ams_beta_is_zero) {
328:       PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,NULL));
329:     } else if (jac->beta_Poisson) {
330:       hm = (Mat_HYPRE*)(jac->beta_Poisson->data);
331:       PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
332:       PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,parcsr));
333:     }
334:     if (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1])) {
335:       PetscInt           i;
336:       HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3];
337:       if (jac->ND_PiFull) {
338:         hm = (Mat_HYPRE*)(jac->ND_PiFull->data);
339:         PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&nd_parcsrfull)));
340:       } else {
341:         nd_parcsrfull = NULL;
342:       }
343:       for (i=0;i<3;++i) {
344:         if (jac->ND_Pi[i]) {
345:           hm = (Mat_HYPRE*)(jac->ND_Pi[i]->data);
346:           PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&nd_parcsr[i])));
347:         } else {
348:           nd_parcsr[i] = NULL;
349:         }
350:       }
351:       PetscStackCallStandard(HYPRE_AMSSetInterpolations,(jac->hsolver,nd_parcsrfull,nd_parcsr[0],nd_parcsr[1],nd_parcsr[2]));
352:     }
353:   }
354:   /* special case for ADS */
355:   if (jac->setup == HYPRE_ADSSetup) {
356:     Mat_HYPRE          *hm;
357:     HYPRE_ParCSRMatrix parcsr;
358:     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])))) {
359:       SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs either the coordinate vectors via PCSetCoordinates() or the interpolation matrices via PCHYPRESetInterpolations");
360:     }
361:     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");
362:     if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient");
363:     if (!jac->C) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs the discrete curl operator via PCHYPRESetDiscreteGradient");
364:     if (jac->coords[0]) {
365:       HYPRE_ParVector coords[3];
366:       coords[0] = NULL;
367:       coords[1] = NULL;
368:       coords[2] = NULL;
369:       if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0]->ij,(void**)(&coords[0])));
370:       if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1]->ij,(void**)(&coords[1])));
371:       if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2]->ij,(void**)(&coords[2])));
372:       PetscStackCallStandard(HYPRE_ADSSetCoordinateVectors,(jac->hsolver,coords[0],coords[1],coords[2]));
373:     }
374:     hm = (Mat_HYPRE*)(jac->G->data);
375:     PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
376:     PetscStackCallStandard(HYPRE_ADSSetDiscreteGradient,(jac->hsolver,parcsr));
377:     hm = (Mat_HYPRE*)(jac->C->data);
378:     PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
379:     PetscStackCallStandard(HYPRE_ADSSetDiscreteCurl,(jac->hsolver,parcsr));
380:     if ((jac->RT_PiFull || (jac->RT_Pi[0] && jac->RT_Pi[1])) && (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]))) {
381:       PetscInt           i;
382:       HYPRE_ParCSRMatrix rt_parcsrfull, rt_parcsr[3];
383:       HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3];
384:       if (jac->RT_PiFull) {
385:         hm = (Mat_HYPRE*)(jac->RT_PiFull->data);
386:         PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&rt_parcsrfull)));
387:       } else {
388:         rt_parcsrfull = NULL;
389:       }
390:       for (i=0;i<3;++i) {
391:         if (jac->RT_Pi[i]) {
392:           hm = (Mat_HYPRE*)(jac->RT_Pi[i]->data);
393:           PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&rt_parcsr[i])));
394:         } else {
395:           rt_parcsr[i] = NULL;
396:         }
397:       }
398:       if (jac->ND_PiFull) {
399:         hm = (Mat_HYPRE*)(jac->ND_PiFull->data);
400:         PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&nd_parcsrfull)));
401:       } else {
402:         nd_parcsrfull = NULL;
403:       }
404:       for (i=0;i<3;++i) {
405:         if (jac->ND_Pi[i]) {
406:           hm = (Mat_HYPRE*)(jac->ND_Pi[i]->data);
407:           PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&nd_parcsr[i])));
408:         } else {
409:           nd_parcsr[i] = NULL;
410:         }
411:       }
412:       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]));
413:     }
414:   }
415:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
416:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b->ij,(void**)&bv));
417:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x->ij,(void**)&xv));
418:   PetscStackCallStandard(jac->setup,(jac->hsolver,hmat,bv,xv));
419:   return(0);
420: }

422: static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
423: {
424:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
425:   Mat_HYPRE          *hjac = (Mat_HYPRE*)(jac->hpmat->data);
426:   PetscErrorCode     ierr;
427:   HYPRE_ParCSRMatrix hmat;
428:   HYPRE_ParVector    jbv,jxv;
429:   PetscInt           hierr;

432:   PetscCitationsRegister(hypreCitation,&cite);
433:   if (!jac->applyrichardson) {VecSet(x,0.0);}
434:   VecHYPRE_IJVectorPushVecRead(hjac->b,b);
435:   if (jac->applyrichardson) { VecHYPRE_IJVectorPushVec(hjac->x,x); }
436:   else { VecHYPRE_IJVectorPushVecWrite(hjac->x,x); }
437:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
438:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b->ij,(void**)&jbv));
439:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x->ij,(void**)&jxv));
440:   PetscStackCall("Hypre solve",h(*jac->solve)(jac->hsolver,hmat,jbv,jxv);
441:   if (hierr && hierr != HYPRE_ERROR_CONV) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
442:   if (hierr) hypre__global_error = 0;);

444:   if (jac->setup == HYPRE_AMSSetup && jac->ams_beta_is_zero_part) {
445:     PetscStackCallStandard(HYPRE_AMSProjectOutGradients,(jac->hsolver,jxv));
446:   }
447:   VecHYPRE_IJVectorPopVec(hjac->x);
448:   VecHYPRE_IJVectorPopVec(hjac->b);
449:   return(0);
450: }

452: static PetscErrorCode PCReset_HYPRE(PC pc)
453: {
454:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;

458:   MatDestroy(&jac->hpmat);
459:   MatDestroy(&jac->G);
460:   MatDestroy(&jac->C);
461:   MatDestroy(&jac->alpha_Poisson);
462:   MatDestroy(&jac->beta_Poisson);
463:   MatDestroy(&jac->RT_PiFull);
464:   MatDestroy(&jac->RT_Pi[0]);
465:   MatDestroy(&jac->RT_Pi[1]);
466:   MatDestroy(&jac->RT_Pi[2]);
467:   MatDestroy(&jac->ND_PiFull);
468:   MatDestroy(&jac->ND_Pi[0]);
469:   MatDestroy(&jac->ND_Pi[1]);
470:   MatDestroy(&jac->ND_Pi[2]);
471:   VecHYPRE_IJVectorDestroy(&jac->coords[0]);
472:   VecHYPRE_IJVectorDestroy(&jac->coords[1]);
473:   VecHYPRE_IJVectorDestroy(&jac->coords[2]);
474:   VecHYPRE_IJVectorDestroy(&jac->constants[0]);
475:   VecHYPRE_IJVectorDestroy(&jac->constants[1]);
476:   VecHYPRE_IJVectorDestroy(&jac->constants[2]);
477:   PCHYPREResetNearNullSpace_Private(pc);
478:   jac->ams_beta_is_zero = PETSC_FALSE;
479:   jac->dim = 0;
480:   return(0);
481: }

483: static PetscErrorCode PCDestroy_HYPRE(PC pc)
484: {
485:   PC_HYPRE                 *jac = (PC_HYPRE*)pc->data;
486:   PetscErrorCode           ierr;

489:   PCReset_HYPRE(pc);
490:   if (jac->destroy) PetscStackCallStandard(jac->destroy,(jac->hsolver));
491:   PetscFree(jac->hypre_type);
492:   if (jac->comm_hypre != MPI_COMM_NULL) {MPI_Comm_free(&(jac->comm_hypre));}
493:   PetscFree(pc->data);

495:   PetscObjectChangeTypeName((PetscObject)pc,0);
496:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);
497:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);
498:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetCoordinates_C",NULL);
499:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",NULL);
500:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",NULL);
501:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",NULL);
502:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetConstantEdgeVectors_C",NULL);
503:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",NULL);
504:   PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",NULL);
505:   PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",NULL);
506:   return(0);
507: }

509: /* --------------------------------------------------------------------------------------------*/
510: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PetscOptionItems *PetscOptionsObject,PC pc)
511: {
512:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
514:   PetscBool      flag;

517:   PetscOptionsHead(PetscOptionsObject,"HYPRE Pilut Options");
518:   PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
519:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter));
520:   PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
521:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol));
522:   PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
523:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize));
524:   PetscOptionsTail();
525:   return(0);
526: }

528: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
529: {
530:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
532:   PetscBool      iascii;

535:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
536:   if (iascii) {
537:     PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut preconditioning\n");
538:     if (jac->maxiter != PETSC_DEFAULT) {
539:       PetscViewerASCIIPrintf(viewer,"    maximum number of iterations %d\n",jac->maxiter);
540:     } else {
541:       PetscViewerASCIIPrintf(viewer,"    default maximum number of iterations \n");
542:     }
543:     if (jac->tol != PETSC_DEFAULT) {
544:       PetscViewerASCIIPrintf(viewer,"    drop tolerance %g\n",(double)jac->tol);
545:     } else {
546:       PetscViewerASCIIPrintf(viewer,"    default drop tolerance \n");
547:     }
548:     if (jac->factorrowsize != PETSC_DEFAULT) {
549:       PetscViewerASCIIPrintf(viewer,"    factor row size %d\n",jac->factorrowsize);
550:     } else {
551:       PetscViewerASCIIPrintf(viewer,"    default factor row size \n");
552:     }
553:   }
554:   return(0);
555: }

557: /* --------------------------------------------------------------------------------------------*/
558: static PetscErrorCode PCSetFromOptions_HYPRE_Euclid(PetscOptionItems *PetscOptionsObject,PC pc)
559: {
560:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
562:   PetscBool      flag,eu_bj = jac->eu_bj ? PETSC_TRUE : PETSC_FALSE;

565:   PetscOptionsHead(PetscOptionsObject,"HYPRE Euclid Options");
566:   PetscOptionsInt("-pc_hypre_euclid_level","Factorization levels","None",jac->eu_level,&jac->eu_level,&flag);
567:   if (flag) PetscStackCallStandard(HYPRE_EuclidSetLevel,(jac->hsolver,jac->eu_level));

569:   PetscOptionsReal("-pc_hypre_euclid_droptolerance","Drop tolerance for ILU(k) in Euclid","None",jac->eu_droptolerance,&jac->eu_droptolerance,&flag);
570:   if (flag) {
571:     PetscMPIInt size;

573:     MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
574:     if (size > 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"hypre's Euclid does not support a parallel drop tolerance");
575:     PetscStackCallStandard(HYPRE_EuclidSetILUT,(jac->hsolver,jac->eu_droptolerance));
576:   }

578:   PetscOptionsBool("-pc_hypre_euclid_bj", "Use Block Jacobi for ILU in Euclid", "None", eu_bj,&eu_bj,&flag);
579:   if (flag) {
580:     jac->eu_bj = eu_bj ? 1 : 0;
581:     PetscStackCallStandard(HYPRE_EuclidSetBJ,(jac->hsolver,jac->eu_bj));
582:   }
583:   PetscOptionsTail();
584:   return(0);
585: }

587: static PetscErrorCode PCView_HYPRE_Euclid(PC pc,PetscViewer viewer)
588: {
589:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
591:   PetscBool      iascii;

594:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
595:   if (iascii) {
596:     PetscViewerASCIIPrintf(viewer,"  HYPRE Euclid preconditioning\n");
597:     if (jac->eu_level != PETSC_DEFAULT) {
598:       PetscViewerASCIIPrintf(viewer,"    factorization levels %d\n",jac->eu_level);
599:     } else {
600:       PetscViewerASCIIPrintf(viewer,"    default factorization levels \n");
601:     }
602:     PetscViewerASCIIPrintf(viewer,"    drop tolerance %g\n",(double)jac->eu_droptolerance);
603:     PetscViewerASCIIPrintf(viewer,"    use Block-Jacobi? %D\n",jac->eu_bj);
604:   }
605:   return(0);
606: }

608: /* --------------------------------------------------------------------------------------------*/

610: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
611: {
612:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
613:   Mat_HYPRE          *hjac = (Mat_HYPRE*)(jac->hpmat->data);
614:   PetscErrorCode     ierr;
615:   HYPRE_ParCSRMatrix hmat;
616:   HYPRE_ParVector    jbv,jxv;
617:   PetscInt           hierr;

620:   PetscCitationsRegister(hypreCitation,&cite);
621:   VecSet(x,0.0);
622:   VecHYPRE_IJVectorPushVecRead(hjac->x,b);
623:   VecHYPRE_IJVectorPushVecWrite(hjac->b,x);

625:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
626:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b->ij,(void**)&jbv));
627:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x->ij,(void**)&jxv));

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

634:   VecHYPRE_IJVectorPopVec(hjac->x);
635:   VecHYPRE_IJVectorPopVec(hjac->b);
636:   return(0);
637: }

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

642: static const char *HYPREBoomerAMGCycleType[]   = {"","V","W"};
643: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
644: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
645: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
646: static const char *HYPREBoomerAMGSmoothType[]  = {"Schwarz-smoothers","Pilut","ParaSails","Euclid"};
647: static const char *HYPREBoomerAMGRelaxType[]   = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi",
648:                                                   "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi",
649:                                                   "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination",
650:                                                   "" /* 10 */, "" /* 11 */, "" /* 12 */, "l1-Gauss-Seidel" /* nonsymmetric */, "backward-l1-Gauss-Seidel" /* nonsymmetric */,
651:                                                   "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"};
652: static const char *HYPREBoomerAMGInterpType[]  = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
653:                                                   "ext+i-cc", "standard", "standard-wts", "block", "block-wtd", "FF", "FF1",
654:                                                   "ext", "ad-wts", "ext-mm", "ext+i-mm", "ext+e-mm"};
655: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PetscOptionItems *PetscOptionsObject,PC pc)
656: {
657:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
659:   PetscInt       bs,n,indx,level;
660:   PetscBool      flg, tmp_truth;
661:   double         tmpdbl, twodbl[2];
662:   const char     *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};

665:   PetscOptionsHead(PetscOptionsObject,"HYPRE BoomerAMG Options");
666:   PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);
667:   if (flg) {
668:     jac->cycletype = indx+1;
669:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
670:   }
671:   PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
672:   if (flg) {
673:     if (jac->maxlevels < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
674:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
675:   }
676:   PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
677:   if (flg) {
678:     if (jac->maxiter < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
679:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
680:   }
681:   PetscOptionsReal("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);
682:   if (flg) {
683:     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);
684:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
685:   }
686:   bs = 1;
687:   if (pc->pmat) {
688:     MatGetBlockSize(pc->pmat,&bs);
689:   }
690:   PetscOptionsInt("-pc_hypre_boomeramg_numfunctions","Number of functions","HYPRE_BoomerAMGSetNumFunctions",bs,&bs,&flg);
691:   if (flg) {
692:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
693:   }

695:   PetscOptionsReal("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);
696:   if (flg) {
697:     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);
698:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
699:   }

701:   PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);
702:   if (flg) {
703:     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);
704:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
705:   }

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

711:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
712:   }

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

718:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
719:   }

721:   PetscOptionsReal("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
722:   if (flg) {
723:     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);
724:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
725:   }
726:   PetscOptionsReal("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
727:   if (flg) {
728:     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);
729:     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);
730:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
731:   }

733:   /* Grid sweeps */
734:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);
735:   if (flg) {
736:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx));
737:     /* modify the jac structure so we can view the updated options with PC_View */
738:     jac->gridsweeps[0] = indx;
739:     jac->gridsweeps[1] = indx;
740:     /*defaults coarse to 1 */
741:     jac->gridsweeps[2] = 1;
742:   }
743:   PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen","Use a nodal based coarsening 1-6","HYPRE_BoomerAMGSetNodal",jac->nodal_coarsening,&jac->nodal_coarsening,&flg);
744:   if (flg) {
745:     PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,jac->nodal_coarsening));
746:   }
747:   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);
748:   if (flg) {
749:     PetscStackCallStandard(HYPRE_BoomerAMGSetNodalDiag,(jac->hsolver,jac->nodal_coarsening_diag));
750:   }
751:   PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_variant","Variant of algorithm 1-3","HYPRE_BoomerAMGSetInterpVecVariant",jac->vec_interp_variant, &jac->vec_interp_variant,&flg);
752:   if (flg) {
753:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecVariant,(jac->hsolver,jac->vec_interp_variant));
754:   }
755:   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);
756:   if (flg) {
757:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecQMax,(jac->hsolver,jac->vec_interp_qmax));
758:   }
759:   PetscOptionsBool("-pc_hypre_boomeramg_vec_interp_smooth","Whether to smooth the interpolation vectors","HYPRE_BoomerAMGSetSmoothInterpVectors",jac->vec_interp_smooth, &jac->vec_interp_smooth,&flg);
760:   if (flg) {
761:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothInterpVectors,(jac->hsolver,jac->vec_interp_smooth));
762:   }
763:   PetscOptionsInt("-pc_hypre_boomeramg_interp_refine","Preprocess the interpolation matrix through iterative weight refinement","HYPRE_BoomerAMGSetInterpRefine",jac->interp_refine, &jac->interp_refine,&flg);
764:   if (flg) {
765:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpRefine,(jac->hsolver,jac->interp_refine));
766:   }
767:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);
768:   if (flg) {
769:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1));
770:     jac->gridsweeps[0] = indx;
771:   }
772:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);
773:   if (flg) {
774:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2));
775:     jac->gridsweeps[1] = indx;
776:   }
777:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);
778:   if (flg) {
779:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3));
780:     jac->gridsweeps[2] = indx;
781:   }

783:   /* Smooth type */
784:   PetscOptionsEList("-pc_hypre_boomeramg_smooth_type","Enable more complex smoothers","None",HYPREBoomerAMGSmoothType,ALEN(HYPREBoomerAMGSmoothType),HYPREBoomerAMGSmoothType[0],&indx,&flg);
785:   if (flg) {
786:     jac->smoothtype = indx;
787:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,indx+6));
788:     jac->smoothnumlevels = 25;
789:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,25));
790:   }

792:   /* Number of smoothing levels */
793:   PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_levels","Number of levels on which more complex smoothers are used","None",25,&indx,&flg);
794:   if (flg && (jac->smoothtype != -1)) {
795:     jac->smoothnumlevels = indx;
796:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,indx));
797:   }

799:   /* Number of levels for ILU(k) for Euclid */
800:   PetscOptionsInt("-pc_hypre_boomeramg_eu_level","Number of levels for ILU(k) in Euclid smoother","None",0,&indx,&flg);
801:   if (flg && (jac->smoothtype == 3)) {
802:     jac->eu_level = indx;
803:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,indx));
804:   }

806:   /* Filter for ILU(k) for Euclid */
807:   double droptolerance;
808:   PetscOptionsReal("-pc_hypre_boomeramg_eu_droptolerance","Drop tolerance for ILU(k) in Euclid smoother","None",0,&droptolerance,&flg);
809:   if (flg && (jac->smoothtype == 3)) {
810:     jac->eu_droptolerance = droptolerance;
811:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,droptolerance));
812:   }

814:   /* Use Block Jacobi ILUT for Euclid */
815:   PetscOptionsBool("-pc_hypre_boomeramg_eu_bj", "Use Block Jacobi for ILU in Euclid smoother?", "None", PETSC_FALSE, &tmp_truth, &flg);
816:   if (flg && (jac->smoothtype == 3)) {
817:     jac->eu_bj = tmp_truth;
818:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuBJ,(jac->hsolver,jac->eu_bj));
819:   }

821:   /* Relax type */
822:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
823:   if (flg) {
824:     jac->relaxtype[0] = jac->relaxtype[1]  = indx;
825:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx));
826:     /* by default, coarse type set to 9 */
827:     jac->relaxtype[2] = 9;
828:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, 9, 3));
829:   }
830:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
831:   if (flg) {
832:     jac->relaxtype[0] = indx;
833:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1));
834:   }
835:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
836:   if (flg) {
837:     jac->relaxtype[1] = indx;
838:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2));
839:   }
840:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);
841:   if (flg) {
842:     jac->relaxtype[2] = indx;
843:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3));
844:   }

846:   /* Relaxation Weight */
847:   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);
848:   if (flg) {
849:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl));
850:     jac->relaxweight = tmpdbl;
851:   }

853:   n         = 2;
854:   twodbl[0] = twodbl[1] = 1.0;
855:   PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
856:   if (flg) {
857:     if (n == 2) {
858:       indx =  (int)PetscAbsReal(twodbl[1]);
859:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx));
860:     } 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);
861:   }

863:   /* Outer relaxation Weight */
864:   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);
865:   if (flg) {
866:     PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl));
867:     jac->outerrelaxweight = tmpdbl;
868:   }

870:   n         = 2;
871:   twodbl[0] = twodbl[1] = 1.0;
872:   PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
873:   if (flg) {
874:     if (n == 2) {
875:       indx =  (int)PetscAbsReal(twodbl[1]);
876:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx));
877:     } 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);
878:   }

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

883:   if (flg && tmp_truth) {
884:     jac->relaxorder = 0;
885:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
886:   }
887:   PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);
888:   if (flg) {
889:     jac->measuretype = indx;
890:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
891:   }
892:   /* update list length 3/07 */
893:   PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);
894:   if (flg) {
895:     jac->coarsentype = indx;
896:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
897:   }

899:   PetscOptionsInt("-pc_hypre_boomeramg_max_coarse_size", "Maximum size of coarsest grid", "None", jac->maxc, &jac->maxc, &flg);
900:   if (flg) {
901:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxCoarseSize,(jac->hsolver, jac->maxc));
902:   }
903:   PetscOptionsInt("-pc_hypre_boomeramg_min_coarse_size", "Minimum size of coarsest grid", "None", jac->minc, &jac->minc, &flg);
904:   if (flg) {
905:     PetscStackCallStandard(HYPRE_BoomerAMGSetMinCoarseSize,(jac->hsolver, jac->minc));
906:   }

908:   /* AIR */
909: #if PETSC_PKG_HYPRE_VERSION_GE(2,18,0)
910:   PetscOptionsInt("-pc_hypre_boomeramg_restriction_type", "Type of AIR method (distance 1 or 2, 0 means no AIR)", "None", jac->Rtype, &jac->Rtype, NULL);
911:   PetscStackCallStandard(HYPRE_BoomerAMGSetRestriction,(jac->hsolver,jac->Rtype));
912:   if (jac->Rtype) {
913:     jac->interptype = 100; /* no way we can pass this with strings... Set it as default as in MFEM, then users can still customize it back to a different one */

915:     PetscOptionsReal("-pc_hypre_boomeramg_strongthresholdR","Threshold for R","None",jac->Rstrongthreshold,&jac->Rstrongthreshold,NULL);
916:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThresholdR,(jac->hsolver,jac->Rstrongthreshold));

918:     PetscOptionsReal("-pc_hypre_boomeramg_filterthresholdR","Filter threshold for R","None",jac->Rfilterthreshold,&jac->Rfilterthreshold,NULL);
919:     PetscStackCallStandard(HYPRE_BoomerAMGSetFilterThresholdR,(jac->hsolver,jac->Rfilterthreshold));

921:     PetscOptionsReal("-pc_hypre_boomeramg_Adroptol","Defines the drop tolerance for the A-matrices from the 2nd level of AMG","None",jac->Adroptol,&jac->Adroptol,NULL);
922:     PetscStackCallStandard(HYPRE_BoomerAMGSetADropTol,(jac->hsolver,jac->Adroptol));

924:     PetscOptionsInt("-pc_hypre_boomeramg_Adroptype","Drops the entries that are not on the diagonal and smaller than its row norm: type 1: 1-norm, 2: 2-norm, -1: infinity norm","None",jac->Adroptype,&jac->Adroptype,NULL);
925:     PetscStackCallStandard(HYPRE_BoomerAMGSetADropType,(jac->hsolver,jac->Adroptype));
926:   }
927: #endif

929:   /* new 3/07 */
930:   PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);
931:   if (flg || jac->Rtype) {
932:     if (flg) jac->interptype = indx;
933:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
934:   }

936:   PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
937:   if (flg) {
938:     level = 3;
939:     PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);

941:     jac->printstatistics = PETSC_TRUE;
942:     PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level));
943:   }

945:   PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);
946:   if (flg) {
947:     level = 3;
948:     PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);

950:     jac->printstatistics = PETSC_TRUE;
951:     PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level));
952:   }

954:   PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);
955:   if (flg && tmp_truth) {
956:     PetscInt tmp_int;
957:     PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);
958:     if (flg) jac->nodal_relax_levels = tmp_int;
959:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6));
960:     PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1));
961:     PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0));
962:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels));
963:   }

965:   PetscOptionsBool("-pc_hypre_boomeramg_keeptranspose", "Avoid transpose matvecs in preconditioner application", "None", jac->keeptranspose, &jac->keeptranspose, NULL);
966:   PetscStackCallStandard(HYPRE_BoomerAMGSetKeepTranspose,(jac->hsolver,jac->keeptranspose ? 1 : 0));

968:   /* options for ParaSails solvers */
969:   PetscOptionsEList("-pc_hypre_boomeramg_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flg);
970:   if (flg) {
971:     jac->symt = indx;
972:     PetscStackCallStandard(HYPRE_BoomerAMGSetSym,(jac->hsolver,jac->symt));
973:   }

975:   PetscOptionsTail();
976:   return(0);
977: }

979: 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)
980: {
981:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
983:   HYPRE_Int      oits;

986:   PetscCitationsRegister(hypreCitation,&cite);
987:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter));
988:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol));
989:   jac->applyrichardson = PETSC_TRUE;
990:   PCApply_HYPRE(pc,b,y);
991:   jac->applyrichardson = PETSC_FALSE;
992:   PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,&oits));
993:   *outits = oits;
994:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
995:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
996:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
997:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
998:   return(0);
999: }

1001: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
1002: {
1003:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1005:   PetscBool      iascii;

1008:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1009:   if (iascii) {
1010:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG preconditioning\n");
1011:     PetscViewerASCIIPrintf(viewer,"    Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);
1012:     PetscViewerASCIIPrintf(viewer,"    Maximum number of levels %D\n",jac->maxlevels);
1013:     PetscViewerASCIIPrintf(viewer,"    Maximum number of iterations PER hypre call %D\n",jac->maxiter);
1014:     PetscViewerASCIIPrintf(viewer,"    Convergence tolerance PER hypre call %g\n",(double)jac->tol);
1015:     PetscViewerASCIIPrintf(viewer,"    Threshold for strong coupling %g\n",(double)jac->strongthreshold);
1016:     PetscViewerASCIIPrintf(viewer,"    Interpolation truncation factor %g\n",(double)jac->truncfactor);
1017:     PetscViewerASCIIPrintf(viewer,"    Interpolation: max elements per row %D\n",jac->pmax);
1018:     if (jac->interp_refine) {
1019:       PetscViewerASCIIPrintf(viewer,"    Interpolation: number of steps of weighted refinement %D\n",jac->interp_refine);
1020:     }
1021:     PetscViewerASCIIPrintf(viewer,"    Number of levels of aggressive coarsening %D\n",jac->agg_nl);
1022:     PetscViewerASCIIPrintf(viewer,"    Number of paths for aggressive coarsening %D\n",jac->agg_num_paths);

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

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

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

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

1037:     if (jac->relaxorder) {
1038:       PetscViewerASCIIPrintf(viewer,"    Using CF-relaxation\n");
1039:     } else {
1040:       PetscViewerASCIIPrintf(viewer,"    Not using CF-relaxation\n");
1041:     }
1042:     if (jac->smoothtype!=-1) {
1043:       PetscViewerASCIIPrintf(viewer,"    Smooth type          %s\n",HYPREBoomerAMGSmoothType[jac->smoothtype]);
1044:       PetscViewerASCIIPrintf(viewer,"    Smooth num levels    %D\n",jac->smoothnumlevels);
1045:     } else {
1046:       PetscViewerASCIIPrintf(viewer,"    Not using more complex smoothers.\n");
1047:     }
1048:     if (jac->smoothtype==3) {
1049:       PetscViewerASCIIPrintf(viewer,"    Euclid ILU(k) levels %D\n",jac->eu_level);
1050:       PetscViewerASCIIPrintf(viewer,"    Euclid ILU(k) drop tolerance %g\n",(double)jac->eu_droptolerance);
1051:       PetscViewerASCIIPrintf(viewer,"    Euclid ILU use Block-Jacobi? %D\n",jac->eu_bj);
1052:     }
1053:     PetscViewerASCIIPrintf(viewer,"    Measure type        %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
1054:     PetscViewerASCIIPrintf(viewer,"    Coarsen type        %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
1055:     PetscViewerASCIIPrintf(viewer,"    Interpolation type  %s\n",jac->interptype != 100 ? HYPREBoomerAMGInterpType[jac->interptype] : "1pt");
1056:     if (jac->nodal_coarsening) {
1057:       PetscViewerASCIIPrintf(viewer,"    Using nodal coarsening with HYPRE_BOOMERAMGSetNodal() %D\n",jac->nodal_coarsening);
1058:     }
1059:     if (jac->vec_interp_variant) {
1060:       PetscViewerASCIIPrintf(viewer,"    HYPRE_BoomerAMGSetInterpVecVariant() %D\n",jac->vec_interp_variant);
1061:       PetscViewerASCIIPrintf(viewer,"    HYPRE_BoomerAMGSetInterpVecQMax() %D\n",jac->vec_interp_qmax);
1062:       PetscViewerASCIIPrintf(viewer,"    HYPRE_BoomerAMGSetSmoothInterpVectors() %d\n",jac->vec_interp_smooth);
1063:     }
1064:     if (jac->nodal_relax) {
1065:       PetscViewerASCIIPrintf(viewer,"    Using nodal relaxation via Schwarz smoothing on levels %D\n",jac->nodal_relax_levels);
1066:     }

1068:     /* AIR */
1069:     if (jac->Rtype) {
1070:       PetscViewerASCIIPrintf(viewer,"    Using approximate ideal restriction type %D\n",jac->Rtype);
1071:       PetscViewerASCIIPrintf(viewer,"      Threshold for R %g\n",(double)jac->Rstrongthreshold);
1072:       PetscViewerASCIIPrintf(viewer,"      Filter for R %g\n",(double)jac->Rfilterthreshold);
1073:       PetscViewerASCIIPrintf(viewer,"      A drop tolerance %g\n",(double)jac->Adroptol);
1074:       PetscViewerASCIIPrintf(viewer,"      A drop type %D\n",jac->Adroptype);
1075:     }
1076:   }
1077:   return(0);
1078: }

1080: /* --------------------------------------------------------------------------------------------*/
1081: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PetscOptionItems *PetscOptionsObject,PC pc)
1082: {
1083:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1085:   PetscInt       indx;
1086:   PetscBool      flag;
1087:   const char     *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};

1090:   PetscOptionsHead(PetscOptionsObject,"HYPRE ParaSails Options");
1091:   PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);
1092:   PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshold,&jac->threshold,&flag);
1093:   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshold,jac->nlevels));

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

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

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

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

1107:   PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);
1108:   if (flag) {
1109:     jac->symt = indx;
1110:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
1111:   }

1113:   PetscOptionsTail();
1114:   return(0);
1115: }

1117: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
1118: {
1119:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1121:   PetscBool      iascii;
1122:   const char     *symt = 0;

1125:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1126:   if (iascii) {
1127:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails preconditioning\n");
1128:     PetscViewerASCIIPrintf(viewer,"    nlevels %d\n",jac->nlevels);
1129:     PetscViewerASCIIPrintf(viewer,"    threshold %g\n",(double)jac->threshold);
1130:     PetscViewerASCIIPrintf(viewer,"    filter %g\n",(double)jac->filter);
1131:     PetscViewerASCIIPrintf(viewer,"    load balance %g\n",(double)jac->loadbal);
1132:     PetscViewerASCIIPrintf(viewer,"    reuse nonzero structure %s\n",PetscBools[jac->ruse]);
1133:     PetscViewerASCIIPrintf(viewer,"    print info to screen %s\n",PetscBools[jac->logging]);
1134:     if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
1135:     else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
1136:     else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
1137:     else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
1138:     PetscViewerASCIIPrintf(viewer,"    %s\n",symt);
1139:   }
1140:   return(0);
1141: }
1142: /* --------------------------------------------------------------------------------------------*/
1143: static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PetscOptionItems *PetscOptionsObject,PC pc)
1144: {
1145:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1147:   PetscInt       n;
1148:   PetscBool      flag,flag2,flag3,flag4;

1151:   PetscOptionsHead(PetscOptionsObject,"HYPRE AMS Options");
1152:   PetscOptionsInt("-pc_hypre_ams_print_level","Debugging output level for AMS","None",jac->as_print,&jac->as_print,&flag);
1153:   if (flag) PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
1154:   PetscOptionsInt("-pc_hypre_ams_max_iter","Maximum number of AMS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
1155:   if (flag) PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1156:   PetscOptionsInt("-pc_hypre_ams_cycle_type","Cycle type for AMS multigrid","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag);
1157:   if (flag) PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1158:   PetscOptionsReal("-pc_hypre_ams_tol","Error tolerance for AMS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
1159:   if (flag) PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
1160:   PetscOptionsInt("-pc_hypre_ams_relax_type","Relaxation type for AMS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
1161:   PetscOptionsInt("-pc_hypre_ams_relax_times","Number of relaxation steps for AMS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
1162:   PetscOptionsReal("-pc_hypre_ams_relax_weight","Relaxation weight for AMS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
1163:   PetscOptionsReal("-pc_hypre_ams_omega","SSOR coefficient for AMS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
1164:   if (flag || flag2 || flag3 || flag4) {
1165:     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1166:                                                                       jac->as_relax_times,
1167:                                                                       jac->as_relax_weight,
1168:                                                                       jac->as_omega));
1169:   }
1170:   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);
1171:   n = 5;
1172:   PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options","AMG options for vector Poisson","None",jac->as_amg_alpha_opts,&n,&flag2);
1173:   if (flag || flag2) {
1174:     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1175:                                                                      jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1176:                                                                      jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1177:                                                                      jac->as_amg_alpha_theta,
1178:                                                                      jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1179:                                                                      jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1180:   }
1181:   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);
1182:   n = 5;
1183:   PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options","AMG options for scalar Poisson solver","None",jac->as_amg_beta_opts,&n,&flag2);
1184:   if (flag || flag2) {
1185:     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1186:                                                                     jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1187:                                                                     jac->as_amg_beta_opts[2],       /* AMG relax_type */
1188:                                                                     jac->as_amg_beta_theta,
1189:                                                                     jac->as_amg_beta_opts[3],       /* AMG interp_type */
1190:                                                                     jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1191:   }
1192:   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);
1193:   if (flag) { /* override HYPRE's default only if the options is used */
1194:     PetscStackCallStandard(HYPRE_AMSSetProjectionFrequency,(jac->hsolver,jac->ams_proj_freq));
1195:   }
1196:   PetscOptionsTail();
1197:   return(0);
1198: }

1200: static PetscErrorCode PCView_HYPRE_AMS(PC pc,PetscViewer viewer)
1201: {
1202:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1204:   PetscBool      iascii;

1207:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1208:   if (iascii) {
1209:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS preconditioning\n");
1210:     PetscViewerASCIIPrintf(viewer,"    subspace iterations per application %d\n",jac->as_max_iter);
1211:     PetscViewerASCIIPrintf(viewer,"    subspace cycle type %d\n",jac->ams_cycle_type);
1212:     PetscViewerASCIIPrintf(viewer,"    subspace iteration tolerance %g\n",jac->as_tol);
1213:     PetscViewerASCIIPrintf(viewer,"    smoother type %d\n",jac->as_relax_type);
1214:     PetscViewerASCIIPrintf(viewer,"    number of smoothing steps %d\n",jac->as_relax_times);
1215:     PetscViewerASCIIPrintf(viewer,"    smoother weight %g\n",jac->as_relax_weight);
1216:     PetscViewerASCIIPrintf(viewer,"    smoother omega %g\n",jac->as_omega);
1217:     if (jac->alpha_Poisson) {
1218:       PetscViewerASCIIPrintf(viewer,"    vector Poisson solver (passed in by user)\n");
1219:     } else {
1220:       PetscViewerASCIIPrintf(viewer,"    vector Poisson solver (computed) \n");
1221:     }
1222:     PetscViewerASCIIPrintf(viewer,"        boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);
1223:     PetscViewerASCIIPrintf(viewer,"        boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
1224:     PetscViewerASCIIPrintf(viewer,"        boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);
1225:     PetscViewerASCIIPrintf(viewer,"        boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);
1226:     PetscViewerASCIIPrintf(viewer,"        boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
1227:     PetscViewerASCIIPrintf(viewer,"        boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);
1228:     if (!jac->ams_beta_is_zero) {
1229:       if (jac->beta_Poisson) {
1230:         PetscViewerASCIIPrintf(viewer,"    scalar Poisson solver (passed in by user)\n");
1231:       } else {
1232:         PetscViewerASCIIPrintf(viewer,"    scalar Poisson solver (computed) \n");
1233:       }
1234:       PetscViewerASCIIPrintf(viewer,"        boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);
1235:       PetscViewerASCIIPrintf(viewer,"        boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
1236:       PetscViewerASCIIPrintf(viewer,"        boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);
1237:       PetscViewerASCIIPrintf(viewer,"        boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);
1238:       PetscViewerASCIIPrintf(viewer,"        boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
1239:       PetscViewerASCIIPrintf(viewer,"        boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);
1240:       if (jac->ams_beta_is_zero_part) {
1241:         PetscViewerASCIIPrintf(viewer,"        compatible subspace projection frequency %d (-1 HYPRE uses default)\n",jac->ams_proj_freq);
1242:       }
1243:     } else {
1244:       PetscViewerASCIIPrintf(viewer,"    scalar Poisson solver not used (zero-conductivity everywhere) \n");
1245:     }
1246:   }
1247:   return(0);
1248: }

1250: static PetscErrorCode PCSetFromOptions_HYPRE_ADS(PetscOptionItems *PetscOptionsObject,PC pc)
1251: {
1252:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1254:   PetscInt       n;
1255:   PetscBool      flag,flag2,flag3,flag4;

1258:   PetscOptionsHead(PetscOptionsObject,"HYPRE ADS Options");
1259:   PetscOptionsInt("-pc_hypre_ads_print_level","Debugging output level for ADS","None",jac->as_print,&jac->as_print,&flag);
1260:   if (flag) PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
1261:   PetscOptionsInt("-pc_hypre_ads_max_iter","Maximum number of ADS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
1262:   if (flag) PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1263:   PetscOptionsInt("-pc_hypre_ads_cycle_type","Cycle type for ADS multigrid","None",jac->ads_cycle_type,&jac->ads_cycle_type,&flag);
1264:   if (flag) PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ads_cycle_type));
1265:   PetscOptionsReal("-pc_hypre_ads_tol","Error tolerance for ADS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
1266:   if (flag) PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
1267:   PetscOptionsInt("-pc_hypre_ads_relax_type","Relaxation type for ADS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
1268:   PetscOptionsInt("-pc_hypre_ads_relax_times","Number of relaxation steps for ADS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
1269:   PetscOptionsReal("-pc_hypre_ads_relax_weight","Relaxation weight for ADS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
1270:   PetscOptionsReal("-pc_hypre_ads_omega","SSOR coefficient for ADS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
1271:   if (flag || flag2 || flag3 || flag4) {
1272:     PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1273:                                                                       jac->as_relax_times,
1274:                                                                       jac->as_relax_weight,
1275:                                                                       jac->as_omega));
1276:   }
1277:   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);
1278:   n = 5;
1279:   PetscOptionsIntArray("-pc_hypre_ads_ams_options","AMG options for AMS solver inside ADS","None",jac->as_amg_alpha_opts,&n,&flag2);
1280:   PetscOptionsInt("-pc_hypre_ads_ams_cycle_type","Cycle type for AMS solver inside ADS","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag3);
1281:   if (flag || flag2 || flag3) {
1282:     PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type,             /* AMS cycle type */
1283:                                                                 jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1284:                                                                 jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1285:                                                                 jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1286:                                                                 jac->as_amg_alpha_theta,
1287:                                                                 jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1288:                                                                 jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1289:   }
1290:   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);
1291:   n = 5;
1292:   PetscOptionsIntArray("-pc_hypre_ads_amg_options","AMG options for vector AMG solver inside ADS","None",jac->as_amg_beta_opts,&n,&flag2);
1293:   if (flag || flag2) {
1294:     PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1295:                                                                 jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1296:                                                                 jac->as_amg_beta_opts[2],       /* AMG relax_type */
1297:                                                                 jac->as_amg_beta_theta,
1298:                                                                 jac->as_amg_beta_opts[3],       /* AMG interp_type */
1299:                                                                 jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1300:   }
1301:   PetscOptionsTail();
1302:   return(0);
1303: }

1305: static PetscErrorCode PCView_HYPRE_ADS(PC pc,PetscViewer viewer)
1306: {
1307:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1309:   PetscBool      iascii;

1312:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1313:   if (iascii) {
1314:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS preconditioning\n");
1315:     PetscViewerASCIIPrintf(viewer,"    subspace iterations per application %d\n",jac->as_max_iter);
1316:     PetscViewerASCIIPrintf(viewer,"    subspace cycle type %d\n",jac->ads_cycle_type);
1317:     PetscViewerASCIIPrintf(viewer,"    subspace iteration tolerance %g\n",jac->as_tol);
1318:     PetscViewerASCIIPrintf(viewer,"    smoother type %d\n",jac->as_relax_type);
1319:     PetscViewerASCIIPrintf(viewer,"    number of smoothing steps %d\n",jac->as_relax_times);
1320:     PetscViewerASCIIPrintf(viewer,"    smoother weight %g\n",jac->as_relax_weight);
1321:     PetscViewerASCIIPrintf(viewer,"    smoother omega %g\n",jac->as_omega);
1322:     PetscViewerASCIIPrintf(viewer,"    AMS solver using boomerAMG\n");
1323:     PetscViewerASCIIPrintf(viewer,"        subspace cycle type %d\n",jac->ams_cycle_type);
1324:     PetscViewerASCIIPrintf(viewer,"        coarsening type %d\n",jac->as_amg_alpha_opts[0]);
1325:     PetscViewerASCIIPrintf(viewer,"        levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
1326:     PetscViewerASCIIPrintf(viewer,"        relaxation type %d\n",jac->as_amg_alpha_opts[2]);
1327:     PetscViewerASCIIPrintf(viewer,"        interpolation type %d\n",jac->as_amg_alpha_opts[3]);
1328:     PetscViewerASCIIPrintf(viewer,"        max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
1329:     PetscViewerASCIIPrintf(viewer,"        strength threshold %g\n",jac->as_amg_alpha_theta);
1330:     PetscViewerASCIIPrintf(viewer,"    vector Poisson solver using boomerAMG\n");
1331:     PetscViewerASCIIPrintf(viewer,"        coarsening type %d\n",jac->as_amg_beta_opts[0]);
1332:     PetscViewerASCIIPrintf(viewer,"        levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
1333:     PetscViewerASCIIPrintf(viewer,"        relaxation type %d\n",jac->as_amg_beta_opts[2]);
1334:     PetscViewerASCIIPrintf(viewer,"        interpolation type %d\n",jac->as_amg_beta_opts[3]);
1335:     PetscViewerASCIIPrintf(viewer,"        max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
1336:     PetscViewerASCIIPrintf(viewer,"        strength threshold %g\n",jac->as_amg_beta_theta);
1337:   }
1338:   return(0);
1339: }

1341: static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE(PC pc, Mat G)
1342: {
1343:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1344:   PetscBool      ishypre;

1348:   PetscObjectTypeCompare((PetscObject)G,MATHYPRE,&ishypre);
1349:   if (ishypre) {
1350:     PetscObjectReference((PetscObject)G);
1351:     MatDestroy(&jac->G);
1352:     jac->G = G;
1353:   } else {
1354:     MatDestroy(&jac->G);
1355:     MatConvert(G,MATHYPRE,MAT_INITIAL_MATRIX,&jac->G);
1356:   }
1357:   return(0);
1358: }

1360: /*@
1361:  PCHYPRESetDiscreteGradient - Set discrete gradient matrix

1363:    Collective on PC

1365:    Input Parameters:
1366: +  pc - the preconditioning context
1367: -  G - the discrete gradient

1369:    Level: intermediate

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

1375: .seealso:
1376: @*/
1377: PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G)
1378: {

1385:   PetscTryMethod(pc,"PCHYPRESetDiscreteGradient_C",(PC,Mat),(pc,G));
1386:   return(0);
1387: }

1389: static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C)
1390: {
1391:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1392:   PetscBool      ishypre;

1396:   PetscObjectTypeCompare((PetscObject)C,MATHYPRE,&ishypre);
1397:   if (ishypre) {
1398:     PetscObjectReference((PetscObject)C);
1399:     MatDestroy(&jac->C);
1400:     jac->C = C;
1401:   } else {
1402:     MatDestroy(&jac->C);
1403:     MatConvert(C,MATHYPRE,MAT_INITIAL_MATRIX,&jac->C);
1404:   }
1405:   return(0);
1406: }

1408: /*@
1409:  PCHYPRESetDiscreteCurl - Set discrete curl matrix

1411:    Collective on PC

1413:    Input Parameters:
1414: +  pc - the preconditioning context
1415: -  C - the discrete curl

1417:    Level: intermediate

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

1423: .seealso:
1424: @*/
1425: PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C)
1426: {

1433:   PetscTryMethod(pc,"PCHYPRESetDiscreteCurl_C",(PC,Mat),(pc,C));
1434:   return(0);
1435: }

1437: static PetscErrorCode PCHYPRESetInterpolations_HYPRE(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[])
1438: {
1439:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1440:   PetscBool      ishypre;
1442:   PetscInt       i;

1445:   MatDestroy(&jac->RT_PiFull);
1446:   MatDestroy(&jac->ND_PiFull);
1447:   for (i=0;i<3;++i) {
1448:     MatDestroy(&jac->RT_Pi[i]);
1449:     MatDestroy(&jac->ND_Pi[i]);
1450:   }

1452:   jac->dim = dim;
1453:   if (RT_PiFull) {
1454:     PetscObjectTypeCompare((PetscObject)RT_PiFull,MATHYPRE,&ishypre);
1455:     if (ishypre) {
1456:       PetscObjectReference((PetscObject)RT_PiFull);
1457:       jac->RT_PiFull = RT_PiFull;
1458:     } else {
1459:       MatConvert(RT_PiFull,MATHYPRE,MAT_INITIAL_MATRIX,&jac->RT_PiFull);
1460:     }
1461:   }
1462:   if (RT_Pi) {
1463:     for (i=0;i<dim;++i) {
1464:       if (RT_Pi[i]) {
1465:         PetscObjectTypeCompare((PetscObject)RT_Pi[i],MATHYPRE,&ishypre);
1466:         if (ishypre) {
1467:           PetscObjectReference((PetscObject)RT_Pi[i]);
1468:           jac->RT_Pi[i] = RT_Pi[i];
1469:         } else {
1470:           MatConvert(RT_Pi[i],MATHYPRE,MAT_INITIAL_MATRIX,&jac->RT_Pi[i]);
1471:         }
1472:       }
1473:     }
1474:   }
1475:   if (ND_PiFull) {
1476:     PetscObjectTypeCompare((PetscObject)ND_PiFull,MATHYPRE,&ishypre);
1477:     if (ishypre) {
1478:       PetscObjectReference((PetscObject)ND_PiFull);
1479:       jac->ND_PiFull = ND_PiFull;
1480:     } else {
1481:       MatConvert(ND_PiFull,MATHYPRE,MAT_INITIAL_MATRIX,&jac->ND_PiFull);
1482:     }
1483:   }
1484:   if (ND_Pi) {
1485:     for (i=0;i<dim;++i) {
1486:       if (ND_Pi[i]) {
1487:         PetscObjectTypeCompare((PetscObject)ND_Pi[i],MATHYPRE,&ishypre);
1488:         if (ishypre) {
1489:           PetscObjectReference((PetscObject)ND_Pi[i]);
1490:           jac->ND_Pi[i] = ND_Pi[i];
1491:         } else {
1492:           MatConvert(ND_Pi[i],MATHYPRE,MAT_INITIAL_MATRIX,&jac->ND_Pi[i]);
1493:         }
1494:       }
1495:     }
1496:   }

1498:   return(0);
1499: }

1501: /*@
1502:  PCHYPRESetInterpolations - Set interpolation matrices for AMS/ADS preconditioner

1504:    Collective on PC

1506:    Input Parameters:
1507: +  pc - the preconditioning context
1508: -  dim - the dimension of the problem, only used in AMS
1509: -  RT_PiFull - Raviart-Thomas interpolation matrix
1510: -  RT_Pi - x/y/z component of Raviart-Thomas interpolation matrix
1511: -  ND_PiFull - Nedelec interpolation matrix
1512: -  ND_Pi - x/y/z component of Nedelec interpolation matrix

1514:    Notes:
1515:     For AMS, only Nedelec interpolation matrices are needed, the Raviart-Thomas interpolation matrices can be set to NULL.
1516:           For ADS, both type of interpolation matrices are needed.
1517:    Level: intermediate

1519: .seealso:
1520: @*/
1521: PetscErrorCode PCHYPRESetInterpolations(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[])
1522: {
1524:   PetscInt       i;

1528:   if (RT_PiFull) {
1531:   }
1532:   if (RT_Pi) {
1534:     for (i=0;i<dim;++i) {
1535:       if (RT_Pi[i]) {
1538:       }
1539:     }
1540:   }
1541:   if (ND_PiFull) {
1544:   }
1545:   if (ND_Pi) {
1547:     for (i=0;i<dim;++i) {
1548:       if (ND_Pi[i]) {
1551:       }
1552:     }
1553:   }
1554:   PetscTryMethod(pc,"PCHYPRESetInterpolations_C",(PC,PetscInt,Mat,Mat[],Mat,Mat[]),(pc,dim,RT_PiFull,RT_Pi,ND_PiFull,ND_Pi));
1555:   return(0);
1556: }

1558: static PetscErrorCode PCHYPRESetPoissonMatrix_HYPRE(PC pc, Mat A, PetscBool isalpha)
1559: {
1560:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1561:   PetscBool      ishypre;

1565:   PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
1566:   if (ishypre) {
1567:     if (isalpha) {
1568:       PetscObjectReference((PetscObject)A);
1569:       MatDestroy(&jac->alpha_Poisson);
1570:       jac->alpha_Poisson = A;
1571:     } else {
1572:       if (A) {
1573:         PetscObjectReference((PetscObject)A);
1574:       } else {
1575:         jac->ams_beta_is_zero = PETSC_TRUE;
1576:       }
1577:       MatDestroy(&jac->beta_Poisson);
1578:       jac->beta_Poisson = A;
1579:     }
1580:   } else {
1581:     if (isalpha) {
1582:       MatDestroy(&jac->alpha_Poisson);
1583:       MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&jac->alpha_Poisson);
1584:     } else {
1585:       if (A) {
1586:         MatDestroy(&jac->beta_Poisson);
1587:         MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&jac->beta_Poisson);
1588:       } else {
1589:         MatDestroy(&jac->beta_Poisson);
1590:         jac->ams_beta_is_zero = PETSC_TRUE;
1591:       }
1592:     }
1593:   }
1594:   return(0);
1595: }

1597: /*@
1598:  PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix

1600:    Collective on PC

1602:    Input Parameters:
1603: +  pc - the preconditioning context
1604: -  A - the matrix

1606:    Level: intermediate

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

1611: .seealso:
1612: @*/
1613: PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A)
1614: {

1621:   PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_TRUE));
1622:   return(0);
1623: }

1625: /*@
1626:  PCHYPRESetBetaPoissonMatrix - Set Poisson matrix

1628:    Collective on PC

1630:    Input Parameters:
1631: +  pc - the preconditioning context
1632: -  A - the matrix

1634:    Level: intermediate

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

1640: .seealso:
1641: @*/
1642: PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A)
1643: {

1648:   if (A) {
1651:   }
1652:   PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_FALSE));
1653:   return(0);
1654: }

1656: static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE(PC pc,Vec ozz, Vec zoz, Vec zzo)
1657: {
1658:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1659:   PetscErrorCode     ierr;

1662:   /* throw away any vector if already set */
1663:   VecHYPRE_IJVectorDestroy(&jac->constants[0]);
1664:   VecHYPRE_IJVectorDestroy(&jac->constants[1]);
1665:   VecHYPRE_IJVectorDestroy(&jac->constants[2]);
1666:   VecHYPRE_IJVectorCreate(ozz->map,&jac->constants[0]);
1667:   VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);
1668:   VecHYPRE_IJVectorCreate(zoz->map,&jac->constants[1]);
1669:   VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);
1670:   jac->dim = 2;
1671:   if (zzo) {
1672:     VecHYPRE_IJVectorCreate(zzo->map,&jac->constants[2]);
1673:     VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);
1674:     jac->dim++;
1675:   }
1676:   return(0);
1677: }

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

1682:    Collective on PC

1684:    Input Parameters:
1685: +  pc - the preconditioning context
1686: -  ozz - vector representing (1,0,0) (or (1,0) in 2D)
1687: -  zoz - vector representing (0,1,0) (or (0,1) in 2D)
1688: -  zzo - vector representing (0,0,1) (use NULL in 2D)

1690:    Level: intermediate

1692:    Notes:

1694: .seealso:
1695: @*/
1696: PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo)
1697: {

1708:   PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));
1709:   return(0);
1710: }

1712: static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1713: {
1714:   PC_HYPRE        *jac = (PC_HYPRE*)pc->data;
1715:   Vec             tv;
1716:   PetscInt        i;
1717:   PetscErrorCode  ierr;

1720:   /* throw away any coordinate vector if already set */
1721:   VecHYPRE_IJVectorDestroy(&jac->coords[0]);
1722:   VecHYPRE_IJVectorDestroy(&jac->coords[1]);
1723:   VecHYPRE_IJVectorDestroy(&jac->coords[2]);
1724:   jac->dim = dim;

1726:   /* compute IJ vector for coordinates */
1727:   VecCreate(PetscObjectComm((PetscObject)pc),&tv);
1728:   VecSetType(tv,VECSTANDARD);
1729:   VecSetSizes(tv,nloc,PETSC_DECIDE);
1730:   for (i=0;i<dim;i++) {
1731:     PetscScalar *array;
1732:     PetscInt    j;

1734:     VecHYPRE_IJVectorCreate(tv->map,&jac->coords[i]);
1735:     VecGetArrayWrite(tv,&array);
1736:     for (j=0;j<nloc;j++) array[j] = coords[j*dim+i];
1737:     VecRestoreArrayWrite(tv,&array);
1738:     VecHYPRE_IJVectorCopy(tv,jac->coords[i]);
1739:   }
1740:   VecDestroy(&tv);
1741:   return(0);
1742: }

1744: /* ---------------------------------------------------------------------------------*/

1746: static PetscErrorCode  PCHYPREGetType_HYPRE(PC pc,const char *name[])
1747: {
1748:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

1751:   *name = jac->hypre_type;
1752:   return(0);
1753: }

1755: static PetscErrorCode  PCHYPRESetType_HYPRE(PC pc,const char name[])
1756: {
1757:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1759:   PetscBool      flag;

1762:   if (jac->hypre_type) {
1763:     PetscStrcmp(jac->hypre_type,name,&flag);
1764:     if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set");
1765:     return(0);
1766:   } else {
1767:     PetscStrallocpy(name, &jac->hypre_type);
1768:   }

1770:   jac->maxiter         = PETSC_DEFAULT;
1771:   jac->tol             = PETSC_DEFAULT;
1772:   jac->printstatistics = PetscLogPrintInfo;

1774:   PetscStrcmp("pilut",jac->hypre_type,&flag);
1775:   if (flag) {
1776:     MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1777:     PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver));
1778:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
1779:     pc->ops->view           = PCView_HYPRE_Pilut;
1780:     jac->destroy            = HYPRE_ParCSRPilutDestroy;
1781:     jac->setup              = HYPRE_ParCSRPilutSetup;
1782:     jac->solve              = HYPRE_ParCSRPilutSolve;
1783:     jac->factorrowsize      = PETSC_DEFAULT;
1784:     return(0);
1785:   }
1786:   PetscStrcmp("euclid",jac->hypre_type,&flag);
1787:   if (flag) {
1788: #if defined(PETSC_HAVE_64BIT_INDICES)
1789:     SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Hypre Euclid not support with 64 bit indices");
1790: #endif
1791:     MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1792:     PetscStackCallStandard(HYPRE_EuclidCreate,(jac->comm_hypre,&jac->hsolver));
1793:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
1794:     pc->ops->view           = PCView_HYPRE_Euclid;
1795:     jac->destroy            = HYPRE_EuclidDestroy;
1796:     jac->setup              = HYPRE_EuclidSetup;
1797:     jac->solve              = HYPRE_EuclidSolve;
1798:     jac->factorrowsize      = PETSC_DEFAULT;
1799:     jac->eu_level           = PETSC_DEFAULT; /* default */
1800:     return(0);
1801:   }
1802:   PetscStrcmp("parasails",jac->hypre_type,&flag);
1803:   if (flag) {
1804:     MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1805:     PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver));
1806:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
1807:     pc->ops->view           = PCView_HYPRE_ParaSails;
1808:     jac->destroy            = HYPRE_ParaSailsDestroy;
1809:     jac->setup              = HYPRE_ParaSailsSetup;
1810:     jac->solve              = HYPRE_ParaSailsSolve;
1811:     /* initialize */
1812:     jac->nlevels   = 1;
1813:     jac->threshold = .1;
1814:     jac->filter    = .1;
1815:     jac->loadbal   = 0;
1816:     if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE;
1817:     else jac->logging = (int) PETSC_FALSE;

1819:     jac->ruse = (int) PETSC_FALSE;
1820:     jac->symt = 0;
1821:     PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshold,jac->nlevels));
1822:     PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
1823:     PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
1824:     PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
1825:     PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
1826:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
1827:     return(0);
1828:   }
1829:   PetscStrcmp("boomeramg",jac->hypre_type,&flag);
1830:   if (flag) {
1831:     HYPRE_BoomerAMGCreate(&jac->hsolver);
1832:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_BoomerAMG;
1833:     pc->ops->view            = PCView_HYPRE_BoomerAMG;
1834:     pc->ops->applytranspose  = PCApplyTranspose_HYPRE_BoomerAMG;
1835:     pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
1836:     PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",PCGetInterpolations_BoomerAMG);
1837:     PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",PCGetCoarseOperators_BoomerAMG);
1838:     jac->destroy             = HYPRE_BoomerAMGDestroy;
1839:     jac->setup               = HYPRE_BoomerAMGSetup;
1840:     jac->solve               = HYPRE_BoomerAMGSolve;
1841:     jac->applyrichardson     = PETSC_FALSE;
1842:     /* these defaults match the hypre defaults */
1843:     jac->cycletype        = 1;
1844:     jac->maxlevels        = 25;
1845:     jac->maxiter          = 1;
1846:     jac->tol              = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
1847:     jac->truncfactor      = 0.0;
1848:     jac->strongthreshold  = .25;
1849:     jac->maxrowsum        = .9;
1850:     jac->coarsentype      = 6;
1851:     jac->measuretype      = 0;
1852:     jac->gridsweeps[0]    = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
1853:     jac->smoothtype       = -1; /* Not set by default */
1854:     jac->smoothnumlevels  = 25;
1855:     jac->eu_level         = 0;
1856:     jac->eu_droptolerance = 0;
1857:     jac->eu_bj            = 0;
1858:     jac->relaxtype[0]     = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a PC - most likely with CG */
1859:     jac->relaxtype[2]     = 9; /*G.E. */
1860:     jac->relaxweight      = 1.0;
1861:     jac->outerrelaxweight = 1.0;
1862:     jac->relaxorder       = 1;
1863:     jac->interptype       = 0;
1864:     jac->Rtype            = 0;
1865:     jac->Rstrongthreshold = 0.25;
1866:     jac->Rfilterthreshold = 0.0;
1867:     jac->Adroptype        = -1;
1868:     jac->Adroptol         = 0.0;
1869:     jac->agg_nl           = 0;
1870:     jac->agg_interptype   = 4;
1871:     jac->pmax             = 0;
1872:     jac->truncfactor      = 0.0;
1873:     jac->agg_num_paths    = 1;
1874:     jac->maxc             = 9;
1875:     jac->minc             = 1;

1877:     jac->nodal_coarsening      = 0;
1878:     jac->nodal_coarsening_diag = 0;
1879:     jac->vec_interp_variant    = 0;
1880:     jac->vec_interp_qmax       = 0;
1881:     jac->vec_interp_smooth     = PETSC_FALSE;
1882:     jac->interp_refine         = 0;
1883:     jac->nodal_relax           = PETSC_FALSE;
1884:     jac->nodal_relax_levels    = 1;
1885:     jac->rap2                  = 0;

1887:     /* GPU defaults
1888:          from https://hypre.readthedocs.io/en/latest/solvers-boomeramg.html#gpu-supported-options
1889:          and /src/parcsr_ls/par_amg.c */
1890: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1891:     jac->keeptranspose         = PETSC_TRUE;
1892:     jac->mod_rap2              = 1;
1893:     jac->coarsentype           = 8;
1894:     jac->relaxorder            = 0;
1895:     jac->interptype            = 6;
1896:     jac->relaxtype[0]          = 18;
1897:     jac->relaxtype[1]          = 18;
1898:     jac->agg_interptype        = 7;
1899: #else
1900:     jac->keeptranspose         = PETSC_FALSE;
1901:     jac->mod_rap2              = 0;
1902: #endif
1903:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
1904:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
1905:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
1906:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
1907:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
1908:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
1909:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
1910:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
1911:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
1912:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
1913:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
1914:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
1915:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggInterpType,(jac->hsolver,jac->agg_interptype));
1916:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
1917:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
1918:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0]));  /* defaults coarse to 9 */
1919:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /* defaults coarse to 1 */
1920:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxCoarseSize,(jac->hsolver, jac->maxc));
1921:     PetscStackCallStandard(HYPRE_BoomerAMGSetMinCoarseSize,(jac->hsolver, jac->minc));

1923:     /* GPU */
1924: #if PETSC_PKG_HYPRE_VERSION_GE(2,18,0)
1925:     PetscStackCallStandard(HYPRE_BoomerAMGSetKeepTranspose,(jac->hsolver,jac->keeptranspose ? 1 : 0));
1926:     PetscStackCallStandard(HYPRE_BoomerAMGSetRAP2,(jac->hsolver, jac->rap2));
1927:     PetscStackCallStandard(HYPRE_BoomerAMGSetModuleRAP2,(jac->hsolver, jac->mod_rap2));
1928: #endif

1930:     /* AIR */
1931: #if PETSC_PKG_HYPRE_VERSION_GE(2,18,0)
1932:     PetscStackCallStandard(HYPRE_BoomerAMGSetRestriction,(jac->hsolver,jac->Rtype));
1933:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThresholdR,(jac->hsolver,jac->Rstrongthreshold));
1934:     PetscStackCallStandard(HYPRE_BoomerAMGSetFilterThresholdR,(jac->hsolver,jac->Rfilterthreshold));
1935:     PetscStackCallStandard(HYPRE_BoomerAMGSetADropTol,(jac->hsolver,jac->Adroptol));
1936:     PetscStackCallStandard(HYPRE_BoomerAMGSetADropType,(jac->hsolver,jac->Adroptype));
1937: #endif
1938:     return(0);
1939:   }
1940:   PetscStrcmp("ams",jac->hypre_type,&flag);
1941:   if (flag) {
1942:     HYPRE_AMSCreate(&jac->hsolver);
1943:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_AMS;
1944:     pc->ops->view            = PCView_HYPRE_AMS;
1945:     jac->destroy             = HYPRE_AMSDestroy;
1946:     jac->setup               = HYPRE_AMSSetup;
1947:     jac->solve               = HYPRE_AMSSolve;
1948:     jac->coords[0]           = NULL;
1949:     jac->coords[1]           = NULL;
1950:     jac->coords[2]           = NULL;
1951:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */
1952:     jac->as_print           = 0;
1953:     jac->as_max_iter        = 1; /* used as a preconditioner */
1954:     jac->as_tol             = 0.; /* used as a preconditioner */
1955:     jac->ams_cycle_type     = 13;
1956:     /* Smoothing options */
1957:     jac->as_relax_type      = 2;
1958:     jac->as_relax_times     = 1;
1959:     jac->as_relax_weight    = 1.0;
1960:     jac->as_omega           = 1.0;
1961:     /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1962:     jac->as_amg_alpha_opts[0] = 10;
1963:     jac->as_amg_alpha_opts[1] = 1;
1964:     jac->as_amg_alpha_opts[2] = 6;
1965:     jac->as_amg_alpha_opts[3] = 6;
1966:     jac->as_amg_alpha_opts[4] = 4;
1967:     jac->as_amg_alpha_theta   = 0.25;
1968:     /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1969:     jac->as_amg_beta_opts[0] = 10;
1970:     jac->as_amg_beta_opts[1] = 1;
1971:     jac->as_amg_beta_opts[2] = 6;
1972:     jac->as_amg_beta_opts[3] = 6;
1973:     jac->as_amg_beta_opts[4] = 4;
1974:     jac->as_amg_beta_theta   = 0.25;
1975:     PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
1976:     PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1977:     PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1978:     PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
1979:     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1980:                                                                       jac->as_relax_times,
1981:                                                                       jac->as_relax_weight,
1982:                                                                       jac->as_omega));
1983:     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1984:                                                                      jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1985:                                                                      jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1986:                                                                      jac->as_amg_alpha_theta,
1987:                                                                      jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1988:                                                                      jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1989:     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1990:                                                                     jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1991:                                                                     jac->as_amg_beta_opts[2],       /* AMG relax_type */
1992:                                                                     jac->as_amg_beta_theta,
1993:                                                                     jac->as_amg_beta_opts[3],       /* AMG interp_type */
1994:                                                                     jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1995:     /* Zero conductivity */
1996:     jac->ams_beta_is_zero      = PETSC_FALSE;
1997:     jac->ams_beta_is_zero_part = PETSC_FALSE;
1998:     return(0);
1999:   }
2000:   PetscStrcmp("ads",jac->hypre_type,&flag);
2001:   if (flag) {
2002:     HYPRE_ADSCreate(&jac->hsolver);
2003:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_ADS;
2004:     pc->ops->view            = PCView_HYPRE_ADS;
2005:     jac->destroy             = HYPRE_ADSDestroy;
2006:     jac->setup               = HYPRE_ADSSetup;
2007:     jac->solve               = HYPRE_ADSSolve;
2008:     jac->coords[0]           = NULL;
2009:     jac->coords[1]           = NULL;
2010:     jac->coords[2]           = NULL;
2011:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */
2012:     jac->as_print           = 0;
2013:     jac->as_max_iter        = 1; /* used as a preconditioner */
2014:     jac->as_tol             = 0.; /* used as a preconditioner */
2015:     jac->ads_cycle_type     = 13;
2016:     /* Smoothing options */
2017:     jac->as_relax_type      = 2;
2018:     jac->as_relax_times     = 1;
2019:     jac->as_relax_weight    = 1.0;
2020:     jac->as_omega           = 1.0;
2021:     /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */
2022:     jac->ams_cycle_type       = 14;
2023:     jac->as_amg_alpha_opts[0] = 10;
2024:     jac->as_amg_alpha_opts[1] = 1;
2025:     jac->as_amg_alpha_opts[2] = 6;
2026:     jac->as_amg_alpha_opts[3] = 6;
2027:     jac->as_amg_alpha_opts[4] = 4;
2028:     jac->as_amg_alpha_theta   = 0.25;
2029:     /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
2030:     jac->as_amg_beta_opts[0] = 10;
2031:     jac->as_amg_beta_opts[1] = 1;
2032:     jac->as_amg_beta_opts[2] = 6;
2033:     jac->as_amg_beta_opts[3] = 6;
2034:     jac->as_amg_beta_opts[4] = 4;
2035:     jac->as_amg_beta_theta   = 0.25;
2036:     PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
2037:     PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
2038:     PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
2039:     PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
2040:     PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
2041:                                                                       jac->as_relax_times,
2042:                                                                       jac->as_relax_weight,
2043:                                                                       jac->as_omega));
2044:     PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type,             /* AMG coarsen type */
2045:                                                                 jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
2046:                                                                 jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
2047:                                                                 jac->as_amg_alpha_opts[2],       /* AMG relax_type */
2048:                                                                 jac->as_amg_alpha_theta,
2049:                                                                 jac->as_amg_alpha_opts[3],       /* AMG interp_type */
2050:                                                                 jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
2051:     PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
2052:                                                                 jac->as_amg_beta_opts[1],       /* AMG agg_levels */
2053:                                                                 jac->as_amg_beta_opts[2],       /* AMG relax_type */
2054:                                                                 jac->as_amg_beta_theta,
2055:                                                                 jac->as_amg_beta_opts[3],       /* AMG interp_type */
2056:                                                                 jac->as_amg_beta_opts[4]));     /* AMG Pmax */
2057:     return(0);
2058:   }
2059:   PetscFree(jac->hypre_type);

2061:   jac->hypre_type = NULL;
2062:   SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are euclid, pilut, parasails, boomeramg, ams",name);
2063: }

2065: /*
2066:     It only gets here if the HYPRE type has not been set before the call to
2067:    ...SetFromOptions() which actually is most of the time
2068: */
2069: PetscErrorCode PCSetFromOptions_HYPRE(PetscOptionItems *PetscOptionsObject,PC pc)
2070: {
2072:   PetscInt       indx;
2073:   const char     *type[] = {"euclid","pilut","parasails","boomeramg","ams","ads"};
2074:   PetscBool      flg;

2077:   PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");
2078:   PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,ALEN(type),"boomeramg",&indx,&flg);
2079:   if (flg) {
2080:     PCHYPRESetType_HYPRE(pc,type[indx]);
2081:   } else {
2082:     PCHYPRESetType_HYPRE(pc,"boomeramg");
2083:   }
2084:   if (pc->ops->setfromoptions) {
2085:     pc->ops->setfromoptions(PetscOptionsObject,pc);
2086:   }
2087:   PetscOptionsTail();
2088:   return(0);
2089: }

2091: /*@C
2092:      PCHYPRESetType - Sets which hypre preconditioner you wish to use

2094:    Input Parameters:
2095: +     pc - the preconditioner context
2096: -     name - either  euclid, pilut, parasails, boomeramg, ams, ads

2098:    Options Database Keys:
2099:    -pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads

2101:    Level: intermediate

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

2106: @*/
2107: PetscErrorCode  PCHYPRESetType(PC pc,const char name[])
2108: {

2114:   PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));
2115:   return(0);
2116: }

2118: /*@C
2119:      PCHYPREGetType - Gets which hypre preconditioner you are using

2121:    Input Parameter:
2122: .     pc - the preconditioner context

2124:    Output Parameter:
2125: .     name - either  euclid, pilut, parasails, boomeramg, ams, ads

2127:    Level: intermediate

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

2132: @*/
2133: PetscErrorCode  PCHYPREGetType(PC pc,const char *name[])
2134: {

2140:   PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));
2141:   return(0);
2142: }

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

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

2152:    Level: intermediate

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

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

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

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

2177:           MatSetNearNullSpace() - if you provide a near null space to your matrix it is ignored by hypre UNLESS you also use
2178:           the following two options:

2180:    Options Database Keys:
2181: +   -pc_hypre_boomeramg_nodal_coarsen <n> - where n is from 1 to 6 (see HYPRE_BOOMERAMGSetNodal())
2182: -   -pc_hypre_boomeramg_vec_interp_variant <v> - where v is from 1 to 3 (see HYPRE_BoomerAMGSetInterpVecVariant())

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

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

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

2191: M*/

2193: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
2194: {
2195:   PC_HYPRE       *jac;

2199:   PetscNewLog(pc,&jac);

2201:   pc->data                = jac;
2202:   pc->ops->reset          = PCReset_HYPRE;
2203:   pc->ops->destroy        = PCDestroy_HYPRE;
2204:   pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
2205:   pc->ops->setup          = PCSetUp_HYPRE;
2206:   pc->ops->apply          = PCApply_HYPRE;
2207:   jac->comm_hypre         = MPI_COMM_NULL;
2208:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);
2209:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);
2210:   PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);
2211:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);
2212:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",PCHYPRESetDiscreteCurl_HYPRE);
2213:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",PCHYPRESetInterpolations_HYPRE);
2214:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE);
2215:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",PCHYPRESetPoissonMatrix_HYPRE);
2216: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2217: #if defined(HYPRE_USING_HIP)
2218:   PetscHIPInitializeCheck();
2219: #endif
2220: #if defined(HYPRE_USING_CUDA)
2221:   PetscCUDAInitializeCheck();
2222: #endif
2223: #endif
2224:   return(0);
2225: }

2227: /* ---------------------------------------------------------------------------------------------------------------------------------*/

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

2233:   /* keep copy of PFMG options used so may view them */
2234:   PetscInt its;
2235:   double   tol;
2236:   PetscInt relax_type;
2237:   PetscInt rap_type;
2238:   PetscInt num_pre_relax,num_post_relax;
2239:   PetscInt max_levels;
2240: } PC_PFMG;

2242: PetscErrorCode PCDestroy_PFMG(PC pc)
2243: {
2245:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

2248:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
2249:   MPI_Comm_free(&ex->hcomm);
2250:   PetscFree(pc->data);
2251:   return(0);
2252: }

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

2257: PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
2258: {
2260:   PetscBool      iascii;
2261:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

2264:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2265:   if (iascii) {
2266:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG preconditioning\n");
2267:     PetscViewerASCIIPrintf(viewer,"    max iterations %d\n",ex->its);
2268:     PetscViewerASCIIPrintf(viewer,"    tolerance %g\n",ex->tol);
2269:     PetscViewerASCIIPrintf(viewer,"    relax type %s\n",PFMGRelaxType[ex->relax_type]);
2270:     PetscViewerASCIIPrintf(viewer,"    RAP type %s\n",PFMGRAPType[ex->rap_type]);
2271:     PetscViewerASCIIPrintf(viewer,"    number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
2272:     PetscViewerASCIIPrintf(viewer,"    max levels %d\n",ex->max_levels);
2273:   }
2274:   return(0);
2275: }

2277: PetscErrorCode PCSetFromOptions_PFMG(PetscOptionItems *PetscOptionsObject,PC pc)
2278: {
2280:   PC_PFMG        *ex = (PC_PFMG*) pc->data;
2281:   PetscBool      flg = PETSC_FALSE;

2284:   PetscOptionsHead(PetscOptionsObject,"PFMG options");
2285:   PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);
2286:   if (flg) {
2287:     PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,3));
2288:   }
2289:   PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);
2290:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its));
2291:   PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
2292:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax));
2293:   PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
2294:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax));

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

2299:   PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);
2300:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol));
2301:   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);
2302:   PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type));
2303:   PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);
2304:   PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type));
2305:   PetscOptionsTail();
2306:   return(0);
2307: }

2309: PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
2310: {
2311:   PetscErrorCode    ierr;
2312:   PC_PFMG           *ex = (PC_PFMG*) pc->data;
2313:   PetscScalar       *yy;
2314:   const PetscScalar *xx;
2315:   PetscInt          ilower[3],iupper[3];
2316:   HYPRE_Int         hlower[3],hupper[3];
2317:   Mat_HYPREStruct   *mx = (Mat_HYPREStruct*)(pc->pmat->data);

2320:   PetscCitationsRegister(hypreCitation,&cite);
2321:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
2322:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2323:   iupper[0] += ilower[0] - 1;
2324:   iupper[1] += ilower[1] - 1;
2325:   iupper[2] += ilower[2] - 1;
2326:   hlower[0]  = (HYPRE_Int)ilower[0];
2327:   hlower[1]  = (HYPRE_Int)ilower[1];
2328:   hlower[2]  = (HYPRE_Int)ilower[2];
2329:   hupper[0]  = (HYPRE_Int)iupper[0];
2330:   hupper[1]  = (HYPRE_Int)iupper[1];
2331:   hupper[2]  = (HYPRE_Int)iupper[2];

2333:   /* copy x values over to hypre */
2334:   PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
2335:   VecGetArrayRead(x,&xx);
2336:   PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,hlower,hupper,(HYPRE_Complex*)xx));
2337:   VecRestoreArrayRead(x,&xx);
2338:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
2339:   PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx));

2341:   /* copy solution values back to PETSc */
2342:   VecGetArray(y,&yy);
2343:   PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,hlower,hupper,(HYPRE_Complex*)yy));
2344:   VecRestoreArray(y,&yy);
2345:   return(0);
2346: }

2348: 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)
2349: {
2350:   PC_PFMG        *jac = (PC_PFMG*)pc->data;
2352:   HYPRE_Int      oits;

2355:   PetscCitationsRegister(hypreCitation,&cite);
2356:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its));
2357:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol));

2359:   PCApply_PFMG(pc,b,y);
2360:   PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,&oits));
2361:   *outits = oits;
2362:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2363:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
2364:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol));
2365:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its));
2366:   return(0);
2367: }

2369: PetscErrorCode PCSetUp_PFMG(PC pc)
2370: {
2371:   PetscErrorCode  ierr;
2372:   PC_PFMG         *ex = (PC_PFMG*) pc->data;
2373:   Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
2374:   PetscBool       flg;

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

2380:   /* create the hypre solver object and set its information */
2381:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
2382:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
2383:   PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
2384:   PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver));
2385:   return(0);
2386: }

2388: /*MC
2389:      PCPFMG - the hypre PFMG multigrid solver

2391:    Level: advanced

2393:    Options Database:
2394: + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner
2395: . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2396: . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2397: . -pc_pfmg_tol <tol> tolerance of PFMG
2398: . -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
2399: - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin

2401:    Notes:
2402:     This is for CELL-centered descretizations

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

2407: .seealso:  PCMG, MATHYPRESTRUCT
2408: M*/

2410: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
2411: {
2413:   PC_PFMG        *ex;

2416:   PetscNew(&ex); \
2417:   pc->data = ex;

2419:   ex->its            = 1;
2420:   ex->tol            = 1.e-8;
2421:   ex->relax_type     = 1;
2422:   ex->rap_type       = 0;
2423:   ex->num_pre_relax  = 1;
2424:   ex->num_post_relax = 1;
2425:   ex->max_levels     = 0;

2427:   pc->ops->setfromoptions  = PCSetFromOptions_PFMG;
2428:   pc->ops->view            = PCView_PFMG;
2429:   pc->ops->destroy         = PCDestroy_PFMG;
2430:   pc->ops->apply           = PCApply_PFMG;
2431:   pc->ops->applyrichardson = PCApplyRichardson_PFMG;
2432:   pc->ops->setup           = PCSetUp_PFMG;

2434:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
2435:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
2436:   return(0);
2437: }

2439: /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/

2441: /* we know we are working with a HYPRE_SStructMatrix */
2442: typedef struct {
2443:   MPI_Comm            hcomm;       /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
2444:   HYPRE_SStructSolver ss_solver;

2446:   /* keep copy of SYSPFMG options used so may view them */
2447:   PetscInt its;
2448:   double   tol;
2449:   PetscInt relax_type;
2450:   PetscInt num_pre_relax,num_post_relax;
2451: } PC_SysPFMG;

2453: PetscErrorCode PCDestroy_SysPFMG(PC pc)
2454: {
2456:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

2459:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2460:   MPI_Comm_free(&ex->hcomm);
2461:   PetscFree(pc->data);
2462:   return(0);
2463: }

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

2467: PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
2468: {
2470:   PetscBool      iascii;
2471:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

2474:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2475:   if (iascii) {
2476:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG preconditioning\n");
2477:     PetscViewerASCIIPrintf(viewer,"  max iterations %d\n",ex->its);
2478:     PetscViewerASCIIPrintf(viewer,"  tolerance %g\n",ex->tol);
2479:     PetscViewerASCIIPrintf(viewer,"  relax type %s\n",PFMGRelaxType[ex->relax_type]);
2480:     PetscViewerASCIIPrintf(viewer,"  number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
2481:   }
2482:   return(0);
2483: }

2485: PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptionItems *PetscOptionsObject,PC pc)
2486: {
2488:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;
2489:   PetscBool      flg = PETSC_FALSE;

2492:   PetscOptionsHead(PetscOptionsObject,"SysPFMG options");
2493:   PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);
2494:   if (flg) {
2495:     PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,3));
2496:   }
2497:   PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);
2498:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its));
2499:   PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
2500:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax));
2501:   PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
2502:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax));

2504:   PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);
2505:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol));
2506:   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);
2507:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type));
2508:   PetscOptionsTail();
2509:   return(0);
2510: }

2512: PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
2513: {
2514:   PetscErrorCode    ierr;
2515:   PC_SysPFMG        *ex = (PC_SysPFMG*) pc->data;
2516:   PetscScalar       *yy;
2517:   const PetscScalar *xx;
2518:   PetscInt          ilower[3],iupper[3];
2519:   HYPRE_Int         hlower[3],hupper[3];
2520:   Mat_HYPRESStruct  *mx     = (Mat_HYPRESStruct*)(pc->pmat->data);
2521:   PetscInt          ordering= mx->dofs_order;
2522:   PetscInt          nvars   = mx->nvars;
2523:   PetscInt          part    = 0;
2524:   PetscInt          size;
2525:   PetscInt          i;

2528:   PetscCitationsRegister(hypreCitation,&cite);
2529:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
2530:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2531:   iupper[0] += ilower[0] - 1;
2532:   iupper[1] += ilower[1] - 1;
2533:   iupper[2] += ilower[2] - 1;
2534:   hlower[0]  = (HYPRE_Int)ilower[0];
2535:   hlower[1]  = (HYPRE_Int)ilower[1];
2536:   hlower[2]  = (HYPRE_Int)ilower[2];
2537:   hupper[0]  = (HYPRE_Int)iupper[0];
2538:   hupper[1]  = (HYPRE_Int)iupper[1];
2539:   hupper[2]  = (HYPRE_Int)iupper[2];

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

2544:   /* copy x values over to hypre for variable ordering */
2545:   if (ordering) {
2546:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2547:     VecGetArrayRead(x,&xx);
2548:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(xx+(size*i))));
2549:     VecRestoreArrayRead(x,&xx);
2550:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2551:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
2552:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

2554:     /* copy solution values back to PETSc */
2555:     VecGetArray(y,&yy);
2556:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(yy+(size*i))));
2557:     VecRestoreArray(y,&yy);
2558:   } else {      /* nodal ordering must be mapped to variable ordering for sys_pfmg */
2559:     PetscScalar *z;
2560:     PetscInt    j, k;

2562:     PetscMalloc1(nvars*size,&z);
2563:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2564:     VecGetArrayRead(x,&xx);

2566:     /* transform nodal to hypre's variable ordering for sys_pfmg */
2567:     for (i= 0; i< size; i++) {
2568:       k= i*nvars;
2569:       for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
2570:     }
2571:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i))));
2572:     VecRestoreArrayRead(x,&xx);
2573:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2574:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

2576:     /* copy solution values back to PETSc */
2577:     VecGetArray(y,&yy);
2578:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i))));
2579:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
2580:     for (i= 0; i< size; i++) {
2581:       k= i*nvars;
2582:       for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
2583:     }
2584:     VecRestoreArray(y,&yy);
2585:     PetscFree(z);
2586:   }
2587:   return(0);
2588: }

2590: 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)
2591: {
2592:   PC_SysPFMG     *jac = (PC_SysPFMG*)pc->data;
2594:   HYPRE_Int      oits;

2597:   PetscCitationsRegister(hypreCitation,&cite);
2598:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its));
2599:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol));
2600:   PCApply_SysPFMG(pc,b,y);
2601:   PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,&oits));
2602:   *outits = oits;
2603:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2604:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
2605:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol));
2606:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its));
2607:   return(0);
2608: }

2610: PetscErrorCode PCSetUp_SysPFMG(PC pc)
2611: {
2612:   PetscErrorCode   ierr;
2613:   PC_SysPFMG       *ex = (PC_SysPFMG*) pc->data;
2614:   Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
2615:   PetscBool        flg;

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

2621:   /* create the hypre sstruct solver object and set its information */
2622:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2623:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2624:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver));
2625:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2626:   return(0);
2627: }

2629: /*MC
2630:      PCSysPFMG - the hypre SysPFMG multigrid solver

2632:    Level: advanced

2634:    Options Database:
2635: + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner
2636: . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2637: . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2638: . -pc_syspfmg_tol <tol> tolerance of SysPFMG
2639: - -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel

2641:    Notes:
2642:     This is for CELL-centered descretizations

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

2648: .seealso:  PCMG, MATHYPRESSTRUCT
2649: M*/

2651: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
2652: {
2654:   PC_SysPFMG     *ex;

2657:   PetscNew(&ex); \
2658:   pc->data = ex;

2660:   ex->its            = 1;
2661:   ex->tol            = 1.e-8;
2662:   ex->relax_type     = 1;
2663:   ex->num_pre_relax  = 1;
2664:   ex->num_post_relax = 1;

2666:   pc->ops->setfromoptions  = PCSetFromOptions_SysPFMG;
2667:   pc->ops->view            = PCView_SysPFMG;
2668:   pc->ops->destroy         = PCDestroy_SysPFMG;
2669:   pc->ops->apply           = PCApply_SysPFMG;
2670:   pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
2671:   pc->ops->setup           = PCSetUp_SysPFMG;

2673:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
2674:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2675:   return(0);
2676: }