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: PetscErrorCode PCHYPREGetSolver(PC pc, HYPRE_Solver *hsolver)
148: {
149:   PC_HYPRE *jac = (PC_HYPRE *)pc->data;

151:   PetscFunctionBegin;
152:   *hsolver = jac->hsolver;
153:   PetscFunctionReturn(PETSC_SUCCESS);
154: }

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

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

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

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

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

220:   PetscFunctionBegin;
221:   for (i = 0; i < jac->n_hmnull; i++) PetscCall(VecHYPRE_IJVectorDestroy(&jac->hmnull[i]));
222:   PetscCall(PetscFree(jac->hmnull));
223:   PetscCall(PetscFree(jac->phmnull));
224:   PetscCall(VecDestroy(&jac->hmnull_constant));
225:   jac->n_hmnull = 0;
226:   PetscFunctionReturn(PETSC_SUCCESS);
227: }

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

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

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

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

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

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

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

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

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

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

460:   PetscFunctionBegin;
461:   PetscCall(MatDestroy(&jac->hpmat));
462:   PetscCall(MatDestroy(&jac->G));
463:   PetscCall(MatDestroy(&jac->C));
464:   PetscCall(MatDestroy(&jac->alpha_Poisson));
465:   PetscCall(MatDestroy(&jac->beta_Poisson));
466:   PetscCall(MatDestroy(&jac->RT_PiFull));
467:   PetscCall(MatDestroy(&jac->RT_Pi[0]));
468:   PetscCall(MatDestroy(&jac->RT_Pi[1]));
469:   PetscCall(MatDestroy(&jac->RT_Pi[2]));
470:   PetscCall(MatDestroy(&jac->ND_PiFull));
471:   PetscCall(MatDestroy(&jac->ND_Pi[0]));
472:   PetscCall(MatDestroy(&jac->ND_Pi[1]));
473:   PetscCall(MatDestroy(&jac->ND_Pi[2]));
474:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[0]));
475:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[1]));
476:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[2]));
477:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[0]));
478:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[1]));
479:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[2]));
480:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->interior));
481:   PetscCall(PCHYPREResetNearNullSpace_Private(pc));
482:   jac->ams_beta_is_zero      = PETSC_FALSE;
483:   jac->ams_beta_is_zero_part = PETSC_FALSE;
484:   jac->dim                   = 0;
485:   PetscFunctionReturn(PETSC_SUCCESS);
486: }

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

492:   PetscFunctionBegin;
493:   PetscCall(PCReset_HYPRE(pc));
494:   if (jac->destroy) PetscCallExternal(jac->destroy, jac->hsolver);
495:   PetscCall(PetscFree(jac->hypre_type));
496: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
497:   PetscCall(PetscFree(jac->spgemm_type));
498: #endif
499:   if (jac->comm_hypre != MPI_COMM_NULL) PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
500:   PetscCall(PetscFree(pc->data));

502:   PetscCall(PetscObjectChangeTypeName((PetscObject)pc, 0));
503:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetType_C", NULL));
504:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetType_C", NULL));
505:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteGradient_C", NULL));
506:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteCurl_C", NULL));
507:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetInterpolations_C", NULL));
508:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetConstantEdgeVectors_C", NULL));
509:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetPoissonMatrix_C", NULL));
510:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetEdgeConstantVectors_C", NULL));
511:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREAMSSetInteriorNodes_C", NULL));
512:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", NULL));
513:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", NULL));
514:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinSetMatProductAlgorithm_C", NULL));
515:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinGetMatProductAlgorithm_C", NULL));
516:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
517:   PetscFunctionReturn(PETSC_SUCCESS);
518: }

520: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PC pc, PetscOptionItems *PetscOptionsObject)
521: {
522:   PC_HYPRE *jac = (PC_HYPRE *)pc->data;
523:   PetscBool flag;

525:   PetscFunctionBegin;
526:   PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE Pilut Options");
527:   PetscCall(PetscOptionsInt("-pc_hypre_pilut_maxiter", "Number of iterations", "None", jac->maxiter, &jac->maxiter, &flag));
528:   if (flag) PetscCallExternal(HYPRE_ParCSRPilutSetMaxIter, jac->hsolver, jac->maxiter);
529:   PetscCall(PetscOptionsReal("-pc_hypre_pilut_tol", "Drop tolerance", "None", jac->tol, &jac->tol, &flag));
530:   if (flag) PetscCallExternal(HYPRE_ParCSRPilutSetDropTolerance, jac->hsolver, jac->tol);
531:   PetscCall(PetscOptionsInt("-pc_hypre_pilut_factorrowsize", "FactorRowSize", "None", jac->factorrowsize, &jac->factorrowsize, &flag));
532:   if (flag) PetscCallExternal(HYPRE_ParCSRPilutSetFactorRowSize, jac->hsolver, jac->factorrowsize);
533:   PetscOptionsHeadEnd();
534:   PetscFunctionReturn(PETSC_SUCCESS);
535: }

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

542:   PetscFunctionBegin;
543:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
544:   if (iascii) {
545:     PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE Pilut preconditioning\n"));
546:     if (jac->maxiter != PETSC_DEFAULT) {
547:       PetscCall(PetscViewerASCIIPrintf(viewer, "    maximum number of iterations %" PetscInt_FMT "\n", jac->maxiter));
548:     } else {
549:       PetscCall(PetscViewerASCIIPrintf(viewer, "    default maximum number of iterations \n"));
550:     }
551:     if (jac->tol != PETSC_DEFAULT) {
552:       PetscCall(PetscViewerASCIIPrintf(viewer, "    drop tolerance %g\n", (double)jac->tol));
553:     } else {
554:       PetscCall(PetscViewerASCIIPrintf(viewer, "    default drop tolerance \n"));
555:     }
556:     if (jac->factorrowsize != PETSC_DEFAULT) {
557:       PetscCall(PetscViewerASCIIPrintf(viewer, "    factor row size %" PetscInt_FMT "\n", jac->factorrowsize));
558:     } else {
559:       PetscCall(PetscViewerASCIIPrintf(viewer, "    default factor row size \n"));
560:     }
561:   }
562:   PetscFunctionReturn(PETSC_SUCCESS);
563: }

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

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

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

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

584:   PetscCall(PetscOptionsBool("-pc_hypre_euclid_bj", "Use Block Jacobi for ILU in Euclid", "None", eu_bj, &eu_bj, &flag));
585:   if (flag) {
586:     jac->eu_bj = eu_bj ? 1 : 0;
587:     PetscCallExternal(HYPRE_EuclidSetBJ, jac->hsolver, jac->eu_bj);
588:   }
589:   PetscOptionsHeadEnd();
590:   PetscFunctionReturn(PETSC_SUCCESS);
591: }

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

598:   PetscFunctionBegin;
599:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
600:   if (iascii) {
601:     PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE Euclid preconditioning\n"));
602:     if (jac->eu_level != PETSC_DEFAULT) {
603:       PetscCall(PetscViewerASCIIPrintf(viewer, "    factorization levels %" PetscInt_FMT "\n", jac->eu_level));
604:     } else {
605:       PetscCall(PetscViewerASCIIPrintf(viewer, "    default factorization levels \n"));
606:     }
607:     PetscCall(PetscViewerASCIIPrintf(viewer, "    drop tolerance %g\n", (double)jac->eu_droptolerance));
608:     PetscCall(PetscViewerASCIIPrintf(viewer, "    use Block-Jacobi? %" PetscInt_FMT "\n", jac->eu_bj));
609:   }
610:   PetscFunctionReturn(PETSC_SUCCESS);
611: }

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

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

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

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

640:   PetscCall(VecHYPRE_IJVectorPopVec(hjac->x));
641:   PetscCall(VecHYPRE_IJVectorPopVec(hjac->b));
642:   PetscFunctionReturn(PETSC_SUCCESS);
643: }

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

650: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
651:   PetscFunctionBegin;
652:   if (jac->spgemm_type) {
653:     PetscCall(PetscStrcmp(jac->spgemm_type, name, &flag));
654:     PetscCheck(flag, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Cannot reset the HYPRE SpGEMM (really we can)");
655:     PetscFunctionReturn(PETSC_SUCCESS);
656:   } else {
657:     PetscCall(PetscStrallocpy(name, &jac->spgemm_type));
658:   }
659:   PetscCall(PetscStrcmp("cusparse", jac->spgemm_type, &flag));
660:   if (flag) {
661:     PetscCallExternal(HYPRE_SetSpGemmUseCusparse, 1);
662:     PetscFunctionReturn(PETSC_SUCCESS);
663:   }
664:   PetscCall(PetscStrcmp("hypre", jac->spgemm_type, &flag));
665:   if (flag) {
666:     PetscCallExternal(HYPRE_SetSpGemmUseCusparse, 0);
667:     PetscFunctionReturn(PETSC_SUCCESS);
668:   }
669:   jac->spgemm_type = NULL;
670:   SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown HYPRE SpGEM type %s; Choices are cusparse, hypre", name);
671: #endif
672: }

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

678:   PetscFunctionBegin;
680: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
681:   *spgemm = jac->spgemm_type;
682: #endif
683:   PetscFunctionReturn(PETSC_SUCCESS);
684: }

686: static const char *HYPREBoomerAMGCycleType[]   = {"", "V", "W"};
687: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP", "Ruge-Stueben", "", "modifiedRuge-Stueben", "", "", "Falgout", "", "PMIS", "", "HMIS"};
688: static const char *HYPREBoomerAMGMeasureType[] = {"local", "global"};
689: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
690: static const char *HYPREBoomerAMGSmoothType[] = {"Schwarz-smoothers", "Pilut", "ParaSails", "Euclid"};
691: 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"};
692: 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"};
693: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PC pc, PetscOptionItems *PetscOptionsObject)
694: {
695:   PC_HYPRE   *jac = (PC_HYPRE *)pc->data;
696:   PetscInt    bs, n, indx, level;
697:   PetscBool   flg, tmp_truth;
698:   double      tmpdbl, twodbl[2];
699:   const char *symtlist[]           = {"nonsymmetric", "SPD", "nonsymmetric,SPD"};
700:   const char *PCHYPRESpgemmTypes[] = {"cusparse", "hypre"};

702:   PetscFunctionBegin;
703:   PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE BoomerAMG Options");
704:   PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_cycle_type", "Cycle type", "None", HYPREBoomerAMGCycleType + 1, 2, HYPREBoomerAMGCycleType[jac->cycletype], &indx, &flg));
705:   if (flg) {
706:     jac->cycletype = indx + 1;
707:     PetscCallExternal(HYPRE_BoomerAMGSetCycleType, jac->hsolver, jac->cycletype);
708:   }
709:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_max_levels", "Number of levels (of grids) allowed", "None", jac->maxlevels, &jac->maxlevels, &flg));
710:   if (flg) {
711:     PetscCheck(jac->maxlevels >= 2, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Number of levels %" PetscInt_FMT " must be at least two", jac->maxlevels);
712:     PetscCallExternal(HYPRE_BoomerAMGSetMaxLevels, jac->hsolver, jac->maxlevels);
713:   }
714:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_max_iter", "Maximum iterations used PER hypre call", "None", jac->maxiter, &jac->maxiter, &flg));
715:   if (flg) {
716:     PetscCheck(jac->maxiter >= 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Number of iterations %" PetscInt_FMT " must be at least one", jac->maxiter);
717:     PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, jac->maxiter);
718:   }
719:   PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_tol", "Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)", "None", jac->tol, &jac->tol, &flg));
720:   if (flg) {
721:     PetscCheck(jac->tol >= 0.0, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_OUTOFRANGE, "Tolerance %g must be greater than or equal to zero", (double)jac->tol);
722:     PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, jac->tol);
723:   }
724:   bs = 1;
725:   if (pc->pmat) PetscCall(MatGetBlockSize(pc->pmat, &bs));
726:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_numfunctions", "Number of functions", "HYPRE_BoomerAMGSetNumFunctions", bs, &bs, &flg));
727:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetNumFunctions, jac->hsolver, bs);

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

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

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

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

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

