Actual source code: hypre.c

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

  5: #include <petscpkg_version.h>
  6: #include <petsc/private/pcimpl.h>
  7: /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */
  8: #include <petsc/private/matimpl.h>
  9: #include <petsc/private/vecimpl.h>
 10: #include <../src/vec/vec/impls/hypre/vhyp.h>
 11: #include <../src/mat/impls/hypre/mhypre.h>
 12: #include <../src/dm/impls/da/hypre/mhyp.h>
 13: #include <_hypre_parcsr_ls.h>
 14: #include <petscmathypre.h>

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

 20: static PetscBool  cite            = PETSC_FALSE;
 21: static const char hypreCitation[] = "@manual{hypre-web-page,\n  title  = {{\\sl hypre}: High Performance Preconditioners},\n  organization = {Lawrence Livermore National Laboratory},\n  note  = "
 22:                                     "{\\url{https://www.llnl.gov/casc/hypre}}\n}\n";

 24: /*
 25:    Private context (data structure) for the  preconditioner.
 26: */
 27: typedef struct {
 28:   HYPRE_Solver hsolver;
 29:   Mat          hpmat; /* MatHYPRE */

 31:   HYPRE_Int (*destroy)(HYPRE_Solver);
 32:   HYPRE_Int (*solve)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector);
 33:   HYPRE_Int (*setup)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector);

 35:   MPI_Comm comm_hypre;
 36:   char    *hypre_type;

 38:   /* options for Pilut and BoomerAMG*/
 39:   PetscInt  maxiter;
 40:   PetscReal tol;

 42:   /* options for Pilut */
 43:   PetscInt factorrowsize;

 45:   /* options for ParaSails */
 46:   PetscInt  nlevels;
 47:   PetscReal threshold;
 48:   PetscReal filter;
 49:   PetscReal loadbal;
 50:   PetscInt  logging;
 51:   PetscInt  ruse;
 52:   PetscInt  symt;

 54:   /* options for BoomerAMG */
 55:   PetscBool printstatistics;

 57:   /* options for BoomerAMG */
 58:   PetscInt  cycletype;
 59:   PetscInt  maxlevels;
 60:   PetscReal strongthreshold;
 61:   PetscReal maxrowsum;
 62:   PetscInt  gridsweeps[3];
 63:   PetscInt  coarsentype;
 64:   PetscInt  measuretype;
 65:   PetscInt  smoothtype;
 66:   PetscInt  smoothnumlevels;
 67:   PetscInt  eu_level;         /* Number of levels for ILU(k) in Euclid */
 68:   PetscReal eu_droptolerance; /* Drop tolerance for ILU(k) in Euclid */
 69:   PetscInt  eu_bj;            /* Defines use of Block Jacobi ILU in Euclid */
 70:   PetscInt  relaxtype[3];
 71:   PetscReal relaxweight;
 72:   PetscReal outerrelaxweight;
 73:   PetscInt  relaxorder;
 74:   PetscReal truncfactor;
 75:   PetscBool applyrichardson;
 76:   PetscInt  pmax;
 77:   PetscInt  interptype;
 78:   PetscInt  maxc;
 79:   PetscInt  minc;
 80: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
 81:   char *spgemm_type; // this is a global hypre parameter but is closely associated with BoomerAMG
 82: #endif
 83:   /* GPU */
 84:   PetscBool keeptranspose;
 85:   PetscInt  rap2;
 86:   PetscInt  mod_rap2;

 88:   /* AIR */
 89:   PetscInt  Rtype;
 90:   PetscReal Rstrongthreshold;
 91:   PetscReal Rfilterthreshold;
 92:   PetscInt  Adroptype;
 93:   PetscReal Adroptol;

 95:   PetscInt  agg_nl;
 96:   PetscInt  agg_interptype;
 97:   PetscInt  agg_num_paths;
 98:   PetscBool nodal_relax;
 99:   PetscInt  nodal_relax_levels;

101:   PetscInt  nodal_coarsening;
102:   PetscInt  nodal_coarsening_diag;
103:   PetscInt  vec_interp_variant;
104:   PetscInt  vec_interp_qmax;
105:   PetscBool vec_interp_smooth;
106:   PetscInt  interp_refine;

108:   /* NearNullSpace support */
109:   VecHYPRE_IJVector *hmnull;
110:   HYPRE_ParVector   *phmnull;
111:   PetscInt           n_hmnull;
112:   Vec                hmnull_constant;

114:   /* options for AS (Auxiliary Space preconditioners) */
115:   PetscInt  as_print;
116:   PetscInt  as_max_iter;
117:   PetscReal as_tol;
118:   PetscInt  as_relax_type;
119:   PetscInt  as_relax_times;
120:   PetscReal as_relax_weight;
121:   PetscReal as_omega;
122:   PetscInt  as_amg_alpha_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for vector Poisson (AMS) or Curl problem (ADS) */
123:   PetscReal as_amg_alpha_theta;   /* AMG strength for vector Poisson (AMS) or Curl problem (ADS) */
124:   PetscInt  as_amg_beta_opts[5];  /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for scalar Poisson (AMS) or vector Poisson (ADS) */
125:   PetscReal as_amg_beta_theta;    /* AMG strength for scalar Poisson (AMS) or vector Poisson (ADS)  */
126:   PetscInt  ams_cycle_type;
127:   PetscInt  ads_cycle_type;

129:   /* additional data */
130:   Mat G;             /* MatHYPRE */
131:   Mat C;             /* MatHYPRE */
132:   Mat alpha_Poisson; /* MatHYPRE */
133:   Mat beta_Poisson;  /* MatHYPRE */

135:   /* extra information for AMS */
136:   PetscInt          dim; /* geometrical dimension */
137:   VecHYPRE_IJVector coords[3];
138:   VecHYPRE_IJVector constants[3];
139:   VecHYPRE_IJVector interior;
140:   Mat               RT_PiFull, RT_Pi[3];
141:   Mat               ND_PiFull, ND_Pi[3];
142:   PetscBool         ams_beta_is_zero;
143:   PetscBool         ams_beta_is_zero_part;
144:   PetscInt          ams_proj_freq;
145: } PC_HYPRE;

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

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

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

189:   PetscFunctionBegin;
190:   PetscCall(PetscStrcmp(jac->hypre_type, "boomeramg", &same));
191:   PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_NOTSAMETYPE, "Hypre type is not BoomerAMG ");
192:   num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData *)jac->hsolver);
193:   PetscCall(PetscMalloc1(num_levels, &mattmp));
194:   P_array = hypre_ParAMGDataPArray((hypre_ParAMGData *)jac->hsolver);
195:   for (l = 1; l < num_levels; l++) {
196:     PetscCall(MatCreateFromParCSR(P_array[num_levels - 1 - l], MATAIJ, PETSC_OWN_POINTER, &mattmp[l - 1]));
197:     /* We want to own the data, and HYPRE can not touch this matrix any more */
198:     P_array[num_levels - 1 - l] = NULL;
199:   }
200:   *nlevels        = num_levels;
201:   *interpolations = mattmp;
202:   PetscFunctionReturn(PETSC_SUCCESS);
203: }

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

211:   PetscFunctionBegin;
212:   for (i = 0; i < jac->n_hmnull; i++) PetscCall(VecHYPRE_IJVectorDestroy(&jac->hmnull[i]));
213:   PetscCall(PetscFree(jac->hmnull));
214:   PetscCall(PetscFree(jac->phmnull));
215:   PetscCall(VecDestroy(&jac->hmnull_constant));
216:   jac->n_hmnull = 0;
217:   PetscFunctionReturn(PETSC_SUCCESS);
218: }

220: static PetscErrorCode PCSetUp_HYPRE(PC pc)
221: {
222:   PC_HYPRE          *jac = (PC_HYPRE *)pc->data;
223:   Mat_HYPRE         *hjac;
224:   HYPRE_ParCSRMatrix hmat;
225:   HYPRE_ParVector    bv, xv;
226:   PetscBool          ishypre;

228:   PetscFunctionBegin;
229:   /* default type is boomerAMG */
230:   if (!jac->hypre_type) PetscCall(PCHYPRESetType(pc, "boomeramg"));

232:   /* get hypre matrix */
233:   if (pc->flag == DIFFERENT_NONZERO_PATTERN) PetscCall(MatDestroy(&jac->hpmat));
234:   PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRE, &ishypre));
235:   if (!ishypre) {
236:     /* Temporary fix since we do not support MAT_REUSE_MATRIX with HYPRE device */
237: #if defined(PETSC_HAVE_HYPRE_DEVICE)
238:     PetscBool iscuda, iship, iskokkos;

240:     PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &iscuda, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
241:     PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &iship, MATSEQAIJHIPSPARSE, MATMPIAIJHIPSPARSE, ""));
242:     PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &iskokkos, MATSEQAIJKOKKOS, MATMPIAIJKOKKOS, ""));
243:     if (iscuda || iship || iskokkos) PetscCall(MatDestroy(&jac->hpmat));
244: #endif
245:     PetscCall(MatConvert(pc->pmat, MATHYPRE, jac->hpmat ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, &jac->hpmat));
246:   } else {
247:     PetscCall(PetscObjectReference((PetscObject)pc->pmat));
248:     PetscCall(MatDestroy(&jac->hpmat));
249:     jac->hpmat = pc->pmat;
250:   }

252:   /* allow debug */
253:   PetscCall(MatViewFromOptions(jac->hpmat, NULL, "-pc_hypre_mat_view"));
254:   hjac = (Mat_HYPRE *)jac->hpmat->data;

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

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

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

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

435:   PetscFunctionBegin;
436:   PetscCall(PetscCitationsRegister(hypreCitation, &cite));
437:   if (!jac->applyrichardson) PetscCall(VecSet(x, 0.0));
438:   PetscCall(VecHYPRE_IJVectorPushVecRead(hjac->b, b));
439:   if (jac->applyrichardson) PetscCall(VecHYPRE_IJVectorPushVec(hjac->x, x));
440:   else PetscCall(VecHYPRE_IJVectorPushVecWrite(hjac->x, x));
441:   PetscCallExternal(HYPRE_IJMatrixGetObject, hjac->ij, (void **)&hmat);
442:   PetscCallExternal(HYPRE_IJVectorGetObject, hjac->b->ij, (void **)&jbv);
443:   PetscCallExternal(HYPRE_IJVectorGetObject, hjac->x->ij, (void **)&jxv);
444:   PetscStackCallExternalVoid(
445:     "Hypre solve", do {
446:       HYPRE_Int hierr = (*jac->solve)(jac->hsolver, hmat, jbv, jxv);
447:       if (hierr) {
448:         PetscCheck(hierr == HYPRE_ERROR_CONV, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in HYPRE solver, error code %d", (int)hierr);
449:         HYPRE_ClearAllErrors();
450:       }
451:     } while (0));

453:   if (jac->setup == HYPRE_AMSSetup && jac->ams_beta_is_zero_part) PetscCallExternal(HYPRE_AMSProjectOutGradients, jac->hsolver, jxv);
454:   PetscCall(VecHYPRE_IJVectorPopVec(hjac->x));
455:   PetscCall(VecHYPRE_IJVectorPopVec(hjac->b));
456:   PetscFunctionReturn(PETSC_SUCCESS);
457: }

459: static PetscErrorCode PCMatApply_HYPRE_BoomerAMG(PC pc, Mat B, Mat X)
460: {
461:   PC_HYPRE           *jac  = (PC_HYPRE *)pc->data;
462:   Mat_HYPRE          *hjac = (Mat_HYPRE *)jac->hpmat->data;
463:   hypre_ParCSRMatrix *par_matrix;
464:   HYPRE_ParVector     hb, hx;
465:   const PetscScalar  *b;
466:   PetscScalar        *x;
467:   PetscInt            m, N, lda;
468:   hypre_Vector       *x_local;
469:   PetscMemType        type;

471:   PetscFunctionBegin;
472:   PetscCall(PetscCitationsRegister(hypreCitation, &cite));
473:   PetscCallExternal(HYPRE_IJMatrixGetObject, hjac->ij, (void **)&par_matrix);
474:   PetscCall(MatGetLocalSize(B, &m, NULL));
475:   PetscCall(MatGetSize(B, NULL, &N));
476:   PetscCallExternal(HYPRE_ParMultiVectorCreate, hypre_ParCSRMatrixComm(par_matrix), hypre_ParCSRMatrixGlobalNumRows(par_matrix), hypre_ParCSRMatrixRowStarts(par_matrix), N, &hb);
477:   PetscCallExternal(HYPRE_ParMultiVectorCreate, hypre_ParCSRMatrixComm(par_matrix), hypre_ParCSRMatrixGlobalNumRows(par_matrix), hypre_ParCSRMatrixRowStarts(par_matrix), N, &hx);
478:   PetscCall(MatZeroEntries(X));
479:   PetscCall(MatDenseGetArrayReadAndMemType(B, &b, &type));
480:   PetscCall(MatDenseGetLDA(B, &lda));
481:   PetscCheck(lda == m, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot use a LDA different than the number of local rows: % " PetscInt_FMT " != % " PetscInt_FMT, lda, m);
482:   PetscCall(MatDenseGetLDA(X, &lda));
483:   PetscCheck(lda == m, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot use a LDA different than the number of local rows: % " PetscInt_FMT " != % " PetscInt_FMT, lda, m);
484:   x_local = hypre_ParVectorLocalVector(hb);
485:   PetscCallExternal(hypre_SeqVectorSetDataOwner, x_local, 0);
486:   hypre_VectorData(x_local) = (HYPRE_Complex *)b;
487:   PetscCall(MatDenseGetArrayWriteAndMemType(X, &x, NULL));
488:   x_local = hypre_ParVectorLocalVector(hx);
489:   PetscCallExternal(hypre_SeqVectorSetDataOwner, x_local, 0);
490:   hypre_VectorData(x_local) = (HYPRE_Complex *)x;
491:   PetscCallExternal(hypre_ParVectorInitialize_v2, hb, type == PETSC_MEMTYPE_HOST ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE);
492:   PetscCallExternal(hypre_ParVectorInitialize_v2, hx, type == PETSC_MEMTYPE_HOST ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE);
493:   PetscStackCallExternalVoid(
494:     "Hypre solve", do {
495:       HYPRE_Int hierr = (*jac->solve)(jac->hsolver, par_matrix, hb, hx);
496:       if (hierr) {
497:         PetscCheck(hierr == HYPRE_ERROR_CONV, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in HYPRE solver, error code %d", (int)hierr);
498:         HYPRE_ClearAllErrors();
499:       }
500:     } while (0));
501:   PetscCallExternal(HYPRE_ParVectorDestroy, hb);
502:   PetscCallExternal(HYPRE_ParVectorDestroy, hx);
503:   PetscCall(MatDenseRestoreArrayReadAndMemType(B, &b));
504:   PetscCall(MatDenseRestoreArrayWriteAndMemType(X, &x));
505:   PetscFunctionReturn(PETSC_SUCCESS);
506: }

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

512:   PetscFunctionBegin;
513:   PetscCall(MatDestroy(&jac->hpmat));
514:   PetscCall(MatDestroy(&jac->G));
515:   PetscCall(MatDestroy(&jac->C));
516:   PetscCall(MatDestroy(&jac->alpha_Poisson));
517:   PetscCall(MatDestroy(&jac->beta_Poisson));
518:   PetscCall(MatDestroy(&jac->RT_PiFull));
519:   PetscCall(MatDestroy(&jac->RT_Pi[0]));
520:   PetscCall(MatDestroy(&jac->RT_Pi[1]));
521:   PetscCall(MatDestroy(&jac->RT_Pi[2]));
522:   PetscCall(MatDestroy(&jac->ND_PiFull));
523:   PetscCall(MatDestroy(&jac->ND_Pi[0]));
524:   PetscCall(MatDestroy(&jac->ND_Pi[1]));
525:   PetscCall(MatDestroy(&jac->ND_Pi[2]));
526:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[0]));
527:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[1]));
528:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[2]));
529:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[0]));
530:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[1]));
531:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[2]));
532:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->interior));
533:   PetscCall(PCHYPREResetNearNullSpace_Private(pc));
534:   jac->ams_beta_is_zero      = PETSC_FALSE;
535:   jac->ams_beta_is_zero_part = PETSC_FALSE;
536:   jac->dim                   = 0;
537:   PetscFunctionReturn(PETSC_SUCCESS);
538: }

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

