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