762:   /* Grid sweeps */
763:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all", "Number of sweeps for the up and down grid levels", "None", jac->gridsweeps[0], &indx, &flg));
764:   if (flg) {
765:     PetscCallExternal(HYPRE_BoomerAMGSetNumSweeps, jac->hsolver, indx);
766:     /* modify the jac structure so we can view the updated options with PC_View */
767:     jac->gridsweeps[0] = indx;
768:     jac->gridsweeps[1] = indx;
769:     /*defaults coarse to 1 */
770:     jac->gridsweeps[2] = 1;
771:   }
772:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen", "Use a nodal based coarsening 1-6", "HYPRE_BoomerAMGSetNodal", jac->nodal_coarsening, &jac->nodal_coarsening, &flg));
773:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetNodal, jac->hsolver, jac->nodal_coarsening);
774:   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));
775:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetNodalDiag, jac->hsolver, jac->nodal_coarsening_diag);
776:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_variant", "Variant of algorithm 1-3", "HYPRE_BoomerAMGSetInterpVecVariant", jac->vec_interp_variant, &jac->vec_interp_variant, &flg));
777:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetInterpVecVariant, jac->hsolver, jac->vec_interp_variant);
778:   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));
779:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetInterpVecQMax, jac->hsolver, jac->vec_interp_qmax);
780:   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));
781:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetSmoothInterpVectors, jac->hsolver, jac->vec_interp_smooth);
782:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_interp_refine", "Preprocess the interpolation matrix through iterative weight refinement", "HYPRE_BoomerAMGSetInterpRefine", jac->interp_refine, &jac->interp_refine, &flg));
783:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetInterpRefine, jac->hsolver, jac->interp_refine);
784:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down", "Number of sweeps for the down cycles", "None", jac->gridsweeps[0], &indx, &flg));
785:   if (flg) {
786:     PetscCallExternal(HYPRE_BoomerAMGSetCycleNumSweeps, jac->hsolver, indx, 1);
787:     jac->gridsweeps[0] = indx;
788:   }
789:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up", "Number of sweeps for the up cycles", "None", jac->gridsweeps[1], &indx, &flg));
790:   if (flg) {
791:     PetscCallExternal(HYPRE_BoomerAMGSetCycleNumSweeps, jac->hsolver, indx, 2);
792:     jac->gridsweeps[1] = indx;
793:   }
794:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse", "Number of sweeps for the coarse level", "None", jac->gridsweeps[2], &indx, &flg));
795:   if (flg) {
796:     PetscCallExternal(HYPRE_BoomerAMGSetCycleNumSweeps, jac->hsolver, indx, 3);
797:     jac->gridsweeps[2] = indx;
798:   }

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

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

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

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

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

838:   /* Relax type */
839:   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));
840:   if (flg) {
841:     jac->relaxtype[0] = jac->relaxtype[1] = indx;
842:     PetscCallExternal(HYPRE_BoomerAMGSetRelaxType, jac->hsolver, indx);
843:     /* by default, coarse type set to 9 */
844:     jac->relaxtype[2] = 9;
845:     PetscCallExternal(HYPRE_BoomerAMGSetCycleRelaxType, jac->hsolver, 9, 3);
846:   }
847:   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));
848:   if (flg) {
849:     jac->relaxtype[0] = indx;
850:     PetscCallExternal(HYPRE_BoomerAMGSetCycleRelaxType, jac->hsolver, indx, 1);
851:   }
852:   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));
853:   if (flg) {
854:     jac->relaxtype[1] = indx;
855:     PetscCallExternal(HYPRE_BoomerAMGSetCycleRelaxType, jac->hsolver, indx, 2);
856:   }
857:   PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse", "Relax type on coarse grid", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), HYPREBoomerAMGRelaxType[9], &indx, &flg));
858:   if (flg) {
859:     jac->relaxtype[2] = indx;
860:     PetscCallExternal(HYPRE_BoomerAMGSetCycleRelaxType, jac->hsolver, indx, 3);
861:   }

863:   /* Relaxation Weight */
864:   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));
865:   if (flg) {
866:     PetscCallExternal(HYPRE_BoomerAMGSetRelaxWt, jac->hsolver, tmpdbl);
867:     jac->relaxweight = tmpdbl;
868:   }

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

879:   /* Outer relaxation Weight */
880:   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));
881:   if (flg) {
882:     PetscCallExternal(HYPRE_BoomerAMGSetOuterWt, jac->hsolver, tmpdbl);
883:     jac->outerrelaxweight = tmpdbl;
884:   }

886:   n         = 2;
887:   twodbl[0] = twodbl[1] = 1.0;
888:   PetscCall(PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level", "Set the outer relaxation weight for a particular level (weight,level)", "None", twodbl, &n, &flg));
889:   if (flg) {
890:     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);
891:     indx = (int)PetscAbsReal(twodbl[1]);
892:     PetscCallExternal(HYPRE_BoomerAMGSetLevelOuterWt, jac->hsolver, twodbl[0], indx);
893:   }

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

898:   if (flg && tmp_truth) {
899:     jac->relaxorder = 0;
900:     PetscCallExternal(HYPRE_BoomerAMGSetRelaxOrder, jac->hsolver, jac->relaxorder);
901:   }
902:   PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_measure_type", "Measure type", "None", HYPREBoomerAMGMeasureType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGMeasureType), HYPREBoomerAMGMeasureType[0], &indx, &flg));
903:   if (flg) {
904:     jac->measuretype = indx;
905:     PetscCallExternal(HYPRE_BoomerAMGSetMeasureType, jac->hsolver, jac->measuretype);
906:   }
907:   /* update list length 3/07 */
908:   PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type", "Coarsen type", "None", HYPREBoomerAMGCoarsenType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGCoarsenType), HYPREBoomerAMGCoarsenType[6], &indx, &flg));
909:   if (flg) {
910:     jac->coarsentype = indx;
911:     PetscCallExternal(HYPRE_BoomerAMGSetCoarsenType, jac->hsolver, jac->coarsentype);
912:   }

914:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_max_coarse_size", "Maximum size of coarsest grid", "None", jac->maxc, &jac->maxc, &flg));
915:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetMaxCoarseSize, jac->hsolver, jac->maxc);
916:   PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_min_coarse_size", "Minimum size of coarsest grid", "None", jac->minc, &jac->minc, &flg));
917:   if (flg) PetscCallExternal(HYPRE_BoomerAMGSetMinCoarseSize, jac->hsolver, jac->minc);
918: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
919:   // global parameter but is closely associated with BoomerAMG
920:   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));
921:   if (!flg) indx = 0;
922:   PetscCall(PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG(pc, PCHYPRESpgemmTypes[indx]));
923: #endif
924:   /* AIR */
925: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
926:   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));
927:   PetscCallExternal(HYPRE_BoomerAMGSetRestriction, jac->hsolver, jac->Rtype);
928:   if (jac->Rtype) {
929:     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 */

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

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

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

940:     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));
941:     PetscCallExternal(HYPRE_BoomerAMGSetADropType, jac->hsolver, jac->Adroptype);
942:   }
943: #endif

945: #if PETSC_PKG_HYPRE_VERSION_LE(9, 9, 9)
946:   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);
947: #endif

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

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

961:     jac->printstatistics = PETSC_TRUE;
962:     PetscCallExternal(HYPRE_BoomerAMGSetPrintLevel, jac->hsolver, level);
963:   }

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

970:     jac->printstatistics = PETSC_TRUE;
971:     PetscCallExternal(HYPRE_BoomerAMGSetDebugFlag, jac->hsolver, level);
972:   }

974:   PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg));
975:   if (flg && tmp_truth) {
976:     PetscInt tmp_int;
977:     PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", jac->nodal_relax_levels, &tmp_int, &flg));
978:     if (flg) jac->nodal_relax_levels = tmp_int;
979:     PetscCallExternal(HYPRE_BoomerAMGSetSmoothType, jac->hsolver, 6);
980:     PetscCallExternal(HYPRE_BoomerAMGSetDomainType, jac->hsolver, 1);
981:     PetscCallExternal(HYPRE_BoomerAMGSetOverlap, jac->hsolver, 0);
982:     PetscCallExternal(HYPRE_BoomerAMGSetSmoothNumLevels, jac->hsolver, jac->nodal_relax_levels);
983:   }

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

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

995:   PetscOptionsHeadEnd();
996:   PetscFunctionReturn(PETSC_SUCCESS);
997: }

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

1004:   PetscFunctionBegin;
1005:   PetscCall(PetscCitationsRegister(hypreCitation, &cite));
1006:   PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, its * jac->maxiter);
1007:   PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, rtol);
1008:   jac->applyrichardson = PETSC_TRUE;
1009:   PetscCall(PCApply_HYPRE(pc, b, y));
1010:   jac->applyrichardson = PETSC_FALSE;
1011:   PetscCallExternal(HYPRE_BoomerAMGGetNumIterations, jac->hsolver, &oits);
1012:   *outits = oits;
1013:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1014:   else *reason = PCRICHARDSON_CONVERGED_RTOL;
1015:   PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, jac->tol);
1016:   PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, jac->maxiter);
1017:   PetscFunctionReturn(PETSC_SUCCESS);
1018: }

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

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

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

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

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

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

