1 /*
2 Provides an interface to the LLNL package hypre
3 */
4
5 /* Must use hypre 2.0.0 or more recent. */
6
7 #include <petsc/private/pcimpl.h> /*I "petscpc.h" I*/
8 /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */
9 #include <petsc/private/matimpl.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>
15
16 static PetscBool cite = PETSC_FALSE;
17 static const char hypreCitation[] = "@manual{hypre-web-page,\n title = {{\\sl hypre}: High Performance Preconditioners},\n organization = {Lawrence Livermore National Laboratory},\n note = {\\url{https://computation.llnl.gov/projects/hypre-scalable-linear-solvers-multigrid-methods}}\n}\n";
18
19 /*
20 Private context (data structure) for the preconditioner.
21 */
22 typedef struct {
23 HYPRE_Solver hsolver;
24 Mat hpmat; /* MatHYPRE */
25
26 HYPRE_Int (*destroy)(HYPRE_Solver);
27 HYPRE_Int (*solve)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
28 HYPRE_Int (*setup)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
29
30 MPI_Comm comm_hypre;
31 char *hypre_type;
32
33 /* options for Pilut and BoomerAMG*/
34 PetscInt maxiter;
35 PetscReal tol;
36
37 /* options for Pilut */
38 PetscInt factorrowsize;
39
40 /* options for ParaSails */
41 PetscInt nlevels;
42 PetscReal threshold;
43 PetscReal filter;
44 PetscInt sym;
45 PetscReal loadbal;
46 PetscInt logging;
47 PetscInt ruse;
48 PetscInt symt;
49
50 /* options for BoomerAMG */
51 PetscBool printstatistics;
52
53 /* options for BoomerAMG */
54 PetscInt cycletype;
55 PetscInt maxlevels;
56 PetscReal strongthreshold;
57 PetscReal maxrowsum;
58 PetscInt gridsweeps[3];
59 PetscInt coarsentype;
60 PetscInt measuretype;
61 PetscInt smoothtype;
62 PetscInt smoothnumlevels;
63 PetscInt eu_level; /* Number of levels for ILU(k) in Euclid */
64 PetscReal eu_droptolerance; /* Drop tolerance for ILU(k) in Euclid */
65 PetscInt eu_bj; /* Defines use of Block Jacobi ILU in Euclid */
66 PetscInt relaxtype[3];
67 PetscReal relaxweight;
68 PetscReal outerrelaxweight;
69 PetscInt relaxorder;
70 PetscReal truncfactor;
71 PetscBool applyrichardson;
72 PetscInt pmax;
73 PetscInt interptype;
74 PetscInt agg_nl;
75 PetscInt agg_num_paths;
76 PetscBool nodal_relax;
77 PetscInt nodal_relax_levels;
78
79 PetscInt nodal_coarsening;
80 PetscInt nodal_coarsening_diag;
81 PetscInt vec_interp_variant;
82 PetscInt vec_interp_qmax;
83 PetscBool vec_interp_smooth;
84 PetscInt interp_refine;
85
86 HYPRE_IJVector *hmnull;
87 HYPRE_ParVector *phmnull; /* near null space passed to hypre */
88 PetscInt n_hmnull;
89 Vec hmnull_constant;
90 HYPRE_Complex **hmnull_hypre_data_array; /* this is the space in hmnull that was allocated by hypre, it is restored to hypre just before freeing the phmnull vectors */
91
92 /* options for AS (Auxiliary Space preconditioners) */
93 PetscInt as_print;
94 PetscInt as_max_iter;
95 PetscReal as_tol;
96 PetscInt as_relax_type;
97 PetscInt as_relax_times;
98 PetscReal as_relax_weight;
99 PetscReal as_omega;
100 PetscInt as_amg_alpha_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for vector Poisson (AMS) or Curl problem (ADS) */
101 PetscReal as_amg_alpha_theta; /* AMG strength for vector Poisson (AMS) or Curl problem (ADS) */
102 PetscInt as_amg_beta_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for scalar Poisson (AMS) or vector Poisson (ADS) */
103 PetscReal as_amg_beta_theta; /* AMG strength for scalar Poisson (AMS) or vector Poisson (ADS) */
104 PetscInt ams_cycle_type;
105 PetscInt ads_cycle_type;
106
107 /* additional data */
108 Mat G; /* MatHYPRE */
109 Mat C; /* MatHYPRE */
110 Mat alpha_Poisson; /* MatHYPRE */
111 Mat beta_Poisson; /* MatHYPRE */
112
113 /* extra information for AMS */
114 PetscInt dim; /* geometrical dimension */
115 HYPRE_IJVector coords[3];
116 HYPRE_IJVector constants[3];
117 Mat RT_PiFull, RT_Pi[3];
118 Mat ND_PiFull, ND_Pi[3];
119 PetscBool ams_beta_is_zero;
120 PetscBool ams_beta_is_zero_part;
121 PetscInt ams_proj_freq;
122 } PC_HYPRE;
123
PCHYPREGetSolver(PC pc,HYPRE_Solver * hsolver)124 PetscErrorCode PCHYPREGetSolver(PC pc,HYPRE_Solver *hsolver)
125 {
126 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
127
128 PetscFunctionBegin;
129 *hsolver = jac->hsolver;
130 PetscFunctionReturn(0);
131 }
132
133 /*
134 Matrices with AIJ format are created IN PLACE with using (I,J,data) from BoomerAMG. Since the data format in hypre_ParCSRMatrix
135 is different from that used in PETSc, the original hypre_ParCSRMatrix can not be used any more after call this routine.
136 It is used in PCHMG. Other users should avoid using this function.
137 */
PCGetCoarseOperators_BoomerAMG(PC pc,PetscInt * nlevels,Mat * operators[])138 static PetscErrorCode PCGetCoarseOperators_BoomerAMG(PC pc,PetscInt *nlevels,Mat *operators[])
139 {
140 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
141 PetscBool same = PETSC_FALSE;
142 PetscErrorCode ierr;
143 PetscInt num_levels,l;
144 Mat *mattmp;
145 hypre_ParCSRMatrix **A_array;
146
147 PetscFunctionBegin;
148 ierr = PetscStrcmp(jac->hypre_type,"boomeramg",&same);CHKERRQ(ierr);
149 if (!same) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_NOTSAMETYPE,"Hypre type is not BoomerAMG \n");
150 num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData*) (jac->hsolver));
151 ierr = PetscMalloc1(num_levels,&mattmp);CHKERRQ(ierr);
152 A_array = hypre_ParAMGDataAArray((hypre_ParAMGData*) (jac->hsolver));
153 for (l=1; l<num_levels; l++) {
154 ierr = MatCreateFromParCSR(A_array[l],MATAIJ,PETSC_OWN_POINTER, &(mattmp[num_levels-1-l]));CHKERRQ(ierr);
155 /* We want to own the data, and HYPRE can not touch this matrix any more */
156 A_array[l] = NULL;
157 }
158 *nlevels = num_levels;
159 *operators = mattmp;
160 PetscFunctionReturn(0);
161 }
162
163 /*
164 Matrices with AIJ format are created IN PLACE with using (I,J,data) from BoomerAMG. Since the data format in hypre_ParCSRMatrix
165 is different from that used in PETSc, the original hypre_ParCSRMatrix can not be used any more after call this routine.
166 It is used in PCHMG. Other users should avoid using this function.
167 */
PCGetInterpolations_BoomerAMG(PC pc,PetscInt * nlevels,Mat * interpolations[])168 static PetscErrorCode PCGetInterpolations_BoomerAMG(PC pc,PetscInt *nlevels,Mat *interpolations[])
169 {
170 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
171 PetscBool same = PETSC_FALSE;
172 PetscErrorCode ierr;
173 PetscInt num_levels,l;
174 Mat *mattmp;
175 hypre_ParCSRMatrix **P_array;
176
177 PetscFunctionBegin;
178 ierr = PetscStrcmp(jac->hypre_type,"boomeramg",&same);CHKERRQ(ierr);
179 if (!same) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_NOTSAMETYPE,"Hypre type is not BoomerAMG \n");
180 num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData*) (jac->hsolver));
181 ierr = PetscMalloc1(num_levels,&mattmp);CHKERRQ(ierr);
182 P_array = hypre_ParAMGDataPArray((hypre_ParAMGData*) (jac->hsolver));
183 for (l=1; l<num_levels; l++) {
184 ierr = MatCreateFromParCSR(P_array[num_levels-1-l],MATAIJ,PETSC_OWN_POINTER, &(mattmp[l-1]));CHKERRQ(ierr);
185 /* We want to own the data, and HYPRE can not touch this matrix any more */
186 P_array[num_levels-1-l] = NULL;
187 }
188 *nlevels = num_levels;
189 *interpolations = mattmp;
190 PetscFunctionReturn(0);
191 }
192
193 /* Resets (frees) Hypre's representation of the near null space */
PCHYPREResetNearNullSpace_Private(PC pc)194 static PetscErrorCode PCHYPREResetNearNullSpace_Private(PC pc)
195 {
196 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
197 PetscInt i;
198 PetscErrorCode ierr;
199
200 PetscFunctionBegin;
201 for (i=0; i<jac->n_hmnull; i++) {
202 PETSC_UNUSED HYPRE_Complex *harray;
203 VecHYPRE_ParVectorReplacePointer(jac->hmnull[i],jac->hmnull_hypre_data_array[i],harray);
204 PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->hmnull[i]));
205 }
206 ierr = PetscFree(jac->hmnull);CHKERRQ(ierr);
207 ierr = PetscFree(jac->hmnull_hypre_data_array);CHKERRQ(ierr);
208 ierr = PetscFree(jac->phmnull);CHKERRQ(ierr);
209 ierr = VecDestroy(&jac->hmnull_constant);CHKERRQ(ierr);
210 jac->n_hmnull = 0;
211 PetscFunctionReturn(0);
212 }
213
PCSetUp_HYPRE(PC pc)214 static PetscErrorCode PCSetUp_HYPRE(PC pc)
215 {
216 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
217 Mat_HYPRE *hjac;
218 HYPRE_ParCSRMatrix hmat;
219 HYPRE_ParVector bv,xv;
220 PetscBool ishypre;
221 PetscErrorCode ierr;
222
223 PetscFunctionBegin;
224 if (!jac->hypre_type) {
225 ierr = PCHYPRESetType(pc,"boomeramg");CHKERRQ(ierr);
226 }
227
228 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRE,&ishypre);CHKERRQ(ierr);
229 if (!ishypre) {
230 ierr = MatDestroy(&jac->hpmat);CHKERRQ(ierr);
231 ierr = MatConvert(pc->pmat,MATHYPRE,MAT_INITIAL_MATRIX,&jac->hpmat);CHKERRQ(ierr);
232 } else {
233 ierr = PetscObjectReference((PetscObject)pc->pmat);CHKERRQ(ierr);
234 ierr = MatDestroy(&jac->hpmat);CHKERRQ(ierr);
235 jac->hpmat = pc->pmat;
236 }
237 hjac = (Mat_HYPRE*)(jac->hpmat->data);
238
239 /* special case for BoomerAMG */
240 if (jac->setup == HYPRE_BoomerAMGSetup) {
241 MatNullSpace mnull;
242 PetscBool has_const;
243 PetscInt bs,nvec,i;
244 const Vec *vecs;
245 HYPRE_Complex *petscvecarray;
246
247 ierr = MatGetBlockSize(pc->pmat,&bs);CHKERRQ(ierr);
248 if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
249 ierr = MatGetNearNullSpace(pc->mat, &mnull);CHKERRQ(ierr);
250 if (mnull) {
251 ierr = PCHYPREResetNearNullSpace_Private(pc);CHKERRQ(ierr);
252 ierr = MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);CHKERRQ(ierr);
253 ierr = PetscMalloc1(nvec+1,&jac->hmnull);CHKERRQ(ierr);
254 ierr = PetscMalloc1(nvec+1,&jac->hmnull_hypre_data_array);CHKERRQ(ierr);
255 ierr = PetscMalloc1(nvec+1,&jac->phmnull);CHKERRQ(ierr);
256 for (i=0; i<nvec; i++) {
257 ierr = VecHYPRE_IJVectorCreate(vecs[i],&jac->hmnull[i]);CHKERRQ(ierr);
258 ierr = VecGetArrayRead(vecs[i],(const PetscScalar **)&petscvecarray);CHKERRQ(ierr);
259 VecHYPRE_ParVectorReplacePointer(jac->hmnull[i],petscvecarray,jac->hmnull_hypre_data_array[i]);
260 ierr = VecRestoreArrayRead(vecs[i],(const PetscScalar **)&petscvecarray);CHKERRQ(ierr);
261 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[i],(void**)&jac->phmnull[i]));
262 }
263 if (has_const) {
264 ierr = MatCreateVecs(pc->pmat,&jac->hmnull_constant,NULL);CHKERRQ(ierr);
265 ierr = VecSet(jac->hmnull_constant,1);CHKERRQ(ierr);
266 ierr = VecNormalize(jac->hmnull_constant,NULL);
267 ierr = VecHYPRE_IJVectorCreate(jac->hmnull_constant,&jac->hmnull[nvec]);CHKERRQ(ierr);
268 ierr = VecGetArrayRead(jac->hmnull_constant,(const PetscScalar **)&petscvecarray);CHKERRQ(ierr);
269 VecHYPRE_ParVectorReplacePointer(jac->hmnull[nvec],petscvecarray,jac->hmnull_hypre_data_array[nvec]);
270 ierr = VecRestoreArrayRead(jac->hmnull_constant,(const PetscScalar **)&petscvecarray);CHKERRQ(ierr);
271 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[nvec],(void**)&jac->phmnull[nvec]));
272 nvec++;
273 }
274 PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVectors,(jac->hsolver,nvec,jac->phmnull));
275 jac->n_hmnull = nvec;
276 }
277 }
278
279 /* special case for AMS */
280 if (jac->setup == HYPRE_AMSSetup) {
281 Mat_HYPRE *hm;
282 HYPRE_ParCSRMatrix parcsr;
283 if (!jac->coords[0] && !jac->constants[0] && !(jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]))) {
284 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs either the coordinate vectors via PCSetCoordinates() or the edge constant vectors via PCHYPRESetEdgeConstantVectors() or the interpolation matrix via PCHYPRESetInterpolations");
285 }
286 if (jac->dim) {
287 PetscStackCallStandard(HYPRE_AMSSetDimension,(jac->hsolver,jac->dim));
288 }
289 if (jac->constants[0]) {
290 HYPRE_ParVector ozz,zoz,zzo = NULL;
291 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[0],(void**)(&ozz)));
292 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[1],(void**)(&zoz)));
293 if (jac->constants[2]) {
294 PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[2],(void**)(&zzo)));
295 }
296 PetscStackCallStandard(HYPRE_AMSSetEdgeConstantVectors,(jac->hsolver,ozz,zoz,zzo));
297 }
298 if (jac->coords[0]) {
299 HYPRE_ParVector coords[3];
300 coords[0] = NULL;
301 coords[1] = NULL;
302 coords[2] = NULL;
303 if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0],(void**)(&coords[0])));
304 if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1],(void**)(&coords[1])));
305 if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2],(void**)(&coords[2])));
306 PetscStackCallStandard(HYPRE_AMSSetCoordinateVectors,(jac->hsolver,coords[0],coords[1],coords[2]));
307 }
308 if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient");
309 hm = (Mat_HYPRE*)(jac->G->data);
310 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
311 PetscStackCallStandard(HYPRE_AMSSetDiscreteGradient,(jac->hsolver,parcsr));
312 if (jac->alpha_Poisson) {
313 hm = (Mat_HYPRE*)(jac->alpha_Poisson->data);
314 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
315 PetscStackCallStandard(HYPRE_AMSSetAlphaPoissonMatrix,(jac->hsolver,parcsr));
316 }
317 if (jac->ams_beta_is_zero) {
318 PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,NULL));
319 } else if (jac->beta_Poisson) {
320 hm = (Mat_HYPRE*)(jac->beta_Poisson->data);
321 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
322 PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,parcsr));
323 }
324 if (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1])) {
325 PetscInt i;
326 HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3];
327 if (jac->ND_PiFull) {
328 hm = (Mat_HYPRE*)(jac->ND_PiFull->data);
329 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&nd_parcsrfull)));
330 } else {
331 nd_parcsrfull = NULL;
332 }
333 for (i=0;i<3;++i) {
334 if (jac->ND_Pi[i]) {
335 hm = (Mat_HYPRE*)(jac->ND_Pi[i]->data);
336 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&nd_parcsr[i])));
337 } else {
338 nd_parcsr[i] = NULL;
339 }
340 }
341 PetscStackCallStandard(HYPRE_AMSSetInterpolations,(jac->hsolver,nd_parcsrfull,nd_parcsr[0],nd_parcsr[1],nd_parcsr[2]));
342 }
343 }
344 /* special case for ADS */
345 if (jac->setup == HYPRE_ADSSetup) {
346 Mat_HYPRE *hm;
347 HYPRE_ParCSRMatrix parcsr;
348 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])))) {
349 SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs either the coordinate vectors via PCSetCoordinates() or the interpolation matrices via PCHYPRESetInterpolations");
350 }
351 else if (!jac->coords[1] || !jac->coords[2]) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner has been designed for three dimensional problems! For two dimensional problems, use HYPRE AMS instead");
352 if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient");
353 if (!jac->C) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs the discrete curl operator via PCHYPRESetDiscreteGradient");
354 if (jac->coords[0]) {
355 HYPRE_ParVector coords[3];
356 coords[0] = NULL;
357 coords[1] = NULL;
358 coords[2] = NULL;
359 if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0],(void**)(&coords[0])));
360 if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1],(void**)(&coords[1])));
361 if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2],(void**)(&coords[2])));
362 PetscStackCallStandard(HYPRE_ADSSetCoordinateVectors,(jac->hsolver,coords[0],coords[1],coords[2]));
363 }
364 hm = (Mat_HYPRE*)(jac->G->data);
365 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
366 PetscStackCallStandard(HYPRE_ADSSetDiscreteGradient,(jac->hsolver,parcsr));
367 hm = (Mat_HYPRE*)(jac->C->data);
368 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
369 PetscStackCallStandard(HYPRE_ADSSetDiscreteCurl,(jac->hsolver,parcsr));
370 if ((jac->RT_PiFull || (jac->RT_Pi[0] && jac->RT_Pi[1])) && (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]))) {
371 PetscInt i;
372 HYPRE_ParCSRMatrix rt_parcsrfull, rt_parcsr[3];
373 HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3];
374 if (jac->RT_PiFull) {
375 hm = (Mat_HYPRE*)(jac->RT_PiFull->data);
376 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&rt_parcsrfull)));
377 } else {
378 rt_parcsrfull = NULL;
379 }
380 for (i=0;i<3;++i) {
381 if (jac->RT_Pi[i]) {
382 hm = (Mat_HYPRE*)(jac->RT_Pi[i]->data);
383 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&rt_parcsr[i])));
384 } else {
385 rt_parcsr[i] = NULL;
386 }
387 }
388 if (jac->ND_PiFull) {
389 hm = (Mat_HYPRE*)(jac->ND_PiFull->data);
390 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&nd_parcsrfull)));
391 } else {
392 nd_parcsrfull = NULL;
393 }
394 for (i=0;i<3;++i) {
395 if (jac->ND_Pi[i]) {
396 hm = (Mat_HYPRE*)(jac->ND_Pi[i]->data);
397 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&nd_parcsr[i])));
398 } else {
399 nd_parcsr[i] = NULL;
400 }
401 }
402 PetscStackCallStandard(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]));
403 }
404 }
405 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
406 PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b,(void**)&bv));
407 PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x,(void**)&xv));
408 PetscStackCallStandard(jac->setup,(jac->hsolver,hmat,bv,xv));
409 PetscFunctionReturn(0);
410 }
411
PCApply_HYPRE(PC pc,Vec b,Vec x)412 static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
413 {
414 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
415 Mat_HYPRE *hjac = (Mat_HYPRE*)(jac->hpmat->data);
416 PetscErrorCode ierr;
417 HYPRE_ParCSRMatrix hmat;
418 HYPRE_Complex *xv,*sxv;
419 HYPRE_Complex *bv,*sbv;
420 HYPRE_ParVector jbv,jxv;
421 PetscInt hierr;
422
423 PetscFunctionBegin;
424 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr);
425 if (!jac->applyrichardson) {ierr = VecSet(x,0.0);CHKERRQ(ierr);}
426 ierr = VecGetArrayRead(b,(const PetscScalar **)&bv);CHKERRQ(ierr);
427 ierr = VecGetArray(x,(PetscScalar **)&xv);CHKERRQ(ierr);
428 VecHYPRE_ParVectorReplacePointer(hjac->b,bv,sbv);
429 VecHYPRE_ParVectorReplacePointer(hjac->x,xv,sxv);
430
431 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
432 PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b,(void**)&jbv));
433 PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x,(void**)&jxv));
434 PetscStackCall("Hypre solve",hierr = (*jac->solve)(jac->hsolver,hmat,jbv,jxv);
435 if (hierr && hierr != HYPRE_ERROR_CONV) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
436 if (hierr) hypre__global_error = 0;);
437
438 if (jac->setup == HYPRE_AMSSetup && jac->ams_beta_is_zero_part) {
439 PetscStackCallStandard(HYPRE_AMSProjectOutGradients,(jac->hsolver,jxv));
440 }
441 VecHYPRE_ParVectorReplacePointer(hjac->b,sbv,bv);
442 VecHYPRE_ParVectorReplacePointer(hjac->x,sxv,xv);
443 ierr = VecRestoreArray(x,(PetscScalar **)&xv);CHKERRQ(ierr);
444 ierr = VecRestoreArrayRead(b,(const PetscScalar **)&bv);CHKERRQ(ierr);
445 PetscFunctionReturn(0);
446 }
447
PCReset_HYPRE(PC pc)448 static PetscErrorCode PCReset_HYPRE(PC pc)
449 {
450 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
451 PetscErrorCode ierr;
452
453 PetscFunctionBegin;
454 ierr = MatDestroy(&jac->hpmat);CHKERRQ(ierr);
455 ierr = MatDestroy(&jac->G);CHKERRQ(ierr);
456 ierr = MatDestroy(&jac->C);CHKERRQ(ierr);
457 ierr = MatDestroy(&jac->alpha_Poisson);CHKERRQ(ierr);
458 ierr = MatDestroy(&jac->beta_Poisson);CHKERRQ(ierr);
459 ierr = MatDestroy(&jac->RT_PiFull);CHKERRQ(ierr);
460 ierr = MatDestroy(&jac->RT_Pi[0]);CHKERRQ(ierr);
461 ierr = MatDestroy(&jac->RT_Pi[1]);CHKERRQ(ierr);
462 ierr = MatDestroy(&jac->RT_Pi[2]);CHKERRQ(ierr);
463 ierr = MatDestroy(&jac->ND_PiFull);CHKERRQ(ierr);
464 ierr = MatDestroy(&jac->ND_Pi[0]);CHKERRQ(ierr);
465 ierr = MatDestroy(&jac->ND_Pi[1]);CHKERRQ(ierr);
466 ierr = MatDestroy(&jac->ND_Pi[2]);CHKERRQ(ierr);
467 if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0])); jac->coords[0] = NULL;
468 if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1])); jac->coords[1] = NULL;
469 if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2])); jac->coords[2] = NULL;
470 if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0])); jac->constants[0] = NULL;
471 if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1])); jac->constants[1] = NULL;
472 if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2])); jac->constants[2] = NULL;
473 ierr = PCHYPREResetNearNullSpace_Private(pc);CHKERRQ(ierr);
474 jac->ams_beta_is_zero = PETSC_FALSE;
475 jac->dim = 0;
476 PetscFunctionReturn(0);
477 }
478
PCDestroy_HYPRE(PC pc)479 static PetscErrorCode PCDestroy_HYPRE(PC pc)
480 {
481 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
482 PetscErrorCode ierr;
483
484 PetscFunctionBegin;
485 ierr = PCReset_HYPRE(pc);CHKERRQ(ierr);
486 if (jac->destroy) PetscStackCallStandard(jac->destroy,(jac->hsolver));
487 ierr = PetscFree(jac->hypre_type);CHKERRQ(ierr);
488 if (jac->comm_hypre != MPI_COMM_NULL) { ierr = MPI_Comm_free(&(jac->comm_hypre));CHKERRQ(ierr);}
489 ierr = PetscFree(pc->data);CHKERRQ(ierr);
490
491 ierr = PetscObjectChangeTypeName((PetscObject)pc,0);CHKERRQ(ierr);
492 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);CHKERRQ(ierr);
493 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);CHKERRQ(ierr);
494 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetCoordinates_C",NULL);CHKERRQ(ierr);
495 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",NULL);CHKERRQ(ierr);
496 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",NULL);CHKERRQ(ierr);
497 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",NULL);CHKERRQ(ierr);
498 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetConstantEdgeVectors_C",NULL);CHKERRQ(ierr);
499 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",NULL);CHKERRQ(ierr);
500 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",NULL);CHKERRQ(ierr);
501 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",NULL);CHKERRQ(ierr);
502 PetscFunctionReturn(0);
503 }
504
505 /* --------------------------------------------------------------------------------------------*/
PCSetFromOptions_HYPRE_Pilut(PetscOptionItems * PetscOptionsObject,PC pc)506 static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PetscOptionItems *PetscOptionsObject,PC pc)
507 {
508 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
509 PetscErrorCode ierr;
510 PetscBool flag;
511
512 PetscFunctionBegin;
513 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE Pilut Options");CHKERRQ(ierr);
514 ierr = PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);CHKERRQ(ierr);
515 if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter));
516 ierr = PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);CHKERRQ(ierr);
517 if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol));
518 ierr = PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);CHKERRQ(ierr);
519 if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize));
520 ierr = PetscOptionsTail();CHKERRQ(ierr);
521 PetscFunctionReturn(0);
522 }
523
PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)524 static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
525 {
526 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
527 PetscErrorCode ierr;
528 PetscBool iascii;
529
530 PetscFunctionBegin;
531 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
532 if (iascii) {
533 ierr = PetscViewerASCIIPrintf(viewer," HYPRE Pilut preconditioning\n");CHKERRQ(ierr);
534 if (jac->maxiter != PETSC_DEFAULT) {
535 ierr = PetscViewerASCIIPrintf(viewer," maximum number of iterations %d\n",jac->maxiter);CHKERRQ(ierr);
536 } else {
537 ierr = PetscViewerASCIIPrintf(viewer," default maximum number of iterations \n");CHKERRQ(ierr);
538 }
539 if (jac->tol != PETSC_DEFAULT) {
540 ierr = PetscViewerASCIIPrintf(viewer," drop tolerance %g\n",(double)jac->tol);CHKERRQ(ierr);
541 } else {
542 ierr = PetscViewerASCIIPrintf(viewer," default drop tolerance \n");CHKERRQ(ierr);
543 }
544 if (jac->factorrowsize != PETSC_DEFAULT) {
545 ierr = PetscViewerASCIIPrintf(viewer," factor row size %d\n",jac->factorrowsize);CHKERRQ(ierr);
546 } else {
547 ierr = PetscViewerASCIIPrintf(viewer," default factor row size \n");CHKERRQ(ierr);
548 }
549 }
550 PetscFunctionReturn(0);
551 }
552
553 /* --------------------------------------------------------------------------------------------*/
PCSetFromOptions_HYPRE_Euclid(PetscOptionItems * PetscOptionsObject,PC pc)554 static PetscErrorCode PCSetFromOptions_HYPRE_Euclid(PetscOptionItems *PetscOptionsObject,PC pc)
555 {
556 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
557 PetscErrorCode ierr;
558 PetscBool flag;
559
560 PetscFunctionBegin;
561 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE Euclid Options");CHKERRQ(ierr);
562 ierr = PetscOptionsInt("-pc_hypre_euclid_level","Factorization levels","None",jac->eu_level,&jac->eu_level,&flag);CHKERRQ(ierr);
563 if (flag) PetscStackCallStandard(HYPRE_EuclidSetLevel,(jac->hsolver,jac->eu_level));
564 ierr = PetscOptionsTail();CHKERRQ(ierr);
565 PetscFunctionReturn(0);
566 }
567
PCView_HYPRE_Euclid(PC pc,PetscViewer viewer)568 static PetscErrorCode PCView_HYPRE_Euclid(PC pc,PetscViewer viewer)
569 {
570 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
571 PetscErrorCode ierr;
572 PetscBool iascii;
573
574 PetscFunctionBegin;
575 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
576 if (iascii) {
577 ierr = PetscViewerASCIIPrintf(viewer," HYPRE Euclid preconditioning\n");CHKERRQ(ierr);
578 if (jac->eu_level != PETSC_DEFAULT) {
579 ierr = PetscViewerASCIIPrintf(viewer," factorization levels %d\n",jac->eu_level);CHKERRQ(ierr);
580 } else {
581 ierr = PetscViewerASCIIPrintf(viewer," default factorization levels \n");CHKERRQ(ierr);
582 }
583 }
584 PetscFunctionReturn(0);
585 }
586
587 /* --------------------------------------------------------------------------------------------*/
588
PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)589 static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
590 {
591 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
592 Mat_HYPRE *hjac = (Mat_HYPRE*)(jac->hpmat->data);
593 PetscErrorCode ierr;
594 HYPRE_ParCSRMatrix hmat;
595 HYPRE_Complex *xv,*bv;
596 HYPRE_Complex *sbv,*sxv;
597 HYPRE_ParVector jbv,jxv;
598 PetscInt hierr;
599
600 PetscFunctionBegin;
601 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr);
602 ierr = VecSet(x,0.0);CHKERRQ(ierr);
603 ierr = VecGetArrayRead(b,(const PetscScalar**)&bv);CHKERRQ(ierr);
604 ierr = VecGetArray(x,(PetscScalar**)&xv);CHKERRQ(ierr);
605 VecHYPRE_ParVectorReplacePointer(hjac->b,bv,sbv);
606 VecHYPRE_ParVectorReplacePointer(hjac->x,xv,sxv);
607
608 PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
609 PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b,(void**)&jbv));
610 PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x,(void**)&jxv));
611
612 hierr = HYPRE_BoomerAMGSolveT(jac->hsolver,hmat,jbv,jxv);
613 /* error code of 1 in BoomerAMG merely means convergence not achieved */
614 if (hierr && (hierr != 1)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
615 if (hierr) hypre__global_error = 0;
616
617 VecHYPRE_ParVectorReplacePointer(hjac->b,sbv,bv);
618 VecHYPRE_ParVectorReplacePointer(hjac->x,sxv,xv);
619 ierr = VecRestoreArray(x,(PetscScalar**)&xv);CHKERRQ(ierr);
620 ierr = VecRestoreArrayRead(b,(const PetscScalar**)&bv);CHKERRQ(ierr);
621 PetscFunctionReturn(0);
622 }
623
624 /* static array length */
625 #define ALEN(a) (sizeof(a)/sizeof((a)[0]))
626
627 static const char *HYPREBoomerAMGCycleType[] = {"","V","W"};
628 static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
629 static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
630 /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
631 static const char *HYPREBoomerAMGSmoothType[] = {"Schwarz-smoothers","Pilut","ParaSails","Euclid"};
632 static const char *HYPREBoomerAMGRelaxType[] = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi",
633 "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi",
634 "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination",
635 "" /* 10 */, "" /* 11 */, "" /* 12 */, "l1-Gauss-Seidel" /* nonsymmetric */, "backward-l1-Gauss-Seidel" /* nonsymmetric */,
636 "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"};
637 static const char *HYPREBoomerAMGInterpType[] = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
638 "ext+i-cc", "standard", "standard-wts", "block", "block-wtd", "FF", "FF1"};
PCSetFromOptions_HYPRE_BoomerAMG(PetscOptionItems * PetscOptionsObject,PC pc)639 static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PetscOptionItems *PetscOptionsObject,PC pc)
640 {
641 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
642 PetscErrorCode ierr;
643 PetscInt bs,n,indx,level;
644 PetscBool flg, tmp_truth;
645 double tmpdbl, twodbl[2];
646
647 PetscFunctionBegin;
648 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE BoomerAMG Options");CHKERRQ(ierr);
649 ierr = PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);CHKERRQ(ierr);
650 if (flg) {
651 jac->cycletype = indx+1;
652 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
653 }
654 ierr = PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);CHKERRQ(ierr);
655 if (flg) {
656 if (jac->maxlevels < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
657 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
658 }
659 ierr = PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);CHKERRQ(ierr);
660 if (flg) {
661 if (jac->maxiter < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
662 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
663 }
664 ierr = PetscOptionsReal("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);CHKERRQ(ierr);
665 if (flg) {
666 if (jac->tol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance %g must be greater than or equal to zero",(double)jac->tol);
667 PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
668 }
669 bs = 1;
670 if (pc->pmat) {
671 ierr = MatGetBlockSize(pc->pmat,&bs);CHKERRQ(ierr);
672 }
673 ierr = PetscOptionsInt("-pc_hypre_boomeramg_numfunctions","Number of functions","HYPRE_BoomerAMGSetNumFunctions",bs,&bs,&flg);CHKERRQ(ierr);
674 if (flg) {
675 PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
676 }
677
678 ierr = PetscOptionsReal("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);CHKERRQ(ierr);
679 if (flg) {
680 if (jac->truncfactor < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Truncation factor %g must be great than or equal zero",(double)jac->truncfactor);
681 PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
682 }
683
684 ierr = PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);CHKERRQ(ierr);
685 if (flg) {
686 if (jac->pmax < 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"P_max %g must be greater than or equal to zero",(double)jac->pmax);
687 PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
688 }
689
690 ierr = PetscOptionsInt("-pc_hypre_boomeramg_agg_nl","Number of levels of aggressive coarsening","None",jac->agg_nl,&jac->agg_nl,&flg);CHKERRQ(ierr);
691 if (flg) {
692 if (jac->agg_nl < 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %g must be greater than or equal to zero",(double)jac->agg_nl);
693
694 PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
695 }
696
697 ierr = PetscOptionsInt("-pc_hypre_boomeramg_agg_num_paths","Number of paths for aggressive coarsening","None",jac->agg_num_paths,&jac->agg_num_paths,&flg);CHKERRQ(ierr);
698 if (flg) {
699 if (jac->agg_num_paths < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of paths %g must be greater than or equal to 1",(double)jac->agg_num_paths);
700
701 PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
702 }
703
704
705 ierr = PetscOptionsReal("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);CHKERRQ(ierr);
706 if (flg) {
707 if (jac->strongthreshold < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Strong threshold %g must be great than or equal zero",(double)jac->strongthreshold);
708 PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
709 }
710 ierr = PetscOptionsReal("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);CHKERRQ(ierr);
711 if (flg) {
712 if (jac->maxrowsum < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %g must be greater than zero",(double)jac->maxrowsum);
713 if (jac->maxrowsum > 1.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %g must be less than or equal one",(double)jac->maxrowsum);
714 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
715 }
716
717 /* Grid sweeps */
718 ierr = PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);CHKERRQ(ierr);
719 if (flg) {
720 PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx));
721 /* modify the jac structure so we can view the updated options with PC_View */
722 jac->gridsweeps[0] = indx;
723 jac->gridsweeps[1] = indx;
724 /*defaults coarse to 1 */
725 jac->gridsweeps[2] = 1;
726 }
727 ierr = PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen","Use a nodal based coarsening 1-6","HYPRE_BoomerAMGSetNodal",jac->nodal_coarsening,&jac->nodal_coarsening,&flg);CHKERRQ(ierr);
728 if (flg) {
729 PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,jac->nodal_coarsening));
730 }
731 ierr = 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);CHKERRQ(ierr);
732 if (flg) {
733 PetscStackCallStandard(HYPRE_BoomerAMGSetNodalDiag,(jac->hsolver,jac->nodal_coarsening_diag));
734 }
735 ierr = PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_variant","Variant of algorithm 1-3","HYPRE_BoomerAMGSetInterpVecVariant",jac->vec_interp_variant, &jac->vec_interp_variant,&flg);CHKERRQ(ierr);
736 if (flg) {
737 PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecVariant,(jac->hsolver,jac->vec_interp_variant));
738 }
739 ierr = 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);CHKERRQ(ierr);
740 if (flg) {
741 PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecQMax,(jac->hsolver,jac->vec_interp_qmax));
742 }
743 ierr = PetscOptionsBool("-pc_hypre_boomeramg_vec_interp_smooth","Whether to smooth the interpolation vectors","HYPRE_BoomerAMGSetSmoothInterpVectors",jac->vec_interp_smooth, &jac->vec_interp_smooth,&flg);CHKERRQ(ierr);
744 if (flg) {
745 PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothInterpVectors,(jac->hsolver,jac->vec_interp_smooth));
746 }
747 ierr = PetscOptionsInt("-pc_hypre_boomeramg_interp_refine","Preprocess the interpolation matrix through iterative weight refinement","HYPRE_BoomerAMGSetInterpRefine",jac->interp_refine, &jac->interp_refine,&flg);CHKERRQ(ierr);
748 if (flg) {
749 PetscStackCallStandard(HYPRE_BoomerAMGSetInterpRefine,(jac->hsolver,jac->interp_refine));
750 }
751 ierr = PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);CHKERRQ(ierr);
752 if (flg) {
753 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1));
754 jac->gridsweeps[0] = indx;
755 }
756 ierr = PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);CHKERRQ(ierr);
757 if (flg) {
758 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2));
759 jac->gridsweeps[1] = indx;
760 }
761 ierr = PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);CHKERRQ(ierr);
762 if (flg) {
763 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3));
764 jac->gridsweeps[2] = indx;
765 }
766
767 /* Smooth type */
768 ierr = PetscOptionsEList("-pc_hypre_boomeramg_smooth_type","Enable more complex smoothers","None",HYPREBoomerAMGSmoothType,ALEN(HYPREBoomerAMGSmoothType),HYPREBoomerAMGSmoothType[0],&indx,&flg);
769 if (flg) {
770 jac->smoothtype = indx;
771 PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,indx+6));
772 jac->smoothnumlevels = 25;
773 PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,25));
774 }
775
776 /* Number of smoothing levels */
777 ierr = PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_levels","Number of levels on which more complex smoothers are used","None",25,&indx,&flg);CHKERRQ(ierr);
778 if (flg && (jac->smoothtype != -1)) {
779 jac->smoothnumlevels = indx;
780 PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,indx));
781 }
782
783 /* Number of levels for ILU(k) for Euclid */
784 ierr = PetscOptionsInt("-pc_hypre_boomeramg_eu_level","Number of levels for ILU(k) in Euclid smoother","None",0,&indx,&flg);CHKERRQ(ierr);
785 if (flg && (jac->smoothtype == 3)) {
786 jac->eu_level = indx;
787 PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,indx));
788 }
789
790 /* Filter for ILU(k) for Euclid */
791 double droptolerance;
792 ierr = PetscOptionsReal("-pc_hypre_boomeramg_eu_droptolerance","Drop tolerance for ILU(k) in Euclid smoother","None",0,&droptolerance,&flg);CHKERRQ(ierr);
793 if (flg && (jac->smoothtype == 3)) {
794 jac->eu_droptolerance = droptolerance;
795 PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,droptolerance));
796 }
797
798 /* Use Block Jacobi ILUT for Euclid */
799 ierr = PetscOptionsBool("-pc_hypre_boomeramg_eu_bj", "Use Block Jacobi for ILU in Euclid smoother?", "None", PETSC_FALSE, &tmp_truth, &flg);CHKERRQ(ierr);
800 if (flg && (jac->smoothtype == 3)) {
801 jac->eu_bj = tmp_truth;
802 PetscStackCallStandard(HYPRE_BoomerAMGSetEuBJ,(jac->hsolver,jac->eu_bj));
803 }
804
805 /* Relax type */
806 ierr = PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);CHKERRQ(ierr);
807 if (flg) {
808 jac->relaxtype[0] = jac->relaxtype[1] = indx;
809 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx));
810 /* by default, coarse type set to 9 */
811 jac->relaxtype[2] = 9;
812 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, 9, 3));
813 }
814 ierr = PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);CHKERRQ(ierr);
815 if (flg) {
816 jac->relaxtype[0] = indx;
817 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1));
818 }
819 ierr = PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);CHKERRQ(ierr);
820 if (flg) {
821 jac->relaxtype[1] = indx;
822 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2));
823 }
824 ierr = PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);CHKERRQ(ierr);
825 if (flg) {
826 jac->relaxtype[2] = indx;
827 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3));
828 }
829
830 /* Relaxation Weight */
831 ierr = 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);CHKERRQ(ierr);
832 if (flg) {
833 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl));
834 jac->relaxweight = tmpdbl;
835 }
836
837 n = 2;
838 twodbl[0] = twodbl[1] = 1.0;
839 ierr = PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);CHKERRQ(ierr);
840 if (flg) {
841 if (n == 2) {
842 indx = (int)PetscAbsReal(twodbl[1]);
843 PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx));
844 } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relax weight level: you must provide 2 values separated by a comma (and no space), you provided %d",n);
845 }
846
847 /* Outer relaxation Weight */
848 ierr = 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);CHKERRQ(ierr);
849 if (flg) {
850 PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl));
851 jac->outerrelaxweight = tmpdbl;
852 }
853
854 n = 2;
855 twodbl[0] = twodbl[1] = 1.0;
856 ierr = PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);CHKERRQ(ierr);
857 if (flg) {
858 if (n == 2) {
859 indx = (int)PetscAbsReal(twodbl[1]);
860 PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx));
861 } else SETERRQ1(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 %d",n);
862 }
863
864 /* the Relax Order */
865 ierr = PetscOptionsBool("-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", PETSC_FALSE, &tmp_truth, &flg);CHKERRQ(ierr);
866
867 if (flg && tmp_truth) {
868 jac->relaxorder = 0;
869 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
870 }
871 ierr = PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);CHKERRQ(ierr);
872 if (flg) {
873 jac->measuretype = indx;
874 PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
875 }
876 /* update list length 3/07 */
877 ierr = PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);CHKERRQ(ierr);
878 if (flg) {
879 jac->coarsentype = indx;
880 PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
881 }
882
883 /* new 3/07 */
884 ierr = PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);CHKERRQ(ierr);
885 if (flg) {
886 jac->interptype = indx;
887 PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
888 }
889
890 ierr = PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);CHKERRQ(ierr);
891 if (flg) {
892 level = 3;
893 ierr = PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);CHKERRQ(ierr);
894
895 jac->printstatistics = PETSC_TRUE;
896 PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level));
897 }
898
899 ierr = PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);CHKERRQ(ierr);
900 if (flg) {
901 level = 3;
902 ierr = PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);CHKERRQ(ierr);
903
904 jac->printstatistics = PETSC_TRUE;
905 PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level));
906 }
907
908 ierr = PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);CHKERRQ(ierr);
909 if (flg && tmp_truth) {
910 PetscInt tmp_int;
911 ierr = PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);CHKERRQ(ierr);
912 if (flg) jac->nodal_relax_levels = tmp_int;
913 PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6));
914 PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1));
915 PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0));
916 PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels));
917 }
918
919 ierr = PetscOptionsTail();CHKERRQ(ierr);
920 PetscFunctionReturn(0);
921 }
922
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)923 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)
924 {
925 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
926 PetscErrorCode ierr;
927 HYPRE_Int oits;
928
929 PetscFunctionBegin;
930 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr);
931 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter));
932 PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol));
933 jac->applyrichardson = PETSC_TRUE;
934 ierr = PCApply_HYPRE(pc,b,y);CHKERRQ(ierr);
935 jac->applyrichardson = PETSC_FALSE;
936 PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,&oits));
937 *outits = oits;
938 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
939 else *reason = PCRICHARDSON_CONVERGED_RTOL;
940 PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
941 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
942 PetscFunctionReturn(0);
943 }
944
945
PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)946 static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
947 {
948 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
949 PetscErrorCode ierr;
950 PetscBool iascii;
951
952 PetscFunctionBegin;
953 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
954 if (iascii) {
955 ierr = PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG preconditioning\n");CHKERRQ(ierr);
956 ierr = PetscViewerASCIIPrintf(viewer," Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);CHKERRQ(ierr);
957 ierr = PetscViewerASCIIPrintf(viewer," Maximum number of levels %D\n",jac->maxlevels);CHKERRQ(ierr);
958 ierr = PetscViewerASCIIPrintf(viewer," Maximum number of iterations PER hypre call %D\n",jac->maxiter);CHKERRQ(ierr);
959 ierr = PetscViewerASCIIPrintf(viewer," Convergence tolerance PER hypre call %g\n",(double)jac->tol);CHKERRQ(ierr);
960 ierr = PetscViewerASCIIPrintf(viewer," Threshold for strong coupling %g\n",(double)jac->strongthreshold);CHKERRQ(ierr);
961 ierr = PetscViewerASCIIPrintf(viewer," Interpolation truncation factor %g\n",(double)jac->truncfactor);CHKERRQ(ierr);
962 ierr = PetscViewerASCIIPrintf(viewer," Interpolation: max elements per row %D\n",jac->pmax);CHKERRQ(ierr);
963 if (jac->interp_refine) {
964 ierr = PetscViewerASCIIPrintf(viewer," Interpolation: number of steps of weighted refinement %D\n",jac->interp_refine);CHKERRQ(ierr);
965 }
966 ierr = PetscViewerASCIIPrintf(viewer," Number of levels of aggressive coarsening %D\n",jac->agg_nl);CHKERRQ(ierr);
967 ierr = PetscViewerASCIIPrintf(viewer," Number of paths for aggressive coarsening %D\n",jac->agg_num_paths);CHKERRQ(ierr);
968
969 ierr = PetscViewerASCIIPrintf(viewer," Maximum row sums %g\n",(double)jac->maxrowsum);CHKERRQ(ierr);
970
971 ierr = PetscViewerASCIIPrintf(viewer," Sweeps down %D\n",jac->gridsweeps[0]);CHKERRQ(ierr);
972 ierr = PetscViewerASCIIPrintf(viewer," Sweeps up %D\n",jac->gridsweeps[1]);CHKERRQ(ierr);
973 ierr = PetscViewerASCIIPrintf(viewer," Sweeps on coarse %D\n",jac->gridsweeps[2]);CHKERRQ(ierr);
974
975 ierr = PetscViewerASCIIPrintf(viewer," Relax down %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[0]]);CHKERRQ(ierr);
976 ierr = PetscViewerASCIIPrintf(viewer," Relax up %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[1]]);CHKERRQ(ierr);
977 ierr = PetscViewerASCIIPrintf(viewer," Relax on coarse %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[2]]);CHKERRQ(ierr);
978
979 ierr = PetscViewerASCIIPrintf(viewer," Relax weight (all) %g\n",(double)jac->relaxweight);CHKERRQ(ierr);
980 ierr = PetscViewerASCIIPrintf(viewer," Outer relax weight (all) %g\n",(double)jac->outerrelaxweight);CHKERRQ(ierr);
981
982 if (jac->relaxorder) {
983 ierr = PetscViewerASCIIPrintf(viewer," Using CF-relaxation\n");CHKERRQ(ierr);
984 } else {
985 ierr = PetscViewerASCIIPrintf(viewer," Not using CF-relaxation\n");CHKERRQ(ierr);
986 }
987 if (jac->smoothtype!=-1) {
988 ierr = PetscViewerASCIIPrintf(viewer," Smooth type %s\n",HYPREBoomerAMGSmoothType[jac->smoothtype]);CHKERRQ(ierr);
989 ierr = PetscViewerASCIIPrintf(viewer," Smooth num levels %D\n",jac->smoothnumlevels);CHKERRQ(ierr);
990 } else {
991 ierr = PetscViewerASCIIPrintf(viewer," Not using more complex smoothers.\n");CHKERRQ(ierr);
992 }
993 if (jac->smoothtype==3) {
994 ierr = PetscViewerASCIIPrintf(viewer," Euclid ILU(k) levels %D\n",jac->eu_level);CHKERRQ(ierr);
995 ierr = PetscViewerASCIIPrintf(viewer," Euclid ILU(k) drop tolerance %g\n",(double)jac->eu_droptolerance);CHKERRQ(ierr);
996 ierr = PetscViewerASCIIPrintf(viewer," Euclid ILU use Block-Jacobi? %D\n",jac->eu_bj);CHKERRQ(ierr);
997 }
998 ierr = PetscViewerASCIIPrintf(viewer," Measure type %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);CHKERRQ(ierr);
999 ierr = PetscViewerASCIIPrintf(viewer," Coarsen type %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);CHKERRQ(ierr);
1000 ierr = PetscViewerASCIIPrintf(viewer," Interpolation type %s\n",HYPREBoomerAMGInterpType[jac->interptype]);CHKERRQ(ierr);
1001 if (jac->nodal_coarsening) {
1002 ierr = PetscViewerASCIIPrintf(viewer," Using nodal coarsening (with HYPRE_BOOMERAMGSetNodal() %D\n",jac->nodal_coarsening);CHKERRQ(ierr);
1003 }
1004 if (jac->vec_interp_variant) {
1005 ierr = PetscViewerASCIIPrintf(viewer," HYPRE_BoomerAMGSetInterpVecVariant() %D\n",jac->vec_interp_variant);CHKERRQ(ierr);
1006 ierr = PetscViewerASCIIPrintf(viewer," HYPRE_BoomerAMGSetInterpVecQMax() %D\n",jac->vec_interp_qmax);CHKERRQ(ierr);
1007 ierr = PetscViewerASCIIPrintf(viewer," HYPRE_BoomerAMGSetSmoothInterpVectors() %d\n",jac->vec_interp_smooth);CHKERRQ(ierr);
1008 }
1009 if (jac->nodal_relax) {
1010 ierr = PetscViewerASCIIPrintf(viewer," Using nodal relaxation via Schwarz smoothing on levels %D\n",jac->nodal_relax_levels);CHKERRQ(ierr);
1011 }
1012 }
1013 PetscFunctionReturn(0);
1014 }
1015
1016 /* --------------------------------------------------------------------------------------------*/
PCSetFromOptions_HYPRE_ParaSails(PetscOptionItems * PetscOptionsObject,PC pc)1017 static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PetscOptionItems *PetscOptionsObject,PC pc)
1018 {
1019 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1020 PetscErrorCode ierr;
1021 PetscInt indx;
1022 PetscBool flag;
1023 const char *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};
1024
1025 PetscFunctionBegin;
1026 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE ParaSails Options");CHKERRQ(ierr);
1027 ierr = PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);CHKERRQ(ierr);
1028 ierr = PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshold,&jac->threshold,&flag);CHKERRQ(ierr);
1029 if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshold,jac->nlevels));
1030
1031 ierr = PetscOptionsReal("-pc_hypre_parasails_filter","filter","None",jac->filter,&jac->filter,&flag);CHKERRQ(ierr);
1032 if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
1033
1034 ierr = PetscOptionsReal("-pc_hypre_parasails_loadbal","Load balance","None",jac->loadbal,&jac->loadbal,&flag);CHKERRQ(ierr);
1035 if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
1036
1037 ierr = PetscOptionsBool("-pc_hypre_parasails_logging","Print info to screen","None",(PetscBool)jac->logging,(PetscBool*)&jac->logging,&flag);CHKERRQ(ierr);
1038 if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
1039
1040 ierr = PetscOptionsBool("-pc_hypre_parasails_reuse","Reuse nonzero pattern in preconditioner","None",(PetscBool)jac->ruse,(PetscBool*)&jac->ruse,&flag);CHKERRQ(ierr);
1041 if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
1042
1043 ierr = PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);CHKERRQ(ierr);
1044 if (flag) {
1045 jac->symt = indx;
1046 PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
1047 }
1048
1049 ierr = PetscOptionsTail();CHKERRQ(ierr);
1050 PetscFunctionReturn(0);
1051 }
1052
PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)1053 static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
1054 {
1055 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1056 PetscErrorCode ierr;
1057 PetscBool iascii;
1058 const char *symt = 0;
1059
1060 PetscFunctionBegin;
1061 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
1062 if (iascii) {
1063 ierr = PetscViewerASCIIPrintf(viewer," HYPRE ParaSails preconditioning\n");CHKERRQ(ierr);
1064 ierr = PetscViewerASCIIPrintf(viewer," nlevels %d\n",jac->nlevels);CHKERRQ(ierr);
1065 ierr = PetscViewerASCIIPrintf(viewer," threshold %g\n",(double)jac->threshold);CHKERRQ(ierr);
1066 ierr = PetscViewerASCIIPrintf(viewer," filter %g\n",(double)jac->filter);CHKERRQ(ierr);
1067 ierr = PetscViewerASCIIPrintf(viewer," load balance %g\n",(double)jac->loadbal);CHKERRQ(ierr);
1068 ierr = PetscViewerASCIIPrintf(viewer," reuse nonzero structure %s\n",PetscBools[jac->ruse]);CHKERRQ(ierr);
1069 ierr = PetscViewerASCIIPrintf(viewer," print info to screen %s\n",PetscBools[jac->logging]);CHKERRQ(ierr);
1070 if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
1071 else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
1072 else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
1073 else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
1074 ierr = PetscViewerASCIIPrintf(viewer," %s\n",symt);CHKERRQ(ierr);
1075 }
1076 PetscFunctionReturn(0);
1077 }
1078 /* --------------------------------------------------------------------------------------------*/
PCSetFromOptions_HYPRE_AMS(PetscOptionItems * PetscOptionsObject,PC pc)1079 static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PetscOptionItems *PetscOptionsObject,PC pc)
1080 {
1081 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1082 PetscErrorCode ierr;
1083 PetscInt n;
1084 PetscBool flag,flag2,flag3,flag4;
1085
1086 PetscFunctionBegin;
1087 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE AMS Options");CHKERRQ(ierr);
1088 ierr = PetscOptionsInt("-pc_hypre_ams_print_level","Debugging output level for AMS","None",jac->as_print,&jac->as_print,&flag);CHKERRQ(ierr);
1089 if (flag) PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
1090 ierr = PetscOptionsInt("-pc_hypre_ams_max_iter","Maximum number of AMS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);CHKERRQ(ierr);
1091 if (flag) PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1092 ierr = PetscOptionsInt("-pc_hypre_ams_cycle_type","Cycle type for AMS multigrid","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag);CHKERRQ(ierr);
1093 if (flag) PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1094 ierr = PetscOptionsReal("-pc_hypre_ams_tol","Error tolerance for AMS multigrid","None",jac->as_tol,&jac->as_tol,&flag);CHKERRQ(ierr);
1095 if (flag) PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
1096 ierr = PetscOptionsInt("-pc_hypre_ams_relax_type","Relaxation type for AMS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);CHKERRQ(ierr);
1097 ierr = PetscOptionsInt("-pc_hypre_ams_relax_times","Number of relaxation steps for AMS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);CHKERRQ(ierr);
1098 ierr = PetscOptionsReal("-pc_hypre_ams_relax_weight","Relaxation weight for AMS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);CHKERRQ(ierr);
1099 ierr = PetscOptionsReal("-pc_hypre_ams_omega","SSOR coefficient for AMS smoother","None",jac->as_omega,&jac->as_omega,&flag4);CHKERRQ(ierr);
1100 if (flag || flag2 || flag3 || flag4) {
1101 PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1102 jac->as_relax_times,
1103 jac->as_relax_weight,
1104 jac->as_omega));
1105 }
1106 ierr = 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);CHKERRQ(ierr);
1107 n = 5;
1108 ierr = PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options","AMG options for vector Poisson","None",jac->as_amg_alpha_opts,&n,&flag2);CHKERRQ(ierr);
1109 if (flag || flag2) {
1110 PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0], /* AMG coarsen type */
1111 jac->as_amg_alpha_opts[1], /* AMG agg_levels */
1112 jac->as_amg_alpha_opts[2], /* AMG relax_type */
1113 jac->as_amg_alpha_theta,
1114 jac->as_amg_alpha_opts[3], /* AMG interp_type */
1115 jac->as_amg_alpha_opts[4])); /* AMG Pmax */
1116 }
1117 ierr = 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);CHKERRQ(ierr);
1118 n = 5;
1119 ierr = PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options","AMG options for scalar Poisson solver","None",jac->as_amg_beta_opts,&n,&flag2);CHKERRQ(ierr);
1120 if (flag || flag2) {
1121 PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */
1122 jac->as_amg_beta_opts[1], /* AMG agg_levels */
1123 jac->as_amg_beta_opts[2], /* AMG relax_type */
1124 jac->as_amg_beta_theta,
1125 jac->as_amg_beta_opts[3], /* AMG interp_type */
1126 jac->as_amg_beta_opts[4])); /* AMG Pmax */
1127 }
1128 ierr = 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);CHKERRQ(ierr);
1129 if (flag) { /* override HYPRE's default only if the options is used */
1130 PetscStackCallStandard(HYPRE_AMSSetProjectionFrequency,(jac->hsolver,jac->ams_proj_freq));
1131 }
1132 ierr = PetscOptionsTail();CHKERRQ(ierr);
1133 PetscFunctionReturn(0);
1134 }
1135
PCView_HYPRE_AMS(PC pc,PetscViewer viewer)1136 static PetscErrorCode PCView_HYPRE_AMS(PC pc,PetscViewer viewer)
1137 {
1138 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1139 PetscErrorCode ierr;
1140 PetscBool iascii;
1141
1142 PetscFunctionBegin;
1143 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
1144 if (iascii) {
1145 ierr = PetscViewerASCIIPrintf(viewer," HYPRE AMS preconditioning\n");CHKERRQ(ierr);
1146 ierr = PetscViewerASCIIPrintf(viewer," subspace iterations per application %d\n",jac->as_max_iter);CHKERRQ(ierr);
1147 ierr = PetscViewerASCIIPrintf(viewer," subspace cycle type %d\n",jac->ams_cycle_type);CHKERRQ(ierr);
1148 ierr = PetscViewerASCIIPrintf(viewer," subspace iteration tolerance %g\n",jac->as_tol);CHKERRQ(ierr);
1149 ierr = PetscViewerASCIIPrintf(viewer," smoother type %d\n",jac->as_relax_type);CHKERRQ(ierr);
1150 ierr = PetscViewerASCIIPrintf(viewer," number of smoothing steps %d\n",jac->as_relax_times);CHKERRQ(ierr);
1151 ierr = PetscViewerASCIIPrintf(viewer," smoother weight %g\n",jac->as_relax_weight);CHKERRQ(ierr);
1152 ierr = PetscViewerASCIIPrintf(viewer," smoother omega %g\n",jac->as_omega);CHKERRQ(ierr);
1153 if (jac->alpha_Poisson) {
1154 ierr = PetscViewerASCIIPrintf(viewer," vector Poisson solver (passed in by user)\n");CHKERRQ(ierr);
1155 } else {
1156 ierr = PetscViewerASCIIPrintf(viewer," vector Poisson solver (computed) \n");CHKERRQ(ierr);
1157 }
1158 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);CHKERRQ(ierr);
1159 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);CHKERRQ(ierr);
1160 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);CHKERRQ(ierr);
1161 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);CHKERRQ(ierr);
1162 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);CHKERRQ(ierr);
1163 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);CHKERRQ(ierr);
1164 if (!jac->ams_beta_is_zero) {
1165 if (jac->beta_Poisson) {
1166 ierr = PetscViewerASCIIPrintf(viewer," scalar Poisson solver (passed in by user)\n");CHKERRQ(ierr);
1167 } else {
1168 ierr = PetscViewerASCIIPrintf(viewer," scalar Poisson solver (computed) \n");CHKERRQ(ierr);
1169 }
1170 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);CHKERRQ(ierr);
1171 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);CHKERRQ(ierr);
1172 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);CHKERRQ(ierr);
1173 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);CHKERRQ(ierr);
1174 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);CHKERRQ(ierr);
1175 ierr = PetscViewerASCIIPrintf(viewer," boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);CHKERRQ(ierr);
1176 if (jac->ams_beta_is_zero_part) {
1177 ierr = PetscViewerASCIIPrintf(viewer," compatible subspace projection frequency %d (-1 HYPRE uses default)\n",jac->ams_proj_freq);CHKERRQ(ierr);
1178 }
1179 } else {
1180 ierr = PetscViewerASCIIPrintf(viewer," scalar Poisson solver not used (zero-conductivity everywhere) \n");CHKERRQ(ierr);
1181 }
1182 }
1183 PetscFunctionReturn(0);
1184 }
1185
PCSetFromOptions_HYPRE_ADS(PetscOptionItems * PetscOptionsObject,PC pc)1186 static PetscErrorCode PCSetFromOptions_HYPRE_ADS(PetscOptionItems *PetscOptionsObject,PC pc)
1187 {
1188 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1189 PetscErrorCode ierr;
1190 PetscInt n;
1191 PetscBool flag,flag2,flag3,flag4;
1192
1193 PetscFunctionBegin;
1194 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE ADS Options");CHKERRQ(ierr);
1195 ierr = PetscOptionsInt("-pc_hypre_ads_print_level","Debugging output level for ADS","None",jac->as_print,&jac->as_print,&flag);CHKERRQ(ierr);
1196 if (flag) PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
1197 ierr = PetscOptionsInt("-pc_hypre_ads_max_iter","Maximum number of ADS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);CHKERRQ(ierr);
1198 if (flag) PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1199 ierr = PetscOptionsInt("-pc_hypre_ads_cycle_type","Cycle type for ADS multigrid","None",jac->ads_cycle_type,&jac->ads_cycle_type,&flag);CHKERRQ(ierr);
1200 if (flag) PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ads_cycle_type));
1201 ierr = PetscOptionsReal("-pc_hypre_ads_tol","Error tolerance for ADS multigrid","None",jac->as_tol,&jac->as_tol,&flag);CHKERRQ(ierr);
1202 if (flag) PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
1203 ierr = PetscOptionsInt("-pc_hypre_ads_relax_type","Relaxation type for ADS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);CHKERRQ(ierr);
1204 ierr = PetscOptionsInt("-pc_hypre_ads_relax_times","Number of relaxation steps for ADS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);CHKERRQ(ierr);
1205 ierr = PetscOptionsReal("-pc_hypre_ads_relax_weight","Relaxation weight for ADS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);CHKERRQ(ierr);
1206 ierr = PetscOptionsReal("-pc_hypre_ads_omega","SSOR coefficient for ADS smoother","None",jac->as_omega,&jac->as_omega,&flag4);CHKERRQ(ierr);
1207 if (flag || flag2 || flag3 || flag4) {
1208 PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1209 jac->as_relax_times,
1210 jac->as_relax_weight,
1211 jac->as_omega));
1212 }
1213 ierr = 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);CHKERRQ(ierr);
1214 n = 5;
1215 ierr = PetscOptionsIntArray("-pc_hypre_ads_ams_options","AMG options for AMS solver inside ADS","None",jac->as_amg_alpha_opts,&n,&flag2);CHKERRQ(ierr);
1216 ierr = PetscOptionsInt("-pc_hypre_ads_ams_cycle_type","Cycle type for AMS solver inside ADS","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag3);CHKERRQ(ierr);
1217 if (flag || flag2 || flag3) {
1218 PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type, /* AMS cycle type */
1219 jac->as_amg_alpha_opts[0], /* AMG coarsen type */
1220 jac->as_amg_alpha_opts[1], /* AMG agg_levels */
1221 jac->as_amg_alpha_opts[2], /* AMG relax_type */
1222 jac->as_amg_alpha_theta,
1223 jac->as_amg_alpha_opts[3], /* AMG interp_type */
1224 jac->as_amg_alpha_opts[4])); /* AMG Pmax */
1225 }
1226 ierr = 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);CHKERRQ(ierr);
1227 n = 5;
1228 ierr = PetscOptionsIntArray("-pc_hypre_ads_amg_options","AMG options for vector AMG solver inside ADS","None",jac->as_amg_beta_opts,&n,&flag2);CHKERRQ(ierr);
1229 if (flag || flag2) {
1230 PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */
1231 jac->as_amg_beta_opts[1], /* AMG agg_levels */
1232 jac->as_amg_beta_opts[2], /* AMG relax_type */
1233 jac->as_amg_beta_theta,
1234 jac->as_amg_beta_opts[3], /* AMG interp_type */
1235 jac->as_amg_beta_opts[4])); /* AMG Pmax */
1236 }
1237 ierr = PetscOptionsTail();CHKERRQ(ierr);
1238 PetscFunctionReturn(0);
1239 }
1240
PCView_HYPRE_ADS(PC pc,PetscViewer viewer)1241 static PetscErrorCode PCView_HYPRE_ADS(PC pc,PetscViewer viewer)
1242 {
1243 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1244 PetscErrorCode ierr;
1245 PetscBool iascii;
1246
1247 PetscFunctionBegin;
1248 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
1249 if (iascii) {
1250 ierr = PetscViewerASCIIPrintf(viewer," HYPRE ADS preconditioning\n");CHKERRQ(ierr);
1251 ierr = PetscViewerASCIIPrintf(viewer," subspace iterations per application %d\n",jac->as_max_iter);CHKERRQ(ierr);
1252 ierr = PetscViewerASCIIPrintf(viewer," subspace cycle type %d\n",jac->ads_cycle_type);CHKERRQ(ierr);
1253 ierr = PetscViewerASCIIPrintf(viewer," subspace iteration tolerance %g\n",jac->as_tol);CHKERRQ(ierr);
1254 ierr = PetscViewerASCIIPrintf(viewer," smoother type %d\n",jac->as_relax_type);CHKERRQ(ierr);
1255 ierr = PetscViewerASCIIPrintf(viewer," number of smoothing steps %d\n",jac->as_relax_times);CHKERRQ(ierr);
1256 ierr = PetscViewerASCIIPrintf(viewer," smoother weight %g\n",jac->as_relax_weight);CHKERRQ(ierr);
1257 ierr = PetscViewerASCIIPrintf(viewer," smoother omega %g\n",jac->as_omega);CHKERRQ(ierr);
1258 ierr = PetscViewerASCIIPrintf(viewer," AMS solver using boomerAMG\n");CHKERRQ(ierr);
1259 ierr = PetscViewerASCIIPrintf(viewer," subspace cycle type %d\n",jac->ams_cycle_type);CHKERRQ(ierr);
1260 ierr = PetscViewerASCIIPrintf(viewer," coarsening type %d\n",jac->as_amg_alpha_opts[0]);CHKERRQ(ierr);
1261 ierr = PetscViewerASCIIPrintf(viewer," levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);CHKERRQ(ierr);
1262 ierr = PetscViewerASCIIPrintf(viewer," relaxation type %d\n",jac->as_amg_alpha_opts[2]);CHKERRQ(ierr);
1263 ierr = PetscViewerASCIIPrintf(viewer," interpolation type %d\n",jac->as_amg_alpha_opts[3]);CHKERRQ(ierr);
1264 ierr = PetscViewerASCIIPrintf(viewer," max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);CHKERRQ(ierr);
1265 ierr = PetscViewerASCIIPrintf(viewer," strength threshold %g\n",jac->as_amg_alpha_theta);CHKERRQ(ierr);
1266 ierr = PetscViewerASCIIPrintf(viewer," vector Poisson solver using boomerAMG\n");CHKERRQ(ierr);
1267 ierr = PetscViewerASCIIPrintf(viewer," coarsening type %d\n",jac->as_amg_beta_opts[0]);CHKERRQ(ierr);
1268 ierr = PetscViewerASCIIPrintf(viewer," levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);CHKERRQ(ierr);
1269 ierr = PetscViewerASCIIPrintf(viewer," relaxation type %d\n",jac->as_amg_beta_opts[2]);CHKERRQ(ierr);
1270 ierr = PetscViewerASCIIPrintf(viewer," interpolation type %d\n",jac->as_amg_beta_opts[3]);CHKERRQ(ierr);
1271 ierr = PetscViewerASCIIPrintf(viewer," max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);CHKERRQ(ierr);
1272 ierr = PetscViewerASCIIPrintf(viewer," strength threshold %g\n",jac->as_amg_beta_theta);CHKERRQ(ierr);
1273 }
1274 PetscFunctionReturn(0);
1275 }
1276
PCHYPRESetDiscreteGradient_HYPRE(PC pc,Mat G)1277 static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE(PC pc, Mat G)
1278 {
1279 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1280 PetscBool ishypre;
1281 PetscErrorCode ierr;
1282
1283 PetscFunctionBegin;
1284 ierr = PetscObjectTypeCompare((PetscObject)G,MATHYPRE,&ishypre);CHKERRQ(ierr);
1285 if (ishypre) {
1286 ierr = PetscObjectReference((PetscObject)G);CHKERRQ(ierr);
1287 ierr = MatDestroy(&jac->G);CHKERRQ(ierr);
1288 jac->G = G;
1289 } else {
1290 ierr = MatDestroy(&jac->G);CHKERRQ(ierr);
1291 ierr = MatConvert(G,MATHYPRE,MAT_INITIAL_MATRIX,&jac->G);CHKERRQ(ierr);
1292 }
1293 PetscFunctionReturn(0);
1294 }
1295
1296 /*@
1297 PCHYPRESetDiscreteGradient - Set discrete gradient matrix
1298
1299 Collective on PC
1300
1301 Input Parameters:
1302 + pc - the preconditioning context
1303 - G - the discrete gradient
1304
1305 Level: intermediate
1306
1307 Notes:
1308 G should have as many rows as the number of edges and as many columns as the number of vertices in the mesh
1309 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
1310
1311 .seealso:
1312 @*/
PCHYPRESetDiscreteGradient(PC pc,Mat G)1313 PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G)
1314 {
1315 PetscErrorCode ierr;
1316
1317 PetscFunctionBegin;
1318 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1319 PetscValidHeaderSpecific(G,MAT_CLASSID,2);
1320 PetscCheckSameComm(pc,1,G,2);
1321 ierr = PetscTryMethod(pc,"PCHYPRESetDiscreteGradient_C",(PC,Mat),(pc,G));CHKERRQ(ierr);
1322 PetscFunctionReturn(0);
1323 }
1324
PCHYPRESetDiscreteCurl_HYPRE(PC pc,Mat C)1325 static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C)
1326 {
1327 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1328 PetscBool ishypre;
1329 PetscErrorCode ierr;
1330
1331 PetscFunctionBegin;
1332 ierr = PetscObjectTypeCompare((PetscObject)C,MATHYPRE,&ishypre);CHKERRQ(ierr);
1333 if (ishypre) {
1334 ierr = PetscObjectReference((PetscObject)C);CHKERRQ(ierr);
1335 ierr = MatDestroy(&jac->C);CHKERRQ(ierr);
1336 jac->C = C;
1337 } else {
1338 ierr = MatDestroy(&jac->C);CHKERRQ(ierr);
1339 ierr = MatConvert(C,MATHYPRE,MAT_INITIAL_MATRIX,&jac->C);CHKERRQ(ierr);
1340 }
1341 PetscFunctionReturn(0);
1342 }
1343
1344 /*@
1345 PCHYPRESetDiscreteCurl - Set discrete curl matrix
1346
1347 Collective on PC
1348
1349 Input Parameters:
1350 + pc - the preconditioning context
1351 - C - the discrete curl
1352
1353 Level: intermediate
1354
1355 Notes:
1356 C should have as many rows as the number of faces and as many columns as the number of edges in the mesh
1357 Each row of G has as many nonzeros as the number of edges of a face, with column indexes being the global indexes of the corresponding edge: matrix entries are +1 and -1 depending on edge orientation with respect to the face orientation
1358
1359 .seealso:
1360 @*/
PCHYPRESetDiscreteCurl(PC pc,Mat C)1361 PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C)
1362 {
1363 PetscErrorCode ierr;
1364
1365 PetscFunctionBegin;
1366 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1367 PetscValidHeaderSpecific(C,MAT_CLASSID,2);
1368 PetscCheckSameComm(pc,1,C,2);
1369 ierr = PetscTryMethod(pc,"PCHYPRESetDiscreteCurl_C",(PC,Mat),(pc,C));CHKERRQ(ierr);
1370 PetscFunctionReturn(0);
1371 }
1372
PCHYPRESetInterpolations_HYPRE(PC pc,PetscInt dim,Mat RT_PiFull,Mat RT_Pi[],Mat ND_PiFull,Mat ND_Pi[])1373 static PetscErrorCode PCHYPRESetInterpolations_HYPRE(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[])
1374 {
1375 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1376 PetscBool ishypre;
1377 PetscErrorCode ierr;
1378 PetscInt i;
1379 PetscFunctionBegin;
1380
1381 ierr = MatDestroy(&jac->RT_PiFull);CHKERRQ(ierr);
1382 ierr = MatDestroy(&jac->ND_PiFull);CHKERRQ(ierr);
1383 for (i=0;i<3;++i) {
1384 ierr = MatDestroy(&jac->RT_Pi[i]);CHKERRQ(ierr);
1385 ierr = MatDestroy(&jac->ND_Pi[i]);CHKERRQ(ierr);
1386 }
1387
1388 jac->dim = dim;
1389 if (RT_PiFull) {
1390 ierr = PetscObjectTypeCompare((PetscObject)RT_PiFull,MATHYPRE,&ishypre);CHKERRQ(ierr);
1391 if (ishypre) {
1392 ierr = PetscObjectReference((PetscObject)RT_PiFull);CHKERRQ(ierr);
1393 jac->RT_PiFull = RT_PiFull;
1394 } else {
1395 ierr = MatConvert(RT_PiFull,MATHYPRE,MAT_INITIAL_MATRIX,&jac->RT_PiFull);CHKERRQ(ierr);
1396 }
1397 }
1398 if (RT_Pi) {
1399 for (i=0;i<dim;++i) {
1400 if (RT_Pi[i]) {
1401 ierr = PetscObjectTypeCompare((PetscObject)RT_Pi[i],MATHYPRE,&ishypre);CHKERRQ(ierr);
1402 if (ishypre) {
1403 ierr = PetscObjectReference((PetscObject)RT_Pi[i]);CHKERRQ(ierr);
1404 jac->RT_Pi[i] = RT_Pi[i];
1405 } else {
1406 ierr = MatConvert(RT_Pi[i],MATHYPRE,MAT_INITIAL_MATRIX,&jac->RT_Pi[i]);CHKERRQ(ierr);
1407 }
1408 }
1409 }
1410 }
1411 if (ND_PiFull) {
1412 ierr = PetscObjectTypeCompare((PetscObject)ND_PiFull,MATHYPRE,&ishypre);CHKERRQ(ierr);
1413 if (ishypre) {
1414 ierr = PetscObjectReference((PetscObject)ND_PiFull);CHKERRQ(ierr);
1415 jac->ND_PiFull = ND_PiFull;
1416 } else {
1417 ierr = MatConvert(ND_PiFull,MATHYPRE,MAT_INITIAL_MATRIX,&jac->ND_PiFull);CHKERRQ(ierr);
1418 }
1419 }
1420 if (ND_Pi) {
1421 for (i=0;i<dim;++i) {
1422 if (ND_Pi[i]) {
1423 ierr = PetscObjectTypeCompare((PetscObject)ND_Pi[i],MATHYPRE,&ishypre);CHKERRQ(ierr);
1424 if (ishypre) {
1425 ierr = PetscObjectReference((PetscObject)ND_Pi[i]);CHKERRQ(ierr);
1426 jac->ND_Pi[i] = ND_Pi[i];
1427 } else {
1428 ierr = MatConvert(ND_Pi[i],MATHYPRE,MAT_INITIAL_MATRIX,&jac->ND_Pi[i]);CHKERRQ(ierr);
1429 }
1430 }
1431 }
1432 }
1433
1434 PetscFunctionReturn(0);
1435 }
1436
1437 /*@
1438 PCHYPRESetInterpolations - Set interpolation matrices for AMS/ADS preconditioner
1439
1440 Collective on PC
1441
1442 Input Parameters:
1443 + pc - the preconditioning context
1444 - dim - the dimension of the problem, only used in AMS
1445 - RT_PiFull - Raviart-Thomas interpolation matrix
1446 - RT_Pi - x/y/z component of Raviart-Thomas interpolation matrix
1447 - ND_PiFull - Nedelec interpolation matrix
1448 - ND_Pi - x/y/z component of Nedelec interpolation matrix
1449
1450 Notes:
1451 For AMS, only Nedelec interpolation matrices are needed, the Raviart-Thomas interpolation matrices can be set to NULL.
1452 For ADS, both type of interpolation matrices are needed.
1453 Level: intermediate
1454
1455 .seealso:
1456 @*/
PCHYPRESetInterpolations(PC pc,PetscInt dim,Mat RT_PiFull,Mat RT_Pi[],Mat ND_PiFull,Mat ND_Pi[])1457 PetscErrorCode PCHYPRESetInterpolations(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[])
1458 {
1459 PetscErrorCode ierr;
1460 PetscInt i;
1461
1462 PetscFunctionBegin;
1463 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1464 if (RT_PiFull) {
1465 PetscValidHeaderSpecific(RT_PiFull,MAT_CLASSID,3);
1466 PetscCheckSameComm(pc,1,RT_PiFull,3);
1467 }
1468 if (RT_Pi) {
1469 PetscValidPointer(RT_Pi,4);
1470 for (i=0;i<dim;++i) {
1471 if (RT_Pi[i]) {
1472 PetscValidHeaderSpecific(RT_Pi[i],MAT_CLASSID,4);
1473 PetscCheckSameComm(pc,1,RT_Pi[i],4);
1474 }
1475 }
1476 }
1477 if (ND_PiFull) {
1478 PetscValidHeaderSpecific(ND_PiFull,MAT_CLASSID,5);
1479 PetscCheckSameComm(pc,1,ND_PiFull,5);
1480 }
1481 if (ND_Pi) {
1482 PetscValidPointer(ND_Pi,6);
1483 for (i=0;i<dim;++i) {
1484 if (ND_Pi[i]) {
1485 PetscValidHeaderSpecific(ND_Pi[i],MAT_CLASSID,6);
1486 PetscCheckSameComm(pc,1,ND_Pi[i],6);
1487 }
1488 }
1489 }
1490 ierr = PetscTryMethod(pc,"PCHYPRESetInterpolations_C",(PC,PetscInt,Mat,Mat[],Mat,Mat[]),(pc,dim,RT_PiFull,RT_Pi,ND_PiFull,ND_Pi));CHKERRQ(ierr);
1491 PetscFunctionReturn(0);
1492 }
1493
PCHYPRESetPoissonMatrix_HYPRE(PC pc,Mat A,PetscBool isalpha)1494 static PetscErrorCode PCHYPRESetPoissonMatrix_HYPRE(PC pc, Mat A, PetscBool isalpha)
1495 {
1496 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1497 PetscBool ishypre;
1498 PetscErrorCode ierr;
1499
1500 PetscFunctionBegin;
1501 ierr = PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);CHKERRQ(ierr);
1502 if (ishypre) {
1503 if (isalpha) {
1504 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
1505 ierr = MatDestroy(&jac->alpha_Poisson);CHKERRQ(ierr);
1506 jac->alpha_Poisson = A;
1507 } else {
1508 if (A) {
1509 ierr = PetscObjectReference((PetscObject)A);CHKERRQ(ierr);
1510 } else {
1511 jac->ams_beta_is_zero = PETSC_TRUE;
1512 }
1513 ierr = MatDestroy(&jac->beta_Poisson);CHKERRQ(ierr);
1514 jac->beta_Poisson = A;
1515 }
1516 } else {
1517 if (isalpha) {
1518 ierr = MatDestroy(&jac->alpha_Poisson);CHKERRQ(ierr);
1519 ierr = MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&jac->alpha_Poisson);CHKERRQ(ierr);
1520 } else {
1521 if (A) {
1522 ierr = MatDestroy(&jac->beta_Poisson);CHKERRQ(ierr);
1523 ierr = MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&jac->beta_Poisson);CHKERRQ(ierr);
1524 } else {
1525 ierr = MatDestroy(&jac->beta_Poisson);CHKERRQ(ierr);
1526 jac->ams_beta_is_zero = PETSC_TRUE;
1527 }
1528 }
1529 }
1530 PetscFunctionReturn(0);
1531 }
1532
1533 /*@
1534 PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix
1535
1536 Collective on PC
1537
1538 Input Parameters:
1539 + pc - the preconditioning context
1540 - A - the matrix
1541
1542 Level: intermediate
1543
1544 Notes:
1545 A should be obtained by discretizing the vector valued Poisson problem with linear finite elements
1546
1547 .seealso:
1548 @*/
PCHYPRESetAlphaPoissonMatrix(PC pc,Mat A)1549 PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A)
1550 {
1551 PetscErrorCode ierr;
1552
1553 PetscFunctionBegin;
1554 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1555 PetscValidHeaderSpecific(A,MAT_CLASSID,2);
1556 PetscCheckSameComm(pc,1,A,2);
1557 ierr = PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_TRUE));CHKERRQ(ierr);
1558 PetscFunctionReturn(0);
1559 }
1560
1561 /*@
1562 PCHYPRESetBetaPoissonMatrix - Set Poisson matrix
1563
1564 Collective on PC
1565
1566 Input Parameters:
1567 + pc - the preconditioning context
1568 - A - the matrix
1569
1570 Level: intermediate
1571
1572 Notes:
1573 A should be obtained by discretizing the Poisson problem with linear finite elements.
1574 Following HYPRE convention, the scalar Poisson solver of AMS can be turned off by passing NULL.
1575
1576 .seealso:
1577 @*/
PCHYPRESetBetaPoissonMatrix(PC pc,Mat A)1578 PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A)
1579 {
1580 PetscErrorCode ierr;
1581
1582 PetscFunctionBegin;
1583 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1584 if (A) {
1585 PetscValidHeaderSpecific(A,MAT_CLASSID,2);
1586 PetscCheckSameComm(pc,1,A,2);
1587 }
1588 ierr = PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_FALSE));CHKERRQ(ierr);
1589 PetscFunctionReturn(0);
1590 }
1591
PCHYPRESetEdgeConstantVectors_HYPRE(PC pc,Vec ozz,Vec zoz,Vec zzo)1592 static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE(PC pc,Vec ozz, Vec zoz, Vec zzo)
1593 {
1594 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1595 PetscErrorCode ierr;
1596
1597 PetscFunctionBegin;
1598 /* throw away any vector if already set */
1599 if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0]));
1600 if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1]));
1601 if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2]));
1602 jac->constants[0] = NULL;
1603 jac->constants[1] = NULL;
1604 jac->constants[2] = NULL;
1605 ierr = VecHYPRE_IJVectorCreate(ozz,&jac->constants[0]);CHKERRQ(ierr);
1606 ierr = VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);CHKERRQ(ierr);
1607 ierr = VecHYPRE_IJVectorCreate(zoz,&jac->constants[1]);CHKERRQ(ierr);
1608 ierr = VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);CHKERRQ(ierr);
1609 jac->dim = 2;
1610 if (zzo) {
1611 ierr = VecHYPRE_IJVectorCreate(zzo,&jac->constants[2]);CHKERRQ(ierr);
1612 ierr = VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);CHKERRQ(ierr);
1613 jac->dim++;
1614 }
1615 PetscFunctionReturn(0);
1616 }
1617
1618 /*@
1619 PCHYPRESetEdgeConstantVectors - Set the representation of the constant vector fields in edge element basis
1620
1621 Collective on PC
1622
1623 Input Parameters:
1624 + pc - the preconditioning context
1625 - ozz - vector representing (1,0,0) (or (1,0) in 2D)
1626 - zoz - vector representing (0,1,0) (or (0,1) in 2D)
1627 - zzo - vector representing (0,0,1) (use NULL in 2D)
1628
1629 Level: intermediate
1630
1631 Notes:
1632
1633 .seealso:
1634 @*/
PCHYPRESetEdgeConstantVectors(PC pc,Vec ozz,Vec zoz,Vec zzo)1635 PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo)
1636 {
1637 PetscErrorCode ierr;
1638
1639 PetscFunctionBegin;
1640 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
1641 PetscValidHeaderSpecific(ozz,VEC_CLASSID,2);
1642 PetscValidHeaderSpecific(zoz,VEC_CLASSID,3);
1643 if (zzo) PetscValidHeaderSpecific(zzo,VEC_CLASSID,4);
1644 PetscCheckSameComm(pc,1,ozz,2);
1645 PetscCheckSameComm(pc,1,zoz,3);
1646 if (zzo) PetscCheckSameComm(pc,1,zzo,4);
1647 ierr = PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));CHKERRQ(ierr);
1648 PetscFunctionReturn(0);
1649 }
1650
PCSetCoordinates_HYPRE(PC pc,PetscInt dim,PetscInt nloc,PetscReal * coords)1651 static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1652 {
1653 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1654 Vec tv;
1655 PetscInt i;
1656 PetscErrorCode ierr;
1657
1658 PetscFunctionBegin;
1659 /* throw away any coordinate vector if already set */
1660 if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0]));
1661 if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1]));
1662 if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2]));
1663 jac->dim = dim;
1664
1665 /* compute IJ vector for coordinates */
1666 ierr = VecCreate(PetscObjectComm((PetscObject)pc),&tv);CHKERRQ(ierr);
1667 ierr = VecSetType(tv,VECSTANDARD);CHKERRQ(ierr);
1668 ierr = VecSetSizes(tv,nloc,PETSC_DECIDE);CHKERRQ(ierr);
1669 for (i=0;i<dim;i++) {
1670 PetscScalar *array;
1671 PetscInt j;
1672
1673 ierr = VecHYPRE_IJVectorCreate(tv,&jac->coords[i]);CHKERRQ(ierr);
1674 ierr = VecGetArray(tv,&array);CHKERRQ(ierr);
1675 for (j=0;j<nloc;j++) {
1676 array[j] = coords[j*dim+i];
1677 }
1678 PetscStackCallStandard(HYPRE_IJVectorSetValues,(jac->coords[i],nloc,NULL,(HYPRE_Complex*)array));
1679 PetscStackCallStandard(HYPRE_IJVectorAssemble,(jac->coords[i]));
1680 ierr = VecRestoreArray(tv,&array);CHKERRQ(ierr);
1681 }
1682 ierr = VecDestroy(&tv);CHKERRQ(ierr);
1683 PetscFunctionReturn(0);
1684 }
1685
1686 /* ---------------------------------------------------------------------------------*/
1687
PCHYPREGetType_HYPRE(PC pc,const char * name[])1688 static PetscErrorCode PCHYPREGetType_HYPRE(PC pc,const char *name[])
1689 {
1690 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1691
1692 PetscFunctionBegin;
1693 *name = jac->hypre_type;
1694 PetscFunctionReturn(0);
1695 }
1696
PCHYPRESetType_HYPRE(PC pc,const char name[])1697 static PetscErrorCode PCHYPRESetType_HYPRE(PC pc,const char name[])
1698 {
1699 PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1700 PetscErrorCode ierr;
1701 PetscBool flag;
1702
1703 PetscFunctionBegin;
1704 if (jac->hypre_type) {
1705 ierr = PetscStrcmp(jac->hypre_type,name,&flag);CHKERRQ(ierr);
1706 if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set");
1707 PetscFunctionReturn(0);
1708 } else {
1709 ierr = PetscStrallocpy(name, &jac->hypre_type);CHKERRQ(ierr);
1710 }
1711
1712 jac->maxiter = PETSC_DEFAULT;
1713 jac->tol = PETSC_DEFAULT;
1714 jac->printstatistics = PetscLogPrintInfo;
1715
1716 ierr = PetscStrcmp("pilut",jac->hypre_type,&flag);CHKERRQ(ierr);
1717 if (flag) {
1718 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));CHKERRQ(ierr);
1719 PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver));
1720 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
1721 pc->ops->view = PCView_HYPRE_Pilut;
1722 jac->destroy = HYPRE_ParCSRPilutDestroy;
1723 jac->setup = HYPRE_ParCSRPilutSetup;
1724 jac->solve = HYPRE_ParCSRPilutSolve;
1725 jac->factorrowsize = PETSC_DEFAULT;
1726 PetscFunctionReturn(0);
1727 }
1728 ierr = PetscStrcmp("euclid",jac->hypre_type,&flag);CHKERRQ(ierr);
1729 if (flag) {
1730 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));CHKERRQ(ierr);
1731 PetscStackCallStandard(HYPRE_EuclidCreate,(jac->comm_hypre,&jac->hsolver));
1732 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
1733 pc->ops->view = PCView_HYPRE_Euclid;
1734 jac->destroy = HYPRE_EuclidDestroy;
1735 jac->setup = HYPRE_EuclidSetup;
1736 jac->solve = HYPRE_EuclidSolve;
1737 jac->factorrowsize = PETSC_DEFAULT;
1738 jac->eu_level = PETSC_DEFAULT; /* default */
1739 PetscFunctionReturn(0);
1740 }
1741 ierr = PetscStrcmp("parasails",jac->hypre_type,&flag);CHKERRQ(ierr);
1742 if (flag) {
1743 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));CHKERRQ(ierr);
1744 PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver));
1745 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
1746 pc->ops->view = PCView_HYPRE_ParaSails;
1747 jac->destroy = HYPRE_ParaSailsDestroy;
1748 jac->setup = HYPRE_ParaSailsSetup;
1749 jac->solve = HYPRE_ParaSailsSolve;
1750 /* initialize */
1751 jac->nlevels = 1;
1752 jac->threshold = .1;
1753 jac->filter = .1;
1754 jac->loadbal = 0;
1755 if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE;
1756 else jac->logging = (int) PETSC_FALSE;
1757
1758 jac->ruse = (int) PETSC_FALSE;
1759 jac->symt = 0;
1760 PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshold,jac->nlevels));
1761 PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
1762 PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
1763 PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
1764 PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
1765 PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
1766 PetscFunctionReturn(0);
1767 }
1768 ierr = PetscStrcmp("boomeramg",jac->hypre_type,&flag);CHKERRQ(ierr);
1769 if (flag) {
1770 ierr = HYPRE_BoomerAMGCreate(&jac->hsolver);
1771 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_BoomerAMG;
1772 pc->ops->view = PCView_HYPRE_BoomerAMG;
1773 pc->ops->applytranspose = PCApplyTranspose_HYPRE_BoomerAMG;
1774 pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
1775 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",PCGetInterpolations_BoomerAMG);CHKERRQ(ierr);
1776 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",PCGetCoarseOperators_BoomerAMG);CHKERRQ(ierr);
1777 jac->destroy = HYPRE_BoomerAMGDestroy;
1778 jac->setup = HYPRE_BoomerAMGSetup;
1779 jac->solve = HYPRE_BoomerAMGSolve;
1780 jac->applyrichardson = PETSC_FALSE;
1781 /* these defaults match the hypre defaults */
1782 jac->cycletype = 1;
1783 jac->maxlevels = 25;
1784 jac->maxiter = 1;
1785 jac->tol = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
1786 jac->truncfactor = 0.0;
1787 jac->strongthreshold = .25;
1788 jac->maxrowsum = .9;
1789 jac->coarsentype = 6;
1790 jac->measuretype = 0;
1791 jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
1792 jac->smoothtype = -1; /* Not set by default */
1793 jac->smoothnumlevels = 25;
1794 jac->eu_level = 0;
1795 jac->eu_droptolerance = 0;
1796 jac->eu_bj = 0;
1797 jac->relaxtype[0] = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */
1798 jac->relaxtype[2] = 9; /*G.E. */
1799 jac->relaxweight = 1.0;
1800 jac->outerrelaxweight = 1.0;
1801 jac->relaxorder = 1;
1802 jac->interptype = 0;
1803 jac->agg_nl = 0;
1804 jac->pmax = 0;
1805 jac->truncfactor = 0.0;
1806 jac->agg_num_paths = 1;
1807
1808 jac->nodal_coarsening = 0;
1809 jac->nodal_coarsening_diag = 0;
1810 jac->vec_interp_variant = 0;
1811 jac->vec_interp_qmax = 0;
1812 jac->vec_interp_smooth = PETSC_FALSE;
1813 jac->interp_refine = 0;
1814 jac->nodal_relax = PETSC_FALSE;
1815 jac->nodal_relax_levels = 1;
1816 PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
1817 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
1818 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
1819 PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
1820 PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
1821 PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
1822 PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
1823 PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
1824 PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
1825 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
1826 PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
1827 PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
1828 PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
1829 PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
1830 PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0])); /*defaults coarse to 9*/
1831 PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */
1832 PetscFunctionReturn(0);
1833 }
1834 ierr = PetscStrcmp("ams",jac->hypre_type,&flag);CHKERRQ(ierr);
1835 if (flag) {
1836 ierr = HYPRE_AMSCreate(&jac->hsolver);
1837 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_AMS;
1838 pc->ops->view = PCView_HYPRE_AMS;
1839 jac->destroy = HYPRE_AMSDestroy;
1840 jac->setup = HYPRE_AMSSetup;
1841 jac->solve = HYPRE_AMSSolve;
1842 jac->coords[0] = NULL;
1843 jac->coords[1] = NULL;
1844 jac->coords[2] = NULL;
1845 /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */
1846 jac->as_print = 0;
1847 jac->as_max_iter = 1; /* used as a preconditioner */
1848 jac->as_tol = 0.; /* used as a preconditioner */
1849 jac->ams_cycle_type = 13;
1850 /* Smoothing options */
1851 jac->as_relax_type = 2;
1852 jac->as_relax_times = 1;
1853 jac->as_relax_weight = 1.0;
1854 jac->as_omega = 1.0;
1855 /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1856 jac->as_amg_alpha_opts[0] = 10;
1857 jac->as_amg_alpha_opts[1] = 1;
1858 jac->as_amg_alpha_opts[2] = 6;
1859 jac->as_amg_alpha_opts[3] = 6;
1860 jac->as_amg_alpha_opts[4] = 4;
1861 jac->as_amg_alpha_theta = 0.25;
1862 /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1863 jac->as_amg_beta_opts[0] = 10;
1864 jac->as_amg_beta_opts[1] = 1;
1865 jac->as_amg_beta_opts[2] = 6;
1866 jac->as_amg_beta_opts[3] = 6;
1867 jac->as_amg_beta_opts[4] = 4;
1868 jac->as_amg_beta_theta = 0.25;
1869 PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
1870 PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1871 PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1872 PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
1873 PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1874 jac->as_relax_times,
1875 jac->as_relax_weight,
1876 jac->as_omega));
1877 PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0], /* AMG coarsen type */
1878 jac->as_amg_alpha_opts[1], /* AMG agg_levels */
1879 jac->as_amg_alpha_opts[2], /* AMG relax_type */
1880 jac->as_amg_alpha_theta,
1881 jac->as_amg_alpha_opts[3], /* AMG interp_type */
1882 jac->as_amg_alpha_opts[4])); /* AMG Pmax */
1883 PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */
1884 jac->as_amg_beta_opts[1], /* AMG agg_levels */
1885 jac->as_amg_beta_opts[2], /* AMG relax_type */
1886 jac->as_amg_beta_theta,
1887 jac->as_amg_beta_opts[3], /* AMG interp_type */
1888 jac->as_amg_beta_opts[4])); /* AMG Pmax */
1889 /* Zero conductivity */
1890 jac->ams_beta_is_zero = PETSC_FALSE;
1891 jac->ams_beta_is_zero_part = PETSC_FALSE;
1892 PetscFunctionReturn(0);
1893 }
1894 ierr = PetscStrcmp("ads",jac->hypre_type,&flag);CHKERRQ(ierr);
1895 if (flag) {
1896 ierr = HYPRE_ADSCreate(&jac->hsolver);
1897 pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ADS;
1898 pc->ops->view = PCView_HYPRE_ADS;
1899 jac->destroy = HYPRE_ADSDestroy;
1900 jac->setup = HYPRE_ADSSetup;
1901 jac->solve = HYPRE_ADSSolve;
1902 jac->coords[0] = NULL;
1903 jac->coords[1] = NULL;
1904 jac->coords[2] = NULL;
1905 /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */
1906 jac->as_print = 0;
1907 jac->as_max_iter = 1; /* used as a preconditioner */
1908 jac->as_tol = 0.; /* used as a preconditioner */
1909 jac->ads_cycle_type = 13;
1910 /* Smoothing options */
1911 jac->as_relax_type = 2;
1912 jac->as_relax_times = 1;
1913 jac->as_relax_weight = 1.0;
1914 jac->as_omega = 1.0;
1915 /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */
1916 jac->ams_cycle_type = 14;
1917 jac->as_amg_alpha_opts[0] = 10;
1918 jac->as_amg_alpha_opts[1] = 1;
1919 jac->as_amg_alpha_opts[2] = 6;
1920 jac->as_amg_alpha_opts[3] = 6;
1921 jac->as_amg_alpha_opts[4] = 4;
1922 jac->as_amg_alpha_theta = 0.25;
1923 /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1924 jac->as_amg_beta_opts[0] = 10;
1925 jac->as_amg_beta_opts[1] = 1;
1926 jac->as_amg_beta_opts[2] = 6;
1927 jac->as_amg_beta_opts[3] = 6;
1928 jac->as_amg_beta_opts[4] = 4;
1929 jac->as_amg_beta_theta = 0.25;
1930 PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
1931 PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1932 PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1933 PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
1934 PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1935 jac->as_relax_times,
1936 jac->as_relax_weight,
1937 jac->as_omega));
1938 PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type, /* AMG coarsen type */
1939 jac->as_amg_alpha_opts[0], /* AMG coarsen type */
1940 jac->as_amg_alpha_opts[1], /* AMG agg_levels */
1941 jac->as_amg_alpha_opts[2], /* AMG relax_type */
1942 jac->as_amg_alpha_theta,
1943 jac->as_amg_alpha_opts[3], /* AMG interp_type */
1944 jac->as_amg_alpha_opts[4])); /* AMG Pmax */
1945 PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */
1946 jac->as_amg_beta_opts[1], /* AMG agg_levels */
1947 jac->as_amg_beta_opts[2], /* AMG relax_type */
1948 jac->as_amg_beta_theta,
1949 jac->as_amg_beta_opts[3], /* AMG interp_type */
1950 jac->as_amg_beta_opts[4])); /* AMG Pmax */
1951 PetscFunctionReturn(0);
1952 }
1953 ierr = PetscFree(jac->hypre_type);CHKERRQ(ierr);
1954
1955 jac->hypre_type = NULL;
1956 SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are euclid, pilut, parasails, boomeramg, ams",name);
1957 PetscFunctionReturn(0);
1958 }
1959
1960 /*
1961 It only gets here if the HYPRE type has not been set before the call to
1962 ...SetFromOptions() which actually is most of the time
1963 */
PCSetFromOptions_HYPRE(PetscOptionItems * PetscOptionsObject,PC pc)1964 PetscErrorCode PCSetFromOptions_HYPRE(PetscOptionItems *PetscOptionsObject,PC pc)
1965 {
1966 PetscErrorCode ierr;
1967 PetscInt indx;
1968 const char *type[] = {"euclid","pilut","parasails","boomeramg","ams","ads"};
1969 PetscBool flg;
1970
1971 PetscFunctionBegin;
1972 ierr = PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");CHKERRQ(ierr);
1973 ierr = PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,ALEN(type),"boomeramg",&indx,&flg);CHKERRQ(ierr);
1974 if (flg) {
1975 ierr = PCHYPRESetType_HYPRE(pc,type[indx]);CHKERRQ(ierr);
1976 } else {
1977 ierr = PCHYPRESetType_HYPRE(pc,"boomeramg");CHKERRQ(ierr);
1978 }
1979 if (pc->ops->setfromoptions) {
1980 ierr = pc->ops->setfromoptions(PetscOptionsObject,pc);CHKERRQ(ierr);
1981 }
1982 ierr = PetscOptionsTail();CHKERRQ(ierr);
1983 PetscFunctionReturn(0);
1984 }
1985
1986 /*@C
1987 PCHYPRESetType - Sets which hypre preconditioner you wish to use
1988
1989 Input Parameters:
1990 + pc - the preconditioner context
1991 - name - either euclid, pilut, parasails, boomeramg, ams, ads
1992
1993 Options Database Keys:
1994 -pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads
1995
1996 Level: intermediate
1997
1998 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
1999 PCHYPRE
2000
2001 @*/
PCHYPRESetType(PC pc,const char name[])2002 PetscErrorCode PCHYPRESetType(PC pc,const char name[])
2003 {
2004 PetscErrorCode ierr;
2005
2006 PetscFunctionBegin;
2007 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2008 PetscValidCharPointer(name,2);
2009 ierr = PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));CHKERRQ(ierr);
2010 PetscFunctionReturn(0);
2011 }
2012
2013 /*@C
2014 PCHYPREGetType - Gets which hypre preconditioner you are using
2015
2016 Input Parameter:
2017 . pc - the preconditioner context
2018
2019 Output Parameter:
2020 . name - either euclid, pilut, parasails, boomeramg, ams, ads
2021
2022 Level: intermediate
2023
2024 .seealso: PCCreate(), PCHYPRESetType(), PCType (for list of available types), PC,
2025 PCHYPRE
2026
2027 @*/
PCHYPREGetType(PC pc,const char * name[])2028 PetscErrorCode PCHYPREGetType(PC pc,const char *name[])
2029 {
2030 PetscErrorCode ierr;
2031
2032 PetscFunctionBegin;
2033 PetscValidHeaderSpecific(pc,PC_CLASSID,1);
2034 PetscValidPointer(name,2);
2035 ierr = PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));CHKERRQ(ierr);
2036 PetscFunctionReturn(0);
2037 }
2038
2039 /*MC
2040 PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre
2041
2042 Options Database Keys:
2043 + -pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads
2044 - Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX
2045 preconditioner
2046
2047 Level: intermediate
2048
2049 Notes:
2050 Apart from pc_hypre_type (for which there is PCHYPRESetType()),
2051 the many hypre options can ONLY be set via the options database (e.g. the command line
2052 or with PetscOptionsSetValue(), there are no functions to set them)
2053
2054 The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_tol refer to the number of iterations
2055 (V-cycles) and tolerance that boomeramg does EACH time it is called. So for example, if
2056 -pc_hypre_boomeramg_max_iter is set to 2 then 2-V-cycles are being used to define the preconditioner
2057 (-pc_hypre_boomeramg_tol should be set to 0.0 - the default - to strictly use a fixed number of
2058 iterations per hypre call). -ksp_max_it and -ksp_rtol STILL determine the total number of iterations
2059 and tolerance for the Krylov solver. For example, if -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10
2060 then AT MOST twenty V-cycles of boomeramg will be called.
2061
2062 Note that the option -pc_hypre_boomeramg_relax_type_all defaults to symmetric relaxation
2063 (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry.
2064 Otherwise, you may want to use -pc_hypre_boomeramg_relax_type_all SOR/Jacobi.
2065 If you wish to use BoomerAMG WITHOUT a Krylov method use -ksp_type richardson NOT -ksp_type preonly
2066 and use -ksp_max_it to control the number of V-cycles.
2067 (see the PETSc FAQ.html at the PETSc website under the Documentation tab).
2068
2069 2007-02-03 Using HYPRE-1.11.1b, the routine HYPRE_BoomerAMGSolveT and the option
2070 -pc_hypre_parasails_reuse were failing with SIGSEGV. Dalcin L.
2071
2072 MatSetNearNullSpace() - if you provide a near null space to your matrix it is ignored by hypre UNLESS you also use
2073 the following two options:
2074 Options Database Keys:
2075 + -pc_hypre_boomeramg_nodal_coarsen <n> - where n is from 1 to 6 (see HYPRE_BOOMERAMGSetNodal())
2076 - -pc_hypre_boomeramg_vec_interp_variant <v> - where v is from 1 to 3 (see HYPRE_BoomerAMGSetInterpVecVariant())
2077
2078 Depending on the linear system you may see the same or different convergence depending on the values you use.
2079
2080 See PCPFMG for access to the hypre Struct PFMG solver
2081
2082 .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
2083 PCHYPRESetType(), PCPFMG
2084
2085 M*/
2086
PCCreate_HYPRE(PC pc)2087 PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
2088 {
2089 PC_HYPRE *jac;
2090 PetscErrorCode ierr;
2091
2092 PetscFunctionBegin;
2093 ierr = PetscNewLog(pc,&jac);CHKERRQ(ierr);
2094
2095 pc->data = jac;
2096 pc->ops->reset = PCReset_HYPRE;
2097 pc->ops->destroy = PCDestroy_HYPRE;
2098 pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
2099 pc->ops->setup = PCSetUp_HYPRE;
2100 pc->ops->apply = PCApply_HYPRE;
2101 jac->comm_hypre = MPI_COMM_NULL;
2102 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);CHKERRQ(ierr);
2103 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);CHKERRQ(ierr);
2104 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);CHKERRQ(ierr);
2105 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);CHKERRQ(ierr);
2106 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",PCHYPRESetDiscreteCurl_HYPRE);CHKERRQ(ierr);
2107 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",PCHYPRESetInterpolations_HYPRE);CHKERRQ(ierr);
2108 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE);CHKERRQ(ierr);
2109 ierr = PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",PCHYPRESetPoissonMatrix_HYPRE);CHKERRQ(ierr);
2110 PetscFunctionReturn(0);
2111 }
2112
2113 /* ---------------------------------------------------------------------------------------------------------------------------------*/
2114
2115 typedef struct {
2116 MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */
2117 HYPRE_StructSolver hsolver;
2118
2119 /* keep copy of PFMG options used so may view them */
2120 PetscInt its;
2121 double tol;
2122 PetscInt relax_type;
2123 PetscInt rap_type;
2124 PetscInt num_pre_relax,num_post_relax;
2125 PetscInt max_levels;
2126 } PC_PFMG;
2127
PCDestroy_PFMG(PC pc)2128 PetscErrorCode PCDestroy_PFMG(PC pc)
2129 {
2130 PetscErrorCode ierr;
2131 PC_PFMG *ex = (PC_PFMG*) pc->data;
2132
2133 PetscFunctionBegin;
2134 if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
2135 ierr = MPI_Comm_free(&ex->hcomm);CHKERRQ(ierr);
2136 ierr = PetscFree(pc->data);CHKERRQ(ierr);
2137 PetscFunctionReturn(0);
2138 }
2139
2140 static const char *PFMGRelaxType[] = {"Jacobi","Weighted-Jacobi","symmetric-Red/Black-Gauss-Seidel","Red/Black-Gauss-Seidel"};
2141 static const char *PFMGRAPType[] = {"Galerkin","non-Galerkin"};
2142
PCView_PFMG(PC pc,PetscViewer viewer)2143 PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
2144 {
2145 PetscErrorCode ierr;
2146 PetscBool iascii;
2147 PC_PFMG *ex = (PC_PFMG*) pc->data;
2148
2149 PetscFunctionBegin;
2150 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
2151 if (iascii) {
2152 ierr = PetscViewerASCIIPrintf(viewer," HYPRE PFMG preconditioning\n");CHKERRQ(ierr);
2153 ierr = PetscViewerASCIIPrintf(viewer," max iterations %d\n",ex->its);CHKERRQ(ierr);
2154 ierr = PetscViewerASCIIPrintf(viewer," tolerance %g\n",ex->tol);CHKERRQ(ierr);
2155 ierr = PetscViewerASCIIPrintf(viewer," relax type %s\n",PFMGRelaxType[ex->relax_type]);CHKERRQ(ierr);
2156 ierr = PetscViewerASCIIPrintf(viewer," RAP type %s\n",PFMGRAPType[ex->rap_type]);CHKERRQ(ierr);
2157 ierr = PetscViewerASCIIPrintf(viewer," number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);CHKERRQ(ierr);
2158 ierr = PetscViewerASCIIPrintf(viewer," max levels %d\n",ex->max_levels);CHKERRQ(ierr);
2159 }
2160 PetscFunctionReturn(0);
2161 }
2162
PCSetFromOptions_PFMG(PetscOptionItems * PetscOptionsObject,PC pc)2163 PetscErrorCode PCSetFromOptions_PFMG(PetscOptionItems *PetscOptionsObject,PC pc)
2164 {
2165 PetscErrorCode ierr;
2166 PC_PFMG *ex = (PC_PFMG*) pc->data;
2167 PetscBool flg = PETSC_FALSE;
2168
2169 PetscFunctionBegin;
2170 ierr = PetscOptionsHead(PetscOptionsObject,"PFMG options");CHKERRQ(ierr);
2171 ierr = PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);CHKERRQ(ierr);
2172 if (flg) {
2173 PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,3));
2174 }
2175 ierr = PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);CHKERRQ(ierr);
2176 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its));
2177 ierr = PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);CHKERRQ(ierr);
2178 PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax));
2179 ierr = PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);CHKERRQ(ierr);
2180 PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax));
2181
2182 ierr = PetscOptionsInt("-pc_pfmg_max_levels","Max Levels for MG hierarchy","HYPRE_StructPFMGSetMaxLevels",ex->max_levels,&ex->max_levels,NULL);CHKERRQ(ierr);
2183 PetscStackCallStandard(HYPRE_StructPFMGSetMaxLevels,(ex->hsolver,ex->max_levels));
2184
2185 ierr = PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);CHKERRQ(ierr);
2186 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol));
2187 ierr = PetscOptionsEList("-pc_pfmg_relax_type","Relax type for the up and down cycles","HYPRE_StructPFMGSetRelaxType",PFMGRelaxType,ALEN(PFMGRelaxType),PFMGRelaxType[ex->relax_type],&ex->relax_type,NULL);CHKERRQ(ierr);
2188 PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type));
2189 ierr = PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);CHKERRQ(ierr);
2190 PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type));
2191 ierr = PetscOptionsTail();CHKERRQ(ierr);
2192 PetscFunctionReturn(0);
2193 }
2194
PCApply_PFMG(PC pc,Vec x,Vec y)2195 PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
2196 {
2197 PetscErrorCode ierr;
2198 PC_PFMG *ex = (PC_PFMG*) pc->data;
2199 PetscScalar *yy;
2200 const PetscScalar *xx;
2201 PetscInt ilower[3],iupper[3];
2202 HYPRE_Int hlower[3],hupper[3];
2203 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
2204
2205 PetscFunctionBegin;
2206 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr);
2207 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr);
2208 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2209 iupper[0] += ilower[0] - 1;
2210 iupper[1] += ilower[1] - 1;
2211 iupper[2] += ilower[2] - 1;
2212 hlower[0] = (HYPRE_Int)ilower[0];
2213 hlower[1] = (HYPRE_Int)ilower[1];
2214 hlower[2] = (HYPRE_Int)ilower[2];
2215 hupper[0] = (HYPRE_Int)iupper[0];
2216 hupper[1] = (HYPRE_Int)iupper[1];
2217 hupper[2] = (HYPRE_Int)iupper[2];
2218
2219 /* copy x values over to hypre */
2220 PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
2221 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
2222 PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,hlower,hupper,(HYPRE_Complex*)xx));
2223 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
2224 PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
2225 PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
2226
2227 /* copy solution values back to PETSc */
2228 ierr = VecGetArray(y,&yy);CHKERRQ(ierr);
2229 PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,hlower,hupper,(HYPRE_Complex*)yy));
2230 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr);
2231 PetscFunctionReturn(0);
2232 }
2233
PCApplyRichardson_PFMG(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt * outits,PCRichardsonConvergedReason * reason)2234 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)
2235 {
2236 PC_PFMG *jac = (PC_PFMG*)pc->data;
2237 PetscErrorCode ierr;
2238 HYPRE_Int oits;
2239
2240 PetscFunctionBegin;
2241 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr);
2242 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its));
2243 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol));
2244
2245 ierr = PCApply_PFMG(pc,b,y);CHKERRQ(ierr);
2246 PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,&oits));
2247 *outits = oits;
2248 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2249 else *reason = PCRICHARDSON_CONVERGED_RTOL;
2250 PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol));
2251 PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its));
2252 PetscFunctionReturn(0);
2253 }
2254
2255
PCSetUp_PFMG(PC pc)2256 PetscErrorCode PCSetUp_PFMG(PC pc)
2257 {
2258 PetscErrorCode ierr;
2259 PC_PFMG *ex = (PC_PFMG*) pc->data;
2260 Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
2261 PetscBool flg;
2262
2263 PetscFunctionBegin;
2264 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESTRUCT,&flg);CHKERRQ(ierr);
2265 if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESTRUCT with this preconditioner");
2266
2267 /* create the hypre solver object and set its information */
2268 if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
2269 PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
2270 PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
2271 PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver));
2272 PetscFunctionReturn(0);
2273 }
2274
2275 /*MC
2276 PCPFMG - the hypre PFMG multigrid solver
2277
2278 Level: advanced
2279
2280 Options Database:
2281 + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner
2282 . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2283 . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2284 . -pc_pfmg_tol <tol> tolerance of PFMG
2285 . -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
2286 - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin
2287
2288 Notes:
2289 This is for CELL-centered descretizations
2290
2291 This must be used with the MATHYPRESTRUCT matrix type.
2292 This is less general than in hypre, it supports only one block per process defined by a PETSc DMDA.
2293
2294 .seealso: PCMG, MATHYPRESTRUCT
2295 M*/
2296
PCCreate_PFMG(PC pc)2297 PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
2298 {
2299 PetscErrorCode ierr;
2300 PC_PFMG *ex;
2301
2302 PetscFunctionBegin;
2303 ierr = PetscNew(&ex);CHKERRQ(ierr); \
2304 pc->data = ex;
2305
2306 ex->its = 1;
2307 ex->tol = 1.e-8;
2308 ex->relax_type = 1;
2309 ex->rap_type = 0;
2310 ex->num_pre_relax = 1;
2311 ex->num_post_relax = 1;
2312 ex->max_levels = 0;
2313
2314 pc->ops->setfromoptions = PCSetFromOptions_PFMG;
2315 pc->ops->view = PCView_PFMG;
2316 pc->ops->destroy = PCDestroy_PFMG;
2317 pc->ops->apply = PCApply_PFMG;
2318 pc->ops->applyrichardson = PCApplyRichardson_PFMG;
2319 pc->ops->setup = PCSetUp_PFMG;
2320
2321 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));CHKERRQ(ierr);
2322 PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
2323 PetscFunctionReturn(0);
2324 }
2325
2326 /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/
2327
2328 /* we know we are working with a HYPRE_SStructMatrix */
2329 typedef struct {
2330 MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
2331 HYPRE_SStructSolver ss_solver;
2332
2333 /* keep copy of SYSPFMG options used so may view them */
2334 PetscInt its;
2335 double tol;
2336 PetscInt relax_type;
2337 PetscInt num_pre_relax,num_post_relax;
2338 } PC_SysPFMG;
2339
PCDestroy_SysPFMG(PC pc)2340 PetscErrorCode PCDestroy_SysPFMG(PC pc)
2341 {
2342 PetscErrorCode ierr;
2343 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
2344
2345 PetscFunctionBegin;
2346 if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2347 ierr = MPI_Comm_free(&ex->hcomm);CHKERRQ(ierr);
2348 ierr = PetscFree(pc->data);CHKERRQ(ierr);
2349 PetscFunctionReturn(0);
2350 }
2351
2352 static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi","Red/Black-Gauss-Seidel"};
2353
PCView_SysPFMG(PC pc,PetscViewer viewer)2354 PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
2355 {
2356 PetscErrorCode ierr;
2357 PetscBool iascii;
2358 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
2359
2360 PetscFunctionBegin;
2361 ierr = PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);CHKERRQ(ierr);
2362 if (iascii) {
2363 ierr = PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG preconditioning\n");CHKERRQ(ierr);
2364 ierr = PetscViewerASCIIPrintf(viewer," max iterations %d\n",ex->its);CHKERRQ(ierr);
2365 ierr = PetscViewerASCIIPrintf(viewer," tolerance %g\n",ex->tol);CHKERRQ(ierr);
2366 ierr = PetscViewerASCIIPrintf(viewer," relax type %s\n",PFMGRelaxType[ex->relax_type]);CHKERRQ(ierr);
2367 ierr = PetscViewerASCIIPrintf(viewer," number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);CHKERRQ(ierr);
2368 }
2369 PetscFunctionReturn(0);
2370 }
2371
PCSetFromOptions_SysPFMG(PetscOptionItems * PetscOptionsObject,PC pc)2372 PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptionItems *PetscOptionsObject,PC pc)
2373 {
2374 PetscErrorCode ierr;
2375 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
2376 PetscBool flg = PETSC_FALSE;
2377
2378 PetscFunctionBegin;
2379 ierr = PetscOptionsHead(PetscOptionsObject,"SysPFMG options");CHKERRQ(ierr);
2380 ierr = PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);CHKERRQ(ierr);
2381 if (flg) {
2382 PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,3));
2383 }
2384 ierr = PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);CHKERRQ(ierr);
2385 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its));
2386 ierr = PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);CHKERRQ(ierr);
2387 PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax));
2388 ierr = PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);CHKERRQ(ierr);
2389 PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax));
2390
2391 ierr = PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);CHKERRQ(ierr);
2392 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol));
2393 ierr = PetscOptionsEList("-pc_syspfmg_relax_type","Relax type for the up and down cycles","HYPRE_SStructSysPFMGSetRelaxType",SysPFMGRelaxType,ALEN(SysPFMGRelaxType),SysPFMGRelaxType[ex->relax_type],&ex->relax_type,NULL);CHKERRQ(ierr);
2394 PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type));
2395 ierr = PetscOptionsTail();CHKERRQ(ierr);
2396 PetscFunctionReturn(0);
2397 }
2398
PCApply_SysPFMG(PC pc,Vec x,Vec y)2399 PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
2400 {
2401 PetscErrorCode ierr;
2402 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
2403 PetscScalar *yy;
2404 const PetscScalar *xx;
2405 PetscInt ilower[3],iupper[3];
2406 HYPRE_Int hlower[3],hupper[3];
2407 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
2408 PetscInt ordering= mx->dofs_order;
2409 PetscInt nvars = mx->nvars;
2410 PetscInt part = 0;
2411 PetscInt size;
2412 PetscInt i;
2413
2414 PetscFunctionBegin;
2415 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr);
2416 ierr = DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);CHKERRQ(ierr);
2417 /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2418 iupper[0] += ilower[0] - 1;
2419 iupper[1] += ilower[1] - 1;
2420 iupper[2] += ilower[2] - 1;
2421 hlower[0] = (HYPRE_Int)ilower[0];
2422 hlower[1] = (HYPRE_Int)ilower[1];
2423 hlower[2] = (HYPRE_Int)ilower[2];
2424 hupper[0] = (HYPRE_Int)iupper[0];
2425 hupper[1] = (HYPRE_Int)iupper[1];
2426 hupper[2] = (HYPRE_Int)iupper[2];
2427
2428 size = 1;
2429 for (i= 0; i< 3; i++) size *= (iupper[i]-ilower[i]+1);
2430
2431 /* copy x values over to hypre for variable ordering */
2432 if (ordering) {
2433 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2434 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
2435 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(xx+(size*i))));
2436 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
2437 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2438 PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
2439 PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2440
2441 /* copy solution values back to PETSc */
2442 ierr = VecGetArray(y,&yy);CHKERRQ(ierr);
2443 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(yy+(size*i))));
2444 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr);
2445 } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
2446 PetscScalar *z;
2447 PetscInt j, k;
2448
2449 ierr = PetscMalloc1(nvars*size,&z);CHKERRQ(ierr);
2450 PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2451 ierr = VecGetArrayRead(x,&xx);CHKERRQ(ierr);
2452
2453 /* transform nodal to hypre's variable ordering for sys_pfmg */
2454 for (i= 0; i< size; i++) {
2455 k= i*nvars;
2456 for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
2457 }
2458 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i))));
2459 ierr = VecRestoreArrayRead(x,&xx);CHKERRQ(ierr);
2460 PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2461 PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2462
2463 /* copy solution values back to PETSc */
2464 ierr = VecGetArray(y,&yy);CHKERRQ(ierr);
2465 for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i))));
2466 /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
2467 for (i= 0; i< size; i++) {
2468 k= i*nvars;
2469 for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
2470 }
2471 ierr = VecRestoreArray(y,&yy);CHKERRQ(ierr);
2472 ierr = PetscFree(z);CHKERRQ(ierr);
2473 }
2474 PetscFunctionReturn(0);
2475 }
2476
PCApplyRichardson_SysPFMG(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol,PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt * outits,PCRichardsonConvergedReason * reason)2477 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)
2478 {
2479 PC_SysPFMG *jac = (PC_SysPFMG*)pc->data;
2480 PetscErrorCode ierr;
2481 HYPRE_Int oits;
2482
2483 PetscFunctionBegin;
2484 ierr = PetscCitationsRegister(hypreCitation,&cite);CHKERRQ(ierr);
2485 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its));
2486 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol));
2487 ierr = PCApply_SysPFMG(pc,b,y);CHKERRQ(ierr);
2488 PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,&oits));
2489 *outits = oits;
2490 if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2491 else *reason = PCRICHARDSON_CONVERGED_RTOL;
2492 PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol));
2493 PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its));
2494 PetscFunctionReturn(0);
2495 }
2496
PCSetUp_SysPFMG(PC pc)2497 PetscErrorCode PCSetUp_SysPFMG(PC pc)
2498 {
2499 PetscErrorCode ierr;
2500 PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
2501 Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
2502 PetscBool flg;
2503
2504 PetscFunctionBegin;
2505 ierr = PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESSTRUCT,&flg);CHKERRQ(ierr);
2506 if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESSTRUCT with this preconditioner");
2507
2508 /* create the hypre sstruct solver object and set its information */
2509 if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2510 PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2511 PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver));
2512 PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2513 PetscFunctionReturn(0);
2514 }
2515
2516 /*MC
2517 PCSysPFMG - the hypre SysPFMG multigrid solver
2518
2519 Level: advanced
2520
2521 Options Database:
2522 + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner
2523 . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2524 . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2525 . -pc_syspfmg_tol <tol> tolerance of SysPFMG
2526 - -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel
2527
2528 Notes:
2529 This is for CELL-centered descretizations
2530
2531 This must be used with the MATHYPRESSTRUCT matrix type.
2532 This is less general than in hypre, it supports only one part, and one block per process defined by a PETSc DMDA.
2533 Also, only cell-centered variables.
2534
2535 .seealso: PCMG, MATHYPRESSTRUCT
2536 M*/
2537
PCCreate_SysPFMG(PC pc)2538 PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
2539 {
2540 PetscErrorCode ierr;
2541 PC_SysPFMG *ex;
2542
2543 PetscFunctionBegin;
2544 ierr = PetscNew(&ex);CHKERRQ(ierr); \
2545 pc->data = ex;
2546
2547 ex->its = 1;
2548 ex->tol = 1.e-8;
2549 ex->relax_type = 1;
2550 ex->num_pre_relax = 1;
2551 ex->num_post_relax = 1;
2552
2553 pc->ops->setfromoptions = PCSetFromOptions_SysPFMG;
2554 pc->ops->view = PCView_SysPFMG;
2555 pc->ops->destroy = PCDestroy_SysPFMG;
2556 pc->ops->apply = PCApply_SysPFMG;
2557 pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
2558 pc->ops->setup = PCSetUp_SysPFMG;
2559
2560 ierr = MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));CHKERRQ(ierr);
2561 PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2562 PetscFunctionReturn(0);
2563 }
2564