xref: /petsc/src/ksp/pc/impls/hypre/hypre.c (revision 48a46eb9bd028bec07ec0f396b1a3abb43f14558)

/*
   Provides an interface to the LLNL package hypre
*/

#include <petscpkg_version.h>
#include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
/* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */
#include <petsc/private/matimpl.h>
#include <petsc/private/vecimpl.h>
#include <../src/vec/vec/impls/hypre/vhyp.h>
#include <../src/mat/impls/hypre/mhypre.h>
#include <../src/dm/impls/da/hypre/mhyp.h>
#include <_hypre_parcsr_ls.h>
#include <petscmathypre.h>

#if defined(PETSC_HAVE_HYPRE_DEVICE)
#include <petsc/private/deviceimpl.h>
#endif

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

/*
   Private context (data structure) for the  preconditioner.
*/
typedef struct {
  HYPRE_Solver hsolver;
  Mat          hpmat; /* MatHYPRE */

  HYPRE_Int (*destroy)(HYPRE_Solver);
  HYPRE_Int (*solve)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector);
  HYPRE_Int (*setup)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector);

  MPI_Comm comm_hypre;
  char    *hypre_type;

  /* options for Pilut and BoomerAMG*/
  PetscInt  maxiter;
  PetscReal tol;

  /* options for Pilut */
  PetscInt factorrowsize;

  /* options for ParaSails */
  PetscInt  nlevels;
  PetscReal threshold;
  PetscReal filter;
  PetscReal loadbal;
  PetscInt  logging;
  PetscInt  ruse;
  PetscInt  symt;

  /* options for BoomerAMG */
  PetscBool printstatistics;

  /* options for BoomerAMG */
  PetscInt  cycletype;
  PetscInt  maxlevels;
  PetscReal strongthreshold;
  PetscReal maxrowsum;
  PetscInt  gridsweeps[3];
  PetscInt  coarsentype;
  PetscInt  measuretype;
  PetscInt  smoothtype;
  PetscInt  smoothnumlevels;
  PetscInt  eu_level;         /* Number of levels for ILU(k) in Euclid */
  PetscReal eu_droptolerance; /* Drop tolerance for ILU(k) in Euclid */
  PetscInt  eu_bj;            /* Defines use of Block Jacobi ILU in Euclid */
  PetscInt  relaxtype[3];
  PetscReal relaxweight;
  PetscReal outerrelaxweight;
  PetscInt  relaxorder;
  PetscReal truncfactor;
  PetscBool applyrichardson;
  PetscInt  pmax;
  PetscInt  interptype;
  PetscInt  maxc;
  PetscInt  minc;
#if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
  char *spgemm_type; // this is a global hypre parameter but is closely associated with BoomerAMG
#endif
  /* GPU */
  PetscBool keeptranspose;
  PetscInt  rap2;
  PetscInt  mod_rap2;

  /* AIR */
  PetscInt  Rtype;
  PetscReal Rstrongthreshold;
  PetscReal Rfilterthreshold;
  PetscInt  Adroptype;
  PetscReal Adroptol;

  PetscInt  agg_nl;
  PetscInt  agg_interptype;
  PetscInt  agg_num_paths;
  PetscBool nodal_relax;
  PetscInt  nodal_relax_levels;

  PetscInt  nodal_coarsening;
  PetscInt  nodal_coarsening_diag;
  PetscInt  vec_interp_variant;
  PetscInt  vec_interp_qmax;
  PetscBool vec_interp_smooth;
  PetscInt  interp_refine;

  /* NearNullSpace support */
  VecHYPRE_IJVector *hmnull;
  HYPRE_ParVector   *phmnull;
  PetscInt           n_hmnull;
  Vec                hmnull_constant;

  /* options for AS (Auxiliary Space preconditioners) */
  PetscInt  as_print;
  PetscInt  as_max_iter;
  PetscReal as_tol;
  PetscInt  as_relax_type;
  PetscInt  as_relax_times;
  PetscReal as_relax_weight;
  PetscReal as_omega;
  PetscInt  as_amg_alpha_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for vector Poisson (AMS) or Curl problem (ADS) */
  PetscReal as_amg_alpha_theta;   /* AMG strength for vector Poisson (AMS) or Curl problem (ADS) */
  PetscInt  as_amg_beta_opts[5];  /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for scalar Poisson (AMS) or vector Poisson (ADS) */
  PetscReal as_amg_beta_theta;    /* AMG strength for scalar Poisson (AMS) or vector Poisson (ADS)  */
  PetscInt  ams_cycle_type;
  PetscInt  ads_cycle_type;

  /* additional data */
  Mat G;             /* MatHYPRE */
  Mat C;             /* MatHYPRE */
  Mat alpha_Poisson; /* MatHYPRE */
  Mat beta_Poisson;  /* MatHYPRE */

  /* extra information for AMS */
  PetscInt          dim; /* geometrical dimension */
  VecHYPRE_IJVector coords[3];
  VecHYPRE_IJVector constants[3];
  VecHYPRE_IJVector interior;
  Mat               RT_PiFull, RT_Pi[3];
  Mat               ND_PiFull, ND_Pi[3];
  PetscBool         ams_beta_is_zero;
  PetscBool         ams_beta_is_zero_part;
  PetscInt          ams_proj_freq;
} PC_HYPRE;

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

  PetscFunctionBegin;
  *hsolver = jac->hsolver;
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscCall(PetscStrcmp(jac->hypre_type, "boomeramg", &same));
  PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_NOTSAMETYPE, "Hypre type is not BoomerAMG ");
  num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData *)(jac->hsolver));
  PetscCall(PetscMalloc1(num_levels, &mattmp));
  A_array = hypre_ParAMGDataAArray((hypre_ParAMGData *)(jac->hsolver));
  for (l = 1; l < num_levels; l++) {
    PetscCall(MatCreateFromParCSR(A_array[l], MATAIJ, PETSC_OWN_POINTER, &(mattmp[num_levels - 1 - l])));
    /* We want to own the data, and HYPRE can not touch this matrix any more */
    A_array[l] = NULL;
  }
  *nlevels   = num_levels;
  *operators = mattmp;
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscCall(PetscStrcmp(jac->hypre_type, "boomeramg", &same));
  PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_NOTSAMETYPE, "Hypre type is not BoomerAMG ");
  num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData *)(jac->hsolver));
  PetscCall(PetscMalloc1(num_levels, &mattmp));
  P_array = hypre_ParAMGDataPArray((hypre_ParAMGData *)(jac->hsolver));
  for (l = 1; l < num_levels; l++) {
    PetscCall(MatCreateFromParCSR(P_array[num_levels - 1 - l], MATAIJ, PETSC_OWN_POINTER, &(mattmp[l - 1])));
    /* We want to own the data, and HYPRE can not touch this matrix any more */
    P_array[num_levels - 1 - l] = NULL;
  }
  *nlevels        = num_levels;
  *interpolations = mattmp;
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  for (i = 0; i < jac->n_hmnull; i++) PetscCall(VecHYPRE_IJVectorDestroy(&jac->hmnull[i]));
  PetscCall(PetscFree(jac->hmnull));
  PetscCall(PetscFree(jac->phmnull));
  PetscCall(VecDestroy(&jac->hmnull_constant));
  jac->n_hmnull = 0;
  PetscFunctionReturn(0);
}

static PetscErrorCode PCSetUp_HYPRE(PC pc) {
  PC_HYPRE          *jac = (PC_HYPRE *)pc->data;
  Mat_HYPRE         *hjac;
  HYPRE_ParCSRMatrix hmat;
  HYPRE_ParVector    bv, xv;
  PetscBool          ishypre;

  PetscFunctionBegin;
  if (!jac->hypre_type) PetscCall(PCHYPRESetType(pc, "boomeramg"));

  PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRE, &ishypre));
  if (!ishypre) {
    PetscCall(MatDestroy(&jac->hpmat));
    PetscCall(MatConvert(pc->pmat, MATHYPRE, MAT_INITIAL_MATRIX, &jac->hpmat));
    PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)jac->hpmat));
  } else {
    PetscCall(PetscObjectReference((PetscObject)pc->pmat));
    PetscCall(MatDestroy(&jac->hpmat));
    jac->hpmat = pc->pmat;
  }
  /* allow debug */
  PetscCall(MatViewFromOptions(jac->hpmat, NULL, "-pc_hypre_mat_view"));
  hjac = (Mat_HYPRE *)(jac->hpmat->data);

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

    PetscCall(MatGetBlockSize(pc->pmat, &bs));
    if (bs > 1) PetscCallExternal(HYPRE_BoomerAMGSetNumFunctions, jac->hsolver, bs);
    PetscCall(MatGetNearNullSpace(pc->mat, &mnull));
    if (mnull) {
      PetscCall(PCHYPREResetNearNullSpace_Private(pc));
      PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs));
      PetscCall(PetscMalloc1(nvec + 1, &jac->hmnull));
      PetscCall(PetscMalloc1(nvec + 1, &jac->phmnull));
      for (i = 0; i < nvec; i++) {
        PetscCall(VecHYPRE_IJVectorCreate(vecs[i]->map, &jac->hmnull[i]));
        PetscCall(VecHYPRE_IJVectorCopy(vecs[i], jac->hmnull[i]));
        PetscCallExternal(HYPRE_IJVectorGetObject, jac->hmnull[i]->ij, (void **)&jac->phmnull[i]);
      }
      if (has_const) {
        PetscCall(MatCreateVecs(pc->pmat, &jac->hmnull_constant, NULL));
        PetscCall(PetscLogObjectParent((PetscObject)pc, (PetscObject)jac->hmnull_constant));
        PetscCall(VecSet(jac->hmnull_constant, 1));
        PetscCall(VecNormalize(jac->hmnull_constant, NULL));
        PetscCall(VecHYPRE_IJVectorCreate(jac->hmnull_constant->map, &jac->hmnull[nvec]));
        PetscCall(VecHYPRE_IJVectorCopy(jac->hmnull_constant, jac->hmnull[nvec]));
        PetscCallExternal(HYPRE_IJVectorGetObject, jac->hmnull[nvec]->ij, (void **)&jac->phmnull[nvec]);
        nvec++;
      }
      PetscCallExternal(HYPRE_BoomerAMGSetInterpVectors, jac->hsolver, nvec, jac->phmnull);
      jac->n_hmnull = nvec;
    }
  }

  /* special case for AMS */
  if (jac->setup == HYPRE_AMSSetup) {
    Mat_HYPRE         *hm;
    HYPRE_ParCSRMatrix parcsr;
    if (!jac->coords[0] && !jac->constants[0] && !(jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]))) {
      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");
    }
    if (jac->dim) PetscCallExternal(HYPRE_AMSSetDimension, jac->hsolver, jac->dim);
    if (jac->constants[0]) {
      HYPRE_ParVector ozz, zoz, zzo = NULL;
      PetscCallExternal(HYPRE_IJVectorGetObject, jac->constants[0]->ij, (void **)(&ozz));
      PetscCallExternal(HYPRE_IJVectorGetObject, jac->constants[1]->ij, (void **)(&zoz));
      if (jac->constants[2]) PetscCallExternal(HYPRE_IJVectorGetObject, jac->constants[2]->ij, (void **)(&zzo));
      PetscCallExternal(HYPRE_AMSSetEdgeConstantVectors, jac->hsolver, ozz, zoz, zzo);
    }
    if (jac->coords[0]) {
      HYPRE_ParVector coords[3];
      coords[0] = NULL;
      coords[1] = NULL;
      coords[2] = NULL;
      if (jac->coords[0]) PetscCallExternal(HYPRE_IJVectorGetObject, jac->coords[0]->ij, (void **)(&coords[0]));
      if (jac->coords[1]) PetscCallExternal(HYPRE_IJVectorGetObject, jac->coords[1]->ij, (void **)(&coords[1]));
      if (jac->coords[2]) PetscCallExternal(HYPRE_IJVectorGetObject, jac->coords[2]->ij, (void **)(&coords[2]));
      PetscCallExternal(HYPRE_AMSSetCoordinateVectors, jac->hsolver, coords[0], coords[1], coords[2]);
    }
    PetscCheck(jac->G, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE AMS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient");
    hm = (Mat_HYPRE *)(jac->G->data);
    PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&parcsr));
    PetscCallExternal(HYPRE_AMSSetDiscreteGradient, jac->hsolver, parcsr);
    if (jac->alpha_Poisson) {
      hm = (Mat_HYPRE *)(jac->alpha_Poisson->data);
      PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&parcsr));
      PetscCallExternal(HYPRE_AMSSetAlphaPoissonMatrix, jac->hsolver, parcsr);
    }
    if (jac->ams_beta_is_zero) {
      PetscCallExternal(HYPRE_AMSSetBetaPoissonMatrix, jac->hsolver, NULL);
    } else if (jac->beta_Poisson) {
      hm = (Mat_HYPRE *)(jac->beta_Poisson->data);
      PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&parcsr));
      PetscCallExternal(HYPRE_AMSSetBetaPoissonMatrix, jac->hsolver, parcsr);
    } else if (jac->ams_beta_is_zero_part) {
      if (jac->interior) {
        HYPRE_ParVector interior = NULL;
        PetscCallExternal(HYPRE_IJVectorGetObject, jac->interior->ij, (void **)(&interior));
        PetscCallExternal(HYPRE_AMSSetInteriorNodes, jac->hsolver, interior);
      } else {
        jac->ams_beta_is_zero_part = PETSC_FALSE;
      }
    }
    if (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1])) {
      PetscInt           i;
      HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3];
      if (jac->ND_PiFull) {
        hm = (Mat_HYPRE *)(jac->ND_PiFull->data);
        PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&nd_parcsrfull));
      } else {
        nd_parcsrfull = NULL;
      }
      for (i = 0; i < 3; ++i) {
        if (jac->ND_Pi[i]) {
          hm = (Mat_HYPRE *)(jac->ND_Pi[i]->data);
          PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&nd_parcsr[i]));
        } else {
          nd_parcsr[i] = NULL;
        }
      }
      PetscCallExternal(HYPRE_AMSSetInterpolations, jac->hsolver, nd_parcsrfull, nd_parcsr[0], nd_parcsr[1], nd_parcsr[2]);
    }
  }
  /* special case for ADS */
  if (jac->setup == HYPRE_ADSSetup) {
    Mat_HYPRE         *hm;
    HYPRE_ParCSRMatrix parcsr;
    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])))) {
      SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE ADS preconditioner needs either the coordinate vectors via PCSetCoordinates() or the interpolation matrices via PCHYPRESetInterpolations");
    } else PetscCheck(jac->coords[1] && jac->coords[2], PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE ADS preconditioner has been designed for three dimensional problems! For two dimensional problems, use HYPRE AMS instead");
    PetscCheck(jac->G, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE ADS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient");
    PetscCheck(jac->C, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE ADS preconditioner needs the discrete curl operator via PCHYPRESetDiscreteGradient");
    if (jac->coords[0]) {
      HYPRE_ParVector coords[3];
      coords[0] = NULL;
      coords[1] = NULL;
      coords[2] = NULL;
      if (jac->coords[0]) PetscCallExternal(HYPRE_IJVectorGetObject, jac->coords[0]->ij, (void **)(&coords[0]));
      if (jac->coords[1]) PetscCallExternal(HYPRE_IJVectorGetObject, jac->coords[1]->ij, (void **)(&coords[1]));
      if (jac->coords[2]) PetscCallExternal(HYPRE_IJVectorGetObject, jac->coords[2]->ij, (void **)(&coords[2]));
      PetscCallExternal(HYPRE_ADSSetCoordinateVectors, jac->hsolver, coords[0], coords[1], coords[2]);
    }
    hm = (Mat_HYPRE *)(jac->G->data);
    PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&parcsr));
    PetscCallExternal(HYPRE_ADSSetDiscreteGradient, jac->hsolver, parcsr);
    hm = (Mat_HYPRE *)(jac->C->data);
    PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&parcsr));
    PetscCallExternal(HYPRE_ADSSetDiscreteCurl, jac->hsolver, parcsr);
    if ((jac->RT_PiFull || (jac->RT_Pi[0] && jac->RT_Pi[1])) && (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]))) {
      PetscInt           i;
      HYPRE_ParCSRMatrix rt_parcsrfull, rt_parcsr[3];
      HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3];
      if (jac->RT_PiFull) {
        hm = (Mat_HYPRE *)(jac->RT_PiFull->data);
        PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&rt_parcsrfull));
      } else {
        rt_parcsrfull = NULL;
      }
      for (i = 0; i < 3; ++i) {
        if (jac->RT_Pi[i]) {
          hm = (Mat_HYPRE *)(jac->RT_Pi[i]->data);
          PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&rt_parcsr[i]));
        } else {
          rt_parcsr[i] = NULL;
        }
      }
      if (jac->ND_PiFull) {
        hm = (Mat_HYPRE *)(jac->ND_PiFull->data);
        PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&nd_parcsrfull));
      } else {
        nd_parcsrfull = NULL;
      }
      for (i = 0; i < 3; ++i) {
        if (jac->ND_Pi[i]) {
          hm = (Mat_HYPRE *)(jac->ND_Pi[i]->data);
          PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&nd_parcsr[i]));
        } else {
          nd_parcsr[i] = NULL;
        }
      }
      PetscCallExternal(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]);
    }
  }
  PetscCallExternal(HYPRE_IJMatrixGetObject, hjac->ij, (void **)&hmat);
  PetscCallExternal(HYPRE_IJVectorGetObject, hjac->b->ij, (void **)&bv);
  PetscCallExternal(HYPRE_IJVectorGetObject, hjac->x->ij, (void **)&xv);
  PetscCallExternal(jac->setup, jac->hsolver, hmat, bv, xv);
  PetscFunctionReturn(0);
}

