Actual source code: hypre.c

petsc-dev 2014-08-19
Report Typos and Errors
  2: /*
  3:    Provides an interface to the LLNL package hypre
  4: */

  6: /* Must use hypre 2.0.0 or more recent. */

  8: #include <petsc-private/pcimpl.h>          /*I "petscpc.h" I*/
  9: #include <../src/dm/impls/da/hypre/mhyp.h>

 11: static PetscBool cite = PETSC_FALSE;
 12: static const char hypreCitation[] = "@manual{hypre-web-page,\n  title  = {{\\sl hypre}: High Performance Preconditioners},\n  organization = {Lawrence Livermore National Laboratory},\n  note  = {\\url{http://www.llnl.gov/CASC/hypre/}}\n}\n";

 14: /*
 15:    Private context (data structure) for the  preconditioner.
 16: */
 17: typedef struct {
 18:   HYPRE_Solver   hsolver;
 19:   HYPRE_IJMatrix ij;
 20:   HYPRE_IJVector b,x;

 22:   HYPRE_Int (*destroy)(HYPRE_Solver);
 23:   HYPRE_Int (*solve)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
 24:   HYPRE_Int (*setup)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);

 26:   MPI_Comm comm_hypre;
 27:   char     *hypre_type;

 29:   /* options for Pilut and BoomerAMG*/
 30:   PetscInt maxiter;
 31:   double   tol;

 33:   /* options for Pilut */
 34:   PetscInt factorrowsize;

 36:   /* options for ParaSails */
 37:   PetscInt nlevels;
 38:   double   threshhold;
 39:   double   filter;
 40:   PetscInt sym;
 41:   double   loadbal;
 42:   PetscInt logging;
 43:   PetscInt ruse;
 44:   PetscInt symt;

 46:   /* options for Euclid */
 47:   PetscBool bjilu;
 48:   PetscInt  levels;

 50:   /* options for Euclid and BoomerAMG */
 51:   PetscBool printstatistics;

 53:   /* options for BoomerAMG */
 54:   PetscInt  cycletype;
 55:   PetscInt  maxlevels;
 56:   double    strongthreshold;
 57:   double    maxrowsum;
 58:   PetscInt  gridsweeps[3];
 59:   PetscInt  coarsentype;
 60:   PetscInt  measuretype;
 61:   PetscInt  relaxtype[3];
 62:   double    relaxweight;
 63:   double    outerrelaxweight;
 64:   PetscInt  relaxorder;
 65:   double    truncfactor;
 66:   PetscBool applyrichardson;
 67:   PetscInt  pmax;
 68:   PetscInt  interptype;
 69:   PetscInt  agg_nl;
 70:   PetscInt  agg_num_paths;
 71:   PetscInt  nodal_coarsen;
 72:   PetscBool nodal_relax;
 73:   PetscInt  nodal_relax_levels;
 74: } PC_HYPRE;

 78: PetscErrorCode PCHYPREGetSolver(PC pc,HYPRE_Solver *hsolver)
 79: {
 80:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

 83:   *hsolver = jac->hsolver;
 84:   return(0);
 85: }

 89: static PetscErrorCode PCSetUp_HYPRE(PC pc)
 90: {
 91:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
 92:   PetscErrorCode     ierr;
 93:   HYPRE_ParCSRMatrix hmat;
 94:   HYPRE_ParVector    bv,xv;
 95:   PetscInt           bs;

 98:   if (!jac->hypre_type) {
 99:     PCHYPRESetType(pc,"boomeramg");
100:   }

102:   if (pc->setupcalled) {
103:     /* always destroy the old matrix and create a new memory;
104:        hope this does not churn the memory too much. The problem
105:        is I do not know if it is possible to put the matrix back to
106:        its initial state so that we can directly copy the values
107:        the second time through. */
108:     PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij));
109:     jac->ij = 0;
110:   }

112:   if (!jac->ij) { /* create the matrix the first time through */
113:     MatHYPRE_IJMatrixCreate(pc->pmat,&jac->ij);
114:   }
115:   if (!jac->b) { /* create the vectors the first time through */
116:     Vec x,b;
117:     MatGetVecs(pc->pmat,&x,&b);
118:     VecHYPRE_IJVectorCreate(x,&jac->x);
119:     VecHYPRE_IJVectorCreate(b,&jac->b);
120:     VecDestroy(&x);
121:     VecDestroy(&b);
122:   }

124:   /* special case for BoomerAMG */
125:   if (jac->setup == HYPRE_BoomerAMGSetup) {
126:     MatGetBlockSize(pc->pmat,&bs);
127:     if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
128:   };
129:   MatHYPRE_IJMatrixCopy(pc->pmat,jac->ij);
130:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
131:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&bv));
132:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&xv));
133:   PetscStackCall("HYPRE_SetupXXX",(*jac->setup)(jac->hsolver,hmat,bv,xv););
134:   return(0);
135: }

137: /*
138:     Replaces the address where the HYPRE vector points to its data with the address of
139:   PETSc's data. Saves the old address so it can be reset when we are finished with it.
140:   Allows use to get the data into a HYPRE vector without the cost of memcopies
141: */
142: #define HYPREReplacePointer(b,newvalue,savedvalue) { \
143:     hypre_ParVector *par_vector   = (hypre_ParVector*)hypre_IJVectorObject(((hypre_IJVector*)b)); \
144:     hypre_Vector    *local_vector = hypre_ParVectorLocalVector(par_vector); \
145:     savedvalue         = local_vector->data; \
146:     local_vector->data = newvalue;          \
147: }

