Actual source code: hypre.c
1: /*
2: Provides an interface to the LLNL package hypre
3: */
5: #include <petscpkg_version.h>
6: #include <petsc/private/pcimpl.h>
7: /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */
8: #include <petsc/private/matimpl.h>
9: #include <petsc/private/vecimpl.h>
10: #include <../src/vec/vec/impls/hypre/vhyp.h>
11: #include <../src/mat/impls/hypre/mhypre.h>
12: #include <../src/dm/impls/da/hypre/mhyp.h>
13: #include <_hypre_parcsr_ls.h>
14: #include <petscmathypre.h>
16: #if defined(PETSC_HAVE_HYPRE_DEVICE)
17: #include <petsc/private/deviceimpl.h>
18: #endif
20: static PetscBool cite = PETSC_FALSE;
21: static const char hypreCitation[] = "@manual{hypre-web-page,\n title = {{\\sl hypre}: High Performance Preconditioners},\n organization = {Lawrence Livermore National Laboratory},\n note = "
22: "{\\url{https://www.llnl.gov/casc/hypre}}\n}\n";
24: /*
25: Private context (data structure) for the preconditioner.
26: */
27: typedef struct {
28: HYPRE_Solver hsolver;
29: Mat hpmat; /* MatHYPRE */
31: HYPRE_Int (*destroy)(HYPRE_Solver);
32: HYPRE_Int (*solve)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector);
33: HYPRE_Int (*setup)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector);
35: MPI_Comm comm_hypre;
36: char *hypre_type;
38: /* options for Pilut and BoomerAMG*/
39: PetscInt maxiter;
40: PetscReal tol;
42: /* options for Pilut */
43: PetscInt factorrowsize;
45: /* options for ParaSails */
46: PetscInt nlevels;
47: PetscReal threshold;
48: PetscReal filter;
49: PetscReal loadbal;
50: PetscInt logging;
51: PetscInt ruse;
52: PetscInt symt;
54: /* options for BoomerAMG */
55: PetscBool printstatistics;
57: /* options for BoomerAMG */
58: PetscInt cycletype;
59: PetscInt maxlevels;
60: PetscReal strongthreshold;
61: PetscReal maxrowsum;
62: PetscInt gridsweeps[3];
63: 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(PetscBool, 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, l;
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 (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, l;
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 (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;
262: PetscBool flag;
264: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
265: PetscFunctionBegin;
266: if (jac->spgemm_type) {
267: PetscCall(PetscStrcmp(jac->spgemm_type, name, &flag));
268: PetscCheck(flag, PetscObjectComm((PetscObject)pc), PETSC_ERR_ORDER, "PETSc support for resetting the HYPRE SpGEMM is not implemented");
269: PetscFunctionReturn(PETSC_SUCCESS);
270: } else jac->spgemm_type = name;
272: PetscCall(PetscStrcmp("cusparse", jac->spgemm_type, &flag));
273: if (flag) {
274: PetscCallExternal(HYPRE_SetSpGemmUseCusparse, 1);
275: PetscFunctionReturn(PETSC_SUCCESS);
276: }
277: PetscCall(PetscStrcmp("hypre", jac->spgemm_type, &flag));
278: if (flag) {
279: PetscCallExternal(HYPRE_SetSpGemmUseCusparse, 0);
280: PetscFunctionReturn(PETSC_SUCCESS);
281: }
282: jac->spgemm_type = NULL;
283: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown HYPRE SpGEMM type %s; Choices are cusparse, hypre", name);
284: #endif
285: }
287: static PetscErrorCode PCSetUp_HYPRE(PC pc)
288: {
289: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
290: Mat_HYPRE *hjac;
291: HYPRE_ParCSRMatrix hmat;
292: HYPRE_ParVector bv, xv;
293: PetscBool ishypre;
295: PetscFunctionBegin;
296: /* default type is boomerAMG */
297: if (!jac->hypre_type) PetscCall(PCHYPRESetType(pc, "boomeramg"));
299: /* get hypre matrix */
300: if (pc->flag == DIFFERENT_NONZERO_PATTERN) PetscCall(MatDestroy(&jac->hpmat));
301: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRE, &ishypre));
302: if (!ishypre) {
303: #if defined(PETSC_HAVE_HYPRE_DEVICE) && PETSC_PKG_HYPRE_VERSION_LE(2, 30, 0)
304: /* Temporary fix since we do not support MAT_REUSE_MATRIX with HYPRE device */
305: PetscBool iscuda, iship, iskokkos;
307: PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &iscuda, MATSEQAIJCUSPARSE, MATMPIAIJCUSPARSE, ""));
308: PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &iship, MATSEQAIJHIPSPARSE, MATMPIAIJHIPSPARSE, ""));
309: PetscCall(PetscObjectTypeCompareAny((PetscObject)pc->pmat, &iskokkos, MATSEQAIJKOKKOS, MATMPIAIJKOKKOS, ""));
310: if (iscuda || iship || iskokkos) PetscCall(MatDestroy(&jac->hpmat));
311: #endif
312: PetscCall(MatConvert(pc->pmat, MATHYPRE, jac->hpmat ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX, &jac->hpmat));
313: } else {
314: PetscCall(PetscObjectReference((PetscObject)pc->pmat));
315: PetscCall(MatDestroy(&jac->hpmat));
316: jac->hpmat = pc->pmat;
317: }
319: /* allow debug */
320: PetscCall(MatViewFromOptions(jac->hpmat, NULL, "-pc_hypre_mat_view"));
321: hjac = (Mat_HYPRE *)jac->hpmat->data;
323: /* special case for BoomerAMG */
324: if (jac->setup == HYPRE_BoomerAMGSetup) {
325: MatNullSpace mnull;
326: PetscBool has_const;
327: PetscInt bs, nvec, i;
328: PetscMemType memtype;
329: const Vec *vecs;
331: PetscCall(MatGetCurrentMemType(jac->hpmat, &memtype));
332: if (PetscMemTypeDevice(memtype)) {
333: /* GPU defaults
334: From https://hypre.readthedocs.io/en/latest/solvers-boomeramg.html#gpu-supported-options
335: and /src/parcsr_ls/par_amg.c
336: First handle options which users have interfaces for changing */
337: PetscObjectParameterSetDefault(jac, coarsentype, 8);
338: PetscObjectParameterSetDefault(jac, relaxorder, 0);
339: PetscObjectParameterSetDefault(jac, interptype, 6);
340: PetscObjectParameterSetDefault(jac, relaxtype[0], 18);
341: PetscObjectParameterSetDefault(jac, relaxtype[1], 18);
342: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
343: PetscObjectParameterSetDefault(jac, spgemm_type, HYPRESpgemmTypes[0]);
344: #endif
345: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
346: PetscObjectParameterSetDefault(jac, keeptranspose, PETSC_TRUE);
347: PetscObjectParameterSetDefault(jac, mod_rap2, 1);
348: #endif
349: PetscObjectParameterSetDefault(jac, agg_interptype, 7);
350: } else {
351: PetscObjectParameterSetDefault(jac, coarsentype, 6);
352: PetscObjectParameterSetDefault(jac, relaxorder, 1);
353: PetscObjectParameterSetDefault(jac, interptype, 0);
354: PetscObjectParameterSetDefault(jac, relaxtype[0], 6);
355: PetscObjectParameterSetDefault(jac, relaxtype[1], 6); /* Defaults to SYMMETRIC since in PETSc we are using a PC - most likely with CG */
356: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
357: PetscObjectParameterSetDefault(jac, spgemm_type, "hypre");
358: #endif
359: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
360: PetscObjectParameterSetDefault(jac, keeptranspose, PETSC_FALSE);
361: PetscObjectParameterSetDefault(jac, mod_rap2, 0);
362: #endif
363: PetscObjectParameterSetDefault(jac, agg_interptype, 4);
364: }
365: PetscObjectParameterSetDefault(jac, relaxtype[2], 9); /*G.E. */
367: PetscCallExternal(HYPRE_BoomerAMGSetCycleType, jac->hsolver, jac->cycletype);
368: PetscCallExternal(HYPRE_BoomerAMGSetMaxLevels, jac->hsolver, jac->maxlevels);
369: PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, jac->maxiter);
370: PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, jac->tol);
371: PetscCallExternal(HYPRE_BoomerAMGSetTruncFactor, jac->hsolver, jac->truncfactor);
372: PetscCallExternal(HYPRE_BoomerAMGSetStrongThreshold, jac->hsolver, jac->strongthreshold);
373: PetscCallExternal(HYPRE_BoomerAMGSetMaxRowSum, jac->hsolver, jac->maxrowsum);
374: PetscCallExternal(HYPRE_BoomerAMGSetMeasureType, jac->hsolver, jac->measuretype);
375: PetscCallExternal(HYPRE_BoomerAMGSetAggNumLevels, jac->hsolver, jac->agg_nl);
376: PetscCallExternal(HYPRE_BoomerAMGSetPMaxElmts, jac->hsolver, jac->pmax);
377: PetscCallExternal(HYPRE_BoomerAMGSetNumPaths, jac->hsolver, jac->agg_num_paths);
378: PetscCallExternal(HYPRE_BoomerAMGSetCycleNumSweeps, jac->hsolver, jac->gridsweeps[0], 1);
379: PetscCallExternal(HYPRE_BoomerAMGSetCycleNumSweeps, jac->hsolver, jac->gridsweeps[1], 2);
380: PetscCallExternal(HYPRE_BoomerAMGSetCycleNumSweeps, jac->hsolver, jac->gridsweeps[2], 3);
381: PetscCallExternal(HYPRE_BoomerAMGSetMaxCoarseSize, jac->hsolver, jac->maxc);
382: PetscCallExternal(HYPRE_BoomerAMGSetMinCoarseSize, jac->hsolver, jac->minc);
383: PetscCallExternal(HYPRE_BoomerAMGSetCoarsenType, jac->hsolver, jac->coarsentype);
384: PetscCallExternal(HYPRE_BoomerAMGSetRelaxOrder, jac->hsolver, jac->relaxorder);
385: PetscCallExternal(HYPRE_BoomerAMGSetInterpType, jac->hsolver, jac->interptype);
386: PetscCallExternal(HYPRE_BoomerAMGSetRelaxType, jac->hsolver, jac->relaxtype[0]);
387: /* GPU */
388: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
389: PetscCall(PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG(pc, jac->spgemm_type));
390: #endif
391: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
392: PetscCallExternal(HYPRE_BoomerAMGSetKeepTranspose, jac->hsolver, jac->keeptranspose ? 1 : 0);
393: PetscCallExternal(HYPRE_BoomerAMGSetRAP2, jac->hsolver, jac->rap2);
394: PetscCallExternal(HYPRE_BoomerAMGSetModuleRAP2, jac->hsolver, jac->mod_rap2);
395: #endif
396: PetscCallExternal(HYPRE_BoomerAMGSetAggInterpType, jac->hsolver, jac->agg_interptype);
398: /* AIR */
399: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
400: PetscCallExternal(HYPRE_BoomerAMGSetRestriction, jac->hsolver, jac->Rtype);
401: PetscCallExternal(HYPRE_BoomerAMGSetStrongThresholdR, jac->hsolver, jac->Rstrongthreshold);
402: PetscCallExternal(HYPRE_BoomerAMGSetFilterThresholdR, jac->hsolver, jac->Rfilterthreshold);
403: PetscCallExternal(HYPRE_BoomerAMGSetADropTol, jac->hsolver, jac->Adroptol);
404: PetscCallExternal(HYPRE_BoomerAMGSetADropType, jac->hsolver, jac->Adroptype);
405: #endif
407: PetscCall(MatGetBlockSize(pc->pmat, &bs));
408: if (bs > 1) PetscCallExternal(HYPRE_BoomerAMGSetNumFunctions, jac->hsolver, bs);
409: PetscCall(MatGetNearNullSpace(pc->mat, &mnull));
410: if (mnull) {
411: PetscCall(PCHYPREResetNearNullSpace_Private(pc));
412: PetscCall(MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs));
413: PetscCall(PetscMalloc1(nvec + 1, &jac->hmnull));
414: PetscCall(PetscMalloc1(nvec + 1, &jac->phmnull));
415: for (i = 0; i < nvec; i++) {
416: PetscCall(VecHYPRE_IJVectorCreate(vecs[i]->map, &jac->hmnull[i]));
417: PetscCall(VecHYPRE_IJVectorCopy(vecs[i], jac->hmnull[i]));
418: PetscCallExternal(HYPRE_IJVectorGetObject, jac->hmnull[i]->ij, (void **)&jac->phmnull[i]);
419: }
420: if (has_const) {
421: PetscCall(MatCreateVecs(pc->pmat, &jac->hmnull_constant, NULL));
422: PetscCall(VecSet(jac->hmnull_constant, 1));
423: PetscCall(VecNormalize(jac->hmnull_constant, NULL));
424: PetscCall(VecHYPRE_IJVectorCreate(jac->hmnull_constant->map, &jac->hmnull[nvec]));
425: PetscCall(VecHYPRE_IJVectorCopy(jac->hmnull_constant, jac->hmnull[nvec]));
426: PetscCallExternal(HYPRE_IJVectorGetObject, jac->hmnull[nvec]->ij, (void **)&jac->phmnull[nvec]);
427: nvec++;
428: }
429: PetscCallExternal(HYPRE_BoomerAMGSetInterpVectors, jac->hsolver, nvec, jac->phmnull);
430: jac->n_hmnull = nvec;
431: }
432: }
434: /* special case for AMS */
435: if (jac->setup == HYPRE_AMSSetup) {
436: Mat_HYPRE *hm;
437: HYPRE_ParCSRMatrix parcsr;
438: 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()");
439: if (jac->dim) PetscCallExternal(HYPRE_AMSSetDimension, jac->hsolver, jac->dim);
440: if (jac->constants[0]) {
441: HYPRE_ParVector ozz, zoz, zzo = NULL;
442: PetscCallExternal(HYPRE_IJVectorGetObject, jac->constants[0]->ij, (void **)(&ozz));
443: PetscCallExternal(HYPRE_IJVectorGetObject, jac->constants[1]->ij, (void **)(&zoz));
444: if (jac->constants[2]) PetscCallExternal(HYPRE_IJVectorGetObject, jac->constants[2]->ij, (void **)(&zzo));
445: PetscCallExternal(HYPRE_AMSSetEdgeConstantVectors, jac->hsolver, ozz, zoz, zzo);
446: }
447: if (jac->coords[0]) {
448: HYPRE_ParVector coords[3];
449: coords[0] = NULL;
450: coords[1] = NULL;
451: coords[2] = NULL;
452: if (jac->coords[0]) PetscCallExternal(HYPRE_IJVectorGetObject, jac->coords[0]->ij, (void **)(&coords[0]));
453: if (jac->coords[1]) PetscCallExternal(HYPRE_IJVectorGetObject, jac->coords[1]->ij, (void **)(&coords[1]));
454: if (jac->coords[2]) PetscCallExternal(HYPRE_IJVectorGetObject, jac->coords[2]->ij, (void **)(&coords[2]));
455: PetscCallExternal(HYPRE_AMSSetCoordinateVectors, jac->hsolver, coords[0], coords[1], coords[2]);
456: }
457: PetscCheck(jac->G, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE AMS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient");
458: hm = (Mat_HYPRE *)jac->G->data;
459: PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&parcsr));
460: PetscCallExternal(HYPRE_AMSSetDiscreteGradient, jac->hsolver, parcsr);
461: if (jac->alpha_Poisson) {
462: hm = (Mat_HYPRE *)jac->alpha_Poisson->data;
463: PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&parcsr));
464: PetscCallExternal(HYPRE_AMSSetAlphaPoissonMatrix, jac->hsolver, parcsr);
465: }
466: if (jac->ams_beta_is_zero) {
467: PetscCallExternal(HYPRE_AMSSetBetaPoissonMatrix, jac->hsolver, NULL);
468: } else if (jac->beta_Poisson) {
469: hm = (Mat_HYPRE *)jac->beta_Poisson->data;
470: PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&parcsr));
471: PetscCallExternal(HYPRE_AMSSetBetaPoissonMatrix, jac->hsolver, parcsr);
472: } else if (jac->ams_beta_is_zero_part) {
473: if (jac->interior) {
474: HYPRE_ParVector interior = NULL;
475: PetscCallExternal(HYPRE_IJVectorGetObject, jac->interior->ij, (void **)(&interior));
476: PetscCallExternal(HYPRE_AMSSetInteriorNodes, jac->hsolver, interior);
477: } else {
478: jac->ams_beta_is_zero_part = PETSC_FALSE;
479: }
480: }
481: if (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1])) {
482: PetscInt i;
483: HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3];
484: if (jac->ND_PiFull) {
485: hm = (Mat_HYPRE *)jac->ND_PiFull->data;
486: PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&nd_parcsrfull));
487: } else {
488: nd_parcsrfull = NULL;
489: }
490: for (i = 0; i < 3; ++i) {
491: if (jac->ND_Pi[i]) {
492: hm = (Mat_HYPRE *)jac->ND_Pi[i]->data;
493: PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&nd_parcsr[i]));
494: } else {
495: nd_parcsr[i] = NULL;
496: }
497: }
498: PetscCallExternal(HYPRE_AMSSetInterpolations, jac->hsolver, nd_parcsrfull, nd_parcsr[0], nd_parcsr[1], nd_parcsr[2]);
499: }
500: }
501: /* special case for ADS */
502: if (jac->setup == HYPRE_ADSSetup) {
503: Mat_HYPRE *hm;
504: HYPRE_ParCSRMatrix parcsr;
505: 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])))) {
506: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE ADS preconditioner needs either the coordinate vectors via PCSetCoordinates() or the interpolation matrices via PCHYPRESetInterpolations");
507: } 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");
508: PetscCheck(jac->G, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE ADS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient");
509: PetscCheck(jac->C, PetscObjectComm((PetscObject)pc), PETSC_ERR_USER, "HYPRE ADS preconditioner needs the discrete curl operator via PCHYPRESetDiscreteGradient");
510: if (jac->coords[0]) {
511: HYPRE_ParVector coords[3];
512: coords[0] = NULL;
513: coords[1] = NULL;
514: coords[2] = NULL;
515: if (jac->coords[0]) PetscCallExternal(HYPRE_IJVectorGetObject, jac->coords[0]->ij, (void **)(&coords[0]));
516: if (jac->coords[1]) PetscCallExternal(HYPRE_IJVectorGetObject, jac->coords[1]->ij, (void **)(&coords[1]));
517: if (jac->coords[2]) PetscCallExternal(HYPRE_IJVectorGetObject, jac->coords[2]->ij, (void **)(&coords[2]));
518: PetscCallExternal(HYPRE_ADSSetCoordinateVectors, jac->hsolver, coords[0], coords[1], coords[2]);
519: }
520: hm = (Mat_HYPRE *)jac->G->data;
521: PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&parcsr));
522: PetscCallExternal(HYPRE_ADSSetDiscreteGradient, jac->hsolver, parcsr);
523: hm = (Mat_HYPRE *)jac->C->data;
524: PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&parcsr));
525: PetscCallExternal(HYPRE_ADSSetDiscreteCurl, jac->hsolver, parcsr);
526: if ((jac->RT_PiFull || (jac->RT_Pi[0] && jac->RT_Pi[1])) && (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]))) {
527: PetscInt i;
528: HYPRE_ParCSRMatrix rt_parcsrfull, rt_parcsr[3];
529: HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3];
530: if (jac->RT_PiFull) {
531: hm = (Mat_HYPRE *)jac->RT_PiFull->data;
532: PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&rt_parcsrfull));
533: } else {
534: rt_parcsrfull = NULL;
535: }
536: for (i = 0; i < 3; ++i) {
537: if (jac->RT_Pi[i]) {
538: hm = (Mat_HYPRE *)jac->RT_Pi[i]->data;
539: PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&rt_parcsr[i]));
540: } else {
541: rt_parcsr[i] = NULL;
542: }
543: }
544: if (jac->ND_PiFull) {
545: hm = (Mat_HYPRE *)jac->ND_PiFull->data;
546: PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&nd_parcsrfull));
547: } else {
548: nd_parcsrfull = NULL;
549: }
550: for (i = 0; i < 3; ++i) {
551: if (jac->ND_Pi[i]) {
552: hm = (Mat_HYPRE *)jac->ND_Pi[i]->data;
553: PetscCallExternal(HYPRE_IJMatrixGetObject, hm->ij, (void **)(&nd_parcsr[i]));
554: } else {
555: nd_parcsr[i] = NULL;
556: }
557: }
558: PetscCallExternal(HYPRE_ADSSetInterpolations, jac->hsolver, rt_parcsrfull, rt_parcsr[0], rt_parcsr[1], rt_parcsr[2], nd_parcsrfull, nd_parcsr[0], nd_parcsr[1], nd_parcsr[2]);
559: }
560: }
561: PetscCallExternal(HYPRE_IJMatrixGetObject, hjac->ij, (void **)&hmat);
562: PetscCallExternal(HYPRE_IJVectorGetObject, hjac->b->ij, (void **)&bv);
563: PetscCallExternal(HYPRE_IJVectorGetObject, hjac->x->ij, (void **)&xv);
564: PetscCall(PetscFPTrapPush(PETSC_FP_TRAP_OFF));
565: PetscCallExternal(jac->setup, jac->hsolver, hmat, bv, xv);
566: PetscCall(PetscFPTrapPop());
567: PetscFunctionReturn(PETSC_SUCCESS);
568: }
570: static PetscErrorCode PCApply_HYPRE(PC pc, Vec b, Vec x)
571: {
572: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
573: Mat_HYPRE *hjac = (Mat_HYPRE *)jac->hpmat->data;
574: HYPRE_ParCSRMatrix hmat;
575: HYPRE_ParVector jbv, jxv;
577: PetscFunctionBegin;
578: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
579: if (!jac->applyrichardson) PetscCall(VecSet(x, 0.0));
580: PetscCall(VecHYPRE_IJVectorPushVecRead(hjac->b, b));
581: if (jac->applyrichardson) PetscCall(VecHYPRE_IJVectorPushVec(hjac->x, x));
582: else PetscCall(VecHYPRE_IJVectorPushVecWrite(hjac->x, x));
583: PetscCallExternal(HYPRE_IJMatrixGetObject, hjac->ij, (void **)&hmat);
584: PetscCallExternal(HYPRE_IJVectorGetObject, hjac->b->ij, (void **)&jbv);
585: PetscCallExternal(HYPRE_IJVectorGetObject, hjac->x->ij, (void **)&jxv);
586: PetscStackCallExternalVoid(
587: "Hypre solve", do {
588: HYPRE_Int hierr = (*jac->solve)(jac->hsolver, hmat, jbv, jxv);
589: if (hierr) {
590: PetscCheck(hierr == HYPRE_ERROR_CONV, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in HYPRE solver, error code %d", (int)hierr);
591: HYPRE_ClearAllErrors();
592: }
593: } while (0));
595: if (jac->setup == HYPRE_AMSSetup && jac->ams_beta_is_zero_part) PetscCallExternal(HYPRE_AMSProjectOutGradients, jac->hsolver, jxv);
596: PetscCall(VecHYPRE_IJVectorPopVec(hjac->x));
597: PetscCall(VecHYPRE_IJVectorPopVec(hjac->b));
598: PetscFunctionReturn(PETSC_SUCCESS);
599: }
601: static PetscErrorCode PCMatApply_HYPRE_BoomerAMG(PC pc, Mat B, Mat X)
602: {
603: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
604: Mat_HYPRE *hjac = (Mat_HYPRE *)jac->hpmat->data;
605: hypre_ParCSRMatrix *par_matrix;
606: HYPRE_ParVector hb, hx;
607: const PetscScalar *b;
608: PetscScalar *x;
609: PetscInt m, N, lda;
610: hypre_Vector *x_local;
611: PetscMemType type;
613: PetscFunctionBegin;
614: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
615: PetscCallExternal(HYPRE_IJMatrixGetObject, hjac->ij, (void **)&par_matrix);
616: PetscCall(MatGetLocalSize(B, &m, NULL));
617: PetscCall(MatGetSize(B, NULL, &N));
618: PetscCallExternal(HYPRE_ParMultiVectorCreate, hypre_ParCSRMatrixComm(par_matrix), hypre_ParCSRMatrixGlobalNumRows(par_matrix), hypre_ParCSRMatrixRowStarts(par_matrix), N, &hb);
619: PetscCallExternal(HYPRE_ParMultiVectorCreate, hypre_ParCSRMatrixComm(par_matrix), hypre_ParCSRMatrixGlobalNumRows(par_matrix), hypre_ParCSRMatrixRowStarts(par_matrix), N, &hx);
620: PetscCall(MatZeroEntries(X));
621: PetscCall(MatDenseGetArrayReadAndMemType(B, &b, &type));
622: PetscCall(MatDenseGetLDA(B, &lda));
623: 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);
624: PetscCall(MatDenseGetLDA(X, &lda));
625: 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);
626: x_local = hypre_ParVectorLocalVector(hb);
627: PetscCallExternal(hypre_SeqVectorSetDataOwner, x_local, 0);
628: hypre_VectorData(x_local) = (HYPRE_Complex *)b;
629: PetscCall(MatDenseGetArrayWriteAndMemType(X, &x, NULL));
630: x_local = hypre_ParVectorLocalVector(hx);
631: PetscCallExternal(hypre_SeqVectorSetDataOwner, x_local, 0);
632: hypre_VectorData(x_local) = (HYPRE_Complex *)x;
633: PetscCallExternal(hypre_ParVectorInitialize_v2, hb, type == PETSC_MEMTYPE_HOST ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE);
634: PetscCallExternal(hypre_ParVectorInitialize_v2, hx, type == PETSC_MEMTYPE_HOST ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE);
635: PetscStackCallExternalVoid(
636: "Hypre solve", do {
637: HYPRE_Int hierr = (*jac->solve)(jac->hsolver, par_matrix, hb, hx);
638: if (hierr) {
639: PetscCheck(hierr == HYPRE_ERROR_CONV, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in HYPRE solver, error code %d", (int)hierr);
640: HYPRE_ClearAllErrors();
641: }
642: } while (0));
643: PetscCallExternal(HYPRE_ParVectorDestroy, hb);
644: PetscCallExternal(HYPRE_ParVectorDestroy, hx);
645: PetscCall(MatDenseRestoreArrayReadAndMemType(B, &b));
646: PetscCall(MatDenseRestoreArrayWriteAndMemType(X, &x));
647: PetscFunctionReturn(PETSC_SUCCESS);
648: }
650: static PetscErrorCode PCReset_HYPRE(PC pc)
651: {
652: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
654: PetscFunctionBegin;
655: PetscCall(MatDestroy(&jac->hpmat));
656: PetscCall(MatDestroy(&jac->G));
657: PetscCall(MatDestroy(&jac->C));
658: PetscCall(MatDestroy(&jac->alpha_Poisson));
659: PetscCall(MatDestroy(&jac->beta_Poisson));
660: PetscCall(MatDestroy(&jac->RT_PiFull));
661: PetscCall(MatDestroy(&jac->RT_Pi[0]));
662: PetscCall(MatDestroy(&jac->RT_Pi[1]));
663: PetscCall(MatDestroy(&jac->RT_Pi[2]));
664: PetscCall(MatDestroy(&jac->ND_PiFull));
665: PetscCall(MatDestroy(&jac->ND_Pi[0]));
666: PetscCall(MatDestroy(&jac->ND_Pi[1]));
667: PetscCall(MatDestroy(&jac->ND_Pi[2]));
668: PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[0]));
669: PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[1]));
670: PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[2]));
671: PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[0]));
672: PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[1]));
673: PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[2]));
674: PetscCall(VecHYPRE_IJVectorDestroy(&jac->interior));
675: PetscCall(PCHYPREResetNearNullSpace_Private(pc));
676: jac->ams_beta_is_zero = PETSC_FALSE;
677: jac->ams_beta_is_zero_part = PETSC_FALSE;
678: jac->dim = 0;
679: PetscFunctionReturn(PETSC_SUCCESS);
680: }
682: static PetscErrorCode PCDestroy_HYPRE(PC pc)
683: {
684: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
686: PetscFunctionBegin;
687: PetscCall(PCReset_HYPRE(pc));
688: if (jac->destroy) PetscCallExternal(jac->destroy, jac->hsolver);
689: PetscCall(PetscFree(jac->hypre_type));
690: if (jac->comm_hypre != MPI_COMM_NULL) PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
691: PetscCall(PetscFree(pc->data));
693: PetscCall(PetscObjectChangeTypeName((PetscObject)pc, 0));
694: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetType_C", NULL));
695: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetType_C", NULL));
696: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteGradient_C", NULL));
697: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteCurl_C", NULL));
698: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetInterpolations_C", NULL));
699: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetConstantEdgeVectors_C", NULL));
700: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetPoissonMatrix_C", NULL));
701: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetEdgeConstantVectors_C", NULL));
702: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREAMSSetInteriorNodes_C", NULL));
703: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", NULL));
704: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", NULL));
705: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetCFMarkers_C", NULL));
706: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinSetMatProductAlgorithm_C", NULL));
707: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinGetMatProductAlgorithm_C", NULL));
708: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", NULL));
709: PetscFunctionReturn(PETSC_SUCCESS);
710: }
712: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PC pc, PetscOptionItems PetscOptionsObject)
713: {
714: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
715: PetscBool flag;
717: PetscFunctionBegin;
718: PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE Pilut Options");
719: PetscCall(PetscOptionsInt("-pc_hypre_pilut_maxiter", "Number of iterations", "None", jac->maxiter, &jac->maxiter, &flag));
720: if (flag) PetscCallExternal(HYPRE_ParCSRPilutSetMaxIter, jac->hsolver, jac->maxiter);
721: PetscCall(PetscOptionsReal("-pc_hypre_pilut_tol", "Drop tolerance", "None", jac->tol, &jac->tol, &flag));
722: if (flag) PetscCallExternal(HYPRE_ParCSRPilutSetDropTolerance, jac->hsolver, jac->tol);
723: PetscCall(PetscOptionsInt("-pc_hypre_pilut_factorrowsize", "FactorRowSize", "None", jac->factorrowsize, &jac->factorrowsize, &flag));
724: if (flag) PetscCallExternal(HYPRE_ParCSRPilutSetFactorRowSize, jac->hsolver, jac->factorrowsize);
725: PetscOptionsHeadEnd();
726: PetscFunctionReturn(PETSC_SUCCESS);
727: }
729: static PetscErrorCode PCView_HYPRE_Pilut(PC pc, PetscViewer viewer)
730: {
731: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
732: PetscBool isascii;
734: PetscFunctionBegin;
735: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
736: if (isascii) {
737: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE Pilut preconditioning\n"));
738: if (jac->maxiter != PETSC_DEFAULT) {
739: PetscCall(PetscViewerASCIIPrintf(viewer, " maximum number of iterations %" PetscInt_FMT "\n", jac->maxiter));
740: } else {
741: PetscCall(PetscViewerASCIIPrintf(viewer, " default maximum number of iterations \n"));
742: }
743: if (jac->tol != PETSC_DEFAULT) {
744: PetscCall(PetscViewerASCIIPrintf(viewer, " drop tolerance %g\n", (double)jac->tol));
745: } else {
746: PetscCall(PetscViewerASCIIPrintf(viewer, " default drop tolerance \n"));
747: }
748: if (jac->factorrowsize != PETSC_DEFAULT) {
749: PetscCall(PetscViewerASCIIPrintf(viewer, " factor row size %" PetscInt_FMT "\n", jac->factorrowsize));
750: } else {
751: PetscCall(PetscViewerASCIIPrintf(viewer, " default factor row size \n"));
752: }
753: }
754: PetscFunctionReturn(PETSC_SUCCESS);
755: }
757: static const char *HYPREILUType[] = {
758: "Block-Jacobi-ILUk", "Block-Jacobi-ILUT", "", "", "", "", "", "", "", "", /* 0-9 */
759: "GMRES-ILUk", "GMRES-ILUT", "", "", "", "", "", "", "", "", /* 10-19 */
760: "NSH-ILUk", "NSH-ILUT", "", "", "", "", "", "", "", "", /* 20-29 */
761: "RAS-ILUk", "RAS-ILUT", "", "", "", "", "", "", "", "", /* 30-39 */
762: "ddPQ-GMRES-ILUk", "ddPQ-GMRES-ILUT", "", "", "", "", "", "", "", "", /* 40-49 */
763: "GMRES-ILU0" /* 50 */
764: };
766: static const char *HYPREILUIterSetup[] = {"default", "async-in-place", "async-explicit", "sync-explicit", "semisync-explicit"};
768: static PetscErrorCode PCSetFromOptions_HYPRE_ILU(PC pc, PetscOptionItems PetscOptionsObject)
769: {
770: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
771: PetscBool flg;
772: PetscInt indx;
773: PetscReal tmpdbl;
774: PetscBool tmp_truth;
776: PetscFunctionBegin;
777: PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE ILU Options");
779: /* ILU: ILU Type */
780: PetscCall(PetscOptionsEList("-pc_hypre_ilu_type", "Choose ILU Type", "None", HYPREILUType, PETSC_STATIC_ARRAY_LENGTH(HYPREILUType), HYPREILUType[0], &indx, &flg));
781: if (flg) PetscCallExternal(HYPRE_ILUSetType, jac->hsolver, indx);
783: /* ILU: ILU iterative setup type*/
784: PetscCall(PetscOptionsEList("-pc_hypre_ilu_iterative_setup_type", "Set ILU iterative setup type", "None", HYPREILUIterSetup, PETSC_STATIC_ARRAY_LENGTH(HYPREILUIterSetup), HYPREILUIterSetup[0], &indx, &flg));
785: if (flg) PetscCallExternal(HYPRE_ILUSetIterativeSetupType, jac->hsolver, indx);
787: /* ILU: ILU iterative setup option*/
788: PetscCall(PetscOptionsInt("-pc_hypre_ilu_iterative_setup_option", "Set ILU iterative setup option", "None", 0, &indx, &flg));
789: if (flg) PetscCallExternal(HYPRE_ILUSetIterativeSetupOption, jac->hsolver, indx);
791: /* ILU: ILU iterative setup maxiter */
792: PetscCall(PetscOptionsInt("-pc_hypre_ilu_iterative_setup_maxiter", "Set ILU iterative setup maximum iteration count", "None", 0, &indx, &flg));
793: if (flg) PetscCallExternal(HYPRE_ILUSetIterativeSetupMaxIter, jac->hsolver, indx);
795: /* ILU: ILU iterative setup tolerance */
796: PetscCall(PetscOptionsReal("-pc_hypre_ilu_iterative_setup_tolerance", "Set ILU iterative setup tolerance", "None", 0, &tmpdbl, &flg));
797: if (flg) PetscCallExternal(HYPRE_ILUSetIterativeSetupTolerance, jac->hsolver, tmpdbl);
799: /* ILU: ILU Print Level */
800: PetscCall(PetscOptionsInt("-pc_hypre_ilu_print_level", "Set ILU print level", "None", 0, &indx, &flg));
801: if (flg) PetscCallExternal(HYPRE_ILUSetPrintLevel, jac->hsolver, indx);
803: /* ILU: Logging */
804: PetscCall(PetscOptionsInt("-pc_hypre_ilu_logging", "Set ILU logging level", "None", 0, &indx, &flg));
805: if (flg) PetscCallExternal(HYPRE_ILUSetLogging, jac->hsolver, indx);
807: /* ILU: ILU Level */
808: PetscCall(PetscOptionsInt("-pc_hypre_ilu_level", "Set ILU level", "None", 0, &indx, &flg));
809: if (flg) PetscCallExternal(HYPRE_ILUSetLevelOfFill, jac->hsolver, indx);
811: /* ILU: ILU Max NNZ per row */
812: PetscCall(PetscOptionsInt("-pc_hypre_ilu_max_nnz_per_row", "Set maximum NNZ per row", "None", 0, &indx, &flg));
813: if (flg) PetscCallExternal(HYPRE_ILUSetMaxNnzPerRow, jac->hsolver, indx);
815: /* ILU: tolerance */
816: PetscCall(PetscOptionsReal("-pc_hypre_ilu_tol", "Tolerance for ILU", "None", 0, &tmpdbl, &flg));
817: if (flg) PetscCallExternal(HYPRE_ILUSetTol, jac->hsolver, tmpdbl);
819: /* ILU: maximum iteration count */
820: PetscCall(PetscOptionsInt("-pc_hypre_ilu_maxiter", "Set ILU max iterations", "None", 0, &indx, &flg));
821: if (flg) PetscCallExternal(HYPRE_ILUSetMaxIter, jac->hsolver, indx);
823: /* ILU: drop threshold */
824: PetscCall(PetscOptionsReal("-pc_hypre_ilu_drop_threshold", "Drop threshold for ILU", "None", 0, &tmpdbl, &flg));
825: if (flg) PetscCallExternal(HYPRE_ILUSetDropThreshold, jac->hsolver, tmpdbl);
827: /* ILU: Triangular Solve */
828: PetscCall(PetscOptionsBool("-pc_hypre_ilu_tri_solve", "Enable triangular solve", "None", PETSC_FALSE, &tmp_truth, &flg));
829: if (flg) PetscCallExternal(HYPRE_ILUSetTriSolve, jac->hsolver, tmp_truth);
831: /* ILU: Lower Jacobi iteration */
832: PetscCall(PetscOptionsInt("-pc_hypre_ilu_lower_jacobi_iters", "Set lower Jacobi iteration count", "None", 0, &indx, &flg));
833: if (flg) PetscCallExternal(HYPRE_ILUSetLowerJacobiIters, jac->hsolver, indx);
835: /* ILU: Upper Jacobi iteration */
836: PetscCall(PetscOptionsInt("-pc_hypre_ilu_upper_jacobi_iters", "Set upper Jacobi iteration count", "None", 0, &indx, &flg));
837: if (flg) PetscCallExternal(HYPRE_ILUSetUpperJacobiIters, jac->hsolver, indx);
839: /* ILU: local reordering */
840: PetscCall(PetscOptionsBool("-pc_hypre_ilu_local_reordering", "Enable local reordering", "None", PETSC_FALSE, &tmp_truth, &flg));
841: if (flg) PetscCallExternal(HYPRE_ILUSetLocalReordering, jac->hsolver, tmp_truth);
843: PetscOptionsHeadEnd();
844: PetscFunctionReturn(PETSC_SUCCESS);
845: }
847: static PetscErrorCode PCView_HYPRE_ILU(PC pc, PetscViewer viewer)
848: {
849: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
850: hypre_ParILUData *ilu_data = (hypre_ParILUData *)jac->hsolver;
851: PetscBool isascii;
852: PetscInt indx;
853: PetscReal tmpdbl;
854: PetscReal *tmpdbl3;
856: PetscFunctionBegin;
857: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
858: if (isascii) {
859: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE ILU preconditioning\n"));
860: PetscStackCallExternalVoid("hypre_ParILUDataIluType", indx = hypre_ParILUDataIluType(ilu_data));
861: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU type %s (%" PetscInt_FMT ")\n", HYPREILUType[indx], indx));
862: PetscStackCallExternalVoid("hypre_ParILUDataLfil", indx = hypre_ParILUDataLfil(ilu_data));
863: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU level %" PetscInt_FMT "\n", indx));
864: PetscStackCallExternalVoid("hypre_ParILUDataMaxIter", indx = hypre_ParILUDataMaxIter(ilu_data));
865: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU max iterations %" PetscInt_FMT "\n", indx));
866: PetscStackCallExternalVoid("hypre_ParILUDataMaxRowNnz", indx = hypre_ParILUDataMaxRowNnz(ilu_data));
867: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU max NNZ per row %" PetscInt_FMT "\n", indx));
868: PetscStackCallExternalVoid("hypre_ParILUDataTriSolve", indx = hypre_ParILUDataTriSolve(ilu_data));
869: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU triangular solve %" PetscInt_FMT "\n", indx));
870: PetscStackCallExternalVoid("hypre_ParILUDataTol", tmpdbl = hypre_ParILUDataTol(ilu_data));
871: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU tolerance %e\n", tmpdbl));
872: PetscStackCallExternalVoid("hypre_ParILUDataDroptol", tmpdbl3 = hypre_ParILUDataDroptol(ilu_data));
873: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU drop tolerance %e / %e / %e\n", tmpdbl3[0], tmpdbl3[1], tmpdbl3[2]));
874: PetscStackCallExternalVoid("hypre_ParILUDataReorderingType", indx = hypre_ParILUDataReorderingType(ilu_data));
875: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU local reordering %" PetscInt_FMT "\n", indx));
876: PetscStackCallExternalVoid("hypre_ParILUDataLowerJacobiIters", indx = hypre_ParILUDataLowerJacobiIters(ilu_data));
877: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU lower Jacobi iterations %" PetscInt_FMT "\n", indx));
878: PetscStackCallExternalVoid("hypre_ParILUDataUpperJacobiIters", indx = hypre_ParILUDataUpperJacobiIters(ilu_data));
879: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU upper Jacobi iterations %" PetscInt_FMT "\n", indx));
880: PetscStackCallExternalVoid("hypre_ParILUDataPrintLevel", indx = hypre_ParILUDataPrintLevel(ilu_data));
881: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU print level %" PetscInt_FMT "\n", indx));
882: PetscStackCallExternalVoid("hypre_ParILUDataLogging", indx = hypre_ParILUDataLogging(ilu_data));
883: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU logging level %" PetscInt_FMT "\n", indx));
884: PetscStackCallExternalVoid("hypre_ParILUDataIterativeSetupType", indx = hypre_ParILUDataIterativeSetupType(ilu_data));
885: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup type %s (%" PetscInt_FMT ")\n", HYPREILUIterSetup[indx], indx));
886: PetscStackCallExternalVoid("hypre_ParILUDataIterativeSetupOption", indx = hypre_ParILUDataIterativeSetupOption(ilu_data));
887: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup option %" PetscInt_FMT "\n", indx));
888: PetscStackCallExternalVoid("hypre_ParILUDataIterativeSetupMaxIter", indx = hypre_ParILUDataIterativeSetupMaxIter(ilu_data));
889: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup max iterations %" PetscInt_FMT "\n", indx));
890: PetscStackCallExternalVoid("hypre_ParILUDataIterativeSetupTolerance", tmpdbl = hypre_ParILUDataIterativeSetupTolerance(ilu_data));
891: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup tolerance %e\n", tmpdbl));
892: }
893: PetscFunctionReturn(PETSC_SUCCESS);
894: }
896: static PetscErrorCode PCSetFromOptions_HYPRE_Euclid(PC pc, PetscOptionItems PetscOptionsObject)
897: {
898: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
899: PetscBool flag, eu_bj = jac->eu_bj ? PETSC_TRUE : PETSC_FALSE;
901: PetscFunctionBegin;
902: PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE Euclid Options");
903: PetscCall(PetscOptionsInt("-pc_hypre_euclid_level", "Factorization levels", "None", jac->eu_level, &jac->eu_level, &flag));
904: if (flag) PetscCallExternal(HYPRE_EuclidSetLevel, jac->hsolver, jac->eu_level);
906: PetscCall(PetscOptionsReal("-pc_hypre_euclid_droptolerance", "Drop tolerance for ILU(k) in Euclid", "None", jac->eu_droptolerance, &jac->eu_droptolerance, &flag));
907: if (flag) {
908: PetscMPIInt size;
910: PetscCallMPI(MPI_Comm_size(PetscObjectComm((PetscObject)pc), &size));
911: PetscCheck(size == 1, PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "hypre's Euclid does not support a parallel drop tolerance");
912: PetscCallExternal(HYPRE_EuclidSetILUT, jac->hsolver, jac->eu_droptolerance);
913: }
915: PetscCall(PetscOptionsBool("-pc_hypre_euclid_bj", "Use Block Jacobi for ILU in Euclid", "None", eu_bj, &eu_bj, &flag));
916: if (flag) {
917: jac->eu_bj = eu_bj ? 1 : 0;
918: PetscCallExternal(HYPRE_EuclidSetBJ, jac->hsolver, jac->eu_bj);
919: }
920: PetscOptionsHeadEnd();
921: PetscFunctionReturn(PETSC_SUCCESS);
922: }
924: static PetscErrorCode PCView_HYPRE_Euclid(PC pc, PetscViewer viewer)
925: {
926: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
927: PetscBool isascii;
929: PetscFunctionBegin;
930: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
931: if (isascii) {
932: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE Euclid preconditioning\n"));
933: if (jac->eu_level != PETSC_DEFAULT) {
934: PetscCall(PetscViewerASCIIPrintf(viewer, " factorization levels %" PetscInt_FMT "\n", jac->eu_level));
935: } else {
936: PetscCall(PetscViewerASCIIPrintf(viewer, " default factorization levels \n"));
937: }
938: PetscCall(PetscViewerASCIIPrintf(viewer, " drop tolerance %g\n", (double)jac->eu_droptolerance));
939: PetscCall(PetscViewerASCIIPrintf(viewer, " use Block-Jacobi? %" PetscInt_FMT "\n", jac->eu_bj));
940: }
941: PetscFunctionReturn(PETSC_SUCCESS);
942: }
944: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc, Vec b, Vec x)
945: {
946: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
947: Mat_HYPRE *hjac = (Mat_HYPRE *)jac->hpmat->data;
948: HYPRE_ParCSRMatrix hmat;
949: HYPRE_ParVector jbv, jxv;
951: PetscFunctionBegin;
952: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
953: PetscCall(VecSet(x, 0.0));
954: PetscCall(VecHYPRE_IJVectorPushVecRead(hjac->b, b));
955: PetscCall(VecHYPRE_IJVectorPushVecWrite(hjac->x, x));
957: PetscCallExternal(HYPRE_IJMatrixGetObject, hjac->ij, (void **)&hmat);
958: PetscCallExternal(HYPRE_IJVectorGetObject, hjac->b->ij, (void **)&jbv);
959: PetscCallExternal(HYPRE_IJVectorGetObject, hjac->x->ij, (void **)&jxv);
961: PetscStackCallExternalVoid(
962: "Hypre Transpose solve", do {
963: HYPRE_Int hierr = HYPRE_BoomerAMGSolveT(jac->hsolver, hmat, jbv, jxv);
964: if (hierr) {
965: /* error code of 1 in BoomerAMG merely means convergence not achieved */
966: PetscCheck(hierr == 1, PETSC_COMM_SELF, PETSC_ERR_LIB, "Error in HYPRE solver, error code %d", (int)hierr);
967: HYPRE_ClearAllErrors();
968: }
969: } while (0));
971: PetscCall(VecHYPRE_IJVectorPopVec(hjac->x));
972: PetscCall(VecHYPRE_IJVectorPopVec(hjac->b));
973: PetscFunctionReturn(PETSC_SUCCESS);
974: }
976: static PetscErrorCode PCMGGalerkinGetMatProductAlgorithm_HYPRE_BoomerAMG(PC pc, const char *spgemm[])
977: {
978: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
980: PetscFunctionBegin;
982: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
983: *spgemm = jac->spgemm_type;
984: #endif
985: PetscFunctionReturn(PETSC_SUCCESS);
986: }
988: static const char *HYPREBoomerAMGCycleType[] = {"", "V", "W"};
989: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP", "Ruge-Stueben", "", "modifiedRuge-Stueben", "", "", "Falgout", "", "PMIS", "", "HMIS"};
990: static const char *HYPREBoomerAMGMeasureType[] = {"local", "global"};
991: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
992: static const char *HYPREBoomerAMGSmoothType[] = {"ILU", "Schwarz-smoothers", "Pilut", "ParaSails", "Euclid"};
993: 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"};
994: 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"};
996: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PC pc, PetscOptionItems PetscOptionsObject)
997: {
998: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
999: PetscInt bs, n, indx, level;
1000: PetscBool flg, tmp_truth;
1001: PetscReal tmpdbl, twodbl[2];
1002: const char *symtlist[] = {"nonsymmetric", "SPD", "nonsymmetric,SPD"};
1004: PetscFunctionBegin;
1005: PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE BoomerAMG Options");
1006: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_cycle_type", "Cycle type", "None", HYPREBoomerAMGCycleType + 1, 2, HYPREBoomerAMGCycleType[jac->cycletype], &indx, &flg));
1007: if (flg) {
1008: jac->cycletype = indx + 1;
1009: PetscCallExternal(HYPRE_BoomerAMGSetCycleType, jac->hsolver, jac->cycletype);
1010: }
1011: PetscCall(PetscOptionsBoundedInt("-pc_hypre_boomeramg_max_levels", "Number of levels (of grids) allowed", "None", jac->maxlevels, &jac->maxlevels, &flg, 2));
1012: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetMaxLevels, jac->hsolver, jac->maxlevels);
1013: PetscCall(PetscOptionsBoundedInt("-pc_hypre_boomeramg_max_iter", "Maximum iterations used PER hypre call", "None", jac->maxiter, &jac->maxiter, &flg, 1));
1014: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, jac->maxiter);
1015: 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));
1016: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, jac->tol);
1017: bs = 1;
1018: if (pc->pmat) PetscCall(MatGetBlockSize(pc->pmat, &bs));
1019: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_numfunctions", "Number of functions", "HYPRE_BoomerAMGSetNumFunctions", bs, &bs, &flg));
1020: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetNumFunctions, jac->hsolver, bs);
1022: PetscCall(PetscOptionsBoundedReal("-pc_hypre_boomeramg_truncfactor", "Truncation factor for interpolation (0=no truncation)", "None", jac->truncfactor, &jac->truncfactor, &flg, 0.0));
1023: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetTruncFactor, jac->hsolver, jac->truncfactor);
1025: PetscCall(PetscOptionsBoundedInt("-pc_hypre_boomeramg_P_max", "Max elements per row for interpolation operator (0=unlimited)", "None", jac->pmax, &jac->pmax, &flg, 0));
1026: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetPMaxElmts, jac->hsolver, jac->pmax);
1028: PetscCall(PetscOptionsRangeInt("-pc_hypre_boomeramg_agg_nl", "Number of levels of aggressive coarsening", "None", jac->agg_nl, &jac->agg_nl, &flg, 0, jac->maxlevels));
1029: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetAggNumLevels, jac->hsolver, jac->agg_nl);
1031: 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));
1032: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetNumPaths, jac->hsolver, jac->agg_num_paths);
1034: PetscCall(PetscOptionsBoundedReal("-pc_hypre_boomeramg_strong_threshold", "Threshold for being strongly connected", "None", jac->strongthreshold, &jac->strongthreshold, &flg, 0.0));
1035: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetStrongThreshold, jac->hsolver, jac->strongthreshold);
1036: PetscCall(PetscOptionsRangeReal("-pc_hypre_boomeramg_max_row_sum", "Maximum row sum", "None", jac->maxrowsum, &jac->maxrowsum, &flg, 0.0, 1.0));
1037: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetMaxRowSum, jac->hsolver, jac->maxrowsum);
1039: /* Grid sweeps */
1040: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all", "Number of sweeps for the up and down grid levels", "None", jac->gridsweeps[0], &indx, &flg));
1041: if (flg) {
1042: /* modify the jac structure so we can view the updated options with PC_View */
1043: jac->gridsweeps[0] = indx;
1044: jac->gridsweeps[1] = indx;
1045: /*defaults coarse to 1 */
1046: jac->gridsweeps[2] = 1;
1047: }
1048: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen", "Use a nodal based coarsening 1-6", "HYPRE_BoomerAMGSetNodal", jac->nodal_coarsening, &jac->nodal_coarsening, &flg));
1049: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetNodal, jac->hsolver, jac->nodal_coarsening);
1050: 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));
1051: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetNodalDiag, jac->hsolver, jac->nodal_coarsening_diag);
1052: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_variant", "Variant of algorithm 1-3", "HYPRE_BoomerAMGSetInterpVecVariant", jac->vec_interp_variant, &jac->vec_interp_variant, &flg));
1053: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetInterpVecVariant, jac->hsolver, jac->vec_interp_variant);
1054: 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));
1055: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetInterpVecQMax, jac->hsolver, jac->vec_interp_qmax);
1056: 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));
1057: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetSmoothInterpVectors, jac->hsolver, jac->vec_interp_smooth);
1058: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_interp_refine", "Preprocess the interpolation matrix through iterative weight refinement", "HYPRE_BoomerAMGSetInterpRefine", jac->interp_refine, &jac->interp_refine, &flg));
1059: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetInterpRefine, jac->hsolver, jac->interp_refine);
1060: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down", "Number of sweeps for the down cycles", "None", jac->gridsweeps[0], &indx, &flg));
1061: if (flg) {
1062: PetscCallExternal(HYPRE_BoomerAMGSetCycleNumSweeps, jac->hsolver, indx, 1);
1063: jac->gridsweeps[0] = indx;
1064: }
1065: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up", "Number of sweeps for the up cycles", "None", jac->gridsweeps[1], &indx, &flg));
1066: if (flg) {
1067: PetscCallExternal(HYPRE_BoomerAMGSetCycleNumSweeps, jac->hsolver, indx, 2);
1068: jac->gridsweeps[1] = indx;
1069: }
1070: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse", "Number of sweeps for the coarse level", "None", jac->gridsweeps[2], &indx, &flg));
1071: if (flg) {
1072: PetscCallExternal(HYPRE_BoomerAMGSetCycleNumSweeps, jac->hsolver, indx, 3);
1073: jac->gridsweeps[2] = indx;
1074: }
1076: /* Smooth type */
1077: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_smooth_type", "Enable more complex smoothers", "None", HYPREBoomerAMGSmoothType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGSmoothType), HYPREBoomerAMGSmoothType[0], &indx, &flg));
1078: if (flg) {
1079: jac->smoothtype = indx;
1080: PetscCallExternal(HYPRE_BoomerAMGSetSmoothType, jac->hsolver, indx + 5);
1081: jac->smoothnumlevels = 25;
1082: PetscCallExternal(HYPRE_BoomerAMGSetSmoothNumLevels, jac->hsolver, 25);
1083: }
1085: /* Number of smoothing levels */
1086: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_levels", "Number of levels on which more complex smoothers are used", "None", 25, &indx, &flg));
1087: if (flg && (jac->smoothtype != -1)) {
1088: jac->smoothnumlevels = indx;
1089: PetscCallExternal(HYPRE_BoomerAMGSetSmoothNumLevels, jac->hsolver, indx);
1090: }
1092: /* Smooth num sweeps */
1093: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_sweeps", "Set number of smoother sweeps", "None", 1, &indx, &flg));
1094: if (flg && indx > 0) {
1095: jac->smoothsweeps = indx;
1096: PetscCallExternal(HYPRE_BoomerAMGSetSmoothNumSweeps, jac->hsolver, indx);
1097: }
1099: /* ILU: ILU Type */
1100: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_ilu_type", "Choose ILU Type", "None", HYPREILUType, PETSC_STATIC_ARRAY_LENGTH(HYPREILUType), HYPREILUType[0], &indx, &flg));
1101: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetILUType, jac->hsolver, indx);
1103: /* ILU: ILU iterative setup type*/
1104: 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));
1105: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetILUIterSetupType, jac->hsolver, indx);
1107: /* ILU: ILU iterative setup option*/
1108: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_iterative_setup_option", "Set ILU iterative setup option", "None", 0, &indx, &flg));
1109: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetILUIterSetupOption, jac->hsolver, indx);
1111: /* ILU: ILU iterative setup maxiter */
1112: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_iterative_setup_maxiter", "Set ILU iterative setup maximum iteration count", "None", 0, &indx, &flg));
1113: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetILUIterSetupMaxIter, jac->hsolver, indx);
1115: /* ILU: ILU iterative setup tolerance */
1116: PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_ilu_iterative_setup_tolerance", "Set ILU iterative setup tolerance", "None", 0, &tmpdbl, &flg));
1117: if (flg) PetscCallExternal(hypre_BoomerAMGSetILUIterSetupTolerance, jac->hsolver, tmpdbl);
1119: /* ILU: ILU Print Level */
1120: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_print_level", "Set ILU print level", "None", 0, &indx, &flg));
1121: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetPrintLevel, jac->hsolver, indx);
1123: /* ILU: Logging */
1124: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_logging", "Set ILU logging level", "None", 0, &indx, &flg));
1125: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetLogging, jac->hsolver, indx);
1127: /* ILU: ILU Level */
1128: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_level", "Set ILU level", "None", 0, &indx, &flg));
1129: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetILULevel, jac->hsolver, indx);
1131: /* ILU: ILU Max NNZ per row */
1132: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_max_nnz_per_row", "Set maximum NNZ per row", "None", 0, &indx, &flg));
1133: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetILUMaxRowNnz, jac->hsolver, indx);
1135: /* ILU: maximum iteration count */
1136: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_maxiter", "Set ILU max iterations", "None", 0, &indx, &flg));
1137: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetILUMaxIter, jac->hsolver, indx);
1139: /* ILU: drop threshold */
1140: PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_ilu_drop_tol", "Drop tolerance for ILU", "None", 0, &tmpdbl, &flg));
1141: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetILUDroptol, jac->hsolver, tmpdbl);
1143: /* ILU: Triangular Solve */
1144: PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_ilu_tri_solve", "Enable triangular solve", "None", PETSC_FALSE, &tmp_truth, &flg));
1145: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetILUTriSolve, jac->hsolver, tmp_truth);
1147: /* ILU: Lower Jacobi iteration */
1148: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_lower_jacobi_iters", "Set lower Jacobi iteration count", "None", 0, &indx, &flg));
1149: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetILULowerJacobiIters, jac->hsolver, indx);
1151: /* ILU: Upper Jacobi iteration */
1152: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_ilu_upper_jacobi_iters", "Set upper Jacobi iteration count", "None", 0, &indx, &flg));
1153: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetILUUpperJacobiIters, jac->hsolver, indx);
1155: /* ILU: local reordering */
1156: PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_ilu_local_reordering", "Enable local reordering", "None", PETSC_FALSE, &tmp_truth, &flg));
1157: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetILULocalReordering, jac->hsolver, tmp_truth);
1159: /* Number of levels for ILU(k) for Euclid */
1160: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_eu_level", "Number of levels for ILU(k) in Euclid smoother", "None", 0, &indx, &flg));
1161: if (flg && (jac->smoothtype == 4)) {
1162: jac->eu_level = indx;
1163: PetscCallExternal(HYPRE_BoomerAMGSetEuLevel, jac->hsolver, indx);
1164: }
1166: /* Filter for ILU(k) for Euclid */
1167: PetscReal droptolerance;
1168: PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_eu_droptolerance", "Drop tolerance for ILU(k) in Euclid smoother", "None", 0, &droptolerance, &flg));
1169: if (flg && (jac->smoothtype == 4)) {
1170: jac->eu_droptolerance = droptolerance;
1171: PetscCallExternal(HYPRE_BoomerAMGSetEuLevel, jac->hsolver, droptolerance);
1172: }
1174: /* Use Block Jacobi ILUT for Euclid */
1175: PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_eu_bj", "Use Block Jacobi for ILU in Euclid smoother?", "None", PETSC_FALSE, &tmp_truth, &flg));
1176: if (flg && (jac->smoothtype == 4)) {
1177: jac->eu_bj = tmp_truth;
1178: PetscCallExternal(HYPRE_BoomerAMGSetEuBJ, jac->hsolver, jac->eu_bj);
1179: }
1181: /* Relax type */
1182: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all", "Relax type for the up and down cycles", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), HYPREBoomerAMGRelaxType[6], &indx, &flg));
1183: if (flg) {
1184: jac->relaxtype[0] = jac->relaxtype[1] = indx;
1185: PetscCallExternal(HYPRE_BoomerAMGSetRelaxType, jac->hsolver, indx);
1186: /* by default, coarse type set to 9 */
1187: jac->relaxtype[2] = 9;
1188: PetscCallExternal(HYPRE_BoomerAMGSetCycleRelaxType, jac->hsolver, 9, 3);
1189: }
1190: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down", "Relax type for the down cycles", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), HYPREBoomerAMGRelaxType[6], &indx, &flg));
1191: if (flg) {
1192: jac->relaxtype[0] = indx;
1193: PetscCallExternal(HYPRE_BoomerAMGSetCycleRelaxType, jac->hsolver, indx, 1);
1194: }
1195: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up", "Relax type for the up cycles", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), HYPREBoomerAMGRelaxType[6], &indx, &flg));
1196: if (flg) {
1197: jac->relaxtype[1] = indx;
1198: PetscCallExternal(HYPRE_BoomerAMGSetCycleRelaxType, jac->hsolver, indx, 2);
1199: }
1200: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse", "Relax type on coarse grid", "None", HYPREBoomerAMGRelaxType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGRelaxType), HYPREBoomerAMGRelaxType[9], &indx, &flg));
1201: if (flg) {
1202: jac->relaxtype[2] = indx;
1203: PetscCallExternal(HYPRE_BoomerAMGSetCycleRelaxType, jac->hsolver, indx, 3);
1204: }
1206: /* Relaxation Weight */
1207: 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));
1208: if (flg) {
1209: PetscCallExternal(HYPRE_BoomerAMGSetRelaxWt, jac->hsolver, tmpdbl);
1210: jac->relaxweight = tmpdbl;
1211: }
1213: n = 2;
1214: twodbl[0] = twodbl[1] = 1.0;
1215: PetscCall(PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level", "Set the relaxation weight for a particular level (weight,level)", "None", twodbl, &n, &flg));
1216: if (flg) {
1217: 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);
1218: indx = (int)PetscAbsReal(twodbl[1]);
1219: PetscCallExternal(HYPRE_BoomerAMGSetLevelRelaxWt, jac->hsolver, twodbl[0], indx);
1220: }
1222: /* Outer relaxation Weight */
1223: 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));
1224: if (flg) {
1225: PetscCallExternal(HYPRE_BoomerAMGSetOuterWt, jac->hsolver, tmpdbl);
1226: jac->outerrelaxweight = tmpdbl;
1227: }
1229: n = 2;
1230: twodbl[0] = twodbl[1] = 1.0;
1231: PetscCall(PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level", "Set the outer relaxation weight for a particular level (weight,level)", "None", twodbl, &n, &flg));
1232: if (flg) {
1233: 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);
1234: indx = (int)PetscAbsReal(twodbl[1]);
1235: PetscCallExternal(HYPRE_BoomerAMGSetLevelOuterWt, jac->hsolver, twodbl[0], indx);
1236: }
1238: /* the Relax Order */
1239: PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", PETSC_FALSE, &tmp_truth, &flg));
1241: if (flg && tmp_truth) {
1242: jac->relaxorder = 0;
1243: PetscCallExternal(HYPRE_BoomerAMGSetRelaxOrder, jac->hsolver, jac->relaxorder);
1244: }
1245: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_measure_type", "Measure type", "None", HYPREBoomerAMGMeasureType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGMeasureType), HYPREBoomerAMGMeasureType[0], &indx, &flg));
1246: if (flg) {
1247: jac->measuretype = indx;
1248: PetscCallExternal(HYPRE_BoomerAMGSetMeasureType, jac->hsolver, jac->measuretype);
1249: }
1250: /* update list length 3/07 */
1251: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type", "Coarsen type", "None", HYPREBoomerAMGCoarsenType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGCoarsenType), HYPREBoomerAMGCoarsenType[6], &indx, &flg));
1252: if (flg) {
1253: jac->coarsentype = indx;
1254: PetscCallExternal(HYPRE_BoomerAMGSetCoarsenType, jac->hsolver, jac->coarsentype);
1255: }
1257: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_max_coarse_size", "Maximum size of coarsest grid", "None", jac->maxc, &jac->maxc, &flg));
1258: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetMaxCoarseSize, jac->hsolver, jac->maxc);
1259: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_min_coarse_size", "Minimum size of coarsest grid", "None", jac->minc, &jac->minc, &flg));
1260: if (flg) PetscCallExternal(HYPRE_BoomerAMGSetMinCoarseSize, jac->hsolver, jac->minc);
1261: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
1262: // global parameter but is closely associated with BoomerAMG
1263: 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), HYPRESpgemmTypes[0], &indx, &flg));
1264: if (flg) PetscCall(PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG(pc, HYPRESpgemmTypes[indx]));
1265: #endif
1266: /* AIR */
1267: #if PETSC_PKG_HYPRE_VERSION_GE(2, 18, 0)
1268: 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));
1269: PetscCallExternal(HYPRE_BoomerAMGSetRestriction, jac->hsolver, jac->Rtype);
1270: if (jac->Rtype) {
1271: HYPRE_Int **grid_relax_points = hypre_TAlloc(HYPRE_Int *, 4, HYPRE_MEMORY_HOST);
1272: char *prerelax[256];
1273: char *postrelax[256];
1274: char stringF[2] = "F", stringC[2] = "C", stringA[2] = "A";
1275: PetscInt ns_down = 256, ns_up = 256;
1276: PetscBool matchF, matchC, matchA;
1278: 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 */
1280: PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_strongthresholdR", "Threshold for R", "None", jac->Rstrongthreshold, &jac->Rstrongthreshold, NULL));
1281: PetscCallExternal(HYPRE_BoomerAMGSetStrongThresholdR, jac->hsolver, jac->Rstrongthreshold);
1283: PetscCall(PetscOptionsReal("-pc_hypre_boomeramg_filterthresholdR", "Filter threshold for R", "None", jac->Rfilterthreshold, &jac->Rfilterthreshold, NULL));
1284: PetscCallExternal(HYPRE_BoomerAMGSetFilterThresholdR, jac->hsolver, jac->Rfilterthreshold);
1286: 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));
1287: PetscCallExternal(HYPRE_BoomerAMGSetADropTol, jac->hsolver, jac->Adroptol);
1289: 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));
1290: PetscCallExternal(HYPRE_BoomerAMGSetADropType, jac->hsolver, jac->Adroptype);
1291: PetscCall(PetscOptionsStringArray("-pc_hypre_boomeramg_prerelax", "Defines prerelax scheme", "None", prerelax, &ns_down, NULL));
1292: PetscCall(PetscOptionsStringArray("-pc_hypre_boomeramg_postrelax", "Defines postrelax scheme", "None", postrelax, &ns_up, NULL));
1293: 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");
1294: 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");
1296: grid_relax_points[0] = NULL;
1297: grid_relax_points[1] = hypre_TAlloc(HYPRE_Int, ns_down, HYPRE_MEMORY_HOST);
1298: grid_relax_points[2] = hypre_TAlloc(HYPRE_Int, ns_up, HYPRE_MEMORY_HOST);
1299: grid_relax_points[3] = hypre_TAlloc(HYPRE_Int, jac->gridsweeps[2], HYPRE_MEMORY_HOST);
1300: grid_relax_points[3][0] = 0;
1302: // set down relax scheme
1303: for (PetscInt i = 0; i < ns_down; i++) {
1304: PetscCall(PetscStrcasecmp(prerelax[i], stringF, &matchF));
1305: PetscCall(PetscStrcasecmp(prerelax[i], stringC, &matchC));
1306: PetscCall(PetscStrcasecmp(prerelax[i], stringA, &matchA));
1307: PetscCheck(matchF || matchC || matchA, PetscObjectComm((PetscObject)jac), PETSC_ERR_ARG_WRONG, "Valid argument options for -pc_hypre_boomeramg_prerelax are C, F, and A");
1308: if (matchF) grid_relax_points[1][i] = -1;
1309: else if (matchC) grid_relax_points[1][i] = 1;
1310: else if (matchA) grid_relax_points[1][i] = 0;
1311: }
1313: // set up relax scheme
1314: for (PetscInt i = 0; i < ns_up; i++) {
1315: PetscCall(PetscStrcasecmp(postrelax[i], stringF, &matchF));
1316: PetscCall(PetscStrcasecmp(postrelax[i], stringC, &matchC));
1317: PetscCall(PetscStrcasecmp(postrelax[i], stringA, &matchA));
1318: PetscCheck(matchF || matchC || matchA, PetscObjectComm((PetscObject)jac), PETSC_ERR_ARG_WRONG, "Valid argument options for -pc_hypre_boomeramg_postrelax are C, F, and A");
1319: if (matchF) grid_relax_points[2][i] = -1;
1320: else if (matchC) grid_relax_points[2][i] = 1;
1321: else if (matchA) grid_relax_points[2][i] = 0;
1322: }
1324: // set coarse relax scheme
1325: for (PetscInt i = 0; i < jac->gridsweeps[2]; i++) grid_relax_points[3][i] = 0;
1327: // Pass relax schemes to hypre
1328: PetscCallExternal(HYPRE_BoomerAMGSetGridRelaxPoints, jac->hsolver, grid_relax_points);
1330: // cleanup memory
1331: for (PetscInt i = 0; i < ns_down; i++) PetscCall(PetscFree(prerelax[i]));
1332: for (PetscInt i = 0; i < ns_up; i++) PetscCall(PetscFree(postrelax[i]));
1333: }
1334: #endif
1336: #if PETSC_PKG_HYPRE_VERSION_LE(9, 9, 9)
1337: 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);
1338: #endif
1340: /* new 3/07 */
1341: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_interp_type", "Interpolation type", "None", HYPREBoomerAMGInterpType, PETSC_STATIC_ARRAY_LENGTH(HYPREBoomerAMGInterpType), HYPREBoomerAMGInterpType[0], &indx, &flg));
1342: if (flg || jac->Rtype) {
1343: if (flg) jac->interptype = indx;
1344: PetscCallExternal(HYPRE_BoomerAMGSetInterpType, jac->hsolver, jac->interptype);
1345: }
1347: PetscCall(PetscOptionsName("-pc_hypre_boomeramg_print_statistics", "Print statistics", "None", &flg));
1348: if (flg) {
1349: level = 3;
1350: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_print_statistics", "Print statistics", "None", level, &level, NULL));
1352: jac->printstatistics = PETSC_TRUE;
1353: PetscCallExternal(HYPRE_BoomerAMGSetPrintLevel, jac->hsolver, level);
1354: }
1356: PetscCall(PetscOptionsName("-pc_hypre_boomeramg_print_debug", "Print debug information", "None", &flg));
1357: if (flg) {
1358: level = 3;
1359: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_print_debug", "Print debug information", "None", level, &level, NULL));
1361: jac->printstatistics = PETSC_TRUE;
1362: PetscCallExternal(HYPRE_BoomerAMGSetDebugFlag, jac->hsolver, level);
1363: }
1365: PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg));
1366: if (flg && tmp_truth) {
1367: PetscInt tmp_int;
1368: PetscCall(PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", jac->nodal_relax_levels, &tmp_int, &flg));
1369: if (flg) jac->nodal_relax_levels = tmp_int;
1370: PetscCallExternal(HYPRE_BoomerAMGSetSmoothType, jac->hsolver, 6);
1371: PetscCallExternal(HYPRE_BoomerAMGSetDomainType, jac->hsolver, 1);
1372: PetscCallExternal(HYPRE_BoomerAMGSetOverlap, jac->hsolver, 0);
1373: PetscCallExternal(HYPRE_BoomerAMGSetSmoothNumLevels, jac->hsolver, jac->nodal_relax_levels);
1374: }
1376: PetscCall(PetscOptionsBool("-pc_hypre_boomeramg_keeptranspose", "Avoid transpose matvecs in preconditioner application", "None", jac->keeptranspose, &jac->keeptranspose, NULL));
1377: PetscCallExternal(HYPRE_BoomerAMGSetKeepTranspose, jac->hsolver, jac->keeptranspose ? 1 : 0);
1379: /* options for ParaSails solvers */
1380: PetscCall(PetscOptionsEList("-pc_hypre_boomeramg_parasails_sym", "Symmetry of matrix and preconditioner", "None", symtlist, PETSC_STATIC_ARRAY_LENGTH(symtlist), symtlist[0], &indx, &flg));
1381: if (flg) {
1382: jac->symt = indx;
1383: PetscCallExternal(HYPRE_BoomerAMGSetSym, jac->hsolver, jac->symt);
1384: }
1386: PetscOptionsHeadEnd();
1387: PetscFunctionReturn(PETSC_SUCCESS);
1388: }
1390: 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)
1391: {
1392: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1393: HYPRE_Int oits;
1395: PetscFunctionBegin;
1396: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
1397: PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, its * jac->maxiter);
1398: PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, rtol);
1399: jac->applyrichardson = PETSC_TRUE;
1400: PetscCall(PCApply_HYPRE(pc, b, y));
1401: jac->applyrichardson = PETSC_FALSE;
1402: PetscCallExternal(HYPRE_BoomerAMGGetNumIterations, jac->hsolver, &oits);
1403: *outits = oits;
1404: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1405: else *reason = PCRICHARDSON_CONVERGED_RTOL;
1406: PetscCallExternal(HYPRE_BoomerAMGSetTol, jac->hsolver, jac->tol);
1407: PetscCallExternal(HYPRE_BoomerAMGSetMaxIter, jac->hsolver, jac->maxiter);
1408: PetscFunctionReturn(PETSC_SUCCESS);
1409: }
1411: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc, PetscViewer viewer)
1412: {
1413: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1414: hypre_ParAMGData *amg_data = (hypre_ParAMGData *)jac->hsolver;
1415: PetscBool isascii;
1416: PetscInt indx;
1417: PetscReal val;
1419: PetscFunctionBegin;
1420: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1421: if (isascii) {
1422: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE BoomerAMG preconditioning\n"));
1423: PetscCall(PetscViewerASCIIPrintf(viewer, " Cycle type %s\n", HYPREBoomerAMGCycleType[jac->cycletype]));
1424: PetscCall(PetscViewerASCIIPrintf(viewer, " Maximum number of levels %" PetscInt_FMT "\n", jac->maxlevels));
1425: PetscCall(PetscViewerASCIIPrintf(viewer, " Maximum number of iterations PER hypre call %" PetscInt_FMT "\n", jac->maxiter));
1426: PetscCall(PetscViewerASCIIPrintf(viewer, " Convergence tolerance PER hypre call %g\n", (double)jac->tol));
1427: PetscCall(PetscViewerASCIIPrintf(viewer, " Threshold for strong coupling %g\n", (double)jac->strongthreshold));
1428: PetscCall(PetscViewerASCIIPrintf(viewer, " Interpolation truncation factor %g\n", (double)jac->truncfactor));
1429: PetscCall(PetscViewerASCIIPrintf(viewer, " Interpolation: max elements per row %" PetscInt_FMT "\n", jac->pmax));
1430: if (jac->interp_refine) PetscCall(PetscViewerASCIIPrintf(viewer, " Interpolation: number of steps of weighted refinement %" PetscInt_FMT "\n", jac->interp_refine));
1431: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of levels of aggressive coarsening %" PetscInt_FMT "\n", jac->agg_nl));
1432: PetscCall(PetscViewerASCIIPrintf(viewer, " Number of paths for aggressive coarsening %" PetscInt_FMT "\n", jac->agg_num_paths));
1434: PetscCall(PetscViewerASCIIPrintf(viewer, " Maximum row sums %g\n", (double)jac->maxrowsum));
1436: PetscCall(PetscViewerASCIIPrintf(viewer, " Sweeps down %" PetscInt_FMT "\n", jac->gridsweeps[0]));
1437: PetscCall(PetscViewerASCIIPrintf(viewer, " Sweeps up %" PetscInt_FMT "\n", jac->gridsweeps[1]));
1438: PetscCall(PetscViewerASCIIPrintf(viewer, " Sweeps on coarse %" PetscInt_FMT "\n", jac->gridsweeps[2]));
1440: PetscCall(PetscViewerASCIIPrintf(viewer, " Relax down %s\n", HYPREBoomerAMGRelaxType[jac->relaxtype[0]]));
1441: PetscCall(PetscViewerASCIIPrintf(viewer, " Relax up %s\n", HYPREBoomerAMGRelaxType[jac->relaxtype[1]]));
1442: PetscCall(PetscViewerASCIIPrintf(viewer, " Relax on coarse %s\n", HYPREBoomerAMGRelaxType[jac->relaxtype[2]]));
1444: PetscCall(PetscViewerASCIIPrintf(viewer, " Relax weight (all) %g\n", (double)jac->relaxweight));
1445: PetscCall(PetscViewerASCIIPrintf(viewer, " Outer relax weight (all) %g\n", (double)jac->outerrelaxweight));
1447: PetscCall(PetscViewerASCIIPrintf(viewer, " Maximum size of coarsest grid %" PetscInt_FMT "\n", jac->maxc));
1448: PetscCall(PetscViewerASCIIPrintf(viewer, " Minimum size of coarsest grid %" PetscInt_FMT "\n", jac->minc));
1450: if (jac->relaxorder) {
1451: PetscCall(PetscViewerASCIIPrintf(viewer, " Using CF-relaxation\n"));
1452: } else {
1453: PetscCall(PetscViewerASCIIPrintf(viewer, " Not using CF-relaxation\n"));
1454: }
1455: if (jac->smoothtype != -1) {
1456: PetscCall(PetscViewerASCIIPrintf(viewer, " Smooth type %s\n", HYPREBoomerAMGSmoothType[jac->smoothtype]));
1457: PetscCall(PetscViewerASCIIPrintf(viewer, " Smooth num levels %" PetscInt_FMT "\n", jac->smoothnumlevels));
1458: PetscCall(PetscViewerASCIIPrintf(viewer, " Smooth num sweeps %" PetscInt_FMT "\n", jac->smoothsweeps));
1459: if (jac->smoothtype == 0) {
1460: PetscStackCallExternalVoid("hypre_ParAMGDataILUType", indx = hypre_ParAMGDataILUType(amg_data));
1461: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU type %s (%" PetscInt_FMT ")\n", HYPREILUType[indx], indx));
1462: PetscStackCallExternalVoid("hypre_ParAMGDataILULevel", indx = hypre_ParAMGDataILULevel(amg_data));
1463: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU level %" PetscInt_FMT "\n", indx));
1464: PetscStackCallExternalVoid("hypre_ParAMGDataILUMaxIter", indx = hypre_ParAMGDataILUMaxIter(amg_data));
1465: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU max iterations %" PetscInt_FMT "\n", indx));
1466: PetscStackCallExternalVoid("hypre_ParAMGDataILUMaxRowNnz", indx = hypre_ParAMGDataILUMaxRowNnz(amg_data));
1467: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU max NNZ per row %" PetscInt_FMT "\n", indx));
1468: PetscStackCallExternalVoid("hypre_ParAMGDataILUTriSolve", indx = hypre_ParAMGDataILUTriSolve(amg_data));
1469: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU triangular solve %" PetscInt_FMT "\n", indx));
1470: PetscStackCallExternalVoid("hypre_ParAMGDataTol", val = hypre_ParAMGDataTol(amg_data));
1471: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU tolerance %e\n", val));
1472: PetscStackCallExternalVoid("hypre_ParAMGDataILUDroptol", val = hypre_ParAMGDataILUDroptol(amg_data));
1473: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU drop tolerance %e\n", val));
1474: PetscStackCallExternalVoid("hypre_ParAMGDataILULocalReordering", indx = hypre_ParAMGDataILULocalReordering(amg_data));
1475: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU local reordering %" PetscInt_FMT "\n", indx));
1476: PetscStackCallExternalVoid("hypre_ParAMGDataILULowerJacobiIters", indx = hypre_ParAMGDataILULowerJacobiIters(amg_data));
1477: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU lower Jacobi iterations %" PetscInt_FMT "\n", indx));
1478: PetscStackCallExternalVoid("hypre_ParAMGDataILUUpperJacobiIters", indx = hypre_ParAMGDataILUUpperJacobiIters(amg_data));
1479: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU upper Jacobi iterations %" PetscInt_FMT "\n", indx));
1480: PetscStackCallExternalVoid("hypre_ParAMGDataPrintLevel", indx = hypre_ParAMGDataPrintLevel(amg_data));
1481: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU print level %" PetscInt_FMT "\n", indx));
1482: PetscStackCallExternalVoid("hypre_ParAMGDataLogging", indx = hypre_ParAMGDataLogging(amg_data));
1483: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU logging level %" PetscInt_FMT "\n", indx));
1484: PetscStackCallExternalVoid("hypre_ParAMGDataILUIterSetupType", indx = hypre_ParAMGDataILUIterSetupType(amg_data));
1485: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup type %s (%" PetscInt_FMT ")\n", HYPREILUIterSetup[indx], indx));
1486: PetscStackCallExternalVoid("hypre_ParAMGDataILUIterSetupOption", indx = hypre_ParAMGDataILUIterSetupOption(amg_data));
1487: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup option %" PetscInt_FMT "\n", indx));
1488: PetscStackCallExternalVoid("hypre_ParAMGDataILUIterSetupMaxIter", indx = hypre_ParAMGDataILUIterSetupMaxIter(amg_data));
1489: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup max iterations %" PetscInt_FMT "\n", indx));
1490: PetscStackCallExternalVoid("hypre_ParAMGDataILUIterSetupTolerance", val = hypre_ParAMGDataILUIterSetupTolerance(amg_data));
1491: PetscCall(PetscViewerASCIIPrintf(viewer, " ILU iterative setup tolerance %e\n", val));
1492: }
1493: } else {
1494: PetscCall(PetscViewerASCIIPrintf(viewer, " Not using more complex smoothers.\n"));
1495: }
1496: if (jac->smoothtype == 3) {
1497: PetscCall(PetscViewerASCIIPrintf(viewer, " Euclid ILU(k) levels %" PetscInt_FMT "\n", jac->eu_level));
1498: PetscCall(PetscViewerASCIIPrintf(viewer, " Euclid ILU(k) drop tolerance %g\n", (double)jac->eu_droptolerance));
1499: PetscCall(PetscViewerASCIIPrintf(viewer, " Euclid ILU use Block-Jacobi? %" PetscInt_FMT "\n", jac->eu_bj));
1500: }
1501: PetscCall(PetscViewerASCIIPrintf(viewer, " Measure type %s\n", HYPREBoomerAMGMeasureType[jac->measuretype]));
1502: PetscCall(PetscViewerASCIIPrintf(viewer, " Coarsen type %s\n", HYPREBoomerAMGCoarsenType[jac->coarsentype]));
1503: PetscCall(PetscViewerASCIIPrintf(viewer, " Interpolation type %s\n", jac->interptype != 100 ? HYPREBoomerAMGInterpType[jac->interptype] : "1pt"));
1504: if (jac->nodal_coarsening) PetscCall(PetscViewerASCIIPrintf(viewer, " Using nodal coarsening with HYPRE_BOOMERAMGSetNodal() %" PetscInt_FMT "\n", jac->nodal_coarsening));
1505: if (jac->vec_interp_variant) {
1506: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE_BoomerAMGSetInterpVecVariant() %" PetscInt_FMT "\n", jac->vec_interp_variant));
1507: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE_BoomerAMGSetInterpVecQMax() %" PetscInt_FMT "\n", jac->vec_interp_qmax));
1508: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE_BoomerAMGSetSmoothInterpVectors() %d\n", jac->vec_interp_smooth));
1509: }
1510: if (jac->nodal_relax) PetscCall(PetscViewerASCIIPrintf(viewer, " Using nodal relaxation via Schwarz smoothing on levels %" PetscInt_FMT "\n", jac->nodal_relax_levels));
1511: #if PETSC_PKG_HYPRE_VERSION_GE(2, 23, 0)
1512: PetscCall(PetscViewerASCIIPrintf(viewer, " SpGEMM type %s\n", jac->spgemm_type));
1513: #else
1514: PetscCall(PetscViewerASCIIPrintf(viewer, " SpGEMM type %s\n", "hypre"));
1515: #endif
1516: /* AIR */
1517: if (jac->Rtype) {
1518: PetscCall(PetscViewerASCIIPrintf(viewer, " Using approximate ideal restriction type %" PetscInt_FMT "\n", jac->Rtype));
1519: PetscCall(PetscViewerASCIIPrintf(viewer, " Threshold for R %g\n", (double)jac->Rstrongthreshold));
1520: PetscCall(PetscViewerASCIIPrintf(viewer, " Filter for R %g\n", (double)jac->Rfilterthreshold));
1521: PetscCall(PetscViewerASCIIPrintf(viewer, " A drop tolerance %g\n", (double)jac->Adroptol));
1522: PetscCall(PetscViewerASCIIPrintf(viewer, " A drop type %" PetscInt_FMT "\n", jac->Adroptype));
1523: }
1524: }
1525: PetscFunctionReturn(PETSC_SUCCESS);
1526: }
1528: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PC pc, PetscOptionItems PetscOptionsObject)
1529: {
1530: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1531: PetscInt indx;
1532: PetscBool flag;
1533: const char *symtlist[] = {"nonsymmetric", "SPD", "nonsymmetric,SPD"};
1535: PetscFunctionBegin;
1536: PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE ParaSails Options");
1537: PetscCall(PetscOptionsInt("-pc_hypre_parasails_nlevels", "Number of number of levels", "None", jac->nlevels, &jac->nlevels, 0));
1538: PetscCall(PetscOptionsReal("-pc_hypre_parasails_thresh", "Threshold", "None", jac->threshold, &jac->threshold, &flag));
1539: if (flag) PetscCallExternal(HYPRE_ParaSailsSetParams, jac->hsolver, jac->threshold, jac->nlevels);
1541: PetscCall(PetscOptionsReal("-pc_hypre_parasails_filter", "filter", "None", jac->filter, &jac->filter, &flag));
1542: if (flag) PetscCallExternal(HYPRE_ParaSailsSetFilter, jac->hsolver, jac->filter);
1544: PetscCall(PetscOptionsReal("-pc_hypre_parasails_loadbal", "Load balance", "None", jac->loadbal, &jac->loadbal, &flag));
1545: if (flag) PetscCallExternal(HYPRE_ParaSailsSetLoadbal, jac->hsolver, jac->loadbal);
1547: PetscCall(PetscOptionsBool("-pc_hypre_parasails_logging", "Print info to screen", "None", (PetscBool)jac->logging, (PetscBool *)&jac->logging, &flag));
1548: if (flag) PetscCallExternal(HYPRE_ParaSailsSetLogging, jac->hsolver, jac->logging);
1550: PetscCall(PetscOptionsBool("-pc_hypre_parasails_reuse", "Reuse nonzero pattern in preconditioner", "None", (PetscBool)jac->ruse, (PetscBool *)&jac->ruse, &flag));
1551: if (flag) PetscCallExternal(HYPRE_ParaSailsSetReuse, jac->hsolver, jac->ruse);
1553: PetscCall(PetscOptionsEList("-pc_hypre_parasails_sym", "Symmetry of matrix and preconditioner", "None", symtlist, PETSC_STATIC_ARRAY_LENGTH(symtlist), symtlist[0], &indx, &flag));
1554: if (flag) {
1555: jac->symt = indx;
1556: PetscCallExternal(HYPRE_ParaSailsSetSym, jac->hsolver, jac->symt);
1557: }
1559: PetscOptionsHeadEnd();
1560: PetscFunctionReturn(PETSC_SUCCESS);
1561: }
1563: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc, PetscViewer viewer)
1564: {
1565: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1566: PetscBool isascii;
1567: const char *symt = 0;
1569: PetscFunctionBegin;
1570: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1571: if (isascii) {
1572: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE ParaSails preconditioning\n"));
1573: PetscCall(PetscViewerASCIIPrintf(viewer, " nlevels %" PetscInt_FMT "\n", jac->nlevels));
1574: PetscCall(PetscViewerASCIIPrintf(viewer, " threshold %g\n", (double)jac->threshold));
1575: PetscCall(PetscViewerASCIIPrintf(viewer, " filter %g\n", (double)jac->filter));
1576: PetscCall(PetscViewerASCIIPrintf(viewer, " load balance %g\n", (double)jac->loadbal));
1577: PetscCall(PetscViewerASCIIPrintf(viewer, " reuse nonzero structure %s\n", PetscBools[jac->ruse]));
1578: PetscCall(PetscViewerASCIIPrintf(viewer, " print info to screen %s\n", PetscBools[jac->logging]));
1579: if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
1580: else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
1581: else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
1582: else SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_WRONG, "Unknown HYPRE ParaSails symmetric option %" PetscInt_FMT, jac->symt);
1583: PetscCall(PetscViewerASCIIPrintf(viewer, " %s\n", symt));
1584: }
1585: PetscFunctionReturn(PETSC_SUCCESS);
1586: }
1588: static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PC pc, PetscOptionItems PetscOptionsObject)
1589: {
1590: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1591: PetscInt n;
1592: PetscBool flag, flag2, flag3, flag4;
1594: PetscFunctionBegin;
1595: PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE AMS Options");
1596: PetscCall(PetscOptionsInt("-pc_hypre_ams_print_level", "Debugging output level for AMS", "None", jac->as_print, &jac->as_print, &flag));
1597: if (flag) PetscCallExternal(HYPRE_AMSSetPrintLevel, jac->hsolver, jac->as_print);
1598: 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));
1599: if (flag) PetscCallExternal(HYPRE_AMSSetMaxIter, jac->hsolver, jac->as_max_iter);
1600: PetscCall(PetscOptionsInt("-pc_hypre_ams_cycle_type", "Cycle type for AMS multigrid", "None", jac->ams_cycle_type, &jac->ams_cycle_type, &flag));
1601: if (flag) PetscCallExternal(HYPRE_AMSSetCycleType, jac->hsolver, jac->ams_cycle_type);
1602: PetscCall(PetscOptionsReal("-pc_hypre_ams_tol", "Error tolerance for AMS multigrid", "None", jac->as_tol, &jac->as_tol, &flag));
1603: if (flag) PetscCallExternal(HYPRE_AMSSetTol, jac->hsolver, jac->as_tol);
1604: PetscCall(PetscOptionsInt("-pc_hypre_ams_relax_type", "Relaxation type for AMS smoother", "None", jac->as_relax_type, &jac->as_relax_type, &flag));
1605: PetscCall(PetscOptionsInt("-pc_hypre_ams_relax_times", "Number of relaxation steps for AMS smoother", "None", jac->as_relax_times, &jac->as_relax_times, &flag2));
1606: PetscCall(PetscOptionsReal("-pc_hypre_ams_relax_weight", "Relaxation weight for AMS smoother", "None", jac->as_relax_weight, &jac->as_relax_weight, &flag3));
1607: PetscCall(PetscOptionsReal("-pc_hypre_ams_omega", "SSOR coefficient for AMS smoother", "None", jac->as_omega, &jac->as_omega, &flag4));
1608: if (flag || flag2 || flag3 || flag4) PetscCallExternal(HYPRE_AMSSetSmoothingOptions, jac->hsolver, jac->as_relax_type, jac->as_relax_times, jac->as_relax_weight, jac->as_omega);
1609: 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));
1610: n = 5;
1611: PetscCall(PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options", "AMG options for vector Poisson", "None", jac->as_amg_alpha_opts, &n, &flag2));
1612: if (flag || flag2) {
1613: PetscCallExternal(HYPRE_AMSSetAlphaAMGOptions, jac->hsolver, jac->as_amg_alpha_opts[0], /* AMG coarsen type */
1614: jac->as_amg_alpha_opts[1], /* AMG agg_levels */
1615: jac->as_amg_alpha_opts[2], /* AMG relax_type */
1616: jac->as_amg_alpha_theta, jac->as_amg_alpha_opts[3], /* AMG interp_type */
1617: jac->as_amg_alpha_opts[4]); /* AMG Pmax */
1618: }
1619: 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));
1620: n = 5;
1621: PetscCall(PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options", "AMG options for scalar Poisson solver", "None", jac->as_amg_beta_opts, &n, &flag2));
1622: if (flag || flag2) {
1623: PetscCallExternal(HYPRE_AMSSetBetaAMGOptions, jac->hsolver, jac->as_amg_beta_opts[0], /* AMG coarsen type */
1624: jac->as_amg_beta_opts[1], /* AMG agg_levels */
1625: jac->as_amg_beta_opts[2], /* AMG relax_type */
1626: jac->as_amg_beta_theta, jac->as_amg_beta_opts[3], /* AMG interp_type */
1627: jac->as_amg_beta_opts[4]); /* AMG Pmax */
1628: }
1629: 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));
1630: if (flag) { /* override HYPRE's default only if the options is used */
1631: PetscCallExternal(HYPRE_AMSSetProjectionFrequency, jac->hsolver, jac->ams_proj_freq);
1632: }
1633: PetscOptionsHeadEnd();
1634: PetscFunctionReturn(PETSC_SUCCESS);
1635: }
1637: static PetscErrorCode PCView_HYPRE_AMS(PC pc, PetscViewer viewer)
1638: {
1639: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1640: PetscBool isascii;
1642: PetscFunctionBegin;
1643: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1644: if (isascii) {
1645: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE AMS preconditioning\n"));
1646: PetscCall(PetscViewerASCIIPrintf(viewer, " subspace iterations per application %" PetscInt_FMT "\n", jac->as_max_iter));
1647: PetscCall(PetscViewerASCIIPrintf(viewer, " subspace cycle type %" PetscInt_FMT "\n", jac->ams_cycle_type));
1648: PetscCall(PetscViewerASCIIPrintf(viewer, " subspace iteration tolerance %g\n", (double)jac->as_tol));
1649: PetscCall(PetscViewerASCIIPrintf(viewer, " smoother type %" PetscInt_FMT "\n", jac->as_relax_type));
1650: PetscCall(PetscViewerASCIIPrintf(viewer, " number of smoothing steps %" PetscInt_FMT "\n", jac->as_relax_times));
1651: PetscCall(PetscViewerASCIIPrintf(viewer, " smoother weight %g\n", (double)jac->as_relax_weight));
1652: PetscCall(PetscViewerASCIIPrintf(viewer, " smoother omega %g\n", (double)jac->as_omega));
1653: if (jac->alpha_Poisson) {
1654: PetscCall(PetscViewerASCIIPrintf(viewer, " vector Poisson solver (passed in by user)\n"));
1655: } else {
1656: PetscCall(PetscViewerASCIIPrintf(viewer, " vector Poisson solver (computed) \n"));
1657: }
1658: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG coarsening type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[0]));
1659: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[1]));
1660: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG relaxation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[2]));
1661: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG interpolation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[3]));
1662: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[4]));
1663: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG strength threshold %g\n", (double)jac->as_amg_alpha_theta));
1664: if (!jac->ams_beta_is_zero) {
1665: if (jac->beta_Poisson) {
1666: PetscCall(PetscViewerASCIIPrintf(viewer, " scalar Poisson solver (passed in by user)\n"));
1667: } else {
1668: PetscCall(PetscViewerASCIIPrintf(viewer, " scalar Poisson solver (computed) \n"));
1669: }
1670: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG coarsening type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[0]));
1671: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_beta_opts[1]));
1672: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG relaxation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[2]));
1673: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG interpolation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[3]));
1674: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_beta_opts[4]));
1675: PetscCall(PetscViewerASCIIPrintf(viewer, " boomerAMG strength threshold %g\n", (double)jac->as_amg_beta_theta));
1676: 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));
1677: } else {
1678: PetscCall(PetscViewerASCIIPrintf(viewer, " scalar Poisson solver not used (zero-conductivity everywhere) \n"));
1679: }
1680: }
1681: PetscFunctionReturn(PETSC_SUCCESS);
1682: }
1684: static PetscErrorCode PCSetFromOptions_HYPRE_ADS(PC pc, PetscOptionItems PetscOptionsObject)
1685: {
1686: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1687: PetscInt n;
1688: PetscBool flag, flag2, flag3, flag4;
1690: PetscFunctionBegin;
1691: PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE ADS Options");
1692: PetscCall(PetscOptionsInt("-pc_hypre_ads_print_level", "Debugging output level for ADS", "None", jac->as_print, &jac->as_print, &flag));
1693: if (flag) PetscCallExternal(HYPRE_ADSSetPrintLevel, jac->hsolver, jac->as_print);
1694: 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));
1695: if (flag) PetscCallExternal(HYPRE_ADSSetMaxIter, jac->hsolver, jac->as_max_iter);
1696: PetscCall(PetscOptionsInt("-pc_hypre_ads_cycle_type", "Cycle type for ADS multigrid", "None", jac->ads_cycle_type, &jac->ads_cycle_type, &flag));
1697: if (flag) PetscCallExternal(HYPRE_ADSSetCycleType, jac->hsolver, jac->ads_cycle_type);
1698: PetscCall(PetscOptionsReal("-pc_hypre_ads_tol", "Error tolerance for ADS multigrid", "None", jac->as_tol, &jac->as_tol, &flag));
1699: if (flag) PetscCallExternal(HYPRE_ADSSetTol, jac->hsolver, jac->as_tol);
1700: PetscCall(PetscOptionsInt("-pc_hypre_ads_relax_type", "Relaxation type for ADS smoother", "None", jac->as_relax_type, &jac->as_relax_type, &flag));
1701: PetscCall(PetscOptionsInt("-pc_hypre_ads_relax_times", "Number of relaxation steps for ADS smoother", "None", jac->as_relax_times, &jac->as_relax_times, &flag2));
1702: PetscCall(PetscOptionsReal("-pc_hypre_ads_relax_weight", "Relaxation weight for ADS smoother", "None", jac->as_relax_weight, &jac->as_relax_weight, &flag3));
1703: PetscCall(PetscOptionsReal("-pc_hypre_ads_omega", "SSOR coefficient for ADS smoother", "None", jac->as_omega, &jac->as_omega, &flag4));
1704: if (flag || flag2 || flag3 || flag4) PetscCallExternal(HYPRE_ADSSetSmoothingOptions, jac->hsolver, jac->as_relax_type, jac->as_relax_times, jac->as_relax_weight, jac->as_omega);
1705: 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));
1706: n = 5;
1707: PetscCall(PetscOptionsIntArray("-pc_hypre_ads_ams_options", "AMG options for AMS solver inside ADS", "None", jac->as_amg_alpha_opts, &n, &flag2));
1708: 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));
1709: if (flag || flag2 || flag3) {
1710: PetscCallExternal(HYPRE_ADSSetAMSOptions, jac->hsolver, jac->ams_cycle_type, /* AMS cycle type */
1711: jac->as_amg_alpha_opts[0], /* AMG coarsen type */
1712: jac->as_amg_alpha_opts[1], /* AMG agg_levels */
1713: jac->as_amg_alpha_opts[2], /* AMG relax_type */
1714: jac->as_amg_alpha_theta, jac->as_amg_alpha_opts[3], /* AMG interp_type */
1715: jac->as_amg_alpha_opts[4]); /* AMG Pmax */
1716: }
1717: 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));
1718: n = 5;
1719: PetscCall(PetscOptionsIntArray("-pc_hypre_ads_amg_options", "AMG options for vector AMG solver inside ADS", "None", jac->as_amg_beta_opts, &n, &flag2));
1720: if (flag || flag2) {
1721: PetscCallExternal(HYPRE_ADSSetAMGOptions, jac->hsolver, jac->as_amg_beta_opts[0], /* AMG coarsen type */
1722: jac->as_amg_beta_opts[1], /* AMG agg_levels */
1723: jac->as_amg_beta_opts[2], /* AMG relax_type */
1724: jac->as_amg_beta_theta, jac->as_amg_beta_opts[3], /* AMG interp_type */
1725: jac->as_amg_beta_opts[4]); /* AMG Pmax */
1726: }
1727: PetscOptionsHeadEnd();
1728: PetscFunctionReturn(PETSC_SUCCESS);
1729: }
1731: static PetscErrorCode PCView_HYPRE_ADS(PC pc, PetscViewer viewer)
1732: {
1733: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1734: PetscBool isascii;
1736: PetscFunctionBegin;
1737: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
1738: if (isascii) {
1739: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE ADS preconditioning\n"));
1740: PetscCall(PetscViewerASCIIPrintf(viewer, " subspace iterations per application %" PetscInt_FMT "\n", jac->as_max_iter));
1741: PetscCall(PetscViewerASCIIPrintf(viewer, " subspace cycle type %" PetscInt_FMT "\n", jac->ads_cycle_type));
1742: PetscCall(PetscViewerASCIIPrintf(viewer, " subspace iteration tolerance %g\n", (double)jac->as_tol));
1743: PetscCall(PetscViewerASCIIPrintf(viewer, " smoother type %" PetscInt_FMT "\n", jac->as_relax_type));
1744: PetscCall(PetscViewerASCIIPrintf(viewer, " number of smoothing steps %" PetscInt_FMT "\n", jac->as_relax_times));
1745: PetscCall(PetscViewerASCIIPrintf(viewer, " smoother weight %g\n", (double)jac->as_relax_weight));
1746: PetscCall(PetscViewerASCIIPrintf(viewer, " smoother omega %g\n", (double)jac->as_omega));
1747: PetscCall(PetscViewerASCIIPrintf(viewer, " AMS solver using boomerAMG\n"));
1748: PetscCall(PetscViewerASCIIPrintf(viewer, " subspace cycle type %" PetscInt_FMT "\n", jac->ams_cycle_type));
1749: PetscCall(PetscViewerASCIIPrintf(viewer, " coarsening type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[0]));
1750: PetscCall(PetscViewerASCIIPrintf(viewer, " levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[1]));
1751: PetscCall(PetscViewerASCIIPrintf(viewer, " relaxation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[2]));
1752: PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation type %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[3]));
1753: PetscCall(PetscViewerASCIIPrintf(viewer, " max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_alpha_opts[4]));
1754: PetscCall(PetscViewerASCIIPrintf(viewer, " strength threshold %g\n", (double)jac->as_amg_alpha_theta));
1755: PetscCall(PetscViewerASCIIPrintf(viewer, " vector Poisson solver using boomerAMG\n"));
1756: PetscCall(PetscViewerASCIIPrintf(viewer, " coarsening type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[0]));
1757: PetscCall(PetscViewerASCIIPrintf(viewer, " levels of aggressive coarsening %" PetscInt_FMT "\n", jac->as_amg_beta_opts[1]));
1758: PetscCall(PetscViewerASCIIPrintf(viewer, " relaxation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[2]));
1759: PetscCall(PetscViewerASCIIPrintf(viewer, " interpolation type %" PetscInt_FMT "\n", jac->as_amg_beta_opts[3]));
1760: PetscCall(PetscViewerASCIIPrintf(viewer, " max nonzero elements in interpolation rows %" PetscInt_FMT "\n", jac->as_amg_beta_opts[4]));
1761: PetscCall(PetscViewerASCIIPrintf(viewer, " strength threshold %g\n", (double)jac->as_amg_beta_theta));
1762: }
1763: PetscFunctionReturn(PETSC_SUCCESS);
1764: }
1766: static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE(PC pc, Mat G)
1767: {
1768: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1769: PetscBool ishypre;
1771: PetscFunctionBegin;
1772: PetscCall(PetscObjectTypeCompare((PetscObject)G, MATHYPRE, &ishypre));
1773: if (ishypre) {
1774: PetscCall(PetscObjectReference((PetscObject)G));
1775: PetscCall(MatDestroy(&jac->G));
1776: jac->G = G;
1777: } else {
1778: PetscCall(MatDestroy(&jac->G));
1779: PetscCall(MatConvert(G, MATHYPRE, MAT_INITIAL_MATRIX, &jac->G));
1780: }
1781: PetscFunctionReturn(PETSC_SUCCESS);
1782: }
1784: /*@
1785: PCHYPRESetDiscreteGradient - Set discrete gradient matrix for `PCHYPRE` type of ams or ads
1787: Collective
1789: Input Parameters:
1790: + pc - the preconditioning context
1791: - G - the discrete gradient
1793: Level: intermediate
1795: Notes:
1796: G should have as many rows as the number of edges and as many columns as the number of vertices in the mesh
1798: Each row of G has 2 nonzeros, with column indexes being the global indexes of edge's endpoints: matrix entries are +1 and -1 depending on edge orientation
1800: Developer Notes:
1801: This automatically converts the matrix to `MATHYPRE` if it is not already of that type
1803: .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteCurl()`
1804: @*/
1805: PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G)
1806: {
1807: PetscFunctionBegin;
1810: PetscCheckSameComm(pc, 1, G, 2);
1811: PetscTryMethod(pc, "PCHYPRESetDiscreteGradient_C", (PC, Mat), (pc, G));
1812: PetscFunctionReturn(PETSC_SUCCESS);
1813: }
1815: static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C)
1816: {
1817: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1818: PetscBool ishypre;
1820: PetscFunctionBegin;
1821: PetscCall(PetscObjectTypeCompare((PetscObject)C, MATHYPRE, &ishypre));
1822: if (ishypre) {
1823: PetscCall(PetscObjectReference((PetscObject)C));
1824: PetscCall(MatDestroy(&jac->C));
1825: jac->C = C;
1826: } else {
1827: PetscCall(MatDestroy(&jac->C));
1828: PetscCall(MatConvert(C, MATHYPRE, MAT_INITIAL_MATRIX, &jac->C));
1829: }
1830: PetscFunctionReturn(PETSC_SUCCESS);
1831: }
1833: /*@
1834: PCHYPRESetDiscreteCurl - Set discrete curl matrix for `PCHYPRE` type of ads
1836: Collective
1838: Input Parameters:
1839: + pc - the preconditioning context
1840: - C - the discrete curl
1842: Level: intermediate
1844: Notes:
1845: C should have as many rows as the number of faces and as many columns as the number of edges in the mesh
1847: Each row of C has as many nonzeros as the number of edges of a face, with column indexes being the global indexes of the corresponding edge: matrix entries are +1 and -1 depending on edge orientation with respect to the face orientation
1849: Developer Notes:
1850: This automatically converts the matrix to `MATHYPRE` if it is not already of that type
1852: If this is only for `PCHYPRE` type of ads it should be called `PCHYPREADSSetDiscreteCurl()`
1854: .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`
1855: @*/
1856: PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C)
1857: {
1858: PetscFunctionBegin;
1861: PetscCheckSameComm(pc, 1, C, 2);
1862: PetscTryMethod(pc, "PCHYPRESetDiscreteCurl_C", (PC, Mat), (pc, C));
1863: PetscFunctionReturn(PETSC_SUCCESS);
1864: }
1866: static PetscErrorCode PCHYPRESetInterpolations_HYPRE(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[])
1867: {
1868: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1869: PetscBool ishypre;
1870: PetscInt i;
1872: PetscFunctionBegin;
1873: PetscCall(MatDestroy(&jac->RT_PiFull));
1874: PetscCall(MatDestroy(&jac->ND_PiFull));
1875: for (i = 0; i < 3; ++i) {
1876: PetscCall(MatDestroy(&jac->RT_Pi[i]));
1877: PetscCall(MatDestroy(&jac->ND_Pi[i]));
1878: }
1880: jac->dim = dim;
1881: if (RT_PiFull) {
1882: PetscCall(PetscObjectTypeCompare((PetscObject)RT_PiFull, MATHYPRE, &ishypre));
1883: if (ishypre) {
1884: PetscCall(PetscObjectReference((PetscObject)RT_PiFull));
1885: jac->RT_PiFull = RT_PiFull;
1886: } else {
1887: PetscCall(MatConvert(RT_PiFull, MATHYPRE, MAT_INITIAL_MATRIX, &jac->RT_PiFull));
1888: }
1889: }
1890: if (RT_Pi) {
1891: for (i = 0; i < dim; ++i) {
1892: if (RT_Pi[i]) {
1893: PetscCall(PetscObjectTypeCompare((PetscObject)RT_Pi[i], MATHYPRE, &ishypre));
1894: if (ishypre) {
1895: PetscCall(PetscObjectReference((PetscObject)RT_Pi[i]));
1896: jac->RT_Pi[i] = RT_Pi[i];
1897: } else {
1898: PetscCall(MatConvert(RT_Pi[i], MATHYPRE, MAT_INITIAL_MATRIX, &jac->RT_Pi[i]));
1899: }
1900: }
1901: }
1902: }
1903: if (ND_PiFull) {
1904: PetscCall(PetscObjectTypeCompare((PetscObject)ND_PiFull, MATHYPRE, &ishypre));
1905: if (ishypre) {
1906: PetscCall(PetscObjectReference((PetscObject)ND_PiFull));
1907: jac->ND_PiFull = ND_PiFull;
1908: } else {
1909: PetscCall(MatConvert(ND_PiFull, MATHYPRE, MAT_INITIAL_MATRIX, &jac->ND_PiFull));
1910: }
1911: }
1912: if (ND_Pi) {
1913: for (i = 0; i < dim; ++i) {
1914: if (ND_Pi[i]) {
1915: PetscCall(PetscObjectTypeCompare((PetscObject)ND_Pi[i], MATHYPRE, &ishypre));
1916: if (ishypre) {
1917: PetscCall(PetscObjectReference((PetscObject)ND_Pi[i]));
1918: jac->ND_Pi[i] = ND_Pi[i];
1919: } else {
1920: PetscCall(MatConvert(ND_Pi[i], MATHYPRE, MAT_INITIAL_MATRIX, &jac->ND_Pi[i]));
1921: }
1922: }
1923: }
1924: }
1925: PetscFunctionReturn(PETSC_SUCCESS);
1926: }
1928: /*@
1929: PCHYPRESetInterpolations - Set interpolation matrices for `PCHYPRE` type of ams or ads
1931: Collective
1933: Input Parameters:
1934: + pc - the preconditioning context
1935: . dim - the dimension of the problem, only used in AMS
1936: . RT_PiFull - Raviart-Thomas interpolation matrix
1937: . RT_Pi - x/y/z component of Raviart-Thomas interpolation matrix
1938: . ND_PiFull - Nedelec interpolation matrix
1939: - ND_Pi - x/y/z component of Nedelec interpolation matrix
1941: Level: intermediate
1943: Notes:
1944: For AMS, only Nedelec interpolation matrices are needed, the Raviart-Thomas interpolation matrices can be set to NULL.
1946: For ADS, both type of interpolation matrices are needed.
1948: Developer Notes:
1949: This automatically converts the matrix to `MATHYPRE` if it is not already of that type
1951: .seealso: [](ch_ksp), `PCHYPRE`
1952: @*/
1953: PetscErrorCode PCHYPRESetInterpolations(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[])
1954: {
1955: PetscInt i;
1957: PetscFunctionBegin;
1959: if (RT_PiFull) {
1961: PetscCheckSameComm(pc, 1, RT_PiFull, 3);
1962: }
1963: if (RT_Pi) {
1964: PetscAssertPointer(RT_Pi, 4);
1965: for (i = 0; i < dim; ++i) {
1966: if (RT_Pi[i]) {
1968: PetscCheckSameComm(pc, 1, RT_Pi[i], 4);
1969: }
1970: }
1971: }
1972: if (ND_PiFull) {
1974: PetscCheckSameComm(pc, 1, ND_PiFull, 5);
1975: }
1976: if (ND_Pi) {
1977: PetscAssertPointer(ND_Pi, 6);
1978: for (i = 0; i < dim; ++i) {
1979: if (ND_Pi[i]) {
1981: PetscCheckSameComm(pc, 1, ND_Pi[i], 6);
1982: }
1983: }
1984: }
1985: PetscTryMethod(pc, "PCHYPRESetInterpolations_C", (PC, PetscInt, Mat, Mat[], Mat, Mat[]), (pc, dim, RT_PiFull, RT_Pi, ND_PiFull, ND_Pi));
1986: PetscFunctionReturn(PETSC_SUCCESS);
1987: }
1989: static PetscErrorCode PCHYPRESetPoissonMatrix_HYPRE(PC pc, Mat A, PetscBool isalpha)
1990: {
1991: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
1992: PetscBool ishypre;
1994: PetscFunctionBegin;
1995: PetscCall(PetscObjectTypeCompare((PetscObject)A, MATHYPRE, &ishypre));
1996: if (ishypre) {
1997: if (isalpha) {
1998: PetscCall(PetscObjectReference((PetscObject)A));
1999: PetscCall(MatDestroy(&jac->alpha_Poisson));
2000: jac->alpha_Poisson = A;
2001: } else {
2002: if (A) {
2003: PetscCall(PetscObjectReference((PetscObject)A));
2004: } else {
2005: jac->ams_beta_is_zero = PETSC_TRUE;
2006: }
2007: PetscCall(MatDestroy(&jac->beta_Poisson));
2008: jac->beta_Poisson = A;
2009: }
2010: } else {
2011: if (isalpha) {
2012: PetscCall(MatDestroy(&jac->alpha_Poisson));
2013: PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &jac->alpha_Poisson));
2014: } else {
2015: if (A) {
2016: PetscCall(MatDestroy(&jac->beta_Poisson));
2017: PetscCall(MatConvert(A, MATHYPRE, MAT_INITIAL_MATRIX, &jac->beta_Poisson));
2018: } else {
2019: PetscCall(MatDestroy(&jac->beta_Poisson));
2020: jac->ams_beta_is_zero = PETSC_TRUE;
2021: }
2022: }
2023: }
2024: PetscFunctionReturn(PETSC_SUCCESS);
2025: }
2027: /*@
2028: PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix for `PCHYPRE` of type ams
2030: Collective
2032: Input Parameters:
2033: + pc - the preconditioning context
2034: - A - the matrix
2036: Level: intermediate
2038: Note:
2039: A should be obtained by discretizing the vector valued Poisson problem with linear finite elements
2041: Developer Notes:
2042: This automatically converts the matrix to `MATHYPRE` if it is not already of that type
2044: If this is only for `PCHYPRE` type of ams it should be called `PCHYPREAMSSetAlphaPoissonMatrix()`
2046: .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetBetaPoissonMatrix()`
2047: @*/
2048: PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A)
2049: {
2050: PetscFunctionBegin;
2053: PetscCheckSameComm(pc, 1, A, 2);
2054: PetscTryMethod(pc, "PCHYPRESetPoissonMatrix_C", (PC, Mat, PetscBool), (pc, A, PETSC_TRUE));
2055: PetscFunctionReturn(PETSC_SUCCESS);
2056: }
2058: /*@
2059: PCHYPRESetBetaPoissonMatrix - Set Poisson matrix for `PCHYPRE` of type ams
2061: Collective
2063: Input Parameters:
2064: + pc - the preconditioning context
2065: - A - the matrix, or NULL to turn it off
2067: Level: intermediate
2069: Note:
2070: A should be obtained by discretizing the Poisson problem with linear finite elements.
2072: Developer Notes:
2073: This automatically converts the matrix to `MATHYPRE` if it is not already of that type
2075: If this is only for `PCHYPRE` type of ams it should be called `PCHYPREAMSPCHYPRESetBetaPoissonMatrix()`
2077: .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetAlphaPoissonMatrix()`
2078: @*/
2079: PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A)
2080: {
2081: PetscFunctionBegin;
2083: if (A) {
2085: PetscCheckSameComm(pc, 1, A, 2);
2086: }
2087: PetscTryMethod(pc, "PCHYPRESetPoissonMatrix_C", (PC, Mat, PetscBool), (pc, A, PETSC_FALSE));
2088: PetscFunctionReturn(PETSC_SUCCESS);
2089: }
2091: static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE(PC pc, Vec ozz, Vec zoz, Vec zzo)
2092: {
2093: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
2095: PetscFunctionBegin;
2096: /* throw away any vector if already set */
2097: PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[0]));
2098: PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[1]));
2099: PetscCall(VecHYPRE_IJVectorDestroy(&jac->constants[2]));
2100: PetscCall(VecHYPRE_IJVectorCreate(ozz->map, &jac->constants[0]));
2101: PetscCall(VecHYPRE_IJVectorCopy(ozz, jac->constants[0]));
2102: PetscCall(VecHYPRE_IJVectorCreate(zoz->map, &jac->constants[1]));
2103: PetscCall(VecHYPRE_IJVectorCopy(zoz, jac->constants[1]));
2104: jac->dim = 2;
2105: if (zzo) {
2106: PetscCall(VecHYPRE_IJVectorCreate(zzo->map, &jac->constants[2]));
2107: PetscCall(VecHYPRE_IJVectorCopy(zzo, jac->constants[2]));
2108: jac->dim++;
2109: }
2110: PetscFunctionReturn(PETSC_SUCCESS);
2111: }
2113: /*@
2114: PCHYPRESetEdgeConstantVectors - Set the representation of the constant vector fields in the edge element basis for `PCHYPRE` of type ams
2116: Collective
2118: Input Parameters:
2119: + pc - the preconditioning context
2120: . ozz - vector representing (1,0,0) (or (1,0) in 2D)
2121: . zoz - vector representing (0,1,0) (or (0,1) in 2D)
2122: - zzo - vector representing (0,0,1) (use NULL in 2D)
2124: Level: intermediate
2126: Developer Notes:
2127: If this is only for `PCHYPRE` type of ams it should be called `PCHYPREAMSSetEdgeConstantVectors()`
2129: .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetAlphaPoissonMatrix()`
2130: @*/
2131: PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo)
2132: {
2133: PetscFunctionBegin;
2138: PetscCheckSameComm(pc, 1, ozz, 2);
2139: PetscCheckSameComm(pc, 1, zoz, 3);
2140: if (zzo) PetscCheckSameComm(pc, 1, zzo, 4);
2141: PetscTryMethod(pc, "PCHYPRESetEdgeConstantVectors_C", (PC, Vec, Vec, Vec), (pc, ozz, zoz, zzo));
2142: PetscFunctionReturn(PETSC_SUCCESS);
2143: }
2145: static PetscErrorCode PCHYPREAMSSetInteriorNodes_HYPRE(PC pc, Vec interior)
2146: {
2147: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
2149: PetscFunctionBegin;
2150: PetscCall(VecHYPRE_IJVectorDestroy(&jac->interior));
2151: PetscCall(VecHYPRE_IJVectorCreate(interior->map, &jac->interior));
2152: PetscCall(VecHYPRE_IJVectorCopy(interior, jac->interior));
2153: jac->ams_beta_is_zero_part = PETSC_TRUE;
2154: PetscFunctionReturn(PETSC_SUCCESS);
2155: }
2157: /*@
2158: PCHYPREAMSSetInteriorNodes - Set the list of interior nodes to a zero-conductivity region for `PCHYPRE` of type ams
2160: Collective
2162: Input Parameters:
2163: + pc - the preconditioning context
2164: - interior - vector. node is interior if its entry in the array is 1.0.
2166: Level: intermediate
2168: Note:
2169: This calls `HYPRE_AMSSetInteriorNodes()`
2171: .seealso: [](ch_ksp), `PCHYPRE`, `PCHYPRESetDiscreteGradient()`, `PCHYPRESetDiscreteCurl()`, `PCHYPRESetAlphaPoissonMatrix()`
2172: @*/
2173: PetscErrorCode PCHYPREAMSSetInteriorNodes(PC pc, Vec interior)
2174: {
2175: PetscFunctionBegin;
2178: PetscCheckSameComm(pc, 1, interior, 2);
2179: PetscTryMethod(pc, "PCHYPREAMSSetInteriorNodes_C", (PC, Vec), (pc, interior));
2180: PetscFunctionReturn(PETSC_SUCCESS);
2181: }
2183: static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
2184: {
2185: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
2186: Vec tv;
2187: PetscInt i;
2189: PetscFunctionBegin;
2190: /* throw away any coordinate vector if already set */
2191: PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[0]));
2192: PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[1]));
2193: PetscCall(VecHYPRE_IJVectorDestroy(&jac->coords[2]));
2194: jac->dim = dim;
2196: /* compute IJ vector for coordinates */
2197: PetscCall(VecCreate(PetscObjectComm((PetscObject)pc), &tv));
2198: PetscCall(VecSetType(tv, VECSTANDARD));
2199: PetscCall(VecSetSizes(tv, nloc, PETSC_DECIDE));
2200: for (i = 0; i < dim; i++) {
2201: PetscScalar *array;
2202: PetscInt j;
2204: PetscCall(VecHYPRE_IJVectorCreate(tv->map, &jac->coords[i]));
2205: PetscCall(VecGetArrayWrite(tv, &array));
2206: for (j = 0; j < nloc; j++) array[j] = coords[j * dim + i];
2207: PetscCall(VecRestoreArrayWrite(tv, &array));
2208: PetscCall(VecHYPRE_IJVectorCopy(tv, jac->coords[i]));
2209: }
2210: PetscCall(VecDestroy(&tv));
2211: PetscFunctionReturn(PETSC_SUCCESS);
2212: }
2214: static PetscErrorCode PCHYPREGetType_HYPRE(PC pc, const char *name[])
2215: {
2216: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
2218: PetscFunctionBegin;
2219: *name = jac->hypre_type;
2220: PetscFunctionReturn(PETSC_SUCCESS);
2221: }
2223: static PetscErrorCode PCHYPRESetType_HYPRE(PC pc, const char name[])
2224: {
2225: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
2226: PetscBool flag;
2228: PetscFunctionBegin;
2229: if (jac->hypre_type) {
2230: PetscCall(PetscStrcmp(jac->hypre_type, name, &flag));
2231: if (flag) PetscFunctionReturn(PETSC_SUCCESS);
2232: }
2234: PetscCall(PCReset_HYPRE(pc));
2235: PetscCall(PetscFree(jac->hypre_type));
2236: PetscCall(PetscStrallocpy(name, &jac->hypre_type));
2238: jac->maxiter = PETSC_DEFAULT;
2239: jac->tol = PETSC_DEFAULT;
2240: jac->printstatistics = PetscLogPrintInfo;
2242: PetscCall(PetscStrcmp("ilu", jac->hypre_type, &flag));
2243: if (flag) {
2244: PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
2245: PetscCallExternal(HYPRE_ILUCreate, &jac->hsolver);
2246: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ILU;
2247: pc->ops->view = PCView_HYPRE_ILU;
2248: jac->destroy = HYPRE_ILUDestroy;
2249: jac->setup = HYPRE_ILUSetup;
2250: jac->solve = HYPRE_ILUSolve;
2251: jac->factorrowsize = PETSC_DEFAULT;
2252: PetscFunctionReturn(PETSC_SUCCESS);
2253: }
2255: PetscCall(PetscStrcmp("pilut", jac->hypre_type, &flag));
2256: if (flag) {
2257: PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
2258: PetscCallExternal(HYPRE_ParCSRPilutCreate, jac->comm_hypre, &jac->hsolver);
2259: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
2260: pc->ops->view = PCView_HYPRE_Pilut;
2261: jac->destroy = HYPRE_ParCSRPilutDestroy;
2262: jac->setup = HYPRE_ParCSRPilutSetup;
2263: jac->solve = HYPRE_ParCSRPilutSolve;
2264: jac->factorrowsize = PETSC_DEFAULT;
2265: PetscFunctionReturn(PETSC_SUCCESS);
2266: }
2267: PetscCall(PetscStrcmp("euclid", jac->hypre_type, &flag));
2268: if (flag) {
2269: #if defined(PETSC_USE_64BIT_INDICES)
2270: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_SUP, "Hypre Euclid does not support 64-bit indices");
2271: #endif
2272: PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
2273: PetscCallExternal(HYPRE_EuclidCreate, jac->comm_hypre, &jac->hsolver);
2274: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
2275: pc->ops->view = PCView_HYPRE_Euclid;
2276: jac->destroy = HYPRE_EuclidDestroy;
2277: jac->setup = HYPRE_EuclidSetup;
2278: jac->solve = HYPRE_EuclidSolve;
2279: jac->factorrowsize = PETSC_DEFAULT;
2280: jac->eu_level = PETSC_DEFAULT; /* default */
2281: PetscFunctionReturn(PETSC_SUCCESS);
2282: }
2283: PetscCall(PetscStrcmp("parasails", jac->hypre_type, &flag));
2284: if (flag) {
2285: PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &jac->comm_hypre));
2286: PetscCallExternal(HYPRE_ParaSailsCreate, jac->comm_hypre, &jac->hsolver);
2287: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
2288: pc->ops->view = PCView_HYPRE_ParaSails;
2289: jac->destroy = HYPRE_ParaSailsDestroy;
2290: jac->setup = HYPRE_ParaSailsSetup;
2291: jac->solve = HYPRE_ParaSailsSolve;
2292: /* initialize */
2293: jac->nlevels = 1;
2294: jac->threshold = .1;
2295: jac->filter = .1;
2296: jac->loadbal = 0;
2297: if (PetscLogPrintInfo) jac->logging = (int)PETSC_TRUE;
2298: else jac->logging = (int)PETSC_FALSE;
2300: jac->ruse = (int)PETSC_FALSE;
2301: jac->symt = 0;
2302: PetscCallExternal(HYPRE_ParaSailsSetParams, jac->hsolver, jac->threshold, jac->nlevels);
2303: PetscCallExternal(HYPRE_ParaSailsSetFilter, jac->hsolver, jac->filter);
2304: PetscCallExternal(HYPRE_ParaSailsSetLoadbal, jac->hsolver, jac->loadbal);
2305: PetscCallExternal(HYPRE_ParaSailsSetLogging, jac->hsolver, jac->logging);
2306: PetscCallExternal(HYPRE_ParaSailsSetReuse, jac->hsolver, jac->ruse);
2307: PetscCallExternal(HYPRE_ParaSailsSetSym, jac->hsolver, jac->symt);
2308: PetscFunctionReturn(PETSC_SUCCESS);
2309: }
2310: PetscCall(PetscStrcmp("boomeramg", jac->hypre_type, &flag));
2311: if (flag) {
2312: PetscCallExternal(HYPRE_BoomerAMGCreate, &jac->hsolver);
2313: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_BoomerAMG;
2314: pc->ops->view = PCView_HYPRE_BoomerAMG;
2315: pc->ops->applytranspose = PCApplyTranspose_HYPRE_BoomerAMG;
2316: pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
2317: pc->ops->matapply = PCMatApply_HYPRE_BoomerAMG;
2318: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetInterpolations_C", PCGetInterpolations_BoomerAMG));
2319: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCGetCoarseOperators_C", PCGetCoarseOperators_BoomerAMG));
2320: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetCFMarkers_C", PCHYPREGetCFMarkers_BoomerAMG));
2321: jac->destroy = HYPRE_BoomerAMGDestroy;
2322: jac->setup = HYPRE_BoomerAMGSetup;
2323: jac->solve = HYPRE_BoomerAMGSolve;
2324: jac->applyrichardson = PETSC_FALSE;
2325: /* these defaults match the hypre defaults */
2326: jac->cycletype = 1;
2327: jac->maxlevels = 25;
2328: jac->maxiter = 1;
2329: jac->tol = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
2330: jac->truncfactor = 0.0;
2331: jac->strongthreshold = .25;
2332: jac->maxrowsum = .9;
2333: jac->measuretype = 0;
2334: jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
2335: jac->smoothtype = -1; /* Not set by default */
2336: jac->smoothnumlevels = 25;
2337: jac->eu_level = 0;
2338: jac->eu_droptolerance = 0;
2339: jac->eu_bj = 0;
2340: jac->relaxweight = 1.0;
2341: jac->outerrelaxweight = 1.0;
2342: jac->Rtype = 0;
2343: jac->Rstrongthreshold = 0.25;
2344: jac->Rfilterthreshold = 0.0;
2345: jac->Adroptype = -1;
2346: jac->Adroptol = 0.0;
2347: jac->agg_nl = 0;
2348: jac->pmax = 0;
2349: jac->truncfactor = 0.0;
2350: jac->agg_num_paths = 1;
2351: jac->maxc = 9;
2352: jac->minc = 1;
2353: jac->nodal_coarsening = 0;
2354: jac->nodal_coarsening_diag = 0;
2355: jac->vec_interp_variant = 0;
2356: jac->vec_interp_qmax = 0;
2357: jac->vec_interp_smooth = PETSC_FALSE;
2358: jac->interp_refine = 0;
2359: jac->nodal_relax = PETSC_FALSE;
2360: jac->nodal_relax_levels = 1;
2361: jac->rap2 = 0;
2362: PetscObjectParameterSetDefault(jac, relaxorder, -1); /* Initialize with invalid value so we can recognize user input */
2363: PetscFunctionReturn(PETSC_SUCCESS);
2364: }
2365: PetscCall(PetscStrcmp("ams", jac->hypre_type, &flag));
2366: if (flag) {
2367: PetscCallExternal(HYPRE_AMSCreate, &jac->hsolver);
2368: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_AMS;
2369: pc->ops->view = PCView_HYPRE_AMS;
2370: jac->destroy = HYPRE_AMSDestroy;
2371: jac->setup = HYPRE_AMSSetup;
2372: jac->solve = HYPRE_AMSSolve;
2373: jac->coords[0] = NULL;
2374: jac->coords[1] = NULL;
2375: jac->coords[2] = NULL;
2376: jac->interior = NULL;
2377: /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */
2378: jac->as_print = 0;
2379: jac->as_max_iter = 1; /* used as a preconditioner */
2380: jac->as_tol = 0.; /* used as a preconditioner */
2381: jac->ams_cycle_type = 13;
2382: /* Smoothing options */
2383: jac->as_relax_type = 2;
2384: jac->as_relax_times = 1;
2385: jac->as_relax_weight = 1.0;
2386: jac->as_omega = 1.0;
2387: /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
2388: jac->as_amg_alpha_opts[0] = 10;
2389: jac->as_amg_alpha_opts[1] = 1;
2390: jac->as_amg_alpha_opts[2] = 6;
2391: jac->as_amg_alpha_opts[3] = 6;
2392: jac->as_amg_alpha_opts[4] = 4;
2393: jac->as_amg_alpha_theta = 0.25;
2394: /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
2395: jac->as_amg_beta_opts[0] = 10;
2396: jac->as_amg_beta_opts[1] = 1;
2397: jac->as_amg_beta_opts[2] = 6;
2398: jac->as_amg_beta_opts[3] = 6;
2399: jac->as_amg_beta_opts[4] = 4;
2400: jac->as_amg_beta_theta = 0.25;
2401: PetscCallExternal(HYPRE_AMSSetPrintLevel, jac->hsolver, jac->as_print);
2402: PetscCallExternal(HYPRE_AMSSetMaxIter, jac->hsolver, jac->as_max_iter);
2403: PetscCallExternal(HYPRE_AMSSetCycleType, jac->hsolver, jac->ams_cycle_type);
2404: PetscCallExternal(HYPRE_AMSSetTol, jac->hsolver, jac->as_tol);
2405: PetscCallExternal(HYPRE_AMSSetSmoothingOptions, jac->hsolver, jac->as_relax_type, jac->as_relax_times, jac->as_relax_weight, jac->as_omega);
2406: PetscCallExternal(HYPRE_AMSSetAlphaAMGOptions, jac->hsolver, jac->as_amg_alpha_opts[0], /* AMG coarsen type */
2407: jac->as_amg_alpha_opts[1], /* AMG agg_levels */
2408: jac->as_amg_alpha_opts[2], /* AMG relax_type */
2409: jac->as_amg_alpha_theta, jac->as_amg_alpha_opts[3], /* AMG interp_type */
2410: jac->as_amg_alpha_opts[4]); /* AMG Pmax */
2411: PetscCallExternal(HYPRE_AMSSetBetaAMGOptions, jac->hsolver, jac->as_amg_beta_opts[0], /* AMG coarsen type */
2412: jac->as_amg_beta_opts[1], /* AMG agg_levels */
2413: jac->as_amg_beta_opts[2], /* AMG relax_type */
2414: jac->as_amg_beta_theta, jac->as_amg_beta_opts[3], /* AMG interp_type */
2415: jac->as_amg_beta_opts[4]); /* AMG Pmax */
2416: /* Zero conductivity */
2417: jac->ams_beta_is_zero = PETSC_FALSE;
2418: jac->ams_beta_is_zero_part = PETSC_FALSE;
2419: PetscFunctionReturn(PETSC_SUCCESS);
2420: }
2421: PetscCall(PetscStrcmp("ads", jac->hypre_type, &flag));
2422: if (flag) {
2423: PetscCallExternal(HYPRE_ADSCreate, &jac->hsolver);
2424: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ADS;
2425: pc->ops->view = PCView_HYPRE_ADS;
2426: jac->destroy = HYPRE_ADSDestroy;
2427: jac->setup = HYPRE_ADSSetup;
2428: jac->solve = HYPRE_ADSSolve;
2429: jac->coords[0] = NULL;
2430: jac->coords[1] = NULL;
2431: jac->coords[2] = NULL;
2432: /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */
2433: jac->as_print = 0;
2434: jac->as_max_iter = 1; /* used as a preconditioner */
2435: jac->as_tol = 0.; /* used as a preconditioner */
2436: jac->ads_cycle_type = 13;
2437: /* Smoothing options */
2438: jac->as_relax_type = 2;
2439: jac->as_relax_times = 1;
2440: jac->as_relax_weight = 1.0;
2441: jac->as_omega = 1.0;
2442: /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */
2443: jac->ams_cycle_type = 14;
2444: jac->as_amg_alpha_opts[0] = 10;
2445: jac->as_amg_alpha_opts[1] = 1;
2446: jac->as_amg_alpha_opts[2] = 6;
2447: jac->as_amg_alpha_opts[3] = 6;
2448: jac->as_amg_alpha_opts[4] = 4;
2449: jac->as_amg_alpha_theta = 0.25;
2450: /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
2451: jac->as_amg_beta_opts[0] = 10;
2452: jac->as_amg_beta_opts[1] = 1;
2453: jac->as_amg_beta_opts[2] = 6;
2454: jac->as_amg_beta_opts[3] = 6;
2455: jac->as_amg_beta_opts[4] = 4;
2456: jac->as_amg_beta_theta = 0.25;
2457: PetscCallExternal(HYPRE_ADSSetPrintLevel, jac->hsolver, jac->as_print);
2458: PetscCallExternal(HYPRE_ADSSetMaxIter, jac->hsolver, jac->as_max_iter);
2459: PetscCallExternal(HYPRE_ADSSetCycleType, jac->hsolver, jac->ams_cycle_type);
2460: PetscCallExternal(HYPRE_ADSSetTol, jac->hsolver, jac->as_tol);
2461: PetscCallExternal(HYPRE_ADSSetSmoothingOptions, jac->hsolver, jac->as_relax_type, jac->as_relax_times, jac->as_relax_weight, jac->as_omega);
2462: PetscCallExternal(HYPRE_ADSSetAMSOptions, jac->hsolver, jac->ams_cycle_type, /* AMG coarsen type */
2463: jac->as_amg_alpha_opts[0], /* AMG coarsen type */
2464: jac->as_amg_alpha_opts[1], /* AMG agg_levels */
2465: jac->as_amg_alpha_opts[2], /* AMG relax_type */
2466: jac->as_amg_alpha_theta, jac->as_amg_alpha_opts[3], /* AMG interp_type */
2467: jac->as_amg_alpha_opts[4]); /* AMG Pmax */
2468: PetscCallExternal(HYPRE_ADSSetAMGOptions, jac->hsolver, jac->as_amg_beta_opts[0], /* AMG coarsen type */
2469: jac->as_amg_beta_opts[1], /* AMG agg_levels */
2470: jac->as_amg_beta_opts[2], /* AMG relax_type */
2471: jac->as_amg_beta_theta, jac->as_amg_beta_opts[3], /* AMG interp_type */
2472: jac->as_amg_beta_opts[4]); /* AMG Pmax */
2473: PetscFunctionReturn(PETSC_SUCCESS);
2474: }
2475: PetscCall(PetscFree(jac->hypre_type));
2477: jac->hypre_type = NULL;
2478: SETERRQ(PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_UNKNOWN_TYPE, "Unknown HYPRE preconditioner %s; Choices are euclid, ilu, pilut, parasails, boomeramg, ams, ads", name);
2479: }
2481: /*
2482: It only gets here if the HYPRE type has not been set before the call to
2483: ...SetFromOptions() which actually is most of the time
2484: */
2485: static PetscErrorCode PCSetFromOptions_HYPRE(PC pc, PetscOptionItems PetscOptionsObject)
2486: {
2487: PetscInt indx;
2488: const char *type[] = {"ilu", "euclid", "pilut", "parasails", "boomeramg", "ams", "ads"};
2489: PetscBool flg;
2490: PC_HYPRE *jac = (PC_HYPRE *)pc->data;
2492: PetscFunctionBegin;
2493: PetscOptionsHeadBegin(PetscOptionsObject, "HYPRE preconditioner options");
2494: PetscCall(PetscOptionsEList("-pc_hypre_type", "HYPRE preconditioner type", "PCHYPRESetType", type, PETSC_STATIC_ARRAY_LENGTH(type), "boomeramg", &indx, &flg));
2495: if (flg) PetscCall(PCHYPRESetType_HYPRE(pc, type[indx]));
2496: /*
2497: Set the type if it was never set.
2498: */
2499: if (!jac->hypre_type) PetscCall(PCHYPRESetType_HYPRE(pc, "boomeramg"));
2500: PetscTryTypeMethod(pc, setfromoptions, PetscOptionsObject);
2501: PetscOptionsHeadEnd();
2502: PetscFunctionReturn(PETSC_SUCCESS);
2503: }
2505: /*@
2506: PCHYPRESetType - Sets which hypre preconditioner you wish to use
2508: Input Parameters:
2509: + pc - the preconditioner context
2510: - name - either euclid, ilu, pilut, parasails, boomeramg, ams, ads
2512: Options Database Key:
2513: . pc_hypre_type - One of euclid, ilu, pilut, parasails, boomeramg, ams, ads
2515: Level: intermediate
2517: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHYPRE`
2518: @*/
2519: PetscErrorCode PCHYPRESetType(PC pc, const char name[])
2520: {
2521: PetscFunctionBegin;
2523: PetscAssertPointer(name, 2);
2524: PetscTryMethod(pc, "PCHYPRESetType_C", (PC, const char[]), (pc, name));
2525: PetscFunctionReturn(PETSC_SUCCESS);
2526: }
2528: /*@C
2529: PCHYPREGetCFMarkers - Gets CF marker arrays for all levels (except the finest level)
2531: Logically Collective
2533: Input Parameter:
2534: . pc - the preconditioner context
2536: Output Parameters:
2537: + n_per_level - the number of nodes per level (size of `num_levels`)
2538: - CFMarkers - the Coarse/Fine Boolean arrays (size of `num_levels` - 1)
2540: Note:
2541: Caller is responsible for memory management of `n_per_level` and `CFMarkers` pointers. That is they should free them with `PetscFree()` when no longer needed.
2543: Level: advanced
2545: .seealso: [](ch_ksp), `PC`, `PCMG`, `PCMGGetRestriction()`, `PCMGSetInterpolation()`, `PCMGGetRScale()`, `PCMGGetInterpolation()`, `PCGetInterpolations()`
2546: @*/
2547: PetscErrorCode PCHYPREGetCFMarkers(PC pc, PetscInt *n_per_level[], PetscBT *CFMarkers[])
2548: {
2549: PetscFunctionBegin;
2551: PetscAssertPointer(n_per_level, 2);
2552: PetscAssertPointer(CFMarkers, 3);
2553: PetscUseMethod(pc, "PCHYPREGetCFMarkers_C", (PC, PetscInt *[], PetscBT *[]), (pc, n_per_level, CFMarkers));
2554: PetscFunctionReturn(PETSC_SUCCESS);
2555: }
2557: /*@
2558: PCHYPREGetType - Gets which hypre preconditioner you are using
2560: Input Parameter:
2561: . pc - the preconditioner context
2563: Output Parameter:
2564: . name - either euclid, ilu, pilut, parasails, boomeramg, ams, ads
2566: Level: intermediate
2568: .seealso: [](ch_ksp), `PCCreate()`, `PCHYPRESetType()`, `PCType`, `PC`, `PCHYPRE`
2569: @*/
2570: PetscErrorCode PCHYPREGetType(PC pc, const char *name[])
2571: {
2572: PetscFunctionBegin;
2574: PetscAssertPointer(name, 2);
2575: PetscTryMethod(pc, "PCHYPREGetType_C", (PC, const char *[]), (pc, name));
2576: PetscFunctionReturn(PETSC_SUCCESS);
2577: }
2579: /*@
2580: PCMGGalerkinSetMatProductAlgorithm - Set type of SpGEMM for hypre to use on GPUs
2582: Logically Collective
2584: Input Parameters:
2585: + pc - the hypre context
2586: - name - one of 'cusparse', 'hypre'
2588: Options Database Key:
2589: . -pc_mg_galerkin_mat_product_algorithm <cusparse,hypre> - Type of SpGEMM to use in hypre
2591: Level: intermediate
2593: Developer Notes:
2594: How the name starts with `PCMG`, should it not be `PCHYPREBoomerAMG`?
2596: .seealso: [](ch_ksp), `PCHYPRE`, `PCMGGalerkinGetMatProductAlgorithm()`
2597: @*/
2598: PetscErrorCode PCMGGalerkinSetMatProductAlgorithm(PC pc, const char name[])
2599: {
2600: PetscFunctionBegin;
2602: PetscTryMethod(pc, "PCMGGalerkinSetMatProductAlgorithm_C", (PC, const char[]), (pc, name));
2603: PetscFunctionReturn(PETSC_SUCCESS);
2604: }
2606: /*@
2607: PCMGGalerkinGetMatProductAlgorithm - Get type of SpGEMM for hypre to use on GPUs
2609: Not Collective
2611: Input Parameter:
2612: . pc - the multigrid context
2614: Output Parameter:
2615: . name - one of 'cusparse', 'hypre'
2617: Level: intermediate
2619: .seealso: [](ch_ksp), `PCHYPRE`, `PCMGGalerkinSetMatProductAlgorithm()`
2620: @*/
2621: PetscErrorCode PCMGGalerkinGetMatProductAlgorithm(PC pc, const char *name[])
2622: {
2623: PetscFunctionBegin;
2625: PetscTryMethod(pc, "PCMGGalerkinGetMatProductAlgorithm_C", (PC, const char *[]), (pc, name));
2626: PetscFunctionReturn(PETSC_SUCCESS);
2627: }
2629: /*MC
2630: PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre as PETSc `PC`
2632: Options Database Keys:
2633: + -pc_hypre_type - One of `euclid`, `pilut`, `parasails`, `boomeramg`, `ams`, or `ads`
2634: . -pc_hypre_boomeramg_nodal_coarsen <n> - where n is from 1 to 6 (see `HYPRE_BoomerAMGSetNodal()`)
2635: . -pc_hypre_boomeramg_vec_interp_variant <v> - where v is from 1 to 3 (see `HYPRE_BoomerAMGSetInterpVecVariant()`)
2636: - Many others - run with `-pc_type hypre` `-pc_hypre_type XXX` `-help` to see options for the XXX preconditioner
2638: Level: intermediate
2640: Notes:
2641: Apart from `-pc_hypre_type` (for which there is `PCHYPRESetType()`),
2642: the many hypre options can ONLY be set via the options database (e.g. the command line
2643: or with `PetscOptionsSetValue()`, there are no functions to set them)
2645: The options `-pc_hypre_boomeramg_max_iter` and `-pc_hypre_boomeramg_tol` refer to the number of iterations
2646: (V-cycles) and tolerance that boomerAMG does EACH time it is called. So for example, if
2647: `-pc_hypre_boomeramg_max_iter` is set to 2 then 2-V-cycles are being used to define the preconditioner
2648: (`-pc_hypre_boomeramg_tol` should be set to 0.0 - the default - to strictly use a fixed number of
2649: iterations per hypre call). `-ksp_max_it` and `-ksp_rtol` STILL determine the total number of iterations
2650: and tolerance for the Krylov solver. For example, if `-pc_hypre_boomeramg_max_iter` is 2 and `-ksp_max_it` is 10
2651: then AT MOST twenty V-cycles of boomeramg will be used.
2653: Note that the option `-pc_hypre_boomeramg_relax_type_all` defaults to symmetric relaxation
2654: (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry.
2655: Otherwise, you may want to use `-pc_hypre_boomeramg_relax_type_all SOR/Jacobi`.
2657: `MatSetNearNullSpace()` - if you provide a near null space to your matrix it is ignored by hypre UNLESS you also use
2658: the following two options: `-pc_hypre_boomeramg_nodal_coarsen <n> -pc_hypre_boomeramg_vec_interp_variant <v>`
2660: See `PCPFMG`, `PCSMG`, and `PCSYSPFMG` for access to hypre's other (nonalgebraic) multigrid solvers
2662: For `PCHYPRE` type of `ams` or `ads` auxiliary data must be provided to the preconditioner with `PCHYPRESetDiscreteGradient()`,
2663: `PCHYPRESetDiscreteCurl()`, `PCHYPRESetInterpolations()`, `PCHYPRESetAlphaPoissonMatrix()`, `PCHYPRESetBetaPoissonMatrix()`, `PCHYPRESetEdgeConstantVectors()`,
2664: `PCHYPREAMSSetInteriorNodes()`
2666: 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
2667: since it is usually slower, one should use a `KSPType` of `KSPRICHARDSON`
2668: (or equivalently `-ksp_type richardson`) to achieve this. Using `KSPPREONLY` will not work since it only applies a single cycle of multigrid.
2670: PETSc provides its own geometric and algebraic multigrid solvers `PCMG` and `PCGAMG`, also see `PCHMG` which is useful for certain multicomponent problems
2672: GPU Notes:
2673: To configure hypre BoomerAMG so that it can utilize NVIDIA GPUs run ./configure --download-hypre --with-cuda
2674: Then pass `VECCUDA` vectors and `MATAIJCUSPARSE` matrices to the solvers and PETSc will automatically utilize hypre's GPU solvers.
2676: To configure hypre BoomerAMG so that it can utilize AMD GPUs run ./configure --download-hypre --with-hip
2677: Then pass `VECHIP` vectors to the solvers and PETSc will automatically utilize hypre's GPU solvers.
2679: .seealso: [](ch_ksp), `PCCreate()`, `PCSetType()`, `PCType`, `PC`, `PCHYPRESetType()`, `PCPFMG`, `PCGAMG`, `PCSYSPFMG`, `PCSMG`, `PCHYPRESetDiscreteGradient()`,
2680: `PCHYPRESetDiscreteCurl()`, `PCHYPRESetInterpolations()`, `PCHYPRESetAlphaPoissonMatrix()`, `PCHYPRESetBetaPoissonMatrix()`, `PCHYPRESetEdgeConstantVectors()`,
2681: PCHYPREAMSSetInteriorNodes()
2682: M*/
2684: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
2685: {
2686: PC_HYPRE *jac;
2688: PetscFunctionBegin;
2689: PetscCall(PetscNew(&jac));
2691: pc->data = jac;
2692: pc->ops->reset = PCReset_HYPRE;
2693: pc->ops->destroy = PCDestroy_HYPRE;
2694: pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
2695: pc->ops->setup = PCSetUp_HYPRE;
2696: pc->ops->apply = PCApply_HYPRE;
2697: jac->hypre_type = NULL;
2698: jac->comm_hypre = MPI_COMM_NULL;
2699: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetType_C", PCHYPRESetType_HYPRE));
2700: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREGetType_C", PCHYPREGetType_HYPRE));
2701: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCSetCoordinates_C", PCSetCoordinates_HYPRE));
2702: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteGradient_C", PCHYPRESetDiscreteGradient_HYPRE));
2703: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetDiscreteCurl_C", PCHYPRESetDiscreteCurl_HYPRE));
2704: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetInterpolations_C", PCHYPRESetInterpolations_HYPRE));
2705: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetEdgeConstantVectors_C", PCHYPRESetEdgeConstantVectors_HYPRE));
2706: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPREAMSSetInteriorNodes_C", PCHYPREAMSSetInteriorNodes_HYPRE));
2707: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCHYPRESetPoissonMatrix_C", PCHYPRESetPoissonMatrix_HYPRE));
2708: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinSetMatProductAlgorithm_C", PCMGGalerkinSetMatProductAlgorithm_HYPRE_BoomerAMG));
2709: PetscCall(PetscObjectComposeFunction((PetscObject)pc, "PCMGGalerkinGetMatProductAlgorithm_C", PCMGGalerkinGetMatProductAlgorithm_HYPRE_BoomerAMG));
2710: #if defined(PETSC_HAVE_HYPRE_DEVICE)
2711: #if defined(HYPRE_USING_HIP)
2712: PetscCall(PetscDeviceInitialize(PETSC_DEVICE_HIP));
2713: #endif
2714: #if defined(HYPRE_USING_CUDA)
2715: PetscCall(PetscDeviceInitialize(PETSC_DEVICE_CUDA));
2716: #endif
2717: #endif
2718: PetscHYPREInitialize();
2719: PetscFunctionReturn(PETSC_SUCCESS);
2720: }
2722: typedef struct {
2723: MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */
2724: HYPRE_StructSolver hsolver;
2726: /* keep copy of PFMG options used so may view them */
2727: PetscInt its;
2728: PetscReal tol;
2729: PetscInt relax_type;
2730: PetscInt rap_type;
2731: PetscInt num_pre_relax, num_post_relax;
2732: PetscInt max_levels;
2733: PetscInt skip_relax;
2734: PetscBool print_statistics;
2735: } PC_PFMG;
2737: static PetscErrorCode PCDestroy_PFMG(PC pc)
2738: {
2739: PC_PFMG *ex = (PC_PFMG *)pc->data;
2741: PetscFunctionBegin;
2742: if (ex->hsolver) PetscCallExternal(HYPRE_StructPFMGDestroy, ex->hsolver);
2743: PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2744: PetscCall(PetscFree(pc->data));
2745: PetscFunctionReturn(PETSC_SUCCESS);
2746: }
2748: static const char *PFMGRelaxType[] = {"Jacobi", "Weighted-Jacobi", "symmetric-Red/Black-Gauss-Seidel", "Red/Black-Gauss-Seidel"};
2749: static const char *PFMGRAPType[] = {"Galerkin", "non-Galerkin"};
2751: static PetscErrorCode PCView_PFMG(PC pc, PetscViewer viewer)
2752: {
2753: PetscBool isascii;
2754: PC_PFMG *ex = (PC_PFMG *)pc->data;
2756: PetscFunctionBegin;
2757: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2758: if (isascii) {
2759: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE PFMG preconditioning\n"));
2760: PetscCall(PetscViewerASCIIPrintf(viewer, " max iterations %" PetscInt_FMT "\n", ex->its));
2761: PetscCall(PetscViewerASCIIPrintf(viewer, " tolerance %g\n", ex->tol));
2762: PetscCall(PetscViewerASCIIPrintf(viewer, " relax type %s\n", PFMGRelaxType[ex->relax_type]));
2763: PetscCall(PetscViewerASCIIPrintf(viewer, " RAP type %s\n", PFMGRAPType[ex->rap_type]));
2764: PetscCall(PetscViewerASCIIPrintf(viewer, " number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax));
2765: PetscCall(PetscViewerASCIIPrintf(viewer, " max levels %" PetscInt_FMT "\n", ex->max_levels));
2766: PetscCall(PetscViewerASCIIPrintf(viewer, " skip relax %" PetscInt_FMT "\n", ex->skip_relax));
2767: }
2768: PetscFunctionReturn(PETSC_SUCCESS);
2769: }
2771: static PetscErrorCode PCSetFromOptions_PFMG(PC pc, PetscOptionItems PetscOptionsObject)
2772: {
2773: PC_PFMG *ex = (PC_PFMG *)pc->data;
2775: PetscFunctionBegin;
2776: PetscOptionsHeadBegin(PetscOptionsObject, "PFMG options");
2777: PetscCall(PetscOptionsBool("-pc_pfmg_print_statistics", "Print statistics", "HYPRE_StructPFMGSetPrintLevel", ex->print_statistics, &ex->print_statistics, NULL));
2778: PetscCall(PetscOptionsInt("-pc_pfmg_its", "Number of iterations of PFMG to use as preconditioner", "HYPRE_StructPFMGSetMaxIter", ex->its, &ex->its, NULL));
2779: PetscCallExternal(HYPRE_StructPFMGSetMaxIter, ex->hsolver, ex->its);
2780: 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));
2781: PetscCallExternal(HYPRE_StructPFMGSetNumPreRelax, ex->hsolver, ex->num_pre_relax);
2782: 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));
2783: PetscCallExternal(HYPRE_StructPFMGSetNumPostRelax, ex->hsolver, ex->num_post_relax);
2785: PetscCall(PetscOptionsInt("-pc_pfmg_max_levels", "Max Levels for MG hierarchy", "HYPRE_StructPFMGSetMaxLevels", ex->max_levels, &ex->max_levels, NULL));
2786: PetscCallExternal(HYPRE_StructPFMGSetMaxLevels, ex->hsolver, ex->max_levels);
2788: PetscCall(PetscOptionsReal("-pc_pfmg_tol", "Tolerance of PFMG", "HYPRE_StructPFMGSetTol", ex->tol, &ex->tol, NULL));
2789: PetscCallExternal(HYPRE_StructPFMGSetTol, ex->hsolver, ex->tol);
2790: 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));
2791: PetscCallExternal(HYPRE_StructPFMGSetRelaxType, ex->hsolver, ex->relax_type);
2792: PetscCall(PetscOptionsEList("-pc_pfmg_rap_type", "RAP type", "HYPRE_StructPFMGSetRAPType", PFMGRAPType, PETSC_STATIC_ARRAY_LENGTH(PFMGRAPType), PFMGRAPType[ex->rap_type], &ex->rap_type, NULL));
2793: PetscCallExternal(HYPRE_StructPFMGSetRAPType, ex->hsolver, ex->rap_type);
2794: 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));
2795: PetscCallExternal(HYPRE_StructPFMGSetSkipRelax, ex->hsolver, ex->skip_relax);
2796: PetscOptionsHeadEnd();
2797: PetscFunctionReturn(PETSC_SUCCESS);
2798: }
2800: static PetscErrorCode PCApply_PFMG(PC pc, Vec x, Vec y)
2801: {
2802: PC_PFMG *ex = (PC_PFMG *)pc->data;
2803: PetscScalar *yy;
2804: const PetscScalar *xx;
2805: PetscInt ilower[3], iupper[3];
2806: HYPRE_Int hlower[3], hupper[3];
2807: Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data;
2809: PetscFunctionBegin;
2810: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
2811: PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
2812: /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2813: iupper[0] += ilower[0] - 1;
2814: iupper[1] += ilower[1] - 1;
2815: iupper[2] += ilower[2] - 1;
2816: hlower[0] = (HYPRE_Int)ilower[0];
2817: hlower[1] = (HYPRE_Int)ilower[1];
2818: hlower[2] = (HYPRE_Int)ilower[2];
2819: hupper[0] = (HYPRE_Int)iupper[0];
2820: hupper[1] = (HYPRE_Int)iupper[1];
2821: hupper[2] = (HYPRE_Int)iupper[2];
2823: /* copy x values over to hypre */
2824: PetscCallExternal(HYPRE_StructVectorSetConstantValues, mx->hb, 0.0);
2825: PetscCall(VecGetArrayRead(x, &xx));
2826: PetscCallExternal(HYPRE_StructVectorSetBoxValues, mx->hb, hlower, hupper, (HYPRE_Complex *)xx);
2827: PetscCall(VecRestoreArrayRead(x, &xx));
2828: PetscCallExternal(HYPRE_StructVectorAssemble, mx->hb);
2829: PetscCallExternal(HYPRE_StructPFMGSolve, ex->hsolver, mx->hmat, mx->hb, mx->hx);
2831: /* copy solution values back to PETSc */
2832: PetscCall(VecGetArray(y, &yy));
2833: PetscCallExternal(HYPRE_StructVectorGetBoxValues, mx->hx, hlower, hupper, (HYPRE_Complex *)yy);
2834: PetscCall(VecRestoreArray(y, &yy));
2835: PetscFunctionReturn(PETSC_SUCCESS);
2836: }
2838: 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)
2839: {
2840: PC_PFMG *jac = (PC_PFMG *)pc->data;
2841: HYPRE_Int oits;
2843: PetscFunctionBegin;
2844: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
2845: PetscCallExternal(HYPRE_StructPFMGSetMaxIter, jac->hsolver, its * jac->its);
2846: PetscCallExternal(HYPRE_StructPFMGSetTol, jac->hsolver, rtol);
2848: PetscCall(PCApply_PFMG(pc, b, y));
2849: PetscCallExternal(HYPRE_StructPFMGGetNumIterations, jac->hsolver, &oits);
2850: *outits = oits;
2851: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2852: else *reason = PCRICHARDSON_CONVERGED_RTOL;
2853: PetscCallExternal(HYPRE_StructPFMGSetTol, jac->hsolver, jac->tol);
2854: PetscCallExternal(HYPRE_StructPFMGSetMaxIter, jac->hsolver, jac->its);
2855: PetscFunctionReturn(PETSC_SUCCESS);
2856: }
2858: static PetscErrorCode PCSetUp_PFMG(PC pc)
2859: {
2860: PC_PFMG *ex = (PC_PFMG *)pc->data;
2861: Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data;
2862: PetscBool flg;
2864: PetscFunctionBegin;
2865: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRESTRUCT, &flg));
2866: PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use MATHYPRESTRUCT with this preconditioner");
2868: /* create the hypre solver object and set its information */
2869: if (ex->hsolver) PetscCallExternal(HYPRE_StructPFMGDestroy, ex->hsolver);
2870: PetscCallExternal(HYPRE_StructPFMGCreate, ex->hcomm, &ex->hsolver);
2872: // Print Hypre statistics about the solve process
2873: if (ex->print_statistics) PetscCallExternal(HYPRE_StructPFMGSetPrintLevel, ex->hsolver, 3);
2875: // The hypre options must be repeated here because the StructPFMG was destroyed and recreated
2876: PetscCallExternal(HYPRE_StructPFMGSetMaxIter, ex->hsolver, ex->its);
2877: PetscCallExternal(HYPRE_StructPFMGSetNumPreRelax, ex->hsolver, ex->num_pre_relax);
2878: PetscCallExternal(HYPRE_StructPFMGSetNumPostRelax, ex->hsolver, ex->num_post_relax);
2879: PetscCallExternal(HYPRE_StructPFMGSetMaxLevels, ex->hsolver, ex->max_levels);
2880: PetscCallExternal(HYPRE_StructPFMGSetTol, ex->hsolver, ex->tol);
2881: PetscCallExternal(HYPRE_StructPFMGSetRelaxType, ex->hsolver, ex->relax_type);
2882: PetscCallExternal(HYPRE_StructPFMGSetRAPType, ex->hsolver, ex->rap_type);
2884: PetscCallExternal(HYPRE_StructPFMGSetup, ex->hsolver, mx->hmat, mx->hb, mx->hx);
2885: PetscCallExternal(HYPRE_StructPFMGSetZeroGuess, ex->hsolver);
2886: PetscFunctionReturn(PETSC_SUCCESS);
2887: }
2889: /*MC
2890: PCPFMG - the hypre PFMG multigrid solver
2892: Options Database Keys:
2893: + -pc_pfmg_its <its> - number of iterations of PFMG to use as preconditioner
2894: . -pc_pfmg_num_pre_relax <steps> - number of smoothing steps before coarse grid solve
2895: . -pc_pfmg_num_post_relax <steps> - number of smoothing steps after coarse grid solve
2896: . -pc_pfmg_tol <tol> - tolerance of PFMG
2897: . -pc_pfmg_relax_type - relaxation type for the up and down cycles, one of Jacobi,Weighted-Jacobi,symmetric-Red/Black-Gauss-Seidel,Red/Black-Gauss-Seidel
2898: . -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin
2899: - -pc_pfmg_skip_relax - skip relaxation on certain grids for isotropic problems. This can greatly improve efficiency by eliminating unnecessary relaxations
2900: when the underlying problem is isotropic, one of 0,1
2902: Level: advanced
2904: Notes:
2905: This is for CELL-centered descretizations
2907: See `PCSYSPFMG` for a version suitable for systems of PDEs, and `PCSMG`
2909: See `PCHYPRE` for hypre's BoomerAMG algebraic multigrid solver
2911: This must be used with the `MATHYPRESTRUCT` matrix type.
2913: This provides only some of the functionality of PFMG, it supports only one block per process defined by a PETSc `DMDA`.
2915: .seealso: [](ch_ksp), `PCMG`, `MATHYPRESTRUCT`, `PCHYPRE`, `PCGAMG`, `PCSYSPFMG`, `PCSMG`
2916: M*/
2918: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
2919: {
2920: PC_PFMG *ex;
2922: PetscFunctionBegin;
2923: PetscCall(PetscNew(&ex));
2924: pc->data = ex;
2926: ex->its = 1;
2927: ex->tol = 1.e-8;
2928: ex->relax_type = 1;
2929: ex->rap_type = 0;
2930: ex->num_pre_relax = 1;
2931: ex->num_post_relax = 1;
2932: ex->max_levels = 0;
2933: ex->skip_relax = 0;
2934: ex->print_statistics = PETSC_FALSE;
2936: pc->ops->setfromoptions = PCSetFromOptions_PFMG;
2937: pc->ops->view = PCView_PFMG;
2938: pc->ops->destroy = PCDestroy_PFMG;
2939: pc->ops->apply = PCApply_PFMG;
2940: pc->ops->applyrichardson = PCApplyRichardson_PFMG;
2941: pc->ops->setup = PCSetUp_PFMG;
2943: PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2944: PetscHYPREInitialize();
2945: PetscCallExternal(HYPRE_StructPFMGCreate, ex->hcomm, &ex->hsolver);
2946: PetscFunctionReturn(PETSC_SUCCESS);
2947: }
2949: /* we know we are working with a HYPRE_SStructMatrix */
2950: typedef struct {
2951: MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
2952: HYPRE_SStructSolver ss_solver;
2954: /* keep copy of SYSPFMG options used so may view them */
2955: PetscInt its;
2956: PetscReal tol;
2957: PetscInt relax_type;
2958: PetscInt num_pre_relax, num_post_relax;
2959: } PC_SysPFMG;
2961: static PetscErrorCode PCDestroy_SysPFMG(PC pc)
2962: {
2963: PC_SysPFMG *ex = (PC_SysPFMG *)pc->data;
2965: PetscFunctionBegin;
2966: if (ex->ss_solver) PetscCallExternal(HYPRE_SStructSysPFMGDestroy, ex->ss_solver);
2967: PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
2968: PetscCall(PetscFree(pc->data));
2969: PetscFunctionReturn(PETSC_SUCCESS);
2970: }
2972: static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi", "Red/Black-Gauss-Seidel"};
2974: static PetscErrorCode PCView_SysPFMG(PC pc, PetscViewer viewer)
2975: {
2976: PetscBool isascii;
2977: PC_SysPFMG *ex = (PC_SysPFMG *)pc->data;
2979: PetscFunctionBegin;
2980: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
2981: if (isascii) {
2982: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE SysPFMG preconditioning\n"));
2983: PetscCall(PetscViewerASCIIPrintf(viewer, " max iterations %" PetscInt_FMT "\n", ex->its));
2984: PetscCall(PetscViewerASCIIPrintf(viewer, " tolerance %g\n", ex->tol));
2985: PetscCall(PetscViewerASCIIPrintf(viewer, " relax type %s\n", PFMGRelaxType[ex->relax_type]));
2986: PetscCall(PetscViewerASCIIPrintf(viewer, " number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax));
2987: }
2988: PetscFunctionReturn(PETSC_SUCCESS);
2989: }
2991: static PetscErrorCode PCSetFromOptions_SysPFMG(PC pc, PetscOptionItems PetscOptionsObject)
2992: {
2993: PC_SysPFMG *ex = (PC_SysPFMG *)pc->data;
2994: PetscBool flg = PETSC_FALSE;
2996: PetscFunctionBegin;
2997: PetscOptionsHeadBegin(PetscOptionsObject, "SysPFMG options");
2998: PetscCall(PetscOptionsBool("-pc_syspfmg_print_statistics", "Print statistics", "HYPRE_SStructSysPFMGSetPrintLevel", flg, &flg, NULL));
2999: if (flg) PetscCallExternal(HYPRE_SStructSysPFMGSetPrintLevel, ex->ss_solver, 3);
3000: PetscCall(PetscOptionsInt("-pc_syspfmg_its", "Number of iterations of SysPFMG to use as preconditioner", "HYPRE_SStructSysPFMGSetMaxIter", ex->its, &ex->its, NULL));
3001: PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter, ex->ss_solver, ex->its);
3002: 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));
3003: PetscCallExternal(HYPRE_SStructSysPFMGSetNumPreRelax, ex->ss_solver, ex->num_pre_relax);
3004: 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));
3005: PetscCallExternal(HYPRE_SStructSysPFMGSetNumPostRelax, ex->ss_solver, ex->num_post_relax);
3007: PetscCall(PetscOptionsReal("-pc_syspfmg_tol", "Tolerance of SysPFMG", "HYPRE_SStructSysPFMGSetTol", ex->tol, &ex->tol, NULL));
3008: PetscCallExternal(HYPRE_SStructSysPFMGSetTol, ex->ss_solver, ex->tol);
3009: 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));
3010: PetscCallExternal(HYPRE_SStructSysPFMGSetRelaxType, ex->ss_solver, ex->relax_type);
3011: PetscOptionsHeadEnd();
3012: PetscFunctionReturn(PETSC_SUCCESS);
3013: }
3015: static PetscErrorCode PCApply_SysPFMG(PC pc, Vec x, Vec y)
3016: {
3017: PC_SysPFMG *ex = (PC_SysPFMG *)pc->data;
3018: PetscScalar *yy;
3019: const PetscScalar *xx;
3020: PetscInt ilower[3], iupper[3];
3021: HYPRE_Int hlower[3], hupper[3];
3022: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)pc->pmat->data;
3023: PetscInt ordering = mx->dofs_order;
3024: PetscInt nvars = mx->nvars;
3025: PetscInt part = 0;
3026: PetscInt size;
3027: PetscInt i;
3029: PetscFunctionBegin;
3030: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
3031: PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
3032: /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
3033: iupper[0] += ilower[0] - 1;
3034: iupper[1] += ilower[1] - 1;
3035: iupper[2] += ilower[2] - 1;
3036: hlower[0] = (HYPRE_Int)ilower[0];
3037: hlower[1] = (HYPRE_Int)ilower[1];
3038: hlower[2] = (HYPRE_Int)ilower[2];
3039: hupper[0] = (HYPRE_Int)iupper[0];
3040: hupper[1] = (HYPRE_Int)iupper[1];
3041: hupper[2] = (HYPRE_Int)iupper[2];
3043: size = 1;
3044: for (i = 0; i < 3; i++) size *= (iupper[i] - ilower[i] + 1);
3046: /* copy x values over to hypre for variable ordering */
3047: if (ordering) {
3048: PetscCallExternal(HYPRE_SStructVectorSetConstantValues, mx->ss_b, 0.0);
3049: PetscCall(VecGetArrayRead(x, &xx));
3050: for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorSetBoxValues, mx->ss_b, part, hlower, hupper, i, (HYPRE_Complex *)(xx + (size * i)));
3051: PetscCall(VecRestoreArrayRead(x, &xx));
3052: PetscCallExternal(HYPRE_SStructVectorAssemble, mx->ss_b);
3053: PetscCallExternal(HYPRE_SStructMatrixMatvec, 1.0, mx->ss_mat, mx->ss_b, 0.0, mx->ss_x);
3054: PetscCallExternal(HYPRE_SStructSysPFMGSolve, ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x);
3056: /* copy solution values back to PETSc */
3057: PetscCall(VecGetArray(y, &yy));
3058: for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues, mx->ss_x, part, hlower, hupper, i, (HYPRE_Complex *)(yy + (size * i)));
3059: PetscCall(VecRestoreArray(y, &yy));
3060: } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
3061: PetscScalar *z;
3062: PetscInt j, k;
3064: PetscCall(PetscMalloc1(nvars * size, &z));
3065: PetscCallExternal(HYPRE_SStructVectorSetConstantValues, mx->ss_b, 0.0);
3066: PetscCall(VecGetArrayRead(x, &xx));
3068: /* transform nodal to hypre's variable ordering for sys_pfmg */
3069: for (i = 0; i < size; i++) {
3070: k = i * nvars;
3071: for (j = 0; j < nvars; j++) z[j * size + i] = xx[k + j];
3072: }
3073: for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorSetBoxValues, mx->ss_b, part, hlower, hupper, i, (HYPRE_Complex *)(z + (size * i)));
3074: PetscCall(VecRestoreArrayRead(x, &xx));
3075: PetscCallExternal(HYPRE_SStructVectorAssemble, mx->ss_b);
3076: PetscCallExternal(HYPRE_SStructSysPFMGSolve, ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x);
3078: /* copy solution values back to PETSc */
3079: PetscCall(VecGetArray(y, &yy));
3080: for (i = 0; i < nvars; i++) PetscCallExternal(HYPRE_SStructVectorGetBoxValues, mx->ss_x, part, hlower, hupper, i, (HYPRE_Complex *)(z + (size * i)));
3081: /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
3082: for (i = 0; i < size; i++) {
3083: k = i * nvars;
3084: for (j = 0; j < nvars; j++) yy[k + j] = z[j * size + i];
3085: }
3086: PetscCall(VecRestoreArray(y, &yy));
3087: PetscCall(PetscFree(z));
3088: }
3089: PetscFunctionReturn(PETSC_SUCCESS);
3090: }
3092: 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)
3093: {
3094: PC_SysPFMG *jac = (PC_SysPFMG *)pc->data;
3095: HYPRE_Int oits;
3097: PetscFunctionBegin;
3098: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
3099: PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter, jac->ss_solver, its * jac->its);
3100: PetscCallExternal(HYPRE_SStructSysPFMGSetTol, jac->ss_solver, rtol);
3101: PetscCall(PCApply_SysPFMG(pc, b, y));
3102: PetscCallExternal(HYPRE_SStructSysPFMGGetNumIterations, jac->ss_solver, &oits);
3103: *outits = oits;
3104: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
3105: else *reason = PCRICHARDSON_CONVERGED_RTOL;
3106: PetscCallExternal(HYPRE_SStructSysPFMGSetTol, jac->ss_solver, jac->tol);
3107: PetscCallExternal(HYPRE_SStructSysPFMGSetMaxIter, jac->ss_solver, jac->its);
3108: PetscFunctionReturn(PETSC_SUCCESS);
3109: }
3111: static PetscErrorCode PCSetUp_SysPFMG(PC pc)
3112: {
3113: PC_SysPFMG *ex = (PC_SysPFMG *)pc->data;
3114: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)pc->pmat->data;
3115: PetscBool flg;
3117: PetscFunctionBegin;
3118: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRESSTRUCT, &flg));
3119: PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use MATHYPRESSTRUCT with this preconditioner");
3121: /* create the hypre sstruct solver object and set its information */
3122: if (ex->ss_solver) PetscCallExternal(HYPRE_SStructSysPFMGDestroy, ex->ss_solver);
3123: PetscCallExternal(HYPRE_SStructSysPFMGCreate, ex->hcomm, &ex->ss_solver);
3124: PetscCallExternal(HYPRE_SStructSysPFMGSetZeroGuess, ex->ss_solver);
3125: PetscCallExternal(HYPRE_SStructSysPFMGSetup, ex->ss_solver, mx->ss_mat, mx->ss_b, mx->ss_x);
3126: PetscFunctionReturn(PETSC_SUCCESS);
3127: }
3129: /*MC
3130: PCSYSPFMG - the hypre SysPFMG multigrid solver
3132: Level: advanced
3134: Options Database Keys:
3135: + -pc_syspfmg_its <its> - number of iterations of SysPFMG to use as preconditioner
3136: . -pc_syspfmg_num_pre_relax <steps> - number of smoothing steps before coarse grid
3137: . -pc_syspfmg_num_post_relax <steps> - number of smoothing steps after coarse grid
3138: . -pc_syspfmg_tol <tol> - tolerance of SysPFMG
3139: - -pc_syspfmg_relax_type <Weighted-Jacobi,Red/Black-Gauss-Seidel> - relaxation type for the up and down cycles
3141: Notes:
3142: See `PCPFMG` for hypre's PFMG that works for a scalar PDE and `PCSMG`
3144: See `PCHYPRE` for hypre's BoomerAMG algebraic multigrid solver
3146: This is for CELL-centered descretizations
3148: This must be used with the `MATHYPRESSTRUCT` matrix type.
3150: 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`.
3152: .seealso: [](ch_ksp), `PCMG`, `MATHYPRESSTRUCT`, `PCPFMG`, `PCHYPRE`, `PCGAMG`, `PCSMG`
3153: M*/
3155: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
3156: {
3157: PC_SysPFMG *ex;
3159: PetscFunctionBegin;
3160: PetscCall(PetscNew(&ex));
3161: pc->data = ex;
3163: ex->its = 1;
3164: ex->tol = 1.e-8;
3165: ex->relax_type = 1;
3166: ex->num_pre_relax = 1;
3167: ex->num_post_relax = 1;
3169: pc->ops->setfromoptions = PCSetFromOptions_SysPFMG;
3170: pc->ops->view = PCView_SysPFMG;
3171: pc->ops->destroy = PCDestroy_SysPFMG;
3172: pc->ops->apply = PCApply_SysPFMG;
3173: pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
3174: pc->ops->setup = PCSetUp_SysPFMG;
3176: PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
3177: PetscHYPREInitialize();
3178: PetscCallExternal(HYPRE_SStructSysPFMGCreate, ex->hcomm, &ex->ss_solver);
3179: PetscFunctionReturn(PETSC_SUCCESS);
3180: }
3182: /* PC SMG */
3183: typedef struct {
3184: MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */
3185: HYPRE_StructSolver hsolver;
3186: PetscInt its; /* keep copy of SMG options used so may view them */
3187: PetscReal tol;
3188: PetscBool print_statistics;
3189: PetscInt num_pre_relax, num_post_relax;
3190: } PC_SMG;
3192: static PetscErrorCode PCDestroy_SMG(PC pc)
3193: {
3194: PC_SMG *ex = (PC_SMG *)pc->data;
3196: PetscFunctionBegin;
3197: if (ex->hsolver) PetscCallExternal(HYPRE_StructSMGDestroy, ex->hsolver);
3198: PetscCall(PetscCommRestoreComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
3199: PetscCall(PetscFree(pc->data));
3200: PetscFunctionReturn(PETSC_SUCCESS);
3201: }
3203: static PetscErrorCode PCView_SMG(PC pc, PetscViewer viewer)
3204: {
3205: PetscBool isascii;
3206: PC_SMG *ex = (PC_SMG *)pc->data;
3208: PetscFunctionBegin;
3209: PetscCall(PetscObjectTypeCompare((PetscObject)viewer, PETSCVIEWERASCII, &isascii));
3210: if (isascii) {
3211: PetscCall(PetscViewerASCIIPrintf(viewer, " HYPRE SMG preconditioning\n"));
3212: PetscCall(PetscViewerASCIIPrintf(viewer, " max iterations %" PetscInt_FMT "\n", ex->its));
3213: PetscCall(PetscViewerASCIIPrintf(viewer, " tolerance %g\n", ex->tol));
3214: PetscCall(PetscViewerASCIIPrintf(viewer, " number pre-relax %" PetscInt_FMT " post-relax %" PetscInt_FMT "\n", ex->num_pre_relax, ex->num_post_relax));
3215: }
3216: PetscFunctionReturn(PETSC_SUCCESS);
3217: }
3219: static PetscErrorCode PCSetFromOptions_SMG(PC pc, PetscOptionItems PetscOptionsObject)
3220: {
3221: PC_SMG *ex = (PC_SMG *)pc->data;
3223: PetscFunctionBegin;
3224: PetscOptionsHeadBegin(PetscOptionsObject, "SMG options");
3226: PetscCall(PetscOptionsInt("-pc_smg_its", "Number of iterations of SMG to use as preconditioner", "HYPRE_StructSMGSetMaxIter", ex->its, &ex->its, NULL));
3227: 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));
3228: 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));
3229: PetscCall(PetscOptionsReal("-pc_smg_tol", "Tolerance of SMG", "HYPRE_StructSMGSetTol", ex->tol, &ex->tol, NULL));
3231: PetscOptionsHeadEnd();
3232: PetscFunctionReturn(PETSC_SUCCESS);
3233: }
3235: static PetscErrorCode PCApply_SMG(PC pc, Vec x, Vec y)
3236: {
3237: PC_SMG *ex = (PC_SMG *)pc->data;
3238: PetscScalar *yy;
3239: const PetscScalar *xx;
3240: PetscInt ilower[3], iupper[3];
3241: HYPRE_Int hlower[3], hupper[3];
3242: Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data;
3244: PetscFunctionBegin;
3245: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
3246: PetscCall(DMDAGetCorners(mx->da, &ilower[0], &ilower[1], &ilower[2], &iupper[0], &iupper[1], &iupper[2]));
3247: /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
3248: iupper[0] += ilower[0] - 1;
3249: iupper[1] += ilower[1] - 1;
3250: iupper[2] += ilower[2] - 1;
3251: hlower[0] = (HYPRE_Int)ilower[0];
3252: hlower[1] = (HYPRE_Int)ilower[1];
3253: hlower[2] = (HYPRE_Int)ilower[2];
3254: hupper[0] = (HYPRE_Int)iupper[0];
3255: hupper[1] = (HYPRE_Int)iupper[1];
3256: hupper[2] = (HYPRE_Int)iupper[2];
3258: /* copy x values over to hypre */
3259: PetscCallExternal(HYPRE_StructVectorSetConstantValues, mx->hb, 0.0);
3260: PetscCall(VecGetArrayRead(x, &xx));
3261: PetscCallExternal(HYPRE_StructVectorSetBoxValues, mx->hb, hlower, hupper, (HYPRE_Complex *)xx);
3262: PetscCall(VecRestoreArrayRead(x, &xx));
3263: PetscCallExternal(HYPRE_StructVectorAssemble, mx->hb);
3264: PetscCallExternal(HYPRE_StructSMGSolve, ex->hsolver, mx->hmat, mx->hb, mx->hx);
3266: /* copy solution values back to PETSc */
3267: PetscCall(VecGetArray(y, &yy));
3268: PetscCallExternal(HYPRE_StructVectorGetBoxValues, mx->hx, hlower, hupper, (HYPRE_Complex *)yy);
3269: PetscCall(VecRestoreArray(y, &yy));
3270: PetscFunctionReturn(PETSC_SUCCESS);
3271: }
3273: 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)
3274: {
3275: PC_SMG *jac = (PC_SMG *)pc->data;
3276: HYPRE_Int oits;
3278: PetscFunctionBegin;
3279: PetscCall(PetscCitationsRegister(hypreCitation, &cite));
3280: PetscCallExternal(HYPRE_StructSMGSetMaxIter, jac->hsolver, its * jac->its);
3281: PetscCallExternal(HYPRE_StructSMGSetTol, jac->hsolver, rtol);
3283: PetscCall(PCApply_SMG(pc, b, y));
3284: PetscCallExternal(HYPRE_StructSMGGetNumIterations, jac->hsolver, &oits);
3285: *outits = oits;
3286: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
3287: else *reason = PCRICHARDSON_CONVERGED_RTOL;
3288: PetscCallExternal(HYPRE_StructSMGSetTol, jac->hsolver, jac->tol);
3289: PetscCallExternal(HYPRE_StructSMGSetMaxIter, jac->hsolver, jac->its);
3290: PetscFunctionReturn(PETSC_SUCCESS);
3291: }
3293: static PetscErrorCode PCSetUp_SMG(PC pc)
3294: {
3295: PetscInt i, dim;
3296: PC_SMG *ex = (PC_SMG *)pc->data;
3297: Mat_HYPREStruct *mx = (Mat_HYPREStruct *)pc->pmat->data;
3298: PetscBool flg;
3299: DMBoundaryType p[3];
3300: PetscInt M[3];
3302: PetscFunctionBegin;
3303: PetscCall(PetscObjectTypeCompare((PetscObject)pc->pmat, MATHYPRESTRUCT, &flg));
3304: PetscCheck(flg, PetscObjectComm((PetscObject)pc), PETSC_ERR_ARG_INCOMP, "Must use MATHYPRESTRUCT with this preconditioner");
3306: PetscCall(DMDAGetInfo(mx->da, &dim, &M[0], &M[1], &M[2], 0, 0, 0, 0, 0, &p[0], &p[1], &p[2], 0));
3307: // Check if power of 2 in periodic directions
3308: for (i = 0; i < dim; i++) {
3309: 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]);
3310: }
3312: /* create the hypre solver object and set its information */
3313: if (ex->hsolver) PetscCallExternal(HYPRE_StructSMGDestroy, ex->hsolver);
3314: PetscCallExternal(HYPRE_StructSMGCreate, ex->hcomm, &ex->hsolver);
3315: // The hypre options must be set here and not in SetFromOptions because it is created here!
3316: PetscCallExternal(HYPRE_StructSMGSetMaxIter, ex->hsolver, ex->its);
3317: PetscCallExternal(HYPRE_StructSMGSetNumPreRelax, ex->hsolver, ex->num_pre_relax);
3318: PetscCallExternal(HYPRE_StructSMGSetNumPostRelax, ex->hsolver, ex->num_post_relax);
3319: PetscCallExternal(HYPRE_StructSMGSetTol, ex->hsolver, ex->tol);
3321: PetscCallExternal(HYPRE_StructSMGSetup, ex->hsolver, mx->hmat, mx->hb, mx->hx);
3322: PetscCallExternal(HYPRE_StructSMGSetZeroGuess, ex->hsolver);
3323: PetscFunctionReturn(PETSC_SUCCESS);
3324: }
3326: /*MC
3327: PCSMG - the hypre (structured grid) SMG multigrid solver
3329: Level: advanced
3331: Options Database Keys:
3332: + -pc_smg_its <its> - number of iterations of SMG to use as preconditioner
3333: . -pc_smg_num_pre_relax <steps> - number of smoothing steps before coarse grid
3334: . -pc_smg_num_post_relax <steps> - number of smoothing steps after coarse grid
3335: - -pc_smg_tol <tol> - tolerance of SMG
3337: Notes:
3338: This is for CELL-centered descretizations
3340: This must be used with the `MATHYPRESTRUCT` `MatType`.
3342: This does not provide all the functionality of hypre's SMG solver, it supports only one block per process defined by a PETSc `DMDA`.
3344: See `PCSYSPFMG`, `PCSMG`, `PCPFMG`, and `PCHYPRE` for access to hypre's other preconditioners
3346: .seealso: `PCMG`, `MATHYPRESTRUCT`, `PCPFMG`, `PCSYSPFMG`, `PCHYPRE`, `PCGAMG`
3347: M*/
3349: PETSC_EXTERN PetscErrorCode PCCreate_SMG(PC pc)
3350: {
3351: PC_SMG *ex;
3353: PetscFunctionBegin;
3354: PetscCall(PetscNew(&ex));
3355: pc->data = ex;
3357: ex->its = 1;
3358: ex->tol = 1.e-8;
3359: ex->num_pre_relax = 1;
3360: ex->num_post_relax = 1;
3362: pc->ops->setfromoptions = PCSetFromOptions_SMG;
3363: pc->ops->view = PCView_SMG;
3364: pc->ops->destroy = PCDestroy_SMG;
3365: pc->ops->apply = PCApply_SMG;
3366: pc->ops->applyrichardson = PCApplyRichardson_SMG;
3367: pc->ops->setup = PCSetUp_SMG;
3369: PetscCall(PetscCommGetComm(PetscObjectComm((PetscObject)pc), &ex->hcomm));
3370: PetscHYPREInitialize();
3371: PetscCallExternal(HYPRE_StructSMGCreate, ex->hcomm, &ex->hsolver);
3372: PetscFunctionReturn(PETSC_SUCCESS);
3373: }