static PetscErrorCode PCApply_HYPRE(PC pc, Vec b, Vec x) {
  PC_HYPRE          *jac  = (PC_HYPRE *)pc->data;
  Mat_HYPRE         *hjac = (Mat_HYPRE *)(jac->hpmat->data);
  HYPRE_ParCSRMatrix hmat;
  HYPRE_ParVector    jbv, jxv;

  PetscFunctionBegin;
  PetscCall(PetscCitationsRegister(hypreCitation, &cite));
  if (!jac->applyrichardson) PetscCall(VecSet(x, 0.0));
  PetscCall(VecHYPRE_IJVectorPushVecRead(hjac->b, b));
  if (jac->applyrichardson) PetscCall(VecHYPRE_IJVectorPushVec(hjac->x, x));
  else PetscCall(VecHYPRE_IJVectorPushVecWrite(hjac->x, x));
  PetscCallExternal(HYPRE_IJMatrixGetObject, hjac->ij, (void **)&hmat);
  PetscCallExternal(HYPRE_IJVectorGetObject, hjac->b->ij, (void **)&jbv);
  PetscCallExternal(HYPRE_IJVectorGetObject, hjac->x->ij, (void **)&jxv);
  PetscStackCallExternalVoid(
    "Hypre solve", do {
      HYPRE_Int hierr = (*jac->solve)(jac->hsolver, hmat, jbv, jxv);
      if (hierr) {
        PetscCheck(hierr == HYPRE_ERROR_CONV, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in HYPRE solver, error code %d", (int)hierr);
        hypre__global_error = 0;
      }
    } while (0));

  if (jac->setup == HYPRE_AMSSetup && jac->ams_beta_is_zero_part) PetscCallExternal(HYPRE_AMSProjectOutGradients, jac->hsolver, jxv);
  PetscCall(VecHYPRE_IJVectorPopVec(hjac->x));
  PetscCall(VecHYPRE_IJVectorPopVec(hjac->b));
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscCall(MatDestroy(&jac->hpmat));
  PetscCall(MatDestroy(&jac->G));
  PetscCall(MatDestroy(&jac->C));
  PetscCall(MatDestroy(&jac->alpha_Poisson));
  PetscCall(MatDestroy(&jac->beta_Poisson));
  PetscCall(MatDestroy(&jac->RT_PiFull));
  PetscCall(MatDestroy(&jac->RT_Pi[0]));
  PetscCall(MatDestroy(&jac->RT_Pi[1]));
  PetscCall(MatDestroy(&jac->RT_Pi[2]));
  PetscCall(MatDestroy(&jac->ND_PiFull));
  PetscCall(MatDestroy(&jac->ND_Pi[0]));
  PetscCall(MatDestroy(&jac->ND_Pi[1]));
  PetscCall(MatDestroy(&jac->ND_Pi[2]));
  PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[0]));
  PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[1]));
  PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[2]));
  PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[0]));
  PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[1]));
  PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[2]));
  PetscCall(VecHYPRE_IJVectorDestroy(&jac->interior));
  PetscCall(PCHYPREResetNearNullSpace_Private(pc));
  jac->ams_beta_is_zero      = PETSC_FALSE;
  jac->ams_beta_is_zero_part = PETSC_FALSE;
  jac->dim                   = 0;
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscCall(PCReset_HYPRE(pc));
  if (jac->destroy) PetscCallExternal(jac->destroy, jac->hsolver);
  PetscCall(PetscFree(jac->hypre_type));
#if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
  PetscCall(PetscFree(jac->spgemm_type));
#endif
  if (jac->comm_hypre != MPI_COMM_NULL) PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
  PetscCall(PetscFree(pc->data));

  PetscCall(PetscObjectChangeTypeName((PetscObject)pc, 0));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetType_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetType_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteGradient_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteCurl_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetInterpolations_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetConstantEdgeVectors_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetPoissonMatrix_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetEdgeConstantVectors_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREAMSSetInteriorNodes_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinSetMatProductAlgorithm_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinGetMatProductAlgorithm_C", NULL));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
  PetscFunctionReturn(0);
}

/* --------------------------------------------------------------------------------------------*/
static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PC pc, PetscOptionItems *PetscOptionsObject) {
  PC_HYPRE *jac = (PC_HYPRE *)pc->data;
  PetscBool flag;

  PetscFunctionBegin;
  PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE Pilut Options");
  PetscCall(PetscOptionsInt("-pc_hypre_pilut_maxiter", "Number of iterations", "None", jac->maxiter, &jac->maxiter, &flag));
  if (flag) PetscCallExternal(HYPRE_ParCSRPilutSetMaxIter, jac->hsolver, jac->maxiter);
  PetscCall(PetscOptionsReal("-pc_hypre_pilut_tol", "Drop tolerance", "None", jac->tol, &jac->tol, &flag));
  if (flag) PetscCallExternal(HYPRE_ParCSRPilutSetDropTolerance, jac->hsolver, jac->tol);
  PetscCall(PetscOptionsInt("-pc_hypre_pilut_factorrowsize", "FactorRowSize", "None", jac->factorrowsize, &jac->factorrowsize, &flag));
  if (flag) PetscCallExternal(HYPRE_ParCSRPilutSetFactorRowSize, jac->hsolver, jac->factorrowsize);
  PetscOptionsHeadEnd();
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
  if (iascii) {
    PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE Pilut preconditioning\n"));
    if (jac->maxiter != PETSC_DEFAULT) {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    maximum number of iterations %" PetscInt_FMT "\n", jac->maxiter));
    } else {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    default maximum number of iterations \n"));
    }
    if (jac->tol != PETSC_DEFAULT) {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    drop tolerance %g\n", (double)jac->tol));
    } else {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    default drop tolerance \n"));
    }
    if (jac->factorrowsize != PETSC_DEFAULT) {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    factor row size %" PetscInt_FMT "\n", jac->factorrowsize));
    } else {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    default factor row size \n"));
    }
  }
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE Euclid Options");
  PetscCall(PetscOptionsInt("-pc_hypre_euclid_level", "Factorization levels", "None", jac->eu_level, &jac->eu_level, &flag));
  if (flag) PetscCallExternal(HYPRE_EuclidSetLevel, jac->hsolver, jac->eu_level);

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

    PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
    PetscCheck(size == 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "hypre's Euclid does not support a parallel drop tolerance");
    PetscCallExternal(HYPRE_EuclidSetILUT, jac->hsolver, jac->eu_droptolerance);
  }

  PetscCall(PetscOptionsBool("-pc_hypre_euclid_bj", "Use Block Jacobi for ILU in Euclid", "None", eu_bj, &eu_bj, &flag));
  if (flag) {
    jac->eu_bj = eu_bj ? 1 : 0;
    PetscCallExternal(HYPRE_EuclidSetBJ, jac->hsolver, jac->eu_bj);
  }
  PetscOptionsHeadEnd();
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
  if (iascii) {
    PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE Euclid preconditioning\n"));
    if (jac->eu_level != PETSC_DEFAULT) {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    factorization levels %" PetscInt_FMT "\n", jac->eu_level));
    } else {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    default factorization levels \n"));
    }
    PetscCall(PetscViewerASCIIPrintf(viewer, "    drop tolerance %g\n", (double)jac->eu_droptolerance));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    use Block-Jacobi? %" PetscInt_FMT "\n", jac->eu_bj));
  }
  PetscFunctionReturn(0);
}

/* --------------------------------------------------------------------------------------------*/