151: static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
152: {
153:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
154:   PetscErrorCode     ierr;
155:   HYPRE_ParCSRMatrix hmat;
156:   PetscScalar        *bv,*xv;
157:   HYPRE_ParVector    jbv,jxv;
158:   PetscScalar        *sbv,*sxv;
159:   PetscInt           hierr;

162:   PetscCitationsRegister(hypreCitation,&cite);
163:   if (!jac->applyrichardson) {VecSet(x,0.0);}
164:   VecGetArray(b,&bv);
165:   VecGetArray(x,&xv);
166:   HYPREReplacePointer(jac->b,bv,sbv);
167:   HYPREReplacePointer(jac->x,xv,sxv);

169:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
170:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
171:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));
172:   PetscStackCall("Hypre solve",h(*jac->solve)(jac->hsolver,hmat,jbv,jxv);
173:                                if (hierr && hierr != HYPRE_ERROR_CONV) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
174:                                if (hierr) hypre__global_error = 0;);

176:   HYPREReplacePointer(jac->b,sbv,bv);
177:   HYPREReplacePointer(jac->x,sxv,xv);
178:   VecRestoreArray(x,&xv);
179:   VecRestoreArray(b,&bv);
180:   return(0);
181: }

185: static PetscErrorCode PCDestroy_HYPRE(PC pc)
186: {
187:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;

191:   if (jac->ij) PetscStackCallStandard(HYPRE_IJMatrixDestroy,(jac->ij));
192:   if (jac->b) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->b));
193:   if (jac->x) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->x));
194:   if (jac->destroy) PetscStackCall("HYPRE_DestroyXXX",(*jac->destroy)(jac->hsolver););
195:   PetscFree(jac->hypre_type);
196:   if (jac->comm_hypre != MPI_COMM_NULL) { MPI_Comm_free(&(jac->comm_hypre));}
197:   PetscFree(pc->data);

199:   PetscObjectChangeTypeName((PetscObject)pc,0);
200:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);
201:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);
202:   return(0);
203: }

205: /* --------------------------------------------------------------------------------------------*/
208: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PC pc)
209: {
210:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
212:   PetscBool      flag;

215:   PetscOptionsHead("HYPRE Pilut Options");
216:   PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
217:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter));
218:   PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
219:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol));
220:   PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
221:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize));
222:   PetscOptionsTail();
223:   return(0);
224: }

228: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
229: {
230:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
232:   PetscBool      iascii;

235:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
236:   if (iascii) {
237:     PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut preconditioning\n");
238:     if (jac->maxiter != PETSC_DEFAULT) {
239:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: maximum number of iterations %d\n",jac->maxiter);
240:     } else {
241:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default maximum number of iterations \n");
242:     }
243:     if (jac->tol != PETSC_DEFAULT) {
244:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: drop tolerance %g\n",(double)jac->tol);
245:     } else {
246:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default drop tolerance \n");
247:     }
248:     if (jac->factorrowsize != PETSC_DEFAULT) {
249:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: factor row size %d\n",jac->factorrowsize);
250:     } else {
251:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default factor row size \n");
252:     }
253:   }
254:   return(0);
255: }

257: /* --------------------------------------------------------------------------------------------*/
260: static PetscErrorCode PCSetFromOptions_HYPRE_Euclid(PC pc)
261: {
262:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
264:   PetscBool      flag;
265:   char           *args[8],levels[16];
266:   PetscInt       cnt = 0;

269:   PetscOptionsHead("HYPRE Euclid Options");
270:   PetscOptionsInt("-pc_hypre_euclid_levels","Number of levels of fill ILU(k)","None",jac->levels,&jac->levels,&flag);
271:   if (flag) {
272:     if (jac->levels < 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be nonegative",jac->levels);
273:     PetscSNPrintf(levels,sizeof(levels),"%D",jac->levels);
274:     args[cnt++] = (char*)"-level"; args[cnt++] = levels;
275:   }
276:   PetscOptionsBool("-pc_hypre_euclid_bj","Use block Jacobi ILU(k)","None",jac->bjilu,&jac->bjilu,NULL);
277:   if (jac->bjilu) {
278:     args[cnt++] =(char*) "-bj"; args[cnt++] = (char*)"1";
279:   }

281:   PetscOptionsBool("-pc_hypre_euclid_print_statistics","Print statistics","None",jac->printstatistics,&jac->printstatistics,NULL);
282:   if (jac->printstatistics) {
283:     args[cnt++] = (char*)"-eu_stats"; args[cnt++] = (char*)"1";
284:     args[cnt++] = (char*)"-eu_mem"; args[cnt++] = (char*)"1";
285:   }
286:   PetscOptionsTail();
287:   if (cnt) PetscStackCallStandard(HYPRE_EuclidSetParams,(jac->hsolver,cnt,args));
288:   return(0);
289: }

293: static PetscErrorCode PCView_HYPRE_Euclid(PC pc,PetscViewer viewer)
294: {
295:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
297:   PetscBool      iascii;

300:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
301:   if (iascii) {
302:     PetscViewerASCIIPrintf(viewer,"  HYPRE Euclid preconditioning\n");
303:     PetscViewerASCIIPrintf(viewer,"  HYPRE Euclid: number of levels %d\n",jac->levels);
304:     if (jac->bjilu) {
305:       PetscViewerASCIIPrintf(viewer,"  HYPRE Euclid: Using block Jacobi ILU instead of parallel ILU\n");
306:     }
307:   }
308:   return(0);
309: }

311: /* --------------------------------------------------------------------------------------------*/

315: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
316: {
317:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
318:   PetscErrorCode     ierr;
319:   HYPRE_ParCSRMatrix hmat;
320:   PetscScalar        *bv,*xv;
321:   HYPRE_ParVector    jbv,jxv;
322:   PetscScalar        *sbv,*sxv;
323:   PetscInt           hierr;

326:   PetscCitationsRegister(hypreCitation,&cite);
327:   VecSet(x,0.0);
328:   VecGetArray(b,&bv);
329:   VecGetArray(x,&xv);
330:   HYPREReplacePointer(jac->b,bv,sbv);
331:   HYPREReplacePointer(jac->x,xv,sxv);

333:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(jac->ij,(void**)&hmat));
334:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->b,(void**)&jbv));
335:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->x,(void**)&jxv));