1053:     if (jac->relaxorder) {
1054:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Using CF-relaxation\n"));
1055:     } else {
1056:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Not using CF-relaxation\n"));
1057:     }
1058:     if (jac->smoothtype != -1) {
1059:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Smooth type          %s\n", HYPREBoomerAMGSmoothType[jac->smoothtype]));
1060:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Smooth num levels    %" PetscInt_FMT "\n", jac->smoothnumlevels));
1061:     } else {
1062:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Not using more complex smoothers.\n"));
1063:     }
1064:     if (jac->smoothtype == 3) {
1065:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Euclid ILU(k) levels %" PetscInt_FMT "\n", jac->eu_level));
1066:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Euclid ILU(k) drop tolerance %g\n", (double)jac->eu_droptolerance));
1067:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Euclid ILU use Block-Jacobi? %" PetscInt_FMT "\n", jac->eu_bj));
1068:     }
1069:     PetscCall(PetscViewerASCIIPrintf(viewer, "    Measure type        %s\n", HYPREBoomerAMGMeasureType[jac->measuretype]));
1070:     PetscCall(PetscViewerASCIIPrintf(viewer, "    Coarsen type        %s\n", HYPREBoomerAMGCoarsenType[jac->coarsentype]));
1071:     PetscCall(PetscViewerASCIIPrintf(viewer, "    Interpolation type  %s\n", jac->interptype != 100 ? HYPREBoomerAMGInterpType[jac->interptype] : "1pt"));
1072:     if (jac->nodal_coarsening) PetscCall(PetscViewerASCIIPrintf(viewer, "    Using nodal coarsening with HYPRE_BOOMERAMGSetNodal() %" PetscInt_FMT "\n", jac->nodal_coarsening));
1073:     if (jac->vec_interp_variant) {
1074:       PetscCall(PetscViewerASCIIPrintf(viewer, "    HYPRE_BoomerAMGSetInterpVecVariant() %" PetscInt_FMT "\n", jac->vec_interp_variant));
1075:       PetscCall(PetscViewerASCIIPrintf(viewer, "    HYPRE_BoomerAMGSetInterpVecQMax() %" PetscInt_FMT "\n", jac->vec_interp_qmax));
1076:       PetscCall(PetscViewerASCIIPrintf(viewer, "    HYPRE_BoomerAMGSetSmoothInterpVectors() %d\n", jac->vec_interp_smooth));
1077:     }
1078:     if (jac->nodal_relax) PetscCall(PetscViewerASCIIPrintf(viewer, "    Using nodal relaxation via Schwarz smoothing on levels %" PetscInt_FMT "\n", jac->nodal_relax_levels));
1079: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
1080:     PetscCall(PetscViewerASCIIPrintf(viewer, "    SpGEMM type         %s\n", jac->spgemm_type));
1081: #endif
1082:     /* AIR */
1083:     if (jac->Rtype) {
1084:       PetscCall(PetscViewerASCIIPrintf(viewer, "    Using approximate ideal restriction type %" PetscInt_FMT "\n", jac->Rtype));
1085:       PetscCall(PetscViewerASCIIPrintf(viewer, "      Threshold for R %g\n", (double)jac->Rstrongthreshold));
1086:       PetscCall(PetscViewerASCIIPrintf(viewer, "      Filter for R %g\n", (double)jac->Rfilterthreshold));
1087:       PetscCall(PetscViewerASCIIPrintf(viewer, "      A drop tolerance %g\n", (double)jac->Adroptol));
1088:       PetscCall(PetscViewerASCIIPrintf(viewer, "      A drop type %" PetscInt_FMT "\n", jac->Adroptype));
1089:     }
1090:   }
1091:   PetscFunctionReturn(PETSC_SUCCESS);
1092: }

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

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

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

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

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

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

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

1125:   PetscOptionsHeadEnd();
1126:   PetscFunctionReturn(PETSC_SUCCESS);
1127: }

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

1135:   PetscFunctionBegin;
1136:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1137:   if (iascii) {
1138:     PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE ParaSails preconditioning\n"));
1139:     PetscCall(PetscViewerASCIIPrintf(viewer, "    nlevels %" PetscInt_FMT "\n", jac->nlevels));
1140:     PetscCall(PetscViewerASCIIPrintf(viewer, "    threshold %g\n", (double)jac->threshold));
1141:     PetscCall(PetscViewerASCIIPrintf(viewer, "    filter %g\n", (double)jac->filter));
1142:     PetscCall(PetscViewerASCIIPrintf(viewer, "    load balance %g\n", (double)jac->loadbal));
1143:     PetscCall(PetscViewerASCIIPrintf(viewer, "    reuse nonzero structure %s\n", PetscBools[jac->ruse]));
1144:     PetscCall(PetscViewerASCIIPrintf(viewer, "    print info to screen %s\n", PetscBools[jac->logging]));
1145:     if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
1146:     else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
1147:     else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
1148:     else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Unknown HYPRE ParaSails symmetric option %" PetscInt_FMT, jac->symt);
1149:     PetscCall(PetscViewerASCIIPrintf(viewer, "    %s\n", symt));
1150:   }
1151:   PetscFunctionReturn(PETSC_SUCCESS);
1152: }

