Actual source code: hypre.c
petsc-3.3-p7 2013-05-11
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);
22:
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;
76: static PetscErrorCode PCSetUp_HYPRE(PC pc)
77: {
78: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
79: PetscErrorCode ierr;
80: HYPRE_ParCSRMatrix hmat;
81: HYPRE_ParVector bv,xv;
82: PetscInt bs;
85: if (!jac->hypre_type) {
86: PCHYPRESetType(pc,"boomeramg");
87: }
89: if (pc->setupcalled) {
90: /* always destroy the old matrix and create a new memory;
91: hope this does not churn the memory too much. The problem
92: is I do not know if it is possible to put the matrix back to
93: its initial state so that we can directly copy the values
94: the second time through. */
95: PetscStackCallHypre(0,HYPRE_IJMatrixDestroy,(jac->ij));
96: jac->ij = 0;
97: }
99: if (!jac->ij) { /* create the matrix the first time through */
100: MatHYPRE_IJMatrixCreate(pc->pmat,&jac->ij);
101: }
102: if (!jac->b) { /* create the vectors the first time through */
103: Vec x,b;
104: MatGetVecs(pc->pmat,&x,&b);
105: VecHYPRE_IJVectorCreate(x,&jac->x);
106: VecHYPRE_IJVectorCreate(b,&jac->b);
107: VecDestroy(&x);
108: VecDestroy(&b);
109: }
111: /* special case for BoomerAMG */
112: if (jac->setup == HYPRE_BoomerAMGSetup) {
113: MatGetBlockSize(pc->pmat,&bs);
114: if (bs > 1) {
115: PetscStackCallHypre(0,HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
116: }
117: };
118: MatHYPRE_IJMatrixCopy(pc->pmat,jac->ij);
119: PetscStackCallHypre(0,HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
120: PetscStackCallHypre(0,HYPRE_IJVectorGetObject,(jac->b,(void**)&bv));
121: PetscStackCallHypre(0,HYPRE_IJVectorGetObject,(jac->x,(void**)&xv));
122: PetscStackCallHypre("HYPRE_SetupXXX",(*jac->setup),(jac->hsolver,hmat,bv,xv));
123: return(0);
124: }
126: /*
127: Replaces the address where the HYPRE vector points to its data with the address of
128: PETSc's data. Saves the old address so it can be reset when we are finished with it.
129: Allows use to get the data into a HYPRE vector without the cost of memcopies
130: */
131: #define HYPREReplacePointer(b,newvalue,savedvalue) {\
132: hypre_ParVector *par_vector = (hypre_ParVector *)hypre_IJVectorObject(((hypre_IJVector*)b));\
133: hypre_Vector *local_vector = hypre_ParVectorLocalVector(par_vector);\
134: savedvalue = local_vector->data;\
135: local_vector->data = newvalue;}
139: static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
140: {
141: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
142: PetscErrorCode ierr;
143: HYPRE_ParCSRMatrix hmat;
144: PetscScalar *bv,*xv;
145: HYPRE_ParVector jbv,jxv;
146: PetscScalar *sbv,*sxv;
147: PetscInt hierr;
150: if (!jac->applyrichardson) {VecSet(x,0.0);}
151: VecGetArray(b,&bv);
152: VecGetArray(x,&xv);
153: HYPREReplacePointer(jac->b,bv,sbv);
154: HYPREReplacePointer(jac->x,xv,sxv);
156: PetscStackCallHypre(0,HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
157: PetscStackCallHypre(0,HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
158: PetscStackCallHypre(0,HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));
159: h(*jac->solve)(jac->hsolver,hmat,jbv,jxv);
160: if (hierr && hierr != HYPRE_ERROR_CONV) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
161: if (hierr) hypre__global_error = 0;
162:
163: HYPREReplacePointer(jac->b,sbv,bv);
164: HYPREReplacePointer(jac->x,sxv,xv);
165: VecRestoreArray(x,&xv);
166: VecRestoreArray(b,&bv);
167: return(0);
168: }
172: static PetscErrorCode PCDestroy_HYPRE(PC pc)
173: {
174: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
178: if (jac->ij) PetscStackCallHypre(0,HYPRE_IJMatrixDestroy,(jac->ij));
179: if (jac->b) PetscStackCallHypre(0,HYPRE_IJVectorDestroy,(jac->b));
180: if (jac->x) PetscStackCallHypre(0,HYPRE_IJVectorDestroy,(jac->x));
181: if (jac->destroy) PetscStackCallHypre("HYPRE_DistroyXXX",(*jac->destroy),(jac->hsolver));
182: PetscFree(jac->hypre_type);
183: if (jac->comm_hypre != MPI_COMM_NULL) { MPI_Comm_free(&(jac->comm_hypre));}
184: PetscFree(pc->data);
186: PetscObjectChangeTypeName((PetscObject)pc,0);
187: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCHYPRESetType_C","",PETSC_NULL);
188: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCHYPREGetType_C","",PETSC_NULL);
189: return(0);
190: }
192: /* --------------------------------------------------------------------------------------------*/
195: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PC pc)
196: {
197: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
199: PetscBool flag;
202: PetscOptionsHead("HYPRE Pilut Options");
203: PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
204: if (flag) PetscStackCallHypre(0,HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter));
205: PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
206: if (flag) PetscStackCallHypre(0,HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol));
207: PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
208: if (flag) PetscStackCallHypre(0,HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize));
209: PetscOptionsTail();
210: return(0);
211: }
215: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
216: {
217: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
219: PetscBool iascii;
222: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
223: if (iascii) {
224: PetscViewerASCIIPrintf(viewer," HYPRE Pilut preconditioning\n");
225: if (jac->maxiter != PETSC_DEFAULT) {
226: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: maximum number of iterations %d\n",jac->maxiter);
227: } else {
228: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default maximum number of iterations \n");
229: }
230: if (jac->tol != PETSC_DEFAULT) {
231: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: drop tolerance %G\n",jac->tol);
232: } else {
233: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default drop tolerance \n");
234: }
235: if (jac->factorrowsize != PETSC_DEFAULT) {
236: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: factor row size %d\n",jac->factorrowsize);
237: } else {
238: PetscViewerASCIIPrintf(viewer," HYPRE Pilut: default factor row size \n");
239: }
240: }
241: return(0);
242: }
244: /* --------------------------------------------------------------------------------------------*/
247: static PetscErrorCode PCSetFromOptions_HYPRE_Euclid(PC pc)
248: {
249: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
251: PetscBool flag;
252: char *args[8],levels[16];
253: PetscInt cnt = 0;
256: PetscOptionsHead("HYPRE Euclid Options");
257: PetscOptionsInt("-pc_hypre_euclid_levels","Number of levels of fill ILU(k)","None",jac->levels,&jac->levels,&flag);
258: if (flag) {
259: if (jac->levels < 0) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be nonegative",jac->levels);
260: PetscSNPrintf(levels,sizeof levels,"%D",jac->levels);
261: args[cnt++] = (char*)"-level"; args[cnt++] = levels;
262: }
263: PetscOptionsBool("-pc_hypre_euclid_bj","Use block Jacobi ILU(k)","None",jac->bjilu,&jac->bjilu,PETSC_NULL);
264: if (jac->bjilu) {
265: args[cnt++] =(char*) "-bj"; args[cnt++] = (char*)"1";
266: }
267:
268: PetscOptionsBool("-pc_hypre_euclid_print_statistics","Print statistics","None",jac->printstatistics,&jac->printstatistics,PETSC_NULL);
269: if (jac->printstatistics) {
270: args[cnt++] = (char*)"-eu_stats"; args[cnt++] = (char*)"1";
271: args[cnt++] = (char*)"-eu_mem"; args[cnt++] = (char*)"1";
272: }
273: PetscOptionsTail();
274: if (cnt) PetscStackCallHypre(0,HYPRE_EuclidSetParams,(jac->hsolver,cnt,args));
275: return(0);
276: }
280: static PetscErrorCode PCView_HYPRE_Euclid(PC pc,PetscViewer viewer)
281: {
282: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
284: PetscBool iascii;
287: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
288: if (iascii) {
289: PetscViewerASCIIPrintf(viewer," HYPRE Euclid preconditioning\n");
290: PetscViewerASCIIPrintf(viewer," HYPRE Euclid: number of levels %d\n",jac->levels);
291: if (jac->bjilu) {
292: PetscViewerASCIIPrintf(viewer," HYPRE Euclid: Using block Jacobi ILU instead of parallel ILU\n");
293: }
294: }
295: return(0);
296: }
298: /* --------------------------------------------------------------------------------------------*/
302: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
303: {
304: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
305: PetscErrorCode ierr;
306: HYPRE_ParCSRMatrix hmat;
307: PetscScalar *bv,*xv;
308: HYPRE_ParVector jbv,jxv;
309: PetscScalar *sbv,*sxv;
310: PetscInt hierr;
313: VecSet(x,0.0);
314: VecGetArray(b,&bv);
315: VecGetArray(x,&xv);
316: HYPREReplacePointer(jac->b,bv,sbv);
317: HYPREReplacePointer(jac->x,xv,sxv);
319: PetscStackCallHypre(0,HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
320: PetscStackCallHypre(0,HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
321: PetscStackCallHypre(0,HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));
322:
323: hHYPRE_BoomerAMGSolveT(jac->hsolver,hmat,jbv,jxv);
324: /* error code of 1 in BoomerAMG merely means convergence not achieved */
325: if (hierr && (hierr != 1)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
326: if (hierr) hypre__global_error = 0;
327:
328: HYPREReplacePointer(jac->b,sbv,bv);
329: HYPREReplacePointer(jac->x,sxv,xv);
330: VecRestoreArray(x,&xv);
331: VecRestoreArray(b,&bv);
332: return(0);
333: }
335: static const char *HYPREBoomerAMGCycleType[] = {"","V","W"};
336: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
337: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
338: static const char *HYPREBoomerAMGRelaxType[] = {"Jacobi","sequential-Gauss-Seidel","","SOR/Jacobi","backward-SOR/Jacobi","","symmetric-SOR/Jacobi",
339: "","","Gaussian-elimination"};
340: static const char *HYPREBoomerAMGInterpType[] = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
341: "ext+i-cc", "standard", "standard-wts", "", "", "FF", "FF1"};
344: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PC pc)
345: {
346: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
348: PetscInt n,indx,level;
349: PetscBool flg, tmp_truth;
350: double tmpdbl, twodbl[2];
353: PetscOptionsHead("HYPRE BoomerAMG Options");
354: PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);
355: if (flg) {
356: jac->cycletype = indx+1;
357: PetscStackCallHypre(0,HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
358: }
359: PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
360: if (flg) {
361: if (jac->maxlevels < 2) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
362: PetscStackCallHypre(0,HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
363: }
364: PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
365: if (flg) {
366: if (jac->maxiter < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
367: PetscStackCallHypre(0,HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
368: }
369: PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);
370: if (flg) {
371: if (jac->tol < 0.0) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Tolerance %G must be greater than or equal to zero",jac->tol);
372: PetscStackCallHypre(0,HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
373: }
375: PetscOptionsScalar("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);
376: if (flg) {
377: if (jac->truncfactor < 0.0) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Truncation factor %G must be great than or equal zero",jac->truncfactor);
378: PetscStackCallHypre(0,HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
379: }
381: PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator ( 0=unlimited )","None",jac->pmax,&jac->pmax,&flg);
382: if (flg) {
383: if (jac->pmax < 0) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"P_max %G must be greater than or equal to zero",jac->pmax);
384: PetscStackCallHypre(0,HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
385: }
387: PetscOptionsInt("-pc_hypre_boomeramg_agg_nl","Number of levels of aggressive coarsening","None",jac->agg_nl,&jac->agg_nl,&flg);
388: if (flg) {
389: if (jac->agg_nl < 0) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %G must be greater than or equal to zero",jac->agg_nl);
391: PetscStackCallHypre(0,HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
392: }
395: PetscOptionsInt("-pc_hypre_boomeramg_agg_num_paths","Number of paths for aggressive coarsening","None",jac->agg_num_paths,&jac->agg_num_paths,&flg);
396: if (flg) {
397: if (jac->agg_num_paths < 1) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Number of paths %G must be greater than or equal to 1",jac->agg_num_paths);
399: PetscStackCallHypre(0,HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
400: }
403: PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
404: if (flg) {
405: if (jac->strongthreshold < 0.0) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Strong threshold %G must be great than or equal zero",jac->strongthreshold);
406: PetscStackCallHypre(0,HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
407: }
408: PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
409: if (flg) {
410: if (jac->maxrowsum < 0.0) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %G must be greater than zero",jac->maxrowsum);
411: if (jac->maxrowsum > 1.0) SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %G must be less than or equal one",jac->maxrowsum);
412: PetscStackCallHypre(0,HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
413: }
415: /* Grid sweeps */
416: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);
417: if (flg) {
418: PetscStackCallHypre(0,HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx));
419: /* modify the jac structure so we can view the updated options with PC_View */
420: jac->gridsweeps[0] = indx;
421: jac->gridsweeps[1] = indx;
422: /*defaults coarse to 1 */
423: jac->gridsweeps[2] = 1;
424: }
426: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx ,&flg);
427: if (flg) {
428: PetscStackCallHypre(0,HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1));
429: jac->gridsweeps[0] = indx;
430: }
431: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);
432: if (flg) {
433: PetscStackCallHypre(0,HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2));
434: jac->gridsweeps[1] = indx;
435: }
436: PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);
437: if (flg) {
438: PetscStackCallHypre(0,HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3));
439: jac->gridsweeps[2] = indx;
440: }
442: /* Relax type */
443: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[6],&indx,&flg);
444: if (flg) {
445: jac->relaxtype[0] = jac->relaxtype[1] = indx;
446: PetscStackCallHypre(0,HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx));
447: /* by default, coarse type set to 9 */
448: jac->relaxtype[2] = 9;
449:
450: }
451: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[6],&indx,&flg);
452: if (flg) {
453: jac->relaxtype[0] = indx;
454: PetscStackCallHypre(0,HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1));
455: }
456: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[6],&indx,&flg);
457: if (flg) {
458: jac->relaxtype[1] = indx;
459: PetscStackCallHypre(0,HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2));
460: }
461: PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[9],&indx,&flg);
462: if (flg) {
463: jac->relaxtype[2] = indx;
464: PetscStackCallHypre(0,HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3));
465: }
467: /* Relaxation Weight */
468: 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);
469: if (flg) {
470: PetscStackCallHypre(0,HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl));
471: jac->relaxweight = tmpdbl;
472: }
474: n=2;
475: twodbl[0] = twodbl[1] = 1.0;
476: PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
477: if (flg) {
478: if (n == 2) {
479: indx = (int)PetscAbsReal(twodbl[1]);
480: PetscStackCallHypre(0,HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx));
481: } else {
482: SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Relax weight level: you must provide 2 values separated by a comma (and no space), you provided %d",n);
483: }
484: }
486: /* Outer relaxation Weight */
487: 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);
488: if (flg) {
489: PetscStackCallHypre(0,HYPRE_BoomerAMGSetOuterWt,( jac->hsolver, tmpdbl));
490: jac->outerrelaxweight = tmpdbl;
491: }
493: n=2;
494: twodbl[0] = twodbl[1] = 1.0;
495: PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
496: if (flg) {
497: if (n == 2) {
498: indx = (int)PetscAbsReal(twodbl[1]);
499: PetscStackCallHypre(0,HYPRE_BoomerAMGSetLevelOuterWt,( jac->hsolver, twodbl[0], indx));
500: } else {
501: SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_OUTOFRANGE,"Relax weight outer level: You must provide 2 values separated by a comma (and no space), you provided %d",n);
502: }
503: }
505: /* the Relax Order */
506: PetscOptionsBool( "-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", PETSC_FALSE, &tmp_truth, &flg);
508: if (flg) {
509: jac->relaxorder = 0;
510: PetscStackCallHypre(0,HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
511: }
512: PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,2,HYPREBoomerAMGMeasureType[0],&indx,&flg);
513: if (flg) {
514: jac->measuretype = indx;
515: PetscStackCallHypre(0,HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
516: }
517: /* update list length 3/07 */
518: PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,11,HYPREBoomerAMGCoarsenType[6],&indx,&flg);
519: if (flg) {
520: jac->coarsentype = indx;
521: PetscStackCallHypre(0,HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
522: }
523:
524: /* new 3/07 */
525: PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,14,HYPREBoomerAMGInterpType[0],&indx,&flg);
526: if (flg) {
527: jac->interptype = indx;
528: PetscStackCallHypre(0,HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
529: }
531: PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
532: if (flg) {
533: level = 3;
534: PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,PETSC_NULL);
535: jac->printstatistics = PETSC_TRUE;
536: PetscStackCallHypre(0,HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level));
537: }
539: PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);
540: if (flg) {
541: level = 3;
542: PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,PETSC_NULL);
543: jac->printstatistics = PETSC_TRUE;
544: PetscStackCallHypre(0,HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level));
545: }
547: PetscOptionsBool( "-pc_hypre_boomeramg_nodal_coarsen", "HYPRE_BoomerAMGSetNodal()", "None", PETSC_FALSE, &tmp_truth, &flg);
548: if (flg && tmp_truth) {
549: jac->nodal_coarsen = 1;
550: PetscStackCallHypre(0,HYPRE_BoomerAMGSetNodal,(jac->hsolver,1));
551: }
553: PetscOptionsBool( "-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);
554: if (flg && tmp_truth) {
555: PetscInt tmp_int;
556: PetscOptionsInt( "-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);
557: if (flg) jac->nodal_relax_levels = tmp_int;
558: PetscStackCallHypre(0,HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6));
559: PetscStackCallHypre(0,HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1));
560: PetscStackCallHypre(0,HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0));
561: PetscStackCallHypre(0,HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels));
562: }
564: PetscOptionsTail();
565: return(0);
566: }
570: 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)
571: {
572: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
574: PetscInt oits;
577: PetscStackCallHypre(0,HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter));
578: PetscStackCallHypre(0,HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol));
579: jac->applyrichardson = PETSC_TRUE;
580: PCApply_HYPRE(pc,b,y);
581: jac->applyrichardson = PETSC_FALSE;
582: PetscStackCallHypre(0,HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,&oits));
583: *outits = oits;
584: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
585: else *reason = PCRICHARDSON_CONVERGED_RTOL;
586: PetscStackCallHypre(0,HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
587: PetscStackCallHypre(0,HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
588: return(0);
589: }
594: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
595: {
596: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
598: PetscBool iascii;
601: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
602: if (iascii) {
603: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG preconditioning\n");
604: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);
605: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum number of levels %d\n",jac->maxlevels);
606: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum number of iterations PER hypre call %d\n",jac->maxiter);
607: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Convergence tolerance PER hypre call %G\n",jac->tol);
608: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Threshold for strong coupling %G\n",jac->strongthreshold);
609: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Interpolation truncation factor %G\n",jac->truncfactor);
610: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Interpolation: max elements per row %d\n",jac->pmax);
611: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Number of levels of aggressive coarsening %d\n",jac->agg_nl);
612: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Number of paths for aggressive coarsening %d\n",jac->agg_num_paths);
613:
614: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Maximum row sums %G\n",jac->maxrowsum);
616: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps down %d\n",jac->gridsweeps[0]);
617: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps up %d\n",jac->gridsweeps[1]);
618: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Sweeps on coarse %d\n",jac->gridsweeps[2]);
620: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax down %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[0]]);
621: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax up %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[1]]);
622: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax on coarse %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[2]]);
624: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Relax weight (all) %G\n",jac->relaxweight);
625: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Outer relax weight (all) %G\n",jac->outerrelaxweight);
627: if (jac->relaxorder) {
628: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using CF-relaxation\n");
629: } else {
630: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Not using CF-relaxation\n");
631: }
632: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Measure type %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
633: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Coarsen type %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
634: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Interpolation type %s\n",HYPREBoomerAMGInterpType[jac->interptype]);
635: if (jac->nodal_coarsen) {
636: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal coarsening (with HYPRE_BOOMERAMGSetNodal())\n");
637: }
638: if (jac->nodal_relax) {
639: PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal relaxation via Schwarz smoothing on levels %d\n",jac->nodal_relax_levels);
640: }
641: }
642: return(0);
643: }
645: /* --------------------------------------------------------------------------------------------*/
648: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PC pc)
649: {
650: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
652: PetscInt indx;
653: PetscBool flag;
654: const char *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};
657: PetscOptionsHead("HYPRE ParaSails Options");
658: PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);
659: PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshhold,&jac->threshhold,&flag);
660: if (flag) {
661: PetscStackCallHypre(0,HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
662: }
664: PetscOptionsReal("-pc_hypre_parasails_filter","filter","None",jac->filter,&jac->filter,&flag);
665: if (flag) {
666: PetscStackCallHypre(0,HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
667: }
669: PetscOptionsReal("-pc_hypre_parasails_loadbal","Load balance","None",jac->loadbal,&jac->loadbal,&flag);
670: if (flag) {
671: PetscStackCallHypre(0,HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
672: }
674: PetscOptionsBool("-pc_hypre_parasails_logging","Print info to screen","None",(PetscBool)jac->logging,(PetscBool*)&jac->logging,&flag);
675: if (flag) {
676: PetscStackCallHypre(0,HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
677: }
679: PetscOptionsBool("-pc_hypre_parasails_reuse","Reuse nonzero pattern in preconditioner","None",(PetscBool)jac->ruse,(PetscBool*)&jac->ruse,&flag);
680: if (flag) {
681: PetscStackCallHypre(0,HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
682: }
684: PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,3,symtlist[0],&indx,&flag);
685: if (flag) {
686: jac->symt = indx;
687: PetscStackCallHypre(0,HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
688: }
690: PetscOptionsTail();
691: return(0);
692: }
696: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
697: {
698: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
700: PetscBool iascii;
701: const char *symt = 0;;
704: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
705: if (iascii) {
706: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails preconditioning\n");
707: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: nlevels %d\n",jac->nlevels);
708: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: threshold %G\n",jac->threshhold);
709: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: filter %G\n",jac->filter);
710: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: load balance %G\n",jac->loadbal);
711: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: reuse nonzero structure %s\n",PetscBools[jac->ruse]);
712: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: print info to screen %s\n",PetscBools[jac->logging]);
713: if (!jac->symt) {
714: symt = "nonsymmetric matrix and preconditioner";
715: } else if (jac->symt == 1) {
716: symt = "SPD matrix and preconditioner";
717: } else if (jac->symt == 2) {
718: symt = "nonsymmetric matrix but SPD preconditioner";
719: } else {
720: SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
721: }
722: PetscViewerASCIIPrintf(viewer," HYPRE ParaSails: %s\n",symt);
723: }
724: return(0);
725: }
726: /* ---------------------------------------------------------------------------------*/
728: EXTERN_C_BEGIN
731: PetscErrorCode PCHYPREGetType_HYPRE(PC pc,const char *name[])
732: {
733: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
736: *name = jac->hypre_type;
737: return(0);
738: }
739: EXTERN_C_END
741: EXTERN_C_BEGIN
744: PetscErrorCode PCHYPRESetType_HYPRE(PC pc,const char name[])
745: {
746: PC_HYPRE *jac = (PC_HYPRE*)pc->data;
748: PetscBool flag;
751: if (jac->hypre_type) {
752: PetscStrcmp(jac->hypre_type,name,&flag);
753: if (!flag) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set");
754: return(0);
755: } else {
756: PetscStrallocpy(name, &jac->hypre_type);
757: }
758:
759: jac->maxiter = PETSC_DEFAULT;
760: jac->tol = PETSC_DEFAULT;
761: jac->printstatistics = PetscLogPrintInfo;
763: PetscStrcmp("pilut",jac->hypre_type,&flag);
764: if (flag) {
765: PetscStackCallHypre(0,HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver));
766: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
767: pc->ops->view = PCView_HYPRE_Pilut;
768: jac->destroy = HYPRE_ParCSRPilutDestroy;
769: jac->setup = HYPRE_ParCSRPilutSetup;
770: jac->solve = HYPRE_ParCSRPilutSolve;
771: jac->factorrowsize = PETSC_DEFAULT;
772: return(0);
773: }
774: PetscStrcmp("parasails",jac->hypre_type,&flag);
775: if (flag) {
776: PetscStackCallHypre(0,HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver));
777: pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
778: pc->ops->view = PCView_HYPRE_ParaSails;
779: jac->destroy = HYPRE_ParaSailsDestroy;
780: jac->setup = HYPRE_ParaSailsSetup;
781: jac->solve = HYPRE_ParaSailsSolve;
782: /* initialize */
783: jac->nlevels = 1;
784: jac->threshhold = .1;
785: jac->filter = .1;
786: jac->loadbal = 0;
787: if (PetscLogPrintInfo) {
788: jac->logging = (int) PETSC_TRUE;
789: } else {
790: jac->logging = (int) PETSC_FALSE;
791: }
792: jac->ruse = (int) PETSC_FALSE;
793: jac->symt = 0;
794: PetscStackCallHypre(0,HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
795: PetscStackCallHypre(0,HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
796: PetscStackCallHypre(0,HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
797: PetscStackCallHypre(0,HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
798: PetscStackCallHypre(0,HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
799: PetscStackCallHypre(0,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: PetscStackCallHypre(0,HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
852: PetscStackCallHypre(0,HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
853: PetscStackCallHypre(0,HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
854: PetscStackCallHypre(0,HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
855: PetscStackCallHypre(0,HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
856: PetscStackCallHypre(0,HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
857: PetscStackCallHypre(0,HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
858: PetscStackCallHypre(0,HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
859: PetscStackCallHypre(0,HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
860: PetscStackCallHypre(0,HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
861: PetscStackCallHypre(0,HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
862: PetscStackCallHypre(0,HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
863: PetscStackCallHypre(0,HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
864: PetscStackCallHypre(0,HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
865: PetscStackCallHypre(0,HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0])); /*defaults coarse to 9*/
866: PetscStackCallHypre(0,HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */
867: return(0);
868: }
869: PetscFree(jac->hypre_type);
870: jac->hypre_type = PETSC_NULL;
871: SETERRQ1(((PetscObject)pc)->comm,PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are pilut, parasails, euclid, boomeramg",name);
872: return(0);
873: }
874: EXTERN_C_END
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
915:
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
968:
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*/
1000: EXTERN_C_BEGIN
1003: PetscErrorCode PCCreate_HYPRE(PC pc)
1004: {
1005: PC_HYPRE *jac;
1009: PetscNewLog(pc,PC_HYPRE,&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 = PETSC_NULL;
1017: /* duplicate communicator for hypre */
1018: MPI_Comm_dup(((PetscObject)pc)->comm,&(jac->comm_hypre));
1019: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCHYPRESetType_C","PCHYPRESetType_HYPRE",PCHYPRESetType_HYPRE);
1020: PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCHYPREGetType_C","PCHYPREGetType_HYPRE",PCHYPREGetType_HYPRE);
1021: return(0);
1022: }
1023: EXTERN_C_END
1025: /* ---------------------------------------------------------------------------------------------------------------------------------*/
1027: /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */
1028: #include <petsc-private/matimpl.h>
1030: typedef struct {
1031: MPI_Comm hcomm; /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */
1032: HYPRE_StructSolver hsolver;
1034: /* keep copy of PFMG options used so may view them */
1035: PetscInt its;
1036: double tol;
1037: PetscInt relax_type;
1038: PetscInt rap_type;
1039: PetscInt num_pre_relax,num_post_relax;
1040: PetscInt max_levels;
1041: } PC_PFMG;
1045: PetscErrorCode PCDestroy_PFMG(PC pc)
1046: {
1048: PC_PFMG *ex = (PC_PFMG*) pc->data;
1051: if (ex->hsolver) {PetscStackCallHypre(0,HYPRE_StructPFMGDestroy,(ex->hsolver));}
1052: MPI_Comm_free(&ex->hcomm);
1053: PetscFree(pc->data);
1054: return(0);
1055: }
1057: static const char *PFMGRelaxType[] = {"Jacobi","Weighted-Jacobi","symmetric-Red/Black-Gauss-Seidel","Red/Black-Gauss-Seidel"};
1058: static const char *PFMGRAPType[] = {"Galerkin","non-Galerkin"};
1062: PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
1063: {
1065: PetscBool iascii;
1066: PC_PFMG *ex = (PC_PFMG*) pc->data;
1069: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1070: if (iascii) {
1071: PetscViewerASCIIPrintf(viewer," HYPRE PFMG preconditioning\n");
1072: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: max iterations %d\n",ex->its);
1073: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: tolerance %g\n",ex->tol);
1074: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
1075: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: RAP type %s\n",PFMGRAPType[ex->rap_type]);
1076: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
1077: PetscViewerASCIIPrintf(viewer," HYPRE PFMG: max levels %d\n",ex->max_levels);
1078: }
1079: return(0);
1080: }
1085: PetscErrorCode PCSetFromOptions_PFMG(PC pc)
1086: {
1088: PC_PFMG *ex = (PC_PFMG*) pc->data;
1089: PetscBool flg = PETSC_FALSE;
1092: PetscOptionsHead("PFMG options");
1093: PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,PETSC_NULL);
1094: if (flg) {
1095: int level=3;
1096: PetscStackCallHypre(0,HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level));
1097: }
1098: PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,PETSC_NULL);
1099: PetscStackCallHypre(0,HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its));
1100: PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,PETSC_NULL);
1101: PetscStackCallHypre(0,HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax));
1102: PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,PETSC_NULL);
1103: PetscStackCallHypre(0,HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax));
1105: PetscOptionsInt("-pc_pfmg_max_levels","Max Levels for MG hierarchy","HYPRE_StructPFMGSetMaxLevels",ex->max_levels,&ex->max_levels,PETSC_NULL);
1106: PetscStackCallHypre(0,HYPRE_StructPFMGSetMaxLevels,(ex->hsolver,ex->max_levels));
1108: PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,PETSC_NULL);
1109: PetscStackCallHypre(0,HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol));
1110: PetscOptionsEList("-pc_pfmg_relax_type","Relax type for the up and down cycles","HYPRE_StructPFMGSetRelaxType",PFMGRelaxType,4,PFMGRelaxType[ex->relax_type],&ex->relax_type,PETSC_NULL);
1111: PetscStackCallHypre(0,HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type));
1112: PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,2,PFMGRAPType[ex->rap_type],&ex->rap_type,PETSC_NULL);
1113: PetscStackCallHypre(0,HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type));
1114: PetscOptionsTail();
1115: return(0);
1116: }
1120: PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
1121: {
1122: PetscErrorCode ierr;
1123: PC_PFMG *ex = (PC_PFMG*) pc->data;
1124: PetscScalar *xx,*yy;
1125: PetscInt ilower[3],iupper[3];
1126: Mat_HYPREStruct *mx = (Mat_HYPREStruct *)(pc->pmat->data);
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: PetscStackCallHypre(0,HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
1136: VecGetArray(x,&xx);
1137: PetscStackCallHypre(0,HYPRE_StructVectorSetBoxValues,(mx->hb,ilower,iupper,xx));
1138: VecRestoreArray(x,&xx);
1139: PetscStackCallHypre(0,HYPRE_StructVectorAssemble,(mx->hb));
1140: PetscStackCallHypre(0,HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
1142: /* copy solution values back to PETSc */
1143: VecGetArray(y,&yy);
1144: PetscStackCallHypre(0,HYPRE_StructVectorGetBoxValues,(mx->hx,ilower,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: PetscStackCallHypre(0,HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its));
1159: PetscStackCallHypre(0,HYPRE_StructPFMGSetTol,(jac->hsolver,rtol));
1161: PCApply_PFMG(pc,b,y);
1162: PetscStackCallHypre(0,HYPRE_StructPFMGGetNumIterations,(jac->hsolver,&oits));
1163: *outits = oits;
1164: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1165: else *reason = PCRICHARDSON_CONVERGED_RTOL;
1166: PetscStackCallHypre(0,HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol));
1167: PetscStackCallHypre(0,HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its));
1168: return(0);
1169: }
1174: PetscErrorCode PCSetUp_PFMG(PC pc)
1175: {
1176: PetscErrorCode ierr;
1177: PC_PFMG *ex = (PC_PFMG*) pc->data;
1178: Mat_HYPREStruct *mx = (Mat_HYPREStruct *)(pc->pmat->data);
1179: PetscBool flg;
1182: PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESTRUCT,&flg);
1183: if (!flg) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESTRUCT with this preconditioner");
1185: /* create the hypre solver object and set its information */
1186: if (ex->hsolver) {
1187: PetscStackCallHypre(0,HYPRE_StructPFMGDestroy,(ex->hsolver));
1188: }
1189: PetscStackCallHypre(0,HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1190: PCSetFromOptions_PFMG(pc);
1191: PetscStackCallHypre(0,HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
1192: PetscStackCallHypre(0,HYPRE_StructPFMGSetZeroGuess,(ex->hsolver));
1193: return(0);
1194: }
1197: /*MC
1198: PCPFMG - the hypre PFMG multigrid solver
1200: Level: advanced
1202: Options Database:
1203: + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner
1204: . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
1205: . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid
1206: . -pc_pfmg_tol <tol> tolerance of PFMG
1207: . -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
1208: - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin
1210: Notes: This is for CELL-centered descretizations
1212: This must be used with the MATHYPRESTRUCT matrix type.
1213: This is less general than in hypre, it supports only one block per process defined by a PETSc DMDA.
1215: .seealso: PCMG, MATHYPRESTRUCT
1216: M*/
1218: EXTERN_C_BEGIN
1221: PetscErrorCode PCCreate_PFMG(PC pc)
1222: {
1224: PC_PFMG *ex;
1227: PetscNew(PC_PFMG,&ex);\
1228: pc->data = ex;
1230: ex->its = 1;
1231: ex->tol = 1.e-8;
1232: ex->relax_type = 1;
1233: ex->rap_type = 0;
1234: ex->num_pre_relax = 1;
1235: ex->num_post_relax = 1;
1236: ex->max_levels = 0;
1238: pc->ops->setfromoptions = PCSetFromOptions_PFMG;
1239: pc->ops->view = PCView_PFMG;
1240: pc->ops->destroy = PCDestroy_PFMG;
1241: pc->ops->apply = PCApply_PFMG;
1242: pc->ops->applyrichardson = PCApplyRichardson_PFMG;
1243: pc->ops->setup = PCSetUp_PFMG;
1244: MPI_Comm_dup(((PetscObject)pc)->comm,&(ex->hcomm));
1245: PetscStackCallHypre(0,HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1246: return(0);
1247: }
1248: EXTERN_C_END
1250: /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/
1252: /* we know we are working with a HYPRE_SStructMatrix */
1253: typedef struct {
1254: MPI_Comm hcomm; /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
1255: HYPRE_SStructSolver ss_solver;
1257: /* keep copy of SYSPFMG options used so may view them */
1258: PetscInt its;
1259: double tol;
1260: PetscInt relax_type;
1261: PetscInt num_pre_relax,num_post_relax;
1262: } PC_SysPFMG;
1266: PetscErrorCode PCDestroy_SysPFMG(PC pc)
1267: {
1269: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
1272: if (ex->ss_solver) {PetscStackCallHypre(0,HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));}
1273: MPI_Comm_free(&ex->hcomm);
1274: PetscFree(pc->data);
1275: return(0);
1276: }
1278: static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi","Red/Black-Gauss-Seidel"};
1282: PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
1283: {
1285: PetscBool iascii;
1286: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
1289: PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1290: if (iascii) {
1291: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG preconditioning\n");
1292: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: max iterations %d\n",ex->its);
1293: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: tolerance %g\n",ex->tol);
1294: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
1295: PetscViewerASCIIPrintf(viewer," HYPRE SysPFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
1296: }
1297: return(0);
1298: }
1303: PetscErrorCode PCSetFromOptions_SysPFMG(PC pc)
1304: {
1306: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
1307: PetscBool flg = PETSC_FALSE;
1310: PetscOptionsHead("SysPFMG options");
1311: PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,PETSC_NULL);
1312: if (flg) {
1313: int level=3;
1314: PetscStackCallHypre(0,HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level));
1315: }
1316: PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,PETSC_NULL);
1317: PetscStackCallHypre(0,HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its));
1318: PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,PETSC_NULL);
1319: PetscStackCallHypre(0,HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax));
1320: PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,PETSC_NULL);
1321: PetscStackCallHypre(0,HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax));
1323: PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,PETSC_NULL);
1324: PetscStackCallHypre(0,HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol));
1325: PetscOptionsEList("-pc_syspfmg_relax_type","Relax type for the up and down cycles","HYPRE_SStructSysPFMGSetRelaxType",SysPFMGRelaxType,4,SysPFMGRelaxType[ex->relax_type],&ex->relax_type,PETSC_NULL);
1326: PetscStackCallHypre(0,HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type));
1327: PetscOptionsTail();
1328: return(0);
1329: }
1333: PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
1334: {
1335: PetscErrorCode ierr;
1336: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
1337: PetscScalar *xx,*yy;
1338: PetscInt ilower[3],iupper[3];
1339: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)(pc->pmat->data);
1340: PetscInt ordering= mx->dofs_order;
1341: PetscInt nvars= mx->nvars;
1342: PetscInt part= 0;
1343: PetscInt size;
1344: PetscInt i;
1347: DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
1348: iupper[0] += ilower[0] - 1;
1349: iupper[1] += ilower[1] - 1;
1350: iupper[2] += ilower[2] - 1;
1352: size= 1;
1353: for (i= 0; i< 3; i++) {
1354: size*= (iupper[i]-ilower[i]+1);
1355: }
1356: /* copy x values over to hypre for variable ordering */
1357: if (ordering) {
1358: PetscStackCallHypre(0,HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
1359: VecGetArray(x,&xx);
1360: for (i= 0; i< nvars; i++) {
1361: PetscStackCallHypre(0,HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,ilower,iupper,i,xx+(size*i)));
1362: }
1363: VecRestoreArray(x,&xx);
1364: PetscStackCallHypre(0,HYPRE_SStructVectorAssemble,(mx->ss_b));
1365: PetscStackCallHypre(0,HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
1366: PetscStackCallHypre(0,HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
1368: /* copy solution values back to PETSc */
1369: VecGetArray(y,&yy);
1370: for (i= 0; i< nvars; i++) {
1371: PetscStackCallHypre(0,HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,ilower,iupper,i,yy+(size*i)));
1372: }
1373: VecRestoreArray(y,&yy);
1374: } else { /* nodal ordering must be mapped to variable ordering for sys_pfmg */
1375: PetscScalar *z;
1376: PetscInt j, k;
1378: PetscMalloc(nvars*size*sizeof(PetscScalar),&z);
1379: PetscStackCallHypre(0,HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
1380: VecGetArray(x,&xx);
1382: /* transform nodal to hypre's variable ordering for sys_pfmg */
1383: for (i= 0; i< size; i++) {
1384: k= i*nvars;
1385: for (j= 0; j< nvars; j++) {
1386: z[j*size+i]= xx[k+j];
1387: }
1388: }
1389: for (i= 0; i< nvars; i++) {
1390: PetscStackCallHypre(0,HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,ilower,iupper,i,z+(size*i)));
1391: }
1392: VecRestoreArray(x,&xx);
1393: PetscStackCallHypre(0,HYPRE_SStructVectorAssemble,(mx->ss_b));
1394: PetscStackCallHypre(0,HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
1396: /* copy solution values back to PETSc */
1397: VecGetArray(y,&yy);
1398: for (i= 0; i< nvars; i++) {
1399: PetscStackCallHypre(0,HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,ilower,iupper,i,z+(size*i)));
1400: }
1401: /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
1402: for (i= 0; i< size; i++) {
1403: k= i*nvars;
1404: for (j= 0; j< nvars; j++) {
1405: yy[k+j]= z[j*size+i];
1406: }
1407: }
1408: VecRestoreArray(y,&yy);
1409: PetscFree(z);
1410: }
1411: return(0);
1412: }
1416: 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)
1417: {
1418: PC_SysPFMG *jac = (PC_SysPFMG*)pc->data;
1420: PetscInt oits;
1423: PetscStackCallHypre(0,HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its));
1424: PetscStackCallHypre(0,HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol));
1425: PCApply_SysPFMG(pc,b,y);
1426: PetscStackCallHypre(0,HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,&oits));
1427: *outits = oits;
1428: if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1429: else *reason = PCRICHARDSON_CONVERGED_RTOL;
1430: PetscStackCallHypre(0,HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol));
1431: PetscStackCallHypre(0,HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its));
1432: return(0);
1433: }
1438: PetscErrorCode PCSetUp_SysPFMG(PC pc)
1439: {
1440: PetscErrorCode ierr;
1441: PC_SysPFMG *ex = (PC_SysPFMG*) pc->data;
1442: Mat_HYPRESStruct *mx = (Mat_HYPRESStruct *)(pc->pmat->data);
1443: PetscBool flg;
1446: PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESSTRUCT,&flg);
1447: if (!flg) SETERRQ(((PetscObject)pc)->comm,PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESSTRUCT with this preconditioner");
1449: /* create the hypre sstruct solver object and set its information */
1450: if (ex->ss_solver) {
1451: PetscStackCallHypre(0,HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
1452: }
1453: PetscStackCallHypre(0,HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
1454: PCSetFromOptions_SysPFMG(pc);
1455: PetscStackCallHypre(0,HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver));
1456: PetscStackCallHypre(0,HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
1457: return(0);
1458: }
1461: /*MC
1462: PCSysPFMG - the hypre SysPFMG multigrid solver
1464: Level: advanced
1466: Options Database:
1467: + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner
1468: . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
1469: . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid
1470: . -pc_syspfmg_tol <tol> tolerance of SysPFMG
1471: . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel
1473: Notes: This is for CELL-centered descretizations
1475: This must be used with the MATHYPRESSTRUCT matrix type.
1476: This is less general than in hypre, it supports only one part, and one block per process defined by a PETSc DMDA.
1477: Also, only cell-centered variables.
1479: .seealso: PCMG, MATHYPRESSTRUCT
1480: M*/
1482: EXTERN_C_BEGIN
1485: PetscErrorCode PCCreate_SysPFMG(PC pc)
1486: {
1487: PetscErrorCode ierr;
1488: PC_SysPFMG *ex;
1491: PetscNew(PC_SysPFMG,&ex);\
1492: pc->data = ex;
1494: ex->its = 1;
1495: ex->tol = 1.e-8;
1496: ex->relax_type = 1;
1497: ex->num_pre_relax = 1;
1498: ex->num_post_relax = 1;
1500: pc->ops->setfromoptions = PCSetFromOptions_SysPFMG;
1501: pc->ops->view = PCView_SysPFMG;
1502: pc->ops->destroy = PCDestroy_SysPFMG;
1503: pc->ops->apply = PCApply_SysPFMG;
1504: pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
1505: pc->ops->setup = PCSetUp_SysPFMG;
1506: MPI_Comm_dup(((PetscObject)pc)->comm,&(ex->hcomm));
1507: PetscStackCallHypre(0,HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
1508: return(0);
1509: }
1510: EXTERN_C_END