Actual source code: hyppilut.c

  1: #define PETSCKSP_DLL

  3: /*

  5: */

 7:  #include private/pcimpl.h
  9: #include "HYPRE.h"
 10: #include "IJ_mv.h"
 11: #include "parcsr_ls.h"

 14: EXTERN PetscErrorCode MatHYPRE_IJMatrixCreate(Mat,HYPRE_IJMatrix*);
 15: EXTERN PetscErrorCode MatHYPRE_IJMatrixCopy(Mat,HYPRE_IJMatrix);
 16: EXTERN PetscErrorCode VecHYPRE_IJVectorCreate(Vec,HYPRE_IJVector*);

 18: /* 
 19:    Private context (data structure) for the  preconditioner.  
 20: */
 21: typedef struct {
 22:   HYPRE_Solver       hsolver;
 23:   HYPRE_IJMatrix     ij;
 24:   HYPRE_IJVector     b,x;

 26:   PetscErrorCode     (*destroy)(HYPRE_Solver);
 27:   PetscErrorCode     (*solve)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
 28:   PetscErrorCode     (*setup)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
 29: 
 30:   MPI_Comm           comm_hypre;
 31:   char              *hypre_type;

 33:   /* options for pilut and BoomerAMG*/
 34:   int                maxiter;
 35:   double             tol;
 36:   PetscTruth         applyrichardson;

 38:   /* options for pilut */
 39:   int                factorrowsize;

 41:   /* options for parasails */
 42:   int                nlevels;
 43:   double             threshhold;
 44:   double             filter;
 45:   int                sym;
 46:   double             loadbal;
 47:   int                logging;
 48:   int                ruse;
 49:   int                symt;

 51:   /* options for euclid */
 52:   PetscTruth         bjilu;
 53:   int                levels;

 55:   /* options for euclid and BoomerAMG */
 56:   PetscTruth         printstatistics;

 58:   /* options for BoomerAMG */
 59:   int                maxlevels;
 60:   double             strongthreshold;
 61:   double             maxrowsum;
 62:   int                gridsweeps[4];
 63:   int                coarsentype;
 64:   int                measuretype;
 65:   int                relaxtype[4];
 66:   double             relaxweight;
 67:   double             outerrelaxweight;
 68:   int                relaxorder;
 69:   int                **gridrelaxpoints;
 70:   double             truncfactor;
 71: } PC_HYPRE;


 76: static PetscErrorCode PCSetUp_HYPRE(PC pc)
 77: {
 78:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
 79:   PetscErrorCode     ierr;
 80:   HYPRE_ParCSRMatrix hmat;
 81:   HYPRE_ParVector    bv,xv;
 82:   int                hierr;

 85:   if (!jac->ij) { /* create the matrix the first time through */
 86:     MatHYPRE_IJMatrixCreate(pc->pmat,&jac->ij);
 87:   }
 88:   if (!jac->b) {
 89:     Vec vec;
 90:     MatGetVecs(pc->pmat,&vec,0);
 91:     VecHYPRE_IJVectorCreate(vec,&jac->b);
 92:     VecHYPRE_IJVectorCreate(vec,&jac->x);
 93:     VecDestroy(vec);
 94:   }
 95:   MatHYPRE_IJMatrixCopy(pc->pmat,jac->ij);
 96:   HYPRE_IJMatrixGetObject(jac->ij,(void**)&hmat);
 97:   HYPRE_IJVectorGetObject(jac->b,(void**)&bv);
 98:   HYPRE_IJVectorGetObject(jac->x,(void**)&xv);
 99:   h(*jac->setup)(jac->hsolver,hmat,bv,xv);
100:   if (hierr) SETERRQ1(PETSC_ERR_LIB,"Error in HYPRE setup, error code %d",hierr);
101:   return(0);
102: }

104: /*
105:     Replaces the address where the HYPRE vector points to its data with the address of
106:   PETSc's data. Saves the old address so it can be reset when we are finished with it.
107:   Allows use to get the data into a HYPRE vector without the cost of memcopies 
108: */
109: #define HYPREReplacePointer(b,newvalue,savedvalue) {\
110:    hypre_ParVector *par_vector   = (hypre_ParVector *)hypre_IJVectorObject(((hypre_IJVector*)b));\
111:    hypre_Vector    *local_vector = hypre_ParVectorLocalVector(par_vector);\
112:    savedvalue         = local_vector->data;\
113:    local_vector->data = newvalue;}