544:   PetscFunctionBegin;
545:   PetscCall(PCReset_HYPRE(pc));
546:   if (jac->destroy) PetscCallExternal(jac->destroy, jac->hsolver);
547:   PetscCall(PetscFree(jac->hypre_type));
548: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
549:   PetscCall(PetscFree(jac->spgemm_type));
550: #endif
551:   if (jac->comm_hypre != MPI_COMM_NULL) PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
552:   PetscCall(PetscFree(pc->data));

554:   PetscCall(PetscObjectChangeTypeName((PetscObject)pc, 0));
555:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetType_C", NULL));
556:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetType_C", NULL));
557:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteGradient_C", NULL));
558:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteCurl_C", NULL));
559:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetInterpolations_C", NULL));
560:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetConstantEdgeVectors_C", NULL));
561:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetPoissonMatrix_C", NULL));
562:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetEdgeConstantVectors_C", NULL));
563:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREAMSSetInteriorNodes_C", NULL));
564:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", NULL));
565:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", NULL));
566:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinSetMatProductAlgorithm_C", NULL));
567:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinGetMatProductAlgorithm_C", NULL));
568:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
569:   PetscFunctionReturn(PETSC_SUCCESS);
570: }

572: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PC pc, PetscOptionItems *PetscOptionsObject)
573: {
574:   PC_HYPRE *jac = (PC_HYPRE *)pc->data;
575:   PetscBool flag;

577:   PetscFunctionBegin;
578:   PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE Pilut Options");
579:   PetscCall(PetscOptionsInt("-pc_hypre_pilut_maxiter", "Number of iterations", "None", jac->maxiter, &jac->maxiter, &flag));
580:   if (flag) PetscCallExternal(HYPRE_ParCSRPilutSetMaxIter, jac->hsolver, jac->maxiter);
581:   PetscCall(PetscOptionsReal("-pc_hypre_pilut_tol", "Drop tolerance", "None", jac->tol, &jac->tol, &flag));
582:   if (flag) PetscCallExternal(HYPRE_ParCSRPilutSetDropTolerance, jac->hsolver, jac->tol);
583:   PetscCall(PetscOptionsInt("-pc_hypre_pilut_factorrowsize", "FactorRowSize", "None", jac->factorrowsize, &jac->factorrowsize, &flag));
584:   if (flag) PetscCallExternal(HYPRE_ParCSRPilutSetFactorRowSize, jac->hsolver, jac->factorrowsize);
585:   PetscOptionsHeadEnd();
586:   PetscFunctionReturn(PETSC_SUCCESS);
587: }

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

594:   PetscFunctionBegin;
595:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
596:   if (iascii) {
597:     PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE Pilut preconditioning\n"));
598:     if (jac->maxiter != PETSC_DEFAULT) {
599:       PetscCall(PetscViewerASCIIPrintf(viewer, "    maximum number of iterations %" PetscInt_FMT "\n", jac->maxiter));
600:     } else {
601:       PetscCall(PetscViewerASCIIPrintf(viewer, "    default maximum number of iterations \n"));
602:     }
603:     if (jac->tol != PETSC_DEFAULT) {
604:       PetscCall(PetscViewerASCIIPrintf(viewer, "    drop tolerance %g\n", (double)jac->tol));
605:     } else {
606:       PetscCall(PetscViewerASCIIPrintf(viewer, "    default drop tolerance \n"));
607:     }
608:     if (jac->factorrowsize != PETSC_DEFAULT) {
609:       PetscCall(PetscViewerASCIIPrintf(viewer, "    factor row size %" PetscInt_FMT "\n", jac->factorrowsize));
610:     } else {
611:       PetscCall(PetscViewerASCIIPrintf(viewer, "    default factor row size \n"));
612:     }
613:   }
614:   PetscFunctionReturn(PETSC_SUCCESS);
615: }

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

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

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

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

636:   PetscCall(PetscOptionsBool("-pc_hypre_euclid_bj", "Use Block Jacobi for ILU in Euclid", "None", eu_bj, &eu_bj, &flag));
637:   if (flag) {
638:     jac->eu_bj = eu_bj ? 1 : 0;
639:     PetscCallExternal(HYPRE_EuclidSetBJ, jac->hsolver, jac->eu_bj);
640:   }
641:   PetscOptionsHeadEnd();
642:   PetscFunctionReturn(PETSC_SUCCESS);
643: }

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

650:   PetscFunctionBegin;
651:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
652:   if (iascii) {
653:     PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE Euclid preconditioning\n"));
654:     if (jac->eu_level != PETSC_DEFAULT) {
655:       PetscCall(PetscViewerASCIIPrintf(viewer, "    factorization levels %" PetscInt_FMT "\n", jac->eu_level));
656:     } else {
657:       PetscCall(PetscViewerASCIIPrintf(viewer, "    default factorization levels \n"));
658:     }
659:     PetscCall(PetscViewerASCIIPrintf(viewer, "    drop tolerance %g\n", (double)jac->eu_droptolerance));
660:     PetscCall(PetscViewerASCIIPrintf(viewer, "    use Block-Jacobi? %" PetscInt_FMT "\n", jac->eu_bj));
661:   }
662:   PetscFunctionReturn(PETSC_SUCCESS);
663: }

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

672:   PetscFunctionBegin;
673:   PetscCall(PetscCitationsRegister(hypreCitation, &cite));
674:   PetscCall(VecSet(x, 0.0));
675:   PetscCall(VecHYPRE_IJVectorPushVecRead(hjac->x, b));
676:   PetscCall(VecHYPRE_IJVectorPushVecWrite(hjac->b, x));

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

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

692:   PetscCall(VecHYPRE_IJVectorPopVec(hjac->x));
693:   PetscCall(VecHYPRE_IJVectorPopVec(hjac->b));
694:   PetscFunctionReturn(PETSC_SUCCESS);
695: }

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

702: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
703:   PetscFunctionBegin;
704:   if (jac->spgemm_type) {
705:     PetscCall(PetscStrcmp(jac->spgemm_type, name, &flag));
706:     PetscCheck(flag, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Cannot reset the HYPRE SpGEMM (really we can)");
707:     PetscFunctionReturn(PETSC_SUCCESS);
708:   } else {
709:     PetscCall(PetscStrallocpy(name, &jac->spgemm_type));
710:   }
711:   PetscCall(PetscStrcmp("cusparse", jac->spgemm_type, &flag));
712:   if (flag) {
713:     PetscCallExternal(HYPRE_SetSpGemmUseCusparse, 1);
714:     PetscFunctionReturn(PETSC_SUCCESS);
715:   }
716:   PetscCall(PetscStrcmp("hypre", jac->spgemm_type, &flag));
717:   if (flag) {
718:     PetscCallExternal(HYPRE_SetSpGemmUseCusparse, 0);
719:     PetscFunctionReturn(PETSC_SUCCESS);
720:   }
721:   jac->spgemm_type = NULL;
722:   SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown HYPRE SpGEMM type %s; Choices are cusparse, hypre", name);
723: #endif
724: }

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

730:   PetscFunctionBegin;
732: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
733:   *spgemm = jac->spgemm_type;
734: #endif
735:   PetscFunctionReturn(PETSC_SUCCESS);
736: }

738: static const char *HYPREBoomerAMGCycleType[]   = {"", "V", "W"};
739: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP", "Ruge-Stueben", "", "modifiedRuge-Stueben", "", "", "Falgout", "", "PMIS", "", "HMIS"};
740: static const char *HYPREBoomerAMGMeasureType[] = {"local", "global"};
741: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
742: static const char *HYPREBoomerAMGSmoothType[] = {"Schwarz-smoothers", "Pilut", "ParaSails", "Euclid"};
743: 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"};
744: 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"};
745: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PC pc, PetscOptionItems *PetscOptionsObject)
746: {
747:   PC_HYPRE   *jac = (PC_HYPRE *)pc->data;
748:   PetscInt    bs, n, indx, level;
749:   PetscBool   flg, tmp_truth;
750:   PetscReal   tmpdbl, twodbl[2];
751:   const char *symtlist[]           = {"nonsymmetric", "SPD", "nonsymmetric,SPD"};
752:   const char *PCHYPRESpgemmTypes[] = {"cusparse", "hypre"};

754:   PetscFunctionBegin;
755:   PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE BoomerAMG Options");
756:   PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_cycle_type", "Cycle type", "None", HYPREBoomerAMGCycleType + 1, 2, HYPREBoomerAMGCycleType[jac->cycletype], &indx, &flg));
757:   if (flg) {
758:     jac->cycletype = indx + 1;
759:     PetscCallExternal(HYPRE_BoomerAMGSetCycleType, jac->hsolver, jac->cycletype);
760:   }
761:   PetscCall(PetscOptionsBoundedInt("-pc_hypre_boomeramg_max_levels", "Number of levels (of grids) allowed", "None", jac->maxlevels, &jac->maxlevels, &flg, 2));
762:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetMaxLevels, jac->hsolver, jac->maxlevels);
763:   PetscCall(PetscOptionsBoundedInt("-pc_hypre_boomeramg_max_iter", "Maximum iterations used PER hypre call", "None", jac->maxiter, &jac->maxiter, &flg, 1));
764:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, jac->maxiter);
765:   PetscCall(PetscOptionsBoundedReal("-pc_hypre_boomeramg_tol", "Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)", "None", jac->tol, &jac->tol, &flg, 0.0));
766:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, jac->tol);
767:   bs = 1;
768:   if (pc->pmat) PetscCall(MatGetBlockSize(pc->pmat, &bs));
769:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_numfunctions", "Number of functions", "HYPRE_BoomerAMGSetNumFunctions", bs, &bs, &flg));
770:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetNumFunctions, jac->hsolver, bs);

772:   PetscCall(PetscOptionsBoundedReal("-pc_hypre_boomeramg_truncfactor", "Truncation factor for interpolation (0=no truncation)", "None", jac->truncfactor, &jac->truncfactor, &flg, 0.0));
773:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetTruncFactor, jac->hsolver, jac->truncfactor);

775:   PetscCall(PetscOptionsBoundedInt("-pc_hypre_boomeramg_P_max", "Max elements per row for interpolation operator (0=unlimited)", "None", jac->pmax, &jac->pmax, &flg, 0));
776:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetPMaxElmts, jac->hsolver, jac->pmax);

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

781:   PetscCall(PetscOptionsBoundedInt("-pc_hypre_boomeramg_agg_num_paths", "Number of paths for aggressive coarsening", "None", jac->agg_num_paths, &jac->agg_num_paths, &flg, 1));
782:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetNumPaths, jac->hsolver, jac->agg_num_paths);

784:   PetscCall(PetscOptionsBoundedReal("-pc_hypre_boomeramg_strong_threshold", "Threshold for being strongly connected", "None", jac->strongthreshold, &jac->strongthreshold, &flg, 0.0));
785:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetStrongThreshold, jac->hsolver, jac->strongthreshold);
786:   PetscCall(PetscOptionsRangeReal("-pc_hypre_boomeramg_max_row_sum", "Maximum row sum", "None", jac->maxrowsum, &jac->maxrowsum, &flg, 0.0, 1.0));
787:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetMaxRowSum, jac->hsolver, jac->maxrowsum);

789:   /* Grid sweeps */
790:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all", "Number of sweeps for the up and down grid levels", "None", jac->gridsweeps[0], &indx, &flg));
791:   if (flg) {
792:     PetscCallExternal(HYPRE_BoomerAMGSetNumSweeps, jac->hsolver, indx);
793:     /* modify the jac structure so we can view the updated options with PC_View */
794:     jac->gridsweeps[0] = indx;
795:     jac->gridsweeps[1] = indx;
796:     /*defaults coarse to 1 */
797:     jac->gridsweeps[2] = 1;
798:   }
799:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen", "Use a nodal based coarsening 1-6", "HYPRE_BoomerAMGSetNodal", jac->nodal_coarsening, &jac->nodal_coarsening, &flg));
800:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetNodal, jac->hsolver, jac->nodal_coarsening);
801:   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));
802:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetNodalDiag, jac->hsolver, jac->nodal_coarsening_diag);
803:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_variant", "Variant of algorithm 1-3", "HYPRE_BoomerAMGSetInterpVecVariant", jac->vec_interp_variant, &jac->vec_interp_variant, &flg));
804:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetInterpVecVariant, jac->hsolver, jac->vec_interp_variant);
805:   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));
806:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetInterpVecQMax, jac->hsolver, jac->vec_interp_qmax);
807:   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));
808:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetSmoothInterpVectors, jac->hsolver, jac->vec_interp_smooth);
809:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_interp_refine", "Preprocess the interpolation matrix through iterative weight refinement", "HYPRE_BoomerAMGSetInterpRefine", jac->interp_refine, &jac->interp_refine, &flg));
810:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetInterpRefine, jac->hsolver, jac->interp_refine);
811:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down", "Number of sweeps for the down cycles", "None", jac->gridsweeps[0], &indx, &flg));
812:   if (flg) {
813:     PetscCallExternal(HYPRE_BoomerAMGSetCycleNumSweeps, jac->hsolver, indx, 1);
814:     jac->gridsweeps[0] = indx;
815:   }
816:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up", "Number of sweeps for the up cycles", "None", jac->gridsweeps[1], &indx, &flg));
817:   if (flg) {
818:     PetscCallExternal(HYPRE_BoomerAMGSetCycleNumSweeps, jac->hsolver, indx, 2);
819:     jac->gridsweeps[1] = indx;
820:   }
821:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse", "Number of sweeps for the coarse level", "None", jac->gridsweeps[2], &indx, &flg));
822:   if (flg) {
823:     PetscCallExternal(HYPRE_BoomerAMGSetCycleNumSweeps, jac->hsolver, indx, 3);
824:     jac->gridsweeps[2] = indx;
825:   }

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

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

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

850:   /* Filter for ILU(k) for Euclid */
851:   PetscReal droptolerance;
852:   PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_eu_droptolerance", "Drop tolerance for ILU(k) in Euclid smoother", "None", 0, &droptolerance, &flg));
853:   if (flg && (jac->smoothtype == 3)) {
854:     jac->eu_droptolerance = droptolerance;
855:     PetscCallExternal(HYPRE_BoomerAMGSetEuLevel, jac->hsolver, droptolerance);
856:   }

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

865:   /* Relax type */
866:   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));
867:   if (flg) {
868:     jac->relaxtype[0] = jac->relaxtype[1] = indx;
869:     PetscCallExternal(HYPRE_BoomerAMGSetRelaxType, jac->hsolver, indx);
870:     /* by default, coarse type set to 9 */
871:     jac->relaxtype[2] = 9;
872:     PetscCallExternal(HYPRE_BoomerAMGSetCycleRelaxType, jac->hsolver, 9, 3);
873:   }
874:   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));
875:   if (flg) {
876:     jac->relaxtype[0] = indx;
877:     PetscCallExternal(HYPRE_BoomerAMGSetCycleRelaxType, jac->hsolver, indx, 1);
878:   }
879:   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));
880:   if (flg) {
881:     jac->relaxtype[1] = indx;
882:     PetscCallExternal(HYPRE_BoomerAMGSetCycleRelaxType, jac->hsolver, indx, 2);
883:   }
884:   PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse", "Relax type on coarse grid", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), HYPREBoomerAMGRelaxType[9], &indx, &flg));
885:   if (flg) {
886:     jac->relaxtype[2] = indx;
887:     PetscCallExternal(HYPRE_BoomerAMGSetCycleRelaxType, jac->hsolver, indx, 3);
888:   }

890:   /* Relaxation Weight */
891:   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));
892:   if (flg) {
893:     PetscCallExternal(HYPRE_BoomerAMGSetRelaxWt, jac->hsolver, tmpdbl);
894:     jac->relaxweight = tmpdbl;
895:   }

897:   n         = 2;
898:   twodbl[0] = twodbl[1] = 1.0;
899:   PetscCall(PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level", "Set the relaxation weight for a particular level (weight,level)", "None", twodbl, &n, &flg));
900:   if (flg) {
901:     PetscCheck(n == 2, 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);
902:     indx = (int)PetscAbsReal(twodbl[1]);
903:     PetscCallExternal(HYPRE_BoomerAMGSetLevelRelaxWt, jac->hsolver, twodbl[0], indx);
904:   }

906:   /* Outer relaxation Weight */
907:   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));
908:   if (flg) {
909:     PetscCallExternal(HYPRE_BoomerAMGSetOuterWt, jac->hsolver, tmpdbl);
910:     jac->outerrelaxweight = tmpdbl;
911:   }