1154: static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PC pc, PetscOptionItems *PetscOptionsObject)
1155: {
1156:   PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1157:   PetscInt  n;
1158:   PetscBool flag, flag2, flag3, flag4;

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

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

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

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

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

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

1302:   PetscFunctionBegin;
1303:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
1304:   if (iascii) {
1305:     PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE ADS preconditioning\n"));
1306:     PetscCall(PetscViewerASCIIPrintf(viewer, "    subspace iterations per application %" PetscInt_FMT "\n", jac->as_max_iter));
1307:     PetscCall(PetscViewerASCIIPrintf(viewer, "    subspace cycle type %" PetscInt_FMT "\n", jac->ads_cycle_type));
1308:     PetscCall(PetscViewerASCIIPrintf(viewer, "    subspace iteration tolerance %g\n", (double)jac->as_tol));
1309:     PetscCall(PetscViewerASCIIPrintf(viewer, "    smoother type %" PetscInt_FMT "\n", jac->as_relax_type));
1310:     PetscCall(PetscViewerASCIIPrintf(viewer, "    number of smoothing steps %" PetscInt_FMT "\n", jac->as_relax_times));
1311:     PetscCall(PetscViewerASCIIPrintf(viewer, "    smoother weight %g\n", (double)jac->as_relax_weight));
1312:     PetscCall(PetscViewerASCIIPrintf(viewer, "    smoother omega %g\n", (double)jac->as_omega));
1313:     PetscCall(PetscViewerASCIIPrintf(viewer, "    AMS solver using boomerAMG\n"));
1314:     PetscCall(PetscViewerASCIIPrintf(viewer, "        subspace cycle type %" PetscInt_FMT "\n", jac->ams_cycle_type));
1315:     PetscCall(PetscViewerASCIIPrintf(viewer, "        coarsening type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[0]));
1316:     PetscCall(PetscViewerASCIIPrintf(viewer, "        levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[1]));
1317:     PetscCall(PetscViewerASCIIPrintf(viewer, "        relaxation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[2]));
1318:     PetscCall(PetscViewerASCIIPrintf(viewer, "        interpolation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[3]));
1319:     PetscCall(PetscViewerASCIIPrintf(viewer, "        max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[4]));
1320:     PetscCall(PetscViewerASCIIPrintf(viewer, "        strength threshold %g\n", (double)jac->as_amg_alpha_theta));
1321:     PetscCall(PetscViewerASCIIPrintf(viewer, "    vector Poisson solver using boomerAMG\n"));
1322:     PetscCall(PetscViewerASCIIPrintf(viewer, "        coarsening type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[0]));
1323:     PetscCall(PetscViewerASCIIPrintf(viewer, "        levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_beta_opts[1]));
1324:     PetscCall(PetscViewerASCIIPrintf(viewer, "        relaxation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[2]));
1325:     PetscCall(PetscViewerASCIIPrintf(viewer, "        interpolation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[3]));
1326:     PetscCall(PetscViewerASCIIPrintf(viewer, "        max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_beta_opts[4]));
1327:     PetscCall(PetscViewerASCIIPrintf(viewer, "        strength threshold %g\n", (double)jac->as_amg_beta_theta));
1328:   }
1329:   PetscFunctionReturn(PETSC_SUCCESS);
1330: }

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

1337:   PetscFunctionBegin;
1338:   PetscCall(PetscObjectTypeCompare((PetscObject)G, MATHYPRE, &ishypre));
1339:   if (ishypre) {
1340:     PetscCall(PetscObjectReference((PetscObject)G));
1341:     PetscCall(MatDestroy(&jac->G));
1342:     jac->G = G;
1343:   } else {
1344:     PetscCall(MatDestroy(&jac->G));
1345:     PetscCall(MatConvert(G, MATHYPRE, MAT_INITIAL_MATRIX, &jac->G));
1346:   }
1347:   PetscFunctionReturn(PETSC_SUCCESS);
1348: }

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

1353:    Collective

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

1359:    Level: intermediate

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

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

1366:    Developer Note:
1367:    This automatically converts the matrix to `MATHYPRE` if it is not already of that type

1369: .seealso: `PCHYPRE`, `PCHYPRESetDiscreteCurl()`
1370: @*/
1371: PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G)
1372: {
1373:   PetscFunctionBegin;
1376:   PetscCheckSameComm(pc, 1, G, 2);
1377:   PetscTryMethod(pc, "PCHYPRESetDiscreteGradient_C", (PC, Mat), (pc, G));
1378:   PetscFunctionReturn(PETSC_SUCCESS);
1379: }

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

1386:   PetscFunctionBegin;
1387:   PetscCall(PetscObjectTypeCompare((PetscObject)C, MATHYPRE, &ishypre));
1388:   if (ishypre) {
1389:     PetscCall(PetscObjectReference((PetscObject)C));
1390:     PetscCall(MatDestroy(&jac->C));
1391:     jac->C = C;
1392:   } else {
1393:     PetscCall(MatDestroy(&jac->C));
1394:     PetscCall(MatConvert(C, MATHYPRE, MAT_INITIAL_MATRIX, &jac->C));
1395:   }
1396:   PetscFunctionReturn(PETSC_SUCCESS);
1397: }

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

1402:    Collective

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

1408:    Level: intermediate

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

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

1415:    Developer Note:
1416:    This automatically converts the matrix to `MATHYPRE` if it is not already of that type

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

1420: .seealso: `PCHYPRE`, `PCHYPRESetDiscreteGradient()`
1421: @*/
1422: PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C)
1423: {
1424:   PetscFunctionBegin;
1427:   PetscCheckSameComm(pc, 1, C, 2);
1428:   PetscTryMethod(pc, "PCHYPRESetDiscreteCurl_C", (PC, Mat), (pc, C));
1429:   PetscFunctionReturn(PETSC_SUCCESS);
1430: }

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

1439:   PetscCall(MatDestroy(&jac->RT_PiFull));
1440:   PetscCall(MatDestroy(&jac->ND_PiFull));
1441:   for (i = 0; i < 3; ++i) {
1442:     PetscCall(MatDestroy(&jac->RT_Pi[i]));
1443:     PetscCall(MatDestroy(&jac->ND_Pi[i]));
1444:   }

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

1492:   PetscFunctionReturn(PETSC_SUCCESS);
1493: }

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

1498:    Collective

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

1508:    Level: intermediate

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

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

1515:    Developer Note:
1516:    This automatically converts the matrix to `MATHYPRE` if it is not already of that type

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

1524:   PetscFunctionBegin;
1526:   if (RT_PiFull) {
1528:     PetscCheckSameComm(pc, 1, RT_PiFull, 3);
1529:   }
1530:   if (RT_Pi) {
1532:     for (i = 0; i < dim; ++i) {
1533:       if (RT_Pi[i]) {
1535:         PetscCheckSameComm(pc, 1, RT_Pi[i], 4);
1536:       }
1537:     }
1538:   }
1539:   if (ND_PiFull) {
1541:     PetscCheckSameComm(pc, 1, ND_PiFull, 5);
1542:   }
1543:   if (ND_Pi) {
1545:     for (i = 0; i < dim; ++i) {
1546:       if (ND_Pi[i]) {
1548:         PetscCheckSameComm(pc, 1, ND_Pi[i], 6);
1549:       }
1550:     }
1551:   }
1552:   PetscTryMethod(pc, "PCHYPRESetInterpolations_C", (PC, PetscInt, Mat, Mat[], Mat, Mat[]), (pc, dim, RT_PiFull, RT_Pi, ND_PiFull, ND_Pi));
1553:   PetscFunctionReturn(PETSC_SUCCESS);
1554: }

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

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

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

1597:    Collective

1599:    Input Parameters:
1600: +  pc - the preconditioning context
1601: -  A - the matrix

1603:    Level: intermediate

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

1608:    Developer Note:
1609:    This automatically converts the matrix to `MATHYPRE` if it is not already of that type

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

1613: .seealso: `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetBetaPoissonMatrix()`
1614: @*/
1615: PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A)
1616: {
1617:   PetscFunctionBegin;
1620:   PetscCheckSameComm(pc, 1, A, 2);
1621:   PetscTryMethod(pc, "PCHYPRESetPoissonMatrix_C", (PC, Mat, PetscBool), (pc, A, PETSC_TRUE));
1622:   PetscFunctionReturn(PETSC_SUCCESS);
1623: }

1625: /*@
1626:    PCHYPRESetBetaPoissonMatrix - Set Poisson matrix for `PCHYPRE` of type ams

1628:    Collective

1630:    Input Parameters:
1631: +  pc - the preconditioning context
1632: -  A - the matrix, or NULL to turn it off

1634:    Level: intermediate

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

1639:    Developer Note:
1640:    This automatically converts the matrix to `MATHYPRE` if it is not already of that type

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

1644: .seealso: `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetAlphaPoissonMatrix()`
1645: @*/
1646: PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A)
1647: {
1648:   PetscFunctionBegin;
1650:   if (A) {
1652:     PetscCheckSameComm(pc, 1, A, 2);
1653:   }
1654:   PetscTryMethod(pc, "PCHYPRESetPoissonMatrix_C", (PC, Mat, PetscBool), (pc, A, PETSC_FALSE));
1655:   PetscFunctionReturn(PETSC_SUCCESS);
1656: }

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

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

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

1683:    Collective

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

1691:    Level: intermediate

1693:    Developer Note:
1694:    If this is only for  `PCHYPRE` type of ams it should be called `PCHYPREAMSSetEdgeConstantVectors()`

1696: .seealso: `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetAlphaPoissonMatrix()`
1697: @*/
1698: PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo)
1699: {
1700:   PetscFunctionBegin;
1705:   PetscCheckSameComm(pc, 1, ozz, 2);
1706:   PetscCheckSameComm(pc, 1, zoz, 3);
1707:   if (zzo) PetscCheckSameComm(pc, 1, zzo, 4);
1708:   PetscTryMethod(pc, "PCHYPRESetEdgeConstantVectors_C", (PC, Vec, Vec, Vec), (pc, ozz, zoz, zzo));
1709:   PetscFunctionReturn(PETSC_SUCCESS);
1710: }

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

1716:   PetscFunctionBegin;
1717:   PetscCall(VecHYPRE_IJVectorDestroy(&jac->interior));
1718:   PetscCall(VecHYPRE_IJVectorCreate(interior->map, &jac->interior));
1719:   PetscCall(VecHYPRE_IJVectorCopy(interior, jac->interior));
1720:   jac->ams_beta_is_zero_part = PETSC_TRUE;
1721:   PetscFunctionReturn(PETSC_SUCCESS);
1722: }

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

1727:    Collective

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

1733:    Level: intermediate

1735:    Note:
1736:    This calls `HYPRE_AMSSetInteriorNodes()`

1738:    Developer Note:
1739:    If this is only for  `PCHYPRE` type of ams it should be called `PCHYPREAMSSetInteriorNodes()`

1741: .seealso: `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetAlphaPoissonMatrix()`
1742: @*/
1743: PetscErrorCode PCHYPREAMSSetInteriorNodes(PC pc, Vec interior)
1744: {
1745:   PetscFunctionBegin;
1748:   PetscCheckSameComm(pc, 1, interior, 2);
1749:   PetscTryMethod(pc, "PCHYPREAMSSetInteriorNodes_C", (PC, Vec), (pc, interior));
1750:   PetscFunctionReturn(PETSC_SUCCESS);
1751: }

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

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

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

1774:     PetscCall(VecHYPRE_IJVectorCreate(tv->map, &jac->coords[i]));
1775:     PetscCall(VecGetArrayWrite(tv, &array));
1776:     for (j = 0; j < nloc; j++) array[j] = coords[j * dim + i];
1777:     PetscCall(VecRestoreArrayWrite(tv, &array));
1778:     PetscCall(VecHYPRE_IJVectorCopy(tv, jac->coords[i]));
1779:   }
1780:   PetscCall(VecDestroy(&tv));
1781:   PetscFunctionReturn(PETSC_SUCCESS);
1782: }

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

1788:   PetscFunctionBegin;
1789:   *name = jac->hypre_type;
1790:   PetscFunctionReturn(PETSC_SUCCESS);
1791: }

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

1798:   PetscFunctionBegin;
1799:   if (jac->hypre_type) {
1800:     PetscCall(PetscStrcmp(jac->hypre_type, name, &flag));
1801:     PetscCheck(flag, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "Cannot reset the HYPRE preconditioner type once it has been set");
1802:     PetscFunctionReturn(PETSC_SUCCESS);
1803:   } else {
1804:     PetscCall(PetscStrallocpy(name, &jac->hypre_type));
1805:   }

1807:   jac->maxiter         = PETSC_DEFAULT;
1808:   jac->tol             = PETSC_DEFAULT;
1809:   jac->printstatistics = PetscLogPrintInfo;

1811:   PetscCall(PetscStrcmp("pilut", jac->hypre_type, &flag));
1812:   if (flag) {
1813:     PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
1814:     PetscCallExternal(HYPRE_ParCSRPilutCreate, jac->comm_hypre, &jac->hsolver);
1815:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
1816:     pc->ops->view           = PCView_HYPRE_Pilut;
1817:     jac->destroy            = HYPRE_ParCSRPilutDestroy;
1818:     jac->setup              = HYPRE_ParCSRPilutSetup;
1819:     jac->solve              = HYPRE_ParCSRPilutSolve;
1820:     jac->factorrowsize      = PETSC_DEFAULT;
1821:     PetscFunctionReturn(PETSC_SUCCESS);
1822:   }
1823:   PetscCall(PetscStrcmp("euclid", jac->hypre_type, &flag));
1824:   if (flag) {
1825: #if defined(PETSC_USE_64BIT_INDICES)
1826:     SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Hypre Euclid does not support 64-bit indices");
1827: #endif
1828:     PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
1829:     PetscCallExternal(HYPRE_EuclidCreate, jac->comm_hypre, &jac->hsolver);
1830:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
1831:     pc->ops->view           = PCView_HYPRE_Euclid;
1832:     jac->destroy            = HYPRE_EuclidDestroy;
1833:     jac->setup              = HYPRE_EuclidSetup;
1834:     jac->solve              = HYPRE_EuclidSolve;
1835:     jac->factorrowsize      = PETSC_DEFAULT;
1836:     jac->eu_level           = PETSC_DEFAULT; /* default */
1837:     PetscFunctionReturn(PETSC_SUCCESS);
1838:   }
1839:   PetscCall(PetscStrcmp("parasails", jac->hypre_type, &flag));
1840:   if (flag) {
1841:     PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
1842:     PetscCallExternal(HYPRE_ParaSailsCreate, jac->comm_hypre, &jac->hsolver);
1843:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
1844:     pc->ops->view           = PCView_HYPRE_ParaSails;
1845:     jac->destroy            = HYPRE_ParaSailsDestroy;
1846:     jac->setup              = HYPRE_ParaSailsSetup;
1847:     jac->solve              = HYPRE_ParaSailsSolve;
1848:     /* initialize */
1849:     jac->nlevels   = 1;
1850:     jac->threshold = .1;
1851:     jac->filter    = .1;
1852:     jac->loadbal   = 0;
1853:     if (PetscLogPrintInfo) jac->logging = (int)PETSC_TRUE;
1854:     else jac->logging = (int)PETSC_FALSE;

1856:     jac->ruse = (int)PETSC_FALSE;
1857:     jac->symt = 0;
1858:     PetscCallExternal(HYPRE_ParaSailsSetParams, jac->hsolver, jac->threshold, jac->nlevels);
1859:     PetscCallExternal(HYPRE_ParaSailsSetFilter, jac->hsolver, jac->filter);
1860:     PetscCallExternal(HYPRE_ParaSailsSetLoadbal, jac->hsolver, jac->loadbal);
1861:     PetscCallExternal(HYPRE_ParaSailsSetLogging, jac->hsolver, jac->logging);
1862:     PetscCallExternal(HYPRE_ParaSailsSetReuse, jac->hsolver, jac->ruse);
1863:     PetscCallExternal(HYPRE_ParaSailsSetSym, jac->hsolver, jac->symt);
1864:     PetscFunctionReturn(PETSC_SUCCESS);
1865:   }
1866:   PetscCall(PetscStrcmp("boomeramg", jac->hypre_type, &flag));
1867:   if (flag) {
1868:     PetscCallExternal(HYPRE_BoomerAMGCreate, &jac->hsolver);
1869:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_BoomerAMG;
1870:     pc->ops->view            = PCView_HYPRE_BoomerAMG;
1871:     pc->ops->applytranspose  = PCApplyTranspose_HYPRE_BoomerAMG;
1872:     pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
1873:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", PCGetInterpolations_BoomerAMG));
1874:     PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", PCGetCoarseOperators_BoomerAMG));
1875:     jac->destroy         = HYPRE_BoomerAMGDestroy;
1876:     jac->setup           = HYPRE_BoomerAMGSetup;
1877:     jac->solve           = HYPRE_BoomerAMGSolve;
1878:     jac->applyrichardson = PETSC_FALSE;
1879:     /* these defaults match the hypre defaults */
1880:     jac->cycletype       = 1;
1881:     jac->maxlevels       = 25;
1882:     jac->maxiter         = 1;
1883:     jac->tol             = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
1884:     jac->truncfactor     = 0.0;
1885:     jac->strongthreshold = .25;
1886:     jac->maxrowsum       = .9;
1887:     jac->coarsentype     = 6;
1888:     jac->measuretype     = 0;
1889:     jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
1890:     jac->smoothtype                                              = -1; /* Not set by default */
1891:     jac->smoothnumlevels                                         = 25;
1892:     jac->eu_level                                                = 0;
1893:     jac->eu_droptolerance                                        = 0;
1894:     jac->eu_bj                                                   = 0;
1895:     jac->relaxtype[0] = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a PC - most likely with CG */
1896:     jac->relaxtype[2]                     = 9; /*G.E. */
1897:     jac->relaxweight                      = 1.0;
1898:     jac->outerrelaxweight                 = 1.0;
1899:     jac->relaxorder                       = 1;
1900:     jac->interptype                       = 0;
1901:     jac->Rtype                            = 0;
1902:     jac->Rstrongthreshold                 = 0.25;
1903:     jac->Rfilterthreshold                 = 0.0;
1904:     jac->Adroptype                        = -1;
1905:     jac->Adroptol                         = 0.0;
1906:     jac->agg_nl                           = 0;
1907:     jac->agg_interptype                   = 4;
1908:     jac->pmax                             = 0;
1909:     jac->truncfactor                      = 0.0;
1910:     jac->agg_num_paths                    = 1;
1911:     jac->maxc                             = 9;
1912:     jac->minc                             = 1;
1913:     jac->nodal_coarsening                 = 0;
1914:     jac->nodal_coarsening_diag            = 0;
1915:     jac->vec_interp_variant               = 0;
1916:     jac->vec_interp_qmax                  = 0;
1917:     jac->vec_interp_smooth                = PETSC_FALSE;
1918:     jac->interp_refine                    = 0;
1919:     jac->nodal_relax                      = PETSC_FALSE;
1920:     jac->nodal_relax_levels               = 1;
1921:     jac->rap2                             = 0;

1923:     /* GPU defaults
1924:          from https://hypre.readthedocs.io/en/latest/solvers-boomeramg.html#gpu-supported-options
1925:          and /src/parcsr_ls/par_amg.c */
1926: #if defined(PETSC_HAVE_HYPRE_DEVICE)
1927:     jac->keeptranspose  = PETSC_TRUE;
1928:     jac->mod_rap2       = 1;
1929:     jac->coarsentype    = 8;
1930:     jac->relaxorder     = 0;
1931:     jac->interptype     = 6;
1932:     jac->relaxtype[0]   = 18;
1933:     jac->relaxtype[1]   = 18;
1934:     jac->agg_interptype = 7;
1935: #else
1936:     jac->keeptranspose = PETSC_FALSE;
1937:     jac->mod_rap2      = 0;
1938: #endif
1939:     PetscCallExternal(HYPRE_BoomerAMGSetCycleType, jac->hsolver, jac->cycletype);
1940:     PetscCallExternal(HYPRE_BoomerAMGSetMaxLevels, jac->hsolver, jac->maxlevels);
1941:     PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, jac->maxiter);
1942:     PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, jac->tol);
1943:     PetscCallExternal(HYPRE_BoomerAMGSetTruncFactor, jac->hsolver, jac->truncfactor);
1944:     PetscCallExternal(HYPRE_BoomerAMGSetStrongThreshold, jac->hsolver, jac->strongthreshold);
1945:     PetscCallExternal(HYPRE_BoomerAMGSetMaxRowSum, jac->hsolver, jac->maxrowsum);
1946:     PetscCallExternal(HYPRE_BoomerAMGSetCoarsenType, jac->hsolver, jac->coarsentype);
1947:     PetscCallExternal(HYPRE_BoomerAMGSetMeasureType, jac->hsolver, jac->measuretype);
1948:     PetscCallExternal(HYPRE_BoomerAMGSetRelaxOrder, jac->hsolver, jac->relaxorder);
1949:     PetscCallExternal(HYPRE_BoomerAMGSetInterpType, jac->hsolver, jac->interptype);
1950:     PetscCallExternal(HYPRE_BoomerAMGSetAggNumLevels, jac->hsolver, jac->agg_nl);
1951:     PetscCallExternal(HYPRE_BoomerAMGSetAggInterpType, jac->hsolver, jac->agg_interptype);
1952:     PetscCallExternal(HYPRE_BoomerAMGSetPMaxElmts, jac->hsolver, jac->pmax);
1953:     PetscCallExternal(HYPRE_BoomerAMGSetNumPaths, jac->hsolver, jac->agg_num_paths);
1954:     PetscCallExternal(HYPRE_BoomerAMGSetRelaxType, jac->hsolver, jac->relaxtype[0]);  /* defaults coarse to 9 */
1955:     PetscCallExternal(HYPRE_BoomerAMGSetNumSweeps, jac->hsolver, jac->gridsweeps[0]); /* defaults coarse to 1 */
1956:     PetscCallExternal(HYPRE_BoomerAMGSetMaxCoarseSize, jac->hsolver, jac->maxc);
1957:     PetscCallExternal(HYPRE_BoomerAMGSetMinCoarseSize, jac->hsolver, jac->minc);
1958:     /* GPU */
1959: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
1960:     PetscCallExternal(HYPRE_BoomerAMGSetKeepTranspose, jac->hsolver, jac->keeptranspose ? 1 : 0);
1961:     PetscCallExternal(HYPRE_BoomerAMGSetRAP2, jac->hsolver, jac->rap2);
1962:     PetscCallExternal(HYPRE_BoomerAMGSetModuleRAP2, jac->hsolver, jac->mod_rap2);
1963: #endif

1965:     /* AIR */
1966: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
1967:     PetscCallExternal(HYPRE_BoomerAMGSetRestriction, jac->hsolver, jac->Rtype);
1968:     PetscCallExternal(HYPRE_BoomerAMGSetStrongThresholdR, jac->hsolver, jac->Rstrongthreshold);
1969:     PetscCallExternal(HYPRE_BoomerAMGSetFilterThresholdR, jac->hsolver, jac->Rfilterthreshold);
1970:     PetscCallExternal(HYPRE_BoomerAMGSetADropTol, jac->hsolver, jac->Adroptol);
1971:     PetscCallExternal(HYPRE_BoomerAMGSetADropType, jac->hsolver, jac->Adroptype);
1972: #endif
1973:     PetscFunctionReturn(PETSC_SUCCESS);
1974:   }
1975:   PetscCall(PetscStrcmp("ams", jac->hypre_type, &flag));
1976:   if (flag) {
1977:     PetscCallExternal(HYPRE_AMSCreate, &jac->hsolver);
1978:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_AMS;
1979:     pc->ops->view           = PCView_HYPRE_AMS;
1980:     jac->destroy            = HYPRE_AMSDestroy;
1981:     jac->setup              = HYPRE_AMSSetup;
1982:     jac->solve              = HYPRE_AMSSolve;
1983:     jac->coords[0]          = NULL;
1984:     jac->coords[1]          = NULL;
1985:     jac->coords[2]          = NULL;
1986:     jac->interior           = NULL;
1987:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */
1988:     jac->as_print       = 0;
1989:     jac->as_max_iter    = 1;  /* used as a preconditioner */
1990:     jac->as_tol         = 0.; /* used as a preconditioner */
1991:     jac->ams_cycle_type = 13;
1992:     /* Smoothing options */
1993:     jac->as_relax_type   = 2;
1994:     jac->as_relax_times  = 1;
1995:     jac->as_relax_weight = 1.0;
1996:     jac->as_omega        = 1.0;
1997:     /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1998:     jac->as_amg_alpha_opts[0] = 10;
1999:     jac->as_amg_alpha_opts[1] = 1;
2000:     jac->as_amg_alpha_opts[2] = 6;
2001:     jac->as_amg_alpha_opts[3] = 6;
2002:     jac->as_amg_alpha_opts[4] = 4;
2003:     jac->as_amg_alpha_theta   = 0.25;
2004:     /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
2005:     jac->as_amg_beta_opts[0] = 10;
2006:     jac->as_amg_beta_opts[1] = 1;
2007:     jac->as_amg_beta_opts[2] = 6;
2008:     jac->as_amg_beta_opts[3] = 6;
2009:     jac->as_amg_beta_opts[4] = 4;
2010:     jac->as_amg_beta_theta   = 0.25;
2011:     PetscCallExternal(HYPRE_AMSSetPrintLevel, jac->hsolver, jac->as_print);
2012:     PetscCallExternal(HYPRE_AMSSetMaxIter, jac->hsolver, jac->as_max_iter);
2013:     PetscCallExternal(HYPRE_AMSSetCycleType, jac->hsolver, jac->ams_cycle_type);
2014:     PetscCallExternal(HYPRE_AMSSetTol, jac->hsolver, jac->as_tol);
2015:     PetscCallExternal(HYPRE_AMSSetSmoothingOptions, jac->hsolver, jac->as_relax_type, jac->as_relax_times, jac->as_relax_weight, jac->as_omega);
2016:     PetscCallExternal(HYPRE_AMSSetAlphaAMGOptions, jac->hsolver, jac->as_amg_alpha_opts[0], /* AMG coarsen type */
2017:                       jac->as_amg_alpha_opts[1],                                            /* AMG agg_levels */
2018:                       jac->as_amg_alpha_opts[2],                                            /* AMG relax_type */
2019:                       jac->as_amg_alpha_theta, jac->as_amg_alpha_opts[3],                   /* AMG interp_type */
2020:                       jac->as_amg_alpha_opts[4]);                                           /* AMG Pmax */
2021:     PetscCallExternal(HYPRE_AMSSetBetaAMGOptions, jac->hsolver, jac->as_amg_beta_opts[0],   /* AMG coarsen type */
2022:                       jac->as_amg_beta_opts[1],                                             /* AMG agg_levels */
2023:                       jac->as_amg_beta_opts[2],                                             /* AMG relax_type */
2024:                       jac->as_amg_beta_theta, jac->as_amg_beta_opts[3],                     /* AMG interp_type */
2025:                       jac->as_amg_beta_opts[4]);                                            /* AMG Pmax */
2026:     /* Zero conductivity */
2027:     jac->ams_beta_is_zero      = PETSC_FALSE;
2028:     jac->ams_beta_is_zero_part = PETSC_FALSE;
2029:     PetscFunctionReturn(PETSC_SUCCESS);
2030:   }
2031:   PetscCall(PetscStrcmp("ads", jac->hypre_type, &flag));
2032:   if (flag) {
2033:     PetscCallExternal(HYPRE_ADSCreate, &jac->hsolver);
2034:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ADS;
2035:     pc->ops->view           = PCView_HYPRE_ADS;
2036:     jac->destroy            = HYPRE_ADSDestroy;
2037:     jac->setup              = HYPRE_ADSSetup;
2038:     jac->solve              = HYPRE_ADSSolve;
2039:     jac->coords[0]          = NULL;
2040:     jac->coords[1]          = NULL;
2041:     jac->coords[2]          = NULL;
2042:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */
2043:     jac->as_print       = 0;
2044:     jac->as_max_iter    = 1;  /* used as a preconditioner */
2045:     jac->as_tol         = 0.; /* used as a preconditioner */
2046:     jac->ads_cycle_type = 13;
2047:     /* Smoothing options */
2048:     jac->as_relax_type   = 2;
2049:     jac->as_relax_times  = 1;
2050:     jac->as_relax_weight = 1.0;
2051:     jac->as_omega        = 1.0;
2052:     /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */
2053:     jac->ams_cycle_type       = 14;
2054:     jac->as_amg_alpha_opts[0] = 10;
2055:     jac->as_amg_alpha_opts[1] = 1;
2056:     jac->as_amg_alpha_opts[2] = 6;
2057:     jac->as_amg_alpha_opts[3] = 6;
2058:     jac->as_amg_alpha_opts[4] = 4;
2059:     jac->as_amg_alpha_theta   = 0.25;
2060:     /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
2061:     jac->as_amg_beta_opts[0] = 10;
2062:     jac->as_amg_beta_opts[1] = 1;
2063:     jac->as_amg_beta_opts[2] = 6;
2064:     jac->as_amg_beta_opts[3] = 6;
2065:     jac->as_amg_beta_opts[4] = 4;
2066:     jac->as_amg_beta_theta   = 0.25;
2067:     PetscCallExternal(HYPRE_ADSSetPrintLevel, jac->hsolver, jac->as_print);
2068:     PetscCallExternal(HYPRE_ADSSetMaxIter, jac->hsolver, jac->as_max_iter);
2069:     PetscCallExternal(HYPRE_ADSSetCycleType, jac->hsolver, jac->ams_cycle_type);
2070:     PetscCallExternal(HYPRE_ADSSetTol, jac->hsolver, jac->as_tol);
2071:     PetscCallExternal(HYPRE_ADSSetSmoothingOptions, jac->hsolver, jac->as_relax_type, jac->as_relax_times, jac->as_relax_weight, jac->as_omega);
2072:     PetscCallExternal(HYPRE_ADSSetAMSOptions, jac->hsolver, jac->ams_cycle_type,      /* AMG coarsen type */
2073:                       jac->as_amg_alpha_opts[0],                                      /* AMG coarsen type */
2074:                       jac->as_amg_alpha_opts[1],                                      /* AMG agg_levels */
2075:                       jac->as_amg_alpha_opts[2],                                      /* AMG relax_type */
2076:                       jac->as_amg_alpha_theta, jac->as_amg_alpha_opts[3],             /* AMG interp_type */
2077:                       jac->as_amg_alpha_opts[4]);                                     /* AMG Pmax */
2078:     PetscCallExternal(HYPRE_ADSSetAMGOptions, jac->hsolver, jac->as_amg_beta_opts[0], /* AMG coarsen type */
2079:                       jac->as_amg_beta_opts[1],                                       /* AMG agg_levels */
2080:                       jac->as_amg_beta_opts[2],                                       /* AMG relax_type */
2081:                       jac->as_amg_beta_theta, jac->as_amg_beta_opts[3],               /* AMG interp_type */
2082:                       jac->as_amg_beta_opts[4]);                                      /* AMG Pmax */
2083:     PetscFunctionReturn(PETSC_SUCCESS);
2084:   }
2085:   PetscCall(PetscFree(jac->hypre_type));

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

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

2101:   PetscFunctionBegin;
2102:   PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE preconditioner options");
2103:   PetscCall(PetscOptionsEList("-pc_hypre_type", "HYPRE preconditioner type", "PCHYPRESetType", type, PETSC_STATIC_ARRAY_LENGTH(type), "boomeramg", &indx, &flg));
2104:   if (flg) {
2105:     PetscCall(PCHYPRESetType_HYPRE(pc, type[indx]));
2106:   } else {
2107:     PetscCall(PCHYPRESetType_HYPRE(pc, "boomeramg"));
2108:   }
2109:   PetscTryTypeMethod(pc, setfromoptions, PetscOptionsObject);
2110:   PetscOptionsHeadEnd();
2111:   PetscFunctionReturn(PETSC_SUCCESS);
2112: }

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

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

2121:    Options Database Key:
2122:    -pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads

2124:    Level: intermediate

2126: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHYPRE`
2127: @*/
2128: PetscErrorCode PCHYPRESetType(PC pc, const char name[])
2129: {
2130:   PetscFunctionBegin;
2133:   PetscTryMethod(pc, "PCHYPRESetType_C", (PC, const char[]), (pc, name));
2134:   PetscFunctionReturn(PETSC_SUCCESS);
2135: }

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

2140:    Input Parameter:
2141: .     pc - the preconditioner context

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

2146:    Level: intermediate

2148: .seealso: `PCCreate()`, `PCHYPRESetType()`, `PCType`, `PC`, `PCHYPRE`
2149: @*/
2150: PetscErrorCode PCHYPREGetType(PC pc, const char *name[])
2151: {
2152:   PetscFunctionBegin;
2155:   PetscTryMethod(pc, "PCHYPREGetType_C", (PC, const char *[]), (pc, name));
2156:   PetscFunctionReturn(PETSC_SUCCESS);
2157: }

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

2162:    Logically Collective

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

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

2171:    Level: intermediate

2173:    Developer Note:
2174:    How the name starts with `PCMG`, should it not be `PCHYPREBoomerAMG`?

2176: .seealso: `PCHYPRE`, `PCMGGalerkinGetMatProductAlgorithm()`
2177: @*/
2178: PetscErrorCode PCMGGalerkinSetMatProductAlgorithm(PC pc, const char name[])
2179: {
2180:   PetscFunctionBegin;
2182:   PetscTryMethod(pc, "PCMGGalerkinSetMatProductAlgorithm_C", (PC, const char[]), (pc, name));
2183:   PetscFunctionReturn(PETSC_SUCCESS);
2184: }

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

2189:    Not Collective

2191:    Input Parameter:
2192: .  pc - the multigrid context

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

2197:    Level: intermediate

2199: .seealso: `PCHYPRE`, ``PCMGGalerkinSetMatProductAlgorithm()`
2200: @*/
2201: PetscErrorCode PCMGGalerkinGetMatProductAlgorithm(PC pc, const char *name[])
2202: {
2203:   PetscFunctionBegin;
2205:   PetscTryMethod(pc, "PCMGGalerkinGetMatProductAlgorithm_C", (PC, const char *[]), (pc, name));
2206:   PetscFunctionReturn(PETSC_SUCCESS);
2207: }

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

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

2218:    Level: intermediate

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

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

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

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

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

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

2246:   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
2247:   since it is usually slower, one should use a `KSPType` of `KSPRICHARDSON`
2248:   (or equivalently `-ksp_type richardson`) to achieve this. Using `KSPPREONLY` will not work since it only applies a single cycle of multigrid.

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

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

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

2259: .seealso: `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHYPRESetType()`, `PCPFMG`, `PCGAMG`, `PCSYSPFMG`, `PCSMG`, `PCHYPRESetDiscreteGradient()`,
2260:           `PCHYPRESetDiscreteCurl()`, `PCHYPRESetInterpolations()`, `PCHYPRESetAlphaPoissonMatrix()`, `PCHYPRESetBetaPoissonMatrix()`, `PCHYPRESetEdgeConstantVectors()`,
2261:           PCHYPREAMSSetInteriorNodes()
2262: M*/

2264: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
2265: {
2266:   PC_HYPRE *jac;

2268:   PetscFunctionBegin;
2269:   PetscCall(PetscNew(&jac));

2271:   pc->data                = jac;
2272:   pc->ops->reset          = PCReset_HYPRE;
2273:   pc->ops->destroy        = PCDestroy_HYPRE;
2274:   pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
2275:   pc->ops->setup          = PCSetUp_HYPRE;
2276:   pc->ops->apply          = PCApply_HYPRE;
2277:   jac->comm_hypre         = MPI_COMM_NULL;
2278:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetType_C", PCHYPRESetType_HYPRE));
2279:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetType_C", PCHYPREGetType_HYPRE));
2280:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_HYPRE));
2281:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteGradient_C", PCHYPRESetDiscreteGradient_HYPRE));
2282:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteCurl_C", PCHYPRESetDiscreteCurl_HYPRE));
2283:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetInterpolations_C", PCHYPRESetInterpolations_HYPRE));
2284:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetEdgeConstantVectors_C", PCHYPRESetEdgeConstantVectors_HYPRE));
2285:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREAMSSetInteriorNodes_C", PCHYPREAMSSetInteriorNodes_HYPRE));
2286:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetPoissonMatrix_C", PCHYPRESetPoissonMatrix_HYPRE));
2287:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinSetMatProductAlgorithm_C", PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG));
2288:   PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinGetMatProductAlgorithm_C", PCMGGalerkinGetMatProductAlgorithm_HYPRE_BoomerAMG));
2289: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2290:   #if defined(HYPRE_USING_HIP)
2291:   PetscCall(PetscDeviceInitialize(PETSC_DEVICE_HIP));
2292:   #endif
2293:   #if defined(HYPRE_USING_CUDA)
2294:   PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA));
2295:   #endif
2296: #endif
2297:   PetscFunctionReturn(PETSC_SUCCESS);
2298: }

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

