Actual source code: ihypre.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 PetscDefined(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: PetscObjectParameterDeclare(PetscInt, coarsentype);
64: PetscInt measuretype;
65: PetscInt smoothtype;
66: PetscInt smoothsweeps;
67: PetscInt smoothnumlevels;
68: PetscInt eu_level; /* Number of levels for ILU(k) in Euclid */
69: PetscReal eu_droptolerance; /* Drop tolerance for ILU(k) in Euclid */
70: PetscInt eu_bj; /* Defines use of Block Jacobi ILU in Euclid */
71: PetscObjectParameterDeclare(PetscInt, relaxtype[3]);
72: PetscReal relaxweight;
73: PetscReal outerrelaxweight;
74: PetscObjectParameterDeclare(PetscInt, relaxorder);
75: PetscReal truncfactor;
76: PetscBool applyrichardson;
77: PetscInt pmax;
78: PetscObjectParameterDeclare(PetscInt, interptype);
79: PetscInt maxc;
80: PetscInt minc;
81: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
82: PetscObjectParameterDeclarePtr(const char, spgemm_type); // this is a global hypre parameter but is closely associated with BoomerAMG
83: #endif
84: /* GPU */
85: PetscObjectParameterDeclare(PetscBool3, keeptranspose);
86: PetscInt rap2;
87: PetscObjectParameterDeclare(PetscInt, mod_rap2);
89: /* AIR */
90: PetscInt Rtype;
91: PetscReal Rstrongthreshold;
92: PetscReal Rfilterthreshold;
93: PetscInt Adroptype;
94: PetscReal Adroptol;
96: PetscInt agg_nl;
97: PetscObjectParameterDeclare(PetscInt, agg_interptype);
98: PetscInt agg_num_paths;
99: PetscBool nodal_relax;
100: PetscInt nodal_relax_levels;
102: PetscInt nodal_coarsening;
103: PetscInt nodal_coarsening_diag;
104: PetscInt vec_interp_variant;
105: PetscInt vec_interp_qmax;
106: PetscBool vec_interp_smooth;
107: PetscInt interp_refine;
109: /* NearNullSpace support */
110: VecHYPRE_IJVector *hmnull;
111: HYPRE_ParVector *phmnull;
112: PetscInt n_hmnull;
113: Vec hmnull_constant;
115: /* options for AS (Auxiliary Space preconditioners) */
116: PetscInt as_print;
117: PetscInt as_max_iter;
118: PetscReal as_tol;
119: PetscInt as_relax_type;
120: PetscInt as_relax_times;
121: PetscReal as_relax_weight;
122: PetscReal as_omega;
123: PetscInt as_amg_alpha_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for vector Poisson (AMS) or Curl problem (ADS) */
124: PetscReal as_amg_alpha_theta; /* AMG strength for vector Poisson (AMS) or Curl problem (ADS) */
125: PetscInt as_amg_beta_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for scalar Poisson (AMS) or vector Poisson (ADS) */
126: PetscReal as_amg_beta_theta; /* AMG strength for scalar Poisson (AMS) or vector Poisson (ADS) */
127: PetscInt ams_cycle_type;
128: PetscInt ads_cycle_type;
130: /* additional data */
131: Mat G; /* MatHYPRE */
132: Mat C; /* MatHYPRE */
133: Mat alpha_Poisson; /* MatHYPRE */
134: Mat beta_Poisson; /* MatHYPRE */
136: /* extra information for AMS */
137: PetscInt dim; /* geometrical dimension */
138: VecHYPRE_IJVector coords[3];
139: VecHYPRE_IJVector constants[3];
140: VecHYPRE_IJVector interior;
141: Mat RT_PiFull, RT_Pi[3];
142: Mat ND_PiFull, ND_Pi[3];
143: PetscBool ams_beta_is_zero;
144: PetscBool ams_beta_is_zero_part;
145: PetscInt ams_proj_freq;
146: } PC_HYPRE;
148: /*
149: Matrices with AIJ format are created IN PLACE with using (I,J,data) from BoomerAMG. Since the data format in hypre_ParCSRMatrix
150: is different from that used in PETSc, the original hypre_ParCSRMatrix can not be used any more after call this routine.
151: It is used in PCHMG. Other users should avoid using this function.
152: */
153: static PetscErrorCode PCGetCoarseOperators_BoomerAMG(PC pc, PetscInt *nlevels, Mat *operators[])
154: {
155: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
156: PetscBool same;
157: PetscInt num_levels;
158: Mat *mattmp;
159: hypre_ParCSRMatrix **A_array;
161: PetscFunctionBegin;
162: PetscCall(PetscStrcmp(jac->hypre_type, "boomeramg", &same));
163: PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_NOTSAMETYPE, "Hypre type is not BoomerAMG");
164: num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData *)jac->hsolver);
165: PetscCall(PetscMalloc1(num_levels, &mattmp));
166: A_array = hypre_ParAMGDataAArray((hypre_ParAMGData *)jac->hsolver);
167: for (PetscInt l = 1; l < num_levels; l++) {
168: PetscCall(MatCreateFromParCSR(A_array[l], MATAIJ, PETSC_OWN_POINTER, &mattmp[num_levels - 1 - l]));
169: /* We want to own the data, and HYPRE can not touch this matrix any more */
170: A_array[l] = NULL;
171: }
172: *nlevels = num_levels;
173: *operators = mattmp;
174: PetscFunctionReturn(PETSC_SUCCESS);
175: }
177: /*
178: Matrices with AIJ format are created IN PLACE with using (I,J,data) from BoomerAMG. Since the data format in hypre_ParCSRMatrix
179: is different from that used in PETSc, the original hypre_ParCSRMatrix can not be used any more after call this routine.
180: It is used in PCHMG. Other users should avoid using this function.
181: */
182: static PetscErrorCode PCGetInterpolations_BoomerAMG(PC pc, PetscInt *nlevels, Mat *interpolations[])
183: {
184: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
185: PetscBool same;
186: PetscInt num_levels;
187: Mat *mattmp;
188: hypre_ParCSRMatrix **P_array;
190: PetscFunctionBegin;
191: PetscCall(PetscStrcmp(jac->hypre_type, "boomeramg", &same));
192: PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_NOTSAMETYPE, "Hypre type is not BoomerAMG");
193: num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData *)jac->hsolver);
194: PetscCall(PetscMalloc1(num_levels, &mattmp));
195: P_array = hypre_ParAMGDataPArray((hypre_ParAMGData *)jac->hsolver);
196: for (PetscInt l = 1; l < num_levels; l++) {
197: PetscCall(MatCreateFromParCSR(P_array[num_levels - 1 - l], MATAIJ, PETSC_OWN_POINTER, &mattmp[l - 1]));
198: /* We want to own the data, and HYPRE can not touch this matrix any more */
199: P_array[num_levels - 1 - l] = NULL;
200: }
201: *nlevels = num_levels;
202: *interpolations = mattmp;
203: PetscFunctionReturn(PETSC_SUCCESS);
204: }
206: /*
207: Boolean Vecs are created IN PLACE with using data from BoomerAMG.
208: */
209: static PetscErrorCode PCHYPREGetCFMarkers_BoomerAMG(PC pc, PetscInt *n_per_level[], PetscBT *CFMarkers[])
210: {
211: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
212: PetscBool same;
213: PetscInt num_levels, fine_nodes = 0, coarse_nodes;
214: PetscInt *n_per_temp;
215: PetscBT *markertmp;
216: hypre_IntArray **CF_marker_array;
218: PetscFunctionBegin;
219: PetscCall(PetscStrcmp(jac->hypre_type, "boomeramg", &same));
220: PetscCheck(same, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_NOTSAMETYPE, "Hypre type is not BoomerAMG");
221: num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData *)jac->hsolver);
222: PetscCall(PetscMalloc1(num_levels, &n_per_temp));
223: PetscCall(PetscMalloc1(num_levels - 1, &markertmp));
224: CF_marker_array = hypre_ParAMGDataCFMarkerArray((hypre_ParAMGData *)jac->hsolver);
225: for (PetscInt l = 0, CFMaxIndex = num_levels - 2; CFMaxIndex >= 0; l++, CFMaxIndex--) {
226: fine_nodes = hypre_IntArraySize(CF_marker_array[CFMaxIndex]);
227: coarse_nodes = 0;
228: PetscCall(PetscBTCreate(fine_nodes, &markertmp[l]));
229: for (PetscInt k = 0; k < fine_nodes; k++) {
230: if (hypre_IntArrayDataI(CF_marker_array[CFMaxIndex], k) > 0) {
231: PetscCall(PetscBTSet(markertmp[l], k));
232: coarse_nodes++;
233: }
234: }
235: n_per_temp[l] = coarse_nodes;
236: }
237: n_per_temp[num_levels - 1] = fine_nodes;
238: *n_per_level = n_per_temp;
239: *CFMarkers = markertmp;
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: /* Resets (frees) Hypre's representation of the near null space */
244: static PetscErrorCode PCHYPREResetNearNullSpace_Private(PC pc)
245: {
246: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
247: PetscInt i;
249: PetscFunctionBegin;
250: for (i = 0; i < jac->n_hmnull; i++) PetscCall(VecHYPRE_IJVectorDestroy(&jac->hmnull[i]));
251: PetscCall(PetscFree(jac->hmnull));
252: PetscCall(PetscFree(jac->phmnull));
253: PetscCall(VecDestroy(&jac->hmnull_constant));
254: jac->n_hmnull = 0;
255: PetscFunctionReturn(PETSC_SUCCESS);
256: }
258: static const char *HYPRESpgemmTypes[] = {"cusparse", "hypre"};
259: static PetscErrorCode PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG(PC pc, const char name[])
260: {
261: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
263: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
264: PetscFunctionBegin;
265: jac->spgemm_type = name;
266: PetscFunctionReturn(PETSC_SUCCESS);
267: #endif
268: }
270: static PetscErrorCode PCSetUp_HYPRE(PC pc)
271: {
272: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
273: Mat_HYPRE *hjac;
274: HYPRE_ParCSRMatrix hmat;
275: HYPRE_ParVector bv, xv;
276: PetscBool ishypre;
278: PetscFunctionBegin;
279: /* default type is boomerAMG */
280: if (!jac->hypre_type) PetscCall(PCHYPRESetType(pc, "boomeramg"));
282: /* get hypre matrix */
283: if (pc->flag == DIFFERENT_NONZERO_PATTERN) PetscCall(MatDestroy(&jac->hpmat));
284: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRE, &ishypre));
285: if (!ishypre) {
286: #if PetscDefined(HAVE_HYPRE_DEVICE) && PETSC_PKG_HYPRE_VERSION_LE(2, 30, 0)
287: /* Temporary fix since we do not support MAT_REUSE_MATRIX with HYPRE device */
288: PetscBool iscuda, iship, iskokkos;
290: PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &iscuda, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
291: PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &iship, MATSEQAIJHIPSPARSE, MATMPIAIJHIPSPARSE, ""));
292: PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &iskokkos, MATSEQAIJKOKKOS, MATMPIAIJKOKKOS, ""));
293: if (iscuda || iship || iskokkos) PetscCall(MatDestroy(&jac->hpmat));
294: #endif
295: PetscCall(MatConvert(pc->pmat, MATHYPRE, jac->hpmat ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, &jac->hpmat));
296: } else {
297: PetscCall(PetscObjectReference((PetscObject)pc->pmat));
298: PetscCall(MatDestroy(&jac->hpmat));
299: jac->hpmat = pc->pmat;
300: }
302: /* allow debug */
303: PetscCall(MatViewFromOptions(jac->hpmat, NULL, "-pc_hypre_mat_view"));
304: hjac = (Mat_HYPRE *)jac->hpmat->data;
306: /* special case for BoomerAMG */
307: if (jac->setup == HYPRE_BoomerAMGSetup) {
308: MatNullSpace mnull;
309: PetscBool has_const;
310: PetscInt bs, nvec, i;
311: PetscMemType memtype;
312: const Vec *vecs;
314: PetscCall(MatGetCurrentMemType(jac->hpmat, &memtype));
315: if (PetscMemTypeDevice(memtype)) {
316: /* GPU defaults
317: From https://hypre.readthedocs.io/en/latest/solvers-boomeramg.html#gpu-supported-options
318: and /src/parcsr_ls/par_amg.c
319: First handle options which users have interfaces for changing */
320: PetscObjectParameterSetDefault(jac, coarsentype, 8);
321: PetscObjectParameterSetDefault(jac, relaxorder, 0);
322: PetscObjectParameterSetDefault(jac, interptype, 6);
323: PetscObjectParameterSetDefault(jac, relaxtype[0], 18);
324: PetscObjectParameterSetDefault(jac, relaxtype[1], 18);
325: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
326: PetscObjectParameterSetDefault(jac, spgemm_type, HYPRESpgemmTypes[0]);
327: #endif
328: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
329: PetscObjectParameterSetDefault(jac, keeptranspose, PETSC_BOOL3_TRUE);
330: PetscObjectParameterSetDefault(jac, mod_rap2, 1);
331: #endif
332: PetscObjectParameterSetDefault(jac, agg_interptype, 7);
333: } else {
334: PetscObjectParameterSetDefault(jac, coarsentype, 6);
335: PetscObjectParameterSetDefault(jac, relaxorder, 1);
336: PetscObjectParameterSetDefault(jac, interptype, 0);
337: PetscObjectParameterSetDefault(jac, relaxtype[0], 6);
338: PetscObjectParameterSetDefault(jac, relaxtype[1], 6); /* Defaults to SYMMETRIC since in PETSc we are using a PC - most likely with CG */
339: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
340: PetscObjectParameterSetDefault(jac, spgemm_type, "hypre");
341: #endif
342: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
343: PetscObjectParameterSetDefault(jac, keeptranspose, PETSC_BOOL3_FALSE);
344: PetscObjectParameterSetDefault(jac, mod_rap2, 0);
345: #endif
346: PetscObjectParameterSetDefault(jac, agg_interptype, 4);
347: }
348: PetscCallHYPRE(HYPRE_BoomerAMGSetCycleType(jac->hsolver, (HYPRE_Int)jac->cycletype));
349: PetscCallHYPRE(HYPRE_BoomerAMGSetMaxLevels(jac->hsolver, (HYPRE_Int)jac->maxlevels));
350: PetscCallHYPRE(HYPRE_BoomerAMGSetMaxIter(jac->hsolver, (HYPRE_Int)jac->maxiter));
351: PetscCallHYPRE(HYPRE_BoomerAMGSetTol(jac->hsolver, jac->tol));
352: PetscCallHYPRE(HYPRE_BoomerAMGSetTruncFactor(jac->hsolver, jac->truncfactor));
353: PetscCallHYPRE(HYPRE_BoomerAMGSetStrongThreshold(jac->hsolver, jac->strongthreshold));
354: PetscCallHYPRE(HYPRE_BoomerAMGSetMaxRowSum(jac->hsolver, jac->maxrowsum));
355: PetscCallHYPRE(HYPRE_BoomerAMGSetMeasureType(jac->hsolver, (HYPRE_Int)jac->measuretype));
356: PetscCallHYPRE(HYPRE_BoomerAMGSetAggNumLevels(jac->hsolver, (HYPRE_Int)jac->agg_nl));
357: PetscCallHYPRE(HYPRE_BoomerAMGSetPMaxElmts(jac->hsolver, (HYPRE_Int)jac->pmax));
358: PetscCallHYPRE(HYPRE_BoomerAMGSetNumPaths(jac->hsolver, (HYPRE_Int)jac->agg_num_paths));
359: PetscCallHYPRE(HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver, (HYPRE_Int)jac->gridsweeps[0], 1));
360: PetscCallHYPRE(HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver, (HYPRE_Int)jac->gridsweeps[1], 2));
361: PetscCallHYPRE(HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver, (HYPRE_Int)jac->gridsweeps[2], 3));
362: PetscCallHYPRE(HYPRE_BoomerAMGSetMaxCoarseSize(jac->hsolver, (HYPRE_Int)jac->maxc));
363: PetscCallHYPRE(HYPRE_BoomerAMGSetMinCoarseSize(jac->hsolver, (HYPRE_Int)jac->minc));
364: PetscCallHYPRE(HYPRE_BoomerAMGSetCoarsenType(jac->hsolver, (HYPRE_Int)jac->coarsentype));
365: PetscCallHYPRE(HYPRE_BoomerAMGSetRelaxOrder(jac->hsolver, (HYPRE_Int)jac->relaxorder));
366: PetscCallHYPRE(HYPRE_BoomerAMGSetInterpType(jac->hsolver, (HYPRE_Int)jac->interptype));
367: PetscCallHYPRE(HYPRE_BoomerAMGSetRelaxType(jac->hsolver, (HYPRE_Int)jac->relaxtype[0]));
368: PetscCallHYPRE(HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, (HYPRE_Int)jac->relaxtype[0], 1));
369: PetscCallHYPRE(HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, (HYPRE_Int)jac->relaxtype[1], 2));
370: PetscCallHYPRE(HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, (HYPRE_Int)jac->relaxtype[2], 3));
371: /* GPU */
372: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
373: {
374: PetscBool flg_cusparse, flg_hypre;
376: PetscCall(PetscStrcmp("cusparse", jac->spgemm_type, &flg_cusparse));
377: PetscCall(PetscStrcmp("hypre", jac->spgemm_type, &flg_hypre));
378: if (flg_cusparse) PetscCallHYPRE(HYPRE_SetSpGemmUseCusparse(1));
379: else if (flg_hypre) PetscCallHYPRE(HYPRE_SetSpGemmUseCusparse(0));
380: else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown HYPRE SpGEMM type %s; Choices are cusparse, hypre", jac->spgemm_type);
381: }
382: #endif
383: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
384: PetscCallHYPRE(HYPRE_BoomerAMGSetKeepTranspose(jac->hsolver, jac->keeptranspose == PETSC_BOOL3_TRUE ? 1 : 0));
385: PetscCallHYPRE(HYPRE_BoomerAMGSetRAP2(jac->hsolver, (HYPRE_Int)jac->rap2));
386: PetscCallHYPRE(HYPRE_BoomerAMGSetModuleRAP2(jac->hsolver, (HYPRE_Int)jac->mod_rap2));
387: #endif
388: PetscCallHYPRE(HYPRE_BoomerAMGSetAggInterpType(jac->hsolver, (HYPRE_Int)jac->agg_interptype));
390: /* AIR */
391: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
392: PetscCallHYPRE(HYPRE_BoomerAMGSetRestriction(jac->hsolver, (HYPRE_Int)jac->Rtype));
393: PetscCallHYPRE(HYPRE_BoomerAMGSetStrongThresholdR(jac->hsolver, jac->Rstrongthreshold));
394: PetscCallHYPRE(HYPRE_BoomerAMGSetFilterThresholdR(jac->hsolver, jac->Rfilterthreshold));
395: PetscCallHYPRE(HYPRE_BoomerAMGSetADropTol(jac->hsolver, jac->Adroptol));
396: PetscCallHYPRE(HYPRE_BoomerAMGSetADropType(jac->hsolver, (HYPRE_Int)jac->Adroptype));
397: #endif
399: PetscCall(MatGetBlockSize(pc->pmat, &bs));
400: if (bs > 1) PetscCallHYPRE(HYPRE_BoomerAMGSetNumFunctions(jac->hsolver, (HYPRE_Int)bs));
401: PetscCall(MatGetNearNullSpace(pc->mat, &mnull));
402: if (mnull) {
403: PetscCall(PCHYPREResetNearNullSpace_Private(pc));
404: PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs));
405: PetscCall(PetscMalloc1(nvec + 1, &jac->hmnull));
406: PetscCall(PetscMalloc1(nvec + 1, &jac->phmnull));
407: for (i = 0; i < nvec; i++) {
408: PetscCall(VecHYPRE_IJVectorCreate(vecs[i]->map, &jac->hmnull[i]));
409: PetscCall(VecHYPRE_IJVectorCopy(vecs[i], jac->hmnull[i]));
410: PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->hmnull[i]->ij, (void **)&jac->phmnull[i]));
411: }
412: if (has_const) {
413: PetscCall(MatCreateVecs(pc->pmat, &jac->hmnull_constant, NULL));
414: PetscCall(VecSet(jac->hmnull_constant, 1));
415: PetscCall(VecNormalize(jac->hmnull_constant, NULL));
416: PetscCall(VecHYPRE_IJVectorCreate(jac->hmnull_constant->map, &jac->hmnull[nvec]));
417: PetscCall(VecHYPRE_IJVectorCopy(jac->hmnull_constant, jac->hmnull[nvec]));
418: PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->hmnull[nvec]->ij, (void **)&jac->phmnull[nvec]));
419: nvec++;
420: }
421: PetscCallHYPRE(HYPRE_BoomerAMGSetInterpVectors(jac->hsolver, (HYPRE_Int)nvec, jac->phmnull));
422: jac->n_hmnull = nvec;
423: }
424: }
426: /* special case for AMS */
427: if (jac->setup == HYPRE_AMSSetup) {
428: Mat_HYPRE *hm;
429: HYPRE_ParCSRMatrix parcsr;
430: PetscCheck(jac->coords[0] || jac->constants[0] || jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]), 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()");
431: if (jac->dim) PetscCallHYPRE(HYPRE_AMSSetDimension(jac->hsolver, (HYPRE_Int)jac->dim));
432: if (jac->constants[0]) {
433: HYPRE_ParVector ozz, zoz, zzo = NULL;
434: PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->constants[0]->ij, (void **)(&ozz)));
435: PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->constants[1]->ij, (void **)(&zoz)));
436: if (jac->constants[2]) PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->constants[2]->ij, (void **)(&zzo)));
437: PetscCallHYPRE(HYPRE_AMSSetEdgeConstantVectors(jac->hsolver, ozz, zoz, zzo));
438: }
439: if (jac->coords[0]) {
440: HYPRE_ParVector coords[3];
441: coords[0] = NULL;
442: coords[1] = NULL;
443: coords[2] = NULL;
444: if (jac->coords[0]) PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->coords[0]->ij, (void **)(&coords[0])));
445: if (jac->coords[1]) PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->coords[1]->ij, (void **)(&coords[1])));
446: if (jac->coords[2]) PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->coords[2]->ij, (void **)(&coords[2])));
447: PetscCallHYPRE(HYPRE_AMSSetCoordinateVectors(jac->hsolver, coords[0], coords[1], coords[2]));
448: }
449: PetscCheck(jac->G, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE AMS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient");
450: hm = (Mat_HYPRE *)jac->G->data;
451: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&parcsr)));
452: PetscCallHYPRE(HYPRE_AMSSetDiscreteGradient(jac->hsolver, parcsr));
453: if (jac->alpha_Poisson) {
454: hm = (Mat_HYPRE *)jac->alpha_Poisson->data;
455: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&parcsr)));
456: PetscCallHYPRE(HYPRE_AMSSetAlphaPoissonMatrix(jac->hsolver, parcsr));
457: }
458: if (jac->ams_beta_is_zero) {
459: PetscCallHYPRE(HYPRE_AMSSetBetaPoissonMatrix(jac->hsolver, NULL));
460: } else if (jac->beta_Poisson) {
461: hm = (Mat_HYPRE *)jac->beta_Poisson->data;
462: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&parcsr)));
463: PetscCallHYPRE(HYPRE_AMSSetBetaPoissonMatrix(jac->hsolver, parcsr));
464: } else if (jac->ams_beta_is_zero_part) {
465: if (jac->interior) {
466: HYPRE_ParVector interior = NULL;
467: PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->interior->ij, (void **)(&interior)));
468: PetscCallHYPRE(HYPRE_AMSSetInteriorNodes(jac->hsolver, interior));
469: } else {
470: jac->ams_beta_is_zero_part = PETSC_FALSE;
471: }
472: }
473: if (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1])) {
474: HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3];
475: if (jac->ND_PiFull) {
476: hm = (Mat_HYPRE *)jac->ND_PiFull->data;
477: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&nd_parcsrfull)));
478: } else {
479: nd_parcsrfull = NULL;
480: }
481: for (PetscInt i = 0; i < 3; ++i) {
482: if (jac->ND_Pi[i]) {
483: hm = (Mat_HYPRE *)jac->ND_Pi[i]->data;
484: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&nd_parcsr[i])));
485: } else {
486: nd_parcsr[i] = NULL;
487: }
488: }
489: PetscCallHYPRE(HYPRE_AMSSetInterpolations(jac->hsolver, nd_parcsrfull, nd_parcsr[0], nd_parcsr[1], nd_parcsr[2]));
490: }
491: }
492: /* special case for ADS */
493: if (jac->setup == HYPRE_ADSSetup) {
494: Mat_HYPRE *hm;
495: HYPRE_ParCSRMatrix parcsr;
496: 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])))) {
497: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE ADS preconditioner needs either the coordinate vectors via PCSetCoordinates() or the interpolation matrices via PCHYPRESetInterpolations");
498: } 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");
499: PetscCheck(jac->G, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE ADS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient");
500: PetscCheck(jac->C, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE ADS preconditioner needs the discrete curl operator via PCHYPRESetDiscreteGradient");
501: if (jac->coords[0]) {
502: HYPRE_ParVector coords[3];
503: coords[0] = NULL;
504: coords[1] = NULL;
505: coords[2] = NULL;
506: if (jac->coords[0]) PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->coords[0]->ij, (void **)(&coords[0])));
507: if (jac->coords[1]) PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->coords[1]->ij, (void **)(&coords[1])));
508: if (jac->coords[2]) PetscCallHYPRE(HYPRE_IJVectorGetObject(jac->coords[2]->ij, (void **)(&coords[2])));
509: PetscCallHYPRE(HYPRE_ADSSetCoordinateVectors(jac->hsolver, coords[0], coords[1], coords[2]));
510: }
511: hm = (Mat_HYPRE *)jac->G->data;
512: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&parcsr)));
513: PetscCallHYPRE(HYPRE_ADSSetDiscreteGradient(jac->hsolver, parcsr));
514: hm = (Mat_HYPRE *)jac->C->data;
515: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&parcsr)));
516: PetscCallHYPRE(HYPRE_ADSSetDiscreteCurl(jac->hsolver, parcsr));
517: if ((jac->RT_PiFull || (jac->RT_Pi[0] && jac->RT_Pi[1])) && (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]))) {
518: HYPRE_ParCSRMatrix rt_parcsrfull, rt_parcsr[3];
519: HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3];
520: if (jac->RT_PiFull) {
521: hm = (Mat_HYPRE *)jac->RT_PiFull->data;
522: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&rt_parcsrfull)));
523: } else {
524: rt_parcsrfull = NULL;
525: }
526: for (PetscInt i = 0; i < 3; ++i) {
527: if (jac->RT_Pi[i]) {
528: hm = (Mat_HYPRE *)jac->RT_Pi[i]->data;
529: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&rt_parcsr[i])));
530: } else {
531: rt_parcsr[i] = NULL;
532: }
533: }
534: if (jac->ND_PiFull) {
535: hm = (Mat_HYPRE *)jac->ND_PiFull->data;
536: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&nd_parcsrfull)));
537: } else {
538: nd_parcsrfull = NULL;
539: }
540: for (PetscInt i = 0; i < 3; ++i) {
541: if (jac->ND_Pi[i]) {
542: hm = (Mat_HYPRE *)jac->ND_Pi[i]->data;
543: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hm->ij, (void **)(&nd_parcsr[i])));
544: } else {
545: nd_parcsr[i] = NULL;
546: }
547: }
548: PetscCallHYPRE(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]));
549: }
550: }
551: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hjac->ij, (void **)&hmat));
552: PetscCallHYPRE(HYPRE_IJVectorGetObject(hjac->b->ij, (void **)&bv));
553: PetscCallHYPRE(HYPRE_IJVectorGetObject(hjac->x->ij, (void **)&xv));
554: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
555: PetscCallHYPRE((*jac->setup)(jac->hsolver, hmat, bv, xv));
556: PetscCall(PetscFPTrapPop());
557: PetscFunctionReturn(PETSC_SUCCESS);
558: }
560: static PetscErrorCode PCApply_HYPRE(PC pc, Vec b, Vec x)
561: {
562: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
563: Mat_HYPRE *hjac = (Mat_HYPRE *)jac->hpmat->data;
564: HYPRE_ParCSRMatrix hmat;
565: HYPRE_ParVector jbv, jxv;
567: PetscFunctionBegin;
568: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
569: if (!jac->applyrichardson) PetscCall(VecSet(x, 0.0));
570: PetscCall(VecHYPRE_IJVectorPushVecRead(hjac->b, b));
571: if (jac->applyrichardson) PetscCall(VecHYPRE_IJVectorPushVec(hjac->x, x));
572: else PetscCall(VecHYPRE_IJVectorPushVecWrite(hjac->x, x));
573: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hjac->ij, (void **)&hmat));
574: PetscCallHYPRE(HYPRE_IJVectorGetObject(hjac->b->ij, (void **)&jbv));
575: PetscCallHYPRE(HYPRE_IJVectorGetObject(hjac->x->ij, (void **)&jxv));
576: PetscCallExternalVoid(
577: "Hypre solve", do {
578: HYPRE_Int hierr = (*jac->solve)(jac->hsolver, hmat, jbv, jxv);
579: if (hierr) {
580: PetscCheck(hierr == HYPRE_ERROR_CONV, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in HYPRE solver, error code %d", (int)hierr);
581: HYPRE_ClearAllErrors();
582: }
583: } while (0));
585: if (jac->setup == HYPRE_AMSSetup && jac->ams_beta_is_zero_part) PetscCallHYPRE(HYPRE_AMSProjectOutGradients(jac->hsolver, jxv));
586: PetscCall(VecHYPRE_IJVectorPopVec(hjac->x));
587: PetscCall(VecHYPRE_IJVectorPopVec(hjac->b));
588: PetscFunctionReturn(PETSC_SUCCESS);
589: }
591: static PetscErrorCode PCMatApply_HYPRE_BoomerAMG(PC pc, Mat B, Mat X)
592: {
593: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
594: Mat_HYPRE *hjac = (Mat_HYPRE *)jac->hpmat->data;
595: hypre_ParCSRMatrix *par_matrix;
596: HYPRE_ParVector hb, hx;
597: const PetscScalar *b;
598: PetscScalar *x;
599: PetscInt m, N, lda;
600: hypre_Vector *x_local;
601: PetscMemType type;
603: PetscFunctionBegin;
604: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
605: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hjac->ij, (void **)&par_matrix));
606: PetscCall(MatGetLocalSize(B, &m, NULL));
607: PetscCall(MatGetSize(B, NULL, &N));
608: PetscCallHYPRE(HYPRE_ParMultiVectorCreate(hypre_ParCSRMatrixComm(par_matrix), hypre_ParCSRMatrixGlobalNumRows(par_matrix), hypre_ParCSRMatrixRowStarts(par_matrix), (HYPRE_Int)N, &hb));
609: PetscCallHYPRE(HYPRE_ParMultiVectorCreate(hypre_ParCSRMatrixComm(par_matrix), hypre_ParCSRMatrixGlobalNumRows(par_matrix), hypre_ParCSRMatrixRowStarts(par_matrix), (HYPRE_Int)N, &hx));
610: PetscCall(MatZeroEntries(X));
611: PetscCall(MatDenseGetArrayReadAndMemType(B, &b, &type));
612: PetscCall(MatDenseGetLDA(B, &lda));
613: PetscCheck(lda == m, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot use a LDA different than the number of local rows: % " PetscInt_FMT " != % " PetscInt_FMT, lda, m);
614: PetscCall(MatDenseGetLDA(X, &lda));
615: PetscCheck(lda == m, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Cannot use a LDA different than the number of local rows: % " PetscInt_FMT " != % " PetscInt_FMT, lda, m);
616: x_local = hypre_ParVectorLocalVector(hb);
617: PetscCallHYPRE(hypre_SeqVectorSetDataOwner(x_local, 0));
618: hypre_VectorData(x_local) = (HYPRE_Complex *)b;
619: PetscCall(MatDenseGetArrayWriteAndMemType(X, &x, NULL));
620: x_local = hypre_ParVectorLocalVector(hx);
621: PetscCallHYPRE(hypre_SeqVectorSetDataOwner(x_local, 0));
622: hypre_VectorData(x_local) = (HYPRE_Complex *)x;
623: PetscCallHYPRE(hypre_ParVectorInitialize_v2(hb, type == PETSC_MEMTYPE_HOST ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE));
624: PetscCallHYPRE(hypre_ParVectorInitialize_v2(hx, type == PETSC_MEMTYPE_HOST ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE));
625: PetscCallExternalVoid(
626: "Hypre solve", do {
627: HYPRE_Int hierr = (*jac->solve)(jac->hsolver, par_matrix, hb, hx);
628: if (hierr) {
629: PetscCheck(hierr == HYPRE_ERROR_CONV, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in HYPRE solver, error code %d", (int)hierr);
630: HYPRE_ClearAllErrors();
631: }
632: } while (0));
633: PetscCallHYPRE(HYPRE_ParVectorDestroy(hb));
634: PetscCallHYPRE(HYPRE_ParVectorDestroy(hx));
635: PetscCall(MatDenseRestoreArrayReadAndMemType(B, &b));
636: PetscCall(MatDenseRestoreArrayWriteAndMemType(X, &x));
637: PetscFunctionReturn(PETSC_SUCCESS);
638: }
640: static PetscErrorCode PCReset_HYPRE(PC pc)
641: {
642: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
644: PetscFunctionBegin;
645: PetscCall(MatDestroy(&jac->hpmat));
646: PetscCall(MatDestroy(&jac->G));
647: PetscCall(MatDestroy(&jac->C));
648: PetscCall(MatDestroy(&jac->alpha_Poisson));
649: PetscCall(MatDestroy(&jac->beta_Poisson));
650: PetscCall(MatDestroy(&jac->RT_PiFull));
651: PetscCall(MatDestroy(&jac->RT_Pi[0]));
652: PetscCall(MatDestroy(&jac->RT_Pi[1]));
653: PetscCall(MatDestroy(&jac->RT_Pi[2]));
654: PetscCall(MatDestroy(&jac->ND_PiFull));
655: PetscCall(MatDestroy(&jac->ND_Pi[0]));
656: PetscCall(MatDestroy(&jac->ND_Pi[1]));
657: PetscCall(MatDestroy(&jac->ND_Pi[2]));
658: PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[0]));
659: PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[1]));
660: PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[2]));
661: PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[0]));
662: PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[1]));
663: PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[2]));
664: PetscCall(VecHYPRE_IJVectorDestroy(&jac->interior));
665: PetscCall(PCHYPREResetNearNullSpace_Private(pc));
666: jac->ams_beta_is_zero = PETSC_FALSE;
667: jac->ams_beta_is_zero_part = PETSC_FALSE;
668: jac->dim = 0;
669: PetscFunctionReturn(PETSC_SUCCESS);
670: }
672: static PetscErrorCode PCDestroy_HYPRE(PC pc)
673: {
674: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
676: PetscFunctionBegin;
677: PetscCall(PCReset_HYPRE(pc));
678: if (jac->destroy) PetscCallHYPRE((*jac->destroy)(jac->hsolver));
679: PetscCall(PetscFree(jac->hypre_type));
680: if (jac->comm_hypre != MPI_COMM_NULL) PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
681: PetscCall(PetscFree(pc->data));
683: PetscCall(PetscObjectChangeTypeName((PetscObject)pc, 0));
684: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetType_C", NULL));
685: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetType_C", NULL));
686: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteGradient_C", NULL));
687: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteCurl_C", NULL));
688: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetInterpolations_C", NULL));
689: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetPoissonMatrix_C", NULL));
690: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetEdgeConstantVectors_C", NULL));
691: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREAMSSetInteriorNodes_C", NULL));
692: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", NULL));
693: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", NULL));
694: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetCFMarkers_C", NULL));
695: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinSetMatProductAlgorithm_C", NULL));
696: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinGetMatProductAlgorithm_C", NULL));
697: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
698: PetscFunctionReturn(PETSC_SUCCESS);
699: }
701: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PC pc, PetscOptionItems PetscOptionsObject)
702: {
703: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
704: PetscBool flag;
706: PetscFunctionBegin;
707: PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE Pilut Options");
708: PetscCall(PetscOptionsInt("-pc_hypre_pilut_maxiter", "Number of iterations", "None", jac->maxiter, &jac->maxiter, &flag));
709: if (flag) PetscCallHYPRE(HYPRE_ParCSRPilutSetMaxIter(jac->hsolver, (HYPRE_Int)jac->maxiter));
710: PetscCall(PetscOptionsReal("-pc_hypre_pilut_tol", "Drop tolerance", "None", jac->tol, &jac->tol, &flag));
711: if (flag) PetscCallHYPRE(HYPRE_ParCSRPilutSetDropTolerance(jac->hsolver, jac->tol));
712: PetscCall(PetscOptionsInt("-pc_hypre_pilut_factorrowsize", "FactorRowSize", "None", jac->factorrowsize, &jac->factorrowsize, &flag));
713: if (flag) PetscCallHYPRE(HYPRE_ParCSRPilutSetFactorRowSize(jac->hsolver, (HYPRE_Int)jac->factorrowsize));
714: PetscOptionsHeadEnd();
715: PetscFunctionReturn(PETSC_SUCCESS);
716: }
718: static PetscErrorCode PCView_HYPRE_Pilut(PC pc, PetscViewer viewer)
719: {
720: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
721: PetscBool isascii;
723: PetscFunctionBegin;
724: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
725: if (isascii) {
726: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE Pilut preconditioning\n"));
727: if (jac->maxiter != PETSC_DEFAULT) {
728: PetscCall(PetscViewerASCIIPrintf(viewer, " maximum number of iterations %" PetscInt_FMT "\n", jac->maxiter));
729: } else {
730: PetscCall(PetscViewerASCIIPrintf(viewer, " default maximum number of iterations \n"));
731: }
732: if (jac->tol != PETSC_DEFAULT) {
733: PetscCall(PetscViewerASCIIPrintf(viewer, " drop tolerance %g\n", (double)jac->tol));
734: } else {
735: PetscCall(PetscViewerASCIIPrintf(viewer, " default drop tolerance \n"));
736: }
737: if (jac->factorrowsize != PETSC_DEFAULT) {
738: PetscCall(PetscViewerASCIIPrintf(viewer, " factor row size %" PetscInt_FMT "\n", jac->factorrowsize));
739: } else {
740: PetscCall(PetscViewerASCIIPrintf(viewer, " default factor row size \n"));
741: }
742: }
743: PetscFunctionReturn(PETSC_SUCCESS);
744: }
746: static const char *HYPREILUType[] = {
747: "Block-Jacobi-ILUk", "Block-Jacobi-ILUT", "", "", "", "", "", "", "", "", /* 0-9 */
748: "GMRES-ILUk", "GMRES-ILUT", "", "", "", "", "", "", "", "", /* 10-19 */
749: "NSH-ILUk", "NSH-ILUT", "", "", "", "", "", "", "", "", /* 20-29 */
750: "RAS-ILUk", "RAS-ILUT", "", "", "", "", "", "", "", "", /* 30-39 */
751: "ddPQ-GMRES-ILUk", "ddPQ-GMRES-ILUT", "", "", "", "", "", "", "", "", /* 40-49 */
752: "GMRES-ILU0" /* 50 */
753: };
755: static const char *HYPREILUIterSetup[] = {"default", "async-in-place", "async-explicit", "sync-explicit", "semisync-explicit"};
757: static PetscErrorCode PCSetFromOptions_HYPRE_ILU(PC pc, PetscOptionItems PetscOptionsObject)
758: {
759: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
760: PetscBool flg;
761: PetscInt indx;
762: PetscReal tmpdbl;
763: PetscBool tmp_truth;
765: PetscFunctionBegin;
766: PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE ILU Options");
768: /* ILU: ILU Type */
769: PetscCall(PetscOptionsEList("-pc_hypre_ilu_type", "Choose ILU Type", "None", HYPREILUType, PETSC_STATIC_ARRAY_LENGTH(HYPREILUType), HYPREILUType[0], &indx, &flg));
770: if (flg) PetscCallHYPRE(HYPRE_ILUSetType(jac->hsolver, (HYPRE_Int)indx));
772: /* ILU: ILU iterative setup type*/
773: PetscCall(PetscOptionsEList("-pc_hypre_ilu_iterative_setup_type", "Set ILU iterative setup type", "None", HYPREILUIterSetup, PETSC_STATIC_ARRAY_LENGTH(HYPREILUIterSetup), HYPREILUIterSetup[0], &indx, &flg));
774: if (flg) PetscCallHYPRE(HYPRE_ILUSetIterativeSetupType(jac->hsolver, (HYPRE_Int)indx));
776: /* ILU: ILU iterative setup option*/
777: PetscCall(PetscOptionsInt("-pc_hypre_ilu_iterative_setup_option", "Set ILU iterative setup option", "None", 0, &indx, &flg));
778: if (flg) PetscCallHYPRE(HYPRE_ILUSetIterativeSetupOption(jac->hsolver, (HYPRE_Int)indx));
780: /* ILU: ILU iterative setup maxiter */
781: PetscCall(PetscOptionsInt("-pc_hypre_ilu_iterative_setup_maxiter", "Set ILU iterative setup maximum iteration count", "None", 0, &indx, &flg));
782: if (flg) PetscCallHYPRE(HYPRE_ILUSetIterativeSetupMaxIter(jac->hsolver, (HYPRE_Int)indx));
784: /* ILU: ILU iterative setup tolerance */
785: PetscCall(PetscOptionsReal("-pc_hypre_ilu_iterative_setup_tolerance", "Set ILU iterative setup tolerance", "None", 0, &tmpdbl, &flg));
786: if (flg) PetscCallHYPRE(HYPRE_ILUSetIterativeSetupTolerance(jac->hsolver, tmpdbl));
788: /* ILU: ILU Print Level */
789: PetscCall(PetscOptionsInt("-pc_hypre_ilu_print_level", "Set ILU print level", "None", 0, &indx, &flg));
790: if (flg) PetscCallHYPRE(HYPRE_ILUSetPrintLevel(jac->hsolver, (HYPRE_Int)indx));
792: /* ILU: Logging */
793: PetscCall(PetscOptionsInt("-pc_hypre_ilu_logging", "Set ILU logging level", "None", 0, &indx, &flg));
794: if (flg) PetscCallHYPRE(HYPRE_ILUSetLogging(jac->hsolver, (HYPRE_Int)indx));
796: /* ILU: ILU Level */
797: PetscCall(PetscOptionsInt("-pc_hypre_ilu_level", "Set ILU level", "None", 0, &indx, &flg));
798: if (flg) PetscCallHYPRE(HYPRE_ILUSetLevelOfFill(jac->hsolver, (HYPRE_Int)indx));
800: /* ILU: ILU Max NNZ per row */
801: PetscCall(PetscOptionsInt("-pc_hypre_ilu_max_nnz_per_row", "Set maximum NNZ per row", "None", 0, &indx, &flg));
802: if (flg) PetscCallHYPRE(HYPRE_ILUSetMaxNnzPerRow(jac->hsolver, (HYPRE_Int)indx));
804: /* ILU: tolerance */
805: PetscCall(PetscOptionsReal("-pc_hypre_ilu_tol", "Tolerance for ILU", "None", 0, &tmpdbl, &flg));
806: if (flg) PetscCallHYPRE(HYPRE_ILUSetTol(jac->hsolver, tmpdbl));
808: /* ILU: maximum iteration count */
809: PetscCall(PetscOptionsInt("-pc_hypre_ilu_maxiter", "Set ILU max iterations", "None", 0, &indx, &flg));
810: if (flg) PetscCallHYPRE(HYPRE_ILUSetMaxIter(jac->hsolver, (HYPRE_Int)indx));
812: /* ILU: drop threshold */
813: PetscCall(PetscOptionsReal("-pc_hypre_ilu_drop_threshold", "Drop threshold for ILU", "None", 0, &tmpdbl, &flg));
814: if (flg) PetscCallHYPRE(HYPRE_ILUSetDropThreshold(jac->hsolver, tmpdbl));
816: /* ILU: Triangular Solve */
817: PetscCall(PetscOptionsBool("-pc_hypre_ilu_tri_solve", "Enable triangular solve", "None", PETSC_FALSE, &tmp_truth, &flg));
818: if (flg) PetscCallHYPRE(HYPRE_ILUSetTriSolve(jac->hsolver, tmp_truth));
820: /* ILU: Lower Jacobi iteration */
821: PetscCall(PetscOptionsInt("-pc_hypre_ilu_lower_jacobi_iters", "Set lower Jacobi iteration count", "None", 0, &indx, &flg));
822: if (flg) PetscCallHYPRE(HYPRE_ILUSetLowerJacobiIters(jac->hsolver, (HYPRE_Int)indx));
824: /* ILU: Upper Jacobi iteration */
825: PetscCall(PetscOptionsInt("-pc_hypre_ilu_upper_jacobi_iters", "Set upper Jacobi iteration count", "None", 0, &indx, &flg));
826: if (flg) PetscCallHYPRE(HYPRE_ILUSetUpperJacobiIters(jac->hsolver, (HYPRE_Int)indx));
828: /* ILU: local reordering */
829: PetscCall(PetscOptionsBool("-pc_hypre_ilu_local_reordering", "Enable local reordering", "None", PETSC_FALSE, &tmp_truth, &flg));
830: if (flg) PetscCallHYPRE(HYPRE_ILUSetLocalReordering(jac->hsolver, tmp_truth));
832: PetscOptionsHeadEnd();
833: PetscFunctionReturn(PETSC_SUCCESS);
834: }
836: static PetscErrorCode PCView_HYPRE_ILU(PC pc, PetscViewer viewer)
837: {
838: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
839: hypre_ParILUData *ilu_data = (hypre_ParILUData *)jac->hsolver;
840: PetscBool isascii;
841: PetscInt indx;
842: PetscReal tmpdbl;
843: PetscReal *tmpdbl3;
845: PetscFunctionBegin;
846: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
847: if (isascii) {
848: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE ILU preconditioning\n"));
849: PetscCallExternalVoid("hypre_ParILUDataIluType", indx = hypre_ParILUDataIluType(ilu_data));
850: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU type %s (%" PetscInt_FMT ")\n", HYPREILUType[indx], indx));
851: PetscCallExternalVoid("hypre_ParILUDataLfil", indx = hypre_ParILUDataLfil(ilu_data));
852: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU level %" PetscInt_FMT "\n", indx));
853: PetscCallExternalVoid("hypre_ParILUDataMaxIter", indx = hypre_ParILUDataMaxIter(ilu_data));
854: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU max iterations %" PetscInt_FMT "\n", indx));
855: PetscCallExternalVoid("hypre_ParILUDataMaxRowNnz", indx = hypre_ParILUDataMaxRowNnz(ilu_data));
856: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU max NNZ per row %" PetscInt_FMT "\n", indx));
857: PetscCallExternalVoid("hypre_ParILUDataTriSolve", indx = hypre_ParILUDataTriSolve(ilu_data));
858: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU triangular solve %" PetscInt_FMT "\n", indx));
859: PetscCallExternalVoid("hypre_ParILUDataTol", tmpdbl = hypre_ParILUDataTol(ilu_data));
860: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU tolerance %e\n", tmpdbl));
861: PetscCallExternalVoid("hypre_ParILUDataDroptol", tmpdbl3 = hypre_ParILUDataDroptol(ilu_data));
862: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU drop tolerance %e / %e / %e\n", tmpdbl3[0], tmpdbl3[1], tmpdbl3[2]));
863: PetscCallExternalVoid("hypre_ParILUDataReorderingType", indx = hypre_ParILUDataReorderingType(ilu_data));
864: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU local reordering %" PetscInt_FMT "\n", indx));
865: PetscCallExternalVoid("hypre_ParILUDataLowerJacobiIters", indx = hypre_ParILUDataLowerJacobiIters(ilu_data));
866: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU lower Jacobi iterations %" PetscInt_FMT "\n", indx));
867: PetscCallExternalVoid("hypre_ParILUDataUpperJacobiIters", indx = hypre_ParILUDataUpperJacobiIters(ilu_data));
868: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU upper Jacobi iterations %" PetscInt_FMT "\n", indx));
869: PetscCallExternalVoid("hypre_ParILUDataPrintLevel", indx = hypre_ParILUDataPrintLevel(ilu_data));
870: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU print level %" PetscInt_FMT "\n", indx));
871: PetscCallExternalVoid("hypre_ParILUDataLogging", indx = hypre_ParILUDataLogging(ilu_data));
872: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU logging level %" PetscInt_FMT "\n", indx));
873: PetscCallExternalVoid("hypre_ParILUDataIterativeSetupType", indx = hypre_ParILUDataIterativeSetupType(ilu_data));
874: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup type %s (%" PetscInt_FMT ")\n", HYPREILUIterSetup[indx], indx));
875: PetscCallExternalVoid("hypre_ParILUDataIterativeSetupOption", indx = hypre_ParILUDataIterativeSetupOption(ilu_data));
876: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup option %" PetscInt_FMT "\n", indx));
877: PetscCallExternalVoid("hypre_ParILUDataIterativeSetupMaxIter", indx = hypre_ParILUDataIterativeSetupMaxIter(ilu_data));
878: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup max iterations %" PetscInt_FMT "\n", indx));
879: PetscCallExternalVoid("hypre_ParILUDataIterativeSetupTolerance", tmpdbl = hypre_ParILUDataIterativeSetupTolerance(ilu_data));
880: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup tolerance %e\n", tmpdbl));
881: }
882: PetscFunctionReturn(PETSC_SUCCESS);
883: }
885: static PetscErrorCode PCSetFromOptions_HYPRE_Euclid(PC pc, PetscOptionItems PetscOptionsObject)
886: {
887: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
888: PetscBool flag, eu_bj = jac->eu_bj ? PETSC_TRUE : PETSC_FALSE;
890: PetscFunctionBegin;
891: PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE Euclid Options");
892: PetscCall(PetscOptionsInt("-pc_hypre_euclid_level", "Factorization levels", "None", jac->eu_level, &jac->eu_level, &flag));
893: if (flag) PetscCallHYPRE(HYPRE_EuclidSetLevel(jac->hsolver, (HYPRE_Int)jac->eu_level));
895: PetscCall(PetscOptionsReal("-pc_hypre_euclid_droptolerance", "Drop tolerance for ILU(k) in Euclid", "None", jac->eu_droptolerance, &jac->eu_droptolerance, &flag));
896: if (flag) {
897: PetscMPIInt size;
899: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
900: PetscCheck(size == 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "hypre's Euclid does not support a parallel drop tolerance");
901: PetscCallHYPRE(HYPRE_EuclidSetILUT(jac->hsolver, jac->eu_droptolerance));
902: }
904: PetscCall(PetscOptionsBool("-pc_hypre_euclid_bj", "Use Block Jacobi for ILU in Euclid", "None", eu_bj, &eu_bj, &flag));
905: if (flag) {
906: jac->eu_bj = eu_bj ? 1 : 0;
907: PetscCallHYPRE(HYPRE_EuclidSetBJ(jac->hsolver, (HYPRE_Int)jac->eu_bj));
908: }
909: PetscOptionsHeadEnd();
910: PetscFunctionReturn(PETSC_SUCCESS);
911: }
913: static PetscErrorCode PCView_HYPRE_Euclid(PC pc, PetscViewer viewer)
914: {
915: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
916: PetscBool isascii;
918: PetscFunctionBegin;
919: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
920: if (isascii) {
921: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE Euclid preconditioning\n"));
922: if (jac->eu_level != PETSC_DEFAULT) {
923: PetscCall(PetscViewerASCIIPrintf(viewer, " factorization levels %" PetscInt_FMT "\n", jac->eu_level));
924: } else {
925: PetscCall(PetscViewerASCIIPrintf(viewer, " default factorization levels \n"));
926: }
927: PetscCall(PetscViewerASCIIPrintf(viewer, " drop tolerance %g\n", (double)jac->eu_droptolerance));
928: PetscCall(PetscViewerASCIIPrintf(viewer, " use Block-Jacobi? %" PetscInt_FMT "\n", jac->eu_bj));
929: }
930: PetscFunctionReturn(PETSC_SUCCESS);
931: }
933: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc, Vec b, Vec x)
934: {
935: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
936: Mat_HYPRE *hjac = (Mat_HYPRE *)jac->hpmat->data;
937: HYPRE_ParCSRMatrix hmat;
938: HYPRE_ParVector jbv, jxv;
940: PetscFunctionBegin;
941: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
942: PetscCall(VecSet(x, 0.0));
943: PetscCall(VecHYPRE_IJVectorPushVecRead(hjac->b, b));
944: PetscCall(VecHYPRE_IJVectorPushVecWrite(hjac->x, x));
946: PetscCallHYPRE(HYPRE_IJMatrixGetObject(hjac->ij, (void **)&hmat));
947: PetscCallHYPRE(HYPRE_IJVectorGetObject(hjac->b->ij, (void **)&jbv));
948: PetscCallHYPRE(HYPRE_IJVectorGetObject(hjac->x->ij, (void **)&jxv));
950: PetscCallExternalVoid(
951: "Hypre Transpose solve", do {
952: HYPRE_Int hierr = HYPRE_BoomerAMGSolveT(jac->hsolver, hmat, jbv, jxv);
953: if (hierr) {
954: /* error code of 1 in BoomerAMG merely means convergence not achieved */
955: PetscCheck(hierr == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in HYPRE solver, error code %d", (int)hierr);
956: HYPRE_ClearAllErrors();
957: }
958: } while (0));
960: PetscCall(VecHYPRE_IJVectorPopVec(hjac->x));
961: PetscCall(VecHYPRE_IJVectorPopVec(hjac->b));
962: PetscFunctionReturn(PETSC_SUCCESS);
963: }
965: static PetscErrorCode PCMGGalerkinGetMatProductAlgorithm_HYPRE_BoomerAMG(PC pc, const char *spgemm[])
966: {
967: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
969: PetscFunctionBegin;
971: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
972: *spgemm = jac->spgemm_type;
973: #endif
974: PetscFunctionReturn(PETSC_SUCCESS);
975: }
977: static const char *HYPREBoomerAMGCycleType[] = {"", "V", "W"};
978: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP", "Ruge-Stueben", "", "modifiedRuge-Stueben", "", "", "Falgout", "", "PMIS", "", "HMIS"};
979: static const char *HYPREBoomerAMGMeasureType[] = {"local", "global"};
980: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
981: static const char *HYPREBoomerAMGSmoothType[] = {"ILU", "Schwarz-smoothers", "Pilut", "ParaSails", "Euclid"};
982: 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"};
983: 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"};
985: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PC pc, PetscOptionItems PetscOptionsObject)
986: {
987: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
988: PetscInt bs, n, indx, level;
989: PetscBool flg, tmp_truth;
990: PetscReal tmpdbl, twodbl[2];
991: const char *symtlist[] = {"nonsymmetric", "SPD", "nonsymmetric,SPD"};
993: PetscFunctionBegin;
994: PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE BoomerAMG Options");
995: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_cycle_type", "Cycle type", "None", HYPREBoomerAMGCycleType + 1, 2, HYPREBoomerAMGCycleType[jac->cycletype], &indx, &flg));
996: if (flg) {
997: jac->cycletype = indx + 1;
998: PetscCallHYPRE(HYPRE_BoomerAMGSetCycleType(jac->hsolver, (HYPRE_Int)jac->cycletype));
999: }
1000: PetscCall(PetscOptionsBoundedInt("-pc_hypre_boomeramg_max_levels", "Number of levels (of grids) allowed", "None", jac->maxlevels, &jac->maxlevels, &flg, 2));
1001: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetMaxLevels(jac->hsolver, (HYPRE_Int)jac->maxlevels));
1002: PetscCall(PetscOptionsBoundedInt("-pc_hypre_boomeramg_max_iter", "Maximum iterations used PER hypre call", "None", jac->maxiter, &jac->maxiter, &flg, 1));
1003: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetMaxIter(jac->hsolver, (HYPRE_Int)jac->maxiter));
1004: PetscCall(PetscOptionsBoundedReal("-pc_hypre_boomeramg_tol", "Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)", "None", jac->tol, &jac->tol, &flg, 0.0));
1005: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetTol(jac->hsolver, jac->tol));
1006: bs = 1;
1007: if (pc->pmat) PetscCall(MatGetBlockSize(pc->pmat, &bs));
1008: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_numfunctions", "Number of functions", "HYPRE_BoomerAMGSetNumFunctions", bs, &bs, &flg));
1009: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetNumFunctions(jac->hsolver, (HYPRE_Int)bs));
1011: PetscCall(PetscOptionsBoundedReal("-pc_hypre_boomeramg_truncfactor", "Truncation factor for interpolation (0=no truncation)", "None", jac->truncfactor, &jac->truncfactor, &flg, 0.0));
1012: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetTruncFactor(jac->hsolver, jac->truncfactor));
1014: PetscCall(PetscOptionsBoundedInt("-pc_hypre_boomeramg_P_max", "Max elements per row for interpolation operator (0=unlimited)", "None", jac->pmax, &jac->pmax, &flg, 0));
1015: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetPMaxElmts(jac->hsolver, (HYPRE_Int)jac->pmax));
1017: PetscCall(PetscOptionsRangeInt("-pc_hypre_boomeramg_agg_nl", "Number of levels of aggressive coarsening", "None", jac->agg_nl, &jac->agg_nl, &flg, 0, jac->maxlevels));
1018: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetAggNumLevels(jac->hsolver, (HYPRE_Int)jac->agg_nl));
1020: PetscCall(PetscOptionsBoundedInt("-pc_hypre_boomeramg_agg_num_paths", "Number of paths for aggressive coarsening", "None", jac->agg_num_paths, &jac->agg_num_paths, &flg, 1));
1021: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetNumPaths(jac->hsolver, (HYPRE_Int)jac->agg_num_paths));
1023: PetscCall(PetscOptionsBoundedReal("-pc_hypre_boomeramg_strong_threshold", "Threshold for being strongly connected", "None", jac->strongthreshold, &jac->strongthreshold, &flg, 0.0));
1024: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetStrongThreshold(jac->hsolver, jac->strongthreshold));
1025: PetscCall(PetscOptionsRangeReal("-pc_hypre_boomeramg_max_row_sum", "Maximum row sum", "None", jac->maxrowsum, &jac->maxrowsum, &flg, 0.0, 1.0));
1026: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetMaxRowSum(jac->hsolver, jac->maxrowsum));
1028: /* Grid sweeps */
1029: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all", "Number of sweeps for the up and down grid levels", "None", jac->gridsweeps[0], &indx, &flg));
1030: if (flg) {
1031: /* modify the jac structure so we can view the updated options with PC_View */
1032: jac->gridsweeps[0] = indx;
1033: jac->gridsweeps[1] = indx;
1034: /*defaults coarse to 1 */
1035: jac->gridsweeps[2] = 1;
1036: }
1037: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen", "Use a nodal based coarsening 1-6", "HYPRE_BoomerAMGSetNodal", jac->nodal_coarsening, &jac->nodal_coarsening, &flg));
1038: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetNodal(jac->hsolver, (HYPRE_Int)jac->nodal_coarsening));
1039: 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));
1040: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetNodalDiag(jac->hsolver, (HYPRE_Int)jac->nodal_coarsening_diag));
1041: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_variant", "Variant of algorithm 1-3", "HYPRE_BoomerAMGSetInterpVecVariant", jac->vec_interp_variant, &jac->vec_interp_variant, &flg));
1042: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetInterpVecVariant(jac->hsolver, (HYPRE_Int)jac->vec_interp_variant));
1043: 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));
1044: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetInterpVecQMax(jac->hsolver, (HYPRE_Int)jac->vec_interp_qmax));
1045: 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));
1046: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetSmoothInterpVectors(jac->hsolver, jac->vec_interp_smooth));
1047: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_interp_refine", "Preprocess the interpolation matrix through iterative weight refinement", "HYPRE_BoomerAMGSetInterpRefine", jac->interp_refine, &jac->interp_refine, &flg));
1048: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetInterpRefine(jac->hsolver, (HYPRE_Int)jac->interp_refine));
1049: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down", "Number of sweeps for the down cycles", "None", jac->gridsweeps[0], &indx, &flg));
1050: if (flg) {
1051: PetscCallHYPRE(HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver, (HYPRE_Int)indx, 1));
1052: jac->gridsweeps[0] = indx;
1053: }
1054: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up", "Number of sweeps for the up cycles", "None", jac->gridsweeps[1], &indx, &flg));
1055: if (flg) {
1056: PetscCallHYPRE(HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver, (HYPRE_Int)indx, 2));
1057: jac->gridsweeps[1] = indx;
1058: }
1059: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse", "Number of sweeps for the coarse level", "None", jac->gridsweeps[2], &indx, &flg));
1060: if (flg) {
1061: PetscCallHYPRE(HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver, (HYPRE_Int)indx, 3));
1062: jac->gridsweeps[2] = indx;
1063: }
1065: /* Smooth type */
1066: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_smooth_type", "Enable more complex smoothers", "None", HYPREBoomerAMGSmoothType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGSmoothType), HYPREBoomerAMGSmoothType[0], &indx, &flg));
1067: if (flg) {
1068: jac->smoothtype = indx;
1069: PetscCallHYPRE(HYPRE_BoomerAMGSetSmoothType(jac->hsolver, (HYPRE_Int)indx + 5));
1070: jac->smoothnumlevels = 25;
1071: PetscCallHYPRE(HYPRE_BoomerAMGSetSmoothNumLevels(jac->hsolver, 25));
1072: }
1074: /* Number of smoothing levels */
1075: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_levels", "Number of levels on which more complex smoothers are used", "None", 25, &indx, &flg));
1076: if (flg && (jac->smoothtype != -1)) {
1077: jac->smoothnumlevels = indx;
1078: PetscCallHYPRE(HYPRE_BoomerAMGSetSmoothNumLevels(jac->hsolver, (HYPRE_Int)indx));
1079: }
1081: /* Smooth num sweeps */
1082: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_sweeps", "Set number of smoother sweeps", "None", 1, &indx, &flg));
1083: if (flg && indx > 0) {
1084: jac->smoothsweeps = indx;
1085: PetscCallHYPRE(HYPRE_BoomerAMGSetSmoothNumSweeps(jac->hsolver, (HYPRE_Int)indx));
1086: }
1088: /* ILU: ILU Type */
1089: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_ilu_type", "Choose ILU Type", "None", HYPREILUType, PETSC_STATIC_ARRAY_LENGTH(HYPREILUType), HYPREILUType[0], &indx, &flg));
1090: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUType(jac->hsolver, (HYPRE_Int)indx));
1092: /* ILU: ILU iterative setup type*/
1093: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_ilu_iterative_setup_type", "Set ILU iterative setup type", "None", HYPREILUIterSetup, PETSC_STATIC_ARRAY_LENGTH(HYPREILUIterSetup), HYPREILUIterSetup[0], &indx, &flg));
1094: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUIterSetupType(jac->hsolver, (HYPRE_Int)indx));
1096: /* ILU: ILU iterative setup option*/
1097: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_iterative_setup_option", "Set ILU iterative setup option", "None", 0, &indx, &flg));
1098: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUIterSetupOption(jac->hsolver, (HYPRE_Int)indx));
1100: /* ILU: ILU iterative setup maxiter */
1101: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_iterative_setup_maxiter", "Set ILU iterative setup maximum iteration count", "None", 0, &indx, &flg));
1102: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUIterSetupMaxIter(jac->hsolver, (HYPRE_Int)indx));
1104: /* ILU: ILU iterative setup tolerance */
1105: PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_ilu_iterative_setup_tolerance", "Set ILU iterative setup tolerance", "None", 0, &tmpdbl, &flg));
1106: if (flg) PetscCallHYPRE(hypre_BoomerAMGSetILUIterSetupTolerance(jac->hsolver, tmpdbl));
1108: /* ILU: ILU Print Level */
1109: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_print_level", "Set ILU print level", "None", 0, &indx, &flg));
1110: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetPrintLevel(jac->hsolver, (HYPRE_Int)indx));
1112: /* ILU: Logging */
1113: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_logging", "Set ILU logging level", "None", 0, &indx, &flg));
1114: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetLogging(jac->hsolver, (HYPRE_Int)indx));
1116: /* ILU: ILU Level */
1117: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_level", "Set ILU level", "None", 0, &indx, &flg));
1118: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILULevel(jac->hsolver, (HYPRE_Int)indx));
1120: /* ILU: ILU Max NNZ per row */
1121: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_max_nnz_per_row", "Set maximum NNZ per row", "None", 0, &indx, &flg));
1122: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUMaxRowNnz(jac->hsolver, (HYPRE_Int)indx));
1124: /* ILU: maximum iteration count */
1125: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_maxiter", "Set ILU max iterations", "None", 0, &indx, &flg));
1126: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUMaxIter(jac->hsolver, (HYPRE_Int)indx));
1128: /* ILU: drop threshold */
1129: PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_ilu_drop_tol", "Drop tolerance for ILU", "None", 0, &tmpdbl, &flg));
1130: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUDroptol(jac->hsolver, tmpdbl));
1132: /* ILU: Triangular Solve */
1133: PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_ilu_tri_solve", "Enable triangular solve", "None", PETSC_FALSE, &tmp_truth, &flg));
1134: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUTriSolve(jac->hsolver, tmp_truth));
1136: /* ILU: Lower Jacobi iteration */
1137: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_lower_jacobi_iters", "Set lower Jacobi iteration count", "None", 0, &indx, &flg));
1138: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILULowerJacobiIters(jac->hsolver, (HYPRE_Int)indx));
1140: /* ILU: Upper Jacobi iteration */
1141: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_upper_jacobi_iters", "Set upper Jacobi iteration count", "None", 0, &indx, &flg));
1142: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILUUpperJacobiIters(jac->hsolver, (HYPRE_Int)indx));
1144: /* ILU: local reordering */
1145: PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_ilu_local_reordering", "Enable local reordering", "None", PETSC_FALSE, &tmp_truth, &flg));
1146: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetILULocalReordering(jac->hsolver, tmp_truth));
1148: /* Number of levels for ILU(k) for Euclid */
1149: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_eu_level", "Number of levels for ILU(k) in Euclid smoother", "None", 0, &indx, &flg));
1150: if (flg && (jac->smoothtype == 4)) {
1151: jac->eu_level = indx;
1152: PetscCallHYPRE(HYPRE_BoomerAMGSetEuLevel(jac->hsolver, (HYPRE_Int)indx));
1153: }
1155: /* Filter for ILU(k) for Euclid */
1156: PetscReal droptolerance;
1157: PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_eu_droptolerance", "Drop tolerance for ILU(k) in Euclid smoother", "None", 0, &droptolerance, &flg));
1158: if (flg && (jac->smoothtype == 4)) {
1159: jac->eu_droptolerance = droptolerance;
1160: PetscCallHYPRE(HYPRE_BoomerAMGSetEuLevel(jac->hsolver, droptolerance));
1161: }
1163: /* Use Block Jacobi ILUT for Euclid */
1164: PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_eu_bj", "Use Block Jacobi for ILU in Euclid smoother?", "None", PETSC_FALSE, &tmp_truth, &flg));
1165: if (flg && (jac->smoothtype == 4)) {
1166: jac->eu_bj = tmp_truth;
1167: PetscCallHYPRE(HYPRE_BoomerAMGSetEuBJ(jac->hsolver, (HYPRE_Int)jac->eu_bj));
1168: }
1170: /* Relax type */
1171: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all", "Relax type for the up and down cycles", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType),
1172: jac->relaxtype[0] < 0 ? "not yet set" : HYPREBoomerAMGRelaxType[jac->relaxtype[0]], &indx, &flg));
1173: if (flg) jac->relaxtype[0] = jac->relaxtype[1] = indx;
1174: PetscCall(
1175: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down", "Relax type for the down cycles", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), jac->relaxtype[0] < 0 ? "not yet set" : HYPREBoomerAMGRelaxType[jac->relaxtype[0]], &indx, &flg));
1176: if (flg) jac->relaxtype[0] = indx;
1177: PetscCall(
1178: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up", "Relax type for the up cycles", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), jac->relaxtype[1] < 0 ? "not yet set" : HYPREBoomerAMGRelaxType[jac->relaxtype[1]], &indx, &flg));
1179: if (flg) jac->relaxtype[1] = indx;
1180: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse", "Relax type on coarse grid", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), HYPREBoomerAMGRelaxType[jac->relaxtype[2]], &indx, &flg));
1181: if (flg) jac->relaxtype[2] = indx;
1183: /* Relaxation Weight */
1184: 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));
1185: if (flg) {
1186: PetscCallHYPRE(HYPRE_BoomerAMGSetRelaxWt(jac->hsolver, tmpdbl));
1187: jac->relaxweight = tmpdbl;
1188: }
1190: n = 2;
1191: twodbl[0] = twodbl[1] = 1.0;
1192: PetscCall(PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level", "Set the relaxation weight for a particular level (weight,level)", "None", twodbl, &n, &flg));
1193: if (flg) {
1194: 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);
1195: indx = (int)PetscAbsReal(twodbl[1]);
1196: PetscCallHYPRE(HYPRE_BoomerAMGSetLevelRelaxWt(jac->hsolver, twodbl[0], (HYPRE_Int)indx));
1197: }
1199: /* Outer relaxation Weight */
1200: 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));
1201: if (flg) {
1202: PetscCallHYPRE(HYPRE_BoomerAMGSetOuterWt(jac->hsolver, tmpdbl));
1203: jac->outerrelaxweight = tmpdbl;
1204: }
1206: n = 2;
1207: twodbl[0] = twodbl[1] = 1.0;
1208: PetscCall(PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level", "Set the outer relaxation weight for a particular level (weight,level)", "None", twodbl, &n, &flg));
1209: if (flg) {
1210: 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);
1211: indx = (int)PetscAbsReal(twodbl[1]);
1212: PetscCallHYPRE(HYPRE_BoomerAMGSetLevelOuterWt(jac->hsolver, twodbl[0], (HYPRE_Int)indx));
1213: }
1215: /* the Relax Order */
1216: PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", PETSC_FALSE, &tmp_truth, &flg));
1217: if (flg) jac->relaxorder = !tmp_truth;
1218: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_measure_type", "Measure type", "None", HYPREBoomerAMGMeasureType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGMeasureType), HYPREBoomerAMGMeasureType[0], &indx, &flg));
1219: if (flg) {
1220: jac->measuretype = indx;
1221: PetscCallHYPRE(HYPRE_BoomerAMGSetMeasureType(jac->hsolver, (HYPRE_Int)jac->measuretype));
1222: }
1223: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type", "Coarsen type", "None", HYPREBoomerAMGCoarsenType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGCoarsenType), jac->coarsentype < 0 ? "unknown" : HYPREBoomerAMGCoarsenType[jac->coarsentype], &indx, &flg));
1224: if (flg) jac->coarsentype = indx;
1226: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_max_coarse_size", "Maximum size of coarsest grid", "None", jac->maxc, &jac->maxc, &flg));
1227: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetMaxCoarseSize(jac->hsolver, (HYPRE_Int)jac->maxc));
1228: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_min_coarse_size", "Minimum size of coarsest grid", "None", jac->minc, &jac->minc, &flg));
1229: if (flg) PetscCallHYPRE(HYPRE_BoomerAMGSetMinCoarseSize(jac->hsolver, (HYPRE_Int)jac->minc));
1230: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
1231: // global parameter but is closely associated with BoomerAMG
1232: PetscCall(PetscOptionsEList("-pc_mg_galerkin_mat_product_algorithm", "Type of SpGEMM to use in hypre (only for now)", "PCMGGalerkinSetMatProductAlgorithm", HYPRESpgemmTypes, PETSC_STATIC_ARRAY_LENGTH(HYPRESpgemmTypes), jac->spgemm_type, &indx, &flg));
1233: if (flg) PetscCall(PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG(pc, HYPRESpgemmTypes[indx]));
1234: #endif
1235: /* AIR */
1236: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
1237: 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));
1238: PetscCallHYPRE(HYPRE_BoomerAMGSetRestriction(jac->hsolver, (HYPRE_Int)jac->Rtype));
1239: if (jac->Rtype) {
1240: HYPRE_Int **grid_relax_points = hypre_TAlloc(HYPRE_Int *, 4, HYPRE_MEMORY_HOST);
1241: char *prerelax[256];
1242: char *postrelax[256];
1243: char stringF[2] = "F", stringC[2] = "C", stringA[2] = "A";
1244: PetscInt ns_down = 256, ns_up = 256;
1245: PetscBool matchF, matchC, matchA;
1247: 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 */
1249: PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_strongthresholdR", "Threshold for R", "None", jac->Rstrongthreshold, &jac->Rstrongthreshold, NULL));
1250: PetscCallHYPRE(HYPRE_BoomerAMGSetStrongThresholdR(jac->hsolver, jac->Rstrongthreshold));
1252: PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_filterthresholdR", "Filter threshold for R", "None", jac->Rfilterthreshold, &jac->Rfilterthreshold, NULL));
1253: PetscCallHYPRE(HYPRE_BoomerAMGSetFilterThresholdR(jac->hsolver, jac->Rfilterthreshold));
1255: 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));
1256: PetscCallHYPRE(HYPRE_BoomerAMGSetADropTol(jac->hsolver, (HYPRE_Int)jac->Adroptol));
1258: 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));
1259: PetscCallHYPRE(HYPRE_BoomerAMGSetADropType(jac->hsolver, (HYPRE_Int)jac->Adroptype));
1260: PetscCall(PetscOptionsStringArray("-pc_hypre_boomeramg_prerelax", "Defines prerelax scheme", "None", prerelax, &ns_down, NULL));
1261: PetscCall(PetscOptionsStringArray("-pc_hypre_boomeramg_postrelax", "Defines postrelax scheme", "None", postrelax, &ns_up, NULL));
1262: PetscCheck(ns_down == jac->gridsweeps[0], PetscObjectComm((PetscObject)jac), PETSC_ERR_ARG_SIZ, "The number of arguments passed to -pc_hypre_boomeramg_prerelax must match the number passed to -pc_hypre_bomeramg_grid_sweeps_down");
1263: PetscCheck(ns_up == jac->gridsweeps[1], PetscObjectComm((PetscObject)jac), PETSC_ERR_ARG_SIZ, "The number of arguments passed to -pc_hypre_boomeramg_postrelax must match the number passed to -pc_hypre_bomeramg_grid_sweeps_up");
1265: grid_relax_points[0] = NULL;
1266: grid_relax_points[1] = hypre_TAlloc(HYPRE_Int, ns_down, HYPRE_MEMORY_HOST);
1267: grid_relax_points[2] = hypre_TAlloc(HYPRE_Int, ns_up, HYPRE_MEMORY_HOST);
1268: grid_relax_points[3] = hypre_TAlloc(HYPRE_Int, jac->gridsweeps[2], HYPRE_MEMORY_HOST);
1269: grid_relax_points[3][0] = 0;
1271: // set down relax scheme
1272: for (PetscInt i = 0; i < ns_down; i++) {
1273: PetscCall(PetscStrcasecmp(prerelax[i], stringF, &matchF));
1274: PetscCall(PetscStrcasecmp(prerelax[i], stringC, &matchC));
1275: PetscCall(PetscStrcasecmp(prerelax[i], stringA, &matchA));
1276: PetscCheck(matchF || matchC || matchA, PetscObjectComm((PetscObject)jac), PETSC_ERR_ARG_WRONG, "Valid argument options for -pc_hypre_boomeramg_prerelax are C, F, and A");
1277: if (matchF) grid_relax_points[1][i] = -1;
1278: else if (matchC) grid_relax_points[1][i] = 1;
1279: else if (matchA) grid_relax_points[1][i] = 0;
1280: }
1282: // set up relax scheme
1283: for (PetscInt i = 0; i < ns_up; i++) {
1284: PetscCall(PetscStrcasecmp(postrelax[i], stringF, &matchF));
1285: PetscCall(PetscStrcasecmp(postrelax[i], stringC, &matchC));
1286: PetscCall(PetscStrcasecmp(postrelax[i], stringA, &matchA));
1287: PetscCheck(matchF || matchC || matchA, PetscObjectComm((PetscObject)jac), PETSC_ERR_ARG_WRONG, "Valid argument options for -pc_hypre_boomeramg_postrelax are C, F, and A");
1288: if (matchF) grid_relax_points[2][i] = -1;
1289: else if (matchC) grid_relax_points[2][i] = 1;
1290: else if (matchA) grid_relax_points[2][i] = 0;
1291: }
1293: // set coarse relax scheme
1294: for (PetscInt i = 0; i < jac->gridsweeps[2]; i++) grid_relax_points[3][i] = 0;
1296: // Pass relax schemes to hypre
1297: PetscCallHYPRE(HYPRE_BoomerAMGSetGridRelaxPoints(jac->hsolver, grid_relax_points));
1299: // cleanup memory
1300: for (PetscInt i = 0; i < ns_down; i++) PetscCall(PetscFree(prerelax[i]));
1301: for (PetscInt i = 0; i < ns_up; i++) PetscCall(PetscFree(postrelax[i]));
1302: }
1303: #endif
1305: #if PETSC_PKG_HYPRE_VERSION_LE(9, 9, 9)
1306: 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);
1307: #endif
1309: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_interp_type", "Interpolation type", "None", HYPREBoomerAMGInterpType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGInterpType),
1310: jac->interptype == 100 ? "1pt" : (jac->interptype < 0 ? "unknown" : HYPREBoomerAMGInterpType[jac->interptype]), &indx, &flg));
1311: if (flg) jac->interptype = indx;
1313: PetscCall(PetscOptionsName("-pc_hypre_boomeramg_print_statistics", "Print statistics", "None", &flg));
1314: if (flg) {
1315: level = 3;
1316: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_print_statistics", "Print statistics", "None", level, &level, NULL));
1318: jac->printstatistics = PETSC_TRUE;
1319: PetscCallHYPRE(HYPRE_BoomerAMGSetPrintLevel(jac->hsolver, (HYPRE_Int)level));
1320: }
1322: PetscCall(PetscOptionsName("-pc_hypre_boomeramg_print_debug", "Print debug information", "None", &flg));
1323: if (flg) {
1324: level = 3;
1325: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_print_debug", "Print debug information", "None", level, &level, NULL));
1327: jac->printstatistics = PETSC_TRUE;
1328: PetscCallHYPRE(HYPRE_BoomerAMGSetDebugFlag(jac->hsolver, (HYPRE_Int)level));
1329: }
1331: PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg));
1332: if (flg && tmp_truth) {
1333: PetscInt tmp_int;
1334: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", (HYPRE_Int)jac->nodal_relax_levels, &tmp_int, &flg));
1335: if (flg) jac->nodal_relax_levels = tmp_int;
1336: PetscCallHYPRE(HYPRE_BoomerAMGSetSmoothType(jac->hsolver, 6));
1337: PetscCallHYPRE(HYPRE_BoomerAMGSetDomainType(jac->hsolver, 1));
1338: PetscCallHYPRE(HYPRE_BoomerAMGSetOverlap(jac->hsolver, 0));
1339: PetscCallHYPRE(HYPRE_BoomerAMGSetSmoothNumLevels(jac->hsolver, (HYPRE_Int)jac->nodal_relax_levels));
1340: }
1342: PetscCall(PetscOptionsBool3("-pc_hypre_boomeramg_keeptranspose", "Avoid transpose matvecs in preconditioner application", "None", jac->keeptranspose, &jac->keeptranspose, NULL));
1344: /* options for ParaSails solvers */
1345: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_parasails_sym", "Symmetry of matrix and preconditioner", "None", symtlist, PETSC_STATIC_ARRAY_LENGTH(symtlist), symtlist[0], &indx, &flg));
1346: if (flg) {
1347: jac->symt = indx;
1348: PetscCallHYPRE(HYPRE_BoomerAMGSetSym(jac->hsolver, (HYPRE_Int)jac->symt));
1349: }
1351: PetscOptionsHeadEnd();
1352: PetscFunctionReturn(PETSC_SUCCESS);
1353: }
1355: 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)
1356: {
1357: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1358: HYPRE_Int oits;
1360: PetscFunctionBegin;
1361: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
1362: PetscCallHYPRE(HYPRE_BoomerAMGSetMaxIter(jac->hsolver, (HYPRE_Int)(its * jac->maxiter)));
1363: PetscCallHYPRE(HYPRE_BoomerAMGSetTol(jac->hsolver, rtol));
1364: jac->applyrichardson = PETSC_TRUE;
1365: PetscCall(PCApply_HYPRE(pc, b, y));
1366: jac->applyrichardson = PETSC_FALSE;
1367: PetscCallHYPRE(HYPRE_BoomerAMGGetNumIterations(jac->hsolver, &oits));
1368: *outits = oits;
1369: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1370: else *reason = PCRICHARDSON_CONVERGED_RTOL;
1371: PetscCallHYPRE(HYPRE_BoomerAMGSetTol(jac->hsolver, jac->tol));
1372: PetscCallHYPRE(HYPRE_BoomerAMGSetMaxIter(jac->hsolver, (HYPRE_Int)jac->maxiter));
1373: PetscFunctionReturn(PETSC_SUCCESS);
1374: }
1376: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc, PetscViewer viewer)
1377: {
1378: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1379: hypre_ParAMGData *amg_data = (hypre_ParAMGData *)jac->hsolver;
1380: PetscBool isascii;
1381: PetscInt indx;
1382: PetscReal val;
1384: PetscFunctionBegin;
1385: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1386: if (isascii) {
1387: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE BoomerAMG preconditioning\n"));
1388: PetscCall(PetscViewerASCIIPrintf(viewer, " Cycle type %s\n", HYPREBoomerAMGCycleType[jac->cycletype]));
1389: PetscCall(PetscViewerASCIIPrintf(viewer, " Maximum number of levels %" PetscInt_FMT "\n", jac->maxlevels));
1390: PetscCall(PetscViewerASCIIPrintf(viewer, " Maximum number of iterations PER hypre call %" PetscInt_FMT "\n", jac->maxiter));
1391: PetscCall(PetscViewerASCIIPrintf(viewer, " Convergence tolerance PER hypre call %g\n", (double)jac->tol));
1392: PetscCall(PetscViewerASCIIPrintf(viewer, " Threshold for strong coupling %g\n", (double)jac->strongthreshold));
1393: PetscCall(PetscViewerASCIIPrintf(viewer, " Interpolation truncation factor %g\n", (double)jac->truncfactor));
1394: PetscCall(PetscViewerASCIIPrintf(viewer, " Interpolation: max elements per row %" PetscInt_FMT "\n", jac->pmax));
1395: if (jac->interp_refine) PetscCall(PetscViewerASCIIPrintf(viewer, " Interpolation: number of steps of weighted refinement %" PetscInt_FMT "\n", jac->interp_refine));
1396: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of levels of aggressive coarsening %" PetscInt_FMT "\n", jac->agg_nl));
1397: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of paths for aggressive coarsening %" PetscInt_FMT "\n", jac->agg_num_paths));
1399: PetscCall(PetscViewerASCIIPrintf(viewer, " Maximum row sums %g\n", (double)jac->maxrowsum));
1401: PetscCall(PetscViewerASCIIPrintf(viewer, " Sweeps down %" PetscInt_FMT "\n", jac->gridsweeps[0]));
1402: PetscCall(PetscViewerASCIIPrintf(viewer, " Sweeps up %" PetscInt_FMT "\n", jac->gridsweeps[1]));
1403: PetscCall(PetscViewerASCIIPrintf(viewer, " Sweeps on coarse %" PetscInt_FMT "\n", jac->gridsweeps[2]));
1405: PetscCall(PetscViewerASCIIPrintf(viewer, " Relax down %s\n", jac->relaxtype[0] < 0 ? "not yet set" : HYPREBoomerAMGRelaxType[jac->relaxtype[0]]));
1406: PetscCall(PetscViewerASCIIPrintf(viewer, " Relax up %s\n", jac->relaxtype[1] < 0 ? "not yet set" : HYPREBoomerAMGRelaxType[jac->relaxtype[1]]));
1407: PetscCall(PetscViewerASCIIPrintf(viewer, " Relax on coarse %s\n", HYPREBoomerAMGRelaxType[jac->relaxtype[2]]));
1409: PetscCall(PetscViewerASCIIPrintf(viewer, " Relax weight (all) %g\n", (double)jac->relaxweight));
1410: PetscCall(PetscViewerASCIIPrintf(viewer, " Outer relax weight (all) %g\n", (double)jac->outerrelaxweight));
1412: PetscCall(PetscViewerASCIIPrintf(viewer, " Maximum size of coarsest grid %" PetscInt_FMT "\n", jac->maxc));
1413: PetscCall(PetscViewerASCIIPrintf(viewer, " Minimum size of coarsest grid %" PetscInt_FMT "\n", jac->minc));
1415: if (jac->relaxorder == PETSC_DECIDE) {
1416: PetscCall(PetscViewerASCIIPrintf(viewer, " CF-relaxation option not yet determined\n"));
1417: } else if (jac->relaxorder) {
1418: PetscCall(PetscViewerASCIIPrintf(viewer, " Using CF-relaxation\n"));
1419: } else {
1420: PetscCall(PetscViewerASCIIPrintf(viewer, " Not using CF-relaxation\n"));
1421: }
1422: if (jac->smoothtype != -1) {
1423: PetscCall(PetscViewerASCIIPrintf(viewer, " Smooth type %s\n", HYPREBoomerAMGSmoothType[jac->smoothtype]));
1424: PetscCall(PetscViewerASCIIPrintf(viewer, " Smooth num levels %" PetscInt_FMT "\n", jac->smoothnumlevels));
1425: PetscCall(PetscViewerASCIIPrintf(viewer, " Smooth num sweeps %" PetscInt_FMT "\n", jac->smoothsweeps));
1426: if (jac->smoothtype == 0) {
1427: PetscCallExternalVoid("hypre_ParAMGDataILUType", indx = hypre_ParAMGDataILUType(amg_data));
1428: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU type %s (%" PetscInt_FMT ")\n", HYPREILUType[indx], indx));
1429: PetscCallExternalVoid("hypre_ParAMGDataILULevel", indx = hypre_ParAMGDataILULevel(amg_data));
1430: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU level %" PetscInt_FMT "\n", indx));
1431: PetscCallExternalVoid("hypre_ParAMGDataILUMaxIter", indx = hypre_ParAMGDataILUMaxIter(amg_data));
1432: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU max iterations %" PetscInt_FMT "\n", indx));
1433: PetscCallExternalVoid("hypre_ParAMGDataILUMaxRowNnz", indx = hypre_ParAMGDataILUMaxRowNnz(amg_data));
1434: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU max NNZ per row %" PetscInt_FMT "\n", indx));
1435: PetscCallExternalVoid("hypre_ParAMGDataILUTriSolve", indx = hypre_ParAMGDataILUTriSolve(amg_data));
1436: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU triangular solve %" PetscInt_FMT "\n", indx));
1437: PetscCallExternalVoid("hypre_ParAMGDataTol", val = hypre_ParAMGDataTol(amg_data));
1438: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU tolerance %e\n", val));
1439: PetscCallExternalVoid("hypre_ParAMGDataILUDroptol", val = hypre_ParAMGDataILUDroptol(amg_data));
1440: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU drop tolerance %e\n", val));
1441: PetscCallExternalVoid("hypre_ParAMGDataILULocalReordering", indx = hypre_ParAMGDataILULocalReordering(amg_data));
1442: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU local reordering %" PetscInt_FMT "\n", indx));
1443: PetscCallExternalVoid("hypre_ParAMGDataILULowerJacobiIters", indx = hypre_ParAMGDataILULowerJacobiIters(amg_data));
1444: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU lower Jacobi iterations %" PetscInt_FMT "\n", indx));
1445: PetscCallExternalVoid("hypre_ParAMGDataILUUpperJacobiIters", indx = hypre_ParAMGDataILUUpperJacobiIters(amg_data));
1446: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU upper Jacobi iterations %" PetscInt_FMT "\n", indx));
1447: PetscCallExternalVoid("hypre_ParAMGDataPrintLevel", indx = hypre_ParAMGDataPrintLevel(amg_data));
1448: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU print level %" PetscInt_FMT "\n", indx));
1449: PetscCallExternalVoid("hypre_ParAMGDataLogging", indx = hypre_ParAMGDataLogging(amg_data));
1450: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU logging level %" PetscInt_FMT "\n", indx));
1451: PetscCallExternalVoid("hypre_ParAMGDataILUIterSetupType", indx = hypre_ParAMGDataILUIterSetupType(amg_data));
1452: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup type %s (%" PetscInt_FMT ")\n", HYPREILUIterSetup[indx], indx));
1453: PetscCallExternalVoid("hypre_ParAMGDataILUIterSetupOption", indx = hypre_ParAMGDataILUIterSetupOption(amg_data));
1454: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup option %" PetscInt_FMT "\n", indx));
1455: PetscCallExternalVoid("hypre_ParAMGDataILUIterSetupMaxIter", indx = hypre_ParAMGDataILUIterSetupMaxIter(amg_data));
1456: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup max iterations %" PetscInt_FMT "\n", indx));
1457: PetscCallExternalVoid("hypre_ParAMGDataILUIterSetupTolerance", val = hypre_ParAMGDataILUIterSetupTolerance(amg_data));
1458: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup tolerance %e\n", val));
1459: }
1460: } else {
1461: PetscCall(PetscViewerASCIIPrintf(viewer, " Not using more complex smoothers.\n"));
1462: }
1463: if (jac->smoothtype == 3) {
1464: PetscCall(PetscViewerASCIIPrintf(viewer, " Euclid ILU(k) levels %" PetscInt_FMT "\n", jac->eu_level));
1465: PetscCall(PetscViewerASCIIPrintf(viewer, " Euclid ILU(k) drop tolerance %g\n", (double)jac->eu_droptolerance));
1466: PetscCall(PetscViewerASCIIPrintf(viewer, " Euclid ILU use Block-Jacobi? %" PetscInt_FMT "\n", jac->eu_bj));
1467: }
1468: PetscCall(PetscViewerASCIIPrintf(viewer, " Measure type %s\n", HYPREBoomerAMGMeasureType[jac->measuretype]));
1469: PetscCall(PetscViewerASCIIPrintf(viewer, " Coarsen type %s\n", jac->coarsentype < 0 ? "not yet set" : HYPREBoomerAMGCoarsenType[jac->coarsentype]));
1470: PetscCall(PetscViewerASCIIPrintf(viewer, " Interpolation type %s\n", jac->interptype != 100 ? (jac->interptype < 0 ? "not yet set" : HYPREBoomerAMGInterpType[jac->interptype]) : "1pt"));
1471: if (jac->nodal_coarsening) PetscCall(PetscViewerASCIIPrintf(viewer, " Using nodal coarsening with HYPRE_BOOMERAMGSetNodal() %" PetscInt_FMT "\n", jac->nodal_coarsening));
1472: if (jac->vec_interp_variant) {
1473: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE_BoomerAMGSetInterpVecVariant() %" PetscInt_FMT "\n", jac->vec_interp_variant));
1474: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE_BoomerAMGSetInterpVecQMax() %" PetscInt_FMT "\n", jac->vec_interp_qmax));
1475: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE_BoomerAMGSetSmoothInterpVectors() %d\n", jac->vec_interp_smooth));
1476: }
1477: if (jac->nodal_relax) PetscCall(PetscViewerASCIIPrintf(viewer, " Using nodal relaxation via Schwarz smoothing on levels %" PetscInt_FMT "\n", jac->nodal_relax_levels));
1478: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
1479: PetscCall(PetscViewerASCIIPrintf(viewer, " SpGEMM type %s\n", jac->spgemm_type));
1480: #else
1481: PetscCall(PetscViewerASCIIPrintf(viewer, " SpGEMM type %s\n", "hypre"));
1482: #endif
1483: /* AIR */
1484: if (jac->Rtype) {
1485: PetscCall(PetscViewerASCIIPrintf(viewer, " Using approximate ideal restriction type %" PetscInt_FMT "\n", jac->Rtype));
1486: PetscCall(PetscViewerASCIIPrintf(viewer, " Threshold for R %g\n", (double)jac->Rstrongthreshold));
1487: PetscCall(PetscViewerASCIIPrintf(viewer, " Filter for R %g\n", (double)jac->Rfilterthreshold));
1488: PetscCall(PetscViewerASCIIPrintf(viewer, " A drop tolerance %g\n", (double)jac->Adroptol));
1489: PetscCall(PetscViewerASCIIPrintf(viewer, " A drop type %" PetscInt_FMT "\n", jac->Adroptype));
1490: }
1491: }
1492: PetscFunctionReturn(PETSC_SUCCESS);
1493: }
1495: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PC pc, PetscOptionItems PetscOptionsObject)
1496: {
1497: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1498: PetscInt indx;
1499: PetscBool flag, value;
1500: const char *symtlist[] = {"nonsymmetric", "SPD", "nonsymmetric,SPD"};
1502: PetscFunctionBegin;
1503: PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE ParaSails Options");
1504: PetscCall(PetscOptionsInt("-pc_hypre_parasails_nlevels", "Number of number of levels", "None", jac->nlevels, &jac->nlevels, 0));
1505: PetscCall(PetscOptionsReal("-pc_hypre_parasails_thresh", "Threshold", "None", jac->threshold, &jac->threshold, &flag));
1506: if (flag) PetscCallHYPRE(HYPRE_ParaSailsSetParams(jac->hsolver, jac->threshold, (HYPRE_Int)jac->nlevels));
1508: PetscCall(PetscOptionsReal("-pc_hypre_parasails_filter", "filter", "None", jac->filter, &jac->filter, &flag));
1509: if (flag) PetscCallHYPRE(HYPRE_ParaSailsSetFilter(jac->hsolver, jac->filter));
1511: PetscCall(PetscOptionsReal("-pc_hypre_parasails_loadbal", "Load balance", "None", jac->loadbal, &jac->loadbal, &flag));
1512: if (flag) PetscCallHYPRE(HYPRE_ParaSailsSetLoadbal(jac->hsolver, (HYPRE_Int)jac->loadbal));
1514: PetscCall(PetscOptionsBool("-pc_hypre_parasails_logging", "Print info to screen", "None", (PetscBool)jac->logging, &value, &flag));
1515: if (flag) {
1516: jac->logging = value ? 1 : 0;
1517: PetscCallHYPRE(HYPRE_ParaSailsSetLogging(jac->hsolver, (HYPRE_Int)jac->logging));
1518: }
1520: PetscCall(PetscOptionsBool("-pc_hypre_parasails_reuse", "Reuse nonzero pattern in preconditioner", "None", (PetscBool)jac->ruse, &value, &flag));
1521: if (flag) {
1522: jac->ruse = value ? 1 : 0;
1523: PetscCallHYPRE(HYPRE_ParaSailsSetReuse(jac->hsolver, (HYPRE_Int)jac->ruse));
1524: }
1526: PetscCall(PetscOptionsEList("-pc_hypre_parasails_sym", "Symmetry of matrix and preconditioner", "None", symtlist, PETSC_STATIC_ARRAY_LENGTH(symtlist), symtlist[0], &indx, &flag));
1527: if (flag) {
1528: jac->symt = indx;
1529: PetscCallHYPRE(HYPRE_ParaSailsSetSym(jac->hsolver, (HYPRE_Int)jac->symt));
1530: }
1532: PetscOptionsHeadEnd();
1533: PetscFunctionReturn(PETSC_SUCCESS);
1534: }
1536: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc, PetscViewer viewer)
1537: {
1538: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1539: PetscBool isascii;
1540: const char *symt = 0;
1542: PetscFunctionBegin;
1543: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1544: if (isascii) {
1545: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE ParaSails preconditioning\n"));
1546: PetscCall(PetscViewerASCIIPrintf(viewer, " nlevels %" PetscInt_FMT "\n", jac->nlevels));
1547: PetscCall(PetscViewerASCIIPrintf(viewer, " threshold %g\n", (double)jac->threshold));
1548: PetscCall(PetscViewerASCIIPrintf(viewer, " filter %g\n", (double)jac->filter));
1549: PetscCall(PetscViewerASCIIPrintf(viewer, " load balance %g\n", (double)jac->loadbal));
1550: PetscCall(PetscViewerASCIIPrintf(viewer, " reuse nonzero structure %s\n", PetscBools[jac->ruse]));
1551: PetscCall(PetscViewerASCIIPrintf(viewer, " print info to screen %s\n", PetscBools[jac->logging]));
1552: if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
1553: else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
1554: else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
1555: else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Unknown HYPRE ParaSails symmetric option %" PetscInt_FMT, jac->symt);
1556: PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", symt));
1557: }
1558: PetscFunctionReturn(PETSC_SUCCESS);
1559: }
1561: static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PC pc, PetscOptionItems PetscOptionsObject)
1562: {
1563: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1564: PetscInt n;
1565: PetscBool flag, flag2, flag3, flag4;
1567: PetscFunctionBegin;
1568: PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE AMS Options");
1569: PetscCall(PetscOptionsInt("-pc_hypre_ams_print_level", "Debugging output level for AMS", "None", jac->as_print, &jac->as_print, &flag));
1570: if (flag) PetscCallHYPRE(HYPRE_AMSSetPrintLevel(jac->hsolver, (HYPRE_Int)jac->as_print));
1571: 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));
1572: if (flag) PetscCallHYPRE(HYPRE_AMSSetMaxIter(jac->hsolver, (HYPRE_Int)jac->as_max_iter));
1573: PetscCall(PetscOptionsInt("-pc_hypre_ams_cycle_type", "Cycle type for AMS multigrid", "None", jac->ams_cycle_type, &jac->ams_cycle_type, &flag));
1574: if (flag) PetscCallHYPRE(HYPRE_AMSSetCycleType(jac->hsolver, (HYPRE_Int)jac->ams_cycle_type));
1575: PetscCall(PetscOptionsReal("-pc_hypre_ams_tol", "Error tolerance for AMS multigrid", "None", jac->as_tol, &jac->as_tol, &flag));
1576: if (flag) PetscCallHYPRE(HYPRE_AMSSetTol(jac->hsolver, jac->as_tol));
1577: PetscCall(PetscOptionsInt("-pc_hypre_ams_relax_type", "Relaxation type for AMS smoother", "None", jac->as_relax_type, &jac->as_relax_type, &flag));
1578: PetscCall(PetscOptionsInt("-pc_hypre_ams_relax_times", "Number of relaxation steps for AMS smoother", "None", jac->as_relax_times, &jac->as_relax_times, &flag2));
1579: PetscCall(PetscOptionsReal("-pc_hypre_ams_relax_weight", "Relaxation weight for AMS smoother", "None", jac->as_relax_weight, &jac->as_relax_weight, &flag3));
1580: PetscCall(PetscOptionsReal("-pc_hypre_ams_omega", "SSOR coefficient for AMS smoother", "None", jac->as_omega, &jac->as_omega, &flag4));
1581: if (flag || flag2 || flag3 || flag4) PetscCallHYPRE(HYPRE_AMSSetSmoothingOptions(jac->hsolver, (HYPRE_Int)jac->as_relax_type, (HYPRE_Int)jac->as_relax_times, jac->as_relax_weight, jac->as_omega));
1582: 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));
1583: n = 5;
1584: PetscCall(PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options", "AMG options for vector Poisson", "None", jac->as_amg_alpha_opts, &n, &flag2));
1585: if (flag || flag2) {
1586: PetscCallHYPRE(HYPRE_AMSSetAlphaAMGOptions(jac->hsolver, (HYPRE_Int)jac->as_amg_alpha_opts[0], /* AMG coarsen type */
1587: (HYPRE_Int)jac->as_amg_alpha_opts[1], /* AMG agg_levels */
1588: (HYPRE_Int)jac->as_amg_alpha_opts[2], /* AMG relax_type */
1589: jac->as_amg_alpha_theta, (HYPRE_Int)jac->as_amg_alpha_opts[3], /* AMG interp_type */
1590: (HYPRE_Int)jac->as_amg_alpha_opts[4])); /* AMG Pmax */
1591: }
1592: 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));
1593: n = 5;
1594: PetscCall(PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options", "AMG options for scalar Poisson solver", "None", jac->as_amg_beta_opts, &n, &flag2));
1595: if (flag || flag2) {
1596: PetscCallHYPRE(HYPRE_AMSSetBetaAMGOptions(jac->hsolver, (HYPRE_Int)jac->as_amg_beta_opts[0], /* AMG coarsen type */
1597: (HYPRE_Int)jac->as_amg_beta_opts[1], /* AMG agg_levels */
1598: (HYPRE_Int)jac->as_amg_beta_opts[2], /* AMG relax_type */
1599: jac->as_amg_beta_theta, (HYPRE_Int)jac->as_amg_beta_opts[3], /* AMG interp_type */
1600: (HYPRE_Int)jac->as_amg_beta_opts[4])); /* AMG Pmax */
1601: }
1602: 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));
1603: if (flag) { /* override HYPRE's default only if the options is used */
1604: PetscCallHYPRE(HYPRE_AMSSetProjectionFrequency(jac->hsolver, (HYPRE_Int)jac->ams_proj_freq));
1605: }
1606: PetscOptionsHeadEnd();
1607: PetscFunctionReturn(PETSC_SUCCESS);
1608: }
1610: static PetscErrorCode PCView_HYPRE_AMS(PC pc, PetscViewer viewer)
1611: {
1612: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1613: PetscBool isascii;
1615: PetscFunctionBegin;
1616: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1617: if (isascii) {
1618: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE AMS preconditioning\n"));
1619: PetscCall(PetscViewerASCIIPrintf(viewer, " subspace iterations per application %" PetscInt_FMT "\n", jac->as_max_iter));
1620: PetscCall(PetscViewerASCIIPrintf(viewer, " subspace cycle type %" PetscInt_FMT "\n", jac->ams_cycle_type));
1621: PetscCall(PetscViewerASCIIPrintf(viewer, " subspace iteration tolerance %g\n", (double)jac->as_tol));
1622: PetscCall(PetscViewerASCIIPrintf(viewer, " smoother type %" PetscInt_FMT "\n", jac->as_relax_type));
1623: PetscCall(PetscViewerASCIIPrintf(viewer, " number of smoothing steps %" PetscInt_FMT "\n", jac->as_relax_times));
1624: PetscCall(PetscViewerASCIIPrintf(viewer, " smoother weight %g\n", (double)jac->as_relax_weight));
1625: PetscCall(PetscViewerASCIIPrintf(viewer, " smoother omega %g\n", (double)jac->as_omega));
1626: if (jac->alpha_Poisson) {
1627: PetscCall(PetscViewerASCIIPrintf(viewer, " vector Poisson solver (passed in by user)\n"));
1628: } else {
1629: PetscCall(PetscViewerASCIIPrintf(viewer, " vector Poisson solver (computed) \n"));
1630: }
1631: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG coarsening type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[0]));
1632: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[1]));
1633: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG relaxation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[2]));
1634: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG interpolation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[3]));
1635: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[4]));
1636: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG strength threshold %g\n", (double)jac->as_amg_alpha_theta));
1637: if (!jac->ams_beta_is_zero) {
1638: if (jac->beta_Poisson) {
1639: PetscCall(PetscViewerASCIIPrintf(viewer, " scalar Poisson solver (passed in by user)\n"));
1640: } else {
1641: PetscCall(PetscViewerASCIIPrintf(viewer, " scalar Poisson solver (computed) \n"));
1642: }
1643: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG coarsening type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[0]));
1644: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_beta_opts[1]));
1645: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG relaxation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[2]));
1646: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG interpolation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[3]));
1647: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_beta_opts[4]));
1648: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG strength threshold %g\n", (double)jac->as_amg_beta_theta));
1649: 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));
1650: } else {
1651: PetscCall(PetscViewerASCIIPrintf(viewer, " scalar Poisson solver not used (zero-conductivity everywhere) \n"));
1652: }
1653: }
1654: PetscFunctionReturn(PETSC_SUCCESS);
1655: }
1657: static PetscErrorCode PCSetFromOptions_HYPRE_ADS(PC pc, PetscOptionItems PetscOptionsObject)
1658: {
1659: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1660: PetscInt n;
1661: PetscBool flag, flag2, flag3, flag4;
1663: PetscFunctionBegin;
1664: PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE ADS Options");
1665: PetscCall(PetscOptionsInt("-pc_hypre_ads_print_level", "Debugging output level for ADS", "None", jac->as_print, &jac->as_print, &flag));
1666: if (flag) PetscCallHYPRE(HYPRE_ADSSetPrintLevel(jac->hsolver, (HYPRE_Int)jac->as_print));
1667: 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));
1668: if (flag) PetscCallHYPRE(HYPRE_ADSSetMaxIter(jac->hsolver, (HYPRE_Int)jac->as_max_iter));
1669: PetscCall(PetscOptionsInt("-pc_hypre_ads_cycle_type", "Cycle type for ADS multigrid", "None", jac->ads_cycle_type, &jac->ads_cycle_type, &flag));
1670: if (flag) PetscCallHYPRE(HYPRE_ADSSetCycleType(jac->hsolver, (HYPRE_Int)jac->ads_cycle_type));
1671: PetscCall(PetscOptionsReal("-pc_hypre_ads_tol", "Error tolerance for ADS multigrid", "None", jac->as_tol, &jac->as_tol, &flag));
1672: if (flag) PetscCallHYPRE(HYPRE_ADSSetTol(jac->hsolver, jac->as_tol));
1673: PetscCall(PetscOptionsInt("-pc_hypre_ads_relax_type", "Relaxation type for ADS smoother", "None", jac->as_relax_type, &jac->as_relax_type, &flag));
1674: PetscCall(PetscOptionsInt("-pc_hypre_ads_relax_times", "Number of relaxation steps for ADS smoother", "None", jac->as_relax_times, &jac->as_relax_times, &flag2));
1675: PetscCall(PetscOptionsReal("-pc_hypre_ads_relax_weight", "Relaxation weight for ADS smoother", "None", jac->as_relax_weight, &jac->as_relax_weight, &flag3));
1676: PetscCall(PetscOptionsReal("-pc_hypre_ads_omega", "SSOR coefficient for ADS smoother", "None", jac->as_omega, &jac->as_omega, &flag4));
1677: if (flag || flag2 || flag3 || flag4) PetscCallHYPRE(HYPRE_ADSSetSmoothingOptions(jac->hsolver, (HYPRE_Int)jac->as_relax_type, (HYPRE_Int)jac->as_relax_times, jac->as_relax_weight, jac->as_omega));
1678: 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));
1679: n = 5;
1680: PetscCall(PetscOptionsIntArray("-pc_hypre_ads_ams_options", "AMG options for AMS solver inside ADS", "None", jac->as_amg_alpha_opts, &n, &flag2));
1681: 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));
1682: if (flag || flag2 || flag3) {
1683: PetscCallHYPRE(HYPRE_ADSSetAMSOptions(jac->hsolver, (HYPRE_Int)jac->ams_cycle_type, /* AMS cycle type */
1684: (HYPRE_Int)jac->as_amg_alpha_opts[0], /* AMG coarsen type */
1685: (HYPRE_Int)jac->as_amg_alpha_opts[1], /* AMG agg_levels */
1686: (HYPRE_Int)jac->as_amg_alpha_opts[2], /* AMG relax_type */
1687: jac->as_amg_alpha_theta, (HYPRE_Int)jac->as_amg_alpha_opts[3], /* AMG interp_type */
1688: (HYPRE_Int)jac->as_amg_alpha_opts[4])); /* AMG Pmax */
1689: }
1690: 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));
1691: n = 5;
1692: PetscCall(PetscOptionsIntArray("-pc_hypre_ads_amg_options", "AMG options for vector AMG solver inside ADS", "None", jac->as_amg_beta_opts, &n, &flag2));
1693: if (flag || flag2) {
1694: PetscCallHYPRE(HYPRE_ADSSetAMGOptions(jac->hsolver, (HYPRE_Int)jac->as_amg_beta_opts[0], /* AMG coarsen type */
1695: (HYPRE_Int)jac->as_amg_beta_opts[1], /* AMG agg_levels */
1696: (HYPRE_Int)jac->as_amg_beta_opts[2], /* AMG relax_type */
1697: jac->as_amg_beta_theta, (HYPRE_Int)jac->as_amg_beta_opts[3], /* AMG interp_type */
1698: (HYPRE_Int)jac->as_amg_beta_opts[4])); /* AMG Pmax */
1699: }
1700: PetscOptionsHeadEnd();
1701: PetscFunctionReturn(PETSC_SUCCESS);
1702: }
1704: static PetscErrorCode PCView_HYPRE_ADS(PC pc, PetscViewer viewer)
1705: {
1706: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1707: PetscBool isascii;
1709: PetscFunctionBegin;
1710: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1711: if (isascii) {
1712: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE ADS preconditioning\n"));
1713: PetscCall(PetscViewerASCIIPrintf(viewer, " subspace iterations per application %" PetscInt_FMT "\n", jac->as_max_iter));
1714: PetscCall(PetscViewerASCIIPrintf(viewer, " subspace cycle type %" PetscInt_FMT "\n", jac->ads_cycle_type));
1715: PetscCall(PetscViewerASCIIPrintf(viewer, " subspace iteration tolerance %g\n", (double)jac->as_tol));
1716: PetscCall(PetscViewerASCIIPrintf(viewer, " smoother type %" PetscInt_FMT "\n", jac->as_relax_type));
1717: PetscCall(PetscViewerASCIIPrintf(viewer, " number of smoothing steps %" PetscInt_FMT "\n", jac->as_relax_times));
1718: PetscCall(PetscViewerASCIIPrintf(viewer, " smoother weight %g\n", (double)jac->as_relax_weight));
1719: PetscCall(PetscViewerASCIIPrintf(viewer, " smoother omega %g\n", (double)jac->as_omega));
1720: PetscCall(PetscViewerASCIIPrintf(viewer, " AMS solver using boomerAMG\n"));
1721: PetscCall(PetscViewerASCIIPrintf(viewer, " subspace cycle type %" PetscInt_FMT "\n", jac->ams_cycle_type));
1722: PetscCall(PetscViewerASCIIPrintf(viewer, " coarsening type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[0]));
1723: PetscCall(PetscViewerASCIIPrintf(viewer, " levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[1]));
1724: PetscCall(PetscViewerASCIIPrintf(viewer, " relaxation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[2]));
1725: PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[3]));
1726: PetscCall(PetscViewerASCIIPrintf(viewer, " max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[4]));
1727: PetscCall(PetscViewerASCIIPrintf(viewer, " strength threshold %g\n", (double)jac->as_amg_alpha_theta));
1728: PetscCall(PetscViewerASCIIPrintf(viewer, " vector Poisson solver using boomerAMG\n"));
1729: PetscCall(PetscViewerASCIIPrintf(viewer, " coarsening type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[0]));
1730: PetscCall(PetscViewerASCIIPrintf(viewer, " levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_beta_opts[1]));
1731: PetscCall(PetscViewerASCIIPrintf(viewer, " relaxation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[2]));
1732: PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[3]));
1733: PetscCall(PetscViewerASCIIPrintf(viewer, " max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_beta_opts[4]));
1734: PetscCall(PetscViewerASCIIPrintf(viewer, " strength threshold %g\n", (double)jac->as_amg_beta_theta));
1735: }
1736: PetscFunctionReturn(PETSC_SUCCESS);
1737: }
1739: static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE(PC pc, Mat G)
1740: {
1741: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1742: PetscBool ishypre;
1744: PetscFunctionBegin;
1745: PetscCall(PetscObjectTypeCompare((PetscObject)G, MATHYPRE, &ishypre));
1746: if (ishypre) {
1747: PetscCall(PetscObjectReference((PetscObject)G));
1748: PetscCall(MatDestroy(&jac->G));
1749: jac->G = G;
1750: } else {
1751: PetscCall(MatDestroy(&jac->G));
1752: PetscCall(MatConvert(G, MATHYPRE, MAT_INITIAL_MATRIX, &jac->G));
1753: }
1754: PetscFunctionReturn(PETSC_SUCCESS);
1755: }
1757: static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C)
1758: {
1759: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1760: PetscBool ishypre;
1762: PetscFunctionBegin;
1763: PetscCall(PetscObjectTypeCompare((PetscObject)C, MATHYPRE, &ishypre));
1764: if (ishypre) {
1765: PetscCall(PetscObjectReference((PetscObject)C));
1766: PetscCall(MatDestroy(&jac->C));
1767: jac->C = C;
1768: } else {
1769: PetscCall(MatDestroy(&jac->C));
1770: PetscCall(MatConvert(C, MATHYPRE, MAT_INITIAL_MATRIX, &jac->C));
1771: }
1772: PetscFunctionReturn(PETSC_SUCCESS);
1773: }
1775: static PetscErrorCode PCHYPRESetInterpolations_HYPRE(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[])
1776: {
1777: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1778: PetscBool ishypre;
1780: PetscFunctionBegin;
1781: PetscCall(MatDestroy(&jac->RT_PiFull));
1782: PetscCall(MatDestroy(&jac->ND_PiFull));
1783: for (PetscInt i = 0; i < 3; ++i) {
1784: PetscCall(MatDestroy(&jac->RT_Pi[i]));
1785: PetscCall(MatDestroy(&jac->ND_Pi[i]));
1786: }
1788: jac->dim = dim;
1789: if (RT_PiFull) {
1790: PetscCall(PetscObjectTypeCompare((PetscObject)RT_PiFull, MATHYPRE, &ishypre));
1791: if (ishypre) {
1792: PetscCall(PetscObjectReference((PetscObject)RT_PiFull));
1793: jac->RT_PiFull = RT_PiFull;
1794: } else {
1795: PetscCall(MatConvert(RT_PiFull, MATHYPRE, MAT_INITIAL_MATRIX, &jac->RT_PiFull));
1796: }
1797: }
1798: if (RT_Pi) {
1799: for (PetscInt i = 0; i < dim; ++i) {
1800: if (RT_Pi[i]) {
1801: PetscCall(PetscObjectTypeCompare((PetscObject)RT_Pi[i], MATHYPRE, &ishypre));
1802: if (ishypre) {
1803: PetscCall(PetscObjectReference((PetscObject)RT_Pi[i]));
1804: jac->RT_Pi[i] = RT_Pi[i];
1805: } else {
1806: PetscCall(MatConvert(RT_Pi[i], MATHYPRE, MAT_INITIAL_MATRIX, &jac->RT_Pi[i]));
1807: }
1808: }
1809: }
1810: }
1811: if (ND_PiFull) {
1812: PetscCall(PetscObjectTypeCompare((PetscObject)ND_PiFull, MATHYPRE, &ishypre));
1813: if (ishypre) {
1814: PetscCall(PetscObjectReference((PetscObject)ND_PiFull));
1815: jac->ND_PiFull = ND_PiFull;
1816: } else {
1817: PetscCall(MatConvert(ND_PiFull, MATHYPRE, MAT_INITIAL_MATRIX, &jac->ND_PiFull));
1818: }
1819: }
1820: if (ND_Pi) {
1821: for (PetscInt i = 0; i < dim; ++i) {
1822: if (ND_Pi[i]) {
1823: PetscCall(PetscObjectTypeCompare((PetscObject)ND_Pi[i], MATHYPRE, &ishypre));
1824: if (ishypre) {
1825: PetscCall(PetscObjectReference((PetscObject)ND_Pi[i]));
1826: jac->ND_Pi[i] = ND_Pi[i];
1827: } else {
1828: PetscCall(MatConvert(ND_Pi[i], MATHYPRE, MAT_INITIAL_MATRIX, &jac->ND_Pi[i]));
1829: }
1830: }
1831: }
1832: }
1833: PetscFunctionReturn(PETSC_SUCCESS);
1834: }
1836: static PetscErrorCode PCHYPRESetPoissonMatrix_HYPRE(PC pc, Mat A, PetscBool isalpha)
1837: {
1838: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1839: PetscBool ishypre;
1841: PetscFunctionBegin;
1842: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre));
1843: if (ishypre) {
1844: if (isalpha) {
1845: PetscCall(PetscObjectReference((PetscObject)A));
1846: PetscCall(MatDestroy(&jac->alpha_Poisson));
1847: jac->alpha_Poisson = A;
1848: } else {
1849: if (A) {
1850: PetscCall(PetscObjectReference((PetscObject)A));
1851: } else {
1852: jac->ams_beta_is_zero = PETSC_TRUE;
1853: }
1854: PetscCall(MatDestroy(&jac->beta_Poisson));
1855: jac->beta_Poisson = A;
1856: }
1857: } else {
1858: if (isalpha) {
1859: PetscCall(MatDestroy(&jac->alpha_Poisson));
1860: PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &jac->alpha_Poisson));
1861: } else {
1862: if (A) {
1863: PetscCall(MatDestroy(&jac->beta_Poisson));
1864: PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &jac->beta_Poisson));
1865: } else {
1866: PetscCall(MatDestroy(&jac->beta_Poisson));
1867: jac->ams_beta_is_zero = PETSC_TRUE;
1868: }
1869: }
1870: }
1871: PetscFunctionReturn(PETSC_SUCCESS);
1872: }
1874: static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE(PC pc, Vec ozz, Vec zoz, Vec zzo)
1875: {
1876: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1878: PetscFunctionBegin;
1879: /* throw away any vector if already set */
1880: PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[0]));
1881: PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[1]));
1882: PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[2]));
1883: PetscCall(VecHYPRE_IJVectorCreate(ozz->map, &jac->constants[0]));
1884: PetscCall(VecHYPRE_IJVectorCopy(ozz, jac->constants[0]));
1885: PetscCall(VecHYPRE_IJVectorCreate(zoz->map, &jac->constants[1]));
1886: PetscCall(VecHYPRE_IJVectorCopy(zoz, jac->constants[1]));
1887: jac->dim = 2;
1888: if (zzo) {
1889: PetscCall(VecHYPRE_IJVectorCreate(zzo->map, &jac->constants[2]));
1890: PetscCall(VecHYPRE_IJVectorCopy(zzo, jac->constants[2]));
1891: jac->dim++;
1892: }
1893: PetscFunctionReturn(PETSC_SUCCESS);
1894: }
1896: static PetscErrorCode PCHYPREAMSSetInteriorNodes_HYPRE(PC pc, Vec interior)
1897: {
1898: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1900: PetscFunctionBegin;
1901: PetscCall(VecHYPRE_IJVectorDestroy(&jac->interior));
1902: PetscCall(VecHYPRE_IJVectorCreate(interior->map, &jac->interior));
1903: PetscCall(VecHYPRE_IJVectorCopy(interior, jac->interior));
1904: jac->ams_beta_is_zero_part = PETSC_TRUE;
1905: PetscFunctionReturn(PETSC_SUCCESS);
1906: }
1908: static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1909: {
1910: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1911: Vec tv;
1913: PetscFunctionBegin;
1914: /* throw away any coordinate vector if already set */
1915: PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[0]));
1916: PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[1]));
1917: PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[2]));
1918: jac->dim = dim;
1920: /* compute IJ vector for coordinates */
1921: PetscCall(VecCreate(PetscObjectComm((PetscObject)pc), &tv));
1922: PetscCall(VecSetType(tv, VECSTANDARD));
1923: PetscCall(VecSetSizes(tv, nloc, PETSC_DECIDE));
1924: for (PetscInt i = 0; i < dim; i++) {
1925: PetscScalar *array;
1927: PetscCall(VecHYPRE_IJVectorCreate(tv->map, &jac->coords[i]));
1928: PetscCall(VecGetArrayWrite(tv, &array));
1929: for (PetscInt j = 0; j < nloc; j++) array[j] = coords[j * dim + i];
1930: PetscCall(VecRestoreArrayWrite(tv, &array));
1931: PetscCall(VecHYPRE_IJVectorCopy(tv, jac->coords[i]));
1932: }
1933: PetscCall(VecDestroy(&tv));
1934: PetscFunctionReturn(PETSC_SUCCESS);
1935: }
1937: static PetscErrorCode PCHYPREGetType_HYPRE(PC pc, const char *name[])
1938: {
1939: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1941: PetscFunctionBegin;
1942: *name = jac->hypre_type;
1943: PetscFunctionReturn(PETSC_SUCCESS);
1944: }
1946: static PetscErrorCode PCHYPRESetType_HYPRE(PC pc, const char name[])
1947: {
1948: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1949: PetscBool flag;
1951: PetscFunctionBegin;
1952: if (jac->hypre_type) {
1953: PetscCall(PetscStrcmp(jac->hypre_type, name, &flag));
1954: if (flag) PetscFunctionReturn(PETSC_SUCCESS);
1955: }
1957: PetscCall(PCReset_HYPRE(pc));
1958: PetscCall(PetscFree(jac->hypre_type));
1959: PetscCall(PetscStrallocpy(name, &jac->hypre_type));
1961: jac->maxiter = PETSC_DEFAULT;
1962: jac->tol = PETSC_DEFAULT;
1963: jac->printstatistics = PetscLogPrintInfo;
1965: PetscCall(PetscStrcmp("ilu", jac->hypre_type, &flag));
1966: if (flag) {
1967: PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
1968: PetscCallHYPRE(HYPRE_ILUCreate(&jac->hsolver));
1969: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ILU;
1970: pc->ops->view = PCView_HYPRE_ILU;
1971: jac->destroy = HYPRE_ILUDestroy;
1972: jac->setup = HYPRE_ILUSetup;
1973: jac->solve = HYPRE_ILUSolve;
1974: jac->factorrowsize = PETSC_DEFAULT;
1975: PetscFunctionReturn(PETSC_SUCCESS);
1976: }
1978: PetscCall(PetscStrcmp("pilut", jac->hypre_type, &flag));
1979: if (flag) {
1980: PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
1981: PetscCallHYPRE(HYPRE_ParCSRPilutCreate(jac->comm_hypre, &jac->hsolver));
1982: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
1983: pc->ops->view = PCView_HYPRE_Pilut;
1984: jac->destroy = HYPRE_ParCSRPilutDestroy;
1985: jac->setup = HYPRE_ParCSRPilutSetup;
1986: jac->solve = HYPRE_ParCSRPilutSolve;
1987: jac->factorrowsize = PETSC_DEFAULT;
1988: PetscFunctionReturn(PETSC_SUCCESS);
1989: }
1990: PetscCall(PetscStrcmp("euclid", jac->hypre_type, &flag));
1991: if (flag) {
1992: PetscCheck(!PetscDefined(USE_64BIT_INDICES), PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Hypre Euclid does not support 64-bit indices");
1993: PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
1994: PetscCallHYPRE(HYPRE_EuclidCreate(jac->comm_hypre, &jac->hsolver));
1995: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
1996: pc->ops->view = PCView_HYPRE_Euclid;
1997: jac->destroy = HYPRE_EuclidDestroy;
1998: jac->setup = HYPRE_EuclidSetup;
1999: jac->solve = HYPRE_EuclidSolve;
2000: jac->factorrowsize = PETSC_DEFAULT;
2001: jac->eu_level = PETSC_DEFAULT; /* default */
2002: PetscFunctionReturn(PETSC_SUCCESS);
2003: }
2004: PetscCall(PetscStrcmp("parasails", jac->hypre_type, &flag));
2005: if (flag) {
2006: PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
2007: PetscCallHYPRE(HYPRE_ParaSailsCreate(jac->comm_hypre, &jac->hsolver));
2008: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
2009: pc->ops->view = PCView_HYPRE_ParaSails;
2010: jac->destroy = HYPRE_ParaSailsDestroy;
2011: jac->setup = HYPRE_ParaSailsSetup;
2012: jac->solve = HYPRE_ParaSailsSolve;
2013: /* initialize */
2014: jac->nlevels = 1;
2015: jac->threshold = .1;
2016: jac->filter = .1;
2017: jac->loadbal = 0;
2018: if (PetscLogPrintInfo) jac->logging = (int)PETSC_TRUE;
2019: else jac->logging = (int)PETSC_FALSE;
2021: jac->ruse = (int)PETSC_FALSE;
2022: jac->symt = 0;
2023: PetscCallHYPRE(HYPRE_ParaSailsSetParams(jac->hsolver, jac->threshold, (HYPRE_Int)jac->nlevels));
2024: PetscCallHYPRE(HYPRE_ParaSailsSetFilter(jac->hsolver, jac->filter));
2025: PetscCallHYPRE(HYPRE_ParaSailsSetLoadbal(jac->hsolver, (HYPRE_Int)jac->loadbal));
2026: PetscCallHYPRE(HYPRE_ParaSailsSetLogging(jac->hsolver, (HYPRE_Int)jac->logging));
2027: PetscCallHYPRE(HYPRE_ParaSailsSetReuse(jac->hsolver, (HYPRE_Int)jac->ruse));
2028: PetscCallHYPRE(HYPRE_ParaSailsSetSym(jac->hsolver, (HYPRE_Int)jac->symt));
2029: PetscFunctionReturn(PETSC_SUCCESS);
2030: }
2031: PetscCall(PetscStrcmp("boomeramg", jac->hypre_type, &flag));
2032: if (flag) {
2033: PetscCallHYPRE(HYPRE_BoomerAMGCreate(&jac->hsolver));
2034: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_BoomerAMG;
2035: pc->ops->view = PCView_HYPRE_BoomerAMG;
2036: pc->ops->applytranspose = PCApplyTranspose_HYPRE_BoomerAMG;
2037: pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
2038: pc->ops->matapply = PCMatApply_HYPRE_BoomerAMG;
2039: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", PCGetInterpolations_BoomerAMG));
2040: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", PCGetCoarseOperators_BoomerAMG));
2041: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetCFMarkers_C", PCHYPREGetCFMarkers_BoomerAMG));
2042: jac->destroy = HYPRE_BoomerAMGDestroy;
2043: jac->setup = HYPRE_BoomerAMGSetup;
2044: jac->solve = HYPRE_BoomerAMGSolve;
2045: jac->applyrichardson = PETSC_FALSE;
2046: /* these defaults match the hypre defaults */
2047: jac->cycletype = 1;
2048: jac->maxlevels = 25;
2049: jac->maxiter = 1;
2050: jac->tol = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
2051: jac->truncfactor = 0.0;
2052: jac->strongthreshold = .25;
2053: jac->maxrowsum = .9;
2054: jac->measuretype = 0;
2055: jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
2056: jac->smoothtype = -1; /* Not set by default */
2057: jac->smoothnumlevels = 25;
2058: jac->eu_level = 0;
2059: jac->eu_droptolerance = 0;
2060: jac->eu_bj = 0;
2061: jac->relaxweight = 1.0;
2062: jac->outerrelaxweight = 1.0;
2063: jac->Rtype = 0;
2064: jac->Rstrongthreshold = 0.25;
2065: jac->Rfilterthreshold = 0.0;
2066: jac->Adroptype = -1;
2067: jac->Adroptol = 0.0;
2068: jac->agg_nl = 0;
2069: jac->pmax = 0;
2070: jac->truncfactor = 0.0;
2071: jac->agg_num_paths = 1;
2072: jac->maxc = 9;
2073: jac->minc = 1;
2074: jac->nodal_coarsening = 0;
2075: jac->nodal_coarsening_diag = 0;
2076: jac->vec_interp_variant = 0;
2077: jac->vec_interp_qmax = 0;
2078: jac->vec_interp_smooth = PETSC_FALSE;
2079: jac->interp_refine = 0;
2080: jac->nodal_relax = PETSC_FALSE;
2081: jac->nodal_relax_levels = 1;
2082: jac->rap2 = 0;
2083: PetscObjectParameterSetDefault(jac, relaxtype[2], 9); /* G.E. */
2085: /*
2086: Initialize the following parameters with invalid value so we can recognize user input that sets the parameter.
2087: If there is no user input they are overwritten in PCSetUp_HYPRE() depending on if the matrix is on the CPU or the GPU
2088: */
2089: PetscObjectParameterSetDefault(jac, relaxorder, PETSC_DECIDE);
2090: PetscObjectParameterSetDefault(jac, coarsentype, PETSC_DECIDE);
2091: PetscObjectParameterSetDefault(jac, interptype, PETSC_DECIDE);
2092: PetscObjectParameterSetDefault(jac, relaxtype[0], PETSC_DECIDE);
2093: PetscObjectParameterSetDefault(jac, relaxtype[1], PETSC_DECIDE);
2094: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
2095: PetscObjectParameterSetDefault(jac, spgemm_type, "not yet set");
2096: #endif
2097: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
2098: PetscObjectParameterSetDefault(jac, keeptranspose, PETSC_BOOL3_UNKNOWN);
2099: PetscObjectParameterSetDefault(jac, mod_rap2, PETSC_DECIDE);
2100: #endif
2101: PetscObjectParameterSetDefault(jac, agg_interptype, PETSC_DECIDE);
2102: PetscFunctionReturn(PETSC_SUCCESS);
2103: }
2104: PetscCall(PetscStrcmp("ams", jac->hypre_type, &flag));
2105: if (flag) {
2106: PetscCallHYPRE(HYPRE_AMSCreate(&jac->hsolver));
2107: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_AMS;
2108: pc->ops->view = PCView_HYPRE_AMS;
2109: jac->destroy = HYPRE_AMSDestroy;
2110: jac->setup = HYPRE_AMSSetup;
2111: jac->solve = HYPRE_AMSSolve;
2112: jac->coords[0] = NULL;
2113: jac->coords[1] = NULL;
2114: jac->coords[2] = NULL;
2115: jac->interior = NULL;
2116: /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */
2117: jac->as_print = 0;
2118: jac->as_max_iter = 1; /* used as a preconditioner */
2119: jac->as_tol = 0.; /* used as a preconditioner */
2120: jac->ams_cycle_type = 13;
2121: /* Smoothing options */
2122: jac->as_relax_type = 2;
2123: jac->as_relax_times = 1;
2124: jac->as_relax_weight = 1.0;
2125: jac->as_omega = 1.0;
2126: /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
2127: jac->as_amg_alpha_opts[0] = 10;
2128: jac->as_amg_alpha_opts[1] = 1;
2129: jac->as_amg_alpha_opts[2] = 6;
2130: jac->as_amg_alpha_opts[3] = 6;
2131: jac->as_amg_alpha_opts[4] = 4;
2132: jac->as_amg_alpha_theta = 0.25;
2133: /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
2134: jac->as_amg_beta_opts[0] = 10;
2135: jac->as_amg_beta_opts[1] = 1;
2136: jac->as_amg_beta_opts[2] = 6;
2137: jac->as_amg_beta_opts[3] = 6;
2138: jac->as_amg_beta_opts[4] = 4;
2139: jac->as_amg_beta_theta = 0.25;
2140: PetscCallHYPRE(HYPRE_AMSSetPrintLevel(jac->hsolver, (HYPRE_Int)jac->as_print));
2141: PetscCallHYPRE(HYPRE_AMSSetMaxIter(jac->hsolver, (HYPRE_Int)jac->as_max_iter));
2142: PetscCallHYPRE(HYPRE_AMSSetCycleType(jac->hsolver, (HYPRE_Int)jac->ams_cycle_type));
2143: PetscCallHYPRE(HYPRE_AMSSetTol(jac->hsolver, jac->as_tol));
2144: PetscCallHYPRE(HYPRE_AMSSetSmoothingOptions(jac->hsolver, (HYPRE_Int)jac->as_relax_type, (HYPRE_Int)jac->as_relax_times, jac->as_relax_weight, jac->as_omega));
2145: PetscCallHYPRE(HYPRE_AMSSetAlphaAMGOptions(jac->hsolver, (HYPRE_Int)jac->as_amg_alpha_opts[0], /* AMG coarsen type */
2146: (HYPRE_Int)jac->as_amg_alpha_opts[1], /* AMG agg_levels */
2147: (HYPRE_Int)jac->as_amg_alpha_opts[2], /* AMG relax_type */
2148: jac->as_amg_alpha_theta, (HYPRE_Int)jac->as_amg_alpha_opts[3], /* AMG interp_type */
2149: (HYPRE_Int)jac->as_amg_alpha_opts[4])); /* AMG Pmax */
2150: PetscCallHYPRE(HYPRE_AMSSetBetaAMGOptions(jac->hsolver, (HYPRE_Int)jac->as_amg_beta_opts[0], /* AMG coarsen type */
2151: (HYPRE_Int)jac->as_amg_beta_opts[1], /* AMG agg_levels */
2152: (HYPRE_Int)jac->as_amg_beta_opts[2], /* AMG relax_type */
2153: jac->as_amg_beta_theta, (HYPRE_Int)jac->as_amg_beta_opts[3], /* AMG interp_type */
2154: (HYPRE_Int)jac->as_amg_beta_opts[4])); /* AMG Pmax */
2155: /* Zero conductivity */
2156: jac->ams_beta_is_zero = PETSC_FALSE;
2157: jac->ams_beta_is_zero_part = PETSC_FALSE;
2158: PetscFunctionReturn(PETSC_SUCCESS);
2159: }
2160: PetscCall(PetscStrcmp("ads", jac->hypre_type, &flag));
2161: if (flag) {
2162: PetscCallHYPRE(HYPRE_ADSCreate(&jac->hsolver));
2163: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ADS;
2164: pc->ops->view = PCView_HYPRE_ADS;
2165: jac->destroy = HYPRE_ADSDestroy;
2166: jac->setup = HYPRE_ADSSetup;
2167: jac->solve = HYPRE_ADSSolve;
2168: jac->coords[0] = NULL;
2169: jac->coords[1] = NULL;
2170: jac->coords[2] = NULL;
2171: /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */
2172: jac->as_print = 0;
2173: jac->as_max_iter = 1; /* used as a preconditioner */
2174: jac->as_tol = 0.; /* used as a preconditioner */
2175: jac->ads_cycle_type = 13;
2176: /* Smoothing options */
2177: jac->as_relax_type = 2;
2178: jac->as_relax_times = 1;
2179: jac->as_relax_weight = 1.0;
2180: jac->as_omega = 1.0;
2181: /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */
2182: jac->ams_cycle_type = 14;
2183: jac->as_amg_alpha_opts[0] = 10;
2184: jac->as_amg_alpha_opts[1] = 1;
2185: jac->as_amg_alpha_opts[2] = 6;
2186: jac->as_amg_alpha_opts[3] = 6;
2187: jac->as_amg_alpha_opts[4] = 4;
2188: jac->as_amg_alpha_theta = 0.25;
2189: /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
2190: jac->as_amg_beta_opts[0] = 10;
2191: jac->as_amg_beta_opts[1] = 1;
2192: jac->as_amg_beta_opts[2] = 6;
2193: jac->as_amg_beta_opts[3] = 6;
2194: jac->as_amg_beta_opts[4] = 4;
2195: jac->as_amg_beta_theta = 0.25;
2196: PetscCallHYPRE(HYPRE_ADSSetPrintLevel(jac->hsolver, (HYPRE_Int)jac->as_print));
2197: PetscCallHYPRE(HYPRE_ADSSetMaxIter(jac->hsolver, (HYPRE_Int)jac->as_max_iter));
2198: PetscCallHYPRE(HYPRE_ADSSetCycleType(jac->hsolver, (HYPRE_Int)jac->ams_cycle_type));
2199: PetscCallHYPRE(HYPRE_ADSSetTol(jac->hsolver, jac->as_tol));
2200: PetscCallHYPRE(HYPRE_ADSSetSmoothingOptions(jac->hsolver, (HYPRE_Int)jac->as_relax_type, (HYPRE_Int)jac->as_relax_times, jac->as_relax_weight, jac->as_omega));
2201: PetscCallHYPRE(HYPRE_ADSSetAMSOptions(jac->hsolver, (HYPRE_Int)jac->ams_cycle_type, /* AMG coarsen type */
2202: (HYPRE_Int)jac->as_amg_alpha_opts[0], /* AMG coarsen type */
2203: (HYPRE_Int)jac->as_amg_alpha_opts[1], /* AMG agg_levels */
2204: (HYPRE_Int)jac->as_amg_alpha_opts[2], /* AMG relax_type */
2205: jac->as_amg_alpha_theta, (HYPRE_Int)jac->as_amg_alpha_opts[3], /* AMG interp_type */
2206: (HYPRE_Int)jac->as_amg_alpha_opts[4])); /* AMG Pmax */
2207: PetscCallHYPRE(HYPRE_ADSSetAMGOptions(jac->hsolver, (HYPRE_Int)jac->as_amg_beta_opts[0], /* AMG coarsen type */
2208: (HYPRE_Int)jac->as_amg_beta_opts[1], /* AMG agg_levels */
2209: (HYPRE_Int)jac->as_amg_beta_opts[2], /* AMG relax_type */
2210: jac->as_amg_beta_theta, (HYPRE_Int)jac->as_amg_beta_opts[3], /* AMG interp_type */
2211: (HYPRE_Int)jac->as_amg_beta_opts[4])); /* AMG Pmax */
2212: PetscFunctionReturn(PETSC_SUCCESS);
2213: }
2214: PetscCall(PetscFree(jac->hypre_type));
2216: jac->hypre_type = NULL;
2217: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown HYPRE preconditioner %s; Choices are euclid, ilu, pilut, parasails, boomeramg, ams, ads", name);
2218: }
2220: /*
2221: It only gets here if the HYPRE type has not been set before the call to
2222: ...SetFromOptions() which actually is most of the time
2223: */
2224: static PetscErrorCode PCSetFromOptions_HYPRE(PC pc, PetscOptionItems PetscOptionsObject)
2225: {
2226: PetscInt indx;
2227: const char *type[] = {"ilu", "euclid", "pilut", "parasails", "boomeramg", "ams", "ads"};
2228: PetscBool flg;
2229: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
2231: PetscFunctionBegin;
2232: PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE preconditioner options");
2233: PetscCall(PetscOptionsEList("-pc_hypre_type", "HYPRE preconditioner type", "PCHYPRESetType", type, PETSC_STATIC_ARRAY_LENGTH(type), "boomeramg", &indx, &flg));
2234: if (flg) PetscCall(PCHYPRESetType_HYPRE(pc, type[indx]));
2235: /*
2236: Set the type if it was never set.
2237: */
2238: if (!jac->hypre_type) PetscCall(PCHYPRESetType_HYPRE(pc, "boomeramg"));
2239: PetscTryTypeMethod(pc, setfromoptions, PetscOptionsObject);
2240: PetscOptionsHeadEnd();
2241: PetscFunctionReturn(PETSC_SUCCESS);
2242: }
2244: /*MC
2245: PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre as PETSc `PC`
2247: Options Database Keys:
2248: + -pc_hypre_type - One of `euclid`, `ilu`, `pilut`, `parasails`, `boomeramg`, `ams`, or `ads`
2249: . -pc_hypre_boomeramg_nodal_coarsen n - where `n` is from 1 to 6 (see `HYPRE_BoomerAMGSetNodal()`)
2250: . -pc_hypre_boomeramg_vec_interp_variant v - where `v` is from 1 to 3 (see `HYPRE_BoomerAMGSetInterpVecVariant()`)
2251: - Many others - run with `-pc_type hypre` `-pc_hypre_type XXX` `-help` to see options for the XXX preconditioner
2253: Level: intermediate
2255: Notes:
2256: Apart from `-pc_hypre_type` (for which there is `PCHYPRESetType()`),
2257: the many hypre options can ONLY be set via the options database (e.g. the command line
2258: or with `PetscOptionsSetValue()`, there are no functions to set them)
2260: The options `-pc_hypre_boomeramg_max_iter` and `-pc_hypre_boomeramg_tol` refer to the number of iterations
2261: (V-cycles) and tolerance that boomerAMG does EACH time it is called. So for example, if
2262: `-pc_hypre_boomeramg_max_iter` is set to 2 then 2-V-cycles are being used to define the preconditioner
2263: (`-pc_hypre_boomeramg_tol` should be set to 0.0 - the default - to strictly use a fixed number of
2264: iterations per hypre call). `-ksp_max_it` and `-ksp_rtol` STILL determine the total number of iterations
2265: and tolerance for the Krylov solver. For example, if `-pc_hypre_boomeramg_max_iter` is 2 and `-ksp_max_it` is 10
2266: then AT MOST twenty V-cycles of BoomerAMG will be used.
2268: Note that the option `-pc_hypre_boomeramg_relax_type_all` defaults to symmetric relaxation
2269: (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry.
2270: Otherwise, you may want to use `-pc_hypre_boomeramg_relax_type_all SOR/Jacobi`.
2272: If you provide a near null space to your matrix with `MatSetNearNullSpace()` it is ignored by hypre's BoomerAMG UNLESS you also use
2273: the following two options: `-pc_hypre_boomeramg_nodal_coarsen <n> -pc_hypre_boomeramg_vec_interp_variant <v>`
2275: See `PCPFMG`, `PCSMG`, and `PCSYSPFMG` for access to hypre's other (nonalgebraic) multigrid solvers
2277: For `PCHYPRE` type of `ams` or `ads` auxiliary data must be provided to the preconditioner with `PCHYPRESetDiscreteGradient()`,
2278: `PCHYPRESetDiscreteCurl()`, `PCHYPRESetInterpolations()`, `PCHYPRESetAlphaPoissonMatrix()`, `PCHYPRESetBetaPoissonMatrix()`, `PCHYPRESetEdgeConstantVectors()`,
2279: `PCHYPREAMSSetInteriorNodes()`
2281: 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
2282: since it is usually slower, one should use a `KSPType` of `KSPRICHARDSON`
2283: (or equivalently `-ksp_type richardson`) to achieve this. Using `KSPPREONLY` will not work since it only applies a single cycle of multigrid.
2285: PETSc provides its own geometric and algebraic multigrid solvers `PCMG` and `PCGAMG`, also see `PCHMG` which is useful for certain multicomponent problems.
2287: hypre supports performance logging via the `Caliper` library. With `--download-hypre --download-caliper`, hypre will be automatically configured with the support.
2289: Enabling Caliper logging requires setting the `CALI_CONFIG` environment variable before running your hypre code. For example,
2291: .vb
2292: export CALI_CONFIG=runtime-report,max_column_width=200,calc.inclusive,mpi-report,output=stdout
2293: .ve
2295: Then run a hypre code, and you will see profiling results on `stdout`. See https://software.llnl.gov/Caliper/#guides for more options.
2297: GPU Notes:
2298: To configure hypre BoomerAMG so that it can utilize NVIDIA GPUs run `./configure --download-hypre --with-cuda`
2299: Then pass `VECCUDA` vectors and `MATAIJCUSPARSE` matrices to the solvers and PETSc will automatically utilize hypre's GPU solvers.
2301: To configure hypre BoomerAMG so that it can utilize AMD GPUs run ./configure --download-hypre --with-hip
2302: Then pass `VECHIP` vectors to the solvers and PETSc will automatically utilize hypre's GPU solvers.
2304: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHYPRESetType()`, `PCPFMG`, `PCGAMG`, `PCAIR`, `PCSYSPFMG`, `PCSMG`, `PCHYPRESetDiscreteGradient()`,
2305: `PCHYPRESetDiscreteCurl()`, `PCHYPRESetInterpolations()`, `PCHYPRESetAlphaPoissonMatrix()`, `PCHYPRESetBetaPoissonMatrix()`, `PCHYPRESetEdgeConstantVectors()`,
2306: `PCHYPREAMSSetInteriorNodes()`
2307: M*/
2309: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
2310: {
2311: PC_HYPRE *jac;
2313: PetscFunctionBegin;
2314: PetscCall(PetscNew(&jac));
2316: pc->data = jac;
2317: pc->ops->reset = PCReset_HYPRE;
2318: pc->ops->destroy = PCDestroy_HYPRE;
2319: pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
2320: pc->ops->setup = PCSetUp_HYPRE;
2321: pc->ops->apply = PCApply_HYPRE;
2322: jac->hypre_type = NULL;
2323: jac->comm_hypre = MPI_COMM_NULL;
2324: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetType_C", PCHYPRESetType_HYPRE));
2325: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetType_C", PCHYPREGetType_HYPRE));
2326: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_HYPRE));
2327: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteGradient_C", PCHYPRESetDiscreteGradient_HYPRE));
2328: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteCurl_C", PCHYPRESetDiscreteCurl_HYPRE));
2329: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetInterpolations_C", PCHYPRESetInterpolations_HYPRE));
2330: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetEdgeConstantVectors_C", PCHYPRESetEdgeConstantVectors_HYPRE));
2331: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREAMSSetInteriorNodes_C", PCHYPREAMSSetInteriorNodes_HYPRE));
2332: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetPoissonMatrix_C", PCHYPRESetPoissonMatrix_HYPRE));
2333: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinSetMatProductAlgorithm_C", PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG));
2334: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinGetMatProductAlgorithm_C", PCMGGalerkinGetMatProductAlgorithm_HYPRE_BoomerAMG));
2335: #if PetscDefined(HAVE_HYPRE_DEVICE)
2336: #if defined(HYPRE_USING_HIP)
2337: PetscCall(PetscDeviceInitialize(PETSC_DEVICE_HIP));
2338: #endif
2339: #if defined(HYPRE_USING_CUDA)
2340: PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA));
2341: #endif
2342: #endif
2343: PetscCall(PetscHYPREInitialize());
2344: PetscFunctionReturn(PETSC_SUCCESS);
2345: }
2347: typedef struct {
2348: MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */
2349: HYPRE_StructSolver hsolver;
2351: /* keep copy of PFMG options used so may view them */
2352: PetscInt its;
2353: PetscReal tol;
2354: PetscInt relax_type;
2355: PetscInt rap_type;
2356: PetscInt num_pre_relax, num_post_relax;
2357: PetscInt max_levels;
2358: PetscInt skip_relax;
2359: PetscBool print_statistics;
2360: } PC_PFMG;
2362: static PetscErrorCode PCDestroy_PFMG(PC pc)
2363: {
2364: PC_PFMG *ex = (PC_PFMG *)pc->data;
2366: PetscFunctionBegin;
2367: if (ex->hsolver) PetscCallHYPRE(HYPRE_StructPFMGDestroy(ex->hsolver));
2368: PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2369: PetscCall(PetscFree(pc->data));
2370: PetscFunctionReturn(PETSC_SUCCESS);
2371: }
2373: static const char *PFMGRelaxType[] = {"Jacobi", "Weighted-Jacobi", "symmetric-Red/Black-Gauss-Seidel", "Red/Black-Gauss-Seidel"};
2374: static const char *PFMGRAPType[] = {"Galerkin", "non-Galerkin"};
2376: static PetscErrorCode PCView_PFMG(PC pc, PetscViewer viewer)
2377: {
2378: PetscBool isascii;
2379: PC_PFMG *ex = (PC_PFMG *)pc->data;
2381: PetscFunctionBegin;
2382: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2383: if (isascii) {
2384: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE PFMG preconditioning\n"));
2385: PetscCall(PetscViewerASCIIPrintf(viewer, " max iterations %" PetscInt_FMT "\n", ex->its));
2386: PetscCall(PetscViewerASCIIPrintf(viewer, " tolerance %g\n", ex->tol));
2387: PetscCall(PetscViewerASCIIPrintf(viewer, " relax type %s\n", PFMGRelaxType[ex->relax_type]));
2388: PetscCall(PetscViewerASCIIPrintf(viewer, " RAP type %s\n", PFMGRAPType[ex->rap_type]));
2389: PetscCall(PetscViewerASCIIPrintf(viewer, " number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax));
2390: PetscCall(PetscViewerASCIIPrintf(viewer, " max levels %" PetscInt_FMT "\n", ex->max_levels));
2391: PetscCall(PetscViewerASCIIPrintf(viewer, " skip relax %" PetscInt_FMT "\n", ex->skip_relax));
2392: }
2393: PetscFunctionReturn(PETSC_SUCCESS);
2394: }
2396: static PetscErrorCode PCSetFromOptions_PFMG(PC pc, PetscOptionItems PetscOptionsObject)
2397: {
2398: PC_PFMG *ex = (PC_PFMG *)pc->data;
2400: PetscFunctionBegin;
2401: PetscOptionsHeadBegin(PetscOptionsObject, "PFMG options");
2402: PetscCall(PetscOptionsBool("-pc_pfmg_print_statistics", "Print statistics", "HYPRE_StructPFMGSetPrintLevel", ex->print_statistics, &ex->print_statistics, NULL));
2403: PetscCall(PetscOptionsInt("-pc_pfmg_its", "Number of iterations of PFMG to use as preconditioner", "HYPRE_StructPFMGSetMaxIter", ex->its, &ex->its, NULL));
2404: PetscCallHYPRE(HYPRE_StructPFMGSetMaxIter(ex->hsolver, (HYPRE_Int)ex->its));
2405: 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));
2406: PetscCallHYPRE(HYPRE_StructPFMGSetNumPreRelax(ex->hsolver, (HYPRE_Int)ex->num_pre_relax));
2407: 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));
2408: PetscCallHYPRE(HYPRE_StructPFMGSetNumPostRelax(ex->hsolver, (HYPRE_Int)ex->num_post_relax));
2410: PetscCall(PetscOptionsInt("-pc_pfmg_max_levels", "Max Levels for MG hierarchy", "HYPRE_StructPFMGSetMaxLevels", ex->max_levels, &ex->max_levels, NULL));
2411: PetscCallHYPRE(HYPRE_StructPFMGSetMaxLevels(ex->hsolver, (HYPRE_Int)ex->max_levels));
2413: PetscCall(PetscOptionsReal("-pc_pfmg_tol", "Tolerance of PFMG", "HYPRE_StructPFMGSetTol", ex->tol, &ex->tol, NULL));
2414: PetscCallHYPRE(HYPRE_StructPFMGSetTol(ex->hsolver, ex->tol));
2415: 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));
2416: PetscCallHYPRE(HYPRE_StructPFMGSetRelaxType(ex->hsolver, (HYPRE_Int)ex->relax_type));
2417: PetscCall(PetscOptionsEList("-pc_pfmg_rap_type", "RAP type", "HYPRE_StructPFMGSetRAPType", PFMGRAPType, PETSC_STATIC_ARRAY_LENGTH(PFMGRAPType), PFMGRAPType[ex->rap_type], &ex->rap_type, NULL));
2418: PetscCallHYPRE(HYPRE_StructPFMGSetRAPType(ex->hsolver, (HYPRE_Int)ex->rap_type));
2419: 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));
2420: PetscCallHYPRE(HYPRE_StructPFMGSetSkipRelax(ex->hsolver, (HYPRE_Int)ex->skip_relax));
2421: PetscOptionsHeadEnd();
2422: PetscFunctionReturn(PETSC_SUCCESS);
2423: }
2425: static PetscErrorCode PCApply_PFMG(PC pc, Vec x, Vec y)
2426: {
2427: PC_PFMG *ex = (PC_PFMG *)pc->data;
2428: PetscScalar *yy;
2429: const PetscScalar *xx;
2430: PetscInt ilower[3], iupper[3];
2431: HYPRE_Int hlower[3], hupper[3];
2432: Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data;
2434: PetscFunctionBegin;
2435: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
2436: PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
2437: /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2438: iupper[0] += ilower[0] - 1;
2439: iupper[1] += ilower[1] - 1;
2440: iupper[2] += ilower[2] - 1;
2441: hlower[0] = (HYPRE_Int)ilower[0];
2442: hlower[1] = (HYPRE_Int)ilower[1];
2443: hlower[2] = (HYPRE_Int)ilower[2];
2444: hupper[0] = (HYPRE_Int)iupper[0];
2445: hupper[1] = (HYPRE_Int)iupper[1];
2446: hupper[2] = (HYPRE_Int)iupper[2];
2448: /* copy x values over to hypre */
2449: PetscCallHYPRE(HYPRE_StructVectorSetConstantValues(mx->hb, 0.0));
2450: PetscCall(VecGetArrayRead(x, &xx));
2451: PetscCallHYPRE(HYPRE_StructVectorSetBoxValues(mx->hb, hlower, hupper, (HYPRE_Complex *)xx));
2452: PetscCall(VecRestoreArrayRead(x, &xx));
2453: PetscCallHYPRE(HYPRE_StructVectorAssemble(mx->hb));
2454: PetscCallHYPRE(HYPRE_StructPFMGSolve(ex->hsolver, mx->hmat, mx->hb, mx->hx));
2456: /* copy solution values back to PETSc */
2457: PetscCall(VecGetArray(y, &yy));
2458: PetscCallHYPRE(HYPRE_StructVectorGetBoxValues(mx->hx, hlower, hupper, (HYPRE_Complex *)yy));
2459: PetscCall(VecRestoreArray(y, &yy));
2460: PetscFunctionReturn(PETSC_SUCCESS);
2461: }
2463: 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)
2464: {
2465: PC_PFMG *jac = (PC_PFMG *)pc->data;
2466: HYPRE_Int oits;
2468: PetscFunctionBegin;
2469: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
2470: PetscCallHYPRE(HYPRE_StructPFMGSetMaxIter(jac->hsolver, (HYPRE_Int)(its * jac->its)));
2471: PetscCallHYPRE(HYPRE_StructPFMGSetTol(jac->hsolver, rtol));
2473: PetscCall(PCApply_PFMG(pc, b, y));
2474: PetscCallHYPRE(HYPRE_StructPFMGGetNumIterations(jac->hsolver, &oits));
2475: *outits = oits;
2476: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2477: else *reason = PCRICHARDSON_CONVERGED_RTOL;
2478: PetscCallHYPRE(HYPRE_StructPFMGSetTol(jac->hsolver, jac->tol));
2479: PetscCallHYPRE(HYPRE_StructPFMGSetMaxIter(jac->hsolver, (HYPRE_Int)jac->its));
2480: PetscFunctionReturn(PETSC_SUCCESS);
2481: }
2483: static PetscErrorCode PCSetUp_PFMG(PC pc)
2484: {
2485: PC_PFMG *ex = (PC_PFMG *)pc->data;
2486: Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data;
2487: PetscBool flg;
2489: PetscFunctionBegin;
2490: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRESTRUCT, &flg));
2491: PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use MATHYPRESTRUCT with this preconditioner");
2493: /* create the hypre solver object and set its information */
2494: if (ex->hsolver) PetscCallHYPRE(HYPRE_StructPFMGDestroy(ex->hsolver));
2495: PetscCallHYPRE(HYPRE_StructPFMGCreate(ex->hcomm, &ex->hsolver));
2497: // Print Hypre statistics about the solve process
2498: if (ex->print_statistics) PetscCallHYPRE(HYPRE_StructPFMGSetPrintLevel(ex->hsolver, 3));
2500: // The hypre options must be repeated here because the StructPFMG was destroyed and recreated
2501: PetscCallHYPRE(HYPRE_StructPFMGSetMaxIter(ex->hsolver, (HYPRE_Int)ex->its));
2502: PetscCallHYPRE(HYPRE_StructPFMGSetNumPreRelax(ex->hsolver, (HYPRE_Int)ex->num_pre_relax));
2503: PetscCallHYPRE(HYPRE_StructPFMGSetNumPostRelax(ex->hsolver, (HYPRE_Int)ex->num_post_relax));
2504: PetscCallHYPRE(HYPRE_StructPFMGSetMaxLevels(ex->hsolver, (HYPRE_Int)ex->max_levels));
2505: PetscCallHYPRE(HYPRE_StructPFMGSetTol(ex->hsolver, ex->tol));
2506: PetscCallHYPRE(HYPRE_StructPFMGSetRelaxType(ex->hsolver, (HYPRE_Int)ex->relax_type));
2507: PetscCallHYPRE(HYPRE_StructPFMGSetRAPType(ex->hsolver, (HYPRE_Int)ex->rap_type));
2509: PetscCallHYPRE(HYPRE_StructPFMGSetup(ex->hsolver, mx->hmat, mx->hb, mx->hx));
2510: PetscCallHYPRE(HYPRE_StructPFMGSetZeroGuess(ex->hsolver));
2511: PetscFunctionReturn(PETSC_SUCCESS);
2512: }
2514: /*MC
2515: PCPFMG - the hypre PFMG multigrid solver
2517: Options Database Keys:
2518: + -pc_pfmg_its its - number of iterations of PFMG to use as preconditioner
2519: . -pc_pfmg_num_pre_relax steps - number of smoothing steps before coarse grid solve
2520: . -pc_pfmg_num_post_relax steps - number of smoothing steps after coarse grid solve
2521: . -pc_pfmg_tol tol - tolerance of PFMG
2522: . -pc_pfmg_relax_type type - relaxation type for the up and down cycles, one of Jacobi,Weighted-Jacobi,symmetric-Red/Black-Gauss-Seidel,Red/Black-Gauss-Seidel
2523: . -pc_pfmg_rap_type type - type of coarse matrix generation, one of Galerkin,non-Galerkin
2524: - -pc_pfmg_skip_relax (0|1) - skip relaxation on certain grids for isotropic problems. This can greatly improve efficiency by eliminating unnecessary relaxations
2525: when the underlying problem is isotropic, one of 0,1
2527: Level: advanced
2529: Notes:
2530: This is for CELL-centered descretizations
2532: See `PCSYSPFMG` for a version suitable for systems of PDEs, and `PCSMG`
2534: See `PCHYPRE` for hypre's BoomerAMG algebraic multigrid solver
2536: This must be used with the `MATHYPRESTRUCT` matrix type.
2538: This provides only some of the functionality of PFMG, it supports only one block per process defined by a PETSc `DMDA`.
2540: .seealso: [](ch_ksp), `PCMG`, `MATHYPRESTRUCT`, `PCHYPRE`, `PCGAMG`, `PCSYSPFMG`, `PCSMG`
2541: M*/
2543: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
2544: {
2545: PC_PFMG *ex;
2547: PetscFunctionBegin;
2548: PetscCall(PetscNew(&ex));
2549: pc->data = ex;
2551: ex->its = 1;
2552: ex->tol = 1.e-8;
2553: ex->relax_type = 1;
2554: ex->rap_type = 0;
2555: ex->num_pre_relax = 1;
2556: ex->num_post_relax = 1;
2557: ex->max_levels = 0;
2558: ex->skip_relax = 0;
2559: ex->print_statistics = PETSC_FALSE;
2561: pc->ops->setfromoptions = PCSetFromOptions_PFMG;
2562: pc->ops->view = PCView_PFMG;
2563: pc->ops->destroy = PCDestroy_PFMG;
2564: pc->ops->apply = PCApply_PFMG;
2565: pc->ops->applyrichardson = PCApplyRichardson_PFMG;
2566: pc->ops->setup = PCSetUp_PFMG;
2568: PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2569: PetscCall(PetscHYPREInitialize());
2570: PetscCallHYPRE(HYPRE_StructPFMGCreate(ex->hcomm, &ex->hsolver));
2571: PetscFunctionReturn(PETSC_SUCCESS);
2572: }
2574: /* we know we are working with a HYPRE_SStructMatrix */
2575: typedef struct {
2576: MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
2577: HYPRE_SStructSolver ss_solver;
2579: /* keep copy of SYSPFMG options used so may view them */
2580: PetscInt its;
2581: PetscReal tol;
2582: PetscInt relax_type;
2583: PetscInt num_pre_relax, num_post_relax;
2584: } PC_SysPFMG;
2586: static PetscErrorCode PCDestroy_SysPFMG(PC pc)
2587: {
2588: PC_SysPFMG *ex = (PC_SysPFMG *)pc->data;
2590: PetscFunctionBegin;
2591: if (ex->ss_solver) PetscCallHYPRE(HYPRE_SStructSysPFMGDestroy(ex->ss_solver));
2592: PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2593: PetscCall(PetscFree(pc->data));
2594: PetscFunctionReturn(PETSC_SUCCESS);
2595: }
2597: static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi", "Red/Black-Gauss-Seidel"};
2599: static PetscErrorCode PCView_SysPFMG(PC pc, PetscViewer viewer)
2600: {
2601: PetscBool isascii;
2602: PC_SysPFMG *ex = (PC_SysPFMG *)pc->data;
2604: PetscFunctionBegin;
2605: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2606: if (isascii) {
2607: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE SysPFMG preconditioning\n"));
2608: PetscCall(PetscViewerASCIIPrintf(viewer, " max iterations %" PetscInt_FMT "\n", ex->its));
2609: PetscCall(PetscViewerASCIIPrintf(viewer, " tolerance %g\n", ex->tol));
2610: PetscCall(PetscViewerASCIIPrintf(viewer, " relax type %s\n", PFMGRelaxType[ex->relax_type]));
2611: PetscCall(PetscViewerASCIIPrintf(viewer, " number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax));
2612: }
2613: PetscFunctionReturn(PETSC_SUCCESS);
2614: }
2616: static PetscErrorCode PCSetFromOptions_SysPFMG(PC pc, PetscOptionItems PetscOptionsObject)
2617: {
2618: PC_SysPFMG *ex = (PC_SysPFMG *)pc->data;
2619: PetscBool flg = PETSC_FALSE;
2621: PetscFunctionBegin;
2622: PetscOptionsHeadBegin(PetscOptionsObject, "SysPFMG options");
2623: PetscCall(PetscOptionsBool("-pc_syspfmg_print_statistics", "Print statistics", "HYPRE_SStructSysPFMGSetPrintLevel", flg, &flg, NULL));
2624: if (flg) PetscCallHYPRE(HYPRE_SStructSysPFMGSetPrintLevel(ex->ss_solver, 3));
2625: PetscCall(PetscOptionsInt("-pc_syspfmg_its", "Number of iterations of SysPFMG to use as preconditioner", "HYPRE_SStructSysPFMGSetMaxIter", ex->its, &ex->its, NULL));
2626: PetscCallHYPRE(HYPRE_SStructSysPFMGSetMaxIter(ex->ss_solver, (HYPRE_Int)ex->its));
2627: 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));
2628: PetscCallHYPRE(HYPRE_SStructSysPFMGSetNumPreRelax(ex->ss_solver, (HYPRE_Int)ex->num_pre_relax));
2629: 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));
2630: PetscCallHYPRE(HYPRE_SStructSysPFMGSetNumPostRelax(ex->ss_solver, (HYPRE_Int)ex->num_post_relax));
2632: PetscCall(PetscOptionsReal("-pc_syspfmg_tol", "Tolerance of SysPFMG", "HYPRE_SStructSysPFMGSetTol", ex->tol, &ex->tol, NULL));
2633: PetscCallHYPRE(HYPRE_SStructSysPFMGSetTol(ex->ss_solver, ex->tol));
2634: 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));
2635: PetscCallHYPRE(HYPRE_SStructSysPFMGSetRelaxType(ex->ss_solver, (HYPRE_Int)ex->relax_type));
2636: PetscOptionsHeadEnd();
2637: PetscFunctionReturn(PETSC_SUCCESS);
2638: }
2640: static PetscErrorCode PCApply_SysPFMG(PC pc, Vec x, Vec y)
2641: {
2642: PC_SysPFMG *ex = (PC_SysPFMG *)pc->data;
2643: PetscScalar *yy;
2644: const PetscScalar *xx;
2645: PetscInt ilower[3], iupper[3];
2646: HYPRE_Int hlower[3], hupper[3];
2647: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)pc->pmat->data;
2648: PetscInt ordering = mx->dofs_order;
2649: PetscInt nvars = mx->nvars;
2650: HYPRE_Int part = 0;
2651: PetscInt size;
2653: PetscFunctionBegin;
2654: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
2655: PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
2656: /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2657: iupper[0] += ilower[0] - 1;
2658: iupper[1] += ilower[1] - 1;
2659: iupper[2] += ilower[2] - 1;
2660: hlower[0] = (HYPRE_Int)ilower[0];
2661: hlower[1] = (HYPRE_Int)ilower[1];
2662: hlower[2] = (HYPRE_Int)ilower[2];
2663: hupper[0] = (HYPRE_Int)iupper[0];
2664: hupper[1] = (HYPRE_Int)iupper[1];
2665: hupper[2] = (HYPRE_Int)iupper[2];
2667: size = 1;
2668: for (PetscInt i = 0; i < 3; i++) size *= (iupper[i] - ilower[i] + 1);
2670: /* copy x values over to hypre for variable ordering */
2671: if (ordering) {
2672: PetscCallHYPRE(HYPRE_SStructVectorSetConstantValues(mx->ss_b, 0.0));
2673: PetscCall(VecGetArrayRead(x, &xx));
2674: for (PetscInt i = 0; i < nvars; i++) PetscCallHYPRE(HYPRE_SStructVectorSetBoxValues(mx->ss_b, part, hlower, hupper, (HYPRE_Int)i, (HYPRE_Complex *)(xx + (size * i))));
2675: PetscCall(VecRestoreArrayRead(x, &xx));
2676: PetscCallHYPRE(HYPRE_SStructVectorAssemble(mx->ss_b));
2677: PetscCallHYPRE(HYPRE_SStructMatrixMatvec(1.0, mx->ss_mat, mx->ss_b, 0.0, mx->ss_x));
2678: PetscCallHYPRE(HYPRE_SStructSysPFMGSolve(ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x));
2680: /* copy solution values back to PETSc */
2681: PetscCall(VecGetArray(y, &yy));
2682: for (PetscInt i = 0; i < nvars; i++) PetscCallHYPRE(HYPRE_SStructVectorGetBoxValues(mx->ss_x, part, hlower, hupper, (HYPRE_Int)i, (HYPRE_Complex *)(yy + (size * i))));
2683: PetscCall(VecRestoreArray(y, &yy));
2684: } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
2685: PetscScalar *z;
2686: PetscInt k;
2688: PetscCall(PetscMalloc1(nvars * size, &z));
2689: PetscCallHYPRE(HYPRE_SStructVectorSetConstantValues(mx->ss_b, 0.0));
2690: PetscCall(VecGetArrayRead(x, &xx));
2692: /* transform nodal to hypre's variable ordering for sys_pfmg */
2693: for (PetscInt i = 0; i < size; i++) {
2694: k = i * nvars;
2695: for (PetscInt j = 0; j < nvars; j++) z[j * size + i] = xx[k + j];
2696: }
2697: for (PetscInt i = 0; i < nvars; i++) PetscCallHYPRE(HYPRE_SStructVectorSetBoxValues(mx->ss_b, part, hlower, hupper, (HYPRE_Int)i, (HYPRE_Complex *)(z + (size * i))));
2698: PetscCall(VecRestoreArrayRead(x, &xx));
2699: PetscCallHYPRE(HYPRE_SStructVectorAssemble(mx->ss_b));
2700: PetscCallHYPRE(HYPRE_SStructSysPFMGSolve(ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x));
2702: /* copy solution values back to PETSc */
2703: PetscCall(VecGetArray(y, &yy));
2704: for (PetscInt i = 0; i < nvars; i++) PetscCallHYPRE(HYPRE_SStructVectorGetBoxValues(mx->ss_x, part, hlower, hupper, (HYPRE_Int)i, (HYPRE_Complex *)(z + (size * i))));
2705: /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
2706: for (PetscInt i = 0; i < size; i++) {
2707: k = i * nvars;
2708: for (PetscInt j = 0; j < nvars; j++) yy[k + j] = z[j * size + i];
2709: }
2710: PetscCall(VecRestoreArray(y, &yy));
2711: PetscCall(PetscFree(z));
2712: }
2713: PetscFunctionReturn(PETSC_SUCCESS);
2714: }
2716: 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)
2717: {
2718: PC_SysPFMG *jac = (PC_SysPFMG *)pc->data;
2719: HYPRE_Int oits;
2721: PetscFunctionBegin;
2722: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
2723: PetscCallHYPRE(HYPRE_SStructSysPFMGSetMaxIter(jac->ss_solver, (HYPRE_Int)(its * jac->its)));
2724: PetscCallHYPRE(HYPRE_SStructSysPFMGSetTol(jac->ss_solver, rtol));
2725: PetscCall(PCApply_SysPFMG(pc, b, y));
2726: PetscCallHYPRE(HYPRE_SStructSysPFMGGetNumIterations(jac->ss_solver, &oits));
2727: *outits = oits;
2728: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2729: else *reason = PCRICHARDSON_CONVERGED_RTOL;
2730: PetscCallHYPRE(HYPRE_SStructSysPFMGSetTol(jac->ss_solver, jac->tol));
2731: PetscCallHYPRE(HYPRE_SStructSysPFMGSetMaxIter(jac->ss_solver, (HYPRE_Int)jac->its));
2732: PetscFunctionReturn(PETSC_SUCCESS);
2733: }
2735: static PetscErrorCode PCSetUp_SysPFMG(PC pc)
2736: {
2737: PC_SysPFMG *ex = (PC_SysPFMG *)pc->data;
2738: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)pc->pmat->data;
2739: PetscBool flg;
2741: PetscFunctionBegin;
2742: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRESSTRUCT, &flg));
2743: PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use MATHYPRESSTRUCT with this preconditioner");
2745: /* create the hypre sstruct solver object and set its information */
2746: if (ex->ss_solver) PetscCallHYPRE(HYPRE_SStructSysPFMGDestroy(ex->ss_solver));
2747: PetscCallHYPRE(HYPRE_SStructSysPFMGCreate(ex->hcomm, &ex->ss_solver));
2748: PetscCallHYPRE(HYPRE_SStructSysPFMGSetZeroGuess(ex->ss_solver));
2749: PetscCallHYPRE(HYPRE_SStructSysPFMGSetup(ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x));
2750: PetscFunctionReturn(PETSC_SUCCESS);
2751: }
2753: /*MC
2754: PCSYSPFMG - the hypre SysPFMG multigrid solver
2756: Level: advanced
2758: Options Database Keys:
2759: + -pc_syspfmg_its its - number of iterations of SysPFMG to use as preconditioner
2760: . -pc_syspfmg_num_pre_relax steps - number of smoothing steps before coarse grid
2761: . -pc_syspfmg_num_post_relax steps - number of smoothing steps after coarse grid
2762: . -pc_syspfmg_tol tol - tolerance of SysPFMG
2763: - -pc_syspfmg_relax_type (Weighted-Jacobi|Red/Black-Gauss-Seidel) - relaxation type for the up and down cycles
2765: Notes:
2766: See `PCPFMG` for hypre's PFMG that works for a scalar PDE and `PCSMG`
2768: See `PCHYPRE` for hypre's BoomerAMG algebraic multigrid solver
2770: This is for CELL-centered descretizations
2772: This must be used with the `MATHYPRESSTRUCT` matrix type.
2774: 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`.
2776: .seealso: [](ch_ksp), `PCMG`, `MATHYPRESSTRUCT`, `PCPFMG`, `PCHYPRE`, `PCGAMG`, `PCSMG`
2777: M*/
2779: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
2780: {
2781: PC_SysPFMG *ex;
2783: PetscFunctionBegin;
2784: PetscCall(PetscNew(&ex));
2785: pc->data = ex;
2787: ex->its = 1;
2788: ex->tol = 1.e-8;
2789: ex->relax_type = 1;
2790: ex->num_pre_relax = 1;
2791: ex->num_post_relax = 1;
2793: pc->ops->setfromoptions = PCSetFromOptions_SysPFMG;
2794: pc->ops->view = PCView_SysPFMG;
2795: pc->ops->destroy = PCDestroy_SysPFMG;
2796: pc->ops->apply = PCApply_SysPFMG;
2797: pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
2798: pc->ops->setup = PCSetUp_SysPFMG;
2800: PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2801: PetscCall(PetscHYPREInitialize());
2802: PetscCallHYPRE(HYPRE_SStructSysPFMGCreate(ex->hcomm, &ex->ss_solver));
2803: PetscFunctionReturn(PETSC_SUCCESS);
2804: }
2806: /* PC SMG */
2807: typedef struct {
2808: MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */
2809: HYPRE_StructSolver hsolver;
2810: PetscInt its; /* keep copy of SMG options used so may view them */
2811: PetscReal tol;
2812: PetscBool print_statistics;
2813: PetscInt num_pre_relax, num_post_relax;
2814: } PC_SMG;
2816: static PetscErrorCode PCDestroy_SMG(PC pc)
2817: {
2818: PC_SMG *ex = (PC_SMG *)pc->data;
2820: PetscFunctionBegin;
2821: if (ex->hsolver) PetscCallHYPRE(HYPRE_StructSMGDestroy(ex->hsolver));
2822: PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2823: PetscCall(PetscFree(pc->data));
2824: PetscFunctionReturn(PETSC_SUCCESS);
2825: }
2827: static PetscErrorCode PCView_SMG(PC pc, PetscViewer viewer)
2828: {
2829: PetscBool isascii;
2830: PC_SMG *ex = (PC_SMG *)pc->data;
2832: PetscFunctionBegin;
2833: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2834: if (isascii) {
2835: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE SMG preconditioning\n"));
2836: PetscCall(PetscViewerASCIIPrintf(viewer, " max iterations %" PetscInt_FMT "\n", ex->its));
2837: PetscCall(PetscViewerASCIIPrintf(viewer, " tolerance %g\n", ex->tol));
2838: PetscCall(PetscViewerASCIIPrintf(viewer, " number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax));
2839: }
2840: PetscFunctionReturn(PETSC_SUCCESS);
2841: }
2843: static PetscErrorCode PCSetFromOptions_SMG(PC pc, PetscOptionItems PetscOptionsObject)
2844: {
2845: PC_SMG *ex = (PC_SMG *)pc->data;
2847: PetscFunctionBegin;
2848: PetscOptionsHeadBegin(PetscOptionsObject, "SMG options");
2850: PetscCall(PetscOptionsInt("-pc_smg_its", "Number of iterations of SMG to use as preconditioner", "HYPRE_StructSMGSetMaxIter", ex->its, &ex->its, NULL));
2851: 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));
2852: 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));
2853: PetscCall(PetscOptionsReal("-pc_smg_tol", "Tolerance of SMG", "HYPRE_StructSMGSetTol", ex->tol, &ex->tol, NULL));
2855: PetscOptionsHeadEnd();
2856: PetscFunctionReturn(PETSC_SUCCESS);
2857: }
2859: static PetscErrorCode PCApply_SMG(PC pc, Vec x, Vec y)
2860: {
2861: PC_SMG *ex = (PC_SMG *)pc->data;
2862: PetscScalar *yy;
2863: const PetscScalar *xx;
2864: PetscInt ilower[3], iupper[3];
2865: HYPRE_Int hlower[3], hupper[3];
2866: Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data;
2868: PetscFunctionBegin;
2869: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
2870: PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
2871: /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2872: iupper[0] += ilower[0] - 1;
2873: iupper[1] += ilower[1] - 1;
2874: iupper[2] += ilower[2] - 1;
2875: hlower[0] = (HYPRE_Int)ilower[0];
2876: hlower[1] = (HYPRE_Int)ilower[1];
2877: hlower[2] = (HYPRE_Int)ilower[2];
2878: hupper[0] = (HYPRE_Int)iupper[0];
2879: hupper[1] = (HYPRE_Int)iupper[1];
2880: hupper[2] = (HYPRE_Int)iupper[2];
2882: /* copy x values over to hypre */
2883: PetscCallHYPRE(HYPRE_StructVectorSetConstantValues(mx->hb, 0.0));
2884: PetscCall(VecGetArrayRead(x, &xx));
2885: PetscCallHYPRE(HYPRE_StructVectorSetBoxValues(mx->hb, hlower, hupper, (HYPRE_Complex *)xx));
2886: PetscCall(VecRestoreArrayRead(x, &xx));
2887: PetscCallHYPRE(HYPRE_StructVectorAssemble(mx->hb));
2888: PetscCallHYPRE(HYPRE_StructSMGSolve(ex->hsolver, mx->hmat, mx->hb, mx->hx));
2890: /* copy solution values back to PETSc */
2891: PetscCall(VecGetArray(y, &yy));
2892: PetscCallHYPRE(HYPRE_StructVectorGetBoxValues(mx->hx, hlower, hupper, (HYPRE_Complex *)yy));
2893: PetscCall(VecRestoreArray(y, &yy));
2894: PetscFunctionReturn(PETSC_SUCCESS);
2895: }
2897: 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)
2898: {
2899: PC_SMG *jac = (PC_SMG *)pc->data;
2900: HYPRE_Int oits;
2902: PetscFunctionBegin;
2903: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
2904: PetscCallHYPRE(HYPRE_StructSMGSetMaxIter(jac->hsolver, (HYPRE_Int)(its * jac->its)));
2905: PetscCallHYPRE(HYPRE_StructSMGSetTol(jac->hsolver, rtol));
2907: PetscCall(PCApply_SMG(pc, b, y));
2908: PetscCallHYPRE(HYPRE_StructSMGGetNumIterations(jac->hsolver, &oits));
2909: *outits = oits;
2910: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2911: else *reason = PCRICHARDSON_CONVERGED_RTOL;
2912: PetscCallHYPRE(HYPRE_StructSMGSetTol(jac->hsolver, jac->tol));
2913: PetscCallHYPRE(HYPRE_StructSMGSetMaxIter(jac->hsolver, (HYPRE_Int)jac->its));
2914: PetscFunctionReturn(PETSC_SUCCESS);
2915: }
2917: static PetscErrorCode PCSetUp_SMG(PC pc)
2918: {
2919: PetscInt dim;
2920: PC_SMG *ex = (PC_SMG *)pc->data;
2921: Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data;
2922: PetscBool flg;
2923: DMBoundaryType p[3];
2924: PetscInt M[3];
2926: PetscFunctionBegin;
2927: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRESTRUCT, &flg));
2928: PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use MATHYPRESTRUCT with this preconditioner");
2930: PetscCall(DMDAGetInfo(mx->da, &dim, &M[0], &M[1], &M[2], 0, 0, 0, 0, 0, &p[0], &p[1], &p[2], 0));
2931: // Check if power of 2 in periodic directions
2932: for (PetscInt i = 0; i < dim; i++) {
2933: PetscCheck((M[i] & (M[i] - 1)) == 0 || p[i] != DM_BOUNDARY_PERIODIC, 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]);
2934: }
2936: /* create the hypre solver object and set its information */
2937: if (ex->hsolver) PetscCallHYPRE(HYPRE_StructSMGDestroy(ex->hsolver));
2938: PetscCallHYPRE(HYPRE_StructSMGCreate(ex->hcomm, &ex->hsolver));
2939: // The hypre options must be set here and not in SetFromOptions because it is created here!
2940: PetscCallHYPRE(HYPRE_StructSMGSetMaxIter(ex->hsolver, (HYPRE_Int)ex->its));
2941: PetscCallHYPRE(HYPRE_StructSMGSetNumPreRelax(ex->hsolver, (HYPRE_Int)ex->num_pre_relax));
2942: PetscCallHYPRE(HYPRE_StructSMGSetNumPostRelax(ex->hsolver, (HYPRE_Int)ex->num_post_relax));
2943: PetscCallHYPRE(HYPRE_StructSMGSetTol(ex->hsolver, ex->tol));
2945: PetscCallHYPRE(HYPRE_StructSMGSetup(ex->hsolver, mx->hmat, mx->hb, mx->hx));
2946: PetscCallHYPRE(HYPRE_StructSMGSetZeroGuess(ex->hsolver));
2947: PetscFunctionReturn(PETSC_SUCCESS);
2948: }
2950: /*MC
2951: PCSMG - the hypre (structured grid) SMG multigrid solver
2953: Level: advanced
2955: Options Database Keys:
2956: + -pc_smg_its its - number of iterations of SMG to use as preconditioner
2957: . -pc_smg_num_pre_relax steps - number of smoothing steps before coarse grid
2958: . -pc_smg_num_post_relax steps - number of smoothing steps after coarse grid
2959: - -pc_smg_tol tol - tolerance of SMG
2961: Notes:
2962: This is for CELL-centered descretizations
2964: This must be used with the `MATHYPRESTRUCT` `MatType`.
2966: This does not provide all the functionality of hypre's SMG solver, it supports only one block per process defined by a PETSc `DMDA`.
2968: See `PCSYSPFMG`, `PCSMG`, `PCPFMG`, and `PCHYPRE` for access to hypre's other preconditioners
2970: .seealso: `PCMG`, `MATHYPRESTRUCT`, `PCPFMG`, `PCSYSPFMG`, `PCHYPRE`, `PCGAMG`
2971: M*/
2973: PETSC_EXTERN PetscErrorCode PCCreate_SMG(PC pc)
2974: {
2975: PC_SMG *ex;
2977: PetscFunctionBegin;
2978: PetscCall(PetscNew(&ex));
2979: pc->data = ex;
2981: ex->its = 1;
2982: ex->tol = 1.e-8;
2983: ex->num_pre_relax = 1;
2984: ex->num_post_relax = 1;
2986: pc->ops->setfromoptions = PCSetFromOptions_SMG;
2987: pc->ops->view = PCView_SMG;
2988: pc->ops->destroy = PCDestroy_SMG;
2989: pc->ops->apply = PCApply_SMG;
2990: pc->ops->applyrichardson = PCApplyRichardson_SMG;
2991: pc->ops->setup = PCSetUp_SMG;
2993: PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2994: PetscCall(PetscHYPREInitialize());
2995: PetscCallHYPRE(HYPRE_StructSMGCreate(ex->hcomm, &ex->hsolver));
2996: PetscFunctionReturn(PETSC_SUCCESS);
2997: }