913:   n         = 2;
914:   twodbl[0] = twodbl[1] = 1.0;
915:   PetscCall(PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level", "Set the outer relaxation weight for a particular level (weight,level)", "None", twodbl, &n, &flg));
916:   if (flg) {
917:     PetscCheck(n == 2, 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);
918:     indx = (int)PetscAbsReal(twodbl[1]);
919:     PetscCallExternal(HYPRE_BoomerAMGSetLevelOuterWt, jac->hsolver, twodbl[0], indx);
920:   }

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

925:   if (flg && tmp_truth) {
926:     jac->relaxorder = 0;
927:     PetscCallExternal(HYPRE_BoomerAMGSetRelaxOrder, jac->hsolver, jac->relaxorder);
928:   }
929:   PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_measure_type", "Measure type", "None", HYPREBoomerAMGMeasureType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGMeasureType), HYPREBoomerAMGMeasureType[0], &indx, &flg));
930:   if (flg) {
931:     jac->measuretype = indx;
932:     PetscCallExternal(HYPRE_BoomerAMGSetMeasureType, jac->hsolver, jac->measuretype);
933:   }
934:   /* update list length 3/07 */
935:   PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type", "Coarsen type", "None", HYPREBoomerAMGCoarsenType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGCoarsenType), HYPREBoomerAMGCoarsenType[6], &indx, &flg));
936:   if (flg) {
937:     jac->coarsentype = indx;
938:     PetscCallExternal(HYPRE_BoomerAMGSetCoarsenType, jac->hsolver, jac->coarsentype);
939:   }

941:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_max_coarse_size", "Maximum size of coarsest grid", "None", jac->maxc, &jac->maxc, &flg));
942:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetMaxCoarseSize, jac->hsolver, jac->maxc);
943:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_min_coarse_size", "Minimum size of coarsest grid", "None", jac->minc, &jac->minc, &flg));
944:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetMinCoarseSize, jac->hsolver, jac->minc);
945: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
946:   // global parameter but is closely associated with BoomerAMG
947:   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));
948:   #if defined(PETSC_HAVE_HYPRE_DEVICE)
949:   if (!flg) indx = 0;
950:   PetscCall(PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG(pc, PCHYPRESpgemmTypes[indx]));
951:   #else
952:   PetscCall(PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG(pc, "hypre"));
953:   #endif
954: #endif
955:   /* AIR */
956: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
957:   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));
958:   PetscCallExternal(HYPRE_BoomerAMGSetRestriction, jac->hsolver, jac->Rtype);
959:   if (jac->Rtype) {
960:     HYPRE_Int **grid_relax_points = hypre_TAlloc(HYPRE_Int *, 4, HYPRE_MEMORY_HOST);
961:     char       *prerelax[256];
962:     char       *postrelax[256];
963:     char        stringF[2] = "F", stringC[2] = "C", stringA[2] = "A";
964:     PetscInt    ns_down = 256, ns_up = 256;
965:     PetscBool   matchF, matchC, matchA;

967:     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 */

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

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

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

978:     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));
979:     PetscCallExternal(HYPRE_BoomerAMGSetADropType, jac->hsolver, jac->Adroptype);
980:     PetscCall(PetscOptionsStringArray("-pc_hypre_boomeramg_prerelax", "Defines prerelax scheme", "None", prerelax, &ns_down, NULL));
981:     PetscCall(PetscOptionsStringArray("-pc_hypre_boomeramg_postrelax", "Defines postrelax scheme", "None", postrelax, &ns_up, NULL));
982:     PetscCheck(ns_down == jac->gridsweeps[0], PetscObjectComm((PetscObject)jac), PETSC_ERR_ARG_SIZ, "The number of arguments passed to -pc_hypre_boomeramg_prerelax must match the number passed to -pc_hypre_bomeramg_grid_sweeps_down");
983:     PetscCheck(ns_up == jac->gridsweeps[1], PetscObjectComm((PetscObject)jac), PETSC_ERR_ARG_SIZ, "The number of arguments passed to -pc_hypre_boomeramg_postrelax must match the number passed to -pc_hypre_bomeramg_grid_sweeps_up");

985:     grid_relax_points[0]    = NULL;
986:     grid_relax_points[1]    = hypre_TAlloc(HYPRE_Int, ns_down, HYPRE_MEMORY_HOST);
987:     grid_relax_points[2]    = hypre_TAlloc(HYPRE_Int, ns_up, HYPRE_MEMORY_HOST);
988:     grid_relax_points[3]    = hypre_TAlloc(HYPRE_Int, jac->gridsweeps[2], HYPRE_MEMORY_HOST);
989:     grid_relax_points[3][0] = 0;

991:     // set down relax scheme
992:     for (PetscInt i = 0; i < ns_down; i++) {
993:       PetscCall(PetscStrcasecmp(prerelax[i], stringF, &matchF));
994:       PetscCall(PetscStrcasecmp(prerelax[i], stringC, &matchC));
995:       PetscCall(PetscStrcasecmp(prerelax[i], stringA, &matchA));
996:       PetscCheck(matchF || matchC || matchA, PetscObjectComm((PetscObject)jac), PETSC_ERR_ARG_WRONG, "Valid argument options for -pc_hypre_boomeramg_prerelax are C, F, and A");
997:       if (matchF) grid_relax_points[1][i] = -1;
998:       else if (matchC) grid_relax_points[1][i] = 1;
999:       else if (matchA) grid_relax_points[1][i] = 0;
1000:     }

1002:     // set up relax scheme
1003:     for (PetscInt i = 0; i < ns_up; i++) {
1004:       PetscCall(PetscStrcasecmp(postrelax[i], stringF, &matchF));
1005:       PetscCall(PetscStrcasecmp(postrelax[i], stringC, &matchC));
1006:       PetscCall(PetscStrcasecmp(postrelax[i], stringA, &matchA));
1007:       PetscCheck(matchF || matchC || matchA, PetscObjectComm((PetscObject)jac), PETSC_ERR_ARG_WRONG, "Valid argument options for -pc_hypre_boomeramg_postrelax are C, F, and A");
1008:       if (matchF) grid_relax_points[2][i] = -1;
1009:       else if (matchC) grid_relax_points[2][i] = 1;
1010:       else if (matchA) grid_relax_points[2][i] = 0;
1011:     }

1013:     // set coarse relax scheme
1014:     for (PetscInt i = 0; i < jac->gridsweeps[2]; i++) grid_relax_points[3][i] = 0;

1016:     // Pass relax schemes to hypre
1017:     PetscCallExternal(HYPRE_BoomerAMGSetGridRelaxPoints, jac->hsolver, grid_relax_points);

1019:     // cleanup memory
1020:     for (PetscInt i = 0; i < ns_down; i++) PetscCall(PetscFree(prerelax[i]));
1021:     for (PetscInt i = 0; i < ns_up; i++) PetscCall(PetscFree(postrelax[i]));
1022:   }
1023: #endif

1025: #if PETSC_PKG_HYPRE_VERSION_LE(9, 9, 9)
1026:   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);
1027: #endif

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

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

1041:     jac->printstatistics = PETSC_TRUE;
1042:     PetscCallExternal(HYPRE_BoomerAMGSetPrintLevel, jac->hsolver, level);
1043:   }

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

1050:     jac->printstatistics = PETSC_TRUE;
1051:     PetscCallExternal(HYPRE_BoomerAMGSetDebugFlag, jac->hsolver, level);
1052:   }

1054:   PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg));
1055:   if (flg && tmp_truth) {
1056:     PetscInt tmp_int;
1057:     PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", jac->nodal_relax_levels, &tmp_int, &flg));
1058:     if (flg) jac->nodal_relax_levels = tmp_int;
1059:     PetscCallExternal(HYPRE_BoomerAMGSetSmoothType, jac->hsolver, 6);
1060:     PetscCallExternal(HYPRE_BoomerAMGSetDomainType, jac->hsolver, 1);
1061:     PetscCallExternal(HYPRE_BoomerAMGSetOverlap, jac->hsolver, 0);
1062:     PetscCallExternal(HYPRE_BoomerAMGSetSmoothNumLevels, jac->hsolver, jac->nodal_relax_levels);
1063:   }

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

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

1075:   PetscOptionsHeadEnd();
1076:   PetscFunctionReturn(PETSC_SUCCESS);
1077: }

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

1084:   PetscFunctionBegin;
1085:   PetscCall(PetscCitationsRegister(hypreCitation, &cite));
1086:   PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, its * jac->maxiter);
1087:   PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, rtol);
1088:   jac->applyrichardson = PETSC_TRUE;
1089:   PetscCall(PCApply_HYPRE(pc, b, y));
1090:   jac->applyrichardson = PETSC_FALSE;
1091:   PetscCallExternal(HYPRE_BoomerAMGGetNumIterations, jac->hsolver, &oits);
1092:   *outits = oits;
1093:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1094:   else *reason = PCRICHARDSON_CONVERGED_RTOL;
1095:   PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, jac->tol);
1096:   PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, jac->maxiter);
1097:   PetscFunctionReturn(PETSC_SUCCESS);
1098: }

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

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

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

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

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

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

1133:     PetscCall(PetscViewerASCIIPrintf(viewer, "    Maximum size of coarsest grid %" PetscInt_FMT "\n", jac->maxc));
1134:     PetscCall(PetscViewerASCIIPrintf(viewer, "    Minimum size of coarsest grid %" PetscInt_FMT "\n", jac->minc));

1136:     if (jac->relaxorder) {
1137:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Using CF-relaxation\n"));
1138:     } else {
1139:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Not using CF-relaxation\n"));
1140:     }
1141:     if (jac->smoothtype != -1) {
1142:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Smooth type          %s\n", HYPREBoomerAMGSmoothType[jac->smoothtype]));
1143:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Smooth num levels    %" PetscInt_FMT "\n", jac->smoothnumlevels));
1144:     } else {
1145:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Not using more complex smoothers.\n"));
1146:     }
1147:     if (jac->smoothtype == 3) {
1148:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Euclid ILU(k) levels %" PetscInt_FMT "\n", jac->eu_level));
1149:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Euclid ILU(k) drop tolerance %g\n", (double)jac->eu_droptolerance));
1150:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Euclid ILU use Block-Jacobi? %" PetscInt_FMT "\n", jac->eu_bj));
1151:     }
1152:     PetscCall(PetscViewerASCIIPrintf(viewer, "    Measure type        %s\n", HYPREBoomerAMGMeasureType[jac->measuretype]));
1153:     PetscCall(PetscViewerASCIIPrintf(viewer, "    Coarsen type        %s\n", HYPREBoomerAMGCoarsenType[jac->coarsentype]));
1154:     PetscCall(PetscViewerASCIIPrintf(viewer, "    Interpolation type  %s\n", jac->interptype != 100 ? HYPREBoomerAMGInterpType[jac->interptype] : "1pt"));
1155:     if (jac->nodal_coarsening) PetscCall(PetscViewerASCIIPrintf(viewer, "    Using nodal coarsening with HYPRE_BOOMERAMGSetNodal() %" PetscInt_FMT "\n", jac->nodal_coarsening));
1156:     if (jac->vec_interp_variant) {
1157:       PetscCall(PetscViewerASCIIPrintf(viewer, "    HYPRE_BoomerAMGSetInterpVecVariant() %" PetscInt_FMT "\n", jac->vec_interp_variant));
1158:       PetscCall(PetscViewerASCIIPrintf(viewer, "    HYPRE_BoomerAMGSetInterpVecQMax() %" PetscInt_FMT "\n", jac->vec_interp_qmax));
1159:       PetscCall(PetscViewerASCIIPrintf(viewer, "    HYPRE_BoomerAMGSetSmoothInterpVectors() %d\n", jac->vec_interp_smooth));
1160:     }
1161:     if (jac->nodal_relax) PetscCall(PetscViewerASCIIPrintf(viewer, "    Using nodal relaxation via Schwarz smoothing on levels %" PetscInt_FMT "\n", jac->nodal_relax_levels));
1162: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
1163:     PetscCall(PetscViewerASCIIPrintf(viewer, "    SpGEMM type         %s\n", jac->spgemm_type));
1164: #else
1165:     PetscCall(PetscViewerASCIIPrintf(viewer, "    SpGEMM type         %s\n", "hypre"));
1166: #endif
1167:     /* AIR */
1168:     if (jac->Rtype) {
1169:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Using approximate ideal restriction type %" PetscInt_FMT "\n", jac->Rtype));
1170:       PetscCall(PetscViewerASCIIPrintf(viewer, "      Threshold for R %g\n", (double)jac->Rstrongthreshold));
1171:       PetscCall(PetscViewerASCIIPrintf(viewer, "      Filter for R %g\n", (double)jac->Rfilterthreshold));
1172:       PetscCall(PetscViewerASCIIPrintf(viewer, "      A drop tolerance %g\n", (double)jac->Adroptol));
1173:       PetscCall(PetscViewerASCIIPrintf(viewer, "      A drop type %" PetscInt_FMT "\n", jac->Adroptype));
1174:     }
1175:   }
1176:   PetscFunctionReturn(PETSC_SUCCESS);
1177: }

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

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

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

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

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

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

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

1210:   PetscOptionsHeadEnd();
1211:   PetscFunctionReturn(PETSC_SUCCESS);
1212: }

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

1220:   PetscFunctionBegin;
1221:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1222:   if (iascii) {
1223:     PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE ParaSails preconditioning\n"));
1224:     PetscCall(PetscViewerASCIIPrintf(viewer, "    nlevels %" PetscInt_FMT "\n", jac->nlevels));
1225:     PetscCall(PetscViewerASCIIPrintf(viewer, "    threshold %g\n", (double)jac->threshold));
1226:     PetscCall(PetscViewerASCIIPrintf(viewer, "    filter %g\n", (double)jac->filter));
1227:     PetscCall(PetscViewerASCIIPrintf(viewer, "    load balance %g\n", (double)jac->loadbal));
1228:     PetscCall(PetscViewerASCIIPrintf(viewer, "    reuse nonzero structure %s\n", PetscBools[jac->ruse]));
1229:     PetscCall(PetscViewerASCIIPrintf(viewer, "    print info to screen %s\n", PetscBools[jac->logging]));
1230:     if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
1231:     else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
1232:     else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
1233:     else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Unknown HYPRE ParaSails symmetric option %" PetscInt_FMT, jac->symt);
1234:     PetscCall(PetscViewerASCIIPrintf(viewer, "    %s\n", symt));
1235:   }
1236:   PetscFunctionReturn(PETSC_SUCCESS);
1237: }