2304:   /* keep copy of PFMG options used so may view them */
2305:   PetscInt  its;
2306:   double    tol;
2307:   PetscInt  relax_type;
2308:   PetscInt  rap_type;
2309:   PetscInt  num_pre_relax, num_post_relax;
2310:   PetscInt  max_levels;
2311:   PetscInt  skip_relax;
2312:   PetscBool print_statistics;
2313: } PC_PFMG;

2315: PetscErrorCode PCDestroy_PFMG(PC pc)
2316: {
2317:   PC_PFMG *ex = (PC_PFMG *)pc->data;

2319:   PetscFunctionBegin;
2320:   if (ex->hsolver) PetscCallExternal(HYPRE_StructPFMGDestroy, ex->hsolver);
2321:   PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2322:   PetscCall(PetscFree(pc->data));
2323:   PetscFunctionReturn(PETSC_SUCCESS);
2324: }

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

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

2334:   PetscFunctionBegin;
2335:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
2336:   if (iascii) {
2337:     PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE PFMG preconditioning\n"));
2338:     PetscCall(PetscViewerASCIIPrintf(viewer, "    max iterations %" PetscInt_FMT "\n", ex->its));
2339:     PetscCall(PetscViewerASCIIPrintf(viewer, "    tolerance %g\n", ex->tol));
2340:     PetscCall(PetscViewerASCIIPrintf(viewer, "    relax type %s\n", PFMGRelaxType[ex->relax_type]));
2341:     PetscCall(PetscViewerASCIIPrintf(viewer, "    RAP type %s\n", PFMGRAPType[ex->rap_type]));
2342:     PetscCall(PetscViewerASCIIPrintf(viewer, "    number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax));
2343:     PetscCall(PetscViewerASCIIPrintf(viewer, "    max levels %" PetscInt_FMT "\n", ex->max_levels));
2344:     PetscCall(PetscViewerASCIIPrintf(viewer, "    skip relax %" PetscInt_FMT "\n", ex->skip_relax));
2345:   }
2346:   PetscFunctionReturn(PETSC_SUCCESS);
2347: }

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

2353:   PetscFunctionBegin;
2354:   PetscOptionsHeadBegin(PetscOptionsObject, "PFMG options");
2355:   PetscCall(PetscOptionsBool("-pc_pfmg_print_statistics", "Print statistics", "HYPRE_StructPFMGSetPrintLevel", ex->print_statistics, &ex->print_statistics, NULL));
2356:   PetscCall(PetscOptionsInt("-pc_pfmg_its", "Number of iterations of PFMG to use as preconditioner", "HYPRE_StructPFMGSetMaxIter", ex->its, &ex->its, NULL));
2357:   PetscCallExternal(HYPRE_StructPFMGSetMaxIter, ex->hsolver, ex->its);
2358:   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));
2359:   PetscCallExternal(HYPRE_StructPFMGSetNumPreRelax, ex->hsolver, ex->num_pre_relax);
2360:   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));
2361:   PetscCallExternal(HYPRE_StructPFMGSetNumPostRelax, ex->hsolver, ex->num_post_relax);

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