static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc, Vec b, Vec x) {
  PC_HYPRE          *jac  = (PC_HYPRE *)pc->data;
  Mat_HYPRE         *hjac = (Mat_HYPRE *)(jac->hpmat->data);
  HYPRE_ParCSRMatrix hmat;
  HYPRE_ParVector    jbv, jxv;

  PetscFunctionBegin;
  PetscCall(PetscCitationsRegister(hypreCitation, &cite));
  PetscCall(VecSet(x, 0.0));
  PetscCall(VecHYPRE_IJVectorPushVecRead(hjac->x, b));
  PetscCall(VecHYPRE_IJVectorPushVecWrite(hjac->b, x));

  PetscCallExternal(HYPRE_IJMatrixGetObject, hjac->ij, (void **)&hmat);
  PetscCallExternal(HYPRE_IJVectorGetObject, hjac->b->ij, (void **)&jbv);
  PetscCallExternal(HYPRE_IJVectorGetObject, hjac->x->ij, (void **)&jxv);

  PetscStackCallExternalVoid(
    "Hypre Transpose solve", do {
      HYPRE_Int hierr = HYPRE_BoomerAMGSolveT(jac->hsolver, hmat, jbv, jxv);
      if (hierr) {
        /* error code of 1 in BoomerAMG merely means convergence not achieved */
        PetscCheck(hierr == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in HYPRE solver, error code %d", (int)hierr);
        hypre__global_error = 0;
      }
    } while (0));

  PetscCall(VecHYPRE_IJVectorPopVec(hjac->x));
  PetscCall(VecHYPRE_IJVectorPopVec(hjac->b));
  PetscFunctionReturn(0);
}

static PetscErrorCode PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG(PC pc, const char name[]) {
  PC_HYPRE *jac = (PC_HYPRE *)pc->data;
  PetscBool flag;

#if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
  PetscFunctionBegin;
  if (jac->spgemm_type) {
    PetscCall(PetscStrcmp(jac->spgemm_type, name, &flag));
    PetscCheck(flag, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Cannot reset the HYPRE SpGEMM (really we can)");
    PetscFunctionReturn(0);
  } else {
    PetscCall(PetscStrallocpy(name, &jac->spgemm_type));
  }
  PetscCall(PetscStrcmp("cusparse", jac->spgemm_type, &flag));
  if (flag) {
    PetscCallExternal(HYPRE_SetSpGemmUseCusparse, 1);
    PetscFunctionReturn(0);
  }
  PetscCall(PetscStrcmp("hypre", jac->spgemm_type, &flag));
  if (flag) {
    PetscCallExternal(HYPRE_SetSpGemmUseCusparse, 0);
    PetscFunctionReturn(0);
  }
  jac->spgemm_type = NULL;
  SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown HYPRE SpGEM type %s; Choices are cusparse, hypre", name);
#endif
}

static PetscErrorCode PCMGGalerkinGetMatProductAlgorithm_HYPRE_BoomerAMG(PC pc, const char *spgemm[]) {
  PC_HYPRE *jac = (PC_HYPRE *)pc->data;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
#if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
  *spgemm = jac->spgemm_type;
#endif
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE BoomerAMG Options");
  PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_cycle_type", "Cycle type", "None", HYPREBoomerAMGCycleType + 1, 2, HYPREBoomerAMGCycleType[jac->cycletype], &indx, &flg));
  if (flg) {
    jac->cycletype = indx + 1;
    PetscCallExternal(HYPRE_BoomerAMGSetCycleType, jac->hsolver, jac->cycletype);
  }
  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_max_levels", "Number of levels (of grids) allowed", "None", jac->maxlevels, &jac->maxlevels, &flg));
  if (flg) {
    PetscCheck(jac->maxlevels >= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Number of levels %" PetscInt_FMT " must be at least two", jac->maxlevels);
    PetscCallExternal(HYPRE_BoomerAMGSetMaxLevels, jac->hsolver, jac->maxlevels);
  }
  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_max_iter", "Maximum iterations used PER hypre call", "None", jac->maxiter, &jac->maxiter, &flg));
  if (flg) {
    PetscCheck(jac->maxiter >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Number of iterations %" PetscInt_FMT " must be at least one", jac->maxiter);
    PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, jac->maxiter);
  }
  PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_tol", "Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)", "None", jac->tol, &jac->tol, &flg));
  if (flg) {
    PetscCheck(jac->tol >= 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Tolerance %g must be greater than or equal to zero", (double)jac->tol);
    PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, jac->tol);
  }
  bs = 1;
  if (pc->pmat) PetscCall(MatGetBlockSize(pc->pmat, &bs));
  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_numfunctions", "Number of functions", "HYPRE_BoomerAMGSetNumFunctions", bs, &bs, &flg));
  if (flg) PetscCallExternal(HYPRE_BoomerAMGSetNumFunctions, jac->hsolver, bs);

  PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_truncfactor", "Truncation factor for interpolation (0=no truncation)", "None", jac->truncfactor, &jac->truncfactor, &flg));
  if (flg) {
    PetscCheck(jac->truncfactor >= 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Truncation factor %g must be great than or equal zero", (double)jac->truncfactor);
    PetscCallExternal(HYPRE_BoomerAMGSetTruncFactor, jac->hsolver, jac->truncfactor);
  }

  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_P_max", "Max elements per row for interpolation operator (0=unlimited)", "None", jac->pmax, &jac->pmax, &flg));
  if (flg) {
    PetscCheck(jac->pmax >= 0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "P_max %" PetscInt_FMT " must be greater than or equal to zero", jac->pmax);
    PetscCallExternal(HYPRE_BoomerAMGSetPMaxElmts, jac->hsolver, jac->pmax);
  }

  PetscCall(PetscOptionsRangeInt("-pc_hypre_boomeramg_agg_nl", "Number of levels of aggressive coarsening", "None", jac->agg_nl, &jac->agg_nl, &flg, 0, jac->maxlevels));
  if (flg) PetscCallExternal(HYPRE_BoomerAMGSetAggNumLevels, jac->hsolver, jac->agg_nl);

  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_agg_num_paths", "Number of paths for aggressive coarsening", "None", jac->agg_num_paths, &jac->agg_num_paths, &flg));
  if (flg) {
    PetscCheck(jac->agg_num_paths >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Number of paths %" PetscInt_FMT " must be greater than or equal to 1", jac->agg_num_paths);
    PetscCallExternal(HYPRE_BoomerAMGSetNumPaths, jac->hsolver, jac->agg_num_paths);
  }

  PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_strong_threshold", "Threshold for being strongly connected", "None", jac->strongthreshold, &jac->strongthreshold, &flg));
  if (flg) {
    PetscCheck(jac->strongthreshold >= 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Strong threshold %g must be great than or equal zero", (double)jac->strongthreshold);
    PetscCallExternal(HYPRE_BoomerAMGSetStrongThreshold, jac->hsolver, jac->strongthreshold);
  }
  PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_max_row_sum", "Maximum row sum", "None", jac->maxrowsum, &jac->maxrowsum, &flg));
  if (flg) {
    PetscCheck(jac->maxrowsum >= 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Maximum row sum %g must be greater than zero", (double)jac->maxrowsum);
    PetscCheck(jac->maxrowsum <= 1.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Maximum row sum %g must be less than or equal one", (double)jac->maxrowsum);
    PetscCallExternal(HYPRE_BoomerAMGSetMaxRowSum, jac->hsolver, jac->maxrowsum);
  }

  /* Grid sweeps */
  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all", "Number of sweeps for the up and down grid levels", "None", jac->gridsweeps[0], &indx, &flg));
  if (flg) {
    PetscCallExternal(HYPRE_BoomerAMGSetNumSweeps, jac->hsolver, indx);
    /* modify the jac structure so we can view the updated options with PC_View */
    jac->gridsweeps[0] = indx;
    jac->gridsweeps[1] = indx;
    /*defaults coarse to 1 */
    jac->gridsweeps[2] = 1;
  }
  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen", "Use a nodal based coarsening 1-6", "HYPRE_BoomerAMGSetNodal", jac->nodal_coarsening, &jac->nodal_coarsening, &flg));
  if (flg) PetscCallExternal(HYPRE_BoomerAMGSetNodal, jac->hsolver, jac->nodal_coarsening);
  PetscCall(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));
  if (flg) PetscCallExternal(HYPRE_BoomerAMGSetNodalDiag, jac->hsolver, jac->nodal_coarsening_diag);
  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_variant", "Variant of algorithm 1-3", "HYPRE_BoomerAMGSetInterpVecVariant", jac->vec_interp_variant, &jac->vec_interp_variant, &flg));
  if (flg) PetscCallExternal(HYPRE_BoomerAMGSetInterpVecVariant, jac->hsolver, jac->vec_interp_variant);
  PetscCall(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));
  if (flg) PetscCallExternal(HYPRE_BoomerAMGSetInterpVecQMax, jac->hsolver, jac->vec_interp_qmax);
  PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_vec_interp_smooth", "Whether to smooth the interpolation vectors", "HYPRE_BoomerAMGSetSmoothInterpVectors", jac->vec_interp_smooth, &jac->vec_interp_smooth, &flg));
  if (flg) PetscCallExternal(HYPRE_BoomerAMGSetSmoothInterpVectors, jac->hsolver, jac->vec_interp_smooth);
  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_interp_refine", "Preprocess the interpolation matrix through iterative weight refinement", "HYPRE_BoomerAMGSetInterpRefine", jac->interp_refine, &jac->interp_refine, &flg));
  if (flg) PetscCallExternal(HYPRE_BoomerAMGSetInterpRefine, jac->hsolver, jac->interp_refine);
  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down", "Number of sweeps for the down cycles", "None", jac->gridsweeps[0], &indx, &flg));
  if (flg) {
    PetscCallExternal(HYPRE_BoomerAMGSetCycleNumSweeps, jac->hsolver, indx, 1);
    jac->gridsweeps[0] = indx;
  }
  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up", "Number of sweeps for the up cycles", "None", jac->gridsweeps[1], &indx, &flg));
  if (flg) {
    PetscCallExternal(HYPRE_BoomerAMGSetCycleNumSweeps, jac->hsolver, indx, 2);
    jac->gridsweeps[1] = indx;
  }
  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse", "Number of sweeps for the coarse level", "None", jac->gridsweeps[2], &indx, &flg));
  if (flg) {
    PetscCallExternal(HYPRE_BoomerAMGSetCycleNumSweeps, jac->hsolver, indx, 3);
    jac->gridsweeps[2] = indx;
  }

  /* Smooth type */
  PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_smooth_type", "Enable more complex smoothers", "None", HYPREBoomerAMGSmoothType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGSmoothType), HYPREBoomerAMGSmoothType[0], &indx, &flg));
  if (flg) {
    jac->smoothtype = indx;
    PetscCallExternal(HYPRE_BoomerAMGSetSmoothType, jac->hsolver, indx + 6);
    jac->smoothnumlevels = 25;
    PetscCallExternal(HYPRE_BoomerAMGSetSmoothNumLevels, jac->hsolver, 25);
  }

  /* Number of smoothing levels */
  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_levels", "Number of levels on which more complex smoothers are used", "None", 25, &indx, &flg));
  if (flg && (jac->smoothtype != -1)) {
    jac->smoothnumlevels = indx;
    PetscCallExternal(HYPRE_BoomerAMGSetSmoothNumLevels, jac->hsolver, indx);
  }

  /* Number of levels for ILU(k) for Euclid */
  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_eu_level", "Number of levels for ILU(k) in Euclid smoother", "None", 0, &indx, &flg));
  if (flg && (jac->smoothtype == 3)) {
    jac->eu_level = indx;
    PetscCallExternal(HYPRE_BoomerAMGSetEuLevel, jac->hsolver, indx);
  }

  /* Filter for ILU(k) for Euclid */
  double droptolerance;
  PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_eu_droptolerance", "Drop tolerance for ILU(k) in Euclid smoother", "None", 0, &droptolerance, &flg));
  if (flg && (jac->smoothtype == 3)) {
    jac->eu_droptolerance = droptolerance;
    PetscCallExternal(HYPRE_BoomerAMGSetEuLevel, jac->hsolver, droptolerance);
  }

  /* Use Block Jacobi ILUT for Euclid */
  PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_eu_bj", "Use Block Jacobi for ILU in Euclid smoother?", "None", PETSC_FALSE, &tmp_truth, &flg));
  if (flg && (jac->smoothtype == 3)) {
    jac->eu_bj = tmp_truth;
    PetscCallExternal(HYPRE_BoomerAMGSetEuBJ, jac->hsolver, jac->eu_bj);
  }

  /* Relax type */
  PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all", "Relax type for the up and down cycles", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), HYPREBoomerAMGRelaxType[6], &indx, &flg));
  if (flg) {
    jac->relaxtype[0] = jac->relaxtype[1] = indx;
    PetscCallExternal(HYPRE_BoomerAMGSetRelaxType, jac->hsolver, indx);
    /* by default, coarse type set to 9 */
    jac->relaxtype[2] = 9;
    PetscCallExternal(HYPRE_BoomerAMGSetCycleRelaxType, jac->hsolver, 9, 3);
  }
  PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down", "Relax type for the down cycles", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), HYPREBoomerAMGRelaxType[6], &indx, &flg));
  if (flg) {
    jac->relaxtype[0] = indx;
    PetscCallExternal(HYPRE_BoomerAMGSetCycleRelaxType, jac->hsolver, indx, 1);
  }
  PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up", "Relax type for the up cycles", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), HYPREBoomerAMGRelaxType[6], &indx, &flg));
  if (flg) {
    jac->relaxtype[1] = indx;
    PetscCallExternal(HYPRE_BoomerAMGSetCycleRelaxType, jac->hsolver, indx, 2);
  }
  PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse", "Relax type on coarse grid", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), HYPREBoomerAMGRelaxType[9], &indx, &flg));
  if (flg) {
    jac->relaxtype[2] = indx;
    PetscCallExternal(HYPRE_BoomerAMGSetCycleRelaxType, jac->hsolver, indx, 3);
  }

  /* Relaxation Weight */
  PetscCall(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));
  if (flg) {
    PetscCallExternal(HYPRE_BoomerAMGSetRelaxWt, jac->hsolver, tmpdbl);
    jac->relaxweight = tmpdbl;
  }

  n         = 2;
  twodbl[0] = twodbl[1] = 1.0;
  PetscCall(PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level", "Set the relaxation weight for a particular level (weight,level)", "None", twodbl, &n, &flg));
  if (flg) {
    if (n == 2) {
      indx = (int)PetscAbsReal(twodbl[1]);
      PetscCallExternal(HYPRE_BoomerAMGSetLevelRelaxWt, jac->hsolver, twodbl[0], indx);
    } else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Relax weight level: you must provide 2 values separated by a comma (and no space), you provided %" PetscInt_FMT, n);
  }

  /* Outer relaxation Weight */
  PetscCall(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));
  if (flg) {
    PetscCallExternal(HYPRE_BoomerAMGSetOuterWt, jac->hsolver, tmpdbl);
    jac->outerrelaxweight = tmpdbl;
  }

  n         = 2;
  twodbl[0] = twodbl[1] = 1.0;
  PetscCall(PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level", "Set the outer relaxation weight for a particular level (weight,level)", "None", twodbl, &n, &flg));
  if (flg) {
    if (n == 2) {
      indx = (int)PetscAbsReal(twodbl[1]);
      PetscCallExternal(HYPRE_BoomerAMGSetLevelOuterWt, jac->hsolver, twodbl[0], indx);
    } else SETERRQ(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 %" PetscInt_FMT, n);
  }

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

  if (flg && tmp_truth) {
    jac->relaxorder = 0;
    PetscCallExternal(HYPRE_BoomerAMGSetRelaxOrder, jac->hsolver, jac->relaxorder);
  }
  PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_measure_type", "Measure type", "None", HYPREBoomerAMGMeasureType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGMeasureType), HYPREBoomerAMGMeasureType[0], &indx, &flg));
  if (flg) {
    jac->measuretype = indx;
    PetscCallExternal(HYPRE_BoomerAMGSetMeasureType, jac->hsolver, jac->measuretype);
  }
  /* update list length 3/07 */
  PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type", "Coarsen type", "None", HYPREBoomerAMGCoarsenType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGCoarsenType), HYPREBoomerAMGCoarsenType[6], &indx, &flg));
  if (flg) {
    jac->coarsentype = indx;
    PetscCallExternal(HYPRE_BoomerAMGSetCoarsenType, jac->hsolver, jac->coarsentype);
  }

  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_max_coarse_size", "Maximum size of coarsest grid", "None", jac->maxc, &jac->maxc, &flg));
  if (flg) PetscCallExternal(HYPRE_BoomerAMGSetMaxCoarseSize, jac->hsolver, jac->maxc);
  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_min_coarse_size", "Minimum size of coarsest grid", "None", jac->minc, &jac->minc, &flg));
  if (flg) PetscCallExternal(HYPRE_BoomerAMGSetMinCoarseSize, jac->hsolver, jac->minc);
#if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
  // global parameter but is closely associated with BoomerAMG
  PetscCall(PetscOptionsEList("-pc_mg_galerkin_mat_product_algorithm", "Type of SpGEMM to use in hypre (only for now)", "PCMGGalerkinSetMatProductAlgorithm", PCHYPRESpgemmTypes, PETSC_STATIC_ARRAY_LENGTH(PCHYPRESpgemmTypes), PCHYPRESpgemmTypes[0], &indx, &flg));
  if (!flg) indx = 0;
  PetscCall(PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG(pc, PCHYPRESpgemmTypes[indx]));
#endif
  /* AIR */
#if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
  PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_restriction_type", "Type of AIR method (distance 1 or 2, 0 means no AIR)", "None", jac->Rtype, &jac->Rtype, NULL));
  PetscCallExternal(HYPRE_BoomerAMGSetRestriction, jac->hsolver, jac->Rtype);
  if (jac->Rtype) {
    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 */

    PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_strongthresholdR", "Threshold for R", "None", jac->Rstrongthreshold, &jac->Rstrongthreshold, NULL));
    PetscCallExternal(HYPRE_BoomerAMGSetStrongThresholdR, jac->hsolver, jac->Rstrongthreshold);

    PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_filterthresholdR", "Filter threshold for R", "None", jac->Rfilterthreshold, &jac->Rfilterthreshold, NULL));
    PetscCallExternal(HYPRE_BoomerAMGSetFilterThresholdR, jac->hsolver, jac->Rfilterthreshold);

    PetscCall(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));
    PetscCallExternal(HYPRE_BoomerAMGSetADropTol, jac->hsolver, jac->Adroptol);

    PetscCall(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));
    PetscCallExternal(HYPRE_BoomerAMGSetADropType, jac->hsolver, jac->Adroptype);
  }
#endif

#if PETSC_PKG_HYPRE_VERSION_LE(9, 9, 9)
  PetscCheck(!jac->Rtype || !jac->agg_nl, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "-pc_hypre_boomeramg_restriction_type (%" PetscInt_FMT ") and -pc_hypre_boomeramg_agg_nl (%" PetscInt_FMT ")", jac->Rtype, jac->agg_nl);
#endif

  /* new 3/07 */
  PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_interp_type", "Interpolation type", "None", HYPREBoomerAMGInterpType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGInterpType), HYPREBoomerAMGInterpType[0], &indx, &flg));
  if (flg || jac->Rtype) {
    if (flg) jac->interptype = indx;
    PetscCallExternal(HYPRE_BoomerAMGSetInterpType, jac->hsolver, jac->interptype);
  }

  PetscCall(PetscOptionsName("-pc_hypre_boomeramg_print_statistics", "Print statistics", "None", &flg));
  if (flg) {
    level = 3;
    PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_print_statistics", "Print statistics", "None", level, &level, NULL));

    jac->printstatistics = PETSC_TRUE;
    PetscCallExternal(HYPRE_BoomerAMGSetPrintLevel, jac->hsolver, level);
  }

  PetscCall(PetscOptionsName("-pc_hypre_boomeramg_print_debug", "Print debug information", "None", &flg));
  if (flg) {
    level = 3;
    PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_print_debug", "Print debug information", "None", level, &level, NULL));

    jac->printstatistics = PETSC_TRUE;
    PetscCallExternal(HYPRE_BoomerAMGSetDebugFlag, jac->hsolver, level);
  }

  PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg));
  if (flg && tmp_truth) {
    PetscInt tmp_int;
    PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", jac->nodal_relax_levels, &tmp_int, &flg));
    if (flg) jac->nodal_relax_levels = tmp_int;
    PetscCallExternal(HYPRE_BoomerAMGSetSmoothType, jac->hsolver, 6);
    PetscCallExternal(HYPRE_BoomerAMGSetDomainType, jac->hsolver, 1);
    PetscCallExternal(HYPRE_BoomerAMGSetOverlap, jac->hsolver, 0);
    PetscCallExternal(HYPRE_BoomerAMGSetSmoothNumLevels, jac->hsolver, jac->nodal_relax_levels);
  }

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

  /* options for ParaSails solvers */
  PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_parasails_sym", "Symmetry of matrix and preconditioner", "None", symtlist, PETSC_STATIC_ARRAY_LENGTH(symtlist), symtlist[0], &indx, &flg));
  if (flg) {
    jac->symt = indx;
    PetscCallExternal(HYPRE_BoomerAMGSetSym, jac->hsolver, jac->symt);
  }

  PetscOptionsHeadEnd();
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscCall(PetscCitationsRegister(hypreCitation, &cite));
  PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, its * jac->maxiter);
  PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, rtol);
  jac->applyrichardson = PETSC_TRUE;
  PetscCall(PCApply_HYPRE(pc, b, y));
  jac->applyrichardson = PETSC_FALSE;
  PetscCallExternal(HYPRE_BoomerAMGGetNumIterations, jac->hsolver, &oits);
  *outits = oits;
  if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
  else *reason = PCRICHARDSON_CONVERGED_RTOL;
  PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, jac->tol);
  PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, jac->maxiter);
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
  if (iascii) {
    PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE BoomerAMG preconditioning\n"));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    Cycle type %s\n", HYPREBoomerAMGCycleType[jac->cycletype]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    Maximum number of levels %" PetscInt_FMT "\n", jac->maxlevels));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    Maximum number of iterations PER hypre call %" PetscInt_FMT "\n", jac->maxiter));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    Convergence tolerance PER hypre call %g\n", (double)jac->tol));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    Threshold for strong coupling %g\n", (double)jac->strongthreshold));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    Interpolation truncation factor %g\n", (double)jac->truncfactor));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    Interpolation: max elements per row %" PetscInt_FMT "\n", jac->pmax));
    if (jac->interp_refine) PetscCall(PetscViewerASCIIPrintf(viewer, "    Interpolation: number of steps of weighted refinement %" PetscInt_FMT "\n", jac->interp_refine));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    Number of levels of aggressive coarsening %" PetscInt_FMT "\n", jac->agg_nl));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    Number of paths for aggressive coarsening %" PetscInt_FMT "\n", jac->agg_num_paths));

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

    PetscCall(PetscViewerASCIIPrintf(viewer, "    Sweeps down         %" PetscInt_FMT "\n", jac->gridsweeps[0]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    Sweeps up           %" PetscInt_FMT "\n", jac->gridsweeps[1]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    Sweeps on coarse    %" PetscInt_FMT "\n", jac->gridsweeps[2]));

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

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

    if (jac->relaxorder) {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    Using CF-relaxation\n"));
    } else {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    Not using CF-relaxation\n"));
    }
    if (jac->smoothtype != -1) {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    Smooth type          %s\n", HYPREBoomerAMGSmoothType[jac->smoothtype]));
      PetscCall(PetscViewerASCIIPrintf(viewer, "    Smooth num levels    %" PetscInt_FMT "\n", jac->smoothnumlevels));
    } else {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    Not using more complex smoothers.\n"));
    }
    if (jac->smoothtype == 3) {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    Euclid ILU(k) levels %" PetscInt_FMT "\n", jac->eu_level));
      PetscCall(PetscViewerASCIIPrintf(viewer, "    Euclid ILU(k) drop tolerance %g\n", (double)jac->eu_droptolerance));
      PetscCall(PetscViewerASCIIPrintf(viewer, "    Euclid ILU use Block-Jacobi? %" PetscInt_FMT "\n", jac->eu_bj));
    }
    PetscCall(PetscViewerASCIIPrintf(viewer, "    Measure type        %s\n", HYPREBoomerAMGMeasureType[jac->measuretype]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    Coarsen type        %s\n", HYPREBoomerAMGCoarsenType[jac->coarsentype]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    Interpolation type  %s\n", jac->interptype != 100 ? HYPREBoomerAMGInterpType[jac->interptype] : "1pt"));
    if (jac->nodal_coarsening) PetscCall(PetscViewerASCIIPrintf(viewer, "    Using nodal coarsening with HYPRE_BOOMERAMGSetNodal() %" PetscInt_FMT "\n", jac->nodal_coarsening));
    if (jac->vec_interp_variant) {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    HYPRE_BoomerAMGSetInterpVecVariant() %" PetscInt_FMT "\n", jac->vec_interp_variant));
      PetscCall(PetscViewerASCIIPrintf(viewer, "    HYPRE_BoomerAMGSetInterpVecQMax() %" PetscInt_FMT "\n", jac->vec_interp_qmax));
      PetscCall(PetscViewerASCIIPrintf(viewer, "    HYPRE_BoomerAMGSetSmoothInterpVectors() %d\n", jac->vec_interp_smooth));
    }
    if (jac->nodal_relax) PetscCall(PetscViewerASCIIPrintf(viewer, "    Using nodal relaxation via Schwarz smoothing on levels %" PetscInt_FMT "\n", jac->nodal_relax_levels));