1239: static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PC pc, PetscOptionItems *PetscOptionsObject)
1240: {
1241:   PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1242:   PetscInt  n;
1243:   PetscBool flag, flag2, flag3, flag4;

1245:   PetscFunctionBegin;
1246:   PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE AMS Options");
1247:   PetscCall(PetscOptionsInt("-pc_hypre_ams_print_level", "Debugging output level for AMS", "None", jac->as_print, &jac->as_print, &flag));
1248:   if (flag) PetscCallExternal(HYPRE_AMSSetPrintLevel, jac->hsolver, jac->as_print);
1249:   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));
1250:   if (flag) PetscCallExternal(HYPRE_AMSSetMaxIter, jac->hsolver, jac->as_max_iter);
1251:   PetscCall(PetscOptionsInt("-pc_hypre_ams_cycle_type", "Cycle type for AMS multigrid", "None", jac->ams_cycle_type, &jac->ams_cycle_type, &flag));
1252:   if (flag) PetscCallExternal(HYPRE_AMSSetCycleType, jac->hsolver, jac->ams_cycle_type);
1253:   PetscCall(PetscOptionsReal("-pc_hypre_ams_tol", "Error tolerance for AMS multigrid", "None", jac->as_tol, &jac->as_tol, &flag));
1254:   if (flag) PetscCallExternal(HYPRE_AMSSetTol, jac->hsolver, jac->as_tol);
1255:   PetscCall(PetscOptionsInt("-pc_hypre_ams_relax_type", "Relaxation type for AMS smoother", "None", jac->as_relax_type, &jac->as_relax_type, &flag));
1256:   PetscCall(PetscOptionsInt("-pc_hypre_ams_relax_times", "Number of relaxation steps for AMS smoother", "None", jac->as_relax_times, &jac->as_relax_times, &flag2));
1257:   PetscCall(PetscOptionsReal("-pc_hypre_ams_relax_weight", "Relaxation weight for AMS smoother", "None", jac->as_relax_weight, &jac->as_relax_weight, &flag3));
1258:   PetscCall(PetscOptionsReal("-pc_hypre_ams_omega", "SSOR coefficient for AMS smoother", "None", jac->as_omega, &jac->as_omega, &flag4));
1259:   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);
1260:   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));
1261:   n = 5;
1262:   PetscCall(PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options", "AMG options for vector Poisson", "None", jac->as_amg_alpha_opts, &n, &flag2));
1263:   if (flag || flag2) {
1264:     PetscCallExternal(HYPRE_AMSSetAlphaAMGOptions, jac->hsolver, jac->as_amg_alpha_opts[0], /* AMG coarsen type */
1265:                       jac->as_amg_alpha_opts[1],                                            /* AMG agg_levels */
1266:                       jac->as_amg_alpha_opts[2],                                            /* AMG relax_type */
1267:                       jac->as_amg_alpha_theta, jac->as_amg_alpha_opts[3],                   /* AMG interp_type */
1268:                       jac->as_amg_alpha_opts[4]);                                           /* AMG Pmax */
1269:   }
1270:   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));
1271:   n = 5;
1272:   PetscCall(PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options", "AMG options for scalar Poisson solver", "None", jac->as_amg_beta_opts, &n, &flag2));
1273:   if (flag || flag2) {
1274:     PetscCallExternal(HYPRE_AMSSetBetaAMGOptions, jac->hsolver, jac->as_amg_beta_opts[0], /* AMG coarsen type */
1275:                       jac->as_amg_beta_opts[1],                                           /* AMG agg_levels */
1276:                       jac->as_amg_beta_opts[2],                                           /* AMG relax_type */
1277:                       jac->as_amg_beta_theta, jac->as_amg_beta_opts[3],                   /* AMG interp_type */
1278:                       jac->as_amg_beta_opts[4]);                                          /* AMG Pmax */
1279:   }
1280:   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));
1281:   if (flag) { /* override HYPRE's default only if the options is used */
1282:     PetscCallExternal(HYPRE_AMSSetProjectionFrequency, jac->hsolver, jac->ams_proj_freq);
1283:   }
1284:   PetscOptionsHeadEnd();
1285:   PetscFunctionReturn(PETSC_SUCCESS);
1286: }

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

1293:   PetscFunctionBegin;
1294:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1295:   if (iascii) {
1296:     PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE AMS preconditioning\n"));
1297:     PetscCall(PetscViewerASCIIPrintf(viewer, "    subspace iterations per application %" PetscInt_FMT "\n", jac->as_max_iter));
1298:     PetscCall(PetscViewerASCIIPrintf(viewer, "    subspace cycle type %" PetscInt_FMT "\n", jac->ams_cycle_type));
1299:     PetscCall(PetscViewerASCIIPrintf(viewer, "    subspace iteration tolerance %g\n", (double)jac->as_tol));
1300:     PetscCall(PetscViewerASCIIPrintf(viewer, "    smoother type %" PetscInt_FMT "\n", jac->as_relax_type));
1301:     PetscCall(PetscViewerASCIIPrintf(viewer, "    number of smoothing steps %" PetscInt_FMT "\n", jac->as_relax_times));
1302:     PetscCall(PetscViewerASCIIPrintf(viewer, "    smoother weight %g\n", (double)jac->as_relax_weight));
1303:     PetscCall(PetscViewerASCIIPrintf(viewer, "    smoother omega %g\n", (double)jac->as_omega));
1304:     if (jac->alpha_Poisson) {
1305:       PetscCall(PetscViewerASCIIPrintf(viewer, "    vector Poisson solver (passed in by user)\n"));
1306:     } else {
1307:       PetscCall(PetscViewerASCIIPrintf(viewer, "    vector Poisson solver (computed) \n"));
1308:     }
1309:     PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG coarsening type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[0]));
1310:     PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[1]));
1311:     PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG relaxation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[2]));
1312:     PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG interpolation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[3]));
1313:     PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[4]));
1314:     PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG strength threshold %g\n", (double)jac->as_amg_alpha_theta));
1315:     if (!jac->ams_beta_is_zero) {
1316:       if (jac->beta_Poisson) {
1317:         PetscCall(PetscViewerASCIIPrintf(viewer, "    scalar Poisson solver (passed in by user)\n"));
1318:       } else {
1319:         PetscCall(PetscViewerASCIIPrintf(viewer, "    scalar Poisson solver (computed) \n"));
1320:       }
1321:       PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG coarsening type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[0]));
1322:       PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_beta_opts[1]));
1323:       PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG relaxation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[2]));
1324:       PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG interpolation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[3]));
1325:       PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_beta_opts[4]));
1326:       PetscCall(PetscViewerASCIIPrintf(viewer, "        boomerAMG strength threshold %g\n", (double)jac->as_amg_beta_theta));
1327:       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));
1328:     } else {
1329:       PetscCall(PetscViewerASCIIPrintf(viewer, "    scalar Poisson solver not used (zero-conductivity everywhere) \n"));
1330:     }
1331:   }
1332:   PetscFunctionReturn(PETSC_SUCCESS);
1333: }

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

1341:   PetscFunctionBegin;
1342:   PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE ADS Options");
1343:   PetscCall(PetscOptionsInt("-pc_hypre_ads_print_level", "Debugging output level for ADS", "None", jac->as_print, &jac->as_print, &flag));
1344:   if (flag) PetscCallExternal(HYPRE_ADSSetPrintLevel, jac->hsolver, jac->as_print);
1345:   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));
1346:   if (flag) PetscCallExternal(HYPRE_ADSSetMaxIter, jac->hsolver, jac->as_max_iter);
1347:   PetscCall(PetscOptionsInt("-pc_hypre_ads_cycle_type", "Cycle type for ADS multigrid", "None", jac->ads_cycle_type, &jac->ads_cycle_type, &flag));
1348:   if (flag) PetscCallExternal(HYPRE_ADSSetCycleType, jac->hsolver, jac->ads_cycle_type);
1349:   PetscCall(PetscOptionsReal("-pc_hypre_ads_tol", "Error tolerance for ADS multigrid", "None", jac->as_tol, &jac->as_tol, &flag));
1350:   if (flag) PetscCallExternal(HYPRE_ADSSetTol, jac->hsolver, jac->as_tol);
1351:   PetscCall(PetscOptionsInt("-pc_hypre_ads_relax_type", "Relaxation type for ADS smoother", "None", jac->as_relax_type, &jac->as_relax_type, &flag));
1352:   PetscCall(PetscOptionsInt("-pc_hypre_ads_relax_times", "Number of relaxation steps for ADS smoother", "None", jac->as_relax_times, &jac->as_relax_times, &flag2));
1353:   PetscCall(PetscOptionsReal("-pc_hypre_ads_relax_weight", "Relaxation weight for ADS smoother", "None", jac->as_relax_weight, &jac->as_relax_weight, &flag3));
1354:   PetscCall(PetscOptionsReal("-pc_hypre_ads_omega", "SSOR coefficient for ADS smoother", "None", jac->as_omega, &jac->as_omega, &flag4));
1355:   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);
1356:   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));
1357:   n = 5;
1358:   PetscCall(PetscOptionsIntArray("-pc_hypre_ads_ams_options", "AMG options for AMS solver inside ADS", "None", jac->as_amg_alpha_opts, &n, &flag2));
1359:   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));
1360:   if (flag || flag2 || flag3) {
1361:     PetscCallExternal(HYPRE_ADSSetAMSOptions, jac->hsolver, jac->ams_cycle_type, /* AMS cycle type */
1362:                       jac->as_amg_alpha_opts[0],                                 /* AMG coarsen type */
1363:                       jac->as_amg_alpha_opts[1],                                 /* AMG agg_levels */
1364:                       jac->as_amg_alpha_opts[2],                                 /* AMG relax_type */
1365:                       jac->as_amg_alpha_theta, jac->as_amg_alpha_opts[3],        /* AMG interp_type */
1366:                       jac->as_amg_alpha_opts[4]);                                /* AMG Pmax */
1367:   }
1368:   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));
1369:   n = 5;
1370:   PetscCall(PetscOptionsIntArray("-pc_hypre_ads_amg_options", "AMG options for vector AMG solver inside ADS", "None", jac->as_amg_beta_opts, &n, &flag2));
1371:   if (flag || flag2) {
1372:     PetscCallExternal(HYPRE_ADSSetAMGOptions, jac->hsolver, jac->as_amg_beta_opts[0], /* AMG coarsen type */
1373:                       jac->as_amg_beta_opts[1],                                       /* AMG agg_levels */
1374:                       jac->as_amg_beta_opts[2],                                       /* AMG relax_type */
1375:                       jac->as_amg_beta_theta, jac->as_amg_beta_opts[3],               /* AMG interp_type */
1376:                       jac->as_amg_beta_opts[4]);                                      /* AMG Pmax */
1377:   }
1378:   PetscOptionsHeadEnd();
1379:   PetscFunctionReturn(PETSC_SUCCESS);
1380: }

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