117: static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
118: {
119:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
120:   PetscErrorCode     ierr;
121:   HYPRE_ParCSRMatrix hmat;
122:   PetscScalar        *bv,*xv;
123:   HYPRE_ParVector    jbv,jxv;
124:   PetscScalar        *sbv,*sxv;
125:   int                hierr;

128:   if (!jac->applyrichardson) {VecSet(x,0.0);}
129:   VecGetArray(b,&bv);
130:   VecGetArray(x,&xv);
131:   HYPREReplacePointer(jac->b,bv,sbv);
132:   HYPREReplacePointer(jac->x,xv,sxv);

134:   HYPRE_IJMatrixGetObject(jac->ij,(void**)&hmat);
135:   HYPRE_IJVectorGetObject(jac->b,(void**)&jbv);
136:   HYPRE_IJVectorGetObject(jac->x,(void**)&jxv);
137:   h(*jac->solve)(jac->hsolver,hmat,jbv,jxv);
138:   /* error code of 1 in boomerAMG merely means convergence not achieved */
139:   if (hierr && (hierr != 1 || jac->solve != HYPRE_BoomerAMGSolve)) SETERRQ1(PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
140:   if (hierr) hypre__global_error = 0;
141: 
142:   HYPREReplacePointer(jac->b,sbv,bv);
143:   HYPREReplacePointer(jac->x,sxv,xv);
144:   VecRestoreArray(x,&xv);
145:   VecRestoreArray(b,&bv);
146:   return(0);
147: }

151: static PetscErrorCode PCDestroy_HYPRE(PC pc)
152: {
153:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;

157:   if (jac->ij) { HYPRE_IJMatrixDestroy(jac->ij); }
158:   if (jac->b) { HYPRE_IJVectorDestroy(jac->b); }
159:   if (jac->x) { HYPRE_IJVectorDestroy(jac->x); }
160:   PetscStrfree(jac->hypre_type);
161:   (*jac->destroy)(jac->hsolver);
162:   MPI_Comm_free(&(jac->comm_hypre));
163:   PetscFree(jac);
164:   return(0);
165: }

167: /* --------------------------------------------------------------------------------------------*/
170: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PC pc)
171: {
172:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
174:   PetscTruth     flag;

177:   PetscOptionsHead("HYPRE Pilut Options");
178:   PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
179:   if (flag) {
180:     HYPRE_ParCSRPilutSetMaxIter(jac->hsolver,jac->maxiter);
181:   }
182:   PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
183:   if (flag) {
184:     HYPRE_ParCSRPilutSetDropTolerance(jac->hsolver,jac->tol);
185:   }
186:   PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
187:   if (flag) {
188:     HYPRE_ParCSRPilutSetFactorRowSize(jac->hsolver,jac->factorrowsize);
189:   }
190:   PetscOptionsTail();
191:   return(0);
192: }

196: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
197: {
198:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
200:   PetscTruth     iascii;

203:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
204:   if (iascii) {
205:     PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut preconditioning\n");
206:     if (jac->maxiter != PETSC_DEFAULT) {
207:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: maximum number of iterations %d\n",jac->maxiter);
208:     } else {
209:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default maximum number of iterations \n");
210:     }
211:     if (jac->tol != PETSC_DEFAULT) {
212:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: drop tolerance %G\n",jac->tol);
213:     } else {
214:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default drop tolerance \n");
215:     }
216:     if (jac->factorrowsize != PETSC_DEFAULT) {
217:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: factor row size %d\n",jac->factorrowsize);
218:     } else {
219:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default factor row size \n");
220:     }
221:   }
222:   return(0);
223: }