337:   hHYPRE_BoomerAMGSolveT(jac->hsolver,hmat,jbv,jxv);
338:   /* error code of 1 in BoomerAMG merely means convergence not achieved */
339:   if (hierr && (hierr != 1)) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
340:   if (hierr) hypre__global_error = 0;

342:   HYPREReplacePointer(jac->b,sbv,bv);
343:   HYPREReplacePointer(jac->x,sxv,xv);
344:   VecRestoreArray(x,&xv);
345:   VecRestoreArray(b,&bv);
346:   return(0);
347: }

349: /* static array length */
350: #define ALEN(a) (sizeof(a)/sizeof((a)[0]))

352: static const char *HYPREBoomerAMGCycleType[]   = {"","V","W"};
353: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
354: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
355: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
356: static const char *HYPREBoomerAMGRelaxType[]   = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi",
357:                                                   "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi",
358:                                                   "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination",
359:                                                   "" /* 10 */, "" /* 11 */, "" /* 12 */, "" /* 13 */, "" /* 14 */,
360:                                                   "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"};
361: static const char *HYPREBoomerAMGInterpType[]  = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
362:                                                   "ext+i-cc", "standard", "standard-wts", "", "", "FF", "FF1"};
365: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PC pc)
366: {
367:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
369:   PetscInt       n,indx,level;
370:   PetscBool      flg, tmp_truth;
371:   double         tmpdbl, twodbl[2];

374:   PetscOptionsHead("HYPRE BoomerAMG Options");
375:   PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);
376:   if (flg) {
377:     jac->cycletype = indx+1;
378:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
379:   }
380:   PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
381:   if (flg) {
382:     if (jac->maxlevels < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
383:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
384:   }
385:   PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
386:   if (flg) {
387:     if (jac->maxiter < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
388:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
389:   }
390:   PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);
391:   if (flg) {
392:     if (jac->tol < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Tolerance %g must be greater than or equal to zero",(double)jac->tol);
393:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
394:   }

396:   PetscOptionsScalar("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);
397:   if (flg) {
398:     if (jac->truncfactor < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Truncation factor %g must be great than or equal zero",(double)jac->truncfactor);
399:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
400:   }

402:   PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);
403:   if (flg) {
404:     if (jac->pmax < 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"P_max %g must be greater than or equal to zero",(double)jac->pmax);
405:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
406:   }

408:   PetscOptionsInt("-pc_hypre_boomeramg_agg_nl","Number of levels of aggressive coarsening","None",jac->agg_nl,&jac->agg_nl,&flg);
409:   if (flg) {
410:     if (jac->agg_nl < 0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %g must be greater than or equal to zero",(double)jac->agg_nl);

412:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
413:   }


416:   PetscOptionsInt("-pc_hypre_boomeramg_agg_num_paths","Number of paths for aggressive coarsening","None",jac->agg_num_paths,&jac->agg_num_paths,&flg);
417:   if (flg) {
418:     if (jac->agg_num_paths < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of paths %g must be greater than or equal to 1",(double)jac->agg_num_paths);

420:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
421:   }


424:   PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
425:   if (flg) {
426:     if (jac->strongthreshold < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Strong threshold %g must be great than or equal zero",(double)jac->strongthreshold);
427:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
428:   }
429:   PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
430:   if (flg) {
431:     if (jac->maxrowsum < 0.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %g must be greater than zero",(double)jac->maxrowsum);
432:     if (jac->maxrowsum > 1.0) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Maximum row sum %g must be less than or equal one",(double)jac->maxrowsum);
433:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
434:   }

436:   /* Grid sweeps */
437:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);
438:   if (flg) {
439:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx));
440:     /* modify the jac structure so we can view the updated options with PC_View */
441:     jac->gridsweeps[0] = indx;
442:     jac->gridsweeps[1] = indx;
443:     /*defaults coarse to 1 */
444:     jac->gridsweeps[2] = 1;
445:   }

447:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);
448:   if (flg) {
449:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1));
450:     jac->gridsweeps[0] = indx;
451:   }
452:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);
453:   if (flg) {
454:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2));
455:     jac->gridsweeps[1] = indx;
456:   }
457:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);
458:   if (flg) {
459:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3));
460:     jac->gridsweeps[2] = indx;
461:   }

463:   /* Relax type */
464:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
465:   if (flg) {
466:     jac->relaxtype[0] = jac->relaxtype[1]  = indx;
467:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx));
468:     /* by default, coarse type set to 9 */
469:     jac->relaxtype[2] = 9;

471:   }
472:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
473:   if (flg) {
474:     jac->relaxtype[0] = indx;
475:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1));
476:   }
477:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
478:   if (flg) {
479:     jac->relaxtype[1] = indx;
480:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2));
481:   }
482:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);
483:   if (flg) {
484:     jac->relaxtype[2] = indx;
485:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3));
486:   }

488:   /* Relaxation Weight */
489:   PetscOptionsReal("-pc_hypre_boomeramg_relax_weight_all","Relaxation weight for all levels (0 = hypre estimates, -k = determined with k CG steps)","None",jac->relaxweight, &tmpdbl,&flg);
490:   if (flg) {
491:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl));
492:     jac->relaxweight = tmpdbl;
493:   }

495:   n         = 2;
496:   twodbl[0] = twodbl[1] = 1.0;
497:   PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
498:   if (flg) {
499:     if (n == 2) {
500:       indx =  (int)PetscAbsReal(twodbl[1]);
501:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx));
502:     } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relax weight level: you must provide 2 values separated by a comma (and no space), you provided %d",n);
503:   }

505:   /* Outer relaxation Weight */
506:   PetscOptionsReal("-pc_hypre_boomeramg_outer_relax_weight_all","Outer relaxation weight for all levels (-k = determined with k CG steps)","None",jac->outerrelaxweight, &tmpdbl,&flg);
507:   if (flg) {
508:     PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl));
509:     jac->outerrelaxweight = tmpdbl;
510:   }

