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