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