225: /* --------------------------------------------------------------------------------------------*/
228: static PetscErrorCode PCSetFromOptions_HYPRE_Euclid(PC pc)
229: {
230:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
232:   PetscTruth     flag;
233:   char           *args[2];

236:   PetscOptionsHead("HYPRE Euclid Options");
237:   PetscOptionsInt("-pc_hypre_euclid_levels","Number of levels of fill ILU(k)","None",jac->levels,&jac->levels,&flag);
238:   if (flag) {
239:     char levels[16];
240:     if (jac->levels < 0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be nonegative",jac->levels);
241:     sprintf(levels,"%d",jac->levels);
242:     args[0] = (char*)"-level"; args[1] = levels;
243:     HYPRE_EuclidSetParams(jac->hsolver,2,args);
244:   }
245:   PetscOptionsTruth("-pc_hypre_euclid_bj","Use block Jacobi ILU(k)","None",jac->bjilu,&jac->bjilu,PETSC_NULL);
246:   if (jac->bjilu) {
247:     args[0] =(char*) "-bj"; args[1] = (char*)"1";
248:     HYPRE_EuclidSetParams(jac->hsolver,2,args);
249:   }
250: 
251:   PetscOptionsTruth("-pc_hypre_euclid_print_statistics","Print statistics","None",jac->printstatistics,&jac->printstatistics,PETSC_NULL);
252:   if (jac->printstatistics) {
253:     args[0] = (char*)"-eu_stats"; args[1] = (char*)"1";
254:     HYPRE_EuclidSetParams(jac->hsolver,2,args);
255:     args[0] = (char*)"-eu_mem"; args[1] = (char*)"1";
256:     HYPRE_EuclidSetParams(jac->hsolver,2,args);
257:   }
258:   PetscOptionsTail();
259:   return(0);
260: }

264: static PetscErrorCode PCView_HYPRE_Euclid(PC pc,PetscViewer viewer)
265: {
266:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
268:   PetscTruth     iascii;

271:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
272:   if (iascii) {
273:     PetscViewerASCIIPrintf(viewer,"  HYPRE Euclid preconditioning\n");
274:     PetscViewerASCIIPrintf(viewer,"  HYPRE Euclid: number of levels %d\n",jac->levels);
275:     if (jac->bjilu) {
276:       PetscViewerASCIIPrintf(viewer,"  HYPRE Euclid: Using block Jacobi ILU instead of parallel ILU\n");
277:     }
278:   }
279:   return(0);
280: }

282: /* --------------------------------------------------------------------------------------------*/

284: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout"};
285: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
286: static const char *HYPREBoomerAMGRelaxType[]   = {"Jacobi","sequential-Gauss-Seidel","","SOR/Jacobi","backward-SOR/Jacobi","","symmetric-SOR/Jacobi",
287:                                                   "","","Gaussian-elimination"};
290: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PC pc)
291: {
292:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
294:   int            n,indx;
295:   PetscTruth     flg, tmp_truth;
296:   double         tmpdbl, twodbl[2];

299:   PetscOptionsHead("HYPRE BoomerAMG Options");
300:   PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
301:   if (flg) {
302:     if (jac->maxlevels < 2) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
303:     HYPRE_BoomerAMGSetMaxLevels(jac->hsolver,jac->maxlevels);
304:   }
305:   PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
306:   if (flg) {
307:     if (jac->maxiter < 1) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
308:     HYPRE_BoomerAMGSetMaxIter(jac->hsolver,jac->maxiter);
309:   }
310:   PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call","None",jac->tol,&jac->tol,&flg);
311:   if (flg) {
312:     if (jac->tol < 0.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Tolerance %G must be great than or equal zero",jac->tol);
313:     HYPRE_BoomerAMGSetTol(jac->hsolver,jac->tol);
314:   }

316:   PetscOptionsScalar("-pc_hypre_boomeramg_truncfactor","Truncation factor","None",jac->truncfactor,&jac->truncfactor,&flg);
317:   if (flg) {
318:     if (jac->truncfactor < 0.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Truncation factor %G must be great than or equal zero",jac->truncfactor);
319:     HYPRE_BoomerAMGSetTruncFactor(jac->hsolver,jac->truncfactor);
320:   }

322:   PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
323:   if (flg) {
324:     if (jac->strongthreshold < 0.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Strong threshold %G must be great than or equal zero",jac->strongthreshold);
325:     HYPRE_BoomerAMGSetStrongThreshold(jac->hsolver,jac->strongthreshold);
326:   }
327:   PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
328:   if (flg) {
329:     if (jac->maxrowsum < 0.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %G must be greater than zero",jac->maxrowsum);
330:     if (jac->maxrowsum > 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %G must be less than or equal one",jac->maxrowsum);
331:     HYPRE_BoomerAMGSetMaxRowSum(jac->hsolver,jac->maxrowsum);
332:   }

334:   /* Grid sweeps */
335:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for all grid levels (fine, up, and down)","None", jac->gridsweeps[0], &indx ,&flg);
336:   if (flg) {
337:     HYPRE_BoomerAMGSetNumSweeps(jac->hsolver,indx);
338:     /* modify the jac structure so we can view the updated options with PC_View */
339:     jac->gridsweeps[0] = indx;
340:     jac->gridsweeps[1] = jac->gridsweeps[2] = jac->gridsweeps[0];
341:     jac->gridsweeps[3] = 1;  /*The coarse level is not affected by this function - hypre code sets to 1*/
342:   }
343:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_fine","Number of sweeps for the fine level","None", jac->gridsweeps[0], &indx ,&flg);
344:   if (flg) {
345:     HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver,indx, 0);
346:     jac->gridsweeps[0] = indx;
347:   }
348:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[2], &indx ,&flg);
349:   if (flg) {
350:     HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver,indx, 1);
351:     jac->gridsweeps[1] = indx;
352:   }
353:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None", jac->gridsweeps[1], &indx ,&flg);
354:   if (flg) {
355:     HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver,indx, 2);
356:     jac->gridsweeps[2] = indx;
357:   }
358:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None", jac->gridsweeps[3], &indx ,&flg);
359:   if (flg) {
360:     HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver,indx, 3);
361:     jac->gridsweeps[3] = indx;
362:   }

364:   /* Relax type */
365:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for fine, up, and down cycles (coarse level set to gaussian elimination)","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[3],&indx,&flg);
366:   if (flg) {
367:     jac->relaxtype[0] = jac->relaxtype[1] = jac->relaxtype[2] = indx;
368:     jac->relaxtype[3] = 9; /* hypre code sets coarse grid to 9 (G.E.)*/
369:     hypre_BoomerAMGSetRelaxType(jac->hsolver, indx);
370:   }
371:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_fine","Relax type on fine grid","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[3],&indx,&flg);
372:   if (flg) {
373:     jac->relaxtype[0] = indx;
374:     HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, indx, 0);
375:   }
376:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[3],&indx,&flg);
377:   if (flg) {
378:     jac->relaxtype[1] = indx;
379:     HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, indx, 1);
380:   }
381:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[3],&indx,&flg);
382:   if (flg) {
383:     jac->relaxtype[2] = indx;
384:     HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, indx, 2);
385:   }
386:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[9],&indx,&flg);
387:   if (flg) {
388:     jac->relaxtype[3] = indx;
389:     HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, indx, 3);
390:   }