512:   n         = 2;
513:   twodbl[0] = twodbl[1] = 1.0;
514:   PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
515:   if (flg) {
516:     if (n == 2) {
517:       indx =  (int)PetscAbsReal(twodbl[1]);
518:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx));
519:     } else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Relax weight outer level: You must provide 2 values separated by a comma (and no space), you provided %d",n);
520:   }

522:   /* the Relax Order */
523:   PetscOptionsBool("-pc_hypre_boomeramg_no_CF", "Do not use CF-relaxation", "None", PETSC_FALSE, &tmp_truth, &flg);

525:   if (flg) {
526:     jac->relaxorder = 0;
527:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
528:   }
529:   PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);
530:   if (flg) {
531:     jac->measuretype = indx;
532:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
533:   }
534:   /* update list length 3/07 */
535:   PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);
536:   if (flg) {
537:     jac->coarsentype = indx;
538:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
539:   }

541:   /* new 3/07 */
542:   PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);
543:   if (flg) {
544:     jac->interptype = indx;
545:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
546:   }

548:   PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
549:   if (flg) {
550:     level = 3;
551:     PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);

553:     jac->printstatistics = PETSC_TRUE;
554:     PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level));
555:   }

557:   PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);
558:   if (flg) {
559:     level = 3;
560:     PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);

562:     jac->printstatistics = PETSC_TRUE;
563:     PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level));
564:   }

566:   PetscOptionsBool("-pc_hypre_boomeramg_nodal_coarsen", "HYPRE_BoomerAMGSetNodal()", "None", PETSC_FALSE, &tmp_truth, &flg);
567:   if (flg && tmp_truth) {
568:     jac->nodal_coarsen = 1;
569:     PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,1));
570:   }

572:   PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);
573:   if (flg && tmp_truth) {
574:     PetscInt tmp_int;
575:     PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);
576:     if (flg) jac->nodal_relax_levels = tmp_int;
577:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6));
578:     PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1));
579:     PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0));
580:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels));
581:   }

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

589: static PetscErrorCode PCApplyRichardson_HYPRE_BoomerAMG(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
590: {
591:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
593:   PetscInt       oits;

596:   PetscCitationsRegister(hypreCitation,&cite);
597:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter));
598:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol));
599:   jac->applyrichardson = PETSC_TRUE;
600:   PCApply_HYPRE(pc,b,y);
601:   jac->applyrichardson = PETSC_FALSE;
602:   PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
603:   *outits = oits;
604:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
605:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
606:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
607:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
608:   return(0);
609: }


614: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
615: {
616:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
618:   PetscBool      iascii;

621:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
622:   if (iascii) {
623:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG preconditioning\n");
624:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);
625:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of levels %d\n",jac->maxlevels);
626:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of iterations PER hypre call %d\n",jac->maxiter);
627:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Convergence tolerance PER hypre call %g\n",(double)jac->tol);
628:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Threshold for strong coupling %g\n",(double)jac->strongthreshold);
629:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation truncation factor %g\n",(double)jac->truncfactor);
630:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation: max elements per row %d\n",jac->pmax);
631:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Number of levels of aggressive coarsening %d\n",jac->agg_nl);
632:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Number of paths for aggressive coarsening %d\n",jac->agg_num_paths);

634:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum row sums %g\n",(double)jac->maxrowsum);

636:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps down         %d\n",jac->gridsweeps[0]);
637:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps up           %d\n",jac->gridsweeps[1]);
638:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps on coarse    %d\n",jac->gridsweeps[2]);

640:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax down          %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[0]]);
641:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax up            %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[1]]);
642:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax on coarse     %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[2]]);

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

647:     if (jac->relaxorder) {
648:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Using CF-relaxation\n");
649:     } else {
650:       PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Not using CF-relaxation\n");
651:     }
652:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Measure type        %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
653:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Coarsen type        %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
654:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Interpolation type  %s\n",HYPREBoomerAMGInterpType[jac->interptype]);
655:     if (jac->nodal_coarsen) {
656:       PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal coarsening (with HYPRE_BOOMERAMGSetNodal())\n");
657:     }
658:     if (jac->nodal_relax) {
659:       PetscViewerASCIIPrintf(viewer," HYPRE BoomerAMG: Using nodal relaxation via Schwarz smoothing on levels %d\n",jac->nodal_relax_levels);
660:     }
661:   }
662:   return(0);
663: }

665: /* --------------------------------------------------------------------------------------------*/
668: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PC pc)
669: {
670:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
672:   PetscInt       indx;
673:   PetscBool      flag;
674:   const char     *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};

677:   PetscOptionsHead("HYPRE ParaSails Options");
678:   PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);
679:   PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshhold,&jac->threshhold,&flag);
680:   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));

682:   PetscOptionsReal("-pc_hypre_parasails_filter","filter","None",jac->filter,&jac->filter,&flag);
683:   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));

685:   PetscOptionsReal("-pc_hypre_parasails_loadbal","Load balance","None",jac->loadbal,&jac->loadbal,&flag);
686:   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));

688:   PetscOptionsBool("-pc_hypre_parasails_logging","Print info to screen","None",(PetscBool)jac->logging,(PetscBool*)&jac->logging,&flag);
689:   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));

691:   PetscOptionsBool("-pc_hypre_parasails_reuse","Reuse nonzero pattern in preconditioner","None",(PetscBool)jac->ruse,(PetscBool*)&jac->ruse,&flag);
692:   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));

694:   PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);
695:   if (flag) {
696:     jac->symt = indx;
697:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
698:   }

700:   PetscOptionsTail();
701:   return(0);
702: }

706: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
707: {
708:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
710:   PetscBool      iascii;
711:   const char     *symt = 0;;

714:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
715:   if (iascii) {
716:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails preconditioning\n");
717:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: nlevels %d\n",jac->nlevels);
718:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: threshold %g\n",(double)jac->threshhold);
719:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: filter %g\n",(double)jac->filter);
720:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: load balance %g\n",(double)jac->loadbal);
721:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: reuse nonzero structure %s\n",PetscBools[jac->ruse]);
722:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: print info to screen %s\n",PetscBools[jac->logging]);
723:     if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
724:     else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
725:     else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
726:     else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
727:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: %s\n",symt);
728:   }
729:   return(0);
730: }
731: /* ---------------------------------------------------------------------------------*/

