Actual source code: hypre.c
petsc-3.7.7 2017-09-25
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> /*I "petscpc.h" I*/
9: #include <../src/dm/impls/da/hypre/mhyp.h>
10: #include <_hypre_parcsr_ls.h>
12: static PetscBool cite = PETSC_FALSE;
13: 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";
15: /*
16: Private context (data structure) for the preconditioner.
17: */
18: typedef struct {
19: HYPRE_Solver hsolver;
20: HYPRE_IJMatrix ij;
21: HYPRE_IJVector b,x;
23: HYPRE_Int (*destroy)(HYPRE_Solver);
24: HYPRE_Int (*solve)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
25: HYPRE_Int (*setup)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
26: HYPRE_Int (*setdgrad)(HYPRE_Solver,HYPRE_ParCSRMatrix);
27: HYPRE_Int (*setdcurl)(HYPRE_Solver,HYPRE_ParCSRMatrix);
28: HYPRE_Int (*setcoord)(HYPRE_Solver,HYPRE_ParVector,HYPRE_ParVector,HYPRE_ParVector);
29: HYPRE_Int (*setdim)(HYPRE_Solver,HYPRE_Int);
31: MPI_Comm comm_hypre;
32: char *hypre_type;
34: /* options for Pilut and BoomerAMG*/
35: PetscInt maxiter;
36: double tol;
38: /* options for Pilut */
39: PetscInt factorrowsize;
41: /* options for ParaSails */
42: PetscInt nlevels;
43: double threshhold;
44: double filter;
45: PetscInt sym;
46: double loadbal;
47: PetscInt logging;
48: PetscInt ruse;
49: PetscInt symt;
51: /* options for BoomerAMG */
52: PetscBool printstatistics;
54: /* options for BoomerAMG */
55: PetscInt cycletype;
56: PetscInt maxlevels;
57: double strongthreshold;
58: double maxrowsum;
59: PetscInt gridsweeps[3];
60: PetscInt coarsentype;
61: PetscInt measuretype;
62: PetscInt smoothtype;
63: PetscInt smoothnumlevels;
64: PetscInt eu_level; /* Number of levels for ILU(k) in Euclid */
65: double eu_droptolerance; /* Drop tolerance for ILU(k) in Euclid */
66: PetscInt eu_bj; /* Defines use of Block Jacobi ILU in Euclid */
67: PetscInt relaxtype[3];
68: double relaxweight;
69: double outerrelaxweight;
70: PetscInt relaxorder;
71: double truncfactor;
72: PetscBool applyrichardson;
73: PetscInt pmax;
74: PetscInt interptype;
75: PetscInt agg_nl;
76: PetscInt agg_num_paths;
77: PetscInt nodal_coarsen;
78: PetscBool nodal_relax;
79: PetscInt nodal_relax_levels;
81: PetscInt nodal_coarsening;
82: PetscInt vec_interp_variant;
83: HYPRE_IJVector *hmnull;
84: HYPRE_ParVector *phmnull; /* near null space passed to hypre */
85: PetscInt n_hmnull;
86: Vec hmnull_constant;
87: 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 */
89: /* options for AS (Auxiliary Space preconditioners) */
90: PetscInt as_print;
91: PetscInt as_max_iter;
92: PetscReal as_tol;
93: PetscInt as_relax_type;
94: PetscInt as_relax_times;
95: PetscReal as_relax_weight;
96: PetscReal as_omega;
97: PetscInt as_amg_alpha_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for vector Poisson (AMS) or Curl problem (ADS) */
98: PetscReal as_amg_alpha_theta; /* AMG strength for vector Poisson (AMS) or Curl problem (ADS) */
99: PetscInt as_amg_beta_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for scalar Poisson (AMS) or vector Poisson (ADS) */
100: PetscReal as_amg_beta_theta; /* AMG strength for scalar Poisson (AMS) or vector Poisson (ADS) */
101: PetscInt ams_cycle_type;
102: PetscInt ads_cycle_type;
104: /* additional data */
105: HYPRE_IJVector coords[3];
106: HYPRE_IJVector constants[3];
107: HYPRE_IJMatrix G;
108: HYPRE_IJMatrix C;
109: HYPRE_IJMatrix alpha_Poisson;
110: HYPRE_IJMatrix beta_Poisson;
111: PetscBool ams_beta_is_zero;
112: } PC_HYPRE;
116: PetscErrorCode PCHYPREGetSolver(PC pc,HYPRE_Solver *hsolver)
117: {
118: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
121: *hsolver = jac->hsolver;
122: return(0);
123: }
125: /*
126: Replaces the address where the HYPRE vector points to its data with the address of
127: PETSc's data. Saves the old address so it can be reset when we are finished with it.
128: Allows use to get the data into a HYPRE vector without the cost of memcopies
129: */
130: #define HYPREReplacePointer(b,newvalue,savedvalue) { \
131: hypre_ParVector *par_vector = (hypre_ParVector*)hypre_IJVectorObject(((hypre_IJVector*)b)); \
132: hypre_Vector *local_vector = hypre_ParVectorLocalVector(par_vector); \
133: savedvalue = local_vector->data; \
134: local_vector->data = newvalue; \
135: }
139: static PetscErrorCode PCSetUp_HYPRE(PC pc)
140: {
141: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
142: PetscErrorCode ierr;
143: HYPRE_ParCSRMatrix hmat;
144: HYPRE_ParVector bv,xv;
145: PetscInt bs;
148: if (!jac->hypre_type) {
149: PCHYPRESetType(pc,"boomeramg");
150: }
152: if (pc->setupcalled) {
153: /* always destroy the old matrix and create a new memory;
154: hope this does not churn the memory too much. The problem
155: is I do not know if it is possible to put the matrix back to
156: its initial state so that we can directly copy the values
157: the second time through. */
158: PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij));
159: jac->ij = 0;
160: }
162: if (!jac->ij) { /* create the matrix the first time through */
163: MatHYPRE_IJMatrixCreate(pc->pmat,&jac->ij);
164: }
165: if (!jac->b) { /* create the vectors the first time through */
166: Vec x,b;
167: MatCreateVecs(pc->pmat,&x,&b);
168: VecHYPRE_IJVectorCreate(x,&jac->x);
169: VecHYPRE_IJVectorCreate(b,&jac->b);
170: VecDestroy(&x);
171: VecDestroy(&b);
172: }
174: /* special case for BoomerAMG */
175: if (jac->setup == HYPRE_BoomerAMGSetup) {
176: MatNullSpace mnull;
177: PetscBool has_const;
178: PetscInt nvec,i;
179: const Vec *vecs;
180: PetscScalar *petscvecarray;
181:
182: MatGetBlockSize(pc->pmat,&bs);
183: if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
184: MatGetNearNullSpace(pc->mat, &mnull);
185: if (mnull) {
186: MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);
187: PetscMalloc1(nvec+1,&jac->hmnull);
188: PetscMalloc1(nvec+1,&jac->hmnull_hypre_data_array);
189: PetscMalloc1(nvec+1,&jac->phmnull);
190: for (i=0; i<nvec; i++) {
191: VecHYPRE_IJVectorCreate(vecs[i],&jac->hmnull[i]);
192: VecGetArrayRead(vecs[i],(const PetscScalar **)&petscvecarray);
193: HYPREReplacePointer(jac->hmnull[i],petscvecarray,jac->hmnull_hypre_data_array[i]);
194: VecRestoreArrayRead(vecs[i],(const PetscScalar **)&petscvecarray);
195: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[i],(void**)&jac->phmnull[i]));
196: }
197: if (has_const) {
198: MatCreateVecs(pc->pmat,&jac->hmnull_constant,NULL);
199: VecSet(jac->hmnull_constant,1);
200: VecNormalize(jac->hmnull_constant,NULL);
201: VecHYPRE_IJVectorCreate(jac->hmnull_constant,&jac->hmnull[nvec]);
202: VecGetArrayRead(jac->hmnull_constant,(const PetscScalar **)&petscvecarray);
203: HYPREReplacePointer(jac->hmnull[nvec],petscvecarray,jac->hmnull_hypre_data_array[nvec]);
204: VecRestoreArrayRead(jac->hmnull_constant,(const PetscScalar **)&petscvecarray);
205: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[nvec],(void**)&jac->phmnull[nvec]));
206: nvec++;
207: }
208: PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVectors,(jac->hsolver,nvec,jac->phmnull));
209: jac->n_hmnull = nvec;
210: }
211: }
213: /* special case for AMS */
214: if (jac->setup == HYPRE_AMSSetup) {
215: if (!jac->coords[0] && !jac->constants[0]) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs either coordinate vectors via PCSetCoordinates() or edge constant vectors via PCHYPRESetEdgeConstantVectors()");
216: if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs discrete gradient operator via PCHYPRESetDiscreteGradient");
217: }
218: /* special case for ADS */
219: if (jac->setup == HYPRE_ADSSetup) {
220: if (!jac->coords[0]) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs coordinate vectors via PCSetCoordinates()");
221: 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");
222: if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs discrete gradient operator via PCHYPRESetDiscreteGradient");
223: if (!jac->C) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs discrete curl operator via PCHYPRESetDiscreteGradient");
224: }
225: MatHYPRE_IJMatrixCopy(pc->pmat,jac->ij);
226: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
227: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&bv));
228: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&xv));
229: PetscStackCall("HYPRE_SetupXXX",(*jac->setup)(jac->hsolver,hmat,bv,xv););
230: return(0);
231: }
235: static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
236: {
237: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
238: PetscErrorCode ierr;
239: HYPRE_ParCSRMatrix hmat;
240: PetscScalar *xv;
241: const PetscScalar *bv,*sbv;
242: HYPRE_ParVector jbv,jxv;
243: PetscScalar *sxv;
244: PetscInt hierr;
247: PetscCitationsRegister(hypreCitation,&cite);
248: if (!jac->applyrichardson) {VecSet(x,0.0);}
249: VecGetArrayRead(b,&bv);
250: VecGetArray(x,&xv);
251: HYPREReplacePointer(jac->b,(PetscScalar*)bv,sbv);
252: HYPREReplacePointer(jac->x,xv,sxv);
254: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
255: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
256: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));
257: PetscStackCall("Hypre solve",h(*jac->solve)(jac->hsolver,hmat,jbv,jxv);
258: if (hierr && hierr != HYPRE_ERROR_CONV) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
259: if (hierr) hypre__global_error = 0;);
261: HYPREReplacePointer(jac->b,(PetscScalar*)sbv,bv);
262: HYPREReplacePointer(jac->x,sxv,xv);
263: VecRestoreArray(x,&xv);
264: VecRestoreArrayRead(b,&bv);
265: return(0);
266: }
270: static PetscErrorCode PCReset_HYPRE(PC pc)
271: {
272: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
276: if (jac->ij) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij)); jac->ij = NULL;
277: if (jac->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->b)); jac->b = NULL;
278: if (jac->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->x)); jac->x = NULL;
279: if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0])); jac->coords[0] = NULL;
280: if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1])); jac->coords[1] = NULL;
281: if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2])); jac->coords[2] = NULL;
282: if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0])); jac->constants[0] = NULL;
283: if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1])); jac->constants[1] = NULL;
284: if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2])); jac->constants[2] = NULL;
285: if (jac->G) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->G)); jac->G = NULL;
286: if (jac->C) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->C)); jac->C = NULL;
287: if (jac->alpha_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->alpha_Poisson)); jac->alpha_Poisson = NULL;
288: if (jac->beta_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->beta_Poisson)); jac->beta_Poisson = NULL;
289: if (jac->n_hmnull && jac->hmnull) {
290: PetscInt i;
291: PETSC_UNUSED PetscScalar *petscvecarray;
293: for (i=0; i<jac->n_hmnull; i++) {
294: HYPREReplacePointer(jac->hmnull[i],jac->hmnull_hypre_data_array[i],petscvecarray);
295: PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->hmnull[i]));
296: }
297: PetscFree(jac->hmnull);
298: PetscFree(jac->hmnull_hypre_data_array);
299: PetscFree(jac->phmnull);
300: VecDestroy(&jac->hmnull_constant);
301: }
302: return(0);
303: }
307: static PetscErrorCode PCDestroy_HYPRE(PC pc)
308: {
309: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
310: PetscErrorCode ierr;
313: PCReset_HYPRE(pc);
314: if (jac->destroy) PetscStackCall("HYPRE_DestroyXXX",(*jac->destroy)(jac->hsolver););
315: PetscFree(jac->hypre_type);
316: if (jac->comm_hypre != MPI_COMM_NULL) { MPI_Comm_free(&(jac->comm_hypre));}
317: PetscFree(pc->data);
319: PetscObjectChangeTypeName((PetscObject)pc,0);
320: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);
321: PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);
322: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetCoordinates_C",NULL);
323: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",NULL);
324: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",NULL);
325: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetConstantEdgeVectors_C",NULL);
326: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetAlphaPoissonMatrix_C",NULL);
327: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetBetaPoissonMatrix_C",NULL);
328: return(0);
329: }
331: /* --------------------------------------------------------------------------------------------*/
334: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PetscOptionItems *PetscOptionsObject,PC pc)
335: {
336: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
338: PetscBool flag;
341: PetscOptionsHead(PetscOptionsObject,"HYPRE Pilut Options");
342: PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
343: if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter));
344: PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
345: if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol));
346: PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
347: if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize));
348: PetscOptionsTail();
349: return(0);
350: }
354: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
355: {
356: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
358: PetscBool iascii;
361: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
362: if (iascii) {
363: PetscViewerASCIIPrintf(viewer," HYPRE Pilut preconditioning\n");
364: if (jac->maxiter != PETSC_DEFAULT) {
365: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: maximum number of iterations %d\n",jac->maxiter);
366: } else {
367: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default maximum number of iterations \n");
368: }
369: if (jac->tol != PETSC_DEFAULT) {
370: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: drop tolerance %g\n",(double)jac->tol);
371: } else {
372: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default drop tolerance \n");
373: }
374: if (jac->factorrowsize != PETSC_DEFAULT) {
375: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: factor row size %d\n",jac->factorrowsize);
376: } else {
377: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default factor row size \n");
378: }
379: }
380: return(0);
381: }
383: /* --------------------------------------------------------------------------------------------*/
387: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
388: {
389: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
390: PetscErrorCode ierr;
391: HYPRE_ParCSRMatrix hmat;
392: PetscScalar *xv;
393: const PetscScalar *bv;
394: HYPRE_ParVector jbv,jxv;
395: PetscScalar *sbv,*sxv;
396: PetscInt hierr;
399: PetscCitationsRegister(hypreCitation,&cite);
400: VecSet(x,0.0);
401: VecGetArrayRead(b,&bv);
402: VecGetArray(x,&xv);
403: HYPREReplacePointer(jac->b,(PetscScalar*)bv,sbv);
404: HYPREReplacePointer(jac->x,xv,sxv);
406: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
407: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
408: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));
410: hHYPRE_BoomerAMGSolveT(jac->hsolver,hmat,jbv,jxv);
411: /* error code of 1 in BoomerAMG merely means convergence not achieved */
412: if (hierr && (hierr != 1)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
413: if (hierr) hypre__global_error = 0;
415: HYPREReplacePointer(jac->b,sbv,bv);
416: HYPREReplacePointer(jac->x,sxv,xv);
417: VecRestoreArray(x,&xv);
418: VecRestoreArrayRead(b,&bv);
419: return(0);
420: }
422: /* static array length */
423: #define ALEN(a) (sizeof(a)/sizeof((a)[0]))
425: static const char *HYPREBoomerAMGCycleType[] = {"","V","W"};
426: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
427: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
428: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
429: static const char *HYPREBoomerAMGSmoothType[] = {"Schwarz-smoothers","Pilut","ParaSails","Euclid"};
430: static const char *HYPREBoomerAMGRelaxType[] = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi",
431: "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi",
432: "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination",
433: "" /* 10 */, "" /* 11 */, "" /* 12 */, "" /* 13 */, "" /* 14 */,
434: "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"};
435: static const char *HYPREBoomerAMGInterpType[] = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
436: "ext+i-cc", "standard", "standard-wts", "", "", "FF", "FF1"};
439: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PetscOptionItems *PetscOptionsObject,PC pc)
440: {
441: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
443: PetscInt n,indx,level;
444: PetscBool flg, tmp_truth;
445: double tmpdbl, twodbl[2];
448: PetscOptionsHead(PetscOptionsObject,"HYPRE BoomerAMG Options");
449: PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);
450: if (flg) {
451: jac->cycletype = indx+1;
452: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
453: }
454: PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
455: if (flg) {
456: if (jac->maxlevels < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
457: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
458: }
459: PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
460: if (flg) {
461: if (jac->maxiter < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
462: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
463: }
464: PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);
465: if (flg) {
466: 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);
467: PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
468: }
470: PetscOptionsScalar("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);
471: if (flg) {
472: 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);
473: PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
474: }
476: PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);
477: if (flg) {
478: 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);
479: PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
480: }
482: PetscOptionsInt("-pc_hypre_boomeramg_agg_nl","Number of levels of aggressive coarsening","None",jac->agg_nl,&jac->agg_nl,&flg);
483: if (flg) {
484: 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);
486: PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
487: }
490: PetscOptionsInt("-pc_hypre_boomeramg_agg_num_paths","Number of paths for aggressive coarsening","None",jac->agg_num_paths,&jac->agg_num_paths,&flg);
491: if (flg) {
492: 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);
494: PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
495: }
498: PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
499: if (flg) {
500: 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);
501: PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
502: }
503: PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
504: if (flg) {
505: 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);
506: 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);
507: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
508: }
510: /* Grid sweeps */
511: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);
512: if (flg) {
513: PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx));
514: /* modify the jac structure so we can view the updated options with PC_View */
515: jac->gridsweeps[0] = indx;
516: jac->gridsweeps[1] = indx;
517: /*defaults coarse to 1 */
518: jac->gridsweeps[2] = 1;
519: }
521: PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen","Use a nodal based coarsening 1-6","HYPRE_BoomerAMGSetNodal",jac->nodal_coarsening, &jac->nodal_coarsening,&flg);
522: if (flg) {
523: PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,jac->nodal_coarsening));
524: }
526: PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_variant","Variant of algorithm 1-3","HYPRE_BoomerAMGSetInterpVecVariant",jac->vec_interp_variant, &jac->vec_interp_variant,&flg);
527: if (flg) {
528: PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecVariant,(jac->hsolver,jac->vec_interp_variant));
529: }
531: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);
532: if (flg) {
533: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1));
534: jac->gridsweeps[0] = indx;
535: }
536: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);
537: if (flg) {
538: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2));
539: jac->gridsweeps[1] = indx;
540: }
541: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);
542: if (flg) {
543: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3));
544: jac->gridsweeps[2] = indx;
545: }
547: /* Smooth type */
548: PetscOptionsEList("-pc_hypre_boomeramg_smooth_type","Enable more complex smoothers","None",HYPREBoomerAMGSmoothType,ALEN(HYPREBoomerAMGSmoothType),HYPREBoomerAMGSmoothType[0],&indx,&flg);
549: if (flg) {
550: jac->smoothtype = indx;
551: PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,indx+6));
552: jac->smoothnumlevels = 25;
553: PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,25));
554: }
556: /* Number of smoothing levels */
557: PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_levels","Number of levels on which more complex smoothers are used","None",25,&indx,&flg);
558: if (flg && (jac->smoothtype != -1)) {
559: jac->smoothnumlevels = indx;
560: PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,indx));
561: }
563: /* Number of levels for ILU(k) for Euclid */
564: PetscOptionsInt("-pc_hypre_boomeramg_eu_level","Number of levels for ILU(k) in Euclid smoother","None",0,&indx,&flg);
565: if (flg && (jac->smoothtype == 3)) {
566: jac->eu_level = indx;
567: PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,indx));
568: }
570: /* Filter for ILU(k) for Euclid */
571: double droptolerance;
572: PetscOptionsScalar("-pc_hypre_boomeramg_eu_droptolerance","Drop tolerance for ILU(k) in Euclid smoother","None",0,&droptolerance,&flg);
573: if (flg && (jac->smoothtype == 3)) {
574: jac->eu_droptolerance = droptolerance;
575: PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,droptolerance));
576: }
578: /* Use Block Jacobi ILUT for Euclid */
579: PetscOptionsBool("-pc_hypre_boomeramg_eu_bj", "Use Block Jacobi for ILU in Euclid smoother?", "None", PETSC_FALSE, &tmp_truth, &flg);
580: if (flg && (jac->smoothtype == 3)) {
581: jac->eu_bj = tmp_truth;
582: PetscStackCallStandard(HYPRE_BoomerAMGSetEuBJ,(jac->hsolver,jac->eu_bj));
583: }
585: /* Relax type */
586: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
587: if (flg) {
588: jac->relaxtype[0] = jac->relaxtype[1] = indx;
589: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx));
590: /* by default, coarse type set to 9 */
591: jac->relaxtype[2] = 9;
593: }
594: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
595: if (flg) {
596: jac->relaxtype[0] = indx;
597: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1));
598: }
599: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
600: if (flg) {
601: jac->relaxtype[1] = indx;
602: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2));
603: }
604: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);
605: if (flg) {
606: jac->relaxtype[2] = indx;
607: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3));
608: }
610: /* Relaxation Weight */
611: 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);
612: if (flg) {
613: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl));
614: jac->relaxweight = tmpdbl;
615: }
617: n = 2;
618: twodbl[0] = twodbl[1] = 1.0;
619: PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
620: if (flg) {
621: if (n == 2) {
622: indx = (int)PetscAbsReal(twodbl[1]);
623: PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx));
624: } 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);
625: }
627: /* Outer relaxation Weight */
628: 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);
629: if (flg) {
630: PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl));
631: jac->outerrelaxweight = tmpdbl;
632: }
634: n = 2;
635: twodbl[0] = twodbl[1] = 1.0;
636: PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
637: if (flg) {
638: if (n == 2) {
639: indx = (int)PetscAbsReal(twodbl[1]);
640: PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx));
641: } 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);
642: }
644: /* the Relax Order */
645: PetscOptionsBool("-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", PETSC_FALSE, &tmp_truth, &flg);
647: if (flg && tmp_truth) {
648: jac->relaxorder = 0;
649: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
650: }
651: PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);
652: if (flg) {
653: jac->measuretype = indx;
654: PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
655: }
656: /* update list length 3/07 */
657: PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);
658: if (flg) {
659: jac->coarsentype = indx;
660: PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
661: }
663: /* new 3/07 */
664: PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);
665: if (flg) {
666: jac->interptype = indx;
667: PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
668: }
670: PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
671: if (flg) {
672: level = 3;
673: PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);
675: jac->printstatistics = PETSC_TRUE;
676: PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level));
677: }
679: PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);
680: if (flg) {
681: level = 3;
682: PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);
684: jac->printstatistics = PETSC_TRUE;
685: PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level));
686: }
688: PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);
689: if (flg && tmp_truth) {
690: PetscInt tmp_int;
691: PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);
692: if (flg) jac->nodal_relax_levels = tmp_int;
693: PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6));
694: PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1));
695: PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0));
696: PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels));
697: }
699: PetscOptionsTail();
700: return(0);
701: }
705: 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)
706: {
707: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
709: PetscInt oits;
712: PetscCitationsRegister(hypreCitation,&cite);
713: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter));
714: PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol));
715: jac->applyrichardson = PETSC_TRUE;
716: PCApply_HYPRE(pc,b,y);
717: jac->applyrichardson = PETSC_FALSE;
718: PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
719: *outits = oits;
720: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
721: else *reason = PCRICHARDSON_CONVERGED_RTOL;
722: PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
723: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
724: return(0);
725: }
730: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
731: {
732: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
734: PetscBool iascii;
737: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
738: if (iascii) {
739: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG preconditioning\n");
740: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);
741: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum number of levels %d\n",jac->maxlevels);
742: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum number of iterations PER hypre call %d\n",jac->maxiter);
743: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Convergence tolerance PER hypre call %g\n",(double)jac->tol);
744: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Threshold for strong coupling %g\n",(double)jac->strongthreshold);
745: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Interpolation truncation factor %g\n",(double)jac->truncfactor);
746: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Interpolation: max elements per row %d\n",jac->pmax);
747: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Number of levels of aggressive coarsening %d\n",jac->agg_nl);
748: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Number of paths for aggressive coarsening %d\n",jac->agg_num_paths);
750: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum row sums %g\n",(double)jac->maxrowsum);
752: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps down %d\n",jac->gridsweeps[0]);
753: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps up %d\n",jac->gridsweeps[1]);
754: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps on coarse %d\n",jac->gridsweeps[2]);
756: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax down %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[0]]);
757: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax up %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[1]]);
758: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax on coarse %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[2]]);
760: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax weight (all) %g\n",(double)jac->relaxweight);
761: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Outer relax weight (all) %g\n",(double)jac->outerrelaxweight);
763: if (jac->relaxorder) {
764: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using CF-relaxation\n");
765: } else {
766: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Not using CF-relaxation\n");
767: }
768: if (jac->smoothtype!=-1) {
769: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Smooth type %s\n",HYPREBoomerAMGSmoothType[jac->smoothtype]);
770: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Smooth num levels %d\n",jac->smoothnumlevels);
771: } else {
772: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Not using more complex smoothers.\n");
773: }
774: if (jac->smoothtype==3) {
775: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Euclid ILU(k) levels %d\n",jac->eu_level);
776: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Euclid ILU(k) drop tolerance %g\n",jac->eu_droptolerance);
777: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Euclid ILU use Block-Jacobi? %d\n",jac->eu_bj);
778: }
779: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Measure type %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
780: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Coarsen type %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
781: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Interpolation type %s\n",HYPREBoomerAMGInterpType[jac->interptype]);
782: if (jac->nodal_coarsening) {
783: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal coarsening (with HYPRE_BOOMERAMGSetNodal() %D\n",jac->nodal_coarsening);
784: }
785: if (jac->vec_interp_variant) {
786: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: HYPRE_BoomerAMGSetInterpVecVariant() %D\n",jac->vec_interp_variant);
787: }
788: if (jac->nodal_relax) {
789: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal relaxation via Schwarz smoothing on levels %d\n",jac->nodal_relax_levels);
790: }
791: }
792: return(0);
793: }
795: /* --------------------------------------------------------------------------------------------*/
798: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PetscOptionItems *PetscOptionsObject,PC pc)
799: {
800: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
802: PetscInt indx;
803: PetscBool flag;
804: const char *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};
807: PetscOptionsHead(PetscOptionsObject,"HYPRE ParaSails Options");
808: PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);
809: PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshhold,&jac->threshhold,&flag);
810: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
812: PetscOptionsReal("-pc_hypre_parasails_filter","filter","None",jac->filter,&jac->filter,&flag);
813: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
815: PetscOptionsReal("-pc_hypre_parasails_loadbal","Load balance","None",jac->loadbal,&jac->loadbal,&flag);
816: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
818: PetscOptionsBool("-pc_hypre_parasails_logging","Print info to screen","None",(PetscBool)jac->logging,(PetscBool*)&jac->logging,&flag);
819: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
821: PetscOptionsBool("-pc_hypre_parasails_reuse","Reuse nonzero pattern in preconditioner","None",(PetscBool)jac->ruse,(PetscBool*)&jac->ruse,&flag);
822: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
824: PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);
825: if (flag) {
826: jac->symt = indx;
827: PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
828: }
830: PetscOptionsTail();
831: return(0);
832: }
836: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
837: {
838: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
840: PetscBool iascii;
841: const char *symt = 0;;
844: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
845: if (iascii) {
846: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails preconditioning\n");
847: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: nlevels %d\n",jac->nlevels);
848: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: threshold %g\n",(double)jac->threshhold);
849: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: filter %g\n",(double)jac->filter);
850: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: load balance %g\n",(double)jac->loadbal);
851: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: reuse nonzero structure %s\n",PetscBools[jac->ruse]);
852: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: print info to screen %s\n",PetscBools[jac->logging]);
853: if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
854: else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
855: else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
856: else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
857: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: %s\n",symt);
858: }
859: return(0);
860: }
861: /* --------------------------------------------------------------------------------------------*/
864: static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PetscOptionItems *PetscOptionsObject,PC pc)
865: {
866: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
868: PetscInt n;
869: PetscBool flag,flag2,flag3,flag4;
872: PetscOptionsHead(PetscOptionsObject,"HYPRE AMS Options");
873: PetscOptionsInt("-pc_hypre_ams_print_level","Debugging output level for AMS","None",jac->as_print,&jac->as_print,&flag);
874: if (flag) PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
875: PetscOptionsInt("-pc_hypre_ams_max_iter","Maximum number of AMS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
876: if (flag) PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
877: PetscOptionsInt("-pc_hypre_ams_cycle_type","Cycle type for AMS multigrid","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag);
878: if (flag) PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
879: PetscOptionsReal("-pc_hypre_ams_tol","Error tolerance for AMS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
880: if (flag) PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
881: PetscOptionsInt("-pc_hypre_ams_relax_type","Relaxation type for AMS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
882: PetscOptionsInt("-pc_hypre_ams_relax_times","Number of relaxation steps for AMS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
883: PetscOptionsReal("-pc_hypre_ams_relax_weight","Relaxation weight for AMS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
884: PetscOptionsReal("-pc_hypre_ams_omega","SSOR coefficient for AMS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
885: if (flag || flag2 || flag3 || flag4) {
886: PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
887: jac->as_relax_times,
888: jac->as_relax_weight,
889: jac->as_omega));
890: }
891: 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);
892: n = 5;
893: PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options","AMG options for vector Poisson","None",jac->as_amg_alpha_opts,&n,&flag2);
894: if (flag || flag2) {
895: PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0], /* AMG coarsen type */
896: jac->as_amg_alpha_opts[1], /* AMG agg_levels */
897: jac->as_amg_alpha_opts[2], /* AMG relax_type */
898: jac->as_amg_alpha_theta,
899: jac->as_amg_alpha_opts[3], /* AMG interp_type */
900: jac->as_amg_alpha_opts[4])); /* AMG Pmax */
901: }
902: 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);
903: n = 5;
904: PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options","AMG options for scalar Poisson solver","None",jac->as_amg_beta_opts,&n,&flag2);
905: if (flag || flag2) {
906: PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */
907: jac->as_amg_beta_opts[1], /* AMG agg_levels */
908: jac->as_amg_beta_opts[2], /* AMG relax_type */
909: jac->as_amg_beta_theta,
910: jac->as_amg_beta_opts[3], /* AMG interp_type */
911: jac->as_amg_beta_opts[4])); /* AMG Pmax */
912: }
913: PetscOptionsTail();
914: return(0);
915: }
919: static PetscErrorCode PCView_HYPRE_AMS(PC pc,PetscViewer viewer)
920: {
921: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
923: PetscBool iascii;
926: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
927: if (iascii) {
928: PetscViewerASCIIPrintf(viewer," HYPRE AMS preconditioning\n");
929: PetscViewerASCIIPrintf(viewer," HYPRE AMS: subspace iterations per application %d\n",jac->as_max_iter);
930: PetscViewerASCIIPrintf(viewer," HYPRE AMS: subspace cycle type %d\n",jac->ams_cycle_type);
931: PetscViewerASCIIPrintf(viewer," HYPRE AMS: subspace iteration tolerance %g\n",jac->as_tol);
932: PetscViewerASCIIPrintf(viewer," HYPRE AMS: smoother type %d\n",jac->as_relax_type);
933: PetscViewerASCIIPrintf(viewer," HYPRE AMS: number of smoothing steps %d\n",jac->as_relax_times);
934: PetscViewerASCIIPrintf(viewer," HYPRE AMS: smoother weight %g\n",jac->as_relax_weight);
935: PetscViewerASCIIPrintf(viewer," HYPRE AMS: smoother omega %g\n",jac->as_omega);
936: if (jac->alpha_Poisson) {
937: PetscViewerASCIIPrintf(viewer," HYPRE AMS: vector Poisson solver (passed in by user)\n");
938: } else {
939: PetscViewerASCIIPrintf(viewer," HYPRE AMS: vector Poisson solver (computed) \n");
940: }
941: PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);
942: PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
943: PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);
944: PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);
945: PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
946: PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);
947: if (!jac->ams_beta_is_zero) {
948: if (jac->beta_Poisson) {
949: PetscViewerASCIIPrintf(viewer," HYPRE AMS: scalar Poisson solver (passed in by user)\n");
950: } else {
951: PetscViewerASCIIPrintf(viewer," HYPRE AMS: scalar Poisson solver (computed) \n");
952: }
953: PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);
954: PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
955: PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);
956: PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);
957: PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
958: PetscViewerASCIIPrintf(viewer," HYPRE AMS: boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);
959: }
960: }
961: return(0);
962: }
966: static PetscErrorCode PCSetFromOptions_HYPRE_ADS(PetscOptionItems *PetscOptionsObject,PC pc)
967: {
968: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
970: PetscInt n;
971: PetscBool flag,flag2,flag3,flag4;
974: PetscOptionsHead(PetscOptionsObject,"HYPRE ADS Options");
975: PetscOptionsInt("-pc_hypre_ads_print_level","Debugging output level for ADS","None",jac->as_print,&jac->as_print,&flag);
976: if (flag) PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
977: PetscOptionsInt("-pc_hypre_ads_max_iter","Maximum number of ADS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
978: if (flag) PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
979: PetscOptionsInt("-pc_hypre_ads_cycle_type","Cycle type for ADS multigrid","None",jac->ads_cycle_type,&jac->ads_cycle_type,&flag);
980: if (flag) PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ads_cycle_type));
981: PetscOptionsReal("-pc_hypre_ads_tol","Error tolerance for ADS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
982: if (flag) PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
983: PetscOptionsInt("-pc_hypre_ads_relax_type","Relaxation type for ADS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
984: PetscOptionsInt("-pc_hypre_ads_relax_times","Number of relaxation steps for ADS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
985: PetscOptionsReal("-pc_hypre_ads_relax_weight","Relaxation weight for ADS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
986: PetscOptionsReal("-pc_hypre_ads_omega","SSOR coefficient for ADS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
987: if (flag || flag2 || flag3 || flag4) {
988: PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
989: jac->as_relax_times,
990: jac->as_relax_weight,
991: jac->as_omega));
992: }
993: 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);
994: n = 5;
995: PetscOptionsIntArray("-pc_hypre_ads_ams_options","AMG options for AMS solver inside ADS","None",jac->as_amg_alpha_opts,&n,&flag2);
996: PetscOptionsInt("-pc_hypre_ads_ams_cycle_type","Cycle type for AMS solver inside ADS","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag3);
997: if (flag || flag2 || flag3) {
998: PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type, /* AMS cycle type */
999: jac->as_amg_alpha_opts[0], /* AMG coarsen type */
1000: jac->as_amg_alpha_opts[1], /* AMG agg_levels */
1001: jac->as_amg_alpha_opts[2], /* AMG relax_type */
1002: jac->as_amg_alpha_theta,
1003: jac->as_amg_alpha_opts[3], /* AMG interp_type */
1004: jac->as_amg_alpha_opts[4])); /* AMG Pmax */
1005: }
1006: 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);
1007: n = 5;
1008: PetscOptionsIntArray("-pc_hypre_ads_amg_options","AMG options for vector AMG solver inside ADS","None",jac->as_amg_beta_opts,&n,&flag2);
1009: if (flag || flag2) {
1010: PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */
1011: jac->as_amg_beta_opts[1], /* AMG agg_levels */
1012: jac->as_amg_beta_opts[2], /* AMG relax_type */
1013: jac->as_amg_beta_theta,
1014: jac->as_amg_beta_opts[3], /* AMG interp_type */
1015: jac->as_amg_beta_opts[4])); /* AMG Pmax */
1016: }
1017: PetscOptionsTail();
1018: return(0);
1019: }
1023: static PetscErrorCode PCView_HYPRE_ADS(PC pc,PetscViewer viewer)
1024: {
1025: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1027: PetscBool iascii;
1030: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1031: if (iascii) {
1032: PetscViewerASCIIPrintf(viewer," HYPRE ADS preconditioning\n");
1033: PetscViewerASCIIPrintf(viewer," HYPRE ADS: subspace iterations per application %d\n",jac->as_max_iter);
1034: PetscViewerASCIIPrintf(viewer," HYPRE ADS: subspace cycle type %d\n",jac->ads_cycle_type);
1035: PetscViewerASCIIPrintf(viewer," HYPRE ADS: subspace iteration tolerance %g\n",jac->as_tol);
1036: PetscViewerASCIIPrintf(viewer," HYPRE ADS: smoother type %d\n",jac->as_relax_type);
1037: PetscViewerASCIIPrintf(viewer," HYPRE ADS: number of smoothing steps %d\n",jac->as_relax_times);
1038: PetscViewerASCIIPrintf(viewer," HYPRE ADS: smoother weight %g\n",jac->as_relax_weight);
1039: PetscViewerASCIIPrintf(viewer," HYPRE ADS: smoother omega %g\n",jac->as_omega);
1040: PetscViewerASCIIPrintf(viewer," HYPRE ADS: AMS solver\n");
1041: PetscViewerASCIIPrintf(viewer," HYPRE ADS: subspace cycle type %d\n",jac->ams_cycle_type);
1042: PetscViewerASCIIPrintf(viewer," HYPRE ADS: boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);
1043: PetscViewerASCIIPrintf(viewer," HYPRE ADS: boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
1044: PetscViewerASCIIPrintf(viewer," HYPRE ADS: boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);
1045: PetscViewerASCIIPrintf(viewer," HYPRE ADS: boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);
1046: PetscViewerASCIIPrintf(viewer," HYPRE ADS: boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
1047: PetscViewerASCIIPrintf(viewer," HYPRE ADS: boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);
1048: PetscViewerASCIIPrintf(viewer," HYPRE ADS: vector Poisson solver\n");
1049: PetscViewerASCIIPrintf(viewer," HYPRE ADS: boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);
1050: PetscViewerASCIIPrintf(viewer," HYPRE ADS: boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
1051: PetscViewerASCIIPrintf(viewer," HYPRE ADS: boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);
1052: PetscViewerASCIIPrintf(viewer," HYPRE ADS: boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);
1053: PetscViewerASCIIPrintf(viewer," HYPRE ADS: boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
1054: PetscViewerASCIIPrintf(viewer," HYPRE ADS: boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);
1055: }
1056: return(0);
1057: }
1061: static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE(PC pc, Mat G)
1062: {
1063: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1064: HYPRE_ParCSRMatrix parcsr_G;
1065: PetscErrorCode ierr;
1068: /* throw away any discrete gradient if already set */
1069: if (jac->G) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->G));
1070: MatHYPRE_IJMatrixCreate(G,&jac->G);
1071: MatHYPRE_IJMatrixCopy(G,jac->G);
1072: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->G,(void**)(&parcsr_G)));
1073: PetscStackCall("Hypre set gradient",(*jac->setdgrad)(jac->hsolver,parcsr_G););
1074: return(0);
1075: }
1079: /*@
1080: PCHYPRESetDiscreteGradient - Set discrete gradient matrix
1082: Collective on PC
1084: Input Parameters:
1085: + pc - the preconditioning context
1086: - G - the discrete gradient
1088: Level: intermediate
1090: Notes: G should have as many rows as the number of edges and as many columns as the number of vertices in the mesh
1091: 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
1093: .seealso:
1094: @*/
1095: PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G)
1096: {
1103: PetscTryMethod(pc,"PCHYPRESetDiscreteGradient_C",(PC,Mat),(pc,G));
1104: return(0);
1105: }
1109: static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C)
1110: {
1111: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1112: HYPRE_ParCSRMatrix parcsr_C;
1113: PetscErrorCode ierr;
1116: /* throw away any discrete curl if already set */
1117: if (jac->C) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->C));
1118: MatHYPRE_IJMatrixCreate(C,&jac->C);
1119: MatHYPRE_IJMatrixCopy(C,jac->C);
1120: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->C,(void**)(&parcsr_C)));
1121: PetscStackCall("Hypre set curl",(*jac->setdcurl)(jac->hsolver,parcsr_C););
1122: return(0);
1123: }
1127: /*@
1128: PCHYPRESetDiscreteCurl - Set discrete curl matrix
1130: Collective on PC
1132: Input Parameters:
1133: + pc - the preconditioning context
1134: - C - the discrete curl
1136: Level: intermediate
1138: Notes: C should have as many rows as the number of faces and as many columns as the number of edges in the mesh
1139: 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
1141: .seealso:
1142: @*/
1143: PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C)
1144: {
1151: PetscTryMethod(pc,"PCHYPRESetDiscreteCurl_C",(PC,Mat),(pc,C));
1152: return(0);
1153: }
1157: static PetscErrorCode PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS(PC pc, Mat A)
1158: {
1159: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1160: HYPRE_ParCSRMatrix parcsr_alpha_Poisson;
1161: PetscErrorCode ierr;
1164: /* throw away any matrix if already set */
1165: if (jac->alpha_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->alpha_Poisson));
1166: MatHYPRE_IJMatrixCreate(A,&jac->alpha_Poisson);
1167: MatHYPRE_IJMatrixCopy(A,jac->alpha_Poisson);
1168: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->alpha_Poisson,(void**)(&parcsr_alpha_Poisson)));
1169: PetscStackCallStandard(HYPRE_AMSSetAlphaPoissonMatrix,(jac->hsolver,parcsr_alpha_Poisson));
1170: return(0);
1171: }
1175: /*@
1176: PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix
1178: Collective on PC
1180: Input Parameters:
1181: + pc - the preconditioning context
1182: - A - the matrix
1184: Level: intermediate
1186: Notes: A should be obtained by discretizing the vector valued Poisson problem with linear finite elements
1188: .seealso:
1189: @*/
1190: PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A)
1191: {
1198: PetscTryMethod(pc,"PCHYPRESetAlphaPoissonMatrix_C",(PC,Mat),(pc,A));
1199: return(0);
1200: }
1204: static PetscErrorCode PCHYPRESetBetaPoissonMatrix_HYPRE_AMS(PC pc, Mat A)
1205: {
1206: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1207: HYPRE_ParCSRMatrix parcsr_beta_Poisson;
1208: PetscErrorCode ierr;
1211: if (!A) {
1212: PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,NULL));
1213: jac->ams_beta_is_zero = PETSC_TRUE;
1214: return(0);
1215: }
1216: jac->ams_beta_is_zero = PETSC_FALSE;
1217: /* throw away any matrix if already set */
1218: if (jac->beta_Poisson) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->beta_Poisson));
1219: MatHYPRE_IJMatrixCreate(A,&jac->beta_Poisson);
1220: MatHYPRE_IJMatrixCopy(A,jac->beta_Poisson);
1221: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->beta_Poisson,(void**)(&parcsr_beta_Poisson)));
1222: PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,parcsr_beta_Poisson));
1223: return(0);
1224: }
1228: /*@
1229: PCHYPRESetBetaPoissonMatrix - Set Poisson matrix
1231: Collective on PC
1233: Input Parameters:
1234: + pc - the preconditioning context
1235: - A - the matrix
1237: Level: intermediate
1239: Notes: A should be obtained by discretizing the Poisson problem with linear finite elements.
1240: Following HYPRE convention, the scalar Poisson solver of AMS can be turned off by passing NULL.
1242: .seealso:
1243: @*/
1244: PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A)
1245: {
1250: if (A) {
1253: }
1254: PetscTryMethod(pc,"PCHYPRESetBetaPoissonMatrix_C",(PC,Mat),(pc,A));
1255: return(0);
1256: }
1260: static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE_AMS(PC pc,Vec ozz, Vec zoz, Vec zzo)
1261: {
1262: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1263: HYPRE_ParVector par_ozz,par_zoz,par_zzo;
1264: PetscInt dim;
1265: PetscErrorCode ierr;
1268: /* throw away any vector if already set */
1269: if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0]));
1270: if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1]));
1271: if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2]));
1272: jac->constants[0] = NULL;
1273: jac->constants[1] = NULL;
1274: jac->constants[2] = NULL;
1275: VecHYPRE_IJVectorCreate(ozz,&jac->constants[0]);
1276: VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);
1277: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[0],(void**)(&par_ozz)));
1278: VecHYPRE_IJVectorCreate(zoz,&jac->constants[1]);
1279: VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);
1280: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[1],(void**)(&par_zoz)));
1281: dim = 2;
1282: if (zzo) {
1283: VecHYPRE_IJVectorCreate(zzo,&jac->constants[2]);
1284: VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);
1285: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[2],(void**)(&par_zzo)));
1286: dim++;
1287: }
1288: PetscStackCallStandard(HYPRE_AMSSetEdgeConstantVectors,(jac->hsolver,par_ozz,par_zoz,par_zzo));
1289: PetscStackCallStandard(HYPRE_AMSSetDimension,(jac->hsolver,dim));
1290: return(0);
1291: }
1295: /*@
1296: PCHYPRESetEdgeConstantVectors - Set the representation of the constant vector fields in edge element basis
1298: Collective on PC
1300: Input Parameters:
1301: + pc - the preconditioning context
1302: - ozz - vector representing (1,0,0) (or (1,0) in 2D)
1303: - zoz - vector representing (0,1,0) (or (0,1) in 2D)
1304: - zzo - vector representing (0,0,1) (use NULL in 2D)
1306: Level: intermediate
1308: Notes:
1310: .seealso:
1311: @*/
1312: PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo)
1313: {
1324: PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));
1325: return(0);
1326: }
1330: static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1331: {
1332: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1333: Vec tv;
1334: HYPRE_ParVector par_coords[3];
1335: PetscInt i;
1336: PetscErrorCode ierr;
1339: /* throw away any coordinate vector if already set */
1340: if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0]));
1341: if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1]));
1342: if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2]));
1343: /* set problem's dimension */
1344: if (jac->setdim) {
1345: PetscStackCall("Hypre set dim",(*jac->setdim)(jac->hsolver,dim););
1346: }
1347: /* compute IJ vector for coordinates */
1348: VecCreate(PetscObjectComm((PetscObject)pc),&tv);
1349: VecSetType(tv,VECSTANDARD);
1350: VecSetSizes(tv,nloc,PETSC_DECIDE);
1351: for (i=0;i<dim;i++) {
1352: PetscScalar *array;
1353: PetscInt j;
1355: VecHYPRE_IJVectorCreate(tv,&jac->coords[i]);
1356: VecGetArray(tv,&array);
1357: for (j=0;j<nloc;j++) {
1358: array[j] = coords[j*dim+i];
1359: }
1360: PetscStackCallStandard(HYPRE_IJVectorSetValues,(jac->coords[i],nloc,NULL,array));
1361: PetscStackCallStandard(HYPRE_IJVectorAssemble,(jac->coords[i]));
1362: VecRestoreArray(tv,&array);
1363: }
1364: VecDestroy(&tv);
1365: /* pass parCSR vectors to AMS solver */
1366: par_coords[0] = NULL;
1367: par_coords[1] = NULL;
1368: par_coords[2] = NULL;
1369: if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0],(void**)(&par_coords[0])));
1370: if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1],(void**)(&par_coords[1])));
1371: if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2],(void**)(&par_coords[2])));
1372: PetscStackCall("Hypre set coords",(*jac->setcoord)(jac->hsolver,par_coords[0],par_coords[1],par_coords[2]););
1373: return(0);
1374: }
1376: /* ---------------------------------------------------------------------------------*/
1380: static PetscErrorCode PCHYPREGetType_HYPRE(PC pc,const char *name[])
1381: {
1382: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1385: *name = jac->hypre_type;
1386: return(0);
1387: }
1391: static PetscErrorCode PCHYPRESetType_HYPRE(PC pc,const char name[])
1392: {
1393: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
1395: PetscBool flag;
1398: if (jac->hypre_type) {
1399: PetscStrcmp(jac->hypre_type,name,&flag);
1400: if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set");
1401: return(0);
1402: } else {
1403: PetscStrallocpy(name, &jac->hypre_type);
1404: }
1406: jac->maxiter = PETSC_DEFAULT;
1407: jac->tol = PETSC_DEFAULT;
1408: jac->printstatistics = PetscLogPrintInfo;
1410: PetscStrcmp("pilut",jac->hypre_type,&flag);
1411: if (flag) {
1412: PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver));
1413: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
1414: pc->ops->view = PCView_HYPRE_Pilut;
1415: jac->destroy = HYPRE_ParCSRPilutDestroy;
1416: jac->setup = HYPRE_ParCSRPilutSetup;
1417: jac->solve = HYPRE_ParCSRPilutSolve;
1418: jac->factorrowsize = PETSC_DEFAULT;
1419: return(0);
1420: }
1421: PetscStrcmp("parasails",jac->hypre_type,&flag);
1422: if (flag) {
1423: PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver));
1424: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
1425: pc->ops->view = PCView_HYPRE_ParaSails;
1426: jac->destroy = HYPRE_ParaSailsDestroy;
1427: jac->setup = HYPRE_ParaSailsSetup;
1428: jac->solve = HYPRE_ParaSailsSolve;
1429: /* initialize */
1430: jac->nlevels = 1;
1431: jac->threshhold = .1;
1432: jac->filter = .1;
1433: jac->loadbal = 0;
1434: if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE;
1435: else jac->logging = (int) PETSC_FALSE;
1437: jac->ruse = (int) PETSC_FALSE;
1438: jac->symt = 0;
1439: PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
1440: PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
1441: PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
1442: PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
1443: PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
1444: PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
1445: return(0);
1446: }
1447: PetscStrcmp("boomeramg",jac->hypre_type,&flag);
1448: if (flag) {
1449: HYPRE_BoomerAMGCreate(&jac->hsolver);
1450: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_BoomerAMG;
1451: pc->ops->view = PCView_HYPRE_BoomerAMG;
1452: pc->ops->applytranspose = PCApplyTranspose_HYPRE_BoomerAMG;
1453: pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
1454: jac->destroy = HYPRE_BoomerAMGDestroy;
1455: jac->setup = HYPRE_BoomerAMGSetup;
1456: jac->solve = HYPRE_BoomerAMGSolve;
1457: jac->applyrichardson = PETSC_FALSE;
1458: /* these defaults match the hypre defaults */
1459: jac->cycletype = 1;
1460: jac->maxlevels = 25;
1461: jac->maxiter = 1;
1462: jac->tol = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
1463: jac->truncfactor = 0.0;
1464: jac->strongthreshold = .25;
1465: jac->maxrowsum = .9;
1466: jac->coarsentype = 6;
1467: jac->measuretype = 0;
1468: jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
1469: jac->smoothtype = -1; /* Not set by default */
1470: jac->smoothnumlevels = 25;
1471: jac->eu_level = 0;
1472: jac->eu_droptolerance = 0;
1473: jac->eu_bj = 0;
1474: jac->relaxtype[0] = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */
1475: jac->relaxtype[2] = 9; /*G.E. */
1476: jac->relaxweight = 1.0;
1477: jac->outerrelaxweight = 1.0;
1478: jac->relaxorder = 1;
1479: jac->interptype = 0;
1480: jac->agg_nl = 0;
1481: jac->pmax = 0;
1482: jac->truncfactor = 0.0;
1483: jac->agg_num_paths = 1;
1485: jac->nodal_coarsen = 0;
1486: jac->nodal_relax = PETSC_FALSE;
1487: jac->nodal_relax_levels = 1;
1488: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
1489: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
1490: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
1491: PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
1492: PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
1493: PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
1494: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
1495: PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
1496: PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
1497: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
1498: PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
1499: PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
1500: PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
1501: PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
1502: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0])); /*defaults coarse to 9*/
1503: PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */
1504: return(0);
1505: }
1506: PetscStrcmp("ams",jac->hypre_type,&flag);
1507: if (flag) {
1508: HYPRE_AMSCreate(&jac->hsolver);
1509: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_AMS;
1510: pc->ops->view = PCView_HYPRE_AMS;
1511: jac->destroy = HYPRE_AMSDestroy;
1512: jac->setup = HYPRE_AMSSetup;
1513: jac->solve = HYPRE_AMSSolve;
1514: jac->setdgrad = HYPRE_AMSSetDiscreteGradient;
1515: jac->setcoord = HYPRE_AMSSetCoordinateVectors;
1516: jac->setdim = HYPRE_AMSSetDimension;
1517: jac->coords[0] = NULL;
1518: jac->coords[1] = NULL;
1519: jac->coords[2] = NULL;
1520: jac->G = NULL;
1521: /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */
1522: jac->as_print = 0;
1523: jac->as_max_iter = 1; /* used as a preconditioner */
1524: jac->as_tol = 0.; /* used as a preconditioner */
1525: jac->ams_cycle_type = 13;
1526: /* Smoothing options */
1527: jac->as_relax_type = 2;
1528: jac->as_relax_times = 1;
1529: jac->as_relax_weight = 1.0;
1530: jac->as_omega = 1.0;
1531: /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1532: jac->as_amg_alpha_opts[0] = 10;
1533: jac->as_amg_alpha_opts[1] = 1;
1534: jac->as_amg_alpha_opts[2] = 6;
1535: jac->as_amg_alpha_opts[3] = 6;
1536: jac->as_amg_alpha_opts[4] = 4;
1537: jac->as_amg_alpha_theta = 0.25;
1538: /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1539: jac->ams_beta_is_zero = PETSC_FALSE;
1540: jac->as_amg_beta_opts[0] = 10;
1541: jac->as_amg_beta_opts[1] = 1;
1542: jac->as_amg_beta_opts[2] = 6;
1543: jac->as_amg_beta_opts[3] = 6;
1544: jac->as_amg_beta_opts[4] = 4;
1545: jac->as_amg_beta_theta = 0.25;
1546: PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
1547: PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1548: PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1549: PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
1550: PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1551: jac->as_relax_times,
1552: jac->as_relax_weight,
1553: jac->as_omega));
1554: PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0], /* AMG coarsen type */
1555: jac->as_amg_alpha_opts[1], /* AMG agg_levels */
1556: jac->as_amg_alpha_opts[2], /* AMG relax_type */
1557: jac->as_amg_alpha_theta,
1558: jac->as_amg_alpha_opts[3], /* AMG interp_type */
1559: jac->as_amg_alpha_opts[4])); /* AMG Pmax */
1560: PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */
1561: jac->as_amg_beta_opts[1], /* AMG agg_levels */
1562: jac->as_amg_beta_opts[2], /* AMG relax_type */
1563: jac->as_amg_beta_theta,
1564: jac->as_amg_beta_opts[3], /* AMG interp_type */
1565: jac->as_amg_beta_opts[4])); /* AMG Pmax */
1566: PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);
1567: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);
1568: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE_AMS);
1569: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetAlphaPoissonMatrix_C",PCHYPRESetAlphaPoissonMatrix_HYPRE_AMS);
1570: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetBetaPoissonMatrix_C",PCHYPRESetBetaPoissonMatrix_HYPRE_AMS);
1571: return(0);
1572: }
1573: PetscStrcmp("ads",jac->hypre_type,&flag);
1574: if (flag) {
1575: HYPRE_ADSCreate(&jac->hsolver);
1576: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ADS;
1577: pc->ops->view = PCView_HYPRE_ADS;
1578: jac->destroy = HYPRE_ADSDestroy;
1579: jac->setup = HYPRE_ADSSetup;
1580: jac->solve = HYPRE_ADSSolve;
1581: jac->setdgrad = HYPRE_ADSSetDiscreteGradient;
1582: jac->setdcurl = HYPRE_ADSSetDiscreteCurl;
1583: jac->setcoord = HYPRE_ADSSetCoordinateVectors;
1584: jac->coords[0] = NULL;
1585: jac->coords[1] = NULL;
1586: jac->coords[2] = NULL;
1587: jac->G = NULL;
1588: jac->C = NULL;
1589: /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */
1590: jac->as_print = 0;
1591: jac->as_max_iter = 1; /* used as a preconditioner */
1592: jac->as_tol = 0.; /* used as a preconditioner */
1593: jac->ads_cycle_type = 13;
1594: /* Smoothing options */
1595: jac->as_relax_type = 2;
1596: jac->as_relax_times = 1;
1597: jac->as_relax_weight = 1.0;
1598: jac->as_omega = 1.0;
1599: /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */
1600: jac->ams_cycle_type = 14;
1601: jac->as_amg_alpha_opts[0] = 10;
1602: jac->as_amg_alpha_opts[1] = 1;
1603: jac->as_amg_alpha_opts[2] = 6;
1604: jac->as_amg_alpha_opts[3] = 6;
1605: jac->as_amg_alpha_opts[4] = 4;
1606: jac->as_amg_alpha_theta = 0.25;
1607: /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1608: jac->as_amg_beta_opts[0] = 10;
1609: jac->as_amg_beta_opts[1] = 1;
1610: jac->as_amg_beta_opts[2] = 6;
1611: jac->as_amg_beta_opts[3] = 6;
1612: jac->as_amg_beta_opts[4] = 4;
1613: jac->as_amg_beta_theta = 0.25;
1614: PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
1615: PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1616: PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1617: PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
1618: PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1619: jac->as_relax_times,
1620: jac->as_relax_weight,
1621: jac->as_omega));
1622: PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type, /* AMG coarsen type */
1623: jac->as_amg_alpha_opts[0], /* AMG coarsen type */
1624: jac->as_amg_alpha_opts[1], /* AMG agg_levels */
1625: jac->as_amg_alpha_opts[2], /* AMG relax_type */
1626: jac->as_amg_alpha_theta,
1627: jac->as_amg_alpha_opts[3], /* AMG interp_type */
1628: jac->as_amg_alpha_opts[4])); /* AMG Pmax */
1629: PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0], /* AMG coarsen type */
1630: jac->as_amg_beta_opts[1], /* AMG agg_levels */
1631: jac->as_amg_beta_opts[2], /* AMG relax_type */
1632: jac->as_amg_beta_theta,
1633: jac->as_amg_beta_opts[3], /* AMG interp_type */
1634: jac->as_amg_beta_opts[4])); /* AMG Pmax */
1635: PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);
1636: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);
1637: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",PCHYPRESetDiscreteCurl_HYPRE);
1638: return(0);
1639: }
1640: PetscFree(jac->hypre_type);
1642: jac->hypre_type = NULL;
1643: SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are pilut, parasails, boomeramg, ams",name);
1644: return(0);
1645: }
1647: /*
1648: It only gets here if the HYPRE type has not been set before the call to
1649: ...SetFromOptions() which actually is most of the time
1650: */
1653: static PetscErrorCode PCSetFromOptions_HYPRE(PetscOptionItems *PetscOptionsObject,PC pc)
1654: {
1656: PetscInt indx;
1657: const char *type[] = {"pilut","parasails","boomeramg","ams","ads"};
1658: PetscBool flg;
1661: PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");
1662: PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"boomeramg",&indx,&flg);
1663: if (flg) {
1664: PCHYPRESetType_HYPRE(pc,type[indx]);
1665: } else {
1666: PCHYPRESetType_HYPRE(pc,"boomeramg");
1667: }
1668: if (pc->ops->setfromoptions) {
1669: pc->ops->setfromoptions(PetscOptionsObject,pc);
1670: }
1671: PetscOptionsTail();
1672: return(0);
1673: }
1677: /*@C
1678: PCHYPRESetType - Sets which hypre preconditioner you wish to use
1680: Input Parameters:
1681: + pc - the preconditioner context
1682: - name - either pilut, parasails, boomeramg, ams, ads
1684: Options Database Keys:
1685: -pc_hypre_type - One of pilut, parasails, boomeramg, ams, ads
1687: Level: intermediate
1689: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
1690: PCHYPRE
1692: @*/
1693: PetscErrorCode PCHYPRESetType(PC pc,const char name[])
1694: {
1700: PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));
1701: return(0);
1702: }
1706: /*@C
1707: PCHYPREGetType - Gets which hypre preconditioner you are using
1709: Input Parameter:
1710: . pc - the preconditioner context
1712: Output Parameter:
1713: . name - either pilut, parasails, boomeramg, ams, ads
1715: Level: intermediate
1717: .seealso: PCCreate(), PCHYPRESetType(), PCType (for list of available types), PC,
1718: PCHYPRE
1720: @*/
1721: PetscErrorCode PCHYPREGetType(PC pc,const char *name[])
1722: {
1728: PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));
1729: return(0);
1730: }
1732: /*MC
1733: PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre
1735: Options Database Keys:
1736: + -pc_hypre_type - One of pilut, parasails, boomeramg, ams, ads
1737: - Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX
1738: preconditioner
1740: Level: intermediate
1742: Notes: Apart from pc_hypre_type (for which there is PCHYPRESetType()),
1743: the many hypre options can ONLY be set via the options database (e.g. the command line
1744: or with PetscOptionsSetValue(), there are no functions to set them)
1746: The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_rtol refer to the number of iterations
1747: (V-cycles) and tolerance that boomeramg does EACH time it is called. So for example, if
1748: -pc_hypre_boomeramg_max_iter is set to 2 then 2-V-cycles are being used to define the preconditioner
1749: (-pc_hypre_boomeramg_rtol should be set to 0.0 - the default - to strictly use a fixed number of
1750: iterations per hypre call). -ksp_max_it and -ksp_rtol STILL determine the total number of iterations
1751: and tolerance for the Krylov solver. For example, if -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10
1752: then AT MOST twenty V-cycles of boomeramg will be called.
1754: Note that the option -pc_hypre_boomeramg_relax_type_all defaults to symmetric relaxation
1755: (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry.
1756: Otherwise, you may want to use -pc_hypre_boomeramg_relax_type_all SOR/Jacobi.
1757: If you wish to use BoomerAMG WITHOUT a Krylov method use -ksp_type richardson NOT -ksp_type preonly
1758: and use -ksp_max_it to control the number of V-cycles.
1759: (see the PETSc FAQ.html at the PETSc website under the Documentation tab).
1761: 2007-02-03 Using HYPRE-1.11.1b, the routine HYPRE_BoomerAMGSolveT and the option
1762: -pc_hypre_parasails_reuse were failing with SIGSEGV. Dalcin L.
1764: MatSetNearNullSpace() - if you provide a near null space to your matrix it is ignored by hypre UNLESS you also use
1765: the two options:
1766: + -pc_hypre_boomeramg_nodal_coarsen <n> - where n is from 1 to 6 (see HYPRE_BOOMERAMGSetNodal())
1767: - -pc_hypre_boomeramg_vec_interp_variant <v> where v is from 1 to 3 (see HYPRE_BoomerAMGSetInterpVecVariant())
1769: Depending on the linear system you may see the same or different convergence depending on the values you use.
1771: See PCPFMG for access to the hypre Struct PFMG solver
1773: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
1774: PCHYPRESetType(), PCPFMG
1776: M*/
1780: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
1781: {
1782: PC_HYPRE *jac;
1786: PetscNewLog(pc,&jac);
1788: pc->data = jac;
1789: pc->ops->reset = PCReset_HYPRE;
1790: pc->ops->destroy = PCDestroy_HYPRE;
1791: pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
1792: pc->ops->setup = PCSetUp_HYPRE;
1793: pc->ops->apply = PCApply_HYPRE;
1794: jac->comm_hypre = MPI_COMM_NULL;
1795: jac->hypre_type = NULL;
1796: jac->coords[0] = NULL;
1797: jac->coords[1] = NULL;
1798: jac->coords[2] = NULL;
1799: jac->constants[0] = NULL;
1800: jac->constants[1] = NULL;
1801: jac->constants[2] = NULL;
1802: jac->G = NULL;
1803: jac->C = NULL;
1804: jac->alpha_Poisson = NULL;
1805: jac->beta_Poisson = NULL;
1806: jac->setdim = NULL;
1807: jac->hmnull = NULL;
1808: jac->n_hmnull = 0;
1809: /* duplicate communicator for hypre */
1810: MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1811: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);
1812: PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);
1813: return(0);
1814: }
1816: /* ---------------------------------------------------------------------------------------------------------------------------------*/
1818: /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */
1819: #include <petsc/private/matimpl.h>
1821: typedef struct {
1822: MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */
1823: HYPRE_StructSolver hsolver;
1825: /* keep copy of PFMG options used so may view them */
1826: PetscInt its;
1827: double tol;
1828: PetscInt relax_type;
1829: PetscInt rap_type;
1830: PetscInt num_pre_relax,num_post_relax;
1831: PetscInt max_levels;
1832: } PC_PFMG;
1836: PetscErrorCode PCDestroy_PFMG(PC pc)
1837: {
1839: PC_PFMG *ex = (PC_PFMG*) pc->data;
1842: if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1843: MPI_Comm_free(&ex->hcomm);
1844: PetscFree(pc->data);
1845: return(0);
1846: }
1848: static const char *PFMGRelaxType[] = {"Jacobi","Weighted-Jacobi","symmetric-Red/Black-Gauss-Seidel","Red/Black-Gauss-Seidel"};
1849: static const char *PFMGRAPType[] = {"Galerkin","non-Galerkin"};
1853: PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
1854: {
1856: PetscBool iascii;
1857: PC_PFMG *ex = (PC_PFMG*) pc->data;
1860: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1861: if (iascii) {
1862: PetscViewerASCIIPrintf(viewer," HYPRE PFMG preconditioning\n");
1863: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: max iterations %d\n",ex->its);
1864: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: tolerance %g\n",ex->tol);
1865: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
1866: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: RAP type %s\n",PFMGRAPType[ex->rap_type]);
1867: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
1868: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: max levels %d\n",ex->max_levels);
1869: }
1870: return(0);
1871: }
1876: PetscErrorCode PCSetFromOptions_PFMG(PetscOptionItems *PetscOptionsObject,PC pc)
1877: {
1879: PC_PFMG *ex = (PC_PFMG*) pc->data;
1880: PetscBool flg = PETSC_FALSE;
1883: PetscOptionsHead(PetscOptionsObject,"PFMG options");
1884: PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);
1885: if (flg) {
1886: int level=3;
1887: PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level));
1888: }
1889: PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);
1890: PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its));
1891: PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
1892: PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax));
1893: PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
1894: PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax));
1896: PetscOptionsInt("-pc_pfmg_max_levels","Max Levels for MG hierarchy","HYPRE_StructPFMGSetMaxLevels",ex->max_levels,&ex->max_levels,NULL);
1897: PetscStackCallStandard(HYPRE_StructPFMGSetMaxLevels,(ex->hsolver,ex->max_levels));
1899: PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);
1900: PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol));
1901: 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);
1902: PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type));
1903: PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);
1904: PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type));
1905: PetscOptionsTail();
1906: return(0);
1907: }
1911: PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
1912: {
1913: PetscErrorCode ierr;
1914: PC_PFMG *ex = (PC_PFMG*) pc->data;
1915: PetscScalar *yy;
1916: const PetscScalar *xx;
1917: PetscInt ilower[3],iupper[3];
1918: Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
1921: PetscCitationsRegister(hypreCitation,&cite);
1922: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
1923: iupper[0] += ilower[0] - 1;
1924: iupper[1] += ilower[1] - 1;
1925: iupper[2] += ilower[2] - 1;
1927: /* copy x values over to hypre */
1928: PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
1929: VecGetArrayRead(x,&xx);
1930: PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,(PetscScalar*)xx));
1931: VecRestoreArrayRead(x,&xx);
1932: PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
1933: PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
1935: /* copy solution values back to PETSc */
1936: VecGetArray(y,&yy);
1937: PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy));
1938: VecRestoreArray(y,&yy);
1939: return(0);
1940: }
1944: 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)
1945: {
1946: PC_PFMG *jac = (PC_PFMG*)pc->data;
1948: PetscInt oits;
1951: PetscCitationsRegister(hypreCitation,&cite);
1952: PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its));
1953: PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol));
1955: PCApply_PFMG(pc,b,y);
1956: PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
1957: *outits = oits;
1958: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1959: else *reason = PCRICHARDSON_CONVERGED_RTOL;
1960: PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol));
1961: PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its));
1962: return(0);
1963: }
1968: PetscErrorCode PCSetUp_PFMG(PC pc)
1969: {
1970: PetscErrorCode ierr;
1971: PC_PFMG *ex = (PC_PFMG*) pc->data;
1972: Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
1973: PetscBool flg;
1976: PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESTRUCT,&flg);
1977: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESTRUCT with this preconditioner");
1979: /* create the hypre solver object and set its information */
1980: if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1981: PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1982: PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
1983: PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver));
1984: return(0);
1985: }
1988: /*MC
1989: PCPFMG - the hypre PFMG multigrid solver
1991: Level: advanced
1993: Options Database:
1994: + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner
1995: . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
1996: . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid
1997: . -pc_pfmg_tol <tol> tolerance of PFMG
1998: . -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
1999: - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin
2001: Notes: This is for CELL-centered descretizations
2003: This must be used with the MATHYPRESTRUCT matrix type.
2004: This is less general than in hypre, it supports only one block per process defined by a PETSc DMDA.
2006: .seealso: PCMG, MATHYPRESTRUCT
2007: M*/
2011: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
2012: {
2014: PC_PFMG *ex;
2017: PetscNew(&ex); \
2018: pc->data = ex;
2020: ex->its = 1;
2021: ex->tol = 1.e-8;
2022: ex->relax_type = 1;
2023: ex->rap_type = 0;
2024: ex->num_pre_relax = 1;
2025: ex->num_post_relax = 1;
2026: ex->max_levels = 0;
2028: pc->ops->setfromoptions = PCSetFromOptions_PFMG;
2029: pc->ops->view = PCView_PFMG;
2030: pc->ops->destroy = PCDestroy_PFMG;
2031: pc->ops->apply = PCApply_PFMG;
2032: pc->ops->applyrichardson = PCApplyRichardson_PFMG;
2033: pc->ops->setup = PCSetUp_PFMG;
2035: MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
2036: PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
2037: return(0);
2038: }
2040: /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/
2042: /* we know we are working with a HYPRE_SStructMatrix */
2043: typedef struct {
2044: MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
2045: HYPRE_SStructSolver ss_solver;
2047: /* keep copy of SYSPFMG options used so may view them */
2048: PetscInt its;
2049: double tol;
2050: PetscInt relax_type;
2051: PetscInt num_pre_relax,num_post_relax;
2052: } PC_SysPFMG;
2056: PetscErrorCode PCDestroy_SysPFMG(PC pc)
2057: {
2059: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
2062: if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2063: MPI_Comm_free(&ex->hcomm);
2064: PetscFree(pc->data);
2065: return(0);
2066: }
2068: static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi","Red/Black-Gauss-Seidel"};
2072: PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
2073: {
2075: PetscBool iascii;
2076: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
2079: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2080: if (iascii) {
2081: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG preconditioning\n");
2082: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: max iterations %d\n",ex->its);
2083: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: tolerance %g\n",ex->tol);
2084: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
2085: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
2086: }
2087: return(0);
2088: }
2093: PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptionItems *PetscOptionsObject,PC pc)
2094: {
2096: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
2097: PetscBool flg = PETSC_FALSE;
2100: PetscOptionsHead(PetscOptionsObject,"SysPFMG options");
2101: PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);
2102: if (flg) {
2103: int level=3;
2104: PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level));
2105: }
2106: PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);
2107: PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its));
2108: PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
2109: PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax));
2110: PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
2111: PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax));
2113: PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);
2114: PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol));
2115: PetscOptionsEList("-pc_syspfmg_relax_type","Relax type for the up and down cycles","HYPRE_SStructSysPFMGSetRelaxType",SysPFMGRelaxType,4,SysPFMGRelaxType[ex->relax_type],&ex->relax_type,NULL);
2116: PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type));
2117: PetscOptionsTail();
2118: return(0);
2119: }
2123: PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
2124: {
2125: PetscErrorCode ierr;
2126: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
2127: PetscScalar *yy;
2128: const PetscScalar *xx;
2129: PetscInt ilower[3],iupper[3];
2130: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
2131: PetscInt ordering= mx->dofs_order;
2132: PetscInt nvars = mx->nvars;
2133: PetscInt part = 0;
2134: PetscInt size;
2135: PetscInt i;
2138: PetscCitationsRegister(hypreCitation,&cite);
2139: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
2140: iupper[0] += ilower[0] - 1;
2141: iupper[1] += ilower[1] - 1;
2142: iupper[2] += ilower[2] - 1;
2144: size = 1;
2145: for (i= 0; i< 3; i++) size *= (iupper[i]-ilower[i]+1);
2147: /* copy x values over to hypre for variable ordering */
2148: if (ordering) {
2149: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2150: VecGetArrayRead(x,&xx);
2151: for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,(PetscScalar*)xx+(size*i)));
2152: VecRestoreArrayRead(x,&xx);
2153: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2154: PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
2155: PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2157: /* copy solution values back to PETSc */
2158: VecGetArray(y,&yy);
2159: for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i)));
2160: VecRestoreArray(y,&yy);
2161: } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
2162: PetscScalar *z;
2163: PetscInt j, k;
2165: PetscMalloc1(nvars*size,&z);
2166: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2167: VecGetArrayRead(x,&xx);
2169: /* transform nodal to hypre's variable ordering for sys_pfmg */
2170: for (i= 0; i< size; i++) {
2171: k= i*nvars;
2172: for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
2173: }
2174: for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
2175: VecRestoreArrayRead(x,&xx);
2176: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2177: PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2179: /* copy solution values back to PETSc */
2180: VecGetArray(y,&yy);
2181: for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
2182: /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
2183: for (i= 0; i< size; i++) {
2184: k= i*nvars;
2185: for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
2186: }
2187: VecRestoreArray(y,&yy);
2188: PetscFree(z);
2189: }
2190: return(0);
2191: }
2195: 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)
2196: {
2197: PC_SysPFMG *jac = (PC_SysPFMG*)pc->data;
2199: PetscInt oits;
2202: PetscCitationsRegister(hypreCitation,&cite);
2203: PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its));
2204: PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol));
2205: PCApply_SysPFMG(pc,b,y);
2206: PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,(HYPRE_Int *)&oits));
2207: *outits = oits;
2208: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2209: else *reason = PCRICHARDSON_CONVERGED_RTOL;
2210: PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol));
2211: PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its));
2212: return(0);
2213: }
2218: PetscErrorCode PCSetUp_SysPFMG(PC pc)
2219: {
2220: PetscErrorCode ierr;
2221: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
2222: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
2223: PetscBool flg;
2226: PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESSTRUCT,&flg);
2227: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESSTRUCT with this preconditioner");
2229: /* create the hypre sstruct solver object and set its information */
2230: if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2231: PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2232: PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver));
2233: PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2234: return(0);
2235: }
2238: /*MC
2239: PCSysPFMG - the hypre SysPFMG multigrid solver
2241: Level: advanced
2243: Options Database:
2244: + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner
2245: . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2246: . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2247: . -pc_syspfmg_tol <tol> tolerance of SysPFMG
2248: . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel
2250: Notes: This is for CELL-centered descretizations
2252: This must be used with the MATHYPRESSTRUCT matrix type.
2253: This is less general than in hypre, it supports only one part, and one block per process defined by a PETSc DMDA.
2254: Also, only cell-centered variables.
2256: .seealso: PCMG, MATHYPRESSTRUCT
2257: M*/
2261: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
2262: {
2264: PC_SysPFMG *ex;
2267: PetscNew(&ex); \
2268: pc->data = ex;
2270: ex->its = 1;
2271: ex->tol = 1.e-8;
2272: ex->relax_type = 1;
2273: ex->num_pre_relax = 1;
2274: ex->num_post_relax = 1;
2276: pc->ops->setfromoptions = PCSetFromOptions_SysPFMG;
2277: pc->ops->view = PCView_SysPFMG;
2278: pc->ops->destroy = PCDestroy_SysPFMG;
2279: pc->ops->apply = PCApply_SysPFMG;
2280: pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
2281: pc->ops->setup = PCSetUp_SysPFMG;
2283: MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
2284: PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2285: return(0);
2286: }