Actual source code: hypre.c
petsc-3.5.4 2015-05-23
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: static PetscBool cite = PETSC_FALSE;
12: static const char hypreCitation[] = "@manual{hypre-web-page,\n title = {{\\sl hypre}: High Performance Preconditioners},\n organization = {Lawrence Livermore National Laboratory},\n note = {\\url{http://www.llnl.gov/CASC/hypre/}}\n}\n";
14: /*
15: Private context (data structure) for the preconditioner.
16: */
17: typedef struct {
18: HYPRE_Solver hsolver;
19: HYPRE_IJMatrix ij;
20: HYPRE_IJVector b,x;
22: HYPRE_Int (*destroy)(HYPRE_Solver);
23: HYPRE_Int (*solve)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
24: HYPRE_Int (*setup)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
26: MPI_Comm comm_hypre;
27: char *hypre_type;
29: /* options for Pilut and BoomerAMG*/
30: PetscInt maxiter;
31: double tol;
33: /* options for Pilut */
34: PetscInt factorrowsize;
36: /* options for ParaSails */
37: PetscInt nlevels;
38: double threshhold;
39: double filter;
40: PetscInt sym;
41: double loadbal;
42: PetscInt logging;
43: PetscInt ruse;
44: PetscInt symt;
46: /* options for Euclid */
47: PetscBool bjilu;
48: PetscInt levels;
50: /* options for Euclid and BoomerAMG */
51: PetscBool printstatistics;
53: /* options for BoomerAMG */
54: PetscInt cycletype;
55: PetscInt maxlevels;
56: double strongthreshold;
57: double maxrowsum;
58: PetscInt gridsweeps[3];
59: PetscInt coarsentype;
60: PetscInt measuretype;
61: PetscInt relaxtype[3];
62: double relaxweight;
63: double outerrelaxweight;
64: PetscInt relaxorder;
65: double truncfactor;
66: PetscBool applyrichardson;
67: PetscInt pmax;
68: PetscInt interptype;
69: PetscInt agg_nl;
70: PetscInt agg_num_paths;
71: PetscInt nodal_coarsen;
72: PetscBool nodal_relax;
73: PetscInt nodal_relax_levels;
74: } PC_HYPRE;
78: PetscErrorCode PCHYPREGetSolver(PC pc,HYPRE_Solver *hsolver)
79: {
80: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
83: *hsolver = jac->hsolver;
84: return(0);
85: }
89: static PetscErrorCode PCSetUp_HYPRE(PC pc)
90: {
91: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
92: PetscErrorCode ierr;
93: HYPRE_ParCSRMatrix hmat;
94: HYPRE_ParVector bv,xv;
95: PetscInt bs;
98: if (!jac->hypre_type) {
99: PCHYPRESetType(pc,"boomeramg");
100: }
102: if (pc->setupcalled) {
103: /* always destroy the old matrix and create a new memory;
104: hope this does not churn the memory too much. The problem
105: is I do not know if it is possible to put the matrix back to
106: its initial state so that we can directly copy the values
107: the second time through. */
108: PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij));
109: jac->ij = 0;
110: }
112: if (!jac->ij) { /* create the matrix the first time through */
113: MatHYPRE_IJMatrixCreate(pc->pmat,&jac->ij);
114: }
115: if (!jac->b) { /* create the vectors the first time through */
116: Vec x,b;
117: MatGetVecs(pc->pmat,&x,&b);
118: VecHYPRE_IJVectorCreate(x,&jac->x);
119: VecHYPRE_IJVectorCreate(b,&jac->b);
120: VecDestroy(&x);
121: VecDestroy(&b);
122: }
124: /* special case for BoomerAMG */
125: if (jac->setup == HYPRE_BoomerAMGSetup) {
126: MatGetBlockSize(pc->pmat,&bs);
127: if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
128: };
129: MatHYPRE_IJMatrixCopy(pc->pmat,jac->ij);
130: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
131: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&bv));
132: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&xv));
133: PetscStackCall("HYPRE_SetupXXX",(*jac->setup)(jac->hsolver,hmat,bv,xv););
134: return(0);
135: }
137: /*
138: Replaces the address where the HYPRE vector points to its data with the address of
139: PETSc's data. Saves the old address so it can be reset when we are finished with it.
140: Allows use to get the data into a HYPRE vector without the cost of memcopies
141: */
142: #define HYPREReplacePointer(b,newvalue,savedvalue) { \
143: hypre_ParVector *par_vector = (hypre_ParVector*)hypre_IJVectorObject(((hypre_IJVector*)b)); \
144: hypre_Vector *local_vector = hypre_ParVectorLocalVector(par_vector); \
145: savedvalue = local_vector->data; \
146: local_vector->data = newvalue; \
147: }
151: static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
152: {
153: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
154: PetscErrorCode ierr;
155: HYPRE_ParCSRMatrix hmat;
156: PetscScalar *bv,*xv;
157: HYPRE_ParVector jbv,jxv;
158: PetscScalar *sbv,*sxv;
159: PetscInt hierr;
162: PetscCitationsRegister(hypreCitation,&cite);
163: if (!jac->applyrichardson) {VecSet(x,0.0);}
164: VecGetArray(b,&bv);
165: VecGetArray(x,&xv);
166: HYPREReplacePointer(jac->b,bv,sbv);
167: HYPREReplacePointer(jac->x,xv,sxv);
169: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
170: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
171: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));
172: PetscStackCall("Hypre solve",h(*jac->solve)(jac->hsolver,hmat,jbv,jxv);
173: if (hierr && hierr != HYPRE_ERROR_CONV) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
174: if (hierr) hypre__global_error = 0;);
176: HYPREReplacePointer(jac->b,sbv,bv);
177: HYPREReplacePointer(jac->x,sxv,xv);
178: VecRestoreArray(x,&xv);
179: VecRestoreArray(b,&bv);
180: return(0);
181: }
185: static PetscErrorCode PCDestroy_HYPRE(PC pc)
186: {
187: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
191: if (jac->ij) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij));
192: if (jac->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->b));
193: if (jac->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->x));
194: if (jac->destroy) PetscStackCall("HYPRE_DestroyXXX",(*jac->destroy)(jac->hsolver););
195: PetscFree(jac->hypre_type);
196: if (jac->comm_hypre != MPI_COMM_NULL) { MPI_Comm_free(&(jac->comm_hypre));}
197: PetscFree(pc->data);
199: PetscObjectChangeTypeName((PetscObject)pc,0);
200: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);
201: PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);
202: return(0);
203: }
205: /* --------------------------------------------------------------------------------------------*/
208: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PC pc)
209: {
210: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
212: PetscBool flag;
215: PetscOptionsHead("HYPRE Pilut Options");
216: PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
217: if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter));
218: PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
219: if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol));
220: PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
221: if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize));
222: PetscOptionsTail();
223: return(0);
224: }
228: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
229: {
230: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
232: PetscBool iascii;
235: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
236: if (iascii) {
237: PetscViewerASCIIPrintf(viewer," HYPRE Pilut preconditioning\n");
238: if (jac->maxiter != PETSC_DEFAULT) {
239: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: maximum number of iterations %d\n",jac->maxiter);
240: } else {
241: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default maximum number of iterations \n");
242: }
243: if (jac->tol != PETSC_DEFAULT) {
244: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: drop tolerance %g\n",(double)jac->tol);
245: } else {
246: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default drop tolerance \n");
247: }
248: if (jac->factorrowsize != PETSC_DEFAULT) {
249: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: factor row size %d\n",jac->factorrowsize);
250: } else {
251: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default factor row size \n");
252: }
253: }
254: return(0);
255: }
257: /* --------------------------------------------------------------------------------------------*/
260: static PetscErrorCode PCSetFromOptions_HYPRE_Euclid(PC pc)
261: {
262: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
264: PetscBool flag;
265: char *args[8],levels[16];
266: PetscInt cnt = 0;
269: PetscOptionsHead("HYPRE Euclid Options");
270: PetscOptionsInt("-pc_hypre_euclid_levels","Number of levels of fill ILU(k)","None",jac->levels,&jac->levels,&flag);
271: if (flag) {
272: if (jac->levels < 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be nonegative",jac->levels);
273: PetscSNPrintf(levels,sizeof(levels),"%D",jac->levels);
274: args[cnt++] = (char*)"-level"; args[cnt++] = levels;
275: }
276: PetscOptionsBool("-pc_hypre_euclid_bj","Use block Jacobi ILU(k)","None",jac->bjilu,&jac->bjilu,NULL);
277: if (jac->bjilu) {
278: args[cnt++] =(char*) "-bj"; args[cnt++] = (char*)"1";
279: }
281: PetscOptionsBool("-pc_hypre_euclid_print_statistics","Print statistics","None",jac->printstatistics,&jac->printstatistics,NULL);
282: if (jac->printstatistics) {
283: args[cnt++] = (char*)"-eu_stats"; args[cnt++] = (char*)"1";
284: args[cnt++] = (char*)"-eu_mem"; args[cnt++] = (char*)"1";
285: }
286: PetscOptionsTail();
287: if (cnt) PetscStackCallStandard(HYPRE_EuclidSetParams,(jac->hsolver,cnt,args));
288: return(0);
289: }
293: static PetscErrorCode PCView_HYPRE_Euclid(PC pc,PetscViewer viewer)
294: {
295: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
297: PetscBool iascii;
300: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
301: if (iascii) {
302: PetscViewerASCIIPrintf(viewer," HYPRE Euclid preconditioning\n");
303: PetscViewerASCIIPrintf(viewer," HYPRE Euclid: number of levels %d\n",jac->levels);
304: if (jac->bjilu) {
305: PetscViewerASCIIPrintf(viewer," HYPRE Euclid: Using block Jacobi ILU instead of parallel ILU\n");
306: }
307: }
308: return(0);
309: }
311: /* --------------------------------------------------------------------------------------------*/
315: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
316: {
317: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
318: PetscErrorCode ierr;
319: HYPRE_ParCSRMatrix hmat;
320: PetscScalar *bv,*xv;
321: HYPRE_ParVector jbv,jxv;
322: PetscScalar *sbv,*sxv;
323: PetscInt hierr;
326: PetscCitationsRegister(hypreCitation,&cite);
327: VecSet(x,0.0);
328: VecGetArray(b,&bv);
329: VecGetArray(x,&xv);
330: HYPREReplacePointer(jac->b,bv,sbv);
331: HYPREReplacePointer(jac->x,xv,sxv);
333: PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
334: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
335: PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));
337: hHYPRE_BoomerAMGSolveT(jac->hsolver,hmat,jbv,jxv);
338: /* error code of 1 in BoomerAMG merely means convergence not achieved */
339: if (hierr && (hierr != 1)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
340: if (hierr) hypre__global_error = 0;
342: HYPREReplacePointer(jac->b,sbv,bv);
343: HYPREReplacePointer(jac->x,sxv,xv);
344: VecRestoreArray(x,&xv);
345: VecRestoreArray(b,&bv);
346: return(0);
347: }
349: /* static array length */
350: #define ALEN(a) (sizeof(a)/sizeof((a)[0]))
352: static const char *HYPREBoomerAMGCycleType[] = {"","V","W"};
353: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
354: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
355: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
356: static const char *HYPREBoomerAMGRelaxType[] = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi",
357: "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi",
358: "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination",
359: "" /* 10 */, "" /* 11 */, "" /* 12 */, "" /* 13 */, "" /* 14 */,
360: "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"};
361: static const char *HYPREBoomerAMGInterpType[] = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
362: "ext+i-cc", "standard", "standard-wts", "", "", "FF", "FF1"};
365: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PC pc)
366: {
367: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
369: PetscInt n,indx,level;
370: PetscBool flg, tmp_truth;
371: double tmpdbl, twodbl[2];
374: PetscOptionsHead("HYPRE BoomerAMG Options");
375: PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);
376: if (flg) {
377: jac->cycletype = indx+1;
378: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
379: }
380: PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
381: if (flg) {
382: if (jac->maxlevels < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
383: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
384: }
385: PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
386: if (flg) {
387: if (jac->maxiter < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
388: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
389: }
390: PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);
391: if (flg) {
392: if (jac->tol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance %g must be greater than or equal to zero",(double)jac->tol);
393: PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
394: }
396: PetscOptionsScalar("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);
397: if (flg) {
398: if (jac->truncfactor < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Truncation factor %g must be great than or equal zero",(double)jac->truncfactor);
399: PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
400: }
402: PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);
403: if (flg) {
404: if (jac->pmax < 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"P_max %g must be greater than or equal to zero",(double)jac->pmax);
405: PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
406: }
408: PetscOptionsInt("-pc_hypre_boomeramg_agg_nl","Number of levels of aggressive coarsening","None",jac->agg_nl,&jac->agg_nl,&flg);
409: if (flg) {
410: if (jac->agg_nl < 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %g must be greater than or equal to zero",(double)jac->agg_nl);
412: PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
413: }
416: PetscOptionsInt("-pc_hypre_boomeramg_agg_num_paths","Number of paths for aggressive coarsening","None",jac->agg_num_paths,&jac->agg_num_paths,&flg);
417: if (flg) {
418: if (jac->agg_num_paths < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of paths %g must be greater than or equal to 1",(double)jac->agg_num_paths);
420: PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
421: }
424: PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
425: if (flg) {
426: if (jac->strongthreshold < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Strong threshold %g must be great than or equal zero",(double)jac->strongthreshold);
427: PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
428: }
429: PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
430: if (flg) {
431: if (jac->maxrowsum < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %g must be greater than zero",(double)jac->maxrowsum);
432: if (jac->maxrowsum > 1.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %g must be less than or equal one",(double)jac->maxrowsum);
433: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
434: }
436: /* Grid sweeps */
437: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);
438: if (flg) {
439: PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx));
440: /* modify the jac structure so we can view the updated options with PC_View */
441: jac->gridsweeps[0] = indx;
442: jac->gridsweeps[1] = indx;
443: /*defaults coarse to 1 */
444: jac->gridsweeps[2] = 1;
445: }
447: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);
448: if (flg) {
449: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1));
450: jac->gridsweeps[0] = indx;
451: }
452: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);
453: if (flg) {
454: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2));
455: jac->gridsweeps[1] = indx;
456: }
457: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);
458: if (flg) {
459: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3));
460: jac->gridsweeps[2] = indx;
461: }
463: /* Relax type */
464: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
465: if (flg) {
466: jac->relaxtype[0] = jac->relaxtype[1] = indx;
467: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx));
468: /* by default, coarse type set to 9 */
469: jac->relaxtype[2] = 9;
471: }
472: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
473: if (flg) {
474: jac->relaxtype[0] = indx;
475: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1));
476: }
477: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
478: if (flg) {
479: jac->relaxtype[1] = indx;
480: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2));
481: }
482: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);
483: if (flg) {
484: jac->relaxtype[2] = indx;
485: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3));
486: }
488: /* Relaxation Weight */
489: 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);
490: if (flg) {
491: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl));
492: jac->relaxweight = tmpdbl;
493: }
495: n = 2;
496: twodbl[0] = twodbl[1] = 1.0;
497: PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
498: if (flg) {
499: if (n == 2) {
500: indx = (int)PetscAbsReal(twodbl[1]);
501: PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx));
502: } 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);
503: }
505: /* Outer relaxation Weight */
506: 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);
507: if (flg) {
508: PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl));
509: jac->outerrelaxweight = tmpdbl;
510: }
512: n = 2;
513: twodbl[0] = twodbl[1] = 1.0;
514: PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
515: if (flg) {
516: if (n == 2) {
517: indx = (int)PetscAbsReal(twodbl[1]);
518: PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx));
519: } 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);
520: }
522: /* the Relax Order */
523: PetscOptionsBool("-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", PETSC_FALSE, &tmp_truth, &flg);
525: if (flg) {
526: jac->relaxorder = 0;
527: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
528: }
529: PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);
530: if (flg) {
531: jac->measuretype = indx;
532: PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
533: }
534: /* update list length 3/07 */
535: PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);
536: if (flg) {
537: jac->coarsentype = indx;
538: PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
539: }
541: /* new 3/07 */
542: PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);
543: if (flg) {
544: jac->interptype = indx;
545: PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
546: }
548: PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
549: if (flg) {
550: level = 3;
551: PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);
553: jac->printstatistics = PETSC_TRUE;
554: PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level));
555: }
557: PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);
558: if (flg) {
559: level = 3;
560: PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);
562: jac->printstatistics = PETSC_TRUE;
563: PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level));
564: }
566: PetscOptionsBool("-pc_hypre_boomeramg_nodal_coarsen", "HYPRE_BoomerAMGSetNodal()", "None", PETSC_FALSE, &tmp_truth, &flg);
567: if (flg && tmp_truth) {
568: jac->nodal_coarsen = 1;
569: PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,1));
570: }
572: PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);
573: if (flg && tmp_truth) {
574: PetscInt tmp_int;
575: PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);
576: if (flg) jac->nodal_relax_levels = tmp_int;
577: PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6));
578: PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1));
579: PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0));
580: PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels));
581: }
583: PetscOptionsTail();
584: return(0);
585: }
589: 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)
590: {
591: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
593: PetscInt oits;
596: PetscCitationsRegister(hypreCitation,&cite);
597: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter));
598: PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol));
599: jac->applyrichardson = PETSC_TRUE;
600: PCApply_HYPRE(pc,b,y);
601: jac->applyrichardson = PETSC_FALSE;
602: PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
603: *outits = oits;
604: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
605: else *reason = PCRICHARDSON_CONVERGED_RTOL;
606: PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
607: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
608: return(0);
609: }
614: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
615: {
616: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
618: PetscBool iascii;
621: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
622: if (iascii) {
623: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG preconditioning\n");
624: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);
625: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum number of levels %d\n",jac->maxlevels);
626: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum number of iterations PER hypre call %d\n",jac->maxiter);
627: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Convergence tolerance PER hypre call %g\n",(double)jac->tol);
628: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Threshold for strong coupling %g\n",(double)jac->strongthreshold);
629: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Interpolation truncation factor %g\n",(double)jac->truncfactor);
630: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Interpolation: max elements per row %d\n",jac->pmax);
631: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Number of levels of aggressive coarsening %d\n",jac->agg_nl);
632: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Number of paths for aggressive coarsening %d\n",jac->agg_num_paths);
634: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum row sums %g\n",(double)jac->maxrowsum);
636: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps down %d\n",jac->gridsweeps[0]);
637: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps up %d\n",jac->gridsweeps[1]);
638: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps on coarse %d\n",jac->gridsweeps[2]);
640: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax down %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[0]]);
641: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax up %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[1]]);
642: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax on coarse %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[2]]);
644: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax weight (all) %g\n",(double)jac->relaxweight);
645: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Outer relax weight (all) %g\n",(double)jac->outerrelaxweight);
647: if (jac->relaxorder) {
648: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using CF-relaxation\n");
649: } else {
650: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Not using CF-relaxation\n");
651: }
652: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Measure type %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
653: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Coarsen type %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
654: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Interpolation type %s\n",HYPREBoomerAMGInterpType[jac->interptype]);
655: if (jac->nodal_coarsen) {
656: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal coarsening (with HYPRE_BOOMERAMGSetNodal())\n");
657: }
658: if (jac->nodal_relax) {
659: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal relaxation via Schwarz smoothing on levels %d\n",jac->nodal_relax_levels);
660: }
661: }
662: return(0);
663: }
665: /* --------------------------------------------------------------------------------------------*/
668: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PC pc)
669: {
670: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
672: PetscInt indx;
673: PetscBool flag;
674: const char *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};
677: PetscOptionsHead("HYPRE ParaSails Options");
678: PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);
679: PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshhold,&jac->threshhold,&flag);
680: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
682: PetscOptionsReal("-pc_hypre_parasails_filter","filter","None",jac->filter,&jac->filter,&flag);
683: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
685: PetscOptionsReal("-pc_hypre_parasails_loadbal","Load balance","None",jac->loadbal,&jac->loadbal,&flag);
686: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
688: PetscOptionsBool("-pc_hypre_parasails_logging","Print info to screen","None",(PetscBool)jac->logging,(PetscBool*)&jac->logging,&flag);
689: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
691: PetscOptionsBool("-pc_hypre_parasails_reuse","Reuse nonzero pattern in preconditioner","None",(PetscBool)jac->ruse,(PetscBool*)&jac->ruse,&flag);
692: if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
694: PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);
695: if (flag) {
696: jac->symt = indx;
697: PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
698: }
700: PetscOptionsTail();
701: return(0);
702: }
706: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
707: {
708: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
710: PetscBool iascii;
711: const char *symt = 0;;
714: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
715: if (iascii) {
716: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails preconditioning\n");
717: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: nlevels %d\n",jac->nlevels);
718: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: threshold %g\n",(double)jac->threshhold);
719: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: filter %g\n",(double)jac->filter);
720: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: load balance %g\n",(double)jac->loadbal);
721: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: reuse nonzero structure %s\n",PetscBools[jac->ruse]);
722: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: print info to screen %s\n",PetscBools[jac->logging]);
723: if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
724: else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
725: else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
726: else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
727: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: %s\n",symt);
728: }
729: return(0);
730: }
731: /* ---------------------------------------------------------------------------------*/
735: static PetscErrorCode PCHYPREGetType_HYPRE(PC pc,const char *name[])
736: {
737: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
740: *name = jac->hypre_type;
741: return(0);
742: }
746: static PetscErrorCode PCHYPRESetType_HYPRE(PC pc,const char name[])
747: {
748: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
750: PetscBool flag;
753: if (jac->hypre_type) {
754: PetscStrcmp(jac->hypre_type,name,&flag);
755: if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set");
756: return(0);
757: } else {
758: PetscStrallocpy(name, &jac->hypre_type);
759: }
761: jac->maxiter = PETSC_DEFAULT;
762: jac->tol = PETSC_DEFAULT;
763: jac->printstatistics = PetscLogPrintInfo;
765: PetscStrcmp("pilut",jac->hypre_type,&flag);
766: if (flag) {
767: PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver));
768: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
769: pc->ops->view = PCView_HYPRE_Pilut;
770: jac->destroy = HYPRE_ParCSRPilutDestroy;
771: jac->setup = HYPRE_ParCSRPilutSetup;
772: jac->solve = HYPRE_ParCSRPilutSolve;
773: jac->factorrowsize = PETSC_DEFAULT;
774: return(0);
775: }
776: PetscStrcmp("parasails",jac->hypre_type,&flag);
777: if (flag) {
778: PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver));
779: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
780: pc->ops->view = PCView_HYPRE_ParaSails;
781: jac->destroy = HYPRE_ParaSailsDestroy;
782: jac->setup = HYPRE_ParaSailsSetup;
783: jac->solve = HYPRE_ParaSailsSolve;
784: /* initialize */
785: jac->nlevels = 1;
786: jac->threshhold = .1;
787: jac->filter = .1;
788: jac->loadbal = 0;
789: if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE;
790: else jac->logging = (int) PETSC_FALSE;
792: jac->ruse = (int) PETSC_FALSE;
793: jac->symt = 0;
794: PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
795: PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
796: PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
797: PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
798: PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
799: PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
800: return(0);
801: }
802: PetscStrcmp("euclid",jac->hypre_type,&flag);
803: if (flag) {
804: HYPRE_EuclidCreate(jac->comm_hypre,&jac->hsolver);
805: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
806: pc->ops->view = PCView_HYPRE_Euclid;
807: jac->destroy = HYPRE_EuclidDestroy;
808: jac->setup = HYPRE_EuclidSetup;
809: jac->solve = HYPRE_EuclidSolve;
810: /* initialization */
811: jac->bjilu = PETSC_FALSE;
812: jac->levels = 1;
813: return(0);
814: }
815: PetscStrcmp("boomeramg",jac->hypre_type,&flag);
816: if (flag) {
817: HYPRE_BoomerAMGCreate(&jac->hsolver);
818: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_BoomerAMG;
819: pc->ops->view = PCView_HYPRE_BoomerAMG;
820: pc->ops->applytranspose = PCApplyTranspose_HYPRE_BoomerAMG;
821: pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
822: jac->destroy = HYPRE_BoomerAMGDestroy;
823: jac->setup = HYPRE_BoomerAMGSetup;
824: jac->solve = HYPRE_BoomerAMGSolve;
825: jac->applyrichardson = PETSC_FALSE;
826: /* these defaults match the hypre defaults */
827: jac->cycletype = 1;
828: jac->maxlevels = 25;
829: jac->maxiter = 1;
830: jac->tol = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
831: jac->truncfactor = 0.0;
832: jac->strongthreshold = .25;
833: jac->maxrowsum = .9;
834: jac->coarsentype = 6;
835: jac->measuretype = 0;
836: jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
837: jac->relaxtype[0] = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */
838: jac->relaxtype[2] = 9; /*G.E. */
839: jac->relaxweight = 1.0;
840: jac->outerrelaxweight = 1.0;
841: jac->relaxorder = 1;
842: jac->interptype = 0;
843: jac->agg_nl = 0;
844: jac->pmax = 0;
845: jac->truncfactor = 0.0;
846: jac->agg_num_paths = 1;
848: jac->nodal_coarsen = 0;
849: jac->nodal_relax = PETSC_FALSE;
850: jac->nodal_relax_levels = 1;
851: PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
852: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
853: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
854: PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
855: PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
856: PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
857: PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
858: PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
859: PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
860: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
861: PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
862: PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
863: PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
864: PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
865: PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0])); /*defaults coarse to 9*/
866: PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */
867: return(0);
868: }
869: PetscFree(jac->hypre_type);
871: jac->hypre_type = NULL;
872: SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are pilut, parasails, euclid, boomeramg",name);
873: return(0);
874: }
876: /*
877: It only gets here if the HYPRE type has not been set before the call to
878: ...SetFromOptions() which actually is most of the time
879: */
882: static PetscErrorCode PCSetFromOptions_HYPRE(PC pc)
883: {
885: PetscInt indx;
886: const char *type[] = {"pilut","parasails","boomeramg","euclid"};
887: PetscBool flg;
890: PetscOptionsHead("HYPRE preconditioner options");
891: PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"boomeramg",&indx,&flg);
892: if (flg) {
893: PCHYPRESetType_HYPRE(pc,type[indx]);
894: } else {
895: PCHYPRESetType_HYPRE(pc,"boomeramg");
896: }
897: if (pc->ops->setfromoptions) {
898: pc->ops->setfromoptions(pc);
899: }
900: PetscOptionsTail();
901: return(0);
902: }
906: /*@C
907: PCHYPRESetType - Sets which hypre preconditioner you wish to use
909: Input Parameters:
910: + pc - the preconditioner context
911: - name - either pilut, parasails, boomeramg, euclid
913: Options Database Keys:
914: -pc_hypre_type - One of pilut, parasails, boomeramg, euclid
916: Level: intermediate
918: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
919: PCHYPRE
921: @*/
922: PetscErrorCode PCHYPRESetType(PC pc,const char name[])
923: {
929: PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));
930: return(0);
931: }
935: /*@C
936: PCHYPREGetType - Gets which hypre preconditioner you are using
938: Input Parameter:
939: . pc - the preconditioner context
941: Output Parameter:
942: . name - either pilut, parasails, boomeramg, euclid
944: Level: intermediate
946: .seealso: PCCreate(), PCHYPRESetType(), PCType (for list of available types), PC,
947: PCHYPRE
949: @*/
950: PetscErrorCode PCHYPREGetType(PC pc,const char *name[])
951: {
957: PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));
958: return(0);
959: }
961: /*MC
962: PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre
964: Options Database Keys:
965: + -pc_hypre_type - One of pilut, parasails, boomeramg, euclid
966: - Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX
967: preconditioner
969: Level: intermediate
971: Notes: Apart from pc_hypre_type (for which there is PCHYPRESetType()),
972: the many hypre options can ONLY be set via the options database (e.g. the command line
973: or with PetscOptionsSetValue(), there are no functions to set them)
975: The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_rtol refer to the number of iterations
976: (V-cycles) and tolerance that boomeramg does EACH time it is called. So for example, if
977: -pc_hypre_boomeramg_max_iter is set to 2 then 2-V-cycles are being used to define the preconditioner
978: (-pc_hypre_boomeramg_rtol should be set to 0.0 - the default - to strictly use a fixed number of
979: iterations per hypre call). -ksp_max_it and -ksp_rtol STILL determine the total number of iterations
980: and tolerance for the Krylov solver. For example, if -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10
981: then AT MOST twenty V-cycles of boomeramg will be called.
983: Note that the option -pc_hypre_boomeramg_relax_type_all defaults to symmetric relaxation
984: (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry.
985: Otherwise, you may want to use -pc_hypre_boomeramg_relax_type_all SOR/Jacobi.
986: If you wish to use BoomerAMG WITHOUT a Krylov method use -ksp_type richardson NOT -ksp_type preonly
987: and use -ksp_max_it to control the number of V-cycles.
988: (see the PETSc FAQ.html at the PETSc website under the Documentation tab).
990: 2007-02-03 Using HYPRE-1.11.1b, the routine HYPRE_BoomerAMGSolveT and the option
991: -pc_hypre_parasails_reuse were failing with SIGSEGV. Dalcin L.
993: See PCPFMG for access to the hypre Struct PFMG solver
995: .seealso: PCCreate(), PCSetType(), PCType (for list of available types), PC,
996: PCHYPRESetType(), PCPFMG
998: M*/
1002: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
1003: {
1004: PC_HYPRE *jac;
1008: PetscNewLog(pc,&jac);
1010: pc->data = jac;
1011: pc->ops->destroy = PCDestroy_HYPRE;
1012: pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
1013: pc->ops->setup = PCSetUp_HYPRE;
1014: pc->ops->apply = PCApply_HYPRE;
1015: jac->comm_hypre = MPI_COMM_NULL;
1016: jac->hypre_type = NULL;
1017: /* duplicate communicator for hypre */
1018: MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1019: PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);
1020: PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);
1021: return(0);
1022: }
1024: /* ---------------------------------------------------------------------------------------------------------------------------------*/
1026: /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */
1027: #include <petsc-private/matimpl.h>
1029: typedef struct {
1030: MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */
1031: HYPRE_StructSolver hsolver;
1033: /* keep copy of PFMG options used so may view them */
1034: PetscInt its;
1035: double tol;
1036: PetscInt relax_type;
1037: PetscInt rap_type;
1038: PetscInt num_pre_relax,num_post_relax;
1039: PetscInt max_levels;
1040: } PC_PFMG;
1044: PetscErrorCode PCDestroy_PFMG(PC pc)
1045: {
1047: PC_PFMG *ex = (PC_PFMG*) pc->data;
1050: if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1051: MPI_Comm_free(&ex->hcomm);
1052: PetscFree(pc->data);
1053: return(0);
1054: }
1056: static const char *PFMGRelaxType[] = {"Jacobi","Weighted-Jacobi","symmetric-Red/Black-Gauss-Seidel","Red/Black-Gauss-Seidel"};
1057: static const char *PFMGRAPType[] = {"Galerkin","non-Galerkin"};
1061: PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
1062: {
1064: PetscBool iascii;
1065: PC_PFMG *ex = (PC_PFMG*) pc->data;
1068: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1069: if (iascii) {
1070: PetscViewerASCIIPrintf(viewer," HYPRE PFMG preconditioning\n");
1071: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: max iterations %d\n",ex->its);
1072: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: tolerance %g\n",ex->tol);
1073: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
1074: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: RAP type %s\n",PFMGRAPType[ex->rap_type]);
1075: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
1076: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: max levels %d\n",ex->max_levels);
1077: }
1078: return(0);
1079: }
1084: PetscErrorCode PCSetFromOptions_PFMG(PC pc)
1085: {
1087: PC_PFMG *ex = (PC_PFMG*) pc->data;
1088: PetscBool flg = PETSC_FALSE;
1091: PetscOptionsHead("PFMG options");
1092: PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);
1093: if (flg) {
1094: int level=3;
1095: PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level));
1096: }
1097: PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);
1098: PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its));
1099: PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
1100: PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax));
1101: PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
1102: PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax));
1104: PetscOptionsInt("-pc_pfmg_max_levels","Max Levels for MG hierarchy","HYPRE_StructPFMGSetMaxLevels",ex->max_levels,&ex->max_levels,NULL);
1105: PetscStackCallStandard(HYPRE_StructPFMGSetMaxLevels,(ex->hsolver,ex->max_levels));
1107: PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);
1108: PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol));
1109: 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);
1110: PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type));
1111: PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);
1112: PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type));
1113: PetscOptionsTail();
1114: return(0);
1115: }
1119: PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
1120: {
1121: PetscErrorCode ierr;
1122: PC_PFMG *ex = (PC_PFMG*) pc->data;
1123: PetscScalar *xx,*yy;
1124: PetscInt ilower[3],iupper[3];
1125: Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
1128: PetscCitationsRegister(hypreCitation,&cite);
1129: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
1130: iupper[0] += ilower[0] - 1;
1131: iupper[1] += ilower[1] - 1;
1132: iupper[2] += ilower[2] - 1;
1134: /* copy x values over to hypre */
1135: PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
1136: VecGetArray(x,&xx);
1137: PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,xx));
1138: VecRestoreArray(x,&xx);
1139: PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
1140: PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
1142: /* copy solution values back to PETSc */
1143: VecGetArray(y,&yy);
1144: PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy));
1145: VecRestoreArray(y,&yy);
1146: return(0);
1147: }
1151: 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)
1152: {
1153: PC_PFMG *jac = (PC_PFMG*)pc->data;
1155: PetscInt oits;
1158: PetscCitationsRegister(hypreCitation,&cite);
1159: PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its));
1160: PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol));
1162: PCApply_PFMG(pc,b,y);
1163: PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
1164: *outits = oits;
1165: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1166: else *reason = PCRICHARDSON_CONVERGED_RTOL;
1167: PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol));
1168: PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its));
1169: return(0);
1170: }
1175: PetscErrorCode PCSetUp_PFMG(PC pc)
1176: {
1177: PetscErrorCode ierr;
1178: PC_PFMG *ex = (PC_PFMG*) pc->data;
1179: Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
1180: PetscBool flg;
1183: PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESTRUCT,&flg);
1184: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESTRUCT with this preconditioner");
1186: /* create the hypre solver object and set its information */
1187: if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1188: PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1189: PCSetFromOptions_PFMG(pc);
1190: PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
1191: PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver));
1192: return(0);
1193: }
1196: /*MC
1197: PCPFMG - the hypre PFMG multigrid solver
1199: Level: advanced
1201: Options Database:
1202: + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner
1203: . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
1204: . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid
1205: . -pc_pfmg_tol <tol> tolerance of PFMG
1206: . -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
1207: - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin
1209: Notes: This is for CELL-centered descretizations
1211: This must be used with the MATHYPRESTRUCT matrix type.
1212: This is less general than in hypre, it supports only one block per process defined by a PETSc DMDA.
1214: .seealso: PCMG, MATHYPRESTRUCT
1215: M*/
1219: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
1220: {
1222: PC_PFMG *ex;
1225: PetscNew(&ex); \
1226: pc->data = ex;
1228: ex->its = 1;
1229: ex->tol = 1.e-8;
1230: ex->relax_type = 1;
1231: ex->rap_type = 0;
1232: ex->num_pre_relax = 1;
1233: ex->num_post_relax = 1;
1234: ex->max_levels = 0;
1236: pc->ops->setfromoptions = PCSetFromOptions_PFMG;
1237: pc->ops->view = PCView_PFMG;
1238: pc->ops->destroy = PCDestroy_PFMG;
1239: pc->ops->apply = PCApply_PFMG;
1240: pc->ops->applyrichardson = PCApplyRichardson_PFMG;
1241: pc->ops->setup = PCSetUp_PFMG;
1243: MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
1244: PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1245: return(0);
1246: }
1248: /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/
1250: /* we know we are working with a HYPRE_SStructMatrix */
1251: typedef struct {
1252: MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
1253: HYPRE_SStructSolver ss_solver;
1255: /* keep copy of SYSPFMG options used so may view them */
1256: PetscInt its;
1257: double tol;
1258: PetscInt relax_type;
1259: PetscInt num_pre_relax,num_post_relax;
1260: } PC_SysPFMG;
1264: PetscErrorCode PCDestroy_SysPFMG(PC pc)
1265: {
1267: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
1270: if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
1271: MPI_Comm_free(&ex->hcomm);
1272: PetscFree(pc->data);
1273: return(0);
1274: }
1276: static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi","Red/Black-Gauss-Seidel"};
1280: PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
1281: {
1283: PetscBool iascii;
1284: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
1287: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1288: if (iascii) {
1289: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG preconditioning\n");
1290: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: max iterations %d\n",ex->its);
1291: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: tolerance %g\n",ex->tol);
1292: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
1293: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
1294: }
1295: return(0);
1296: }
1301: PetscErrorCode PCSetFromOptions_SysPFMG(PC pc)
1302: {
1304: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
1305: PetscBool flg = PETSC_FALSE;
1308: PetscOptionsHead("SysPFMG options");
1309: PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);
1310: if (flg) {
1311: int level=3;
1312: PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level));
1313: }
1314: PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);
1315: PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its));
1316: PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
1317: PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax));
1318: PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
1319: PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax));
1321: PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);
1322: PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol));
1323: 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);
1324: PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type));
1325: PetscOptionsTail();
1326: return(0);
1327: }
1331: PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
1332: {
1333: PetscErrorCode ierr;
1334: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
1335: PetscScalar *xx,*yy;
1336: PetscInt ilower[3],iupper[3];
1337: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
1338: PetscInt ordering= mx->dofs_order;
1339: PetscInt nvars = mx->nvars;
1340: PetscInt part = 0;
1341: PetscInt size;
1342: PetscInt i;
1345: PetscCitationsRegister(hypreCitation,&cite);
1346: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
1347: iupper[0] += ilower[0] - 1;
1348: iupper[1] += ilower[1] - 1;
1349: iupper[2] += ilower[2] - 1;
1351: size = 1;
1352: for (i= 0; i< 3; i++) size *= (iupper[i]-ilower[i]+1);
1354: /* copy x values over to hypre for variable ordering */
1355: if (ordering) {
1356: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
1357: VecGetArray(x,&xx);
1358: for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,xx+(size*i)));
1359: VecRestoreArray(x,&xx);
1360: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
1361: PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
1362: PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
1364: /* copy solution values back to PETSc */
1365: VecGetArray(y,&yy);
1366: for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i)));
1367: VecRestoreArray(y,&yy);
1368: } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
1369: PetscScalar *z;
1370: PetscInt j, k;
1372: PetscMalloc1(nvars*size,&z);
1373: PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
1374: VecGetArray(x,&xx);
1376: /* transform nodal to hypre's variable ordering for sys_pfmg */
1377: for (i= 0; i< size; i++) {
1378: k= i*nvars;
1379: for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
1380: }
1381: for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
1382: VecRestoreArray(x,&xx);
1383: PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
1384: PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
1386: /* copy solution values back to PETSc */
1387: VecGetArray(y,&yy);
1388: for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
1389: /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
1390: for (i= 0; i< size; i++) {
1391: k= i*nvars;
1392: for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
1393: }
1394: VecRestoreArray(y,&yy);
1395: PetscFree(z);
1396: }
1397: return(0);
1398: }
1402: 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)
1403: {
1404: PC_SysPFMG *jac = (PC_SysPFMG*)pc->data;
1406: PetscInt oits;
1409: PetscCitationsRegister(hypreCitation,&cite);
1410: PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its));
1411: PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol));
1412: PCApply_SysPFMG(pc,b,y);
1413: PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,(HYPRE_Int *)&oits));
1414: *outits = oits;
1415: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1416: else *reason = PCRICHARDSON_CONVERGED_RTOL;
1417: PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol));
1418: PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its));
1419: return(0);
1420: }
1425: PetscErrorCode PCSetUp_SysPFMG(PC pc)
1426: {
1427: PetscErrorCode ierr;
1428: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
1429: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
1430: PetscBool flg;
1433: PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESSTRUCT,&flg);
1434: if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESSTRUCT with this preconditioner");
1436: /* create the hypre sstruct solver object and set its information */
1437: if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
1438: PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
1439: PCSetFromOptions_SysPFMG(pc);
1440: PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver));
1441: PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
1442: return(0);
1443: }
1446: /*MC
1447: PCSysPFMG - the hypre SysPFMG multigrid solver
1449: Level: advanced
1451: Options Database:
1452: + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner
1453: . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
1454: . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid
1455: . -pc_syspfmg_tol <tol> tolerance of SysPFMG
1456: . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel
1458: Notes: This is for CELL-centered descretizations
1460: This must be used with the MATHYPRESSTRUCT matrix type.
1461: This is less general than in hypre, it supports only one part, and one block per process defined by a PETSc DMDA.
1462: Also, only cell-centered variables.
1464: .seealso: PCMG, MATHYPRESSTRUCT
1465: M*/
1469: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
1470: {
1472: PC_SysPFMG *ex;
1475: PetscNew(&ex); \
1476: pc->data = ex;
1478: ex->its = 1;
1479: ex->tol = 1.e-8;
1480: ex->relax_type = 1;
1481: ex->num_pre_relax = 1;
1482: ex->num_post_relax = 1;
1484: pc->ops->setfromoptions = PCSetFromOptions_SysPFMG;
1485: pc->ops->view = PCView_SysPFMG;
1486: pc->ops->destroy = PCDestroy_SysPFMG;
1487: pc->ops->apply = PCApply_SysPFMG;
1488: pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
1489: pc->ops->setup = PCSetUp_SysPFMG;
1491: MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
1492: PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
1493: return(0);
1494: }