#if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
    PetscCall(PetscViewerASCIIPrintf(viewer, "    SpGEMM type         %s\n", jac->spgemm_type));
#endif
    /* AIR */
    if (jac->Rtype) {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    Using approximate ideal restriction type %" PetscInt_FMT "\n", jac->Rtype));
      PetscCall(PetscViewerASCIIPrintf(viewer, "      Threshold for R %g\n", (double)jac->Rstrongthreshold));
      PetscCall(PetscViewerASCIIPrintf(viewer, "      Filter for R %g\n", (double)jac->Rfilterthreshold));
      PetscCall(PetscViewerASCIIPrintf(viewer, "      A drop tolerance %g\n", (double)jac->Adroptol));
      PetscCall(PetscViewerASCIIPrintf(viewer, "      A drop type %" PetscInt_FMT "\n", jac->Adroptype));
    }
  }
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE ParaSails Options");
  PetscCall(PetscOptionsInt("-pc_hypre_parasails_nlevels", "Number of number of levels", "None", jac->nlevels, &jac->nlevels, 0));
  PetscCall(PetscOptionsReal("-pc_hypre_parasails_thresh", "Threshold", "None", jac->threshold, &jac->threshold, &flag));
  if (flag) PetscCallExternal(HYPRE_ParaSailsSetParams, jac->hsolver, jac->threshold, jac->nlevels);

  PetscCall(PetscOptionsReal("-pc_hypre_parasails_filter", "filter", "None", jac->filter, &jac->filter, &flag));
  if (flag) PetscCallExternal(HYPRE_ParaSailsSetFilter, jac->hsolver, jac->filter);

  PetscCall(PetscOptionsReal("-pc_hypre_parasails_loadbal", "Load balance", "None", jac->loadbal, &jac->loadbal, &flag));
  if (flag) PetscCallExternal(HYPRE_ParaSailsSetLoadbal, jac->hsolver, jac->loadbal);

  PetscCall(PetscOptionsBool("-pc_hypre_parasails_logging", "Print info to screen", "None", (PetscBool)jac->logging, (PetscBool *)&jac->logging, &flag));
  if (flag) PetscCallExternal(HYPRE_ParaSailsSetLogging, jac->hsolver, jac->logging);

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

  PetscCall(PetscOptionsEList("-pc_hypre_parasails_sym", "Symmetry of matrix and preconditioner", "None", symtlist, PETSC_STATIC_ARRAY_LENGTH(symtlist), symtlist[0], &indx, &flag));
  if (flag) {
    jac->symt = indx;
    PetscCallExternal(HYPRE_ParaSailsSetSym, jac->hsolver, jac->symt);
  }

  PetscOptionsHeadEnd();
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
  if (iascii) {
    PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE ParaSails preconditioning\n"));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    nlevels %" PetscInt_FMT "\n", jac->nlevels));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    threshold %g\n", (double)jac->threshold));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    filter %g\n", (double)jac->filter));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    load balance %g\n", (double)jac->loadbal));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    reuse nonzero structure %s\n", PetscBools[jac->ruse]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    print info to screen %s\n", PetscBools[jac->logging]));
    if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
    else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
    else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
    else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Unknown HYPRE ParaSails symmetric option %" PetscInt_FMT, jac->symt);
    PetscCall(PetscViewerASCIIPrintf(viewer, "    %s\n", symt));
  }
  PetscFunctionReturn(0);
}
/* --------------------------------------------------------------------------------------------*/
static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PC pc, PetscOptionItems *PetscOptionsObject) {
  PC_HYPRE *jac = (PC_HYPRE *)pc->data;
  PetscInt  n;
  PetscBool flag, flag2, flag3, flag4;

  PetscFunctionBegin;
  PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE AMS Options");
  PetscCall(PetscOptionsInt("-pc_hypre_ams_print_level", "Debugging output level for AMS", "None", jac->as_print, &jac->as_print, &flag));
  if (flag) PetscCallExternal(HYPRE_AMSSetPrintLevel, jac->hsolver, jac->as_print);
  PetscCall(PetscOptionsInt("-pc_hypre_ams_max_iter", "Maximum number of AMS multigrid iterations within PCApply", "None", jac->as_max_iter, &jac->as_max_iter, &flag));
  if (flag) PetscCallExternal(HYPRE_AMSSetMaxIter, jac->hsolver, jac->as_max_iter);
  PetscCall(PetscOptionsInt("-pc_hypre_ams_cycle_type", "Cycle type for AMS multigrid", "None", jac->ams_cycle_type, &jac->ams_cycle_type, &flag));
  if (flag) PetscCallExternal(HYPRE_AMSSetCycleType, jac->hsolver, jac->ams_cycle_type);
  PetscCall(PetscOptionsReal("-pc_hypre_ams_tol", "Error tolerance for AMS multigrid", "None", jac->as_tol, &jac->as_tol, &flag));
  if (flag) PetscCallExternal(HYPRE_AMSSetTol, jac->hsolver, jac->as_tol);
  PetscCall(PetscOptionsInt("-pc_hypre_ams_relax_type", "Relaxation type for AMS smoother", "None", jac->as_relax_type, &jac->as_relax_type, &flag));
  PetscCall(PetscOptionsInt("-pc_hypre_ams_relax_times", "Number of relaxation steps for AMS smoother", "None", jac->as_relax_times, &jac->as_relax_times, &flag2));
  PetscCall(PetscOptionsReal("-pc_hypre_ams_relax_weight", "Relaxation weight for AMS smoother", "None", jac->as_relax_weight, &jac->as_relax_weight, &flag3));
  PetscCall(PetscOptionsReal("-pc_hypre_ams_omega", "SSOR coefficient for AMS smoother", "None", jac->as_omega, &jac->as_omega, &flag4));
  if (flag || flag2 || flag3 || flag4) PetscCallExternal(HYPRE_AMSSetSmoothingOptions, jac->hsolver, jac->as_relax_type, jac->as_relax_times, jac->as_relax_weight, jac->as_omega);
  PetscCall(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));
  n = 5;
  PetscCall(PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options", "AMG options for vector Poisson", "None", jac->as_amg_alpha_opts, &n, &flag2));
  if (flag || flag2) {
    PetscCallExternal(HYPRE_AMSSetAlphaAMGOptions, jac->hsolver, jac->as_amg_alpha_opts[0], /* AMG coarsen type */
                      jac->as_amg_alpha_opts[1],                                            /* AMG agg_levels */
                      jac->as_amg_alpha_opts[2],                                            /* AMG relax_type */
                      jac->as_amg_alpha_theta, jac->as_amg_alpha_opts[3],                   /* AMG interp_type */
                      jac->as_amg_alpha_opts[4]);                                           /* AMG Pmax */
  }
  PetscCall(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));
  n = 5;
  PetscCall(PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options", "AMG options for scalar Poisson solver", "None", jac->as_amg_beta_opts, &n, &flag2));
  if (flag || flag2) {
    PetscCallExternal(HYPRE_AMSSetBetaAMGOptions, jac->hsolver, jac->as_amg_beta_opts[0], /* AMG coarsen type */
                      jac->as_amg_beta_opts[1],                                           /* AMG agg_levels */
                      jac->as_amg_beta_opts[2],                                           /* AMG relax_type */
                      jac->as_amg_beta_theta, jac->as_amg_beta_opts[3],                   /* AMG interp_type */
                      jac->as_amg_beta_opts[4]);                                          /* AMG Pmax */
  }
  PetscCall(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));
  if (flag) { /* override HYPRE's default only if the options is used */
    PetscCallExternal(HYPRE_AMSSetProjectionFrequency, jac->hsolver, jac->ams_proj_freq);
  }
  PetscOptionsHeadEnd();
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
  if (iascii) {
    PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE AMS preconditioning\n"));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    subspace iterations per application %" PetscInt_FMT "\n", jac->as_max_iter));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    subspace cycle type %" PetscInt_FMT "\n", jac->ams_cycle_type));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    subspace iteration tolerance %g\n", (double)jac->as_tol));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    smoother type %" PetscInt_FMT "\n", jac->as_relax_type));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    number of smoothing steps %" PetscInt_FMT "\n", jac->as_relax_times));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    smoother weight %g\n", (double)jac->as_relax_weight));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    smoother omega %g\n", (double)jac->as_omega));
    if (jac->alpha_Poisson) {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    vector Poisson solver (passed in by user)\n"));
    } else {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    vector Poisson solver (computed) \n"));
    }
    PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG coarsening type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[0]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[1]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG relaxation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[2]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG interpolation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[3]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[4]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG strength threshold %g\n", (double)jac->as_amg_alpha_theta));
    if (!jac->ams_beta_is_zero) {
      if (jac->beta_Poisson) {
        PetscCall(PetscViewerASCIIPrintf(viewer, "    scalar Poisson solver (passed in by user)\n"));
      } else {
        PetscCall(PetscViewerASCIIPrintf(viewer, "    scalar Poisson solver (computed) \n"));
      }
      PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG coarsening type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[0]));
      PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_beta_opts[1]));
      PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG relaxation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[2]));
      PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG interpolation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[3]));
      PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_beta_opts[4]));
      PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG strength threshold %g\n", (double)jac->as_amg_beta_theta));
      if (jac->ams_beta_is_zero_part) PetscCall(PetscViewerASCIIPrintf(viewer, "        compatible subspace projection frequency %" PetscInt_FMT " (-1 HYPRE uses default)\n", jac->ams_proj_freq));
    } else {
      PetscCall(PetscViewerASCIIPrintf(viewer, "    scalar Poisson solver not used (zero-conductivity everywhere) \n"));
    }
  }
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE ADS Options");
  PetscCall(PetscOptionsInt("-pc_hypre_ads_print_level", "Debugging output level for ADS", "None", jac->as_print, &jac->as_print, &flag));
  if (flag) PetscCallExternal(HYPRE_ADSSetPrintLevel, jac->hsolver, jac->as_print);
  PetscCall(PetscOptionsInt("-pc_hypre_ads_max_iter", "Maximum number of ADS multigrid iterations within PCApply", "None", jac->as_max_iter, &jac->as_max_iter, &flag));
  if (flag) PetscCallExternal(HYPRE_ADSSetMaxIter, jac->hsolver, jac->as_max_iter);
  PetscCall(PetscOptionsInt("-pc_hypre_ads_cycle_type", "Cycle type for ADS multigrid", "None", jac->ads_cycle_type, &jac->ads_cycle_type, &flag));
  if (flag) PetscCallExternal(HYPRE_ADSSetCycleType, jac->hsolver, jac->ads_cycle_type);
  PetscCall(PetscOptionsReal("-pc_hypre_ads_tol", "Error tolerance for ADS multigrid", "None", jac->as_tol, &jac->as_tol, &flag));
  if (flag) PetscCallExternal(HYPRE_ADSSetTol, jac->hsolver, jac->as_tol);
  PetscCall(PetscOptionsInt("-pc_hypre_ads_relax_type", "Relaxation type for ADS smoother", "None", jac->as_relax_type, &jac->as_relax_type, &flag));
  PetscCall(PetscOptionsInt("-pc_hypre_ads_relax_times", "Number of relaxation steps for ADS smoother", "None", jac->as_relax_times, &jac->as_relax_times, &flag2));
  PetscCall(PetscOptionsReal("-pc_hypre_ads_relax_weight", "Relaxation weight for ADS smoother", "None", jac->as_relax_weight, &jac->as_relax_weight, &flag3));
  PetscCall(PetscOptionsReal("-pc_hypre_ads_omega", "SSOR coefficient for ADS smoother", "None", jac->as_omega, &jac->as_omega, &flag4));
  if (flag || flag2 || flag3 || flag4) PetscCallExternal(HYPRE_ADSSetSmoothingOptions, jac->hsolver, jac->as_relax_type, jac->as_relax_times, jac->as_relax_weight, jac->as_omega);
  PetscCall(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));
  n = 5;
  PetscCall(PetscOptionsIntArray("-pc_hypre_ads_ams_options", "AMG options for AMS solver inside ADS", "None", jac->as_amg_alpha_opts, &n, &flag2));
  PetscCall(PetscOptionsInt("-pc_hypre_ads_ams_cycle_type", "Cycle type for AMS solver inside ADS", "None", jac->ams_cycle_type, &jac->ams_cycle_type, &flag3));
  if (flag || flag2 || flag3) {
    PetscCallExternal(HYPRE_ADSSetAMSOptions, jac->hsolver, jac->ams_cycle_type, /* AMS cycle type */
                      jac->as_amg_alpha_opts[0],                                 /* AMG coarsen type */
                      jac->as_amg_alpha_opts[1],                                 /* AMG agg_levels */
                      jac->as_amg_alpha_opts[2],                                 /* AMG relax_type */
                      jac->as_amg_alpha_theta, jac->as_amg_alpha_opts[3],        /* AMG interp_type */
                      jac->as_amg_alpha_opts[4]);                                /* AMG Pmax */
  }
  PetscCall(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));
  n = 5;
  PetscCall(PetscOptionsIntArray("-pc_hypre_ads_amg_options", "AMG options for vector AMG solver inside ADS", "None", jac->as_amg_beta_opts, &n, &flag2));
  if (flag || flag2) {
    PetscCallExternal(HYPRE_ADSSetAMGOptions, jac->hsolver, jac->as_amg_beta_opts[0], /* AMG coarsen type */
                      jac->as_amg_beta_opts[1],                                       /* AMG agg_levels */
                      jac->as_amg_beta_opts[2],                                       /* AMG relax_type */
                      jac->as_amg_beta_theta, jac->as_amg_beta_opts[3],               /* AMG interp_type */
                      jac->as_amg_beta_opts[4]);                                      /* AMG Pmax */
  }
  PetscOptionsHeadEnd();
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
  if (iascii) {
    PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE ADS preconditioning\n"));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    subspace iterations per application %" PetscInt_FMT "\n", jac->as_max_iter));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    subspace cycle type %" PetscInt_FMT "\n", jac->ads_cycle_type));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    subspace iteration tolerance %g\n", (double)jac->as_tol));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    smoother type %" PetscInt_FMT "\n", jac->as_relax_type));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    number of smoothing steps %" PetscInt_FMT "\n", jac->as_relax_times));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    smoother weight %g\n", (double)jac->as_relax_weight));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    smoother omega %g\n", (double)jac->as_omega));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    AMS solver using boomerAMG\n"));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        subspace cycle type %" PetscInt_FMT "\n", jac->ams_cycle_type));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        coarsening type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[0]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[1]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        relaxation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[2]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        interpolation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[3]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[4]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        strength threshold %g\n", (double)jac->as_amg_alpha_theta));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    vector Poisson solver using boomerAMG\n"));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        coarsening type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[0]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_beta_opts[1]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        relaxation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[2]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        interpolation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[3]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_beta_opts[4]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "        strength threshold %g\n", (double)jac->as_amg_beta_theta));
  }
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)G, MATHYPRE, &ishypre));
  if (ishypre) {
    PetscCall(PetscObjectReference((PetscObject)G));
    PetscCall(MatDestroy(&jac->G));
    jac->G = G;
  } else {
    PetscCall(MatDestroy(&jac->G));
    PetscCall(MatConvert(G, MATHYPRE, MAT_INITIAL_MATRIX, &jac->G));
  }
  PetscFunctionReturn(0);
}