1387:   PetscFunctionBegin;
1388:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1389:   if (iascii) {
1390:     PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE ADS preconditioning\n"));
1391:     PetscCall(PetscViewerASCIIPrintf(viewer, "    subspace iterations per application %" PetscInt_FMT "\n", jac->as_max_iter));
1392:     PetscCall(PetscViewerASCIIPrintf(viewer, "    subspace cycle type %" PetscInt_FMT "\n", jac->ads_cycle_type));
1393:     PetscCall(PetscViewerASCIIPrintf(viewer, "    subspace iteration tolerance %g\n", (double)jac->as_tol));
1394:     PetscCall(PetscViewerASCIIPrintf(viewer, "    smoother type %" PetscInt_FMT "\n", jac->as_relax_type));
1395:     PetscCall(PetscViewerASCIIPrintf(viewer, "    number of smoothing steps %" PetscInt_FMT "\n", jac->as_relax_times));
1396:     PetscCall(PetscViewerASCIIPrintf(viewer, "    smoother weight %g\n", (double)jac->as_relax_weight));
1397:     PetscCall(PetscViewerASCIIPrintf(viewer, "    smoother omega %g\n", (double)jac->as_omega));
1398:     PetscCall(PetscViewerASCIIPrintf(viewer, "    AMS solver using boomerAMG\n"));
1399:     PetscCall(PetscViewerASCIIPrintf(viewer, "        subspace cycle type %" PetscInt_FMT "\n", jac->ams_cycle_type));
1400:     PetscCall(PetscViewerASCIIPrintf(viewer, "        coarsening type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[0]));
1401:     PetscCall(PetscViewerASCIIPrintf(viewer, "        levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[1]));
1402:     PetscCall(PetscViewerASCIIPrintf(viewer, "        relaxation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[2]));
1403:     PetscCall(PetscViewerASCIIPrintf(viewer, "        interpolation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[3]));
1404:     PetscCall(PetscViewerASCIIPrintf(viewer, "        max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[4]));
1405:     PetscCall(PetscViewerASCIIPrintf(viewer, "        strength threshold %g\n", (double)jac->as_amg_alpha_theta));
1406:     PetscCall(PetscViewerASCIIPrintf(viewer, "    vector Poisson solver using boomerAMG\n"));
1407:     PetscCall(PetscViewerASCIIPrintf(viewer, "        coarsening type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[0]));
1408:     PetscCall(PetscViewerASCIIPrintf(viewer, "        levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_beta_opts[1]));
1409:     PetscCall(PetscViewerASCIIPrintf(viewer, "        relaxation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[2]));
1410:     PetscCall(PetscViewerASCIIPrintf(viewer, "        interpolation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[3]));
1411:     PetscCall(PetscViewerASCIIPrintf(viewer, "        max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_beta_opts[4]));
1412:     PetscCall(PetscViewerASCIIPrintf(viewer, "        strength threshold %g\n", (double)jac->as_amg_beta_theta));
1413:   }
1414:   PetscFunctionReturn(PETSC_SUCCESS);
1415: }

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

1422:   PetscFunctionBegin;
1423:   PetscCall(PetscObjectTypeCompare((PetscObject)G, MATHYPRE, &ishypre));
1424:   if (ishypre) {
1425:     PetscCall(PetscObjectReference((PetscObject)G));
1426:     PetscCall(MatDestroy(&jac->G));
1427:     jac->G = G;
1428:   } else {
1429:     PetscCall(MatDestroy(&jac->G));
1430:     PetscCall(MatConvert(G, MATHYPRE, MAT_INITIAL_MATRIX, &jac->G));
1431:   }
1432:   PetscFunctionReturn(PETSC_SUCCESS);
1433: }

1435: /*@
1436:   PCHYPRESetDiscreteGradient - Set discrete gradient matrix for `PCHYPRE` type of ams or ads

1438:   Collective

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

1444:   Level: intermediate

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

1449:   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

1451:   Developer Notes:
1452:   This automatically converts the matrix to `MATHYPRE` if it is not already of that type

1454: .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteCurl()`
1455: @*/
1456: PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G)
1457: {
1458:   PetscFunctionBegin;
1461:   PetscCheckSameComm(pc, 1, G, 2);
1462:   PetscTryMethod(pc, "PCHYPRESetDiscreteGradient_C", (PC, Mat), (pc, G));
1463:   PetscFunctionReturn(PETSC_SUCCESS);
1464: }

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

1471:   PetscFunctionBegin;
1472:   PetscCall(PetscObjectTypeCompare((PetscObject)C, MATHYPRE, &ishypre));
1473:   if (ishypre) {
1474:     PetscCall(PetscObjectReference((PetscObject)C));
1475:     PetscCall(MatDestroy(&jac->C));
1476:     jac->C = C;
1477:   } else {
1478:     PetscCall(MatDestroy(&jac->C));
1479:     PetscCall(MatConvert(C, MATHYPRE, MAT_INITIAL_MATRIX, &jac->C));
1480:   }
1481:   PetscFunctionReturn(PETSC_SUCCESS);
1482: }

1484: /*@
1485:   PCHYPRESetDiscreteCurl - Set discrete curl matrx for `PCHYPRE` type of ads

1487:   Collective

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

1493:   Level: intermediate

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

1498:   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

1500:   Developer Notes:
1501:   This automatically converts the matrix to `MATHYPRE` if it is not already of that type

1503:   If this is only for  `PCHYPRE` type of ads it should be called `PCHYPREADSSetDiscreteCurl()`

1505: .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`
1506: @*/
1507: PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C)
1508: {
1509:   PetscFunctionBegin;
1512:   PetscCheckSameComm(pc, 1, C, 2);
1513:   PetscTryMethod(pc, "PCHYPRESetDiscreteCurl_C", (PC, Mat), (pc, C));
1514:   PetscFunctionReturn(PETSC_SUCCESS);
1515: }

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

1523:   PetscFunctionBegin;
1524:   PetscCall(MatDestroy(&jac->RT_PiFull));
1525:   PetscCall(MatDestroy(&jac->ND_PiFull));
1526:   for (i = 0; i < 3; ++i) {
1527:     PetscCall(MatDestroy(&jac->RT_Pi[i]));
1528:     PetscCall(MatDestroy(&jac->ND_Pi[i]));
1529:   }

1531:   jac->dim = dim;
1532:   if (RT_PiFull) {
1533:     PetscCall(PetscObjectTypeCompare((PetscObject)RT_PiFull, MATHYPRE, &ishypre));
1534:     if (ishypre) {
1535:       PetscCall(PetscObjectReference((PetscObject)RT_PiFull));
1536:       jac->RT_PiFull = RT_PiFull;
1537:     } else {
1538:       PetscCall(MatConvert(RT_PiFull, MATHYPRE, MAT_INITIAL_MATRIX, &jac->RT_PiFull));
1539:     }
1540:   }
1541:   if (RT_Pi) {
1542:     for (i = 0; i < dim; ++i) {
1543:       if (RT_Pi[i]) {
1544:         PetscCall(PetscObjectTypeCompare((PetscObject)RT_Pi[i], MATHYPRE, &ishypre));
1545:         if (ishypre) {
1546:           PetscCall(PetscObjectReference((PetscObject)RT_Pi[i]));
1547:           jac->RT_Pi[i] = RT_Pi[i];
1548:         } else {
1549:           PetscCall(MatConvert(RT_Pi[i], MATHYPRE, MAT_INITIAL_MATRIX, &jac->RT_Pi[i]));
1550:         }
1551:       }
1552:     }
1553:   }
1554:   if (ND_PiFull) {
1555:     PetscCall(PetscObjectTypeCompare((PetscObject)ND_PiFull, MATHYPRE, &ishypre));
1556:     if (ishypre) {
1557:       PetscCall(PetscObjectReference((PetscObject)ND_PiFull));
1558:       jac->ND_PiFull = ND_PiFull;
1559:     } else {
1560:       PetscCall(MatConvert(ND_PiFull, MATHYPRE, MAT_INITIAL_MATRIX, &jac->ND_PiFull));
1561:     }
1562:   }
1563:   if (ND_Pi) {
1564:     for (i = 0; i < dim; ++i) {
1565:       if (ND_Pi[i]) {
1566:         PetscCall(PetscObjectTypeCompare((PetscObject)ND_Pi[i], MATHYPRE, &ishypre));
1567:         if (ishypre) {
1568:           PetscCall(PetscObjectReference((PetscObject)ND_Pi[i]));
1569:           jac->ND_Pi[i] = ND_Pi[i];
1570:         } else {
1571:           PetscCall(MatConvert(ND_Pi[i], MATHYPRE, MAT_INITIAL_MATRIX, &jac->ND_Pi[i]));
1572:         }
1573:       }
1574:     }
1575:   }
1576:   PetscFunctionReturn(PETSC_SUCCESS);
1577: }

1579: /*@
1580:   PCHYPRESetInterpolations - Set interpolation matrices for `PCHYPRE` type of ams or ads

1582:   Collective

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

1592:   Level: intermediate

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

1597:   For ADS, both type of interpolation matrices are needed.

1599:   Developer Notes:
1600:   This automatically converts the matrix to `MATHYPRE` if it is not already of that type

1602: .seealso: [](ch_ksp), `PCHYPRE`
1603: @*/
1604: PetscErrorCode PCHYPRESetInterpolations(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[])
1605: {
1606:   PetscInt i;

1608:   PetscFunctionBegin;
1610:   if (RT_PiFull) {
1612:     PetscCheckSameComm(pc, 1, RT_PiFull, 3);
1613:   }
1614:   if (RT_Pi) {
1615:     PetscAssertPointer(RT_Pi, 4);
1616:     for (i = 0; i < dim; ++i) {
1617:       if (RT_Pi[i]) {
1619:         PetscCheckSameComm(pc, 1, RT_Pi[i], 4);
1620:       }
1621:     }
1622:   }
1623:   if (ND_PiFull) {
1625:     PetscCheckSameComm(pc, 1, ND_PiFull, 5);
1626:   }
1627:   if (ND_Pi) {
1628:     PetscAssertPointer(ND_Pi, 6);
1629:     for (i = 0; i < dim; ++i) {
1630:       if (ND_Pi[i]) {
1632:         PetscCheckSameComm(pc, 1, ND_Pi[i], 6);
1633:       }
1634:     }
1635:   }
1636:   PetscTryMethod(pc, "PCHYPRESetInterpolations_C", (PC, PetscInt, Mat, Mat[], Mat, Mat[]), (pc, dim, RT_PiFull, RT_Pi, ND_PiFull, ND_Pi));
1637:   PetscFunctionReturn(PETSC_SUCCESS);
1638: }

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

1645:   PetscFunctionBegin;
1646:   PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre));
1647:   if (ishypre) {
1648:     if (isalpha) {
1649:       PetscCall(PetscObjectReference((PetscObject)A));
1650:       PetscCall(MatDestroy(&jac->alpha_Poisson));
1651:       jac->alpha_Poisson = A;
1652:     } else {
1653:       if (A) {
1654:         PetscCall(PetscObjectReference((PetscObject)A));
1655:       } else {
1656:         jac->ams_beta_is_zero = PETSC_TRUE;
1657:       }
1658:       PetscCall(MatDestroy(&jac->beta_Poisson));
1659:       jac->beta_Poisson = A;
1660:     }
1661:   } else {
1662:     if (isalpha) {
1663:       PetscCall(MatDestroy(&jac->alpha_Poisson));
1664:       PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &jac->alpha_Poisson));
1665:     } else {
1666:       if (A) {
1667:         PetscCall(MatDestroy(&jac->beta_Poisson));
1668:         PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &jac->beta_Poisson));
1669:       } else {
1670:         PetscCall(MatDestroy(&jac->beta_Poisson));
1671:         jac->ams_beta_is_zero = PETSC_TRUE;
1672:       }
1673:     }
1674:   }
1675:   PetscFunctionReturn(PETSC_SUCCESS);
1676: }

1678: /*@
1679:   PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix for `PCHYPRE` of type ams

1681:   Collective

1683:   Input Parameters:
1684: + pc - the preconditioning context
1685: - A  - the matrix

1687:   Level: intermediate

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

1692:   Developer Notes:
1693:   This automatically converts the matrix to `MATHYPRE` if it is not already of that type

1695:   If this is only for  `PCHYPRE` type of ams it should be called `PCHYPREAMSSetAlphaPoissonMatrix()`

1697: .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetBetaPoissonMatrix()`
1698: @*/
1699: PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A)
1700: {
1701:   PetscFunctionBegin;
1704:   PetscCheckSameComm(pc, 1, A, 2);
1705:   PetscTryMethod(pc, "PCHYPRESetPoissonMatrix_C", (PC, Mat, PetscBool), (pc, A, PETSC_TRUE));
1706:   PetscFunctionReturn(PETSC_SUCCESS);
1707: }

1709: /*@
1710:   PCHYPRESetBetaPoissonMatrix - Set Poisson matrix for `PCHYPRE` of type ams

1712:   Collective

1714:   Input Parameters:
1715: + pc - the preconditioning context
1716: - A  - the matrix, or NULL to turn it off

1718:   Level: intermediate

1720:   Note:
1721:   A should be obtained by discretizing the Poisson problem with linear finite elements.

1723:   Developer Notes:
1724:   This automatically converts the matrix to `MATHYPRE` if it is not already of that type

1726:   If this is only for  `PCHYPRE` type of ams it should be called `PCHYPREAMSPCHYPRESetBetaPoissonMatrix()`

1728: .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetAlphaPoissonMatrix()`
1729: @*/
1730: PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A)
1731: {
1732:   PetscFunctionBegin;
1734:   if (A) {
1736:     PetscCheckSameComm(pc, 1, A, 2);
1737:   }
1738:   PetscTryMethod(pc, "PCHYPRESetPoissonMatrix_C", (PC, Mat, PetscBool), (pc, A, PETSC_FALSE));
1739:   PetscFunctionReturn(PETSC_SUCCESS);
1740: }

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

1746:   PetscFunctionBegin;
1747:   /* throw away any vector if already set */
1748:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[0]));
1749:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[1]));
1750:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[2]));
1751:   PetscCall(VecHYPRE_IJVectorCreate(ozz->map, &jac->constants[0]));
1752:   PetscCall(VecHYPRE_IJVectorCopy(ozz, jac->constants[0]));
1753:   PetscCall(VecHYPRE_IJVectorCreate(zoz->map, &jac->constants[1]));
1754:   PetscCall(VecHYPRE_IJVectorCopy(zoz, jac->constants[1]));
1755:   jac->dim = 2;
1756:   if (zzo) {
1757:     PetscCall(VecHYPRE_IJVectorCreate(zzo->map, &jac->constants[2]));
1758:     PetscCall(VecHYPRE_IJVectorCopy(zzo, jac->constants[2]));
1759:     jac->dim++;
1760:   }
1761:   PetscFunctionReturn(PETSC_SUCCESS);
1762: }

1764: /*@
1765:   PCHYPRESetEdgeConstantVectors - Set the representation of the constant vector fields in the edge element basis for `PCHYPRE` of type ams

1767:   Collective

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

1775:   Level: intermediate

1777:   Developer Notes:
1778:   If this is only for  `PCHYPRE` type of ams it should be called `PCHYPREAMSSetEdgeConstantVectors()`

1780: .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetAlphaPoissonMatrix()`
1781: @*/
1782: PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo)
1783: {
1784:   PetscFunctionBegin;
1789:   PetscCheckSameComm(pc, 1, ozz, 2);
1790:   PetscCheckSameComm(pc, 1, zoz, 3);
1791:   if (zzo) PetscCheckSameComm(pc, 1, zzo, 4);
1792:   PetscTryMethod(pc, "PCHYPRESetEdgeConstantVectors_C", (PC, Vec, Vec, Vec), (pc, ozz, zoz, zzo));
1793:   PetscFunctionReturn(PETSC_SUCCESS);
1794: }

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

1800:   PetscFunctionBegin;
1801:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->interior));
1802:   PetscCall(VecHYPRE_IJVectorCreate(interior->map, &jac->interior));
1803:   PetscCall(VecHYPRE_IJVectorCopy(interior, jac->interior));
1804:   jac->ams_beta_is_zero_part = PETSC_TRUE;
1805:   PetscFunctionReturn(PETSC_SUCCESS);
1806: }

1808: /*@
1809:   PCHYPREAMSSetInteriorNodes - Set the list of interior nodes to a zero-conductivity region for `PCHYPRE` of type ams

1811:   Collective

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

1817:   Level: intermediate

1819:   Note:
1820:   This calls `HYPRE_AMSSetInteriorNodes()`

1822:   Developer Notes:
1823:   If this is only for  `PCHYPRE` type of ams it should be called `PCHYPREAMSSetInteriorNodes()`

1825: .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetAlphaPoissonMatrix()`
1826: @*/
1827: PetscErrorCode PCHYPREAMSSetInteriorNodes(PC pc, Vec interior)
1828: {
1829:   PetscFunctionBegin;
1832:   PetscCheckSameComm(pc, 1, interior, 2);
1833:   PetscTryMethod(pc, "PCHYPREAMSSetInteriorNodes_C", (PC, Vec), (pc, interior));
1834:   PetscFunctionReturn(PETSC_SUCCESS);
1835: }

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

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

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

1858:     PetscCall(VecHYPRE_IJVectorCreate(tv->map, &jac->coords[i]));
1859:     PetscCall(VecGetArrayWrite(tv, &array));
1860:     for (j = 0; j < nloc; j++) array[j] = coords[j * dim + i];
1861:     PetscCall(VecRestoreArrayWrite(tv, &array));
1862:     PetscCall(VecHYPRE_IJVectorCopy(tv, jac->coords[i]));
1863:   }
1864:   PetscCall(VecDestroy(&tv));
1865:   PetscFunctionReturn(PETSC_SUCCESS);
1866: }

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

1872:   PetscFunctionBegin;
1873:   *name = jac->hypre_type;
1874:   PetscFunctionReturn(PETSC_SUCCESS);
1875: }

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