392:   /* Relaxation Weight */
393:   PetscOptionsReal("-pc_hypre_boomeramg_relax_weight_all","Relaxation weight for all levels (0 = hypre estimates, -k = determined with k CG steps)","None",jac->relaxweight, &tmpdbl ,&flg);
394:   if (flg) {
395:     hypre_BoomerAMGSetRelaxWt( jac->hsolver, tmpdbl);
396:     jac->relaxweight = tmpdbl;
397:   }

399:   n=2;
400:   twodbl[0] = twodbl[1] = 1.0;
401:   PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
402:   if (flg) {
403:     if (n == 2) {
404:       indx =  (int)PetscAbsReal(twodbl[1]);
405:       hypre_BoomerAMGSetLevelRelaxWt( jac->hsolver, twodbl[0], indx);
406:     } else {
407:       SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Relax weight level: you must provide 2 values separated by a comma (and no space), you provided %d",n);
408:     }
409:   }

411:   /* Outer relaxation Weight */
412:   PetscOptionsReal("-pc_hypre_boomeramg_outer_relax_weight_all","Outer relaxation weight for all levels ( -k = determined with k CG steps)","None",jac->outerrelaxweight, &tmpdbl ,&flg);
413:   if (flg) {
414:     hypre_BoomerAMGSetOuterWt( jac->hsolver, tmpdbl);
415:     jac->outerrelaxweight = tmpdbl;
416:   }

418:   n=2;
419:   twodbl[0] = twodbl[1] = 1.0;
420:   PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
421:   if (flg) {
422:     if (n == 2) {
423:       indx =  (int)PetscAbsReal(twodbl[1]);
424:       hypre_BoomerAMGSetLevelOuterWt( jac->hsolver, twodbl[0], indx);
425:     } else {
426:       SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Relax weight outer level: You must provide 2 values separated by a comma (and no space), you provided %d",n);
427:     }
428:   }