/*@
 PCHYPRESetDiscreteGradient - Set discrete gradient matrix

   Collective on PC

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

   Level: intermediate

   Notes:
    G should have as many rows as the number of edges and as many columns as the number of vertices in the mesh

    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

.seealso: `PCHYPRESetDiscreteCurl()`
@*/
PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G) {
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidHeaderSpecific(G, MAT_CLASSID, 2);
  PetscCheckSameComm(pc, 1, G, 2);
  PetscTryMethod(pc, "PCHYPRESetDiscreteGradient_C", (PC, Mat), (pc, G));
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)C, MATHYPRE, &ishypre));
  if (ishypre) {
    PetscCall(PetscObjectReference((PetscObject)C));
    PetscCall(MatDestroy(&jac->C));
    jac->C = C;
  } else {
    PetscCall(MatDestroy(&jac->C));
    PetscCall(MatConvert(C, MATHYPRE, MAT_INITIAL_MATRIX, &jac->C));
  }
  PetscFunctionReturn(0);
}

/*@
 PCHYPRESetDiscreteCurl - Set discrete curl matrix

   Collective on PC

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

   Level: intermediate

   Notes:
    C should have as many rows as the number of faces and as many columns as the number of edges in the mesh

    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

.seealso: `PCHYPRESetDiscreteGradient()`
@*/
PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C) {
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidHeaderSpecific(C, MAT_CLASSID, 2);
  PetscCheckSameComm(pc, 1, C, 2);
  PetscTryMethod(pc, "PCHYPRESetDiscreteCurl_C", (PC, Mat), (pc, C));
  PetscFunctionReturn(0);
}

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

  PetscCall(MatDestroy(&jac->RT_PiFull));
  PetscCall(MatDestroy(&jac->ND_PiFull));
  for (i = 0; i < 3; ++i) {
    PetscCall(MatDestroy(&jac->RT_Pi[i]));
    PetscCall(MatDestroy(&jac->ND_Pi[i]));
  }

  jac->dim = dim;
  if (RT_PiFull) {
    PetscCall(PetscObjectTypeCompare((PetscObject)RT_PiFull, MATHYPRE, &ishypre));
    if (ishypre) {
      PetscCall(PetscObjectReference((PetscObject)RT_PiFull));
      jac->RT_PiFull = RT_PiFull;
    } else {
      PetscCall(MatConvert(RT_PiFull, MATHYPRE, MAT_INITIAL_MATRIX, &jac->RT_PiFull));
    }
  }
  if (RT_Pi) {
    for (i = 0; i < dim; ++i) {
      if (RT_Pi[i]) {
        PetscCall(PetscObjectTypeCompare((PetscObject)RT_Pi[i], MATHYPRE, &ishypre));
        if (ishypre) {
          PetscCall(PetscObjectReference((PetscObject)RT_Pi[i]));
          jac->RT_Pi[i] = RT_Pi[i];
        } else {
          PetscCall(MatConvert(RT_Pi[i], MATHYPRE, MAT_INITIAL_MATRIX, &jac->RT_Pi[i]));
        }
      }
    }
  }
  if (ND_PiFull) {
    PetscCall(PetscObjectTypeCompare((PetscObject)ND_PiFull, MATHYPRE, &ishypre));
    if (ishypre) {
      PetscCall(PetscObjectReference((PetscObject)ND_PiFull));
      jac->ND_PiFull = ND_PiFull;
    } else {
      PetscCall(MatConvert(ND_PiFull, MATHYPRE, MAT_INITIAL_MATRIX, &jac->ND_PiFull));
    }
  }
  if (ND_Pi) {
    for (i = 0; i < dim; ++i) {
      if (ND_Pi[i]) {
        PetscCall(PetscObjectTypeCompare((PetscObject)ND_Pi[i], MATHYPRE, &ishypre));
        if (ishypre) {
          PetscCall(PetscObjectReference((PetscObject)ND_Pi[i]));
          jac->ND_Pi[i] = ND_Pi[i];
        } else {
          PetscCall(MatConvert(ND_Pi[i], MATHYPRE, MAT_INITIAL_MATRIX, &jac->ND_Pi[i]));
        }
      }
    }
  }

  PetscFunctionReturn(0);
}

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

   Collective on PC

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

   Notes:
    For AMS, only Nedelec interpolation matrices are needed, the Raviart-Thomas interpolation matrices can be set to NULL.

    For ADS, both type of interpolation matrices are needed.

   Level: intermediate

@*/
PetscErrorCode PCHYPRESetInterpolations(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[]) {
  PetscInt i;

  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  if (RT_PiFull) {
    PetscValidHeaderSpecific(RT_PiFull, MAT_CLASSID, 3);
    PetscCheckSameComm(pc, 1, RT_PiFull, 3);
  }
  if (RT_Pi) {
    PetscValidPointer(RT_Pi, 4);
    for (i = 0; i < dim; ++i) {
      if (RT_Pi[i]) {
        PetscValidHeaderSpecific(RT_Pi[i], MAT_CLASSID, 4);
        PetscCheckSameComm(pc, 1, RT_Pi[i], 4);
      }
    }
  }
  if (ND_PiFull) {
    PetscValidHeaderSpecific(ND_PiFull, MAT_CLASSID, 5);
    PetscCheckSameComm(pc, 1, ND_PiFull, 5);
  }
  if (ND_Pi) {
    PetscValidPointer(ND_Pi, 6);
    for (i = 0; i < dim; ++i) {
      if (ND_Pi[i]) {
        PetscValidHeaderSpecific(ND_Pi[i], MAT_CLASSID, 6);
        PetscCheckSameComm(pc, 1, ND_Pi[i], 6);
      }
    }
  }
  PetscTryMethod(pc, "PCHYPRESetInterpolations_C", (PC, PetscInt, Mat, Mat[], Mat, Mat[]), (pc, dim, RT_PiFull, RT_Pi, ND_PiFull, ND_Pi));
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre));
  if (ishypre) {
    if (isalpha) {
      PetscCall(PetscObjectReference((PetscObject)A));
      PetscCall(MatDestroy(&jac->alpha_Poisson));
      jac->alpha_Poisson = A;
    } else {
      if (A) {
        PetscCall(PetscObjectReference((PetscObject)A));
      } else {
        jac->ams_beta_is_zero = PETSC_TRUE;
      }
      PetscCall(MatDestroy(&jac->beta_Poisson));
      jac->beta_Poisson = A;
    }
  } else {
    if (isalpha) {
      PetscCall(MatDestroy(&jac->alpha_Poisson));
      PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &jac->alpha_Poisson));
    } else {
      if (A) {
        PetscCall(MatDestroy(&jac->beta_Poisson));
        PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &jac->beta_Poisson));
      } else {
        PetscCall(MatDestroy(&jac->beta_Poisson));
        jac->ams_beta_is_zero = PETSC_TRUE;
      }
    }
  }
  PetscFunctionReturn(0);
}

/*@
 PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix

   Collective on PC

   Input Parameters:
+  pc - the preconditioning context
-  A - the matrix

   Level: intermediate

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

.seealso: `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetBetaPoissonMatrix()`
@*/
PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A) {
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidHeaderSpecific(A, MAT_CLASSID, 2);
  PetscCheckSameComm(pc, 1, A, 2);
  PetscTryMethod(pc, "PCHYPRESetPoissonMatrix_C", (PC, Mat, PetscBool), (pc, A, PETSC_TRUE));
  PetscFunctionReturn(0);
}

/*@
 PCHYPRESetBetaPoissonMatrix - Set Poisson matrix

   Collective on PC

   Input Parameters:
+  pc - the preconditioning context
-  A - the matrix

   Level: intermediate

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

.seealso: `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetAlphaPoissonMatrix()`
@*/
PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A) {
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  if (A) {
    PetscValidHeaderSpecific(A, MAT_CLASSID, 2);
    PetscCheckSameComm(pc, 1, A, 2);
  }
  PetscTryMethod(pc, "PCHYPRESetPoissonMatrix_C", (PC, Mat, PetscBool), (pc, A, PETSC_FALSE));
  PetscFunctionReturn(0);
}

static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE(PC pc, Vec ozz, Vec zoz, Vec zzo) {
  PC_HYPRE *jac = (PC_HYPRE *)pc->data;

  PetscFunctionBegin;
  /* throw away any vector if already set */
  PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[0]));
  PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[1]));
  PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[2]));
  PetscCall(VecHYPRE_IJVectorCreate(ozz->map, &jac->constants[0]));
  PetscCall(VecHYPRE_IJVectorCopy(ozz, jac->constants[0]));
  PetscCall(VecHYPRE_IJVectorCreate(zoz->map, &jac->constants[1]));
  PetscCall(VecHYPRE_IJVectorCopy(zoz, jac->constants[1]));
  jac->dim = 2;
  if (zzo) {
    PetscCall(VecHYPRE_IJVectorCreate(zzo->map, &jac->constants[2]));
    PetscCall(VecHYPRE_IJVectorCopy(zzo, jac->constants[2]));
    jac->dim++;
  }
  PetscFunctionReturn(0);
}

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

   Collective on PC

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

   Level: intermediate

@*/
PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo) {
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidHeaderSpecific(ozz, VEC_CLASSID, 2);
  PetscValidHeaderSpecific(zoz, VEC_CLASSID, 3);
  if (zzo) PetscValidHeaderSpecific(zzo, VEC_CLASSID, 4);
  PetscCheckSameComm(pc, 1, ozz, 2);
  PetscCheckSameComm(pc, 1, zoz, 3);
  if (zzo) PetscCheckSameComm(pc, 1, zzo, 4);
  PetscTryMethod(pc, "PCHYPRESetEdgeConstantVectors_C", (PC, Vec, Vec, Vec), (pc, ozz, zoz, zzo));
  PetscFunctionReturn(0);
}

static PetscErrorCode PCHYPREAMSSetInteriorNodes_HYPRE(PC pc, Vec interior) {
  PC_HYPRE *jac = (PC_HYPRE *)pc->data;

  PetscFunctionBegin;
  PetscCall(VecHYPRE_IJVectorDestroy(&jac->interior));
  PetscCall(VecHYPRE_IJVectorCreate(interior->map, &jac->interior));
  PetscCall(VecHYPRE_IJVectorCopy(interior, jac->interior));
  jac->ams_beta_is_zero_part = PETSC_TRUE;
  PetscFunctionReturn(0);
}

/*@
 PCHYPREAMSSetInteriorNodes - Set the list of interior nodes to a zero-conductivity region.

   Collective on PC

   Input Parameters:
+  pc - the preconditioning context
-  interior - vector. node is interior if its entry in the array is 1.0.

   Level: intermediate

   Note:
   This calls HYPRE_AMSSetInteriorNodes()
.seealso:
@*/
PetscErrorCode PCHYPREAMSSetInteriorNodes(PC pc, Vec interior) {
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidHeaderSpecific(interior, VEC_CLASSID, 2);
  PetscCheckSameComm(pc, 1, interior, 2);
  PetscTryMethod(pc, "PCHYPREAMSSetInteriorNodes_C", (PC, Vec), (pc, interior));
  PetscFunctionReturn(0);
}

static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords) {
  PC_HYPRE *jac = (PC_HYPRE *)pc->data;
  Vec       tv;
  PetscInt  i;

  PetscFunctionBegin;
  /* throw away any coordinate vector if already set */
  PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[0]));
  PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[1]));
  PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[2]));
  jac->dim = dim;

  /* compute IJ vector for coordinates */
  PetscCall(VecCreate(PetscObjectComm((PetscObject)pc), &tv));
  PetscCall(VecSetType(tv, VECSTANDARD));
  PetscCall(VecSetSizes(tv, nloc, PETSC_DECIDE));
  for (i = 0; i < dim; i++) {
    PetscScalar *array;
    PetscInt     j;

    PetscCall(VecHYPRE_IJVectorCreate(tv->map, &jac->coords[i]));
    PetscCall(VecGetArrayWrite(tv, &array));
    for (j = 0; j < nloc; j++) array[j] = coords[j * dim + i];
    PetscCall(VecRestoreArrayWrite(tv, &array));
    PetscCall(VecHYPRE_IJVectorCopy(tv, jac->coords[i]));
  }
  PetscCall(VecDestroy(&tv));
  PetscFunctionReturn(0);
}