1882:   PetscFunctionBegin;
1883:   if (jac->hypre_type) {
1884:     PetscCall(PetscStrcmp(jac->hypre_type, name, &flag));
1885:     PetscCheck(flag, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Cannot reset the HYPRE preconditioner type once it has been set");
1886:     PetscFunctionReturn(PETSC_SUCCESS);
1887:   } else {
1888:     PetscCall(PetscStrallocpy(name, &jac->hypre_type));
1889:   }

1891:   jac->maxiter         = PETSC_DEFAULT;
1892:   jac->tol             = PETSC_DEFAULT;
1893:   jac->printstatistics = PetscLogPrintInfo;

1895:   PetscCall(PetscStrcmp("pilut", jac->hypre_type, &flag));
1896:   if (flag) {
1897:     PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
1898:     PetscCallExternal(HYPRE_ParCSRPilutCreate, jac->comm_hypre, &jac->hsolver);
1899:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
1900:     pc->ops->view           = PCView_HYPRE_Pilut;
1901:     jac->destroy            = HYPRE_ParCSRPilutDestroy;
1902:     jac->setup              = HYPRE_ParCSRPilutSetup;
1903:     jac->solve              = HYPRE_ParCSRPilutSolve;
1904:     jac->factorrowsize      = PETSC_DEFAULT;
1905:     PetscFunctionReturn(PETSC_SUCCESS);
1906:   }
1907:   PetscCall(PetscStrcmp("euclid", jac->hypre_type, &flag));
1908:   if (flag) {
1909: #if defined(PETSC_USE_64BIT_INDICES)
1910:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Hypre Euclid does not support 64-bit indices");
1911: #endif
1912:     PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
1913:     PetscCallExternal(HYPRE_EuclidCreate, jac->comm_hypre, &jac->hsolver);
1914:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
1915:     pc->ops->view           = PCView_HYPRE_Euclid;
1916:     jac->destroy            = HYPRE_EuclidDestroy;
1917:     jac->setup              = HYPRE_EuclidSetup;
1918:     jac->solve              = HYPRE_EuclidSolve;
1919:     jac->factorrowsize      = PETSC_DEFAULT;
1920:     jac->eu_level           = PETSC_DEFAULT; /* default */
1921:     PetscFunctionReturn(PETSC_SUCCESS);
1922:   }
1923:   PetscCall(PetscStrcmp("parasails", jac->hypre_type, &flag));
1924:   if (flag) {
1925:     PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
1926:     PetscCallExternal(HYPRE_ParaSailsCreate, jac->comm_hypre, &jac->hsolver);
1927:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
1928:     pc->ops->view           = PCView_HYPRE_ParaSails;
1929:     jac->destroy            = HYPRE_ParaSailsDestroy;
1930:     jac->setup              = HYPRE_ParaSailsSetup;
1931:     jac->solve              = HYPRE_ParaSailsSolve;
1932:     /* initialize */
1933:     jac->nlevels   = 1;
1934:     jac->threshold = .1;
1935:     jac->filter    = .1;
1936:     jac->loadbal   = 0;
1937:     if (PetscLogPrintInfo) jac->logging = (int)PETSC_TRUE;
1938:     else jac->logging = (int)PETSC_FALSE;

1940:     jac->ruse = (int)PETSC_FALSE;
1941:     jac->symt = 0;
1942:     PetscCallExternal(HYPRE_ParaSailsSetParams, jac->hsolver, jac->threshold, jac->nlevels);
1943:     PetscCallExternal(HYPRE_ParaSailsSetFilter, jac->hsolver, jac->filter);
1944:     PetscCallExternal(HYPRE_ParaSailsSetLoadbal, jac->hsolver, jac->loadbal);
1945:     PetscCallExternal(HYPRE_ParaSailsSetLogging, jac->hsolver, jac->logging);
1946:     PetscCallExternal(HYPRE_ParaSailsSetReuse, jac->hsolver, jac->ruse);
1947:     PetscCallExternal(HYPRE_ParaSailsSetSym, jac->hsolver, jac->symt);
1948:     PetscFunctionReturn(PETSC_SUCCESS);
1949:   }
1950:   PetscCall(PetscStrcmp("boomeramg", jac->hypre_type, &flag));
1951:   if (flag) {
1952:     PetscCallExternal(HYPRE_BoomerAMGCreate, &jac->hsolver);
1953:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_BoomerAMG;
1954:     pc->ops->view            = PCView_HYPRE_BoomerAMG;
1955:     pc->ops->applytranspose  = PCApplyTranspose_HYPRE_BoomerAMG;
1956:     pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
1957:     pc->ops->matapply        = PCMatApply_HYPRE_BoomerAMG;
1958:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", PCGetInterpolations_BoomerAMG));
1959:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", PCGetCoarseOperators_BoomerAMG));
1960:     jac->destroy         = HYPRE_BoomerAMGDestroy;
1961:     jac->setup           = HYPRE_BoomerAMGSetup;
1962:     jac->solve           = HYPRE_BoomerAMGSolve;
1963:     jac->applyrichardson = PETSC_FALSE;
1964:     /* these defaults match the hypre defaults */
1965:     jac->cycletype       = 1;
1966:     jac->maxlevels       = 25;
1967:     jac->maxiter         = 1;
1968:     jac->tol             = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
1969:     jac->truncfactor     = 0.0;
1970:     jac->strongthreshold = .25;
1971:     jac->maxrowsum       = .9;
1972:     jac->coarsentype     = 6;
1973:     jac->measuretype     = 0;
1974:     jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
1975:     jac->smoothtype                                              = -1; /* Not set by default */
1976:     jac->smoothnumlevels                                         = 25;
1977:     jac->eu_level                                                = 0;
1978:     jac->eu_droptolerance                                        = 0;
1979:     jac->eu_bj                                                   = 0;
1980:     jac->relaxtype[0] = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a PC - most likely with CG */
1981:     jac->relaxtype[2]                     = 9; /*G.E. */
1982:     jac->relaxweight                      = 1.0;
1983:     jac->outerrelaxweight                 = 1.0;
1984:     jac->relaxorder                       = 1;
1985:     jac->interptype                       = 0;
1986:     jac->Rtype                            = 0;
1987:     jac->Rstrongthreshold                 = 0.25;
1988:     jac->Rfilterthreshold                 = 0.0;
1989:     jac->Adroptype                        = -1;
1990:     jac->Adroptol                         = 0.0;
1991:     jac->agg_nl                           = 0;
1992:     jac->agg_interptype                   = 4;
1993:     jac->pmax                             = 0;
1994:     jac->truncfactor                      = 0.0;
1995:     jac->agg_num_paths                    = 1;
1996:     jac->maxc                             = 9;
1997:     jac->minc                             = 1;
1998:     jac->nodal_coarsening                 = 0;
1999:     jac->nodal_coarsening_diag            = 0;
2000:     jac->vec_interp_variant               = 0;
2001:     jac->vec_interp_qmax                  = 0;
2002:     jac->vec_interp_smooth                = PETSC_FALSE;
2003:     jac->interp_refine                    = 0;
2004:     jac->nodal_relax                      = PETSC_FALSE;
2005:     jac->nodal_relax_levels               = 1;
2006:     jac->rap2                             = 0;

2008:     /* GPU defaults
2009:          from https://hypre.readthedocs.io/en/latest/solvers-boomeramg.html#gpu-supported-options
2010:          and /src/parcsr_ls/par_amg.c */
2011: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2012:     jac->keeptranspose  = PETSC_TRUE;
2013:     jac->mod_rap2       = 1;
2014:     jac->coarsentype    = 8;
2015:     jac->relaxorder     = 0;
2016:     jac->interptype     = 6;
2017:     jac->relaxtype[0]   = 18;
2018:     jac->relaxtype[1]   = 18;
2019:     jac->agg_interptype = 7;
2020: #else
2021:     jac->keeptranspose = PETSC_FALSE;
2022:     jac->mod_rap2      = 0;
2023: #endif
2024:     PetscCallExternal(HYPRE_BoomerAMGSetCycleType, jac->hsolver, jac->cycletype);
2025:     PetscCallExternal(HYPRE_BoomerAMGSetMaxLevels, jac->hsolver, jac->maxlevels);
2026:     PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, jac->maxiter);
2027:     PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, jac->tol);
2028:     PetscCallExternal(HYPRE_BoomerAMGSetTruncFactor, jac->hsolver, jac->truncfactor);
2029:     PetscCallExternal(HYPRE_BoomerAMGSetStrongThreshold, jac->hsolver, jac->strongthreshold);
2030:     PetscCallExternal(HYPRE_BoomerAMGSetMaxRowSum, jac->hsolver, jac->maxrowsum);
2031:     PetscCallExternal(HYPRE_BoomerAMGSetCoarsenType, jac->hsolver, jac->coarsentype);
2032:     PetscCallExternal(HYPRE_BoomerAMGSetMeasureType, jac->hsolver, jac->measuretype);
2033:     PetscCallExternal(HYPRE_BoomerAMGSetRelaxOrder, jac->hsolver, jac->relaxorder);
2034:     PetscCallExternal(HYPRE_BoomerAMGSetInterpType, jac->hsolver, jac->interptype);
2035:     PetscCallExternal(HYPRE_BoomerAMGSetAggNumLevels, jac->hsolver, jac->agg_nl);
2036:     PetscCallExternal(HYPRE_BoomerAMGSetAggInterpType, jac->hsolver, jac->agg_interptype);
2037:     PetscCallExternal(HYPRE_BoomerAMGSetPMaxElmts, jac->hsolver, jac->pmax);
2038:     PetscCallExternal(HYPRE_BoomerAMGSetNumPaths, jac->hsolver, jac->agg_num_paths);
2039:     PetscCallExternal(HYPRE_BoomerAMGSetRelaxType, jac->hsolver, jac->relaxtype[0]);  /* defaults coarse to 9 */
2040:     PetscCallExternal(HYPRE_BoomerAMGSetNumSweeps, jac->hsolver, jac->gridsweeps[0]); /* defaults coarse to 1 */
2041:     PetscCallExternal(HYPRE_BoomerAMGSetMaxCoarseSize, jac->hsolver, jac->maxc);
2042:     PetscCallExternal(HYPRE_BoomerAMGSetMinCoarseSize, jac->hsolver, jac->minc);
2043:     /* GPU */
2044: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
2045:     PetscCallExternal(HYPRE_BoomerAMGSetKeepTranspose, jac->hsolver, jac->keeptranspose ? 1 : 0);
2046:     PetscCallExternal(HYPRE_BoomerAMGSetRAP2, jac->hsolver, jac->rap2);
2047:     PetscCallExternal(HYPRE_BoomerAMGSetModuleRAP2, jac->hsolver, jac->mod_rap2);
2048: #endif

2050:     /* AIR */
2051: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
2052:     PetscCallExternal(HYPRE_BoomerAMGSetRestriction, jac->hsolver, jac->Rtype);
2053:     PetscCallExternal(HYPRE_BoomerAMGSetStrongThresholdR, jac->hsolver, jac->Rstrongthreshold);
2054:     PetscCallExternal(HYPRE_BoomerAMGSetFilterThresholdR, jac->hsolver, jac->Rfilterthreshold);
2055:     PetscCallExternal(HYPRE_BoomerAMGSetADropTol, jac->hsolver, jac->Adroptol);
2056:     PetscCallExternal(HYPRE_BoomerAMGSetADropType, jac->hsolver, jac->Adroptype);
2057: #endif
2058:     PetscFunctionReturn(PETSC_SUCCESS);
2059:   }
2060:   PetscCall(PetscStrcmp("ams", jac->hypre_type, &flag));
2061:   if (flag) {
2062:     PetscCallExternal(HYPRE_AMSCreate, &jac->hsolver);
2063:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_AMS;
2064:     pc->ops->view           = PCView_HYPRE_AMS;
2065:     jac->destroy            = HYPRE_AMSDestroy;
2066:     jac->setup              = HYPRE_AMSSetup;
2067:     jac->solve              = HYPRE_AMSSolve;
2068:     jac->coords[0]          = NULL;
2069:     jac->coords[1]          = NULL;
2070:     jac->coords[2]          = NULL;
2071:     jac->interior           = NULL;
2072:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */
2073:     jac->as_print       = 0;
2074:     jac->as_max_iter    = 1;  /* used as a preconditioner */
2075:     jac->as_tol         = 0.; /* used as a preconditioner */
2076:     jac->ams_cycle_type = 13;
2077:     /* Smoothing options */
2078:     jac->as_relax_type   = 2;
2079:     jac->as_relax_times  = 1;
2080:     jac->as_relax_weight = 1.0;
2081:     jac->as_omega        = 1.0;
2082:     /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
2083:     jac->as_amg_alpha_opts[0] = 10;
2084:     jac->as_amg_alpha_opts[1] = 1;
2085:     jac->as_amg_alpha_opts[2] = 6;
2086:     jac->as_amg_alpha_opts[3] = 6;
2087:     jac->as_amg_alpha_opts[4] = 4;
2088:     jac->as_amg_alpha_theta   = 0.25;
2089:     /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
2090:     jac->as_amg_beta_opts[0] = 10;
2091:     jac->as_amg_beta_opts[1] = 1;
2092:     jac->as_amg_beta_opts[2] = 6;
2093:     jac->as_amg_beta_opts[3] = 6;
2094:     jac->as_amg_beta_opts[4] = 4;
2095:     jac->as_amg_beta_theta   = 0.25;
2096:     PetscCallExternal(HYPRE_AMSSetPrintLevel, jac->hsolver, jac->as_print);
2097:     PetscCallExternal(HYPRE_AMSSetMaxIter, jac->hsolver, jac->as_max_iter);
2098:     PetscCallExternal(HYPRE_AMSSetCycleType, jac->hsolver, jac->ams_cycle_type);
2099:     PetscCallExternal(HYPRE_AMSSetTol, jac->hsolver, jac->as_tol);
2100:     PetscCallExternal(HYPRE_AMSSetSmoothingOptions, jac->hsolver, jac->as_relax_type, jac->as_relax_times, jac->as_relax_weight, jac->as_omega);
2101:     PetscCallExternal(HYPRE_AMSSetAlphaAMGOptions, jac->hsolver, jac->as_amg_alpha_opts[0], /* AMG coarsen type */
2102:                       jac->as_amg_alpha_opts[1],                                            /* AMG agg_levels */
2103:                       jac->as_amg_alpha_opts[2],                                            /* AMG relax_type */
2104:                       jac->as_amg_alpha_theta, jac->as_amg_alpha_opts[3],                   /* AMG interp_type */
2105:                       jac->as_amg_alpha_opts[4]);                                           /* AMG Pmax */
2106:     PetscCallExternal(HYPRE_AMSSetBetaAMGOptions, jac->hsolver, jac->as_amg_beta_opts[0],   /* AMG coarsen type */
2107:                       jac->as_amg_beta_opts[1],                                             /* AMG agg_levels */
2108:                       jac->as_amg_beta_opts[2],                                             /* AMG relax_type */
2109:                       jac->as_amg_beta_theta, jac->as_amg_beta_opts[3],                     /* AMG interp_type */
2110:                       jac->as_amg_beta_opts[4]);                                            /* AMG Pmax */
2111:     /* Zero conductivity */
2112:     jac->ams_beta_is_zero      = PETSC_FALSE;
2113:     jac->ams_beta_is_zero_part = PETSC_FALSE;
2114:     PetscFunctionReturn(PETSC_SUCCESS);
2115:   }
2116:   PetscCall(PetscStrcmp("ads", jac->hypre_type, &flag));
2117:   if (flag) {
2118:     PetscCallExternal(HYPRE_ADSCreate, &jac->hsolver);
2119:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ADS;
2120:     pc->ops->view           = PCView_HYPRE_ADS;
2121:     jac->destroy            = HYPRE_ADSDestroy;
2122:     jac->setup              = HYPRE_ADSSetup;
2123:     jac->solve              = HYPRE_ADSSolve;
2124:     jac->coords[0]          = NULL;
2125:     jac->coords[1]          = NULL;
2126:     jac->coords[2]          = NULL;
2127:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */
2128:     jac->as_print       = 0;
2129:     jac->as_max_iter    = 1;  /* used as a preconditioner */
2130:     jac->as_tol         = 0.; /* used as a preconditioner */
2131:     jac->ads_cycle_type = 13;
2132:     /* Smoothing options */
2133:     jac->as_relax_type   = 2;
2134:     jac->as_relax_times  = 1;
2135:     jac->as_relax_weight = 1.0;
2136:     jac->as_omega        = 1.0;
2137:     /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */
2138:     jac->ams_cycle_type       = 14;
2139:     jac->as_amg_alpha_opts[0] = 10;
2140:     jac->as_amg_alpha_opts[1] = 1;
2141:     jac->as_amg_alpha_opts[2] = 6;
2142:     jac->as_amg_alpha_opts[3] = 6;
2143:     jac->as_amg_alpha_opts[4] = 4;
2144:     jac->as_amg_alpha_theta   = 0.25;
2145:     /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
2146:     jac->as_amg_beta_opts[0] = 10;
2147:     jac->as_amg_beta_opts[1] = 1;
2148:     jac->as_amg_beta_opts[2] = 6;
2149:     jac->as_amg_beta_opts[3] = 6;
2150:     jac->as_amg_beta_opts[4] = 4;
2151:     jac->as_amg_beta_theta   = 0.25;
2152:     PetscCallExternal(HYPRE_ADSSetPrintLevel, jac->hsolver, jac->as_print);
2153:     PetscCallExternal(HYPRE_ADSSetMaxIter, jac->hsolver, jac->as_max_iter);
2154:     PetscCallExternal(HYPRE_ADSSetCycleType, jac->hsolver, jac->ams_cycle_type);
2155:     PetscCallExternal(HYPRE_ADSSetTol, jac->hsolver, jac->as_tol);
2156:     PetscCallExternal(HYPRE_ADSSetSmoothingOptions, jac->hsolver, jac->as_relax_type, jac->as_relax_times, jac->as_relax_weight, jac->as_omega);
2157:     PetscCallExternal(HYPRE_ADSSetAMSOptions, jac->hsolver, jac->ams_cycle_type,      /* AMG coarsen type */
2158:                       jac->as_amg_alpha_opts[0],                                      /* AMG coarsen type */
2159:                       jac->as_amg_alpha_opts[1],                                      /* AMG agg_levels */
2160:                       jac->as_amg_alpha_opts[2],                                      /* AMG relax_type */
2161:                       jac->as_amg_alpha_theta, jac->as_amg_alpha_opts[3],             /* AMG interp_type */
2162:                       jac->as_amg_alpha_opts[4]);                                     /* AMG Pmax */
2163:     PetscCallExternal(HYPRE_ADSSetAMGOptions, jac->hsolver, jac->as_amg_beta_opts[0], /* AMG coarsen type */
2164:                       jac->as_amg_beta_opts[1],                                       /* AMG agg_levels */
2165:                       jac->as_amg_beta_opts[2],                                       /* AMG relax_type */
2166:                       jac->as_amg_beta_theta, jac->as_amg_beta_opts[3],               /* AMG interp_type */
2167:                       jac->as_amg_beta_opts[4]);                                      /* AMG Pmax */
2168:     PetscFunctionReturn(PETSC_SUCCESS);
2169:   }
2170:   PetscCall(PetscFree(jac->hypre_type));

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

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

2186:   PetscFunctionBegin;
2187:   PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE preconditioner options");
2188:   PetscCall(PetscOptionsEList("-pc_hypre_type", "HYPRE preconditioner type", "PCHYPRESetType", type, PETSC_STATIC_ARRAY_LENGTH(type), "boomeramg", &indx, &flg));
2189:   if (flg) {
2190:     PetscCall(PCHYPRESetType_HYPRE(pc, type[indx]));
2191:   } else {
2192:     PetscCall(PCHYPRESetType_HYPRE(pc, "boomeramg"));
2193:   }
2194:   PetscTryTypeMethod(pc, setfromoptions, PetscOptionsObject);
2195:   PetscOptionsHeadEnd();
2196:   PetscFunctionReturn(PETSC_SUCCESS);
2197: }

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

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

2206:   Options Database Key:
2207: . pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads

2209:   Level: intermediate

