Actual source code: hyppilut.c
1: #define PETSCKSP_DLL
3: /*
5: */
7: #include private/pcimpl.h
9: #include "HYPRE.h"
10: #include "IJ_mv.h"
11: #include "parcsr_ls.h"
14: EXTERN PetscErrorCode MatHYPRE_IJMatrixCreate(Mat,HYPRE_IJMatrix*);
15: EXTERN PetscErrorCode MatHYPRE_IJMatrixCopy(Mat,HYPRE_IJMatrix);
16: EXTERN PetscErrorCode VecHYPRE_IJVectorCreate(Vec,HYPRE_IJVector*);
18: /*
19: Private context (data structure) for the preconditioner.
20: */
21: typedef struct {
22: HYPRE_Solver hsolver;
23: HYPRE_IJMatrix ij;
24: HYPRE_IJVector b,x;
26: PetscErrorCode (*destroy)(HYPRE_Solver);
27: PetscErrorCode (*solve)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
28: PetscErrorCode (*setup)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
29:
30: MPI_Comm comm_hypre;
31: char *hypre_type;
33: /* options for pilut and BoomerAMG*/
34: int maxiter;
35: double tol;
36: PetscTruth applyrichardson;
38: /* options for pilut */
39: int factorrowsize;
41: /* options for parasails */
42: int nlevels;
43: double threshhold;
44: double filter;
45: int sym;
46: double loadbal;
47: int logging;
48: int ruse;
49: int symt;
51: /* options for euclid */
52: PetscTruth bjilu;
53: int levels;
55: /* options for euclid and BoomerAMG */
56: PetscTruth printstatistics;
58: /* options for BoomerAMG */
59: int maxlevels;
60: double strongthreshold;
61: double maxrowsum;
62: int gridsweeps[4];
63: int coarsentype;
64: int measuretype;
65: int relaxtype[4];
66: double relaxweight;
67: double outerrelaxweight;
68: int relaxorder;
69: int **gridrelaxpoints;
70: double truncfactor;
71: } PC_HYPRE;
76: static PetscErrorCode PCSetUp_HYPRE(PC pc)
77: {
78: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
79: PetscErrorCode ierr;
80: HYPRE_ParCSRMatrix hmat;
81: HYPRE_ParVector bv,xv;
82: int hierr;
85: if (!jac->ij) { /* create the matrix the first time through */
86: MatHYPRE_IJMatrixCreate(pc->pmat,&jac->ij);
87: }
88: if (!jac->b) {
89: Vec vec;
90: MatGetVecs(pc->pmat,&vec,0);
91: VecHYPRE_IJVectorCreate(vec,&jac->b);
92: VecHYPRE_IJVectorCreate(vec,&jac->x);
93: VecDestroy(vec);
94: }
95: MatHYPRE_IJMatrixCopy(pc->pmat,jac->ij);
96: HYPRE_IJMatrixGetObject(jac->ij,(void**)&hmat);
97: HYPRE_IJVectorGetObject(jac->b,(void**)&bv);
98: HYPRE_IJVectorGetObject(jac->x,(void**)&xv);
99: h(*jac->setup)(jac->hsolver,hmat,bv,xv);
100: if (hierr) SETERRQ1(PETSC_ERR_LIB,"Error in HYPRE setup, error code %d",hierr);
101: return(0);
102: }
104: /*
105: Replaces the address where the HYPRE vector points to its data with the address of
106: PETSc's data. Saves the old address so it can be reset when we are finished with it.
107: Allows use to get the data into a HYPRE vector without the cost of memcopies
108: */
109: #define HYPREReplacePointer(b,newvalue,savedvalue) {\
110: hypre_ParVector *par_vector = (hypre_ParVector *)hypre_IJVectorObject(((hypre_IJVector*)b));\
111: hypre_Vector *local_vector = hypre_ParVectorLocalVector(par_vector);\
112: savedvalue = local_vector->data;\
113: local_vector->data = newvalue;}
117: static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
118: {
119: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
120: PetscErrorCode ierr;
121: HYPRE_ParCSRMatrix hmat;
122: PetscScalar *bv,*xv;
123: HYPRE_ParVector jbv,jxv;
124: PetscScalar *sbv,*sxv;
125: int hierr;
128: if (!jac->applyrichardson) {VecSet(x,0.0);}
129: VecGetArray(b,&bv);
130: VecGetArray(x,&xv);
131: HYPREReplacePointer(jac->b,bv,sbv);
132: HYPREReplacePointer(jac->x,xv,sxv);
134: HYPRE_IJMatrixGetObject(jac->ij,(void**)&hmat);
135: HYPRE_IJVectorGetObject(jac->b,(void**)&jbv);
136: HYPRE_IJVectorGetObject(jac->x,(void**)&jxv);
137: h(*jac->solve)(jac->hsolver,hmat,jbv,jxv);
138: /* error code of 1 in boomerAMG merely means convergence not achieved */
139: if (hierr && (hierr != 1 || jac->solve != HYPRE_BoomerAMGSolve)) SETERRQ1(PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
140: if (hierr) hypre__global_error = 0;
141:
142: HYPREReplacePointer(jac->b,sbv,bv);
143: HYPREReplacePointer(jac->x,sxv,xv);
144: VecRestoreArray(x,&xv);
145: VecRestoreArray(b,&bv);
146: return(0);
147: }
151: static PetscErrorCode PCDestroy_HYPRE(PC pc)
152: {
153: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
157: if (jac->ij) { HYPRE_IJMatrixDestroy(jac->ij); }
158: if (jac->b) { HYPRE_IJVectorDestroy(jac->b); }
159: if (jac->x) { HYPRE_IJVectorDestroy(jac->x); }
160: PetscStrfree(jac->hypre_type);
161: (*jac->destroy)(jac->hsolver);
162: MPI_Comm_free(&(jac->comm_hypre));
163: PetscFree(jac);
164: return(0);
165: }
167: /* --------------------------------------------------------------------------------------------*/
170: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PC pc)
171: {
172: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
174: PetscTruth flag;
177: PetscOptionsHead("HYPRE Pilut Options");
178: PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
179: if (flag) {
180: HYPRE_ParCSRPilutSetMaxIter(jac->hsolver,jac->maxiter);
181: }
182: PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
183: if (flag) {
184: HYPRE_ParCSRPilutSetDropTolerance(jac->hsolver,jac->tol);
185: }
186: PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
187: if (flag) {
188: HYPRE_ParCSRPilutSetFactorRowSize(jac->hsolver,jac->factorrowsize);
189: }
190: PetscOptionsTail();
191: return(0);
192: }
196: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
197: {
198: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
200: PetscTruth iascii;
203: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
204: if (iascii) {
205: PetscViewerASCIIPrintf(viewer," HYPRE Pilut preconditioning\n");
206: if (jac->maxiter != PETSC_DEFAULT) {
207: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: maximum number of iterations %d\n",jac->maxiter);
208: } else {
209: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default maximum number of iterations \n");
210: }
211: if (jac->tol != PETSC_DEFAULT) {
212: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: drop tolerance %G\n",jac->tol);
213: } else {
214: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default drop tolerance \n");
215: }
216: if (jac->factorrowsize != PETSC_DEFAULT) {
217: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: factor row size %d\n",jac->factorrowsize);
218: } else {
219: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default factor row size \n");
220: }
221: }
222: return(0);
223: }
225: /* --------------------------------------------------------------------------------------------*/
228: static PetscErrorCode PCSetFromOptions_HYPRE_Euclid(PC pc)
229: {
230: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
232: PetscTruth flag;
233: char *args[2];
236: PetscOptionsHead("HYPRE Euclid Options");
237: PetscOptionsInt("-pc_hypre_euclid_levels","Number of levels of fill ILU(k)","None",jac->levels,&jac->levels,&flag);
238: if (flag) {
239: char levels[16];
240: if (jac->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be nonegative",jac->levels);
241: sprintf(levels,"%d",jac->levels);
242: args[0] = (char*)"-level"; args[1] = levels;
243: HYPRE_EuclidSetParams(jac->hsolver,2,args);
244: }
245: PetscOptionsTruth("-pc_hypre_euclid_bj","Use block Jacobi ILU(k)","None",jac->bjilu,&jac->bjilu,PETSC_NULL);
246: if (jac->bjilu) {
247: args[0] =(char*) "-bj"; args[1] = (char*)"1";
248: HYPRE_EuclidSetParams(jac->hsolver,2,args);
249: }
250:
251: PetscOptionsTruth("-pc_hypre_euclid_print_statistics","Print statistics","None",jac->printstatistics,&jac->printstatistics,PETSC_NULL);
252: if (jac->printstatistics) {
253: args[0] = (char*)"-eu_stats"; args[1] = (char*)"1";
254: HYPRE_EuclidSetParams(jac->hsolver,2,args);
255: args[0] = (char*)"-eu_mem"; args[1] = (char*)"1";
256: HYPRE_EuclidSetParams(jac->hsolver,2,args);
257: }
258: PetscOptionsTail();
259: return(0);
260: }
264: static PetscErrorCode PCView_HYPRE_Euclid(PC pc,PetscViewer viewer)
265: {
266: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
268: PetscTruth iascii;
271: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
272: if (iascii) {
273: PetscViewerASCIIPrintf(viewer," HYPRE Euclid preconditioning\n");
274: PetscViewerASCIIPrintf(viewer," HYPRE Euclid: number of levels %d\n",jac->levels);
275: if (jac->bjilu) {
276: PetscViewerASCIIPrintf(viewer," HYPRE Euclid: Using block Jacobi ILU instead of parallel ILU\n");
277: }
278: }
279: return(0);
280: }
282: /* --------------------------------------------------------------------------------------------*/
284: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout"};
285: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
286: static const char *HYPREBoomerAMGRelaxType[] = {"Jacobi","sequential-Gauss-Seidel","","SOR/Jacobi","backward-SOR/Jacobi","","symmetric-SOR/Jacobi",
287: "","","Gaussian-elimination"};
290: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PC pc)
291: {
292: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
294: int n,indx;
295: PetscTruth flg, tmp_truth;
296: double tmpdbl, twodbl[2];
299: PetscOptionsHead("HYPRE BoomerAMG Options");
300: PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
301: if (flg) {
302: if (jac->maxlevels < 2) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
303: HYPRE_BoomerAMGSetMaxLevels(jac->hsolver,jac->maxlevels);
304: }
305: PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
306: if (flg) {
307: if (jac->maxiter < 1) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
308: HYPRE_BoomerAMGSetMaxIter(jac->hsolver,jac->maxiter);
309: }
310: PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call","None",jac->tol,&jac->tol,&flg);
311: if (flg) {
312: if (jac->tol < 0.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Tolerance %G must be great than or equal zero",jac->tol);
313: HYPRE_BoomerAMGSetTol(jac->hsolver,jac->tol);
314: }
316: PetscOptionsScalar("-pc_hypre_boomeramg_truncfactor","Truncation factor","None",jac->truncfactor,&jac->truncfactor,&flg);
317: if (flg) {
318: if (jac->truncfactor < 0.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Truncation factor %G must be great than or equal zero",jac->truncfactor);
319: HYPRE_BoomerAMGSetTruncFactor(jac->hsolver,jac->truncfactor);
320: }
322: PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
323: if (flg) {
324: if (jac->strongthreshold < 0.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Strong threshold %G must be great than or equal zero",jac->strongthreshold);
325: HYPRE_BoomerAMGSetStrongThreshold(jac->hsolver,jac->strongthreshold);
326: }
327: PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
328: if (flg) {
329: if (jac->maxrowsum < 0.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %G must be greater than zero",jac->maxrowsum);
330: if (jac->maxrowsum > 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %G must be less than or equal one",jac->maxrowsum);
331: HYPRE_BoomerAMGSetMaxRowSum(jac->hsolver,jac->maxrowsum);
332: }
334: /* Grid sweeps */
335: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for all grid levels (fine, up, and down)","None", jac->gridsweeps[0], &indx ,&flg);
336: if (flg) {
337: HYPRE_BoomerAMGSetNumSweeps(jac->hsolver,indx);
338: /* modify the jac structure so we can view the updated options with PC_View */
339: jac->gridsweeps[0] = indx;
340: jac->gridsweeps[1] = jac->gridsweeps[2] = jac->gridsweeps[0];
341: jac->gridsweeps[3] = 1; /*The coarse level is not affected by this function - hypre code sets to 1*/
342: }
343: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_fine","Number of sweeps for the fine level","None", jac->gridsweeps[0], &indx ,&flg);
344: if (flg) {
345: HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver,indx, 0);
346: jac->gridsweeps[0] = indx;
347: }
348: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[2], &indx ,&flg);
349: if (flg) {
350: HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver,indx, 1);
351: jac->gridsweeps[1] = indx;
352: }
353: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None", jac->gridsweeps[1], &indx ,&flg);
354: if (flg) {
355: HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver,indx, 2);
356: jac->gridsweeps[2] = indx;
357: }
358: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None", jac->gridsweeps[3], &indx ,&flg);
359: if (flg) {
360: HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver,indx, 3);
361: jac->gridsweeps[3] = indx;
362: }
364: /* Relax type */
365: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for fine, up, and down cycles (coarse level set to gaussian elimination)","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[3],&indx,&flg);
366: if (flg) {
367: jac->relaxtype[0] = jac->relaxtype[1] = jac->relaxtype[2] = indx;
368: jac->relaxtype[3] = 9; /* hypre code sets coarse grid to 9 (G.E.)*/
369: hypre_BoomerAMGSetRelaxType(jac->hsolver, indx);
370: }
371: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_fine","Relax type on fine grid","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[3],&indx,&flg);
372: if (flg) {
373: jac->relaxtype[0] = indx;
374: HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, indx, 0);
375: }
376: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[3],&indx,&flg);
377: if (flg) {
378: jac->relaxtype[1] = indx;
379: HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, indx, 1);
380: }
381: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[3],&indx,&flg);
382: if (flg) {
383: jac->relaxtype[2] = indx;
384: HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, indx, 2);
385: }
386: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[9],&indx,&flg);
387: if (flg) {
388: jac->relaxtype[3] = indx;
389: HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, indx, 3);
390: }
392: /* Relaxation Weight */
393: 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);
394: if (flg) {
395: hypre_BoomerAMGSetRelaxWt( jac->hsolver, tmpdbl);
396: jac->relaxweight = tmpdbl;
397: }
399: n=2;
400: twodbl[0] = twodbl[1] = 1.0;
401: PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
402: if (flg) {
403: if (n == 2) {
404: indx = (int)PetscAbsReal(twodbl[1]);
405: hypre_BoomerAMGSetLevelRelaxWt( jac->hsolver, twodbl[0], indx);
406: } else {
407: SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Relax weight level: you must provide 2 values separated by a comma (and no space), you provided %d",n);
408: }
409: }
411: /* Outer relaxation Weight */
412: 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);
413: if (flg) {
414: hypre_BoomerAMGSetOuterWt( jac->hsolver, tmpdbl);
415: jac->outerrelaxweight = tmpdbl;
416: }
418: n=2;
419: twodbl[0] = twodbl[1] = 1.0;
420: PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
421: if (flg) {
422: if (n == 2) {
423: indx = (int)PetscAbsReal(twodbl[1]);
424: hypre_BoomerAMGSetLevelOuterWt( jac->hsolver, twodbl[0], indx);
425: } else {
426: SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Relax weight outer level: You must provide 2 values separated by a comma (and no space), you provided %d",n);
427: }
428: }
430: /* the Relax Order */
431: /* PetscOptionsName("-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", &flg); */
432: PetscOptionsTruth( "-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", PETSC_FALSE, &tmp_truth, &flg);
434: if (flg) {
435: jac->relaxorder = 0;
436: hypre_BoomerAMGSetRelaxOrder( jac->hsolver, jac->relaxorder);
437: }
438: PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,2,HYPREBoomerAMGMeasureType[0],&indx,&flg);
439: if (flg) {
440: jac->measuretype = indx;
441: HYPRE_BoomerAMGSetMeasureType(jac->hsolver,jac->measuretype);
442: }
443: PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,7,HYPREBoomerAMGCoarsenType[6],&indx,&flg);
444: if (flg) {
445: jac->coarsentype = indx;
446: HYPRE_BoomerAMGSetCoarsenType(jac->hsolver,jac->coarsentype);
447: }
448: PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
449: if (flg) {
450: int level=3;
451: jac->printstatistics = PETSC_TRUE;
452: PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,PETSC_NULL);
453: HYPRE_BoomerAMGSetPrintLevel(jac->hsolver,level);
454: HYPRE_BoomerAMGSetDebugFlag(jac->hsolver,level);
455: }
456: PetscOptionsTail();
457: return(0);
458: }
462: static PetscErrorCode PCApplyRichardson_BoomerAMG(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its)
463: {
464: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
468: HYPRE_BoomerAMGSetMaxIter(jac->hsolver,its*jac->maxiter);
469: HYPRE_BoomerAMGSetTol(jac->hsolver,PetscMin(rtol,jac->tol));
470: jac->applyrichardson = PETSC_TRUE;
471: PCApply_HYPRE(pc,b,y);
472: jac->applyrichardson = PETSC_FALSE;
473: HYPRE_BoomerAMGSetTol(jac->hsolver,jac->tol);
474: HYPRE_BoomerAMGSetMaxIter(jac->hsolver,jac->maxiter);
475: return(0);
476: }
481: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
482: {
483: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
485: PetscTruth iascii;
488: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
489: if (iascii) {
490: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG preconditioning\n");
491: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum number of levels %d\n",jac->maxlevels);
492: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum number of iterations PER hypre call%d\n",jac->maxiter);
493: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Convergence tolerance PER hypre call%G\n",jac->tol);
494: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Threshold for strong coupling %G\n",jac->strongthreshold);
495: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum row sums %G\n",jac->maxrowsum);
497: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps on fine grid %d\n",jac->gridsweeps[0]);
498: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps down %d\n",jac->gridsweeps[1]);
499: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps up %d\n",jac->gridsweeps[2]);
500: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps on coarse %d\n",jac->gridsweeps[3]);
502: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax on fine grid %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[0]]);
503: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax down %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[1]]);
504: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax up %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[2]]);
505: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax on coarse %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[3]]);
507: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax weight (all) %G\n",jac->relaxweight);
508: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Outer relax weight (all) %G\n",jac->outerrelaxweight);
510: if (jac->relaxorder) {
511: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using CF-relaxation\n");
512: } else {
513: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Not using CF-relaxation\n");
514: }
515: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Measure type %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
516: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Coarsen type %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
517: }
518: return(0);
519: }
521: /* --------------------------------------------------------------------------------------------*/
524: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PC pc)
525: {
526: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
528: int indx;
529: PetscTruth flag;
530: const char *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};
533: PetscOptionsHead("HYPRE ParaSails Options");
534: PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);
535: PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshhold,&jac->threshhold,&flag);
536: if (flag) {
537: HYPRE_ParaSailsSetParams(jac->hsolver,jac->threshhold,jac->nlevels);
538: }
540: PetscOptionsReal("-pc_hypre_parasails_filter","filter","None",jac->filter,&jac->filter,&flag);
541: if (flag) {
542: HYPRE_ParaSailsSetFilter(jac->hsolver,jac->filter);
543: }
545: PetscOptionsReal("-pc_hypre_parasails_loadbal","Load balance","None",jac->loadbal,&jac->loadbal,&flag);
546: if (flag) {
547: HYPRE_ParaSailsSetLoadbal(jac->hsolver,jac->loadbal);
548: }
550: PetscOptionsTruth("-pc_hypre_parasails_logging","Print info to screen","None",(PetscTruth)jac->logging,(PetscTruth*)&jac->logging,&flag);
551: if (flag) {
552: HYPRE_ParaSailsSetLogging(jac->hsolver,jac->logging);
553: }
555: PetscOptionsTruth("-pc_hypre_parasails_reuse","Reuse nonzero pattern in preconditioner","None",(PetscTruth)jac->ruse,(PetscTruth*)&jac->ruse,&flag);
556: if (flag) {
557: HYPRE_ParaSailsSetReuse(jac->hsolver,jac->ruse);
558: }
560: PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,3,symtlist[0],&indx,&flag);
561: if (flag) {
562: jac->symt = indx;
563: HYPRE_ParaSailsSetSym(jac->hsolver,jac->symt);
564: }
566: PetscOptionsTail();
567: return(0);
568: }
572: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
573: {
574: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
576: PetscTruth iascii;
577: const char *symt = 0;;
580: PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
581: if (iascii) {
582: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails preconditioning\n");
583: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: nlevels %d\n",jac->nlevels);
584: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: threshold %G\n",jac->threshhold);
585: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: filter %G\n",jac->filter);
586: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: load balance %G\n",jac->loadbal);
587: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: reuse nonzero structure %s\n",PetscTruths[jac->ruse]);
588: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: print info to screen %s\n",PetscTruths[jac->logging]);
589: if (!jac->symt) {
590: symt = "nonsymmetric matrix and preconditioner";
591: } else if (jac->symt == 1) {
592: symt = "SPD matrix and preconditioner";
593: } else if (jac->symt == 2) {
594: symt = "nonsymmetric matrix but SPD preconditioner";
595: } else {
596: SETERRQ1(PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
597: }
598: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: %s\n",symt);
599: }
600: return(0);
601: }
602: /* ---------------------------------------------------------------------------------*/
607: PetscErrorCode PCHYPREGetType_HYPRE(PC pc,const char *name[])
608: {
609: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
612: *name = jac->hypre_type;
613: return(0);
614: }
620: PetscErrorCode PCHYPRESetType_HYPRE(PC pc,const char name[])
621: {
622: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
624: PetscTruth flag;
625: PetscInt bs;
628: if (pc->ops->setup) {
629: SETERRQ(PETSC_ERR_ORDER,"Cannot set the HYPRE preconditioner type once it has been set");
630: }
632: pc->ops->setup = PCSetUp_HYPRE;
633: pc->ops->apply = PCApply_HYPRE;
634: pc->ops->destroy = PCDestroy_HYPRE;
636: jac->maxiter = PETSC_DEFAULT;
637: jac->tol = PETSC_DEFAULT;
638: jac->printstatistics = PetscLogPrintInfo;
640: PetscStrallocpy(name, &jac->hypre_type);
641: PetscStrcmp("pilut",name,&flag);
642: if (flag) {
643: HYPRE_ParCSRPilutCreate(jac->comm_hypre,&jac->hsolver);
644: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
645: pc->ops->view = PCView_HYPRE_Pilut;
646: jac->destroy = HYPRE_ParCSRPilutDestroy;
647: jac->setup = HYPRE_ParCSRPilutSetup;
648: jac->solve = HYPRE_ParCSRPilutSolve;
649: jac->factorrowsize = PETSC_DEFAULT;
650: return(0);
651: }
652: PetscStrcmp("parasails",name,&flag);
653: if (flag) {
654: HYPRE_ParaSailsCreate(jac->comm_hypre,&jac->hsolver);
655: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
656: pc->ops->view = PCView_HYPRE_ParaSails;
657: jac->destroy = HYPRE_ParaSailsDestroy;
658: jac->setup = HYPRE_ParaSailsSetup;
659: jac->solve = HYPRE_ParaSailsSolve;
660: /* initialize */
661: jac->nlevels = 1;
662: jac->threshhold = .1;
663: jac->filter = .1;
664: jac->loadbal = 0;
665: if (PetscLogPrintInfo) {
666: jac->logging = (int) PETSC_TRUE;
667: } else {
668: jac->logging = (int) PETSC_FALSE;
669: }
670: jac->ruse = (int) PETSC_FALSE;
671: jac->symt = 0;
672: HYPRE_ParaSailsSetParams(jac->hsolver,jac->threshhold,jac->nlevels);
673: HYPRE_ParaSailsSetFilter(jac->hsolver,jac->filter);
674: HYPRE_ParaSailsSetLoadbal(jac->hsolver,jac->loadbal);
675: HYPRE_ParaSailsSetLogging(jac->hsolver,jac->logging);
676: HYPRE_ParaSailsSetReuse(jac->hsolver,jac->ruse);
677: HYPRE_ParaSailsSetSym(jac->hsolver,jac->symt);
678: return(0);
679: }
680: PetscStrcmp("euclid",name,&flag);
681: if (flag) {
682: HYPRE_EuclidCreate(jac->comm_hypre,&jac->hsolver);
683: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
684: pc->ops->view = PCView_HYPRE_Euclid;
685: jac->destroy = HYPRE_EuclidDestroy;
686: jac->setup = HYPRE_EuclidSetup;
687: jac->solve = HYPRE_EuclidSolve;
688: /* initialization */
689: jac->bjilu = PETSC_FALSE;
690: jac->levels = 1;
691: return(0);
692: }
693: PetscStrcmp("boomeramg",name,&flag);
694: if (flag) {
695: HYPRE_BoomerAMGCreate(&jac->hsolver);
696: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_BoomerAMG;
697: pc->ops->view = PCView_HYPRE_BoomerAMG;
698: jac->destroy = HYPRE_BoomerAMGDestroy;
699: jac->setup = HYPRE_BoomerAMGSetup;
700: jac->solve = HYPRE_BoomerAMGSolve;
701: pc->ops->applyrichardson = PCApplyRichardson_BoomerAMG;
702: /* these defaults match the hypre defaults */
703: jac->maxlevels = 25;
704: jac->maxiter = 1;
705: jac->tol = 1.e-7;
706: jac->strongthreshold = .25;
707: jac->maxrowsum = .9;
708: jac->coarsentype = 6;
709: jac->measuretype = 0;
710: jac->applyrichardson = PETSC_FALSE;
711: jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = jac->gridsweeps[3] = 1;
712: jac->relaxtype[0] = jac->relaxtype[1] = jac->relaxtype[2] = 3;
713: jac->relaxtype[3] = 9;
714: jac->relaxweight = 1.0;
715: jac->outerrelaxweight = 1.0;
716: jac->relaxorder = 1;
717: HYPRE_BoomerAMGSetMaxLevels(jac->hsolver,jac->maxlevels);
718: HYPRE_BoomerAMGSetMaxIter(jac->hsolver,jac->maxiter);
719: HYPRE_BoomerAMGSetTol(jac->hsolver,jac->tol);
720: HYPRE_BoomerAMGSetTruncFactor(jac->hsolver,jac->truncfactor);
721: HYPRE_BoomerAMGSetStrongThreshold(jac->hsolver,jac->strongthreshold);
722: HYPRE_BoomerAMGSetMaxRowSum(jac->hsolver,jac->maxrowsum);
723: HYPRE_BoomerAMGSetMeasureType(jac->hsolver,jac->measuretype);
724: HYPRE_BoomerAMGSetCoarsenType(jac->hsolver,jac->coarsentype);
725: MatGetBlockSize(pc->pmat,&bs);
726: HYPRE_BoomerAMGSetNumFunctions(jac->hsolver,bs);
727: return(0);
728: }
729: SETERRQ1(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are pilut, parasails, euclid, boomeramg",name);
730: return(0);
731: }
734: /*
735: It only gets here if the HYPRE type has not been set before the call to
736: ...SetFromOptions() which actually is most of the time
737: */
740: static PetscErrorCode PCSetFromOptions_HYPRE(PC pc)
741: {
743: int indx;
744: const char *type[] = {"pilut","parasails","boomeramg","euclid"};
745: PetscTruth flg;
748: PetscOptionsHead("HYPRE preconditioner options");
749: PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"pilut",&indx,&flg);
750: if (PetscOptionsPublishCount) { /* force the default if it was not yet set and user did not set with option */
751: if (!flg && !pc->ops->apply) {
752: flg = PETSC_TRUE;
753: indx = 0;
754: }
755: }
756: if (flg) {
757: PCHYPRESetType_HYPRE(pc,type[indx]);
758: }
759: if (pc->ops->setfromoptions) {
760: pc->ops->setfromoptions(pc);
761: }
762: PetscOptionsTail();
763: return(0);
764: }
768: /*@C
769: PCHYPRESetType - Sets which hypre preconditioner you wish to use
771: Input Parameters:
772: + pc - the preconditioner context
773: - name - either pilut, parasails, boomeramg, euclid
775: Options Database Keys:
776: -pc_hypre_type - One of pilut, parasails, boomeramg, euclid
777:
778: Level: intermediate
780: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
781: PCHYPRE
783: @*/
784: PetscErrorCode PCHYPRESetType(PC pc,const char name[])
785: {
786: PetscErrorCode ierr,(*f)(PC,const char[]);
791: PetscObjectQueryFunction((PetscObject)pc,"PCHYPRESetType_C",(void (**)(void))&f);
792: if (f) {
793: (*f)(pc,name);
794: }
795: return(0);
796: }
800: /*@C
801: PCHYPREGetType - Gets which hypre preconditioner you are using
803: Input Parameter:
804: . pc - the preconditioner context
806: Output Parameter:
807: . name - either pilut, parasails, boomeramg, euclid
809: Level: intermediate
811: .seealso: PCCreate(), PCHYPRESetType(), PCType (for list of available types), PC,
812: PCHYPRE
814: @*/
815: PetscErrorCode PCHYPREGetType(PC pc,const char *name[])
816: {
817: PetscErrorCode ierr,(*f)(PC,const char*[]);
822: PetscObjectQueryFunction((PetscObject)pc,"PCHYPREGetType_C",(void (**)(void))&f);
823: if (f) {
824: (*f)(pc,name);
825: }
826: return(0);
827: }
829: /*MC
830: PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre
832: Options Database Keys:
833: + -pc_hypre_type - One of pilut, parasails, boomeramg, euclid
834: - Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX
835: preconditioner
836:
837: Level: intermediate
839: Notes: Apart from pc_hypre_type (for which there is PCHYPRESetType()),
840: the many hypre options can ONLY be set via the options database (e.g. the command line
841: or with PetscOptionsSetValue(), there are no functions to set them)
843: The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_rtol refer to the number of iterations
844: (V-cycles) that boomeramg does EACH time it is called. So for example, if it is set to 2 then
845: 2-V-cycles are being used to define the preconditioner. -ksp_max_iter and -ksp_rtol STILL determine
846: the total number of iterations and tolerance for the Krylov solver. For example, if
847: -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10 then AT MOST twenty V-cycles of boomeramg
848: will be called.
851: If you wish to use boomerAMG WITHOUT a Krylov method use -ksp_type richardson NOT -ksp_type preonly
852: and use -ksp_max_it to control the number of V-cycles.
853: (see the PETSc FAQ.html at the PETSc website under the Documentation tab).
855: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
856: PCHYPRESetType()
858: M*/
863: PetscErrorCode PCCreate_HYPRE(PC pc)
864: {
865: PC_HYPRE *jac;
869: PetscNew(PC_HYPRE,&jac);
870: pc->data = jac;
871: pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
872: jac->hypre_type = PETSC_NULL;
873: /* Com_dup for hypre */
874: MPI_Comm_dup(pc->comm,&(jac->comm_hypre));
875: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCHYPRESetType_C","PCHYPRESetType_HYPRE",PCHYPRESetType_HYPRE);
876: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCHYPREGetType_C","PCHYPREGetType_HYPRE",PCHYPREGetType_HYPRE);
877: return(0);
878: }