735: static PetscErrorCode  PCHYPREGetType_HYPRE(PC pc,const char *name[])
736: {
737:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

740:   *name = jac->hypre_type;
741:   return(0);
742: }

746: static PetscErrorCode  PCHYPRESetType_HYPRE(PC pc,const char name[])
747: {
748:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
750:   PetscBool      flag;

753:   if (jac->hypre_type) {
754:     PetscStrcmp(jac->hypre_type,name,&flag);
755:     if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set");
756:     return(0);
757:   } else {
758:     PetscStrallocpy(name, &jac->hypre_type);
759:   }

761:   jac->maxiter         = PETSC_DEFAULT;
762:   jac->tol             = PETSC_DEFAULT;
763:   jac->printstatistics = PetscLogPrintInfo;

765:   PetscStrcmp("pilut",jac->hypre_type,&flag);
766:   if (flag) {
767:     PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver));
768:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
769:     pc->ops->view           = PCView_HYPRE_Pilut;
770:     jac->destroy            = HYPRE_ParCSRPilutDestroy;
771:     jac->setup              = HYPRE_ParCSRPilutSetup;
772:     jac->solve              = HYPRE_ParCSRPilutSolve;
773:     jac->factorrowsize      = PETSC_DEFAULT;
774:     return(0);
775:   }
776:   PetscStrcmp("parasails",jac->hypre_type,&flag);
777:   if (flag) {
778:     PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver));
779:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
780:     pc->ops->view           = PCView_HYPRE_ParaSails;
781:     jac->destroy            = HYPRE_ParaSailsDestroy;
782:     jac->setup              = HYPRE_ParaSailsSetup;
783:     jac->solve              = HYPRE_ParaSailsSolve;
784:     /* initialize */
785:     jac->nlevels    = 1;
786:     jac->threshhold = .1;
787:     jac->filter     = .1;
788:     jac->loadbal    = 0;
789:     if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE;
790:     else jac->logging = (int) PETSC_FALSE;

792:     jac->ruse = (int) PETSC_FALSE;
793:     jac->symt = 0;
794:     PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
795:     PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
796:     PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
797:     PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
798:     PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
799:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
800:     return(0);
801:   }
802:   PetscStrcmp("euclid",jac->hypre_type,&flag);
803:   if (flag) {
804:     HYPRE_EuclidCreate(jac->comm_hypre,&jac->hsolver);
805:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
806:     pc->ops->view           = PCView_HYPRE_Euclid;
807:     jac->destroy            = HYPRE_EuclidDestroy;
808:     jac->setup              = HYPRE_EuclidSetup;
809:     jac->solve              = HYPRE_EuclidSolve;
810:     /* initialization */
811:     jac->bjilu              = PETSC_FALSE;
812:     jac->levels             = 1;
813:     return(0);
814:   }
815:   PetscStrcmp("boomeramg",jac->hypre_type,&flag);
816:   if (flag) {
817:     HYPRE_BoomerAMGCreate(&jac->hsolver);
818:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_BoomerAMG;
819:     pc->ops->view            = PCView_HYPRE_BoomerAMG;
820:     pc->ops->applytranspose  = PCApplyTranspose_HYPRE_BoomerAMG;
821:     pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
822:     jac->destroy             = HYPRE_BoomerAMGDestroy;
823:     jac->setup               = HYPRE_BoomerAMGSetup;
824:     jac->solve               = HYPRE_BoomerAMGSolve;
825:     jac->applyrichardson     = PETSC_FALSE;
826:     /* these defaults match the hypre defaults */
827:     jac->cycletype        = 1;
828:     jac->maxlevels        = 25;
829:     jac->maxiter          = 1;
830:     jac->tol              = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
831:     jac->truncfactor      = 0.0;
832:     jac->strongthreshold  = .25;
833:     jac->maxrowsum        = .9;
834:     jac->coarsentype      = 6;
835:     jac->measuretype      = 0;
836:     jac->gridsweeps[0]    = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
837:     jac->relaxtype[0]     = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */
838:     jac->relaxtype[2]     = 9; /*G.E. */
839:     jac->relaxweight      = 1.0;
840:     jac->outerrelaxweight = 1.0;
841:     jac->relaxorder       = 1;
842:     jac->interptype       = 0;
843:     jac->agg_nl           = 0;
844:     jac->pmax             = 0;
845:     jac->truncfactor      = 0.0;
846:     jac->agg_num_paths    = 1;

848:     jac->nodal_coarsen      = 0;
849:     jac->nodal_relax        = PETSC_FALSE;
850:     jac->nodal_relax_levels = 1;
851:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
852:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
853:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
854:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
855:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
856:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
857:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
858:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
859:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
860:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
861:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
862:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
863:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
864:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
865:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0]));  /*defaults coarse to 9*/
866:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */
867:     return(0);
868:   }
869:   PetscFree(jac->hypre_type);

871:   jac->hypre_type = NULL;
872:   SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are pilut, parasails, euclid, boomeramg",name);
873:   return(0);
874: }

876: /*
877:     It only gets here if the HYPRE type has not been set before the call to
878:    ...SetFromOptions() which actually is most of the time
879: */
882: static PetscErrorCode PCSetFromOptions_HYPRE(PC pc)
883: {
885:   PetscInt       indx;
886:   const char     *type[] = {"pilut","parasails","boomeramg","euclid"};
887:   PetscBool      flg;

890:   PetscOptionsHead("HYPRE preconditioner options");
891:   PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"boomeramg",&indx,&flg);
892:   if (flg) {
893:     PCHYPRESetType_HYPRE(pc,type[indx]);
894:   } else {
895:     PCHYPRESetType_HYPRE(pc,"boomeramg");
896:   }
897:   if (pc->ops->setfromoptions) {
898:     pc->ops->setfromoptions(pc);
899:   }
900:   PetscOptionsTail();
901:   return(0);
902: }