2211: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHYPRE`
2212: @*/
2213: PetscErrorCode PCHYPRESetType(PC pc, const char name[])
2214: {
2215:   PetscFunctionBegin;
2217:   PetscAssertPointer(name, 2);
2218:   PetscTryMethod(pc, "PCHYPRESetType_C", (PC, const char[]), (pc, name));
2219:   PetscFunctionReturn(PETSC_SUCCESS);
2220: }

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

2225:   Input Parameter:
2226: . pc - the preconditioner context

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

2231:   Level: intermediate

2233: .seealso: [](ch_ksp), `PCCreate()`, `PCHYPRESetType()`, `PCType`, `PC`, `PCHYPRE`
2234: @*/
2235: PetscErrorCode PCHYPREGetType(PC pc, const char *name[])
2236: {
2237:   PetscFunctionBegin;
2239:   PetscAssertPointer(name, 2);
2240:   PetscTryMethod(pc, "PCHYPREGetType_C", (PC, const char *[]), (pc, name));
2241:   PetscFunctionReturn(PETSC_SUCCESS);
2242: }

2244: /*@C
2245:   PCMGGalerkinSetMatProductAlgorithm - Set type of SpGEMM for hypre to use on GPUs

2247:   Logically Collective

2249:   Input Parameters:
2250: + pc   - the hypre context
2251: - name - one of 'cusparse', 'hypre'

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

2256:   Level: intermediate

2258:   Developer Notes:
2259:   How the name starts with `PCMG`, should it not be `PCHYPREBoomerAMG`?

2261: .seealso: [](ch_ksp), `PCHYPRE`, `PCMGGalerkinGetMatProductAlgorithm()`
2262: @*/
2263: PetscErrorCode PCMGGalerkinSetMatProductAlgorithm(PC pc, const char name[])
2264: {
2265:   PetscFunctionBegin;
2267:   PetscTryMethod(pc, "PCMGGalerkinSetMatProductAlgorithm_C", (PC, const char[]), (pc, name));
2268:   PetscFunctionReturn(PETSC_SUCCESS);
2269: }

2271: /*@C
2272:   PCMGGalerkinGetMatProductAlgorithm - Get type of SpGEMM for hypre to use on GPUs

2274:   Not Collective

2276:   Input Parameter:
2277: . pc - the multigrid context

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

2282:   Level: intermediate

2284: .seealso: [](ch_ksp), `PCHYPRE`, `PCMGGalerkinSetMatProductAlgorithm()`
2285: @*/
2286: PetscErrorCode PCMGGalerkinGetMatProductAlgorithm(PC pc, const char *name[])
2287: {
2288:   PetscFunctionBegin;
2290:   PetscTryMethod(pc, "PCMGGalerkinGetMatProductAlgorithm_C", (PC, const char *[]), (pc, name));
2291:   PetscFunctionReturn(PETSC_SUCCESS);
2292: }

2294: /*MC
2295:      PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre as PETSc `PC`

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

2303:    Level: intermediate

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

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

2318:            Note that the option `-pc_hypre_boomeramg_relax_type_all` defaults to symmetric relaxation
2319:            (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry.
2320:            Otherwise, you may want to use `-pc_hypre_boomeramg_relax_type_all SOR/Jacobi`.

2322:           `MatSetNearNullSpace()` - if you provide a near null space to your matrix it is ignored by hypre UNLESS you also use
2323:           the following two options: `-pc_hypre_boomeramg_nodal_coarsen <n> -pc_hypre_boomeramg_vec_interp_variant <v>`

2325:           See `PCPFMG`, `PCSMG`, and `PCSYSPFMG` for access to hypre's other (nonalgebraic) multigrid solvers

2327:           For `PCHYPRE` type of `ams` or `ads` auxiliary data must be provided to the preconditioner with `PCHYPRESetDiscreteGradient()`,
2328:           `PCHYPRESetDiscreteCurl()`, `PCHYPRESetInterpolations()`, `PCHYPRESetAlphaPoissonMatrix()`, `PCHYPRESetBetaPoissonMatrix()`, `PCHYPRESetEdgeConstantVectors()`,
2329:           `PCHYPREAMSSetInteriorNodes()`

2331:   Sometimes people want to try algebraic multigrid as a "standalone" solver, that is not accelerating it with a Krylov method. Though we generally do not recommend this
2332:   since it is usually slower, one should use a `KSPType` of `KSPRICHARDSON`
2333:   (or equivalently `-ksp_type richardson`) to achieve this. Using `KSPPREONLY` will not work since it only applies a single cycle of multigrid.

2335:    PETSc provides its own geometric and algebraic multigrid solvers `PCMG` and `PCGAMG`, also see `PCHMG` which is useful for certain multicomponent problems

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

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

2344: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHYPRESetType()`, `PCPFMG`, `PCGAMG`, `PCSYSPFMG`, `PCSMG`, `PCHYPRESetDiscreteGradient()`,
2345:           `PCHYPRESetDiscreteCurl()`, `PCHYPRESetInterpolations()`, `PCHYPRESetAlphaPoissonMatrix()`, `PCHYPRESetBetaPoissonMatrix()`, `PCHYPRESetEdgeConstantVectors()`,
2346:           PCHYPREAMSSetInteriorNodes()
2347: M*/

2349: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
2350: {
2351:   PC_HYPRE *jac;

2353:   PetscFunctionBegin;
2354:   PetscCall(PetscNew(&jac));

2356:   pc->data                = jac;
2357:   pc->ops->reset          = PCReset_HYPRE;
2358:   pc->ops->destroy        = PCDestroy_HYPRE;
2359:   pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
2360:   pc->ops->setup          = PCSetUp_HYPRE;
2361:   pc->ops->apply          = PCApply_HYPRE;
2362:   jac->comm_hypre         = MPI_COMM_NULL;
2363:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetType_C", PCHYPRESetType_HYPRE));
2364:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetType_C", PCHYPREGetType_HYPRE));
2365:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_HYPRE));
2366:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteGradient_C", PCHYPRESetDiscreteGradient_HYPRE));
2367:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteCurl_C", PCHYPRESetDiscreteCurl_HYPRE));
2368:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetInterpolations_C", PCHYPRESetInterpolations_HYPRE));
2369:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetEdgeConstantVectors_C", PCHYPRESetEdgeConstantVectors_HYPRE));
2370:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREAMSSetInteriorNodes_C", PCHYPREAMSSetInteriorNodes_HYPRE));
2371:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetPoissonMatrix_C", PCHYPRESetPoissonMatrix_HYPRE));
2372:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinSetMatProductAlgorithm_C", PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG));
2373:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinGetMatProductAlgorithm_C", PCMGGalerkinGetMatProductAlgorithm_HYPRE_BoomerAMG));
2374: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2375:   #if defined(HYPRE_USING_HIP)
2376:   PetscCall(PetscDeviceInitialize(PETSC_DEVICE_HIP));
2377:   #endif
2378:   #if defined(HYPRE_USING_CUDA)
2379:   PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA));
2380:   #endif
2381: #endif
2382:   PetscHYPREInitialize();
2383:   PetscFunctionReturn(PETSC_SUCCESS);
2384: }

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

2390:   /* keep copy of PFMG options used so may view them */
2391:   PetscInt  its;
2392:   PetscReal tol;
2393:   PetscInt  relax_type;
2394:   PetscInt  rap_type;
2395:   PetscInt  num_pre_relax, num_post_relax;
2396:   PetscInt  max_levels;
2397:   PetscInt  skip_relax;
2398:   PetscBool print_statistics;
2399: } PC_PFMG;

2401: static PetscErrorCode PCDestroy_PFMG(PC pc)
2402: {
2403:   PC_PFMG *ex = (PC_PFMG *)pc->data;

2405:   PetscFunctionBegin;
2406:   if (ex->hsolver) PetscCallExternal(HYPRE_StructPFMGDestroy, ex->hsolver);
2407:   PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2408:   PetscCall(PetscFree(pc->data));
2409:   PetscFunctionReturn(PETSC_SUCCESS);
2410: }

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

2415: static PetscErrorCode PCView_PFMG(PC pc, PetscViewer viewer)
2416: {
2417:   PetscBool iascii;
2418:   PC_PFMG  *ex = (PC_PFMG *)pc->data;

2420:   PetscFunctionBegin;
2421:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
2422:   if (iascii) {
2423:     PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE PFMG preconditioning\n"));
2424:     PetscCall(PetscViewerASCIIPrintf(viewer, "    max iterations %" PetscInt_FMT "\n", ex->its));
2425:     PetscCall(PetscViewerASCIIPrintf(viewer, "    tolerance %g\n", ex->tol));
2426:     PetscCall(PetscViewerASCIIPrintf(viewer, "    relax type %s\n", PFMGRelaxType[ex->relax_type]));
2427:     PetscCall(PetscViewerASCIIPrintf(viewer, "    RAP type %s\n", PFMGRAPType[ex->rap_type]));
2428:     PetscCall(PetscViewerASCIIPrintf(viewer, "    number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax));
2429:     PetscCall(PetscViewerASCIIPrintf(viewer, "    max levels %" PetscInt_FMT "\n", ex->max_levels));
2430:     PetscCall(PetscViewerASCIIPrintf(viewer, "    skip relax %" PetscInt_FMT "\n", ex->skip_relax));
2431:   }
2432:   PetscFunctionReturn(PETSC_SUCCESS);
2433: }

2435: static PetscErrorCode PCSetFromOptions_PFMG(PC pc, PetscOptionItems *PetscOptionsObject)
2436: {
2437:   PC_PFMG *ex = (PC_PFMG *)pc->data;

2439:   PetscFunctionBegin;
2440:   PetscOptionsHeadBegin(PetscOptionsObject, "PFMG options");
2441:   PetscCall(PetscOptionsBool("-pc_pfmg_print_statistics", "Print statistics", "HYPRE_StructPFMGSetPrintLevel", ex->print_statistics, &ex->print_statistics, NULL));
2442:   PetscCall(PetscOptionsInt("-pc_pfmg_its", "Number of iterations of PFMG to use as preconditioner", "HYPRE_StructPFMGSetMaxIter", ex->its, &ex->its, NULL));
2443:   PetscCallExternal(HYPRE_StructPFMGSetMaxIter, ex->hsolver, ex->its);
2444:   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));
2445:   PetscCallExternal(HYPRE_StructPFMGSetNumPreRelax, ex->hsolver, ex->num_pre_relax);
2446:   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));
2447:   PetscCallExternal(HYPRE_StructPFMGSetNumPostRelax, ex->hsolver, ex->num_post_relax);

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

2452:   PetscCall(PetscOptionsReal("-pc_pfmg_tol", "Tolerance of PFMG", "HYPRE_StructPFMGSetTol", ex->tol, &ex->tol, NULL));
2453:   PetscCallExternal(HYPRE_StructPFMGSetTol, ex->hsolver, ex->tol);
2454:   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));
2455:   PetscCallExternal(HYPRE_StructPFMGSetRelaxType, ex->hsolver, ex->relax_type);
2456:   PetscCall(PetscOptionsEList("-pc_pfmg_rap_type", "RAP type", "HYPRE_StructPFMGSetRAPType", PFMGRAPType, PETSC_STATIC_ARRAY_LENGTH(PFMGRAPType), PFMGRAPType[ex->rap_type], &ex->rap_type, NULL));
2457:   PetscCallExternal(HYPRE_StructPFMGSetRAPType, ex->hsolver, ex->rap_type);
2458:   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));
2459:   PetscCallExternal(HYPRE_StructPFMGSetSkipRelax, ex->hsolver, ex->skip_relax);
2460:   PetscOptionsHeadEnd();
2461:   PetscFunctionReturn(PETSC_SUCCESS);
2462: }

2464: static PetscErrorCode PCApply_PFMG(PC pc, Vec x, Vec y)
2465: {
2466:   PC_PFMG           *ex = (PC_PFMG *)pc->data;
2467:   PetscScalar       *yy;
2468:   const PetscScalar *xx;
2469:   PetscInt           ilower[3], iupper[3];
2470:   HYPRE_Int          hlower[3], hupper[3];
2471:   Mat_HYPREStruct   *mx = (Mat_HYPREStruct *)pc->pmat->data;

2473:   PetscFunctionBegin;
2474:   PetscCall(PetscCitationsRegister(hypreCitation, &cite));
2475:   PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
2476:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2477:   iupper[0] += ilower[0] - 1;
2478:   iupper[1] += ilower[1] - 1;
2479:   iupper[2] += ilower[2] - 1;
2480:   hlower[0] = (HYPRE_Int)ilower[0];
2481:   hlower[1] = (HYPRE_Int)ilower[1];
2482:   hlower[2] = (HYPRE_Int)ilower[2];
2483:   hupper[0] = (HYPRE_Int)iupper[0];
2484:   hupper[1] = (HYPRE_Int)iupper[1];
2485:   hupper[2] = (HYPRE_Int)iupper[2];

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

2495:   /* copy solution values back to PETSc */
2496:   PetscCall(VecGetArray(y, &yy));
2497:   PetscCallExternal(HYPRE_StructVectorGetBoxValues, mx->hx, hlower, hupper, (HYPRE_Complex *)yy);
2498:   PetscCall(VecRestoreArray(y, &yy));
2499:   PetscFunctionReturn(PETSC_SUCCESS);
2500: }

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

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

2512:   PetscCall(PCApply_PFMG(pc, b, y));
2513:   PetscCallExternal(HYPRE_StructPFMGGetNumIterations, jac->hsolver, &oits);
2514:   *outits = oits;
2515:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2516:   else *reason = PCRICHARDSON_CONVERGED_RTOL;
2517:   PetscCallExternal(HYPRE_StructPFMGSetTol, jac->hsolver, jac->tol);
2518:   PetscCallExternal(HYPRE_StructPFMGSetMaxIter, jac->hsolver, jac->its);
2519:   PetscFunctionReturn(PETSC_SUCCESS);
2520: }

2522: static PetscErrorCode PCSetUp_PFMG(PC pc)
2523: {
2524:   PC_PFMG         *ex = (PC_PFMG *)pc->data;
2525:   Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data;
2526:   PetscBool        flg;

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

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

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

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

2548:   PetscCallExternal(HYPRE_StructPFMGSetup, ex->hsolver, mx->hmat, mx->hb, mx->hx);
2549:   PetscCallExternal(HYPRE_StructPFMGSetZeroGuess, ex->hsolver);
2550:   PetscFunctionReturn(PETSC_SUCCESS);
2551: }

2553: /*MC
2554:      PCPFMG - the hypre PFMG multigrid solver

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

2566:    Level: advanced

2568:    Notes:
2569:    This is for CELL-centered descretizations

2571:    See `PCSYSPFMG` for a version suitable for systems of PDEs, and `PCSMG`

2573:    See `PCHYPRE` for hypre's BoomerAMG algebraic multigrid solver

2575:    This must be used with the `MATHYPRESTRUCT` matrix type.

2577:    This provides only some of the functionality of PFMG, it supports only one block per process defined by a PETSc `DMDA`.

2579: .seealso: [](ch_ksp), `PCMG`, `MATHYPRESTRUCT`, `PCHYPRE`, `PCGAMG`, `PCSYSPFMG`, `PCSMG`
2580: M*/

2582: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
2583: {
2584:   PC_PFMG *ex;

2586:   PetscFunctionBegin;
2587:   PetscCall(PetscNew(&ex));
2588:   pc->data = ex;

2590:   ex->its              = 1;
2591:   ex->tol              = 1.e-8;
2592:   ex->relax_type       = 1;
2593:   ex->rap_type         = 0;
2594:   ex->num_pre_relax    = 1;
2595:   ex->num_post_relax   = 1;
2596:   ex->max_levels       = 0;
2597:   ex->skip_relax       = 0;
2598:   ex->print_statistics = PETSC_FALSE;

2600:   pc->ops->setfromoptions  = PCSetFromOptions_PFMG;
2601:   pc->ops->view            = PCView_PFMG;
2602:   pc->ops->destroy         = PCDestroy_PFMG;
2603:   pc->ops->apply           = PCApply_PFMG;
2604:   pc->ops->applyrichardson = PCApplyRichardson_PFMG;
2605:   pc->ops->setup           = PCSetUp_PFMG;

2607:   PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2608:   PetscHYPREInitialize();
2609:   PetscCallExternal(HYPRE_StructPFMGCreate, ex->hcomm, &ex->hsolver);
2610:   PetscFunctionReturn(PETSC_SUCCESS);
2611: }

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

2618:   /* keep copy of SYSPFMG options used so may view them */
2619:   PetscInt  its;
2620:   PetscReal tol;
2621:   PetscInt  relax_type;
2622:   PetscInt  num_pre_relax, num_post_relax;
2623: } PC_SysPFMG;

2625: static PetscErrorCode PCDestroy_SysPFMG(PC pc)
2626: {
2627:   PC_SysPFMG *ex = (PC_SysPFMG *)pc->data;

2629:   PetscFunctionBegin;
2630:   if (ex->ss_solver) PetscCallExternal(HYPRE_SStructSysPFMGDestroy, ex->ss_solver);
2631:   PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2632:   PetscCall(PetscFree(pc->data));
2633:   PetscFunctionReturn(PETSC_SUCCESS);
2634: }

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

