Actual source code: hypre.c
petsc-3.4.5 2014-06-29
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>
11: /*
12: Private context (data structure) for the preconditioner.
13: */
14: typedef struct {
15: HYPRE_Solver hsolver;
16: HYPRE_IJMatrix ij;
17: HYPRE_IJVector b,x;
19: HYPRE_Int (*destroy)(HYPRE_Solver);
20: HYPRE_Int (*solve)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
21: HYPRE_Int (*setup)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
23: MPI_Comm comm_hypre;
24: char *hypre_type;
26: /* options for Pilut and BoomerAMG*/
27: PetscInt maxiter;
28: double tol;
30: /* options for Pilut */
31: PetscInt factorrowsize;
33: /* options for ParaSails */
34: PetscInt nlevels;
35: double threshhold;
36: double filter;
37: PetscInt sym;
38: double loadbal;
39: PetscInt logging;
40: PetscInt ruse;
41: PetscInt symt;
43: /* options for Euclid */
44: PetscBool bjilu;
45: PetscInt levels;
47: /* options for Euclid and BoomerAMG */
48: PetscBool printstatistics;
50: /* options for BoomerAMG */
51: PetscInt cycletype;
52: PetscInt maxlevels;
53: double strongthreshold;
54: double maxrowsum;
55: PetscInt gridsweeps[3];
56: PetscInt coarsentype;
57: PetscInt measuretype;
58: PetscInt relaxtype[3];
59: double relaxweight;
60: double outerrelaxweight;
61: PetscInt relaxorder;
62: double truncfactor;
63: PetscBool applyrichardson;
64: PetscInt pmax;
65: PetscInt interptype;
66: PetscInt agg_nl;
67: PetscInt agg_num_paths;
68: PetscInt nodal_coarsen;
69: PetscBool nodal_relax;
70: PetscInt nodal_relax_levels;
71: } PC_HYPRE;
75: PetscErrorCode PCHYPREGetSolver(PC pc,HYPRE_Solver *hsolver)
76: {
77: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
80: *hsolver = jac->hsolver;
81: return(0);
82: }
86: static PetscErrorCode PCSetUp_HYPRE(PC pc)
87: {
88: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
89: PetscErrorCode ierr;
90: HYPRE_ParCSRMatrix hmat;
91: HYPRE_ParVector bv,xv;
92: PetscInt bs;
95: if (!jac->hypre_type) {
96: PCHYPRESetType(pc,"boomeramg");
97: }
99: if (pc->setupcalled) {
100: /* always destroy the old matrix and create a new memory;
101: hope this does not churn the memory too much. The problem
102: is I do not know if it is possible to put the matrix back to
103: its initial state so that we can directly copy the values
104: the second time through. */
105: PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij));
106: jac->ij = 0;
107: }
109: if (!jac->ij) { /* create the matrix the first time through */
110: MatHYPRE_IJMatrixCreate(pc->pmat,&jac->ij);
111: }
112: if (!jac->b) { /* create the vectors the first time through */
113: Vec x,b;
114: MatGetVecs(pc->pmat,&x,&b);
115: VecHYPRE_IJVectorCreate(x,&jac->x);
116: VecHYPRE_IJVectorCreate(b,&jac->b);
117: VecDestroy(&x);
118: VecDestroy(&b);
119: }
121: /* special case for BoomerAMG */
122: if (jac->setup == HYPRE_BoomerAMGSetup) {
123: MatGetBlockSize(pc->pmat,&bs);
124: if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
125: };
126: MatHYPRE_IJMatrixCopy(pc->pmat,jac->ij);
127: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
128: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&bv));
129: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&xv));
130: PetscStackCall("HYPRE_SetupXXX",(*jac->setup)(jac->hsolver,hmat,bv,xv););
131: return(0);
132: }
134: /*
135: Replaces the address where the HYPRE vector points to its data with the address of
136: PETSc's data. Saves the old address so it can be reset when we are finished with it.
137: Allows use to get the data into a HYPRE vector without the cost of memcopies
138: */
139: #define HYPREReplacePointer(b,newvalue,savedvalue) { \
140: hypre_ParVector *par_vector = (hypre_ParVector*)hypre_IJVectorObject(((hypre_IJVector*)b)); \
141: hypre_Vector *local_vector = hypre_ParVectorLocalVector(par_vector); \
142: savedvalue = local_vector->data; \
143: local_vector->data = newvalue; \
144: }
148: static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
149: {
150: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
151: PetscErrorCode ierr;
152: HYPRE_ParCSRMatrix hmat;
153: PetscScalar *bv,*xv;
154: HYPRE_ParVector jbv,jxv;
155: PetscScalar *sbv,*sxv;
156: PetscInt hierr;
159: if (!jac->applyrichardson) {VecSet(x,0.0);}
160: VecGetArray(b,&bv);
161: VecGetArray(x,&xv);
162: HYPREReplacePointer(jac->b,bv,sbv);
163: HYPREReplacePointer(jac->x,xv,sxv);
165: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
166: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
167: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));
168: PetscStackCall("Hypre solve",h(*jac->solve)(jac->hsolver,hmat,jbv,jxv);
169: if (hierr && hierr != HYPRE_ERROR_CONV) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
170: if (hierr) hypre__global_error = 0;);
172: HYPREReplacePointer(jac->b,sbv,bv);
173: HYPREReplacePointer(jac->x,sxv,xv);
174: VecRestoreArray(x,&xv);
175: VecRestoreArray(b,&bv);
176: return(0);
177: }
181: static PetscErrorCode PCDestroy_HYPRE(PC pc)
182: {
183: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
187: if (jac->ij) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij));
188: if (jac->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->b));
189: if (jac->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->x));
190: if (jac->destroy) PetscStackCall("HYPRE_DistroyXXX",(*jac->destroy)(jac->hsolver););
191: PetscFree(jac->hypre_type);
192: if (jac->comm_hypre != MPI_COMM_NULL) { MPI_Comm_free(&(jac->comm_hypre));}
193: PetscFree(pc->data);
195: PetscObjectChangeTypeName((PetscObject)pc,0);
196: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);
197: PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);
198: return(0);
199: }
201: /* --------------------------------------------------------------------------------------------*/
204: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PC pc)
205: {
206: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
208: PetscBool flag;
211: PetscOptionsHead("HYPRE Pilut Options");
212: PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
213: if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter));
214: PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
215: if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol));
216: PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
217: if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize));
218: PetscOptionsTail();
219: return(0);
220: }
224: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
225: {
226: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
228: PetscBool iascii;
231: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
232: if (iascii) {
233: PetscViewerASCIIPrintf(viewer," HYPRE Pilut preconditioning\n");
234: if (jac->maxiter != PETSC_DEFAULT) {
235: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: maximum number of iterations %d\n",jac->maxiter);
236: } else {
237: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default maximum number of iterations \n");
238: }
239: if (jac->tol != PETSC_DEFAULT) {
240: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: drop tolerance %G\n",jac->tol);
241: } else {
242: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default drop tolerance \n");
243: }
244: if (jac->factorrowsize != PETSC_DEFAULT) {
245: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: factor row size %d\n",jac->factorrowsize);
246: } else {
247: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default factor row size \n");
248: }
249: }
250: return(0);
251: }
253: /* --------------------------------------------------------------------------------------------*/
256: static PetscErrorCode PCSetFromOptions_HYPRE_Euclid(PC pc)
257: {
258: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
260: PetscBool flag;
261: char *args[8],levels[16];
262: PetscInt cnt = 0;
265: PetscOptionsHead("HYPRE Euclid Options");
266: PetscOptionsInt("-pc_hypre_euclid_levels","Number of levels of fill ILU(k)","None",jac->levels,&jac->levels,&flag);
267: if (flag) {
268: if (jac->levels < 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be nonegative",jac->levels);
269: PetscSNPrintf(levels,sizeof(levels),"%D",jac->levels);
270: args[cnt++] = (char*)"-level"; args[cnt++] = levels;
271: }
272: PetscOptionsBool("-pc_hypre_euclid_bj","Use block Jacobi ILU(k)","None",jac->bjilu,&jac->bjilu,NULL);
273: if (jac->bjilu) {
274: args[cnt++] =(char*) "-bj"; args[cnt++] = (char*)"1";
275: }
277: PetscOptionsBool("-pc_hypre_euclid_print_statistics","Print statistics","None",jac->printstatistics,&jac->printstatistics,NULL);
278: if (jac->printstatistics) {
279: args[cnt++] = (char*)"-eu_stats"; args[cnt++] = (char*)"1";
280: args[cnt++] = (char*)"-eu_mem"; args[cnt++] = (char*)"1";
281: }
282: PetscOptionsTail();
283: if (cnt) PetscStackCallStandard(HYPRE_EuclidSetParams,(jac->hsolver,cnt,args));
284: return(0);
285: }
289: static PetscErrorCode PCView_HYPRE_Euclid(PC pc,PetscViewer viewer)
290: {
291: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
293: PetscBool iascii;
296: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
297: if (iascii) {
298: PetscViewerASCIIPrintf(viewer," HYPRE Euclid preconditioning\n");
299: PetscViewerASCIIPrintf(viewer," HYPRE Euclid: number of levels %d\n",jac->levels);
300: if (jac->bjilu) {
301: PetscViewerASCIIPrintf(viewer," HYPRE Euclid: Using block Jacobi ILU instead of parallel ILU\n");
302: }
303: }
304: return(0);
305: }
307: /* --------------------------------------------------------------------------------------------*/
311: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
312: {
313: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
314: PetscErrorCode ierr;
315: HYPRE_ParCSRMatrix hmat;
316: PetscScalar *bv,*xv;
317: HYPRE_ParVector jbv,jxv;
318: PetscScalar *sbv,*sxv;
319: PetscInt hierr;
322: VecSet(x,0.0);
323: VecGetArray(b,&bv);
324: VecGetArray(x,&xv);
325: HYPREReplacePointer(jac->b,bv,sbv);
326: HYPREReplacePointer(jac->x,xv,sxv);
328: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
329: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
330: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));
332: hHYPRE_BoomerAMGSolveT(jac->hsolver,hmat,jbv,jxv);
333: /* error code of 1 in BoomerAMG merely means convergence not achieved */
334: if (hierr && (hierr != 1)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
335: if (hierr) hypre__global_error = 0;
337: HYPREReplacePointer(jac->b,sbv,bv);
338: HYPREReplacePointer(jac->x,sxv,xv);
339: VecRestoreArray(x,&xv);
340: VecRestoreArray(b,&bv);
341: return(0);
342: }
344: /* static array length */
345: #define ALEN(a) (sizeof(a)/sizeof((a)[0]))
347: static const char *HYPREBoomerAMGCycleType[] = {"","V","W"};
348: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
349: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
350: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
351: static const char *HYPREBoomerAMGRelaxType[] = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi",
352: "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi",
353: "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination",
354: "" /* 10 */, "" /* 11 */, "" /* 12 */, "" /* 13 */, "" /* 14 */,
355: "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"};
356: static const char *HYPREBoomerAMGInterpType[] = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
357: "ext+i-cc", "standard", "standard-wts", "", "", "FF", "FF1"};
360: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PC pc)
361: {
362: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
364: PetscInt n,indx,level;
365: PetscBool flg, tmp_truth;
366: double tmpdbl, twodbl[2];
369: PetscOptionsHead("HYPRE BoomerAMG Options");
370: PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);
371: if (flg) {
372: jac->cycletype = indx+1;
373: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
374: }
375: PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
376: if (flg) {
377: if (jac->maxlevels < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
378: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
379: }
380: PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
381: if (flg) {
382: if (jac->maxiter < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
383: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
384: }
385: PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);
386: if (flg) {
387: if (jac->tol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance %G must be greater than or equal to zero",jac->tol);
388: PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
389: }
391: PetscOptionsScalar("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);
392: if (flg) {
393: if (jac->truncfactor < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Truncation factor %G must be great than or equal zero",jac->truncfactor);
394: PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
395: }
397: PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);
398: if (flg) {
399: if (jac->pmax < 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"P_max %G must be greater than or equal to zero",jac->pmax);
400: PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
401: }
403: PetscOptionsInt("-pc_hypre_boomeramg_agg_nl","Number of levels of aggressive coarsening","None",jac->agg_nl,&jac->agg_nl,&flg);
404: if (flg) {
405: 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",jac->agg_nl);
407: PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
408: }
411: PetscOptionsInt("-pc_hypre_boomeramg_agg_num_paths","Number of paths for aggressive coarsening","None",jac->agg_num_paths,&jac->agg_num_paths,&flg);
412: if (flg) {
413: 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",jac->agg_num_paths);
415: PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
416: }
419: PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
420: if (flg) {
421: if (jac->strongthreshold < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Strong threshold %G must be great than or equal zero",jac->strongthreshold);
422: PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
423: }
424: PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
425: if (flg) {
426: if (jac->maxrowsum < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %G must be greater than zero",jac->maxrowsum);
427: if (jac->maxrowsum > 1.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %G must be less than or equal one",jac->maxrowsum);
428: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
429: }
431: /* Grid sweeps */
432: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);
433: if (flg) {
434: PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx));
435: /* modify the jac structure so we can view the updated options with PC_View */
436: jac->gridsweeps[0] = indx;
437: jac->gridsweeps[1] = indx;
438: /*defaults coarse to 1 */
439: jac->gridsweeps[2] = 1;
440: }
442: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);
443: if (flg) {
444: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1));
445: jac->gridsweeps[0] = indx;
446: }
447: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);
448: if (flg) {
449: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2));
450: jac->gridsweeps[1] = indx;
451: }
452: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);
453: if (flg) {
454: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3));
455: jac->gridsweeps[2] = indx;
456: }
458: /* Relax type */
459: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
460: if (flg) {
461: jac->relaxtype[0] = jac->relaxtype[1] = indx;
462: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx));
463: /* by default, coarse type set to 9 */
464: jac->relaxtype[2] = 9;
466: }
467: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
468: if (flg) {
469: jac->relaxtype[0] = indx;
470: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1));
471: }
472: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
473: if (flg) {
474: jac->relaxtype[1] = indx;
475: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2));
476: }
477: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);
478: if (flg) {
479: jac->relaxtype[2] = indx;
480: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3));
481: }
483: /* Relaxation Weight */
484: 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);
485: if (flg) {
486: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl));
487: jac->relaxweight = tmpdbl;
488: }
490: n = 2;
491: twodbl[0] = twodbl[1] = 1.0;
492: PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
493: if (flg) {
494: if (n == 2) {
495: indx = (int)PetscAbsReal(twodbl[1]);
496: PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx));
497: } 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);
498: }
500: /* Outer relaxation Weight */
501: 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);
502: if (flg) {
503: PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl));
504: jac->outerrelaxweight = tmpdbl;
505: }
507: n = 2;
508: twodbl[0] = twodbl[1] = 1.0;
509: PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
510: if (flg) {
511: if (n == 2) {
512: indx = (int)PetscAbsReal(twodbl[1]);
513: PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx));
514: } 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);
515: }
517: /* the Relax Order */
518: PetscOptionsBool("-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", PETSC_FALSE, &tmp_truth, &flg);
520: if (flg) {
521: jac->relaxorder = 0;
522: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
523: }
524: PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);
525: if (flg) {
526: jac->measuretype = indx;
527: PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
528: }
529: /* update list length 3/07 */
530: PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);
531: if (flg) {
532: jac->coarsentype = indx;
533: PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
534: }
536: /* new 3/07 */
537: PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);
538: if (flg) {
539: jac->interptype = indx;
540: PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
541: }
543: PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
544: if (flg) {
545: level = 3;
546: PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);
548: jac->printstatistics = PETSC_TRUE;
549: PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level));
550: }
552: PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);
553: if (flg) {
554: level = 3;
555: PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);
557: jac->printstatistics = PETSC_TRUE;
558: PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level));
559: }
561: PetscOptionsBool("-pc_hypre_boomeramg_nodal_coarsen", "HYPRE_BoomerAMGSetNodal()", "None", PETSC_FALSE, &tmp_truth, &flg);
562: if (flg && tmp_truth) {
563: jac->nodal_coarsen = 1;
564: PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,1));
565: }
567: PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);
568: if (flg && tmp_truth) {
569: PetscInt tmp_int;
570: PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);
571: if (flg) jac->nodal_relax_levels = tmp_int;
572: PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6));
573: PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1));
574: PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0));
575: PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels));
576: }
578: PetscOptionsTail();
579: return(0);
580: }
584: 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)
585: {
586: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
588: PetscInt oits;
591: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter));
592: PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol));
593: jac->applyrichardson = PETSC_TRUE;
594: PCApply_HYPRE(pc,b,y);
595: jac->applyrichardson = PETSC_FALSE;
596: PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,&oits));
597: *outits = oits;
598: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
599: else *reason = PCRICHARDSON_CONVERGED_RTOL;
600: PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
601: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
602: return(0);
603: }
608: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
609: {
610: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
612: PetscBool iascii;
615: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
616: if (iascii) {
617: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG preconditioning\n");
618: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);
619: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum number of levels %d\n",jac->maxlevels);
620: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum number of iterations PER hypre call %d\n",jac->maxiter);
621: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Convergence tolerance PER hypre call %G\n",jac->tol);
622: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Threshold for strong coupling %G\n",jac->strongthreshold);
623: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Interpolation truncation factor %G\n",jac->truncfactor);
624: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Interpolation: max elements per row %d\n",jac->pmax);
625: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Number of levels of aggressive coarsening %d\n",jac->agg_nl);
626: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Number of paths for aggressive coarsening %d\n",jac->agg_num_paths);
628: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum row sums %G\n",jac->maxrowsum);
630: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps down %d\n",jac->gridsweeps[0]);
631: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps up %d\n",jac->gridsweeps[1]);
632: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps on coarse %d\n",jac->gridsweeps[2]);
634: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax down %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[0]]);
635: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax up %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[1]]);
636: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax on coarse %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[2]]);
638: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax weight (all) %G\n",jac->relaxweight);
639: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Outer relax weight (all) %G\n",jac->outerrelaxweight);
641: if (jac->relaxorder) {
642: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using CF-relaxation\n");
643: } else {
644: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Not using CF-relaxation\n");
645: }
646: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Measure type %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
647: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Coarsen type %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
648: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Interpolation type %s\n",HYPREBoomerAMGInterpType[jac->interptype]);
649: if (jac->nodal_coarsen) {
650: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal coarsening (with HYPRE_BOOMERAMGSetNodal())\n");
651: }
652: if (jac->nodal_relax) {
653: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal relaxation via Schwarz smoothing on levels %d\n",jac->nodal_relax_levels);
654: }
655: }
656: return(0);
657: }
659: /* --------------------------------------------------------------------------------------------*/
662: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PC pc)
663: {
664: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
666: PetscInt indx;
667: PetscBool flag;
668: const char *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};
671: PetscOptionsHead("HYPRE ParaSails Options");
672: PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);
673: PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshhold,&jac->threshhold,&flag);
674: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
676: PetscOptionsReal("-pc_hypre_parasails_filter","filter","None",jac->filter,&jac->filter,&flag);
677: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
679: PetscOptionsReal("-pc_hypre_parasails_loadbal","Load balance","None",jac->loadbal,&jac->loadbal,&flag);
680: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
682: PetscOptionsBool("-pc_hypre_parasails_logging","Print info to screen","None",(PetscBool)jac->logging,(PetscBool*)&jac->logging,&flag);
683: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
685: PetscOptionsBool("-pc_hypre_parasails_reuse","Reuse nonzero pattern in preconditioner","None",(PetscBool)jac->ruse,(PetscBool*)&jac->ruse,&flag);
686: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
688: PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);
689: if (flag) {
690: jac->symt = indx;
691: PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
692: }
694: PetscOptionsTail();
695: return(0);
696: }
700: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
701: {
702: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
704: PetscBool iascii;
705: const char *symt = 0;;
708: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
709: if (iascii) {
710: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails preconditioning\n");
711: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: nlevels %d\n",jac->nlevels);
712: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: threshold %G\n",jac->threshhold);
713: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: filter %G\n",jac->filter);
714: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: load balance %G\n",jac->loadbal);
715: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: reuse nonzero structure %s\n",PetscBools[jac->ruse]);
716: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: print info to screen %s\n",PetscBools[jac->logging]);
717: if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
718: else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
719: else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
720: else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
721: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: %s\n",symt);
722: }
723: return(0);
724: }
725: /* ---------------------------------------------------------------------------------*/
729: static PetscErrorCode PCHYPREGetType_HYPRE(PC pc,const char *name[])
730: {
731: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
734: *name = jac->hypre_type;
735: return(0);
736: }
740: static PetscErrorCode PCHYPRESetType_HYPRE(PC pc,const char name[])
741: {
742: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
744: PetscBool flag;
747: if (jac->hypre_type) {
748: PetscStrcmp(jac->hypre_type,name,&flag);
749: if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set");
750: return(0);
751: } else {
752: PetscStrallocpy(name, &jac->hypre_type);
753: }
755: jac->maxiter = PETSC_DEFAULT;
756: jac->tol = PETSC_DEFAULT;
757: jac->printstatistics = PetscLogPrintInfo;
759: PetscStrcmp("pilut",jac->hypre_type,&flag);
760: if (flag) {
761: PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver));
762: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
763: pc->ops->view = PCView_HYPRE_Pilut;
764: jac->destroy = HYPRE_ParCSRPilutDestroy;
765: jac->setup = HYPRE_ParCSRPilutSetup;
766: jac->solve = HYPRE_ParCSRPilutSolve;
767: jac->factorrowsize = PETSC_DEFAULT;
768: return(0);
769: }
770: PetscStrcmp("parasails",jac->hypre_type,&flag);
771: if (flag) {
772: PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver));
773: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
774: pc->ops->view = PCView_HYPRE_ParaSails;
775: jac->destroy = HYPRE_ParaSailsDestroy;
776: jac->setup = HYPRE_ParaSailsSetup;
777: jac->solve = HYPRE_ParaSailsSolve;
778: /* initialize */
779: jac->nlevels = 1;
780: jac->threshhold = .1;
781: jac->filter = .1;
782: jac->loadbal = 0;
783: if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE;
784: else jac->logging = (int) PETSC_FALSE;
786: jac->ruse = (int) PETSC_FALSE;
787: jac->symt = 0;
788: PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
789: PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
790: PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
791: PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
792: PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
793: PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
794: return(0);
795: }
796: PetscStrcmp("euclid",jac->hypre_type,&flag);
797: if (flag) {
798: HYPRE_EuclidCreate(jac->comm_hypre,&jac->hsolver);
799: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
800: pc->ops->view = PCView_HYPRE_Euclid;
801: jac->destroy = HYPRE_EuclidDestroy;
802: jac->setup = HYPRE_EuclidSetup;
803: jac->solve = HYPRE_EuclidSolve;
804: /* initialization */
805: jac->bjilu = PETSC_FALSE;
806: jac->levels = 1;
807: return(0);
808: }
809: PetscStrcmp("boomeramg",jac->hypre_type,&flag);
810: if (flag) {
811: HYPRE_BoomerAMGCreate(&jac->hsolver);
812: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_BoomerAMG;
813: pc->ops->view = PCView_HYPRE_BoomerAMG;
814: pc->ops->applytranspose = PCApplyTranspose_HYPRE_BoomerAMG;
815: pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
816: jac->destroy = HYPRE_BoomerAMGDestroy;
817: jac->setup = HYPRE_BoomerAMGSetup;
818: jac->solve = HYPRE_BoomerAMGSolve;
819: jac->applyrichardson = PETSC_FALSE;
820: /* these defaults match the hypre defaults */
821: jac->cycletype = 1;
822: jac->maxlevels = 25;
823: jac->maxiter = 1;
824: jac->tol = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
825: jac->truncfactor = 0.0;
826: jac->strongthreshold = .25;
827: jac->maxrowsum = .9;
828: jac->coarsentype = 6;
829: jac->measuretype = 0;
830: jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
831: jac->relaxtype[0] = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */
832: jac->relaxtype[2] = 9; /*G.E. */
833: jac->relaxweight = 1.0;
834: jac->outerrelaxweight = 1.0;
835: jac->relaxorder = 1;
836: jac->interptype = 0;
837: jac->agg_nl = 0;
838: jac->pmax = 0;
839: jac->truncfactor = 0.0;
840: jac->agg_num_paths = 1;
842: jac->nodal_coarsen = 0;
843: jac->nodal_relax = PETSC_FALSE;
844: jac->nodal_relax_levels = 1;
845: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
846: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
847: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
848: PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
849: PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
850: PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
851: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
852: PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
853: PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
854: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
855: PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
856: PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
857: PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
858: PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
859: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0])); /*defaults coarse to 9*/
860: PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */
861: return(0);
862: }
863: PetscFree(jac->hypre_type);
865: jac->hypre_type = NULL;
866: SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are pilut, parasails, euclid, boomeramg",name);
867: return(0);
868: }
870: /*
871: It only gets here if the HYPRE type has not been set before the call to
872: ...SetFromOptions() which actually is most of the time
873: */
876: static PetscErrorCode PCSetFromOptions_HYPRE(PC pc)
877: {
879: PetscInt indx;
880: const char *type[] = {"pilut","parasails","boomeramg","euclid"};
881: PetscBool flg;
884: PetscOptionsHead("HYPRE preconditioner options");
885: PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"boomeramg",&indx,&flg);
886: if (flg) {
887: PCHYPRESetType_HYPRE(pc,type[indx]);
888: } else {
889: PCHYPRESetType_HYPRE(pc,"boomeramg");
890: }
891: if (pc->ops->setfromoptions) {
892: pc->ops->setfromoptions(pc);
893: }
894: PetscOptionsTail();
895: return(0);
896: }
900: /*@C
901: PCHYPRESetType - Sets which hypre preconditioner you wish to use
903: Input Parameters:
904: + pc - the preconditioner context
905: - name - either pilut, parasails, boomeramg, euclid
907: Options Database Keys:
908: -pc_hypre_type - One of pilut, parasails, boomeramg, euclid
910: Level: intermediate
912: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
913: PCHYPRE
915: @*/
916: PetscErrorCode PCHYPRESetType(PC pc,const char name[])
917: {
923: PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));
924: return(0);
925: }
929: /*@C
930: PCHYPREGetType - Gets which hypre preconditioner you are using
932: Input Parameter:
933: . pc - the preconditioner context
935: Output Parameter:
936: . name - either pilut, parasails, boomeramg, euclid
938: Level: intermediate
940: .seealso: PCCreate(), PCHYPRESetType(), PCType (for list of available types), PC,
941: PCHYPRE
943: @*/
944: PetscErrorCode PCHYPREGetType(PC pc,const char *name[])
945: {
951: PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));
952: return(0);
953: }
955: /*MC
956: PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre
958: Options Database Keys:
959: + -pc_hypre_type - One of pilut, parasails, boomeramg, euclid
960: - Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX
961: preconditioner
963: Level: intermediate
965: Notes: Apart from pc_hypre_type (for which there is PCHYPRESetType()),
966: the many hypre options can ONLY be set via the options database (e.g. the command line
967: or with PetscOptionsSetValue(), there are no functions to set them)
969: The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_rtol refer to the number of iterations
970: (V-cycles) and tolerance that boomeramg does EACH time it is called. So for example, if
971: -pc_hypre_boomeramg_max_iter is set to 2 then 2-V-cycles are being used to define the preconditioner
972: (-pc_hypre_boomeramg_rtol should be set to 0.0 - the default - to strictly use a fixed number of
973: iterations per hypre call). -ksp_max_it and -ksp_rtol STILL determine the total number of iterations
974: and tolerance for the Krylov solver. For example, if -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10
975: then AT MOST twenty V-cycles of boomeramg will be called.
977: Note that the option -pc_hypre_boomeramg_relax_type_all defaults to symmetric relaxation
978: (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry.
979: Otherwise, you may want to use -pc_hypre_boomeramg_relax_type_all SOR/Jacobi.
980: If you wish to use BoomerAMG WITHOUT a Krylov method use -ksp_type richardson NOT -ksp_type preonly
981: and use -ksp_max_it to control the number of V-cycles.
982: (see the PETSc FAQ.html at the PETSc website under the Documentation tab).
984: 2007-02-03 Using HYPRE-1.11.1b, the routine HYPRE_BoomerAMGSolveT and the option
985: -pc_hypre_parasails_reuse were failing with SIGSEGV. Dalcin L.
987: See PCPFMG for access to the hypre Struct PFMG solver
989: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
990: PCHYPRESetType(), PCPFMG
992: M*/
996: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
997: {
998: PC_HYPRE *jac;
1002: PetscNewLog(pc,PC_HYPRE,&jac);
1004: pc->data = jac;
1005: pc->ops->destroy = PCDestroy_HYPRE;
1006: pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
1007: pc->ops->setup = PCSetUp_HYPRE;
1008: pc->ops->apply = PCApply_HYPRE;
1009: jac->comm_hypre = MPI_COMM_NULL;
1010: jac->hypre_type = NULL;
1011: /* duplicate communicator for hypre */
1012: MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1013: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);
1014: PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);
1015: return(0);
1016: }
1018: /* ---------------------------------------------------------------------------------------------------------------------------------*/
1020: /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */
1021: #include <petsc-private/matimpl.h>
1023: typedef struct {
1024: MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */
1025: HYPRE_StructSolver hsolver;
1027: /* keep copy of PFMG options used so may view them */
1028: PetscInt its;
1029: double tol;
1030: PetscInt relax_type;
1031: PetscInt rap_type;
1032: PetscInt num_pre_relax,num_post_relax;
1033: PetscInt max_levels;
1034: } PC_PFMG;
1038: PetscErrorCode PCDestroy_PFMG(PC pc)
1039: {
1041: PC_PFMG *ex = (PC_PFMG*) pc->data;
1044: if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1045: MPI_Comm_free(&ex->hcomm);
1046: PetscFree(pc->data);
1047: return(0);
1048: }
1050: static const char *PFMGRelaxType[] = {"Jacobi","Weighted-Jacobi","symmetric-Red/Black-Gauss-Seidel","Red/Black-Gauss-Seidel"};
1051: static const char *PFMGRAPType[] = {"Galerkin","non-Galerkin"};
1055: PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
1056: {
1058: PetscBool iascii;
1059: PC_PFMG *ex = (PC_PFMG*) pc->data;
1062: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1063: if (iascii) {
1064: PetscViewerASCIIPrintf(viewer," HYPRE PFMG preconditioning\n");
1065: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: max iterations %d\n",ex->its);
1066: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: tolerance %g\n",ex->tol);
1067: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
1068: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: RAP type %s\n",PFMGRAPType[ex->rap_type]);
1069: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
1070: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: max levels %d\n",ex->max_levels);
1071: }
1072: return(0);
1073: }
1078: PetscErrorCode PCSetFromOptions_PFMG(PC pc)
1079: {
1081: PC_PFMG *ex = (PC_PFMG*) pc->data;
1082: PetscBool flg = PETSC_FALSE;
1085: PetscOptionsHead("PFMG options");
1086: PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);
1087: if (flg) {
1088: int level=3;
1089: PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level));
1090: }
1091: PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);
1092: PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its));
1093: PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
1094: PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax));
1095: PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
1096: PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax));
1098: PetscOptionsInt("-pc_pfmg_max_levels","Max Levels for MG hierarchy","HYPRE_StructPFMGSetMaxLevels",ex->max_levels,&ex->max_levels,NULL);
1099: PetscStackCallStandard(HYPRE_StructPFMGSetMaxLevels,(ex->hsolver,ex->max_levels));
1101: PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);
1102: PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol));
1103: 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);
1104: PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type));
1105: PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);
1106: PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type));
1107: PetscOptionsTail();
1108: return(0);
1109: }
1113: PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
1114: {
1115: PetscErrorCode ierr;
1116: PC_PFMG *ex = (PC_PFMG*) pc->data;
1117: PetscScalar *xx,*yy;
1118: PetscInt ilower[3],iupper[3];
1119: Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
1122: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
1123: iupper[0] += ilower[0] - 1;
1124: iupper[1] += ilower[1] - 1;
1125: iupper[2] += ilower[2] - 1;
1127: /* copy x values over to hypre */
1128: PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
1129: VecGetArray(x,&xx);
1130: PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,ilower,iupper,xx));
1131: VecRestoreArray(x,&xx);
1132: PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
1133: PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
1135: /* copy solution values back to PETSc */
1136: VecGetArray(y,&yy);
1137: PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,ilower,iupper,yy));
1138: VecRestoreArray(y,&yy);
1139: return(0);
1140: }
1144: 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)
1145: {
1146: PC_PFMG *jac = (PC_PFMG*)pc->data;
1148: PetscInt oits;
1151: PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its));
1152: PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol));
1154: PCApply_PFMG(pc,b,y);
1155: PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,&oits));
1156: *outits = oits;
1157: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1158: else *reason = PCRICHARDSON_CONVERGED_RTOL;
1159: PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol));
1160: PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its));
1161: return(0);
1162: }
1167: PetscErrorCode PCSetUp_PFMG(PC pc)
1168: {
1169: PetscErrorCode ierr;
1170: PC_PFMG *ex = (PC_PFMG*) pc->data;
1171: Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
1172: PetscBool flg;
1175: PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESTRUCT,&flg);
1176: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESTRUCT with this preconditioner");
1178: /* create the hypre solver object and set its information */
1179: if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1180: PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1181: PCSetFromOptions_PFMG(pc);
1182: PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
1183: PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver));
1184: return(0);
1185: }
1188: /*MC
1189: PCPFMG - the hypre PFMG multigrid solver
1191: Level: advanced
1193: Options Database:
1194: + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner
1195: . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
1196: . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid
1197: . -pc_pfmg_tol <tol> tolerance of PFMG
1198: . -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
1199: - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin
1201: Notes: This is for CELL-centered descretizations
1203: This must be used with the MATHYPRESTRUCT matrix type.
1204: This is less general than in hypre, it supports only one block per process defined by a PETSc DMDA.
1206: .seealso: PCMG, MATHYPRESTRUCT
1207: M*/
1211: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
1212: {
1214: PC_PFMG *ex;
1217: PetscNew(PC_PFMG,&ex); \
1218: pc->data = ex;
1220: ex->its = 1;
1221: ex->tol = 1.e-8;
1222: ex->relax_type = 1;
1223: ex->rap_type = 0;
1224: ex->num_pre_relax = 1;
1225: ex->num_post_relax = 1;
1226: ex->max_levels = 0;
1228: pc->ops->setfromoptions = PCSetFromOptions_PFMG;
1229: pc->ops->view = PCView_PFMG;
1230: pc->ops->destroy = PCDestroy_PFMG;
1231: pc->ops->apply = PCApply_PFMG;
1232: pc->ops->applyrichardson = PCApplyRichardson_PFMG;
1233: pc->ops->setup = PCSetUp_PFMG;
1235: MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
1236: PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1237: return(0);
1238: }
1240: /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/
1242: /* we know we are working with a HYPRE_SStructMatrix */
1243: typedef struct {
1244: MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
1245: HYPRE_SStructSolver ss_solver;
1247: /* keep copy of SYSPFMG options used so may view them */
1248: PetscInt its;
1249: double tol;
1250: PetscInt relax_type;
1251: PetscInt num_pre_relax,num_post_relax;
1252: } PC_SysPFMG;
1256: PetscErrorCode PCDestroy_SysPFMG(PC pc)
1257: {
1259: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
1262: if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
1263: MPI_Comm_free(&ex->hcomm);
1264: PetscFree(pc->data);
1265: return(0);
1266: }
1268: static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi","Red/Black-Gauss-Seidel"};
1272: PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
1273: {
1275: PetscBool iascii;
1276: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
1279: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1280: if (iascii) {
1281: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG preconditioning\n");
1282: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: max iterations %d\n",ex->its);
1283: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: tolerance %g\n",ex->tol);
1284: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
1285: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
1286: }
1287: return(0);
1288: }
1293: PetscErrorCode PCSetFromOptions_SysPFMG(PC pc)
1294: {
1296: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
1297: PetscBool flg = PETSC_FALSE;
1300: PetscOptionsHead("SysPFMG options");
1301: PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);
1302: if (flg) {
1303: int level=3;
1304: PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level));
1305: }
1306: PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);
1307: PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its));
1308: PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
1309: PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax));
1310: PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
1311: PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax));
1313: PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);
1314: PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol));
1315: 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);
1316: PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type));
1317: PetscOptionsTail();
1318: return(0);
1319: }
1323: PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
1324: {
1325: PetscErrorCode ierr;
1326: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
1327: PetscScalar *xx,*yy;
1328: PetscInt ilower[3],iupper[3];
1329: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
1330: PetscInt ordering= mx->dofs_order;
1331: PetscInt nvars = mx->nvars;
1332: PetscInt part = 0;
1333: PetscInt size;
1334: PetscInt i;
1337: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
1338: iupper[0] += ilower[0] - 1;
1339: iupper[1] += ilower[1] - 1;
1340: iupper[2] += ilower[2] - 1;
1342: size = 1;
1343: for (i= 0; i< 3; i++) size *= (iupper[i]-ilower[i]+1);
1345: /* copy x values over to hypre for variable ordering */
1346: if (ordering) {
1347: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
1348: VecGetArray(x,&xx);
1349: for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,ilower,iupper,i,xx+(size*i)));
1350: VecRestoreArray(x,&xx);
1351: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
1352: PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
1353: PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
1355: /* copy solution values back to PETSc */
1356: VecGetArray(y,&yy);
1357: for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,ilower,iupper,i,yy+(size*i)));
1358: VecRestoreArray(y,&yy);
1359: } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
1360: PetscScalar *z;
1361: PetscInt j, k;
1363: PetscMalloc(nvars*size*sizeof(PetscScalar),&z);
1364: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
1365: VecGetArray(x,&xx);
1367: /* transform nodal to hypre's variable ordering for sys_pfmg */
1368: for (i= 0; i< size; i++) {
1369: k= i*nvars;
1370: for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
1371: }
1372: for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,ilower,iupper,i,z+(size*i)));
1373: VecRestoreArray(x,&xx);
1374: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
1375: PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
1377: /* copy solution values back to PETSc */
1378: VecGetArray(y,&yy);
1379: for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,ilower,iupper,i,z+(size*i)));
1380: /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
1381: for (i= 0; i< size; i++) {
1382: k= i*nvars;
1383: for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
1384: }
1385: VecRestoreArray(y,&yy);
1386: PetscFree(z);
1387: }
1388: return(0);
1389: }
1393: 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)
1394: {
1395: PC_SysPFMG *jac = (PC_SysPFMG*)pc->data;
1397: PetscInt oits;
1400: PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its));
1401: PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol));
1402: PCApply_SysPFMG(pc,b,y);
1403: PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,&oits));
1404: *outits = oits;
1405: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1406: else *reason = PCRICHARDSON_CONVERGED_RTOL;
1407: PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol));
1408: PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its));
1409: return(0);
1410: }
1415: PetscErrorCode PCSetUp_SysPFMG(PC pc)
1416: {
1417: PetscErrorCode ierr;
1418: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
1419: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
1420: PetscBool flg;
1423: PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESSTRUCT,&flg);
1424: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESSTRUCT with this preconditioner");
1426: /* create the hypre sstruct solver object and set its information */
1427: if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
1428: PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
1429: PCSetFromOptions_SysPFMG(pc);
1430: PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver));
1431: PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
1432: return(0);
1433: }
1436: /*MC
1437: PCSysPFMG - the hypre SysPFMG multigrid solver
1439: Level: advanced
1441: Options Database:
1442: + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner
1443: . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
1444: . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid
1445: . -pc_syspfmg_tol <tol> tolerance of SysPFMG
1446: . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel
1448: Notes: This is for CELL-centered descretizations
1450: This must be used with the MATHYPRESSTRUCT matrix type.
1451: This is less general than in hypre, it supports only one part, and one block per process defined by a PETSc DMDA.
1452: Also, only cell-centered variables.
1454: .seealso: PCMG, MATHYPRESSTRUCT
1455: M*/
1459: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
1460: {
1462: PC_SysPFMG *ex;
1465: PetscNew(PC_SysPFMG,&ex); \
1466: pc->data = ex;
1468: ex->its = 1;
1469: ex->tol = 1.e-8;
1470: ex->relax_type = 1;
1471: ex->num_pre_relax = 1;
1472: ex->num_post_relax = 1;
1474: pc->ops->setfromoptions = PCSetFromOptions_SysPFMG;
1475: pc->ops->view = PCView_SysPFMG;
1476: pc->ops->destroy = PCDestroy_SysPFMG;
1477: pc->ops->apply = PCApply_SysPFMG;
1478: pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
1479: pc->ops->setup = PCSetUp_SysPFMG;
1481: MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
1482: PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
1483: return(0);
1484: }