906: /*@C
907:      PCHYPRESetType - Sets which hypre preconditioner you wish to use

909:    Input Parameters:
910: +     pc - the preconditioner context
911: -     name - either  pilut, parasails, boomeramg, euclid

913:    Options Database Keys:
914:    -pc_hypre_type - One of pilut, parasails, boomeramg, euclid

916:    Level: intermediate

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

921: @*/
922: PetscErrorCode  PCHYPRESetType(PC pc,const char name[])
923: {

929:   PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));
930:   return(0);
931: }

935: /*@C
936:      PCHYPREGetType - Gets which hypre preconditioner you are using

938:    Input Parameter:
939: .     pc - the preconditioner context

941:    Output Parameter:
942: .     name - either  pilut, parasails, boomeramg, euclid

944:    Level: intermediate

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

949: @*/
950: PetscErrorCode  PCHYPREGetType(PC pc,const char *name[])
951: {

957:   PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));
958:   return(0);
959: }

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

964:    Options Database Keys:
965: +   -pc_hypre_type - One of pilut, parasails, boomeramg, euclid
966: -   Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX
967:           preconditioner

969:    Level: intermediate

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

975:           The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_rtol refer to the number of iterations
976:           (V-cycles) and tolerance that boomeramg does EACH time it is called. So for example, if
977:           -pc_hypre_boomeramg_max_iter is set to 2 then 2-V-cycles are being used to define the preconditioner
978:           (-pc_hypre_boomeramg_rtol should be set to 0.0 - the default - to strictly use a fixed number of
979:           iterations per hypre call). -ksp_max_it and -ksp_rtol STILL determine the total number of iterations
980:           and tolerance for the Krylov solver. For example, if -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10
981:           then AT MOST twenty V-cycles of boomeramg will be called.

983:            Note that the option -pc_hypre_boomeramg_relax_type_all defaults to symmetric relaxation
984:            (symmetric-SOR/Jacobi), which is required for Krylov solvers like CG that expect symmetry.
985:            Otherwise, you may want to use -pc_hypre_boomeramg_relax_type_all SOR/Jacobi.
986:           If you wish to use BoomerAMG WITHOUT a Krylov method use -ksp_type richardson NOT -ksp_type preonly
987:           and use -ksp_max_it to control the number of V-cycles.
988:           (see the PETSc FAQ.html at the PETSc website under the Documentation tab).

990:           2007-02-03 Using HYPRE-1.11.1b, the routine HYPRE_BoomerAMGSolveT and the option
991:           -pc_hypre_parasails_reuse were failing with SIGSEGV. Dalcin L.

993:           See PCPFMG for access to the hypre Struct PFMG solver

995: .seealso:  PCCreate(), PCSetType(), PCType (for list of available types), PC,
996:            PCHYPRESetType(), PCPFMG

998: M*/

1002: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
1003: {
1004:   PC_HYPRE       *jac;

1008:   PetscNewLog(pc,&jac);

1010:   pc->data                = jac;
1011:   pc->ops->destroy        = PCDestroy_HYPRE;
1012:   pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
1013:   pc->ops->setup          = PCSetUp_HYPRE;
1014:   pc->ops->apply          = PCApply_HYPRE;
1015:   jac->comm_hypre         = MPI_COMM_NULL;
1016:   jac->hypre_type         = NULL;
1017:   /* duplicate communicator for hypre */
1018:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1019:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);
1020:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);
1021:   return(0);
1022: }

1024: /* ---------------------------------------------------------------------------------------------------------------------------------*/

1026: /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */
1027: #include <petsc-private/matimpl.h>

1029: typedef struct {
1030:   MPI_Comm           hcomm;        /* does not share comm with HYPRE_StructMatrix because need to create solver before getting matrix */
1031:   HYPRE_StructSolver hsolver;

1033:   /* keep copy of PFMG options used so may view them */
1034:   PetscInt its;
1035:   double   tol;
1036:   PetscInt relax_type;
1037:   PetscInt rap_type;
1038:   PetscInt num_pre_relax,num_post_relax;
1039:   PetscInt max_levels;
1040: } PC_PFMG;

1044: PetscErrorCode PCDestroy_PFMG(PC pc)
1045: {
1047:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

1050:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1051:   MPI_Comm_free(&ex->hcomm);
1052:   PetscFree(pc->data);
1053:   return(0);
1054: }

1056: static const char *PFMGRelaxType[] = {"Jacobi","Weighted-Jacobi","symmetric-Red/Black-Gauss-Seidel","Red/Black-Gauss-Seidel"};
1057: static const char *PFMGRAPType[] = {"Galerkin","non-Galerkin"};

1061: PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
1062: {
1064:   PetscBool      iascii;
1065:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

1068:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1069:   if (iascii) {
1070:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG preconditioning\n");
1071:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: max iterations %d\n",ex->its);
1072:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: tolerance %g\n",ex->tol);
1073:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
1074:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: RAP type %s\n",PFMGRAPType[ex->rap_type]);
1075:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
1076:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG: max levels %d\n",ex->max_levels);
1077:   }
1078:   return(0);
1079: }


1084: PetscErrorCode PCSetFromOptions_PFMG(PC pc)
1085: {
1087:   PC_PFMG        *ex = (PC_PFMG*) pc->data;
1088:   PetscBool      flg = PETSC_FALSE;

1091:   PetscOptionsHead("PFMG options");
1092:   PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);
1093:   if (flg) {
1094:     int level=3;
1095:     PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level));
1096:   }
1097:   PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);
1098:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its));
1099:   PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
1100:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax));
1101:   PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
1102:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax));