2638: static PetscErrorCode PCView_SysPFMG(PC pc, PetscViewer viewer)
2639: {
2640:   PetscBool   iascii;
2641:   PC_SysPFMG *ex = (PC_SysPFMG *)pc->data;

2643:   PetscFunctionBegin;
2644:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
2645:   if (iascii) {
2646:     PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE SysPFMG preconditioning\n"));
2647:     PetscCall(PetscViewerASCIIPrintf(viewer, "  max iterations %" PetscInt_FMT "\n", ex->its));
2648:     PetscCall(PetscViewerASCIIPrintf(viewer, "  tolerance %g\n", ex->tol));
2649:     PetscCall(PetscViewerASCIIPrintf(viewer, "  relax type %s\n", PFMGRelaxType[ex->relax_type]));
2650:     PetscCall(PetscViewerASCIIPrintf(viewer, "  number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax));
2651:   }
2652:   PetscFunctionReturn(PETSC_SUCCESS);
2653: }

2655: static PetscErrorCode PCSetFromOptions_SysPFMG(PC pc, PetscOptionItems *PetscOptionsObject)
2656: {
2657:   PC_SysPFMG *ex  = (PC_SysPFMG *)pc->data;
2658:   PetscBool   flg = PETSC_FALSE;

2660:   PetscFunctionBegin;
2661:   PetscOptionsHeadBegin(PetscOptionsObject, "SysPFMG options");
2662:   PetscCall(PetscOptionsBool("-pc_syspfmg_print_statistics", "Print statistics", "HYPRE_SStructSysPFMGSetPrintLevel", flg, &flg, NULL));
2663:   if (flg) PetscCallExternal(HYPRE_SStructSysPFMGSetPrintLevel, ex->ss_solver, 3);
2664:   PetscCall(PetscOptionsInt("-pc_syspfmg_its", "Number of iterations of SysPFMG to use as preconditioner", "HYPRE_SStructSysPFMGSetMaxIter", ex->its, &ex->its, NULL));
2665:   PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter, ex->ss_solver, ex->its);
2666:   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));
2667:   PetscCallExternal(HYPRE_SStructSysPFMGSetNumPreRelax, ex->ss_solver, ex->num_pre_relax);
2668:   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));
2669:   PetscCallExternal(HYPRE_SStructSysPFMGSetNumPostRelax, ex->ss_solver, ex->num_post_relax);

2671:   PetscCall(PetscOptionsReal("-pc_syspfmg_tol", "Tolerance of SysPFMG", "HYPRE_SStructSysPFMGSetTol", ex->tol, &ex->tol, NULL));
2672:   PetscCallExternal(HYPRE_SStructSysPFMGSetTol, ex->ss_solver, ex->tol);
2673:   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));
2674:   PetscCallExternal(HYPRE_SStructSysPFMGSetRelaxType, ex->ss_solver, ex->relax_type);
2675:   PetscOptionsHeadEnd();
2676:   PetscFunctionReturn(PETSC_SUCCESS);
2677: }

2679: static PetscErrorCode PCApply_SysPFMG(PC pc, Vec x, Vec y)
2680: {
2681:   PC_SysPFMG        *ex = (PC_SysPFMG *)pc->data;
2682:   PetscScalar       *yy;
2683:   const PetscScalar *xx;
2684:   PetscInt           ilower[3], iupper[3];
2685:   HYPRE_Int          hlower[3], hupper[3];
2686:   Mat_HYPRESStruct  *mx       = (Mat_HYPRESStruct *)pc->pmat->data;
2687:   PetscInt           ordering = mx->dofs_order;
2688:   PetscInt           nvars    = mx->nvars;
2689:   PetscInt           part     = 0;
2690:   PetscInt           size;
2691:   PetscInt           i;

2693:   PetscFunctionBegin;
2694:   PetscCall(PetscCitationsRegister(hypreCitation, &cite));
2695:   PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
2696:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2697:   iupper[0] += ilower[0] - 1;
2698:   iupper[1] += ilower[1] - 1;
2699:   iupper[2] += ilower[2] - 1;
2700:   hlower[0] = (HYPRE_Int)ilower[0];
2701:   hlower[1] = (HYPRE_Int)ilower[1];
2702:   hlower[2] = (HYPRE_Int)ilower[2];
2703:   hupper[0] = (HYPRE_Int)iupper[0];
2704:   hupper[1] = (HYPRE_Int)iupper[1];
2705:   hupper[2] = (HYPRE_Int)iupper[2];

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

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

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

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

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

2742:     /* copy solution values back to PETSc */
2743:     PetscCall(VecGetArray(y, &yy));
2744:     for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues, mx->ss_x, part, hlower, hupper, i, (HYPRE_Complex *)(z + (size * i)));
2745:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
2746:     for (i = 0; i < size; i++) {
2747:       k = i * nvars;
2748:       for (j = 0; j < nvars; j++) yy[k + j] = z[j * size + i];
2749:     }
2750:     PetscCall(VecRestoreArray(y, &yy));
2751:     PetscCall(PetscFree(z));
2752:   }
2753:   PetscFunctionReturn(PETSC_SUCCESS);
2754: }

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

2761:   PetscFunctionBegin;
2762:   PetscCall(PetscCitationsRegister(hypreCitation, &cite));
2763:   PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter, jac->ss_solver, its * jac->its);
2764:   PetscCallExternal(HYPRE_SStructSysPFMGSetTol, jac->ss_solver, rtol);
2765:   PetscCall(PCApply_SysPFMG(pc, b, y));
2766:   PetscCallExternal(HYPRE_SStructSysPFMGGetNumIterations, jac->ss_solver, &oits);
2767:   *outits = oits;
2768:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2769:   else *reason = PCRICHARDSON_CONVERGED_RTOL;
2770:   PetscCallExternal(HYPRE_SStructSysPFMGSetTol, jac->ss_solver, jac->tol);
2771:   PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter, jac->ss_solver, jac->its);
2772:   PetscFunctionReturn(PETSC_SUCCESS);
2773: }

2775: static PetscErrorCode PCSetUp_SysPFMG(PC pc)
2776: {
2777:   PC_SysPFMG       *ex = (PC_SysPFMG *)pc->data;
2778:   Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)pc->pmat->data;
2779:   PetscBool         flg;

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

2785:   /* create the hypre sstruct solver object and set its information */
2786:   if (ex->ss_solver) PetscCallExternal(HYPRE_SStructSysPFMGDestroy, ex->ss_solver);
2787:   PetscCallExternal(HYPRE_SStructSysPFMGCreate, ex->hcomm, &ex->ss_solver);
2788:   PetscCallExternal(HYPRE_SStructSysPFMGSetZeroGuess, ex->ss_solver);
2789:   PetscCallExternal(HYPRE_SStructSysPFMGSetup, ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x);
2790:   PetscFunctionReturn(PETSC_SUCCESS);
2791: }

2793: /*MC
2794:      PCSYSPFMG - the hypre SysPFMG multigrid solver

2796:    Level: advanced

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

2805:    Notes:
2806:    See `PCPFMG` for hypre's PFMG that works for a scalar PDE and `PCSMG`

2808:    See `PCHYPRE` for hypre's BoomerAMG algebraic multigrid solver

2810:    This is for CELL-centered descretizations

2812:    This must be used with the `MATHYPRESSTRUCT` matrix type.

2814:    This does not give access to all the functionality of hypres SysPFMG, it supports only one part, and one block per process defined by a PETSc `DMDA`.

2816: .seealso: [](ch_ksp), `PCMG`, `MATHYPRESSTRUCT`, `PCPFMG`, `PCHYPRE`, `PCGAMG`, `PCSMG`
2817: M*/

2819: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
2820: {
2821:   PC_SysPFMG *ex;

2823:   PetscFunctionBegin;
2824:   PetscCall(PetscNew(&ex));
2825:   pc->data = ex;

2827:   ex->its            = 1;
2828:   ex->tol            = 1.e-8;
2829:   ex->relax_type     = 1;
2830:   ex->num_pre_relax  = 1;
2831:   ex->num_post_relax = 1;

2833:   pc->ops->setfromoptions  = PCSetFromOptions_SysPFMG;
2834:   pc->ops->view            = PCView_SysPFMG;
2835:   pc->ops->destroy         = PCDestroy_SysPFMG;
2836:   pc->ops->apply           = PCApply_SysPFMG;
2837:   pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
2838:   pc->ops->setup           = PCSetUp_SysPFMG;

2840:   PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2841:   PetscHYPREInitialize();
2842:   PetscCallExternal(HYPRE_SStructSysPFMGCreate, ex->hcomm, &ex->ss_solver);
2843:   PetscFunctionReturn(PETSC_SUCCESS);
2844: }

2846: /* PC SMG */
2847: typedef struct {
2848:   MPI_Comm           hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */
2849:   HYPRE_StructSolver hsolver;
2850:   PetscInt           its; /* keep copy of SMG options used so may view them */
2851:   PetscReal          tol;
2852:   PetscBool          print_statistics;
2853:   PetscInt           num_pre_relax, num_post_relax;
2854: } PC_SMG;

2856: static PetscErrorCode PCDestroy_SMG(PC pc)
2857: {
2858:   PC_SMG *ex = (PC_SMG *)pc->data;

2860:   PetscFunctionBegin;
2861:   if (ex->hsolver) PetscCallExternal(HYPRE_StructSMGDestroy, ex->hsolver);
2862:   PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2863:   PetscCall(PetscFree(pc->data));
2864:   PetscFunctionReturn(PETSC_SUCCESS);
2865: }

2867: static PetscErrorCode PCView_SMG(PC pc, PetscViewer viewer)
2868: {
2869:   PetscBool iascii;
2870:   PC_SMG   *ex = (PC_SMG *)pc->data;

2872:   PetscFunctionBegin;
2873:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
2874:   if (iascii) {
2875:     PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE SMG preconditioning\n"));
2876:     PetscCall(PetscViewerASCIIPrintf(viewer, "    max iterations %" PetscInt_FMT "\n", ex->its));
2877:     PetscCall(PetscViewerASCIIPrintf(viewer, "    tolerance %g\n", ex->tol));
2878:     PetscCall(PetscViewerASCIIPrintf(viewer, "    number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax));
2879:   }
2880:   PetscFunctionReturn(PETSC_SUCCESS);
2881: }

2883: static PetscErrorCode PCSetFromOptions_SMG(PC pc, PetscOptionItems *PetscOptionsObject)
2884: {
2885:   PC_SMG *ex = (PC_SMG *)pc->data;

2887:   PetscFunctionBegin;
2888:   PetscOptionsHeadBegin(PetscOptionsObject, "SMG options");

2890:   PetscCall(PetscOptionsInt("-pc_smg_its", "Number of iterations of SMG to use as preconditioner", "HYPRE_StructSMGSetMaxIter", ex->its, &ex->its, NULL));
2891:   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));
2892:   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));
2893:   PetscCall(PetscOptionsReal("-pc_smg_tol", "Tolerance of SMG", "HYPRE_StructSMGSetTol", ex->tol, &ex->tol, NULL));

2895:   PetscOptionsHeadEnd();
2896:   PetscFunctionReturn(PETSC_SUCCESS);
2897: }

2899: static PetscErrorCode PCApply_SMG(PC pc, Vec x, Vec y)
2900: {
2901:   PC_SMG            *ex = (PC_SMG *)pc->data;
2902:   PetscScalar       *yy;
2903:   const PetscScalar *xx;
2904:   PetscInt           ilower[3], iupper[3];
2905:   HYPRE_Int          hlower[3], hupper[3];
2906:   Mat_HYPREStruct   *mx = (Mat_HYPREStruct *)pc->pmat->data;

2908:   PetscFunctionBegin;
2909:   PetscCall(PetscCitationsRegister(hypreCitation, &cite));
2910:   PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
2911:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2912:   iupper[0] += ilower[0] - 1;
2913:   iupper[1] += ilower[1] - 1;
2914:   iupper[2] += ilower[2] - 1;
2915:   hlower[0] = (HYPRE_Int)ilower[0];
2916:   hlower[1] = (HYPRE_Int)ilower[1];
2917:   hlower[2] = (HYPRE_Int)ilower[2];
2918:   hupper[0] = (HYPRE_Int)iupper[0];
2919:   hupper[1] = (HYPRE_Int)iupper[1];
2920:   hupper[2] = (HYPRE_Int)iupper[2];

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

2930:   /* copy solution values back to PETSc */
2931:   PetscCall(VecGetArray(y, &yy));
2932:   PetscCallExternal(HYPRE_StructVectorGetBoxValues, mx->hx, hlower, hupper, (HYPRE_Complex *)yy);
2933:   PetscCall(VecRestoreArray(y, &yy));
2934:   PetscFunctionReturn(PETSC_SUCCESS);
2935: }

2937: 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)
2938: {
2939:   PC_SMG   *jac = (PC_SMG *)pc->data;
2940:   HYPRE_Int oits;

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

2947:   PetscCall(PCApply_SMG(pc, b, y));
2948:   PetscCallExternal(HYPRE_StructSMGGetNumIterations, jac->hsolver, &oits);
2949:   *outits = oits;
2950:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2951:   else *reason = PCRICHARDSON_CONVERGED_RTOL;
2952:   PetscCallExternal(HYPRE_StructSMGSetTol, jac->hsolver, jac->tol);
2953:   PetscCallExternal(HYPRE_StructSMGSetMaxIter, jac->hsolver, jac->its);
2954:   PetscFunctionReturn(PETSC_SUCCESS);
2955: }

2957: static PetscErrorCode PCSetUp_SMG(PC pc)
2958: {
2959:   PetscInt         i, dim;
2960:   PC_SMG          *ex = (PC_SMG *)pc->data;
2961:   Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data;
2962:   PetscBool        flg;
2963:   DMBoundaryType   p[3];
2964:   PetscInt         M[3];

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

2970:   PetscCall(DMDAGetInfo(mx->da, &dim, &M[0], &M[1], &M[2], 0, 0, 0, 0, 0, &p[0], &p[1], &p[2], 0));
2971:   // Check if power of 2 in periodic directions
2972:   for (i = 0; i < dim; i++) {
2973:     if (((M[i] & (M[i] - 1)) != 0) && (p[i] == DM_BOUNDARY_PERIODIC)) {
2974:       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]);
2975:     }
2976:   }

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

2987:   PetscCallExternal(HYPRE_StructSMGSetup, ex->hsolver, mx->hmat, mx->hb, mx->hx);
2988:   PetscCallExternal(HYPRE_StructSMGSetZeroGuess, ex->hsolver);
2989:   PetscFunctionReturn(PETSC_SUCCESS);
2990: }

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

2995:    Level: advanced

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

3003:    Notes:
3004:    This is for CELL-centered descretizations

3006:    This must be used with the `MATHYPRESTRUCT` `MatType`.

3008:    This does not provide all the functionality of  hypre's SMG solver, it supports only one block per process defined by a PETSc `DMDA`.

3010:    See `PCSYSPFMG`, `PCSMG`, `PCPFMG`, and `PCHYPRE` for access to hypre's other preconditioners

3012: .seealso:  `PCMG`, `MATHYPRESTRUCT`, `PCPFMG`, `PCSYSPFMG`, `PCHYPRE`, `PCGAMG`
3013: M*/

3015: PETSC_EXTERN PetscErrorCode PCCreate_SMG(PC pc)
3016: {
3017:   PC_SMG *ex;

3019:   PetscFunctionBegin;
3020:   PetscCall(PetscNew(&ex));
3021:   pc->data = ex;

3023:   ex->its            = 1;
3024:   ex->tol            = 1.e-8;
3025:   ex->num_pre_relax  = 1;
3026:   ex->num_post_relax = 1;

3028:   pc->ops->setfromoptions  = PCSetFromOptions_SMG;
3029:   pc->ops->view            = PCView_SMG;
3030:   pc->ops->destroy         = PCDestroy_SMG;
3031:   pc->ops->apply           = PCApply_SMG;
3032:   pc->ops->applyrichardson = PCApplyRichardson_SMG;
3033:   pc->ops->setup           = PCSetUp_SMG;

3035:   PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
3036:   PetscHYPREInitialize();
3037:   PetscCallExternal(HYPRE_StructSMGCreate, ex->hcomm, &ex->hsolver);
3038:   PetscFunctionReturn(PETSC_SUCCESS);
3039: }