2366:   PetscCall(PetscOptionsReal("-pc_pfmg_tol", "Tolerance of PFMG", "HYPRE_StructPFMGSetTol", ex->tol, &ex->tol, NULL));
2367:   PetscCallExternal(HYPRE_StructPFMGSetTol, ex->hsolver, ex->tol);
2368:   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));
2369:   PetscCallExternal(HYPRE_StructPFMGSetRelaxType, ex->hsolver, ex->relax_type);
2370:   PetscCall(PetscOptionsEList("-pc_pfmg_rap_type", "RAP type", "HYPRE_StructPFMGSetRAPType", PFMGRAPType, PETSC_STATIC_ARRAY_LENGTH(PFMGRAPType), PFMGRAPType[ex->rap_type], &ex->rap_type, NULL));
2371:   PetscCallExternal(HYPRE_StructPFMGSetRAPType, ex->hsolver, ex->rap_type);
2372:   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));
2373:   PetscCallExternal(HYPRE_StructPFMGSetSkipRelax, ex->hsolver, ex->skip_relax);
2374:   PetscOptionsHeadEnd();
2375:   PetscFunctionReturn(PETSC_SUCCESS);
2376: }

2378: PetscErrorCode PCApply_PFMG(PC pc, Vec x, Vec y)
2379: {
2380:   PC_PFMG           *ex = (PC_PFMG *)pc->data;
2381:   PetscScalar       *yy;
2382:   const PetscScalar *xx;
2383:   PetscInt           ilower[3], iupper[3];
2384:   HYPRE_Int          hlower[3], hupper[3];
2385:   Mat_HYPREStruct   *mx = (Mat_HYPREStruct *)(pc->pmat->data);

2387:   PetscFunctionBegin;
2388:   PetscCall(PetscCitationsRegister(hypreCitation, &cite));
2389:   PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
2390:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2391:   iupper[0] += ilower[0] - 1;
2392:   iupper[1] += ilower[1] - 1;
2393:   iupper[2] += ilower[2] - 1;
2394:   hlower[0] = (HYPRE_Int)ilower[0];
2395:   hlower[1] = (HYPRE_Int)ilower[1];
2396:   hlower[2] = (HYPRE_Int)ilower[2];
2397:   hupper[0] = (HYPRE_Int)iupper[0];
2398:   hupper[1] = (HYPRE_Int)iupper[1];
2399:   hupper[2] = (HYPRE_Int)iupper[2];

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

2409:   /* copy solution values back to PETSc */
2410:   PetscCall(VecGetArray(y, &yy));
2411:   PetscCallExternal(HYPRE_StructVectorGetBoxValues, mx->hx, hlower, hupper, (HYPRE_Complex *)yy);
2412:   PetscCall(VecRestoreArray(y, &yy));
2413:   PetscFunctionReturn(PETSC_SUCCESS);
2414: }

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

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

