Actual source code: hyppilut.c

      1: #define PETSCKSP_DLL

      3: /*

      5: */

     7:  #include src/ksp/pc/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;

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

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

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

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

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

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

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

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

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

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

    128:   if (!jac->applyrichardson) {VecSet(x,zero);}
    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:
    141:   HYPREReplacePointer(jac->b,sbv,bv);
    142:   HYPREReplacePointer(jac->x,sxv,xv);
    143:   VecRestoreArray(x,&xv);
    144:   VecRestoreArray(b,&bv);
    145:   return(0);
    146: }

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

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

    165: /* --------------------------------------------------------------------------------------------*/
    168: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PC pc)
    169: {
    170:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
    172:   PetscTruth     flag;

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

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

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

    223: /* --------------------------------------------------------------------------------------------*/
    226: static PetscErrorCode PCSetFromOptions_HYPRE_Euclid(PC pc)
    227: {
    228:   PC_HYPRE  *jac = (PC_HYPRE*)pc->data;
    230:   PetscTruth flag;
    231:   char       *args[2];

    234:   jac->bjilu              = PETSC_FALSE;
    235:   jac->levels             = 1;

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

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

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

    283: /* --------------------------------------------------------------------------------------------*/

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

    300:   /* these defaults match the hypre defaults */
    301:   jac->maxlevels        = 25;
    302:   jac->maxiter          = 1;
    303:   jac->tol              = 1.e-7;
    304:   jac->strongthreshold  = .25;
    305:   jac->maxrowsum        = .9;
    306:   jac->coarsentype      = 6;
    307:   jac->measuretype      = 0;
    308:   jac->applyrichardson  = PETSC_FALSE;
    309:   jac->gridsweeps[0]    = jac->gridsweeps[1] = jac->gridsweeps[2] = jac->gridsweeps[3]  = 1;
    310:   jac->relaxtype[0]     = jac->relaxtype[1]  = jac->relaxtype[2]  = 3;
    311:   jac->relaxtype[3]     = 9;
    312:   jac->relaxweight      = 1.0;
    313:   jac->outerrelaxweight = 1.0;
    314:   jac->relaxorder       = 1;
    315:
    316:   PetscOptionsHead("HYPRE BoomerAMG Options");
    317:   PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
    318:   if (flg) {
    319:     if (jac->maxlevels < 2) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
    320:   }
    321:   HYPRE_BoomerAMGSetMaxLevels(jac->hsolver,jac->maxlevels);
    322:   PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used","None",jac->maxiter,&jac->maxiter,&flg);
    323:   if (flg) {
    324:     if (jac->maxiter < 1) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
    325:   }
    326:   HYPRE_BoomerAMGSetMaxIter(jac->hsolver,jac->maxiter);
    327:   PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance","None",jac->tol,&jac->tol,&flg);
    328:   if (flg) {
    329:     if (jac->tol < 0.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Tolerance %g must be great than or equal zero",jac->tol);
    330:   }
    331:   HYPRE_BoomerAMGSetTol(jac->hsolver,jac->tol);

    333:   PetscOptionsScalar("-pc_hypre_boomeramg_truncfactor","Truncation factor","None",jac->truncfactor,&jac->truncfactor,&flg);
    334:   if (flg) {
    335:     if (jac->truncfactor < 0.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Truncation factor %g must be great than or equal zero",jac->truncfactor);
    336:   }
    337:   HYPRE_BoomerAMGSetTruncFactor(jac->hsolver,jac->truncfactor);

    339:   PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
    340:   if (flg) {
    341:     if (jac->strongthreshold < 0.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Strong threshold %g must be great than or equal zero",jac->strongthreshold);
    342:   }
    343:   HYPRE_BoomerAMGSetStrongThreshold(jac->hsolver,jac->strongthreshold);
    344:   PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
    345:   if (flg) {
    346:     if (jac->maxrowsum < 0.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %g must be greater than zero",jac->maxrowsum);
    347:     if (jac->maxrowsum > 1.0) SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %g must be less than or equal one",jac->maxrowsum);
    348:   }
    349:   HYPRE_BoomerAMGSetMaxRowSum(jac->hsolver,jac->maxrowsum);

    351:   /* Grid sweeps */
    352:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for all grid levels (fine, up, and down)","None", jac->gridsweeps[0], &indx ,&flg);
    353:   if (flg) {
    354:     HYPRE_BoomerAMGSetNumSweeps(jac->hsolver,indx);
    355:     /* modify the jac structure so we can view the updated options with PC_View */
    356:     jac->gridsweeps[0] = indx;
    357:     jac->gridsweeps[1] = jac->gridsweeps[2] = jac->gridsweeps[0];
    358:     jac->gridsweeps[3] = 1;  /*The coarse level is not affected by this function - hypre code sets to 1*/
    359:   }
    360:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_fine","Number of sweeps for the fine level","None", jac->gridsweeps[0], &indx ,&flg);
    361:   if (flg) {
    362:     HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver,indx, 0);
    363:     jac->gridsweeps[0] = indx;
    364:   }
    365:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[2], &indx ,&flg);
    366:   if (flg) {
    367:     HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver,indx, 1);
    368:     jac->gridsweeps[1] = indx;
    369:   }
    370:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None", jac->gridsweeps[1], &indx ,&flg);
    371:   if (flg) {
    372:     HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver,indx, 2);
    373:     jac->gridsweeps[2] = indx;
    374:   }
    375:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None", jac->gridsweeps[3], &indx ,&flg);
    376:   if (flg) {
    377:     HYPRE_BoomerAMGSetCycleNumSweeps(jac->hsolver,indx, 3);
    378:     jac->gridsweeps[3] = indx;
    379:   }

    381:   /* Relax type */
    382:   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);
    383:   if (flg) {
    384:     jac->relaxtype[0] = jac->relaxtype[1] = jac->relaxtype[2] = indx;
    385:     jac->relaxtype[3] = 9; /* hypre code sets coarse grid to 9 (G.E.)*/
    386:     hypre_BoomerAMGSetRelaxType(jac->hsolver, indx);
    387:   }
    388:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_fine","Relax type on fine grid","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[3],&indx,&flg);
    389:   if (flg) {
    390:     jac->relaxtype[0] = indx;
    391:     HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, indx, 0);
    392:   }
    393:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[3],&indx,&flg);
    394:   if (flg) {
    395:     jac->relaxtype[1] = indx;
    396:     HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, indx, 1);
    397:   }
    398:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[3],&indx,&flg);
    399:   if (flg) {
    400:     jac->relaxtype[2] = indx;
    401:     HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, indx, 2);
    402:   }
    403:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[9],&indx,&flg);
    404:   if (flg) {
    405:     jac->relaxtype[3] = indx;
    406:     HYPRE_BoomerAMGSetCycleRelaxType(jac->hsolver, indx, 3);
    407:   }

    409:   /* Relaxation Weight */
    410:   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);
    411:   if (flg) {
    412:     hypre_BoomerAMGSetRelaxWt( jac->hsolver, tmpdbl);
    413:     jac->relaxweight = tmpdbl;
    414:   }

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

    428:   /* Outer relaxation Weight */
    429:   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);
    430:   if (flg) {
    431:     hypre_BoomerAMGSetOuterWt( jac->hsolver, tmpdbl);
    432:     jac->outerrelaxweight = tmpdbl;
    433:   }

    435:   n=2;
    436:   twodbl[0] = twodbl[1] = 1.0;
    437:   PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
    438:   if (flg) {
    439:     if (n == 2) {
    440:       indx =  (int)PetscAbsReal(twodbl[1]);
    441:       hypre_BoomerAMGSetLevelOuterWt( jac->hsolver, twodbl[0], indx);
    442:     } else {
    443:       SETERRQ1(PETSC_ERR_ARG_OUTOFRANGE,"Relax weight outer level: You must provide 2 values seperated by a comma (and no space), you provided %d",n);
    444:     }
    445:   }

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

    451:   if (flg) {
    452:     jac->relaxorder = 0;
    453:     hypre_BoomerAMGSetRelaxOrder( jac->hsolver, jac->relaxorder);
    454:   }
    455:   PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,2,HYPREBoomerAMGMeasureType[0],&indx,&flg);
    456:   if (flg) {
    457:     jac->measuretype = indx;
    458:   }
    459:   HYPRE_BoomerAMGSetMeasureType(jac->hsolver,jac->measuretype);
    460:   PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,7,HYPREBoomerAMGCoarsenType[6],&indx,&flg);
    461:   if (flg) {
    462:     jac->coarsentype = indx;
    463:   }
    464:   HYPRE_BoomerAMGSetCoarsenType(jac->hsolver,jac->coarsentype);
    465:   PetscOptionsTruth("-pc_hypre_boomeramg_print_statistics","Print statistics","None",jac->printstatistics,&jac->printstatistics,PETSC_NULL);
    466:   if (jac->printstatistics) {
    467:     HYPRE_BoomerAMGSetPrintLevel(jac->hsolver,3);
    468:     HYPRE_BoomerAMGSetDebugFlag(jac->hsolver,3);
    469:   }
    470:   PetscOptionsTail();
    471:   return(0);
    472: }

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

    482:   PetscLogInfo((pc,"PCApplyRichardson_hypre_BoomerAMG: Warning, convergence critera ignored, using %D iterations\n",its));
    483:   HYPRE_BoomerAMGSetMaxIter(jac->hsolver,its);
    484:   HYPRE_BoomerAMGSetTol(jac->hsolver,rtol);
    485:   jac->applyrichardson = PETSC_TRUE;
    486:   PCApply_HYPRE(pc,b,y);
    487:   jac->applyrichardson = PETSC_FALSE;
    488:   HYPRE_BoomerAMGSetTol(jac->hsolver,jac->tol);
    489:   HYPRE_BoomerAMGSetMaxIter(jac->hsolver,jac->maxiter);
    490:   return(0);
    491: }

    496: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
    497: {
    498:   PC_HYPRE    *jac = (PC_HYPRE*)pc->data;
    500:   PetscTruth  iascii;

    503:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    504:   if (iascii) {
    505:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG preconditioning\n");
    506:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of levels %d\n",jac->maxlevels);
    507:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of iterations %d\n",jac->maxiter);
    508:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Convergence tolerance %g\n",jac->tol);
    509:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Threshold for strong coupling %g\n",jac->strongthreshold);
    510:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum row sums %g\n",jac->maxrowsum);

    512:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps on fine grid %d\n",jac->gridsweeps[0]);
    513:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps down         %d\n",jac->gridsweeps[1]);
    514:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps up           %d\n",jac->gridsweeps[2]);
    515:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps on coarse    %d\n",jac->gridsweeps[3]);

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

    522:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax weight  (all)      %g\n",jac->relaxweight);
    523:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Outer relax weight (all) %g\n",jac->outerrelaxweight);

    525:     if (jac->relaxorder) {
    526:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Using CF-relaxation\n");
    527:     } else {
    528:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Not using CF-relaxation\n");
    529:     }
    530:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Measure type        %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
    531:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Coarsen type        %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
    532:   }
    533:   return(0);
    534: }

    536: /* --------------------------------------------------------------------------------------------*/
    539: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PC pc)
    540: {
    541:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
    543:   int            indx;
    544:   PetscTruth     flag;
    545:   const char     *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};

    548:   jac->nlevels     = 1;
    549:   jac->threshhold  = .1;
    550:   jac->filter      = .1;
    551:   jac->loadbal     = 0;
    552:   if (PetscLogPrintInfo) {
    553:     jac->logging     = (int) PETSC_TRUE;
    554:   } else {
    555:     jac->logging     = (int) PETSC_FALSE;
    556:   }
    557:   jac->ruse = (int) PETSC_FALSE;
    558:   jac->symt   = 0;

    560:   PetscOptionsHead("HYPRE ParaSails Options");
    561:   PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);
    562:   PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshhold,&jac->threshhold,0);
    563:   HYPRE_ParaSailsSetParams(jac->hsolver,jac->threshhold,jac->nlevels);

    565:   PetscOptionsReal("-pc_hypre_parasails_filter","filter","None",jac->filter,&jac->filter,0);
    566:   HYPRE_ParaSailsSetFilter(jac->hsolver,jac->filter);

    568:   PetscOptionsReal("-pc_hypre_parasails_loadbal","Load balance","None",jac->loadbal,&jac->loadbal,0);
    569:   HYPRE_ParaSailsSetLoadbal(jac->hsolver,jac->loadbal);

    571:   PetscOptionsTruth("-pc_hypre_parasails_logging","Print info to screen","None",(PetscTruth)jac->logging,(PetscTruth*)&jac->logging,0);
    572:   HYPRE_ParaSailsSetLogging(jac->hsolver,jac->logging);

    574:   PetscOptionsTruth("-pc_hypre_parasails_reuse","Reuse nonzero pattern in preconditioner","None",(PetscTruth)jac->ruse,(PetscTruth*)&jac->ruse,0);
    575:   HYPRE_ParaSailsSetReuse(jac->hsolver,jac->ruse);

    577:   PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,3,symtlist[0],&indx,&flag);
    578:   if (flag) {
    579:     jac->symt = indx;
    580:   }
    581:   HYPRE_ParaSailsSetSym(jac->hsolver,jac->symt);

    583:   PetscOptionsTail();
    584:   return(0);
    585: }

    589: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
    590: {
    591:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
    593:   PetscTruth     iascii;
    594:   const char     *symt = 0;;

    597:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&iascii);
    598:   if (iascii) {
    599:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails preconditioning\n");
    600:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: nlevels %d\n",jac->nlevels);
    601:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: threshold %g\n",jac->threshhold);
    602:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: filter %g\n",jac->filter);
    603:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: load balance %g\n",jac->loadbal);
    604:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: reuse nonzero structure %s\n",PetscTruths[jac->ruse]);
    605:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: print info to screen %s\n",PetscTruths[jac->logging]);
    606:     if (!jac->symt) {
    607:       symt = "nonsymmetric matrix and preconditioner";
    608:     } else if (jac->symt == 1) {
    609:       symt = "SPD matrix and preconditioner";
    610:     } else if (jac->symt == 2) {
    611:       symt = "nonsymmetric matrix but SPD preconditioner";
    612:     } else {
    613:       SETERRQ1(PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
    614:     }
    615:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: %s\n",symt);
    616:   }
    617:   return(0);
    618: }
    619: /* ---------------------------------------------------------------------------------*/

    624: PetscErrorCode PETSCKSP_DLLEXPORT PCHYPRESetType_HYPRE(PC pc,const char name[])
    625: {
    626:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
    628:   PetscTruth     flag;

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

    635:   pc->ops->setup          = PCSetUp_HYPRE;
    636:   pc->ops->apply          = PCApply_HYPRE;
    637:   pc->ops->destroy        = PCDestroy_HYPRE;

    639:   jac->maxiter            = PETSC_DEFAULT;
    640:   jac->tol                = PETSC_DEFAULT;
    641:   jac->printstatistics    = PetscLogPrintInfo;

    643:   PetscStrcmp("pilut",name,&flag);
    644:   if (flag) {
    645:     HYPRE_ParCSRPilutCreate(jac->comm_hypre,&jac->hsolver);
    646:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
    647:     pc->ops->view           = PCView_HYPRE_Pilut;
    648:     jac->destroy            = HYPRE_ParCSRPilutDestroy;
    649:     jac->setup              = HYPRE_ParCSRPilutSetup;
    650:     jac->solve              = HYPRE_ParCSRPilutSolve;
    651:     jac->factorrowsize      = PETSC_DEFAULT;
    652:     return(0);
    653:   }
    654:   PetscStrcmp("parasails",name,&flag);
    655:   if (flag) {
    656:     HYPRE_ParaSailsCreate(jac->comm_hypre,&jac->hsolver);
    657:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
    658:     pc->ops->view           = PCView_HYPRE_ParaSails;
    659:     jac->destroy            = HYPRE_ParaSailsDestroy;
    660:     jac->setup              = HYPRE_ParaSailsSetup;
    661:     jac->solve              = HYPRE_ParaSailsSolve;
    662:     return(0);
    663:   }
    664:   PetscStrcmp("euclid",name,&flag);
    665:   if (flag) {
    666:     HYPRE_EuclidCreate(jac->comm_hypre,&jac->hsolver);
    667:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
    668:     pc->ops->view           = PCView_HYPRE_Euclid;
    669:     jac->destroy            = HYPRE_EuclidDestroy;
    670:     jac->setup              = HYPRE_EuclidSetup;
    671:     jac->solve              = HYPRE_EuclidSolve;
    672:     return(0);
    673:   }
    674:   PetscStrcmp("boomeramg",name,&flag);
    675:   if (flag) {
    676:     HYPRE_BoomerAMGCreate(&jac->hsolver);
    677:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_BoomerAMG;
    678:     pc->ops->view            = PCView_HYPRE_BoomerAMG;
    679:     jac->destroy             = HYPRE_BoomerAMGDestroy;
    680:     jac->setup               = HYPRE_BoomerAMGSetup;
    681:     jac->solve               = HYPRE_BoomerAMGSolve;
    682:     pc->ops->applyrichardson = PCApplyRichardson_BoomerAMG;
    683:     return(0);
    684:   }
    685:   SETERRQ1(PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are pilut, parasails, euclid, boomeramg",name);
    686:   return(0);
    687: }

    690: /*
    691:     It only gets here if the HYPRE type has not been set before the call to
    692:    ...SetFromOptions() which actually is most of the time
    693: */
    696: static PetscErrorCode PCSetFromOptions_HYPRE(PC pc)
    697: {
    699:   int            indx;
    700:   const char     *type[] = {"pilut","parasails","boomeramg","euclid"};
    701:   PetscTruth     flg;

    704:   PetscOptionsHead("HYPRE preconditioner options");
    705:   PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"pilut",&indx,&flg);
    706:   if (PetscOptionsPublishCount) {   /* force the default if it was not yet set and user did not set with option */
    707:     if (!flg && !pc->ops->apply) {
    708:       flg   = PETSC_TRUE;
    709:       indx = 0;
    710:     }
    711:   }
    712:   if (flg) {
    713:     PCHYPRESetType_HYPRE(pc,type[indx]);
    714:   }
    715:   if (pc->ops->setfromoptions) {
    716:     pc->ops->setfromoptions(pc);
    717:   }
    718:   PetscOptionsTail();
    719:   return(0);
    720: }

    724: /*@C
    725:      PCHYPRESetType - Sets which hypre preconditioner you wish to use

    727:    Input Parameters:
    728: +     pc - the preconditioner context
    729: -     name - either  pilut, parasails, boomeramg, euclid

    731:    Options Database Keys:
    732:    -pc_hypre_type - One of pilut, parasails, boomeramg, euclid
    733:
    734:    Level: intermediate

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

    739: @*/
    740: PetscErrorCode PETSCKSP_DLLEXPORT PCHYPRESetType(PC pc,const char name[])
    741: {
    742:   PetscErrorCode ierr,(*f)(PC,const char[]);

    747:   PetscObjectQueryFunction((PetscObject)pc,"PCHYPRESetType_C",(void (**)(void))&f);
    748:   if (f) {
    749:     (*f)(pc,name);
    750:   }
    751:   return(0);
    752: }

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

    757:    Options Database Keys:
    758: +   -pc_hypre_type - One of pilut, parasails, boomerAMG, euclid
    759: -   Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX
    760:           preconditioner
    761:
    762:    Level: intermediate

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

    768:           If you wish to use boomerAMG WITHOUT a Krylov method use -ksp_type richardson NOT -ksp_type preonly
    769:           (see the PETSc FAQ.html at the PETSc website under the Documentation tab).

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

    774: M*/

    779: PetscErrorCode PETSCKSP_DLLEXPORT PCCreate_HYPRE(PC pc)
    780: {
    781:   PC_HYPRE       *jac;

    785:   PetscNew(PC_HYPRE,&jac);
    786:   pc->data                 = jac;
    787:   pc->ops->setfromoptions  = PCSetFromOptions_HYPRE;
    788:   /* Com_dup for hypre */
    789:   MPI_Comm_dup(pc->comm,&(jac->comm_hypre));
    790:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCHYPRESetType_C","PCHYPRESetType_HYPRE",PCHYPRESetType_HYPRE);
    791:   return(0);
    792: }