430:   /* the Relax Order */
431:   /* PetscOptionsName("-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", &flg); */
432:   PetscOptionsTruth( "-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", PETSC_FALSE, &tmp_truth, &flg);

434:   if (flg) {
435:     jac->relaxorder = 0;
436:     hypre_BoomerAMGSetRelaxOrder( jac->hsolver, jac->relaxorder);
437:   }
438:   PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,2,HYPREBoomerAMGMeasureType[0],&indx,&flg);
439:   if (flg) {
440:     jac->measuretype = indx;
441:     HYPRE_BoomerAMGSetMeasureType(jac->hsolver,jac->measuretype);
442:   }
443:   PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,7,HYPREBoomerAMGCoarsenType[6],&indx,&flg);
444:   if (flg) {
445:     jac->coarsentype = indx;
446:     HYPRE_BoomerAMGSetCoarsenType(jac->hsolver,jac->coarsentype);
447:   }
448:   PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
449:   if (flg) {
450:     int level=3;
451:     jac->printstatistics = PETSC_TRUE;
452:     PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,PETSC_NULL);
453:     HYPRE_BoomerAMGSetPrintLevel(jac->hsolver,level);
454:     HYPRE_BoomerAMGSetDebugFlag(jac->hsolver,level);
455:   }
456:   PetscOptionsTail();
457:   return(0);
458: }

462: static PetscErrorCode PCApplyRichardson_BoomerAMG(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its)
463: {
464:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;

468:   HYPRE_BoomerAMGSetMaxIter(jac->hsolver,its*jac->maxiter);
469:   HYPRE_BoomerAMGSetTol(jac->hsolver,PetscMin(rtol,jac->tol));
470:   jac->applyrichardson = PETSC_TRUE;
471:   PCApply_HYPRE(pc,b,y);
472:   jac->applyrichardson = PETSC_FALSE;
473:   HYPRE_BoomerAMGSetTol(jac->hsolver,jac->tol);
474:   HYPRE_BoomerAMGSetMaxIter(jac->hsolver,jac->maxiter);
475:   return(0);
476: }


481: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
482: {
483:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
485:   PetscTruth     iascii;

488:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
489:   if (iascii) {
490:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG preconditioning\n");
491:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of levels %d\n",jac->maxlevels);
492:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of iterations PER hypre call%d\n",jac->maxiter);
493:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Convergence tolerance PER hypre call%G\n",jac->tol);
494:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Threshold for strong coupling %G\n",jac->strongthreshold);
495:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum row sums %G\n",jac->maxrowsum);

497:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps on fine grid %d\n",jac->gridsweeps[0]);
498:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps down         %d\n",jac->gridsweeps[1]);
499:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps up           %d\n",jac->gridsweeps[2]);
500:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps on coarse    %d\n",jac->gridsweeps[3]);

502:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax on fine grid  %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[0]]);
503:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax down          %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[1]]);
504:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax up            %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[2]]);
505:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax on coarse     %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[3]]);

507:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax weight  (all)      %G\n",jac->relaxweight);
508:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Outer relax weight (all) %G\n",jac->outerrelaxweight);

510:     if (jac->relaxorder) {
511:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Using CF-relaxation\n");
512:     } else {
513:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Not using CF-relaxation\n");
514:     }
515:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Measure type        %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
516:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Coarsen type        %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
517:   }
518:   return(0);
519: }

521: /* --------------------------------------------------------------------------------------------*/
524: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PC pc)
525: {
526:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
528:   int            indx;
529:   PetscTruth     flag;
530:   const char     *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};

533:   PetscOptionsHead("HYPRE ParaSails Options");
534:   PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);
535:   PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshhold,&jac->threshhold,&flag);
536:   if (flag) {
537:     HYPRE_ParaSailsSetParams(jac->hsolver,jac->threshhold,jac->nlevels);
538:   }

540:   PetscOptionsReal("-pc_hypre_parasails_filter","filter","None",jac->filter,&jac->filter,&flag);
541:   if (flag) {
542:     HYPRE_ParaSailsSetFilter(jac->hsolver,jac->filter);
543:   }

545:   PetscOptionsReal("-pc_hypre_parasails_loadbal","Load balance","None",jac->loadbal,&jac->loadbal,&flag);
546:   if (flag) {
547:     HYPRE_ParaSailsSetLoadbal(jac->hsolver,jac->loadbal);
548:   }