1104:   PetscOptionsInt("-pc_pfmg_max_levels","Max Levels for MG hierarchy","HYPRE_StructPFMGSetMaxLevels",ex->max_levels,&ex->max_levels,NULL);
1105:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxLevels,(ex->hsolver,ex->max_levels));

1107:   PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);
1108:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol));
1109:   PetscOptionsEList("-pc_pfmg_relax_type","Relax type for the up and down cycles","HYPRE_StructPFMGSetRelaxType",PFMGRelaxType,ALEN(PFMGRelaxType),PFMGRelaxType[ex->relax_type],&ex->relax_type,NULL);
1110:   PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type));
1111:   PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);
1112:   PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type));
1113:   PetscOptionsTail();
1114:   return(0);
1115: }

1119: PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
1120: {
1121:   PetscErrorCode  ierr;
1122:   PC_PFMG         *ex = (PC_PFMG*) pc->data;
1123:   PetscScalar     *xx,*yy;
1124:   PetscInt        ilower[3],iupper[3];
1125:   Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);

1128:   PetscCitationsRegister(hypreCitation,&cite);
1129:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
1130:   iupper[0] += ilower[0] - 1;
1131:   iupper[1] += ilower[1] - 1;
1132:   iupper[2] += ilower[2] - 1;

1134:   /* copy x values over to hypre */
1135:   PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
1136:   VecGetArray(x,&xx);
1137:   PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,xx));
1138:   VecRestoreArray(x,&xx);
1139:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
1140:   PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx));

1142:   /* copy solution values back to PETSc */
1143:   VecGetArray(y,&yy);
1144:   PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,yy));
1145:   VecRestoreArray(y,&yy);
1146:   return(0);
1147: }

1151: static PetscErrorCode PCApplyRichardson_PFMG(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
1152: {
1153:   PC_PFMG        *jac = (PC_PFMG*)pc->data;
1155:   PetscInt       oits;

1158:   PetscCitationsRegister(hypreCitation,&cite);
1159:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its));
1160:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol));

1162:   PCApply_PFMG(pc,b,y);
1163:   PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,(HYPRE_Int *)&oits));
1164:   *outits = oits;
1165:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1166:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
1167:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol));
1168:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its));
1169:   return(0);
1170: }


1175: PetscErrorCode PCSetUp_PFMG(PC pc)
1176: {
1177:   PetscErrorCode  ierr;
1178:   PC_PFMG         *ex = (PC_PFMG*) pc->data;
1179:   Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
1180:   PetscBool       flg;

1183:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESTRUCT,&flg);
1184:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESTRUCT with this preconditioner");

1186:   /* create the hypre solver object and set its information */
1187:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
1188:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1189:   PCSetFromOptions_PFMG(pc);
1190:   PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
1191:   PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver));
1192:   return(0);
1193: }


1196: /*MC
1197:      PCPFMG - the hypre PFMG multigrid solver

1199:    Level: advanced

1201:    Options Database:
1202: + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner
1203: . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
1204: . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid
1205: . -pc_pfmg_tol <tol> tolerance of PFMG
1206: . -pc_pfmg_relax_type -relaxation type for the up and down cycles, one of Jacobi,Weighted-Jacobi,symmetric-Red/Black-Gauss-Seidel,Red/Black-Gauss-Seidel
1207: - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin

1209:    Notes:  This is for CELL-centered descretizations

1211:            This must be used with the MATHYPRESTRUCT matrix type.
1212:            This is less general than in hypre, it supports only one block per process defined by a PETSc DMDA.

1214: .seealso:  PCMG, MATHYPRESTRUCT
1215: M*/

1219: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
1220: {
1222:   PC_PFMG        *ex;

1225:   PetscNew(&ex); \
1226:   pc->data = ex;

1228:   ex->its            = 1;
1229:   ex->tol            = 1.e-8;
1230:   ex->relax_type     = 1;
1231:   ex->rap_type       = 0;
1232:   ex->num_pre_relax  = 1;
1233:   ex->num_post_relax = 1;
1234:   ex->max_levels     = 0;

1236:   pc->ops->setfromoptions  = PCSetFromOptions_PFMG;
1237:   pc->ops->view            = PCView_PFMG;
1238:   pc->ops->destroy         = PCDestroy_PFMG;
1239:   pc->ops->apply           = PCApply_PFMG;
1240:   pc->ops->applyrichardson = PCApplyRichardson_PFMG;
1241:   pc->ops->setup           = PCSetUp_PFMG;

1243:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
1244:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
1245:   return(0);
1246: }

1248: /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/

1250: /* we know we are working with a HYPRE_SStructMatrix */
1251: typedef struct {
1252:   MPI_Comm            hcomm;       /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
1253:   HYPRE_SStructSolver ss_solver;

1255:   /* keep copy of SYSPFMG options used so may view them */
1256:   PetscInt its;
1257:   double   tol;
1258:   PetscInt relax_type;
1259:   PetscInt num_pre_relax,num_post_relax;
1260: } PC_SysPFMG;

1264: PetscErrorCode PCDestroy_SysPFMG(PC pc)
1265: {
1267:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

1270:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
1271:   MPI_Comm_free(&ex->hcomm);
1272:   PetscFree(pc->data);
1273:   return(0);
1274: }

1276: static const char *SysPFMGRelaxType[] = {"Weighted-Jacobi","Red/Black-Gauss-Seidel"};

1280: PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
1281: {
1283:   PetscBool      iascii;
1284:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

1287:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1288:   if (iascii) {
1289:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG preconditioning\n");
1290:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: max iterations %d\n",ex->its);
1291:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: tolerance %g\n",ex->tol);
1292:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: relax type %s\n",PFMGRelaxType[ex->relax_type]);
1293:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG: number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
1294:   }
1295:   return(0);
1296: }


1301: PetscErrorCode PCSetFromOptions_SysPFMG(PC pc)
1302: {
1304:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;
1305:   PetscBool      flg = PETSC_FALSE;

1308:   PetscOptionsHead("SysPFMG options");
1309:   PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);
1310:   if (flg) {
1311:     int level=3;
1312:     PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level));
1313:   }
1314:   PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);
1315:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its));
1316:   PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
1317:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax));
1318:   PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
1319:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax));