2426:   PetscCall(PCApply_PFMG(pc, b, y));
2427:   PetscCallExternal(HYPRE_StructPFMGGetNumIterations, jac->hsolver, &oits);
2428:   *outits = oits;
2429:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2430:   else *reason = PCRICHARDSON_CONVERGED_RTOL;
2431:   PetscCallExternal(HYPRE_StructPFMGSetTol, jac->hsolver, jac->tol);
2432:   PetscCallExternal(HYPRE_StructPFMGSetMaxIter, jac->hsolver, jac->its);
2433:   PetscFunctionReturn(PETSC_SUCCESS);
2434: }

2436: PetscErrorCode PCSetUp_PFMG(PC pc)
2437: {
2438:   PC_PFMG         *ex = (PC_PFMG *)pc->data;
2439:   Mat_HYPREStruct *mx = (Mat_HYPREStruct *)(pc->pmat->data);
2440:   PetscBool        flg;

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

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

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

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

2462:   PetscCallExternal(HYPRE_StructPFMGSetup, ex->hsolver, mx->hmat, mx->hb, mx->hx);
2463:   PetscCallExternal(HYPRE_StructPFMGSetZeroGuess, ex->hsolver);
2464:   PetscFunctionReturn(PETSC_SUCCESS);
2465: }

2467: /*MC
2468:      PCPFMG - the hypre PFMG multigrid solver

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

2480:    Level: advanced

2482:    Notes:
2483:    This is for CELL-centered descretizations

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

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

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

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

2493: .seealso: `PCMG`, `MATHYPRESTRUCT`, `PCHYPRE`, `PCGAMG`, `PCSYSPFMG`, `PCSMG`
2494: M*/

2496: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
2497: {
2498:   PC_PFMG *ex;

2500:   PetscFunctionBegin;
2501:   PetscCall(PetscNew(&ex));
2502:   pc->data = ex;

2504:   ex->its              = 1;
2505:   ex->tol              = 1.e-8;
2506:   ex->relax_type       = 1;
2507:   ex->rap_type         = 0;
2508:   ex->num_pre_relax    = 1;
2509:   ex->num_post_relax   = 1;
2510:   ex->max_levels       = 0;
2511:   ex->skip_relax       = 0;
2512:   ex->print_statistics = PETSC_FALSE;

2514:   pc->ops->setfromoptions  = PCSetFromOptions_PFMG;
2515:   pc->ops->view            = PCView_PFMG;
2516:   pc->ops->destroy         = PCDestroy_PFMG;
2517:   pc->ops->apply           = PCApply_PFMG;
2518:   pc->ops->applyrichardson = PCApplyRichardson_PFMG;
2519:   pc->ops->setup           = PCSetUp_PFMG;

2521:   PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2522:   PetscCallExternal(HYPRE_StructPFMGCreate, ex->hcomm, &ex->hsolver);
2523:   PetscFunctionReturn(PETSC_SUCCESS);
2524: }

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

2531:   /* keep copy of SYSPFMG options used so may view them */
2532:   PetscInt its;
2533:   double   tol;
2534:   PetscInt relax_type;
2535:   PetscInt num_pre_relax, num_post_relax;
2536: } PC_SysPFMG;

2538: PetscErrorCode PCDestroy_SysPFMG(PC pc)
2539: {
2540:   PC_SysPFMG *ex = (PC_SysPFMG *)pc->data;

2542:   PetscFunctionBegin;
2543:   if (ex->ss_solver) PetscCallExternal(HYPRE_SStructSysPFMGDestroy, ex->ss_solver);
2544:   PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2545:   PetscCall(PetscFree(pc->data));
2546:   PetscFunctionReturn(PETSC_SUCCESS);
2547: }

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

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

2556:   PetscFunctionBegin;
2557:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
2558:   if (iascii) {
2559:     PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE SysPFMG preconditioning\n"));
2560:     PetscCall(PetscViewerASCIIPrintf(viewer, "  max iterations %" PetscInt_FMT "\n", ex->its));
2561:     PetscCall(PetscViewerASCIIPrintf(viewer, "  tolerance %g\n", ex->tol));
2562:     PetscCall(PetscViewerASCIIPrintf(viewer, "  relax type %s\n", PFMGRelaxType[ex->relax_type]));
2563:     PetscCall(PetscViewerASCIIPrintf(viewer, "  number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax));
2564:   }
2565:   PetscFunctionReturn(PETSC_SUCCESS);
2566: }

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

2573:   PetscFunctionBegin;
2574:   PetscOptionsHeadBegin(PetscOptionsObject, "SysPFMG options");
2575:   PetscCall(PetscOptionsBool("-pc_syspfmg_print_statistics", "Print statistics", "HYPRE_SStructSysPFMGSetPrintLevel", flg, &flg, NULL));
2576:   if (flg) PetscCallExternal(HYPRE_SStructSysPFMGSetPrintLevel, ex->ss_solver, 3);
2577:   PetscCall(PetscOptionsInt("-pc_syspfmg_its", "Number of iterations of SysPFMG to use as preconditioner", "HYPRE_SStructSysPFMGSetMaxIter", ex->its, &ex->its, NULL));
2578:   PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter, ex->ss_solver, ex->its);
2579:   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));
2580:   PetscCallExternal(HYPRE_SStructSysPFMGSetNumPreRelax, ex->ss_solver, ex->num_pre_relax);
2581:   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));
2582:   PetscCallExternal(HYPRE_SStructSysPFMGSetNumPostRelax, ex->ss_solver, ex->num_post_relax);

2584:   PetscCall(PetscOptionsReal("-pc_syspfmg_tol", "Tolerance of SysPFMG", "HYPRE_SStructSysPFMGSetTol", ex->tol, &ex->tol, NULL));
2585:   PetscCallExternal(HYPRE_SStructSysPFMGSetTol, ex->ss_solver, ex->tol);
2586:   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));
2587:   PetscCallExternal(HYPRE_SStructSysPFMGSetRelaxType, ex->ss_solver, ex->relax_type);
2588:   PetscOptionsHeadEnd();
2589:   PetscFunctionReturn(PETSC_SUCCESS);
2590: }