550:   PetscOptionsTruth("-pc_hypre_parasails_logging","Print info to screen","None",(PetscTruth)jac->logging,(PetscTruth*)&jac->logging,&flag);
551:   if (flag) {
552:     HYPRE_ParaSailsSetLogging(jac->hsolver,jac->logging);
553:   }

555:   PetscOptionsTruth("-pc_hypre_parasails_reuse","Reuse nonzero pattern in preconditioner","None",(PetscTruth)jac->ruse,(PetscTruth*)&jac->ruse,&flag);
556:   if (flag) {
557:     HYPRE_ParaSailsSetReuse(jac->hsolver,jac->ruse);
558:   }

560:   PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,3,symtlist[0],&indx,&flag);
561:   if (flag) {
562:     jac->symt = indx;
563:     HYPRE_ParaSailsSetSym(jac->hsolver,jac->symt);
564:   }

566:   PetscOptionsTail();
567:   return(0);
568: }

572: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
573: {
574:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
576:   PetscTruth     iascii;
577:   const char     *symt = 0;;

580:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
581:   if (iascii) {
582:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails preconditioning\n");
583:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: nlevels %d\n",jac->nlevels);
584:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: threshold %G\n",jac->threshhold);
585:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: filter %G\n",jac->filter);
586:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: load balance %G\n",jac->loadbal);
587:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: reuse nonzero structure %s\n",PetscTruths[jac->ruse]);
588:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: print info to screen %s\n",PetscTruths[jac->logging]);
589:     if (!jac->symt) {
590:       symt = "nonsymmetric matrix and preconditioner";
591:     } else if (jac->symt == 1) {
592:       symt = "SPD matrix and preconditioner";
593:     } else if (jac->symt == 2) {
594:       symt = "nonsymmetric matrix but SPD preconditioner";
595:     } else {
596:       SETERRQ1(PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
597:     }
598:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: %s\n",symt);
599:   }
600:   return(0);
601: }
602: /* ---------------------------------------------------------------------------------*/

607: PetscErrorCode  PCHYPREGetType_HYPRE(PC pc,const char *name[])
608: {
609:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;

612:   *name = jac->hypre_type;
613:   return(0);
614: }

620: PetscErrorCode  PCHYPRESetType_HYPRE(PC pc,const char name[])
621: {
622:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
624:   PetscTruth     flag;
625:   PetscInt       bs;

628:   if (pc->ops->setup) {
629:     SETERRQ(PETSC_ERR_ORDER,"Cannot set the HYPRE preconditioner type once it has been set");
630:   }

632:   pc->ops->setup          = PCSetUp_HYPRE;
633:   pc->ops->apply          = PCApply_HYPRE;
634:   pc->ops->destroy        = PCDestroy_HYPRE;

636:   jac->maxiter            = PETSC_DEFAULT;
637:   jac->tol                = PETSC_DEFAULT;
638:   jac->printstatistics    = PetscLogPrintInfo;

640:   PetscStrallocpy(name, &jac->hypre_type);
641:   PetscStrcmp("pilut",name,&flag);
642:   if (flag) {
643:     HYPRE_ParCSRPilutCreate(jac->comm_hypre,&jac->hsolver);
644:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
645:     pc->ops->view           = PCView_HYPRE_Pilut;
646:     jac->destroy            = HYPRE_ParCSRPilutDestroy;
647:     jac->setup              = HYPRE_ParCSRPilutSetup;
648:     jac->solve              = HYPRE_ParCSRPilutSolve;
649:     jac->factorrowsize      = PETSC_DEFAULT;
650:     return(0);
651:   }
652:   PetscStrcmp("parasails",name,&flag);
653:   if (flag) {
654:     HYPRE_ParaSailsCreate(jac->comm_hypre,&jac->hsolver);
655:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
656:     pc->ops->view           = PCView_HYPRE_ParaSails;
657:     jac->destroy            = HYPRE_ParaSailsDestroy;
658:     jac->setup              = HYPRE_ParaSailsSetup;
659:     jac->solve              = HYPRE_ParaSailsSolve;
660:     /* initialize */
661:     jac->nlevels     = 1;
662:     jac->threshhold  = .1;
663:     jac->filter      = .1;
664:     jac->loadbal     = 0;
665:     if (PetscLogPrintInfo) {
666:       jac->logging     = (int) PETSC_TRUE;
667:     } else {
668:       jac->logging     = (int) PETSC_FALSE;
669:     }
670:     jac->ruse = (int) PETSC_FALSE;
671:     jac->symt   = 0;
672:     HYPRE_ParaSailsSetParams(jac->hsolver,jac->threshhold,jac->nlevels);
673:     HYPRE_ParaSailsSetFilter(jac->hsolver,jac->filter);
674:     HYPRE_ParaSailsSetLoadbal(jac->hsolver,jac->loadbal);
675:     HYPRE_ParaSailsSetLogging(jac->hsolver,jac->logging);
676:     HYPRE_ParaSailsSetReuse(jac->hsolver,jac->ruse);
677:     HYPRE_ParaSailsSetSym(jac->hsolver,jac->symt);
678:     return(0);
679:   }
680:   PetscStrcmp("euclid",name,&flag);
681:   if (flag) {
682:     HYPRE_EuclidCreate(jac->comm_hypre,&jac->hsolver);
683:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
684:     pc->ops->view           = PCView_HYPRE_Euclid;
685:     jac->destroy            = HYPRE_EuclidDestroy;
686:     jac->setup              = HYPRE_EuclidSetup;
687:     jac->solve              = HYPRE_EuclidSolve;
688:     /* initialization */
689:     jac->bjilu              = PETSC_FALSE;
690:     jac->levels             = 1;
691:     return(0);
692:   }
693:   PetscStrcmp("boomeramg",name,&flag);
694:   if (flag) {
695:     HYPRE_BoomerAMGCreate(&jac->hsolver);
696:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_BoomerAMG;
697:     pc->ops->view            = PCView_HYPRE_BoomerAMG;
698:     jac->destroy             = HYPRE_BoomerAMGDestroy;
699:     jac->setup               = HYPRE_BoomerAMGSetup;
700:     jac->solve               = HYPRE_BoomerAMGSolve;
701:     pc->ops->applyrichardson = PCApplyRichardson_BoomerAMG;
702:     /* these defaults match the hypre defaults */
703:     jac->maxlevels        = 25;
704:     jac->maxiter          = 1;
705:     jac->tol              = 1.e-7;
706:     jac->strongthreshold  = .25;
707:     jac->maxrowsum        = .9;
708:     jac->coarsentype      = 6;
709:     jac->measuretype      = 0;
710:     jac->applyrichardson  = PETSC_FALSE;
711:     jac->gridsweeps[0]    = jac->gridsweeps[1] = jac->gridsweeps[2] = jac->gridsweeps[3]  = 1;
712:     jac->relaxtype[0]     = jac->relaxtype[1]  = jac->relaxtype[2]  = 3;
713:     jac->relaxtype[3]     = 9;
714:     jac->relaxweight      = 1.0;
715:     jac->outerrelaxweight = 1.0;
716:     jac->relaxorder       = 1;
717:     HYPRE_BoomerAMGSetMaxLevels(jac->hsolver,jac->maxlevels);
718:     HYPRE_BoomerAMGSetMaxIter(jac->hsolver,jac->maxiter);
719:     HYPRE_BoomerAMGSetTol(jac->hsolver,jac->tol);
720:     HYPRE_BoomerAMGSetTruncFactor(jac->hsolver,jac->truncfactor);
721:     HYPRE_BoomerAMGSetStrongThreshold(jac->hsolver,jac->strongthreshold);
722:     HYPRE_BoomerAMGSetMaxRowSum(jac->hsolver,jac->maxrowsum);
723:     HYPRE_BoomerAMGSetMeasureType(jac->hsolver,jac->measuretype);
724:     HYPRE_BoomerAMGSetCoarsenType(jac->hsolver,jac->coarsentype);
725:     MatGetBlockSize(pc->pmat,&bs);
726:     HYPRE_BoomerAMGSetNumFunctions(jac->hsolver,bs);
727:     return(0);
728:   }
729:   SETERRQ1(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are pilut, parasails, euclid, boomeramg",name);
730:   return(0);
731: }

734: /*
735:     It only gets here if the HYPRE type has not been set before the call to 
736:    ...SetFromOptions() which actually is most of the time
737: */
740: static PetscErrorCode PCSetFromOptions_HYPRE(PC pc)
741: {
743:   int            indx;
744:   const char     *type[] = {"pilut","parasails","boomeramg","euclid"};
745:   PetscTruth     flg;

748:   PetscOptionsHead("HYPRE preconditioner options");
749:   PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"pilut",&indx,&flg);
750:   if (PetscOptionsPublishCount) {   /* force the default if it was not yet set and user did not set with option */
751:     if (!flg && !pc->ops->apply) {
752:       flg   = PETSC_TRUE;
753:       indx = 0;
754:     }
755:   }
756:   if (flg) {
757:     PCHYPRESetType_HYPRE(pc,type[indx]);
758:   }
759:   if (pc->ops->setfromoptions) {
760:     pc->ops->setfromoptions(pc);
761:   }
762:   PetscOptionsTail();
763:   return(0);
764: }

768: /*@C
769:      PCHYPRESetType - Sets which hypre preconditioner you wish to use

771:    Input Parameters:
772: +     pc - the preconditioner context
773: -     name - either  pilut, parasails, boomeramg, euclid

775:    Options Database Keys:
776:    -pc_hypre_type - One of pilut, parasails, boomeramg, euclid
777:  
778:    Level: intermediate

780: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
781:            PCHYPRE

783: @*/
784: PetscErrorCode  PCHYPRESetType(PC pc,const char name[])
785: {
786:   PetscErrorCode ierr,(*f)(PC,const char[]);

791:   PetscObjectQueryFunction((PetscObject)pc,"PCHYPRESetType_C",(void (**)(void))&f);
792:   if (f) {
793:     (*f)(pc,name);
794:   }
795:   return(0);
796: }

800: /*@C
801:      PCHYPREGetType - Gets which hypre preconditioner you are using

803:    Input Parameter:
804: .     pc - the preconditioner context

806:    Output Parameter:
807: .     name - either  pilut, parasails, boomeramg, euclid

809:    Level: intermediate

811: .seealso:  PCCreate(), PCHYPRESetType(), PCType (for list of available types), PC,
812:            PCHYPRE

814: @*/
815: PetscErrorCode  PCHYPREGetType(PC pc,const char *name[])
816: {
817:   PetscErrorCode ierr,(*f)(PC,const char*[]);

822:   PetscObjectQueryFunction((PetscObject)pc,"PCHYPREGetType_C",(void (**)(void))&f);
823:   if (f) {
824:     (*f)(pc,name);
825:   }
826:   return(0);
827: }

829: /*MC
830:      PCHYPRE - Allows you to use the matrix element based preconditioners in the LLNL package hypre

832:    Options Database Keys:
833: +   -pc_hypre_type - One of pilut, parasails, boomeramg, euclid
834: -   Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX
835:           preconditioner
836:  
837:    Level: intermediate

839:    Notes: Apart from pc_hypre_type (for which there is PCHYPRESetType()),
840:           the many hypre options can ONLY be set via the options database (e.g. the command line
841:           or with PetscOptionsSetValue(), there are no functions to set them)

843:           The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_rtol refer to the number of iterations
844:           (V-cycles) that boomeramg does EACH time it is called. So for example, if it is set to 2 then 
845:           2-V-cycles are being used to define the preconditioner. -ksp_max_iter and -ksp_rtol STILL determine
846:           the total number of iterations and tolerance for the Krylov solver. For example, if 
847:           -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10 then AT MOST twenty V-cycles of boomeramg
848:           will be called.


851:           If you wish to use boomerAMG WITHOUT a Krylov method use -ksp_type richardson NOT -ksp_type preonly
852:           and use -ksp_max_it to control the number of V-cycles.
853:           (see the PETSc FAQ.html at the PETSc website under the Documentation tab).

855: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
856:            PCHYPRESetType()

858: M*/

863: PetscErrorCode  PCCreate_HYPRE(PC pc)
864: {
865:   PC_HYPRE       *jac;

869:   PetscNew(PC_HYPRE,&jac);
870:   pc->data                 = jac;
871:   pc->ops->setfromoptions  = PCSetFromOptions_HYPRE;
872:   jac->hypre_type          = PETSC_NULL;
873:   /* Com_dup for hypre */
874:   MPI_Comm_dup(pc->comm,&(jac->comm_hypre));
875:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCHYPRESetType_C","PCHYPRESetType_HYPRE",PCHYPRESetType_HYPRE);
876:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCHYPREGetType_C","PCHYPREGetType_HYPRE",PCHYPREGetType_HYPRE);
877:   return(0);
878: }