/* ---------------------------------------------------------------------------------*/

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

  PetscFunctionBegin;
  *name = jac->hypre_type;
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  if (jac->hypre_type) {
    PetscCall(PetscStrcmp(jac->hypre_type, name, &flag));
    PetscCheck(flag, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Cannot reset the HYPRE preconditioner type once it has been set");
    PetscFunctionReturn(0);
  } else {
    PetscCall(PetscStrallocpy(name, &jac->hypre_type));
  }

  jac->maxiter         = PETSC_DEFAULT;
  jac->tol             = PETSC_DEFAULT;
  jac->printstatistics = PetscLogPrintInfo;

  PetscCall(PetscStrcmp("pilut", jac->hypre_type, &flag));
  if (flag) {
    PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
    PetscCallExternal(HYPRE_ParCSRPilutCreate, jac->comm_hypre, &jac->hsolver);
    pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
    pc->ops->view           = PCView_HYPRE_Pilut;
    jac->destroy            = HYPRE_ParCSRPilutDestroy;
    jac->setup              = HYPRE_ParCSRPilutSetup;
    jac->solve              = HYPRE_ParCSRPilutSolve;
    jac->factorrowsize      = PETSC_DEFAULT;
    PetscFunctionReturn(0);
  }
  PetscCall(PetscStrcmp("euclid", jac->hypre_type, &flag));
  if (flag) {
#if defined(PETSC_USE_64BIT_INDICES)
    SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Hypre Euclid does not support 64 bit indices");
#endif
    PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
    PetscCallExternal(HYPRE_EuclidCreate, jac->comm_hypre, &jac->hsolver);
    pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
    pc->ops->view           = PCView_HYPRE_Euclid;
    jac->destroy            = HYPRE_EuclidDestroy;
    jac->setup              = HYPRE_EuclidSetup;
    jac->solve              = HYPRE_EuclidSolve;
    jac->factorrowsize      = PETSC_DEFAULT;
    jac->eu_level           = PETSC_DEFAULT; /* default */
    PetscFunctionReturn(0);
  }
  PetscCall(PetscStrcmp("parasails", jac->hypre_type, &flag));
  if (flag) {
    PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
    PetscCallExternal(HYPRE_ParaSailsCreate, jac->comm_hypre, &jac->hsolver);
    pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
    pc->ops->view           = PCView_HYPRE_ParaSails;
    jac->destroy            = HYPRE_ParaSailsDestroy;
    jac->setup              = HYPRE_ParaSailsSetup;
    jac->solve              = HYPRE_ParaSailsSolve;
    /* initialize */
    jac->nlevels            = 1;
    jac->threshold          = .1;
    jac->filter             = .1;
    jac->loadbal            = 0;
    if (PetscLogPrintInfo) jac->logging = (int)PETSC_TRUE;
    else jac->logging = (int)PETSC_FALSE;

    jac->ruse = (int)PETSC_FALSE;
    jac->symt = 0;
    PetscCallExternal(HYPRE_ParaSailsSetParams, jac->hsolver, jac->threshold, jac->nlevels);
    PetscCallExternal(HYPRE_ParaSailsSetFilter, jac->hsolver, jac->filter);
    PetscCallExternal(HYPRE_ParaSailsSetLoadbal, jac->hsolver, jac->loadbal);
    PetscCallExternal(HYPRE_ParaSailsSetLogging, jac->hsolver, jac->logging);
    PetscCallExternal(HYPRE_ParaSailsSetReuse, jac->hsolver, jac->ruse);
    PetscCallExternal(HYPRE_ParaSailsSetSym, jac->hsolver, jac->symt);
    PetscFunctionReturn(0);
  }
  PetscCall(PetscStrcmp("boomeramg", jac->hypre_type, &flag));
  if (flag) {
    PetscCallExternal(HYPRE_BoomerAMGCreate, &jac->hsolver);
    pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_BoomerAMG;
    pc->ops->view            = PCView_HYPRE_BoomerAMG;
    pc->ops->applytranspose  = PCApplyTranspose_HYPRE_BoomerAMG;
    pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
    PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", PCGetInterpolations_BoomerAMG));
    PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", PCGetCoarseOperators_BoomerAMG));
    jac->destroy         = HYPRE_BoomerAMGDestroy;
    jac->setup           = HYPRE_BoomerAMGSetup;
    jac->solve           = HYPRE_BoomerAMGSolve;
    jac->applyrichardson = PETSC_FALSE;
    /* these defaults match the hypre defaults */
    jac->cycletype       = 1;
    jac->maxlevels       = 25;
    jac->maxiter         = 1;
    jac->tol             = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
    jac->truncfactor     = 0.0;
    jac->strongthreshold = .25;
    jac->maxrowsum       = .9;
    jac->coarsentype     = 6;
    jac->measuretype     = 0;
    jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
    jac->smoothtype                                              = -1; /* Not set by default */
    jac->smoothnumlevels                                         = 25;
    jac->eu_level                                                = 0;
    jac->eu_droptolerance                                        = 0;
    jac->eu_bj                                                   = 0;
    jac->relaxtype[0] = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a PC - most likely with CG */
    jac->relaxtype[2]                     = 9; /*G.E. */
    jac->relaxweight                      = 1.0;
    jac->outerrelaxweight                 = 1.0;
    jac->relaxorder                       = 1;
    jac->interptype                       = 0;
    jac->Rtype                            = 0;
    jac->Rstrongthreshold                 = 0.25;
    jac->Rfilterthreshold                 = 0.0;
    jac->Adroptype                        = -1;
    jac->Adroptol                         = 0.0;
    jac->agg_nl                           = 0;
    jac->agg_interptype                   = 4;
    jac->pmax                             = 0;
    jac->truncfactor                      = 0.0;
    jac->agg_num_paths                    = 1;
    jac->maxc                             = 9;
    jac->minc                             = 1;
    jac->nodal_coarsening                 = 0;
    jac->nodal_coarsening_diag            = 0;
    jac->vec_interp_variant               = 0;
    jac->vec_interp_qmax                  = 0;
    jac->vec_interp_smooth                = PETSC_FALSE;
    jac->interp_refine                    = 0;
    jac->nodal_relax                      = PETSC_FALSE;
    jac->nodal_relax_levels               = 1;
    jac->rap2                             = 0;

    /* GPU defaults
         from https://hypre.readthedocs.io/en/latest/solvers-boomeramg.html#gpu-supported-options
         and /src/parcsr_ls/par_amg.c */
#if defined(PETSC_HAVE_HYPRE_DEVICE)
    jac->keeptranspose  = PETSC_TRUE;
    jac->mod_rap2       = 1;
    jac->coarsentype    = 8;
    jac->relaxorder     = 0;
    jac->interptype     = 6;
    jac->relaxtype[0]   = 18;
    jac->relaxtype[1]   = 18;
    jac->agg_interptype = 7;
#else
    jac->keeptranspose = PETSC_FALSE;
    jac->mod_rap2      = 0;
#endif
    PetscCallExternal(HYPRE_BoomerAMGSetCycleType, jac->hsolver, jac->cycletype);
    PetscCallExternal(HYPRE_BoomerAMGSetMaxLevels, jac->hsolver, jac->maxlevels);
    PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, jac->maxiter);
    PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, jac->tol);
    PetscCallExternal(HYPRE_BoomerAMGSetTruncFactor, jac->hsolver, jac->truncfactor);
    PetscCallExternal(HYPRE_BoomerAMGSetStrongThreshold, jac->hsolver, jac->strongthreshold);
    PetscCallExternal(HYPRE_BoomerAMGSetMaxRowSum, jac->hsolver, jac->maxrowsum);
    PetscCallExternal(HYPRE_BoomerAMGSetCoarsenType, jac->hsolver, jac->coarsentype);
    PetscCallExternal(HYPRE_BoomerAMGSetMeasureType, jac->hsolver, jac->measuretype);
    PetscCallExternal(HYPRE_BoomerAMGSetRelaxOrder, jac->hsolver, jac->relaxorder);
    PetscCallExternal(HYPRE_BoomerAMGSetInterpType, jac->hsolver, jac->interptype);
    PetscCallExternal(HYPRE_BoomerAMGSetAggNumLevels, jac->hsolver, jac->agg_nl);
    PetscCallExternal(HYPRE_BoomerAMGSetAggInterpType, jac->hsolver, jac->agg_interptype);
    PetscCallExternal(HYPRE_BoomerAMGSetPMaxElmts, jac->hsolver, jac->pmax);
    PetscCallExternal(HYPRE_BoomerAMGSetNumPaths, jac->hsolver, jac->agg_num_paths);
    PetscCallExternal(HYPRE_BoomerAMGSetRelaxType, jac->hsolver, jac->relaxtype[0]);  /* defaults coarse to 9 */
    PetscCallExternal(HYPRE_BoomerAMGSetNumSweeps, jac->hsolver, jac->gridsweeps[0]); /* defaults coarse to 1 */
    PetscCallExternal(HYPRE_BoomerAMGSetMaxCoarseSize, jac->hsolver, jac->maxc);
    PetscCallExternal(HYPRE_BoomerAMGSetMinCoarseSize, jac->hsolver, jac->minc);
    /* GPU */
#if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
    PetscCallExternal(HYPRE_BoomerAMGSetKeepTranspose, jac->hsolver, jac->keeptranspose ? 1 : 0);
    PetscCallExternal(HYPRE_BoomerAMGSetRAP2, jac->hsolver, jac->rap2);
    PetscCallExternal(HYPRE_BoomerAMGSetModuleRAP2, jac->hsolver, jac->mod_rap2);
#endif

    /* AIR */
#if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
    PetscCallExternal(HYPRE_BoomerAMGSetRestriction, jac->hsolver, jac->Rtype);
    PetscCallExternal(HYPRE_BoomerAMGSetStrongThresholdR, jac->hsolver, jac->Rstrongthreshold);
    PetscCallExternal(HYPRE_BoomerAMGSetFilterThresholdR, jac->hsolver, jac->Rfilterthreshold);
    PetscCallExternal(HYPRE_BoomerAMGSetADropTol, jac->hsolver, jac->Adroptol);
    PetscCallExternal(HYPRE_BoomerAMGSetADropType, jac->hsolver, jac->Adroptype);
#endif
    PetscFunctionReturn(0);
  }
  PetscCall(PetscStrcmp("ams", jac->hypre_type, &flag));
  if (flag) {
    PetscCall(HYPRE_AMSCreate(&jac->hsolver));
    pc->ops->setfromoptions   = PCSetFromOptions_HYPRE_AMS;
    pc->ops->view             = PCView_HYPRE_AMS;
    jac->destroy              = HYPRE_AMSDestroy;
    jac->setup                = HYPRE_AMSSetup;
    jac->solve                = HYPRE_AMSSolve;
    jac->coords[0]            = NULL;
    jac->coords[1]            = NULL;
    jac->coords[2]            = NULL;
    jac->interior             = NULL;
    /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */
    jac->as_print             = 0;
    jac->as_max_iter          = 1;  /* used as a preconditioner */
    jac->as_tol               = 0.; /* used as a preconditioner */
    jac->ams_cycle_type       = 13;
    /* Smoothing options */
    jac->as_relax_type        = 2;
    jac->as_relax_times       = 1;
    jac->as_relax_weight      = 1.0;
    jac->as_omega             = 1.0;
    /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
    jac->as_amg_alpha_opts[0] = 10;
    jac->as_amg_alpha_opts[1] = 1;
    jac->as_amg_alpha_opts[2] = 6;
    jac->as_amg_alpha_opts[3] = 6;
    jac->as_amg_alpha_opts[4] = 4;
    jac->as_amg_alpha_theta   = 0.25;
    /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
    jac->as_amg_beta_opts[0]  = 10;
    jac->as_amg_beta_opts[1]  = 1;
    jac->as_amg_beta_opts[2]  = 6;
    jac->as_amg_beta_opts[3]  = 6;
    jac->as_amg_beta_opts[4]  = 4;
    jac->as_amg_beta_theta    = 0.25;
    PetscCallExternal(HYPRE_AMSSetPrintLevel, jac->hsolver, jac->as_print);
    PetscCallExternal(HYPRE_AMSSetMaxIter, jac->hsolver, jac->as_max_iter);
    PetscCallExternal(HYPRE_AMSSetCycleType, jac->hsolver, jac->ams_cycle_type);
    PetscCallExternal(HYPRE_AMSSetTol, jac->hsolver, jac->as_tol);
    PetscCallExternal(HYPRE_AMSSetSmoothingOptions, jac->hsolver, jac->as_relax_type, jac->as_relax_times, jac->as_relax_weight, jac->as_omega);
    PetscCallExternal(HYPRE_AMSSetAlphaAMGOptions, jac->hsolver, jac->as_amg_alpha_opts[0], /* AMG coarsen type */
                      jac->as_amg_alpha_opts[1],                                            /* AMG agg_levels */
                      jac->as_amg_alpha_opts[2],                                            /* AMG relax_type */
                      jac->as_amg_alpha_theta, jac->as_amg_alpha_opts[3],                   /* AMG interp_type */
                      jac->as_amg_alpha_opts[4]);                                           /* AMG Pmax */
    PetscCallExternal(HYPRE_AMSSetBetaAMGOptions, jac->hsolver, jac->as_amg_beta_opts[0],   /* AMG coarsen type */
                      jac->as_amg_beta_opts[1],                                             /* AMG agg_levels */
                      jac->as_amg_beta_opts[2],                                             /* AMG relax_type */
                      jac->as_amg_beta_theta, jac->as_amg_beta_opts[3],                     /* AMG interp_type */
                      jac->as_amg_beta_opts[4]);                                            /* AMG Pmax */
    /* Zero conductivity */
    jac->ams_beta_is_zero      = PETSC_FALSE;
    jac->ams_beta_is_zero_part = PETSC_FALSE;
    PetscFunctionReturn(0);
  }
  PetscCall(PetscStrcmp("ads", jac->hypre_type, &flag));
  if (flag) {
    PetscCall(HYPRE_ADSCreate(&jac->hsolver));
    pc->ops->setfromoptions   = PCSetFromOptions_HYPRE_ADS;
    pc->ops->view             = PCView_HYPRE_ADS;
    jac->destroy              = HYPRE_ADSDestroy;
    jac->setup                = HYPRE_ADSSetup;
    jac->solve                = HYPRE_ADSSolve;
    jac->coords[0]            = NULL;
    jac->coords[1]            = NULL;
    jac->coords[2]            = NULL;
    /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */
    jac->as_print             = 0;
    jac->as_max_iter          = 1;  /* used as a preconditioner */
    jac->as_tol               = 0.; /* used as a preconditioner */
    jac->ads_cycle_type       = 13;
    /* Smoothing options */
    jac->as_relax_type        = 2;
    jac->as_relax_times       = 1;
    jac->as_relax_weight      = 1.0;
    jac->as_omega             = 1.0;
    /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */
    jac->ams_cycle_type       = 14;
    jac->as_amg_alpha_opts[0] = 10;
    jac->as_amg_alpha_opts[1] = 1;
    jac->as_amg_alpha_opts[2] = 6;
    jac->as_amg_alpha_opts[3] = 6;
    jac->as_amg_alpha_opts[4] = 4;
    jac->as_amg_alpha_theta   = 0.25;
    /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
    jac->as_amg_beta_opts[0]  = 10;
    jac->as_amg_beta_opts[1]  = 1;
    jac->as_amg_beta_opts[2]  = 6;
    jac->as_amg_beta_opts[3]  = 6;
    jac->as_amg_beta_opts[4]  = 4;
    jac->as_amg_beta_theta    = 0.25;
    PetscCallExternal(HYPRE_ADSSetPrintLevel, jac->hsolver, jac->as_print);
    PetscCallExternal(HYPRE_ADSSetMaxIter, jac->hsolver, jac->as_max_iter);
    PetscCallExternal(HYPRE_ADSSetCycleType, jac->hsolver, jac->ams_cycle_type);
    PetscCallExternal(HYPRE_ADSSetTol, jac->hsolver, jac->as_tol);
    PetscCallExternal(HYPRE_ADSSetSmoothingOptions, jac->hsolver, jac->as_relax_type, jac->as_relax_times, jac->as_relax_weight, jac->as_omega);
    PetscCallExternal(HYPRE_ADSSetAMSOptions, jac->hsolver, jac->ams_cycle_type,      /* AMG coarsen type */
                      jac->as_amg_alpha_opts[0],                                      /* AMG coarsen type */
                      jac->as_amg_alpha_opts[1],                                      /* AMG agg_levels */
                      jac->as_amg_alpha_opts[2],                                      /* AMG relax_type */
                      jac->as_amg_alpha_theta, jac->as_amg_alpha_opts[3],             /* AMG interp_type */
                      jac->as_amg_alpha_opts[4]);                                     /* AMG Pmax */
    PetscCallExternal(HYPRE_ADSSetAMGOptions, jac->hsolver, jac->as_amg_beta_opts[0], /* AMG coarsen type */
                      jac->as_amg_beta_opts[1],                                       /* AMG agg_levels */
                      jac->as_amg_beta_opts[2],                                       /* AMG relax_type */
                      jac->as_amg_beta_theta, jac->as_amg_beta_opts[3],               /* AMG interp_type */
                      jac->as_amg_beta_opts[4]);                                      /* AMG Pmax */
    PetscFunctionReturn(0);
  }
  PetscCall(PetscFree(jac->hypre_type));

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

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

  PetscFunctionBegin;
  PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE preconditioner options");
  PetscCall(PetscOptionsEList("-pc_hypre_type", "HYPRE preconditioner type", "PCHYPRESetType", type, PETSC_STATIC_ARRAY_LENGTH(type), "boomeramg", &indx, &flg));
  if (flg) {
    PetscCall(PCHYPRESetType_HYPRE(pc, type[indx]));
  } else {
    PetscCall(PCHYPRESetType_HYPRE(pc, "boomeramg"));
  }
  PetscTryTypeMethod(pc, setfromoptions, PetscOptionsObject);
  PetscOptionsHeadEnd();
  PetscFunctionReturn(0);
}

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

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

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

   Level: intermediate

.seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
          `PCHYPRE`

@*/
PetscErrorCode PCHYPRESetType(PC pc, const char name[]) {
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidCharPointer(name, 2);
  PetscTryMethod(pc, "PCHYPRESetType_C", (PC, const char[]), (pc, name));
  PetscFunctionReturn(0);
}

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

   Input Parameter:
.     pc - the preconditioner context

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

   Level: intermediate

.seealso: `PCCreate()`, `PCHYPRESetType()`, `PCType`, `PC`,
          `PCHYPRE`

@*/
PetscErrorCode PCHYPREGetType(PC pc, const char *name[]) {
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscValidPointer(name, 2);
  PetscTryMethod(pc, "PCHYPREGetType_C", (PC, const char *[]), (pc, name));
  PetscFunctionReturn(0);
}

/*@C
   PCMGGalerkinSetMatProductAlgorithm - Set type of SpGEMM for hypre to use

   Logically Collective on PC

   Input Parameters:
+  pc - the hypre context
-  type - one of 'cusparse', 'hypre'

   Options Database Key:
.  -pc_mg_galerkin_mat_product_algorithm <cusparse,hypre> - Type of SpGEMM to use in hypre

   Level: intermediate

.seealso: `PCMGGalerkinGetMatProductAlgorithm()`

@*/
PetscErrorCode PCMGGalerkinSetMatProductAlgorithm(PC pc, const char name[]) {
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscTryMethod(pc, "PCMGGalerkinSetMatProductAlgorithm_C", (PC, const char[]), (pc, name));
  PetscFunctionReturn(0);
}

/*@C
   PCMGGalerkinGetMatProductAlgorithm - Get type of SpGEMM for hypre

   Not Collective

   Input Parameter:
.  pc - the multigrid context

   Output Parameter:
.  name - one of 'cusparse', 'hypre'

   Level: intermediate

.seealso: `PCMGGalerkinSetMatProductAlgorithm()`

@*/
PetscErrorCode PCMGGalerkinGetMatProductAlgorithm(PC pc, const char *name[]) {
  PetscFunctionBegin;
  PetscValidHeaderSpecific(pc, PC_CLASSID, 1);
  PetscTryMethod(pc, "PCMGGalerkinGetMatProductAlgorithm_C", (PC, const char *[]), (pc, name));
  PetscFunctionReturn(0);
}

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

   Options Database Keys:
+   -pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads
.   -pc_hypre_boomeramg_nodal_coarsen <n> - where n is from 1 to 6 (see HYPRE_BOOMERAMGSetNodal())
.   -pc_hypre_boomeramg_vec_interp_variant <v> - where v is from 1 to 3 (see HYPRE_BoomerAMGSetInterpVecVariant())
-   Many others, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX preconditioner

   Level: intermediate

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

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

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

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

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

          See PCPFMG for access to the hypre Struct PFMG solver

   GPU Notes:
     To configure hypre BoomerAMG so that it can utilize NVIDIA GPUs run ./configure --download-hypre --with-cuda
     Then pass VECCUDA vectors and MATAIJCUSPARSE matrices to the solvers and PETSc will automatically utilize hypre's GPU solvers.

     To configure hypre BoomerAMG so that it can utilize AMD GPUs run ./configure --download-hypre --with-hip
     Then pass VECHIP vectors to the solvers and PETSc will automatically utilize hypre's GPU solvers.

.seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`,
          `PCHYPRESetType()`, `PCPFMG`

M*/

PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc) {
  PC_HYPRE *jac;

  PetscFunctionBegin;
  PetscCall(PetscNewLog(pc, &jac));

  pc->data                = jac;
  pc->ops->reset          = PCReset_HYPRE;
  pc->ops->destroy        = PCDestroy_HYPRE;
  pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
  pc->ops->setup          = PCSetUp_HYPRE;
  pc->ops->apply          = PCApply_HYPRE;
  jac->comm_hypre         = MPI_COMM_NULL;
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetType_C", PCHYPRESetType_HYPRE));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetType_C", PCHYPREGetType_HYPRE));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_HYPRE));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteGradient_C", PCHYPRESetDiscreteGradient_HYPRE));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteCurl_C", PCHYPRESetDiscreteCurl_HYPRE));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetInterpolations_C", PCHYPRESetInterpolations_HYPRE));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetEdgeConstantVectors_C", PCHYPRESetEdgeConstantVectors_HYPRE));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREAMSSetInteriorNodes_C", PCHYPREAMSSetInteriorNodes_HYPRE));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetPoissonMatrix_C", PCHYPRESetPoissonMatrix_HYPRE));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinSetMatProductAlgorithm_C", PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG));
  PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinGetMatProductAlgorithm_C", PCMGGalerkinGetMatProductAlgorithm_HYPRE_BoomerAMG));
#if defined(PETSC_HAVE_HYPRE_DEVICE)
#if defined(HYPRE_USING_HIP)
  PetscCall(PetscDeviceInitialize(PETSC_DEVICE_HIP));
#endif
#if defined(HYPRE_USING_CUDA)
  PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA));
#endif
#endif
  PetscFunctionReturn(0);
}

/* ---------------------------------------------------------------------------------------------------------------------------------*/

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

  /* keep copy of PFMG options used so may view them */
  PetscInt  its;
  double    tol;
  PetscInt  relax_type;
  PetscInt  rap_type;
  PetscInt  num_pre_relax, num_post_relax;
  PetscInt  max_levels;
  PetscInt  skip_relax;
  PetscBool print_statistics;
} PC_PFMG;

PetscErrorCode PCDestroy_PFMG(PC pc) {
  PC_PFMG *ex = (PC_PFMG *)pc->data;

  PetscFunctionBegin;
  if (ex->hsolver) PetscCallExternal(HYPRE_StructPFMGDestroy, ex->hsolver);
  PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
  PetscCall(PetscFree(pc->data));
  PetscFunctionReturn(0);
}

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

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

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
  if (iascii) {
    PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE PFMG preconditioning\n"));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    max iterations %" PetscInt_FMT "\n", ex->its));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    tolerance %g\n", ex->tol));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    relax type %s\n", PFMGRelaxType[ex->relax_type]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    RAP type %s\n", PFMGRAPType[ex->rap_type]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    max levels %" PetscInt_FMT "\n", ex->max_levels));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    skip relax %" PetscInt_FMT "\n", ex->skip_relax));
  }
  PetscFunctionReturn(0);
}

PetscErrorCode PCSetFromOptions_PFMG(PC pc, PetscOptionItems *PetscOptionsObject) {
  PC_PFMG *ex = (PC_PFMG *)pc->data;

  PetscFunctionBegin;
  PetscOptionsHeadBegin(PetscOptionsObject, "PFMG options");
  PetscCall(PetscOptionsBool("-pc_pfmg_print_statistics", "Print statistics", "HYPRE_StructPFMGSetPrintLevel", ex->print_statistics, &ex->print_statistics, NULL));
  PetscCall(PetscOptionsInt("-pc_pfmg_its", "Number of iterations of PFMG to use as preconditioner", "HYPRE_StructPFMGSetMaxIter", ex->its, &ex->its, NULL));
  PetscCallExternal(HYPRE_StructPFMGSetMaxIter, ex->hsolver, ex->its);
  PetscCall(PetscOptionsInt("-pc_pfmg_num_pre_relax", "Number of smoothing steps before coarse grid", "HYPRE_StructPFMGSetNumPreRelax", ex->num_pre_relax, &ex->num_pre_relax, NULL));
  PetscCallExternal(HYPRE_StructPFMGSetNumPreRelax, ex->hsolver, ex->num_pre_relax);
  PetscCall(PetscOptionsInt("-pc_pfmg_num_post_relax", "Number of smoothing steps after coarse grid", "HYPRE_StructPFMGSetNumPostRelax", ex->num_post_relax, &ex->num_post_relax, NULL));
  PetscCallExternal(HYPRE_StructPFMGSetNumPostRelax, ex->hsolver, ex->num_post_relax);

  PetscCall(PetscOptionsInt("-pc_pfmg_max_levels", "Max Levels for MG hierarchy", "HYPRE_StructPFMGSetMaxLevels", ex->max_levels, &ex->max_levels, NULL));
  PetscCallExternal(HYPRE_StructPFMGSetMaxLevels, ex->hsolver, ex->max_levels);

  PetscCall(PetscOptionsReal("-pc_pfmg_tol", "Tolerance of PFMG", "HYPRE_StructPFMGSetTol", ex->tol, &ex->tol, NULL));
  PetscCallExternal(HYPRE_StructPFMGSetTol, ex->hsolver, ex->tol);
  PetscCall(PetscOptionsEList("-pc_pfmg_relax_type", "Relax type for the up and down cycles", "HYPRE_StructPFMGSetRelaxType", PFMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(PFMGRelaxType), PFMGRelaxType[ex->relax_type], &ex->relax_type, NULL));
  PetscCallExternal(HYPRE_StructPFMGSetRelaxType, ex->hsolver, ex->relax_type);
  PetscCall(PetscOptionsEList("-pc_pfmg_rap_type", "RAP type", "HYPRE_StructPFMGSetRAPType", PFMGRAPType, PETSC_STATIC_ARRAY_LENGTH(PFMGRAPType), PFMGRAPType[ex->rap_type], &ex->rap_type, NULL));
  PetscCallExternal(HYPRE_StructPFMGSetRAPType, ex->hsolver, ex->rap_type);
  PetscCall(PetscOptionsInt("-pc_pfmg_skip_relax", "Skip relaxation on certain grids for isotropic problems. This can greatly improve efficiency by eliminating unnecessary relaxations when the underlying problem is isotropic", "HYPRE_StructPFMGSetSkipRelax", ex->skip_relax, &ex->skip_relax, NULL));
  PetscCallExternal(HYPRE_StructPFMGSetSkipRelax, ex->hsolver, ex->skip_relax);
  PetscOptionsHeadEnd();
  PetscFunctionReturn(0);
}

PetscErrorCode PCApply_PFMG(PC pc, Vec x, Vec y) {
  PC_PFMG           *ex = (PC_PFMG *)pc->data;
  PetscScalar       *yy;
  const PetscScalar *xx;
  PetscInt           ilower[3], iupper[3];
  HYPRE_Int          hlower[3], hupper[3];
  Mat_HYPREStruct   *mx = (Mat_HYPREStruct *)(pc->pmat->data);

  PetscFunctionBegin;
  PetscCall(PetscCitationsRegister(hypreCitation, &cite));
  PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
  /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
  iupper[0] += ilower[0] - 1;
  iupper[1] += ilower[1] - 1;
  iupper[2] += ilower[2] - 1;
  hlower[0] = (HYPRE_Int)ilower[0];
  hlower[1] = (HYPRE_Int)ilower[1];
  hlower[2] = (HYPRE_Int)ilower[2];
  hupper[0] = (HYPRE_Int)iupper[0];
  hupper[1] = (HYPRE_Int)iupper[1];
  hupper[2] = (HYPRE_Int)iupper[2];

  /* copy x values over to hypre */
  PetscCallExternal(HYPRE_StructVectorSetConstantValues, mx->hb, 0.0);
  PetscCall(VecGetArrayRead(x, &xx));
  PetscCallExternal(HYPRE_StructVectorSetBoxValues, mx->hb, hlower, hupper, (HYPRE_Complex *)xx);
  PetscCall(VecRestoreArrayRead(x, &xx));
  PetscCallExternal(HYPRE_StructVectorAssemble, mx->hb);
  PetscCallExternal(HYPRE_StructPFMGSolve, ex->hsolver, mx->hmat, mx->hb, mx->hx);

  /* copy solution values back to PETSc */
  PetscCall(VecGetArray(y, &yy));
  PetscCallExternal(HYPRE_StructVectorGetBoxValues, mx->hx, hlower, hupper, (HYPRE_Complex *)yy);
  PetscCall(VecRestoreArray(y, &yy));
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscCall(PetscCitationsRegister(hypreCitation, &cite));
  PetscCallExternal(HYPRE_StructPFMGSetMaxIter, jac->hsolver, its * jac->its);
  PetscCallExternal(HYPRE_StructPFMGSetTol, jac->hsolver, rtol);

  PetscCall(PCApply_PFMG(pc, b, y));
  PetscCallExternal(HYPRE_StructPFMGGetNumIterations, jac->hsolver, &oits);
  *outits = oits;
  if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
  else *reason = PCRICHARDSON_CONVERGED_RTOL;
  PetscCallExternal(HYPRE_StructPFMGSetTol, jac->hsolver, jac->tol);
  PetscCallExternal(HYPRE_StructPFMGSetMaxIter, jac->hsolver, jac->its);
  PetscFunctionReturn(0);
}

PetscErrorCode PCSetUp_PFMG(PC pc) {
  PC_PFMG         *ex = (PC_PFMG *)pc->data;
  Mat_HYPREStruct *mx = (Mat_HYPREStruct *)(pc->pmat->data);
  PetscBool        flg;

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRESTRUCT, &flg));
  PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use MATHYPRESTRUCT with this preconditioner");

  /* create the hypre solver object and set its information */
  if (ex->hsolver) PetscCallExternal(HYPRE_StructPFMGDestroy, ex->hsolver);
  PetscCallExternal(HYPRE_StructPFMGCreate, ex->hcomm, &ex->hsolver);

  // Print Hypre statistics about the solve process
  if (ex->print_statistics) PetscCallExternal(HYPRE_StructPFMGSetPrintLevel, ex->hsolver, 3);

  // The hypre options must be repeated here because the StructPFMG was destroyed and recreated
  PetscCallExternal(HYPRE_StructPFMGSetMaxIter, ex->hsolver, ex->its);
  PetscCallExternal(HYPRE_StructPFMGSetNumPreRelax, ex->hsolver, ex->num_pre_relax);
  PetscCallExternal(HYPRE_StructPFMGSetNumPostRelax, ex->hsolver, ex->num_post_relax);
  PetscCallExternal(HYPRE_StructPFMGSetMaxLevels, ex->hsolver, ex->max_levels);
  PetscCallExternal(HYPRE_StructPFMGSetTol, ex->hsolver, ex->tol);
  PetscCallExternal(HYPRE_StructPFMGSetRelaxType, ex->hsolver, ex->relax_type);
  PetscCallExternal(HYPRE_StructPFMGSetRAPType, ex->hsolver, ex->rap_type);

  PetscCallExternal(HYPRE_StructPFMGSetup, ex->hsolver, mx->hmat, mx->hb, mx->hx);
  PetscCallExternal(HYPRE_StructPFMGSetZeroGuess, ex->hsolver);
  PetscFunctionReturn(0);
}

/*MC
     PCPFMG - the hypre PFMG multigrid solver

   Level: advanced

   Options Database:
+ -pc_pfmg_its <its> - number of iterations of PFMG to use as preconditioner
. -pc_pfmg_num_pre_relax <steps> - number of smoothing steps before coarse grid solve
. -pc_pfmg_num_post_relax <steps> - number of smoothing steps after coarse grid solve
. -pc_pfmg_tol <tol> - tolerance of PFMG
. -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
. -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin
- -pc_pfmg_skip_relax - skip relaxation on certain grids for isotropic problems. This can greatly improve efficiency by eliminating unnecessary relaxations when the underlying problem is isotropic, one of 0,1

   Notes:
    This is for CELL-centered descretizations

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

.seealso: `PCMG`, `MATHYPRESTRUCT`
M*/

PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc) {
  PC_PFMG *ex;

  PetscFunctionBegin;
  PetscCall(PetscNew(&ex));
  pc->data = ex;

  ex->its              = 1;
  ex->tol              = 1.e-8;
  ex->relax_type       = 1;
  ex->rap_type         = 0;
  ex->num_pre_relax    = 1;
  ex->num_post_relax   = 1;
  ex->max_levels       = 0;
  ex->skip_relax       = 0;
  ex->print_statistics = PETSC_FALSE;

  pc->ops->setfromoptions  = PCSetFromOptions_PFMG;
  pc->ops->view            = PCView_PFMG;
  pc->ops->destroy         = PCDestroy_PFMG;
  pc->ops->apply           = PCApply_PFMG;
  pc->ops->applyrichardson = PCApplyRichardson_PFMG;
  pc->ops->setup           = PCSetUp_PFMG;

  PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
  PetscCallExternal(HYPRE_StructPFMGCreate, ex->hcomm, &ex->hsolver);
  PetscFunctionReturn(0);
}

/* ---------------------------------------------------------------------------------------------------------------------------------------------------*/

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

  /* keep copy of SYSPFMG options used so may view them */
  PetscInt its;
  double   tol;
  PetscInt relax_type;
  PetscInt num_pre_relax, num_post_relax;
} PC_SysPFMG;

PetscErrorCode PCDestroy_SysPFMG(PC pc) {
  PC_SysPFMG *ex = (PC_SysPFMG *)pc->data;

  PetscFunctionBegin;
  if (ex->ss_solver) PetscCallExternal(HYPRE_SStructSysPFMGDestroy, ex->ss_solver);
  PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
  PetscCall(PetscFree(pc->data));
  PetscFunctionReturn(0);
}

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

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

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
  if (iascii) {
    PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE SysPFMG preconditioning\n"));
    PetscCall(PetscViewerASCIIPrintf(viewer, "  max iterations %" PetscInt_FMT "\n", ex->its));
    PetscCall(PetscViewerASCIIPrintf(viewer, "  tolerance %g\n", ex->tol));
    PetscCall(PetscViewerASCIIPrintf(viewer, "  relax type %s\n", PFMGRelaxType[ex->relax_type]));
    PetscCall(PetscViewerASCIIPrintf(viewer, "  number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax));
  }
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscOptionsHeadBegin(PetscOptionsObject, "SysPFMG options");
  PetscCall(PetscOptionsBool("-pc_syspfmg_print_statistics", "Print statistics", "HYPRE_SStructSysPFMGSetPrintLevel", flg, &flg, NULL));
  if (flg) PetscCallExternal(HYPRE_SStructSysPFMGSetPrintLevel, ex->ss_solver, 3);
  PetscCall(PetscOptionsInt("-pc_syspfmg_its", "Number of iterations of SysPFMG to use as preconditioner", "HYPRE_SStructSysPFMGSetMaxIter", ex->its, &ex->its, NULL));
  PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter, ex->ss_solver, ex->its);
  PetscCall(PetscOptionsInt("-pc_syspfmg_num_pre_relax", "Number of smoothing steps before coarse grid", "HYPRE_SStructSysPFMGSetNumPreRelax", ex->num_pre_relax, &ex->num_pre_relax, NULL));
  PetscCallExternal(HYPRE_SStructSysPFMGSetNumPreRelax, ex->ss_solver, ex->num_pre_relax);
  PetscCall(PetscOptionsInt("-pc_syspfmg_num_post_relax", "Number of smoothing steps after coarse grid", "HYPRE_SStructSysPFMGSetNumPostRelax", ex->num_post_relax, &ex->num_post_relax, NULL));
  PetscCallExternal(HYPRE_SStructSysPFMGSetNumPostRelax, ex->ss_solver, ex->num_post_relax);

  PetscCall(PetscOptionsReal("-pc_syspfmg_tol", "Tolerance of SysPFMG", "HYPRE_SStructSysPFMGSetTol", ex->tol, &ex->tol, NULL));
  PetscCallExternal(HYPRE_SStructSysPFMGSetTol, ex->ss_solver, ex->tol);
  PetscCall(PetscOptionsEList("-pc_syspfmg_relax_type", "Relax type for the up and down cycles", "HYPRE_SStructSysPFMGSetRelaxType", SysPFMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(SysPFMGRelaxType), SysPFMGRelaxType[ex->relax_type], &ex->relax_type, NULL));
  PetscCallExternal(HYPRE_SStructSysPFMGSetRelaxType, ex->ss_solver, ex->relax_type);
  PetscOptionsHeadEnd();
  PetscFunctionReturn(0);
}

PetscErrorCode PCApply_SysPFMG(PC pc, Vec x, Vec y) {
  PC_SysPFMG        *ex = (PC_SysPFMG *)pc->data;
  PetscScalar       *yy;
  const PetscScalar *xx;
  PetscInt           ilower[3], iupper[3];
  HYPRE_Int          hlower[3], hupper[3];
  Mat_HYPRESStruct  *mx       = (Mat_HYPRESStruct *)(pc->pmat->data);
  PetscInt           ordering = mx->dofs_order;
  PetscInt           nvars    = mx->nvars;
  PetscInt           part     = 0;
  PetscInt           size;
  PetscInt           i;

  PetscFunctionBegin;
  PetscCall(PetscCitationsRegister(hypreCitation, &cite));
  PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
  /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
  iupper[0] += ilower[0] - 1;
  iupper[1] += ilower[1] - 1;
  iupper[2] += ilower[2] - 1;
  hlower[0] = (HYPRE_Int)ilower[0];
  hlower[1] = (HYPRE_Int)ilower[1];
  hlower[2] = (HYPRE_Int)ilower[2];
  hupper[0] = (HYPRE_Int)iupper[0];
  hupper[1] = (HYPRE_Int)iupper[1];
  hupper[2] = (HYPRE_Int)iupper[2];

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

  /* copy x values over to hypre for variable ordering */
  if (ordering) {
    PetscCallExternal(HYPRE_SStructVectorSetConstantValues, mx->ss_b, 0.0);
    PetscCall(VecGetArrayRead(x, &xx));
    for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorSetBoxValues, mx->ss_b, part, hlower, hupper, i, (HYPRE_Complex *)(xx + (size * i)));
    PetscCall(VecRestoreArrayRead(x, &xx));
    PetscCallExternal(HYPRE_SStructVectorAssemble, mx->ss_b);
    PetscCallExternal(HYPRE_SStructMatrixMatvec, 1.0, mx->ss_mat, mx->ss_b, 0.0, mx->ss_x);
    PetscCallExternal(HYPRE_SStructSysPFMGSolve, ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x);

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

    PetscCall(PetscMalloc1(nvars * size, &z));
    PetscCallExternal(HYPRE_SStructVectorSetConstantValues, mx->ss_b, 0.0);
    PetscCall(VecGetArrayRead(x, &xx));

    /* transform nodal to hypre's variable ordering for sys_pfmg */
    for (i = 0; i < size; i++) {
      k = i * nvars;
      for (j = 0; j < nvars; j++) z[j * size + i] = xx[k + j];
    }
    for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorSetBoxValues, mx->ss_b, part, hlower, hupper, i, (HYPRE_Complex *)(z + (size * i)));
    PetscCall(VecRestoreArrayRead(x, &xx));
    PetscCallExternal(HYPRE_SStructVectorAssemble, mx->ss_b);
    PetscCallExternal(HYPRE_SStructSysPFMGSolve, ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x);

    /* copy solution values back to PETSc */
    PetscCall(VecGetArray(y, &yy));
    for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues, mx->ss_x, part, hlower, hupper, i, (HYPRE_Complex *)(z + (size * i)));
    /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
    for (i = 0; i < size; i++) {
      k = i * nvars;
      for (j = 0; j < nvars; j++) yy[k + j] = z[j * size + i];
    }
    PetscCall(VecRestoreArray(y, &yy));
    PetscCall(PetscFree(z));
  }
  PetscFunctionReturn(0);
}

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

  PetscFunctionBegin;
  PetscCall(PetscCitationsRegister(hypreCitation, &cite));
  PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter, jac->ss_solver, its * jac->its);
  PetscCallExternal(HYPRE_SStructSysPFMGSetTol, jac->ss_solver, rtol);
  PetscCall(PCApply_SysPFMG(pc, b, y));
  PetscCallExternal(HYPRE_SStructSysPFMGGetNumIterations, jac->ss_solver, &oits);
  *outits = oits;
  if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
  else *reason = PCRICHARDSON_CONVERGED_RTOL;
  PetscCallExternal(HYPRE_SStructSysPFMGSetTol, jac->ss_solver, jac->tol);
  PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter, jac->ss_solver, jac->its);
  PetscFunctionReturn(0);
}

PetscErrorCode PCSetUp_SysPFMG(PC pc) {
  PC_SysPFMG       *ex = (PC_SysPFMG *)pc->data;
  Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)(pc->pmat->data);
  PetscBool         flg;

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRESSTRUCT, &flg));
  PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use MATHYPRESSTRUCT with this preconditioner");

  /* create the hypre sstruct solver object and set its information */
  if (ex->ss_solver) PetscCallExternal(HYPRE_SStructSysPFMGDestroy, ex->ss_solver);
  PetscCallExternal(HYPRE_SStructSysPFMGCreate, ex->hcomm, &ex->ss_solver);
  PetscCallExternal(HYPRE_SStructSysPFMGSetZeroGuess, ex->ss_solver);
  PetscCallExternal(HYPRE_SStructSysPFMGSetup, ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x);
  PetscFunctionReturn(0);
}

/*MC
     PCSysPFMG - the hypre SysPFMG multigrid solver

   Level: advanced

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

   Notes:
    This is for CELL-centered descretizations

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

.seealso: `PCMG`, `MATHYPRESSTRUCT`
M*/

PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc) {
  PC_SysPFMG *ex;

  PetscFunctionBegin;
  PetscCall(PetscNew(&ex));
  pc->data = ex;

  ex->its            = 1;
  ex->tol            = 1.e-8;
  ex->relax_type     = 1;
  ex->num_pre_relax  = 1;
  ex->num_post_relax = 1;

  pc->ops->setfromoptions  = PCSetFromOptions_SysPFMG;
  pc->ops->view            = PCView_SysPFMG;
  pc->ops->destroy         = PCDestroy_SysPFMG;
  pc->ops->apply           = PCApply_SysPFMG;
  pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
  pc->ops->setup           = PCSetUp_SysPFMG;

  PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
  PetscCallExternal(HYPRE_SStructSysPFMGCreate, ex->hcomm, &ex->ss_solver);
  PetscFunctionReturn(0);
}

/* ---------------------------------------------------------------------------------------------------------------------------------*/

// PC SMG
typedef struct {
  MPI_Comm           hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */
  HYPRE_StructSolver hsolver;
  PetscInt           its; /* keep copy of SMG options used so may view them */
  double             tol;
  PetscBool          print_statistics;
  PetscInt           num_pre_relax, num_post_relax;
} PC_SMG;

PetscErrorCode PCDestroy_SMG(PC pc) {
  PC_SMG *ex = (PC_SMG *)pc->data;

  PetscFunctionBegin;
  if (ex->hsolver) PetscCallExternal(HYPRE_StructSMGDestroy, ex->hsolver);
  PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
  PetscCall(PetscFree(pc->data));
  PetscFunctionReturn(0);
}

PetscErrorCode PCView_SMG(PC pc, PetscViewer viewer) {
  PetscBool iascii;
  PC_SMG   *ex = (PC_SMG *)pc->data;

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
  if (iascii) {
    PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE SMG preconditioning\n"));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    max iterations %" PetscInt_FMT "\n", ex->its));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    tolerance %g\n", ex->tol));
    PetscCall(PetscViewerASCIIPrintf(viewer, "    number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax));
  }
  PetscFunctionReturn(0);
}

PetscErrorCode PCSetFromOptions_SMG(PC pc, PetscOptionItems *PetscOptionsObject) {
  PC_SMG *ex = (PC_SMG *)pc->data;

  PetscFunctionBegin;
  PetscOptionsHeadBegin(PetscOptionsObject, "SMG options");

  PetscCall(PetscOptionsInt("-pc_smg_its", "Number of iterations of SMG to use as preconditioner", "HYPRE_StructSMGSetMaxIter", ex->its, &ex->its, NULL));
  PetscCall(PetscOptionsInt("-pc_smg_num_pre_relax", "Number of smoothing steps before coarse grid", "HYPRE_StructSMGSetNumPreRelax", ex->num_pre_relax, &ex->num_pre_relax, NULL));
  PetscCall(PetscOptionsInt("-pc_smg_num_post_relax", "Number of smoothing steps after coarse grid", "HYPRE_StructSMGSetNumPostRelax", ex->num_post_relax, &ex->num_post_relax, NULL));
  PetscCall(PetscOptionsReal("-pc_smg_tol", "Tolerance of SMG", "HYPRE_StructSMGSetTol", ex->tol, &ex->tol, NULL));

  PetscOptionsHeadEnd();
  PetscFunctionReturn(0);
}

PetscErrorCode PCApply_SMG(PC pc, Vec x, Vec y) {
  PC_SMG            *ex = (PC_SMG *)pc->data;
  PetscScalar       *yy;
  const PetscScalar *xx;
  PetscInt           ilower[3], iupper[3];
  HYPRE_Int          hlower[3], hupper[3];
  Mat_HYPREStruct   *mx = (Mat_HYPREStruct *)(pc->pmat->data);

  PetscFunctionBegin;
  PetscCall(PetscCitationsRegister(hypreCitation, &cite));
  PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
  /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
  iupper[0] += ilower[0] - 1;
  iupper[1] += ilower[1] - 1;
  iupper[2] += ilower[2] - 1;
  hlower[0] = (HYPRE_Int)ilower[0];
  hlower[1] = (HYPRE_Int)ilower[1];
  hlower[2] = (HYPRE_Int)ilower[2];
  hupper[0] = (HYPRE_Int)iupper[0];
  hupper[1] = (HYPRE_Int)iupper[1];
  hupper[2] = (HYPRE_Int)iupper[2];

  /* copy x values over to hypre */
  PetscCallExternal(HYPRE_StructVectorSetConstantValues, mx->hb, 0.0);
  PetscCall(VecGetArrayRead(x, &xx));
  PetscCallExternal(HYPRE_StructVectorSetBoxValues, mx->hb, hlower, hupper, (HYPRE_Complex *)xx);
  PetscCall(VecRestoreArrayRead(x, &xx));
  PetscCallExternal(HYPRE_StructVectorAssemble, mx->hb);
  PetscCallExternal(HYPRE_StructSMGSolve, ex->hsolver, mx->hmat, mx->hb, mx->hx);

  /* copy solution values back to PETSc */
  PetscCall(VecGetArray(y, &yy));
  PetscCallExternal(HYPRE_StructVectorGetBoxValues, mx->hx, hlower, hupper, (HYPRE_Complex *)yy);
  PetscCall(VecRestoreArray(y, &yy));
  PetscFunctionReturn(0);
}

static PetscErrorCode PCApplyRichardson_SMG(PC pc, Vec b, Vec y, Vec w, PetscReal rtol, PetscReal abstol, PetscReal dtol, PetscInt its, PetscBool guesszero, PetscInt *outits, PCRichardsonConvergedReason *reason) {
  PC_SMG   *jac = (PC_SMG *)pc->data;
  HYPRE_Int oits;

  PetscFunctionBegin;
  PetscCall(PetscCitationsRegister(hypreCitation, &cite));
  PetscCallExternal(HYPRE_StructSMGSetMaxIter, jac->hsolver, its * jac->its);
  PetscCallExternal(HYPRE_StructSMGSetTol, jac->hsolver, rtol);

  PetscCall(PCApply_SMG(pc, b, y));
  PetscCallExternal(HYPRE_StructSMGGetNumIterations, jac->hsolver, &oits);
  *outits = oits;
  if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
  else *reason = PCRICHARDSON_CONVERGED_RTOL;
  PetscCallExternal(HYPRE_StructSMGSetTol, jac->hsolver, jac->tol);
  PetscCallExternal(HYPRE_StructSMGSetMaxIter, jac->hsolver, jac->its);
  PetscFunctionReturn(0);
}

PetscErrorCode PCSetUp_SMG(PC pc) {
  PetscInt         i, dim;
  PC_SMG          *ex = (PC_SMG *)pc->data;
  Mat_HYPREStruct *mx = (Mat_HYPREStruct *)(pc->pmat->data);
  PetscBool        flg;
  DMBoundaryType   p[3];
  PetscInt         M[3];

  PetscFunctionBegin;
  PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRESTRUCT, &flg));
  PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use MATHYPRESTRUCT with this preconditioner");

  PetscCall(DMDAGetInfo(mx->da, &dim, &M[0], &M[1], &M[2], 0, 0, 0, 0, 0, &p[0], &p[1], &p[2], 0));
  // Check if power of 2 in periodic directions
  for (i = 0; i < dim; i++) {
    if (((M[i] & (M[i] - 1)) != 0) && (p[i] == DM_BOUNDARY_PERIODIC)) {
      SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "With SMG, the number of points in a periodic direction must be a power of 2, but is here %" PetscInt_FMT ".", M[i]);
    }
  }

  /* create the hypre solver object and set its information */
  if (ex->hsolver) PetscCallExternal(HYPRE_StructSMGDestroy, (ex->hsolver));
  PetscCallExternal(HYPRE_StructSMGCreate, ex->hcomm, &ex->hsolver);
  // The hypre options must be set here and not in SetFromOptions because it is created here!
  PetscCallExternal(HYPRE_StructSMGSetMaxIter, ex->hsolver, ex->its);
  PetscCallExternal(HYPRE_StructSMGSetNumPreRelax, ex->hsolver, ex->num_pre_relax);
  PetscCallExternal(HYPRE_StructSMGSetNumPostRelax, ex->hsolver, ex->num_post_relax);
  PetscCallExternal(HYPRE_StructSMGSetTol, ex->hsolver, ex->tol);

  PetscCallExternal(HYPRE_StructSMGSetup, ex->hsolver, mx->hmat, mx->hb, mx->hx);
  PetscCallExternal(HYPRE_StructSMGSetZeroGuess, ex->hsolver);
  PetscFunctionReturn(0);
}

/*MC
     PCSMG - the hypre (structured grid) SMG multigrid solver

   Level: advanced

   Options Database:
+ -pc_smg_its <its> - number of iterations of SMG to use as preconditioner
. -pc_smg_num_pre_relax <steps> - number of smoothing steps before coarse grid
. -pc_smg_num_post_relax <steps> - number of smoothing steps after coarse grid
- -pc_smg_tol <tol> - tolerance of SMG

   Notes:
   This is for CELL-centered descretizations

  This must be used with the `MATHYPRESTRUCT` `MatType`.
  This support is less general than in hypre, it supports only one block per process defined by a PETSc `DMDA`.

.seealso:  PCMG, MATHYPRESTRUCT, PCPFMG
M*/

PETSC_EXTERN PetscErrorCode PCCreate_SMG(PC pc) {
  PC_SMG *ex;

  PetscFunctionBegin;
  PetscCall(PetscNew(&ex));
  pc->data = ex;

  ex->its            = 1;
  ex->tol            = 1.e-8;
  ex->num_pre_relax  = 1;
  ex->num_post_relax = 1;

  pc->ops->setfromoptions  = PCSetFromOptions_SMG;
  pc->ops->view            = PCView_SMG;
  pc->ops->destroy         = PCDestroy_SMG;
  pc->ops->apply           = PCApply_SMG;
  pc->ops->applyrichardson = PCApplyRichardson_SMG;
  pc->ops->setup           = PCSetUp_SMG;

  PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
  PetscCallExternal(HYPRE_StructSMGCreate, ex->hcomm, &ex->hsolver);
  PetscFunctionReturn(0);
}