2592: PetscErrorCode PCApply_SysPFMG(PC pc, Vec x, Vec y)
2593: {
2594:   PC_SysPFMG        *ex = (PC_SysPFMG *)pc->data;
2595:   PetscScalar       *yy;
2596:   const PetscScalar *xx;
2597:   PetscInt           ilower[3], iupper[3];
2598:   HYPRE_Int          hlower[3], hupper[3];
2599:   Mat_HYPRESStruct  *mx       = (Mat_HYPRESStruct *)(pc->pmat->data);
2600:   PetscInt           ordering = mx->dofs_order;
2601:   PetscInt           nvars    = mx->nvars;
2602:   PetscInt           part     = 0;
2603:   PetscInt           size;
2604:   PetscInt           i;

2606:   PetscFunctionBegin;
2607:   PetscCall(PetscCitationsRegister(hypreCitation, &cite));
2608:   PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
2609:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2610:   iupper[0] += ilower[0] - 1;
2611:   iupper[1] += ilower[1] - 1;
2612:   iupper[2] += ilower[2] - 1;
2613:   hlower[0] = (HYPRE_Int)ilower[0];
2614:   hlower[1] = (HYPRE_Int)ilower[1];
2615:   hlower[2] = (HYPRE_Int)ilower[2];
2616:   hupper[0] = (HYPRE_Int)iupper[0];
2617:   hupper[1] = (HYPRE_Int)iupper[1];
2618:   hupper[2] = (HYPRE_Int)iupper[2];

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

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

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

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

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

2655:     /* copy solution values back to PETSc */
2656:     PetscCall(VecGetArray(y, &yy));
2657:     for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues, mx->ss_x, part, hlower, hupper, i, (HYPRE_Complex *)(z + (size * i)));
2658:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
2659:     for (i = 0; i < size; i++) {
2660:       k = i * nvars;
2661:       for (j = 0; j < nvars; j++) yy[k + j] = z[j * size + i];
2662:     }
2663:     PetscCall(VecRestoreArray(y, &yy));
2664:     PetscCall(PetscFree(z));
2665:   }
2666:   PetscFunctionReturn(PETSC_SUCCESS);
2667: }

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

2674:   PetscFunctionBegin;
2675:   PetscCall(PetscCitationsRegister(hypreCitation, &cite));
2676:   PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter, jac->ss_solver, its * jac->its);
2677:   PetscCallExternal(HYPRE_SStructSysPFMGSetTol, jac->ss_solver, rtol);
2678:   PetscCall(PCApply_SysPFMG(pc, b, y));
2679:   PetscCallExternal(HYPRE_SStructSysPFMGGetNumIterations, jac->ss_solver, &oits);
2680:   *outits = oits;
2681:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2682:   else *reason = PCRICHARDSON_CONVERGED_RTOL;
2683:   PetscCallExternal(HYPRE_SStructSysPFMGSetTol, jac->ss_solver, jac->tol);
2684:   PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter, jac->ss_solver, jac->its);
2685:   PetscFunctionReturn(PETSC_SUCCESS);
2686: }

2688: PetscErrorCode PCSetUp_SysPFMG(PC pc)
2689: {
2690:   PC_SysPFMG       *ex = (PC_SysPFMG *)pc->data;
2691:   Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)(pc->pmat->data);
2692:   PetscBool         flg;

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

2698:   /* create the hypre sstruct solver object and set its information */
2699:   if (ex->ss_solver) PetscCallExternal(HYPRE_SStructSysPFMGDestroy, ex->ss_solver);
2700:   PetscCallExternal(HYPRE_SStructSysPFMGCreate, ex->hcomm, &ex->ss_solver);
2701:   PetscCallExternal(HYPRE_SStructSysPFMGSetZeroGuess, ex->ss_solver);
2702:   PetscCallExternal(HYPRE_SStructSysPFMGSetup, ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x);
2703:   PetscFunctionReturn(PETSC_SUCCESS);
2704: }

2706: /*MC
2707:      PCSYSPFMG - the hypre SysPFMG multigrid solver

2709:    Level: advanced

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

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

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

2723:    This is for CELL-centered descretizations

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

2727:    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`.

2729: .seealso: `PCMG`, `MATHYPRESSTRUCT`, `PCPFMG`, `PCHYPRE`, `PCGAMG`, `PCSMG`
2730: M*/

2732: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
2733: {
2734:   PC_SysPFMG *ex;

2736:   PetscFunctionBegin;
2737:   PetscCall(PetscNew(&ex));
2738:   pc->data = ex;

2740:   ex->its            = 1;
2741:   ex->tol            = 1.e-8;
2742:   ex->relax_type     = 1;
2743:   ex->num_pre_relax  = 1;
2744:   ex->num_post_relax = 1;

2746:   pc->ops->setfromoptions  = PCSetFromOptions_SysPFMG;
2747:   pc->ops->view            = PCView_SysPFMG;
2748:   pc->ops->destroy         = PCDestroy_SysPFMG;
2749:   pc->ops->apply           = PCApply_SysPFMG;
2750:   pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
2751:   pc->ops->setup           = PCSetUp_SysPFMG;

2753:   PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2754:   PetscCallExternal(HYPRE_SStructSysPFMGCreate, ex->hcomm, &ex->ss_solver);
2755:   PetscFunctionReturn(PETSC_SUCCESS);
2756: }

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

2768: PetscErrorCode PCDestroy_SMG(PC pc)
2769: {
2770:   PC_SMG *ex = (PC_SMG *)pc->data;

2772:   PetscFunctionBegin;
2773:   if (ex->hsolver) PetscCallExternal(HYPRE_StructSMGDestroy, ex->hsolver);
2774:   PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2775:   PetscCall(PetscFree(pc->data));
2776:   PetscFunctionReturn(PETSC_SUCCESS);
2777: }

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

2784:   PetscFunctionBegin;
2785:   PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &iascii));
2786:   if (iascii) {
2787:     PetscCall(PetscViewerASCIIPrintf(viewer, "  HYPRE SMG preconditioning\n"));
2788:     PetscCall(PetscViewerASCIIPrintf(viewer, "    max iterations %" PetscInt_FMT "\n", ex->its));
2789:     PetscCall(PetscViewerASCIIPrintf(viewer, "    tolerance %g\n", ex->tol));
2790:     PetscCall(PetscViewerASCIIPrintf(viewer, "    number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax));
2791:   }
2792:   PetscFunctionReturn(PETSC_SUCCESS);
2793: }

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

2799:   PetscFunctionBegin;
2800:   PetscOptionsHeadBegin(PetscOptionsObject, "SMG options");

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

2807:   PetscOptionsHeadEnd();
2808:   PetscFunctionReturn(PETSC_SUCCESS);
2809: }

2811: PetscErrorCode PCApply_SMG(PC pc, Vec x, Vec y)
2812: {
2813:   PC_SMG            *ex = (PC_SMG *)pc->data;
2814:   PetscScalar       *yy;
2815:   const PetscScalar *xx;
2816:   PetscInt           ilower[3], iupper[3];
2817:   HYPRE_Int          hlower[3], hupper[3];
2818:   Mat_HYPREStruct   *mx = (Mat_HYPREStruct *)(pc->pmat->data);

2820:   PetscFunctionBegin;
2821:   PetscCall(PetscCitationsRegister(hypreCitation, &cite));
2822:   PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
2823:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2824:   iupper[0] += ilower[0] - 1;
2825:   iupper[1] += ilower[1] - 1;
2826:   iupper[2] += ilower[2] - 1;
2827:   hlower[0] = (HYPRE_Int)ilower[0];
2828:   hlower[1] = (HYPRE_Int)ilower[1];
2829:   hlower[2] = (HYPRE_Int)ilower[2];
2830:   hupper[0] = (HYPRE_Int)iupper[0];
2831:   hupper[1] = (HYPRE_Int)iupper[1];
2832:   hupper[2] = (HYPRE_Int)iupper[2];

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

2842:   /* copy solution values back to PETSc */
2843:   PetscCall(VecGetArray(y, &yy));
2844:   PetscCallExternal(HYPRE_StructVectorGetBoxValues, mx->hx, hlower, hupper, (HYPRE_Complex *)yy);
2845:   PetscCall(VecRestoreArray(y, &yy));
2846:   PetscFunctionReturn(PETSC_SUCCESS);
2847: }

2849: 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)
2850: {
2851:   PC_SMG   *jac = (PC_SMG *)pc->data;
2852:   HYPRE_Int oits;

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

2859:   PetscCall(PCApply_SMG(pc, b, y));
2860:   PetscCallExternal(HYPRE_StructSMGGetNumIterations, jac->hsolver, &oits);
2861:   *outits = oits;
2862:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2863:   else *reason = PCRICHARDSON_CONVERGED_RTOL;
2864:   PetscCallExternal(HYPRE_StructSMGSetTol, jac->hsolver, jac->tol);
2865:   PetscCallExternal(HYPRE_StructSMGSetMaxIter, jac->hsolver, jac->its);
2866:   PetscFunctionReturn(PETSC_SUCCESS);
2867: }

2869: PetscErrorCode PCSetUp_SMG(PC pc)
2870: {
2871:   PetscInt         i, dim;
2872:   PC_SMG          *ex = (PC_SMG *)pc->data;
2873:   Mat_HYPREStruct *mx = (Mat_HYPREStruct *)(pc->pmat->data);
2874:   PetscBool        flg;
2875:   DMBoundaryType   p[3];
2876:   PetscInt         M[3];

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

2882:   PetscCall(DMDAGetInfo(mx->da, &dim, &M[0], &M[1], &M[2], 0, 0, 0, 0, 0, &p[0], &p[1], &p[2], 0));
2883:   // Check if power of 2 in periodic directions
2884:   for (i = 0; i < dim; i++) {
2885:     if (((M[i] & (M[i] - 1)) != 0) && (p[i] == DM_BOUNDARY_PERIODIC)) {
2886:       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]);
2887:     }
2888:   }

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

2899:   PetscCallExternal(HYPRE_StructSMGSetup, ex->hsolver, mx->hmat, mx->hb, mx->hx);
2900:   PetscCallExternal(HYPRE_StructSMGSetZeroGuess, ex->hsolver);
2901:   PetscFunctionReturn(PETSC_SUCCESS);
2902: }

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

2907:    Level: advanced

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

2915:    Notes:
2916:    This is for CELL-centered descretizations

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

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

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

2924: .seealso:  `PCMG`, `MATHYPRESTRUCT`, `PCPFMG`, `PCSYSPFMG`, `PCHYPRE`, `PCGAMG`
2925: M*/

2927: PETSC_EXTERN PetscErrorCode PCCreate_SMG(PC pc)
2928: {
2929:   PC_SMG *ex;

2931:   PetscFunctionBegin;
2932:   PetscCall(PetscNew(&ex));
2933:   pc->data = ex;

2935:   ex->its            = 1;
2936:   ex->tol            = 1.e-8;
2937:   ex->num_pre_relax  = 1;
2938:   ex->num_post_relax = 1;

2940:   pc->ops->setfromoptions  = PCSetFromOptions_SMG;
2941:   pc->ops->view            = PCView_SMG;
2942:   pc->ops->destroy         = PCDestroy_SMG;
2943:   pc->ops->apply           = PCApply_SMG;
2944:   pc->ops->applyrichardson = PCApplyRichardson_SMG;
2945:   pc->ops->setup           = PCSetUp_SMG;

2947:   PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2948:   PetscCallExternal(HYPRE_StructSMGCreate, ex->hcomm, &ex->hsolver);
2949:   PetscFunctionReturn(PETSC_SUCCESS);
2950: }