1321:   PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);
1322:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol));
1323:   PetscOptionsEList("-pc_syspfmg_relax_type","Relax type for the up and down cycles","HYPRE_SStructSysPFMGSetRelaxType",SysPFMGRelaxType,4,SysPFMGRelaxType[ex->relax_type],&ex->relax_type,NULL);
1324:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type));
1325:   PetscOptionsTail();
1326:   return(0);
1327: }

1331: PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
1332: {
1333:   PetscErrorCode   ierr;
1334:   PC_SysPFMG       *ex = (PC_SysPFMG*) pc->data;
1335:   PetscScalar      *xx,*yy;
1336:   PetscInt         ilower[3],iupper[3];
1337:   Mat_HYPRESStruct *mx     = (Mat_HYPRESStruct*)(pc->pmat->data);
1338:   PetscInt         ordering= mx->dofs_order;
1339:   PetscInt         nvars   = mx->nvars;
1340:   PetscInt         part    = 0;
1341:   PetscInt         size;
1342:   PetscInt         i;

1345:   PetscCitationsRegister(hypreCitation,&cite);
1346:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
1347:   iupper[0] += ilower[0] - 1;
1348:   iupper[1] += ilower[1] - 1;
1349:   iupper[2] += ilower[2] - 1;

1351:   size = 1;
1352:   for (i= 0; i< 3; i++) size *= (iupper[i]-ilower[i]+1);

1354:   /* copy x values over to hypre for variable ordering */
1355:   if (ordering) {
1356:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
1357:     VecGetArray(x,&xx);
1358:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,xx+(size*i)));
1359:     VecRestoreArray(x,&xx);
1360:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
1361:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
1362:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

1364:     /* copy solution values back to PETSc */
1365:     VecGetArray(y,&yy);
1366:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,yy+(size*i)));
1367:     VecRestoreArray(y,&yy);
1368:   } else {      /* nodal ordering must be mapped to variable ordering for sys_pfmg */
1369:     PetscScalar *z;
1370:     PetscInt    j, k;

1372:     PetscMalloc1(nvars*size,&z);
1373:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
1374:     VecGetArray(x,&xx);

1376:     /* transform nodal to hypre's variable ordering for sys_pfmg */
1377:     for (i= 0; i< size; i++) {
1378:       k= i*nvars;
1379:       for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
1380:     }
1381:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
1382:     VecRestoreArray(x,&xx);
1383:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
1384:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

1386:     /* copy solution values back to PETSc */
1387:     VecGetArray(y,&yy);
1388:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,(HYPRE_Int *)ilower,(HYPRE_Int *)iupper,i,z+(size*i)));
1389:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
1390:     for (i= 0; i< size; i++) {
1391:       k= i*nvars;
1392:       for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
1393:     }
1394:     VecRestoreArray(y,&yy);
1395:     PetscFree(z);
1396:   }
1397:   return(0);
1398: }

1402: static PetscErrorCode PCApplyRichardson_SysPFMG(PC pc,Vec b,Vec y,Vec w,PetscReal rtol,PetscReal abstol, PetscReal dtol,PetscInt its,PetscBool guesszero,PetscInt *outits,PCRichardsonConvergedReason *reason)
1403: {
1404:   PC_SysPFMG     *jac = (PC_SysPFMG*)pc->data;
1406:   PetscInt       oits;

1409:   PetscCitationsRegister(hypreCitation,&cite);
1410:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its));
1411:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol));
1412:   PCApply_SysPFMG(pc,b,y);
1413:   PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,(HYPRE_Int *)&oits));
1414:   *outits = oits;
1415:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1416:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
1417:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol));
1418:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its));
1419:   return(0);
1420: }


1425: PetscErrorCode PCSetUp_SysPFMG(PC pc)
1426: {
1427:   PetscErrorCode   ierr;
1428:   PC_SysPFMG       *ex = (PC_SysPFMG*) pc->data;
1429:   Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
1430:   PetscBool        flg;

1433:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRESSTRUCT,&flg);
1434:   if (!flg) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_INCOMP,"Must use MATHYPRESSTRUCT with this preconditioner");

1436:   /* create the hypre sstruct solver object and set its information */
1437:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
1438:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
1439:   PCSetFromOptions_SysPFMG(pc);
1440:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver));
1441:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
1442:   return(0);
1443: }


1446: /*MC
1447:      PCSysPFMG - the hypre SysPFMG multigrid solver

1449:    Level: advanced

1451:    Options Database:
1452: + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner
1453: . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
1454: . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid
1455: . -pc_syspfmg_tol <tol> tolerance of SysPFMG
1456: . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel

1458:    Notes:  This is for CELL-centered descretizations

1460:            This must be used with the MATHYPRESSTRUCT matrix type.
1461:            This is less general than in hypre, it supports only one part, and one block per process defined by a PETSc DMDA.
1462:            Also, only cell-centered variables.

1464: .seealso:  PCMG, MATHYPRESSTRUCT
1465: M*/

1469: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
1470: {
1472:   PC_SysPFMG     *ex;

1475:   PetscNew(&ex); \
1476:   pc->data = ex;

1478:   ex->its            = 1;
1479:   ex->tol            = 1.e-8;
1480:   ex->relax_type     = 1;
1481:   ex->num_pre_relax  = 1;
1482:   ex->num_post_relax = 1;

1484:   pc->ops->setfromoptions  = PCSetFromOptions_SysPFMG;
1485:   pc->ops->view            = PCView_SysPFMG;
1486:   pc->ops->destroy         = PCDestroy_SysPFMG;
1487:   pc->ops->apply           = PCApply_SysPFMG;
1488:   pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
1489:   pc->ops->setup           = PCSetUp_SysPFMG;

1491:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
1492:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
1493:   return(0);
1494: }