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