Actual source code: hypre.c

petsc-3.15.0 2021-04-05
Report Typos and Errors
  1: /*
  2:    Provides an interface to the LLNL package hypre
  3: */

  5: #include <petscpkg_version.h>
  6: #include <petsc/private/pcimpl.h>
  7: /* this include is needed ONLY to allow access to the private data inside the Mat object specific to hypre */
  8: #include <petsc/private/matimpl.h>
  9: #include <../src/vec/vec/impls/hypre/vhyp.h>
 10: #include <../src/mat/impls/hypre/mhypre.h>
 11: #include <../src/dm/impls/da/hypre/mhyp.h>
 12: #include <_hypre_parcsr_ls.h>
 13: #include <petscmathypre.h>

 15: static PetscBool cite = PETSC_FALSE;
 16: static const char hypreCitation[] = "@manual{hypre-web-page,\n  title  = {{\\sl hypre}: High Performance Preconditioners},\n  organization = {Lawrence Livermore National Laboratory},\n  note  = {\\url{https://computation.llnl.gov/projects/hypre-scalable-linear-solvers-multigrid-methods}}\n}\n";

 18: /*
 19:    Private context (data structure) for the  preconditioner.
 20: */
 21: typedef struct {
 22:   HYPRE_Solver   hsolver;
 23:   Mat            hpmat; /* MatHYPRE */

 25:   HYPRE_Int (*destroy)(HYPRE_Solver);
 26:   HYPRE_Int (*solve)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
 27:   HYPRE_Int (*setup)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);

 29:   MPI_Comm comm_hypre;
 30:   char     *hypre_type;

 32:   /* options for Pilut and BoomerAMG*/
 33:   PetscInt  maxiter;
 34:   PetscReal tol;

 36:   /* options for Pilut */
 37:   PetscInt factorrowsize;

 39:   /* options for ParaSails */
 40:   PetscInt  nlevels;
 41:   PetscReal threshold;
 42:   PetscReal filter;
 43:   PetscReal loadbal;
 44:   PetscInt  logging;
 45:   PetscInt  ruse;
 46:   PetscInt  symt;

 48:   /* options for BoomerAMG */
 49:   PetscBool printstatistics;

 51:   /* options for BoomerAMG */
 52:   PetscInt  cycletype;
 53:   PetscInt  maxlevels;
 54:   PetscReal strongthreshold;
 55:   PetscReal maxrowsum;
 56:   PetscInt  gridsweeps[3];
 57:   PetscInt  coarsentype;
 58:   PetscInt  measuretype;
 59:   PetscInt  smoothtype;
 60:   PetscInt  smoothnumlevels;
 61:   PetscInt  eu_level;   /* Number of levels for ILU(k) in Euclid */
 62:   PetscReal eu_droptolerance; /* Drop tolerance for ILU(k) in Euclid */
 63:   PetscInt  eu_bj;      /* Defines use of Block Jacobi ILU in Euclid */
 64:   PetscInt  relaxtype[3];
 65:   PetscReal relaxweight;
 66:   PetscReal outerrelaxweight;
 67:   PetscInt  relaxorder;
 68:   PetscReal truncfactor;
 69:   PetscBool applyrichardson;
 70:   PetscInt  pmax;
 71:   PetscInt  interptype;
 72:   PetscInt  maxc;
 73:   PetscInt  minc;

 75:   /* AIR */
 76:   PetscInt  Rtype;
 77:   PetscReal Rstrongthreshold;
 78:   PetscReal Rfilterthreshold;
 79:   PetscInt  Adroptype;
 80:   PetscReal Adroptol;

 82:   PetscInt  agg_nl;
 83:   PetscInt  agg_num_paths;
 84:   PetscBool nodal_relax;
 85:   PetscInt  nodal_relax_levels;
 86:   PetscBool keeptranspose;

 88:   PetscInt  nodal_coarsening;
 89:   PetscInt  nodal_coarsening_diag;
 90:   PetscInt  vec_interp_variant;
 91:   PetscInt  vec_interp_qmax;
 92:   PetscBool vec_interp_smooth;
 93:   PetscInt  interp_refine;

 95:   HYPRE_IJVector  *hmnull;
 96:   HYPRE_ParVector *phmnull;  /* near null space passed to hypre */
 97:   PetscInt        n_hmnull;
 98:   Vec             hmnull_constant;
 99:   HYPRE_Complex   **hmnull_hypre_data_array;   /* this is the space in hmnull that was allocated by hypre, it is restored to hypre just before freeing the phmnull vectors */

101:   /* options for AS (Auxiliary Space preconditioners) */
102:   PetscInt  as_print;
103:   PetscInt  as_max_iter;
104:   PetscReal as_tol;
105:   PetscInt  as_relax_type;
106:   PetscInt  as_relax_times;
107:   PetscReal as_relax_weight;
108:   PetscReal as_omega;
109:   PetscInt  as_amg_alpha_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for vector Poisson (AMS) or Curl problem (ADS) */
110:   PetscReal as_amg_alpha_theta;   /* AMG strength for vector Poisson (AMS) or Curl problem (ADS) */
111:   PetscInt  as_amg_beta_opts[5];  /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for scalar Poisson (AMS) or vector Poisson (ADS) */
112:   PetscReal as_amg_beta_theta;    /* AMG strength for scalar Poisson (AMS) or vector Poisson (ADS)  */
113:   PetscInt  ams_cycle_type;
114:   PetscInt  ads_cycle_type;

116:   /* additional data */
117:   Mat G;             /* MatHYPRE */
118:   Mat C;             /* MatHYPRE */
119:   Mat alpha_Poisson; /* MatHYPRE */
120:   Mat beta_Poisson;  /* MatHYPRE */

122:   /* extra information for AMS */
123:   PetscInt       dim; /* geometrical dimension */
124:   HYPRE_IJVector coords[3];
125:   HYPRE_IJVector constants[3];
126:   Mat            RT_PiFull, RT_Pi[3];
127:   Mat            ND_PiFull, ND_Pi[3];
128:   PetscBool      ams_beta_is_zero;
129:   PetscBool      ams_beta_is_zero_part;
130:   PetscInt       ams_proj_freq;
131: } PC_HYPRE;

133: PetscErrorCode PCHYPREGetSolver(PC pc,HYPRE_Solver *hsolver)
134: {
135:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

138:   *hsolver = jac->hsolver;
139:   return(0);
140: }

142: /*
143:   Matrices with AIJ format are created IN PLACE with using (I,J,data) from BoomerAMG. Since the data format in hypre_ParCSRMatrix
144:   is different from that used in PETSc, the original hypre_ParCSRMatrix can not be used any more after call this routine.
145:   It is used in PCHMG. Other users should avoid using this function.
146: */
147: static PetscErrorCode PCGetCoarseOperators_BoomerAMG(PC pc,PetscInt *nlevels,Mat *operators[])
148: {
149:   PC_HYPRE             *jac  = (PC_HYPRE*)pc->data;
150:   PetscBool            same = PETSC_FALSE;
151:   PetscErrorCode       ierr;
152:   PetscInt             num_levels,l;
153:   Mat                  *mattmp;
154:   hypre_ParCSRMatrix   **A_array;

157:   PetscStrcmp(jac->hypre_type,"boomeramg",&same);
158:   if (!same) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_NOTSAMETYPE,"Hypre type is not BoomerAMG \n");
159:   num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData*) (jac->hsolver));
160:   PetscMalloc1(num_levels,&mattmp);
161:   A_array    = hypre_ParAMGDataAArray((hypre_ParAMGData*) (jac->hsolver));
162:   for (l=1; l<num_levels; l++) {
163:     MatCreateFromParCSR(A_array[l],MATAIJ,PETSC_OWN_POINTER, &(mattmp[num_levels-1-l]));
164:     /* We want to own the data, and HYPRE can not touch this matrix any more */
165:     A_array[l] = NULL;
166:   }
167:   *nlevels = num_levels;
168:   *operators = mattmp;
169:   return(0);
170: }

172: /*
173:   Matrices with AIJ format are created IN PLACE with using (I,J,data) from BoomerAMG. Since the data format in hypre_ParCSRMatrix
174:   is different from that used in PETSc, the original hypre_ParCSRMatrix can not be used any more after call this routine.
175:   It is used in PCHMG. Other users should avoid using this function.
176: */
177: static PetscErrorCode PCGetInterpolations_BoomerAMG(PC pc,PetscInt *nlevels,Mat *interpolations[])
178: {
179:   PC_HYPRE             *jac  = (PC_HYPRE*)pc->data;
180:   PetscBool            same = PETSC_FALSE;
181:   PetscErrorCode       ierr;
182:   PetscInt             num_levels,l;
183:   Mat                  *mattmp;
184:   hypre_ParCSRMatrix   **P_array;

187:   PetscStrcmp(jac->hypre_type,"boomeramg",&same);
188:   if (!same) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_NOTSAMETYPE,"Hypre type is not BoomerAMG \n");
189:   num_levels = hypre_ParAMGDataNumLevels((hypre_ParAMGData*) (jac->hsolver));
190:   PetscMalloc1(num_levels,&mattmp);
191:   P_array  = hypre_ParAMGDataPArray((hypre_ParAMGData*) (jac->hsolver));
192:   for (l=1; l<num_levels; l++) {
193:     MatCreateFromParCSR(P_array[num_levels-1-l],MATAIJ,PETSC_OWN_POINTER, &(mattmp[l-1]));
194:     /* We want to own the data, and HYPRE can not touch this matrix any more */
195:     P_array[num_levels-1-l] = NULL;
196:   }
197:   *nlevels = num_levels;
198:   *interpolations = mattmp;
199:   return(0);
200: }

202: /* Resets (frees) Hypre's representation of the near null space */
203: static PetscErrorCode PCHYPREResetNearNullSpace_Private(PC pc)
204: {
205:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
206:   PetscInt       i;

210:   for (i=0; i<jac->n_hmnull; i++) {
211:     PETSC_UNUSED HYPRE_Complex *harray;
212:     VecHYPRE_ParVectorReplacePointer(jac->hmnull[i],jac->hmnull_hypre_data_array[i],harray);
213:     PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->hmnull[i]));
214:   }
215:   PetscFree(jac->hmnull);
216:   PetscFree(jac->hmnull_hypre_data_array);
217:   PetscFree(jac->phmnull);
218:   VecDestroy(&jac->hmnull_constant);
219:   jac->n_hmnull = 0;
220:   return(0);
221: }

223: static PetscErrorCode PCSetUp_HYPRE(PC pc)
224: {
225:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
226:   Mat_HYPRE          *hjac;
227:   HYPRE_ParCSRMatrix hmat;
228:   HYPRE_ParVector    bv,xv;
229:   PetscBool          ishypre;
230:   PetscErrorCode     ierr;

233:   if (!jac->hypre_type) {
234:     PCHYPRESetType(pc,"boomeramg");
235:   }

237:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRE,&ishypre);
238:   if (!ishypre) {
239:     MatDestroy(&jac->hpmat);
240:     MatConvert(pc->pmat,MATHYPRE,MAT_INITIAL_MATRIX,&jac->hpmat);
241:     PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->hpmat);
242:   } else {
243:     PetscObjectReference((PetscObject)pc->pmat);
244:     MatDestroy(&jac->hpmat);
245:     jac->hpmat = pc->pmat;
246:   }
247:   hjac = (Mat_HYPRE*)(jac->hpmat->data);

249:   /* special case for BoomerAMG */
250:   if (jac->setup == HYPRE_BoomerAMGSetup) {
251:     MatNullSpace    mnull;
252:     PetscBool       has_const;
253:     PetscInt        bs,nvec,i;
254:     const Vec       *vecs;
255:     HYPRE_Complex   *petscvecarray;

257:     MatGetBlockSize(pc->pmat,&bs);
258:     if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
259:     MatGetNearNullSpace(pc->mat, &mnull);
260:     if (mnull) {
261:       PCHYPREResetNearNullSpace_Private(pc);
262:       MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);
263:       PetscMalloc1(nvec+1,&jac->hmnull);
264:       PetscMalloc1(nvec+1,&jac->hmnull_hypre_data_array);
265:       PetscMalloc1(nvec+1,&jac->phmnull);
266:       for (i=0; i<nvec; i++) {
267:         VecHYPRE_IJVectorCreate(vecs[i],&jac->hmnull[i]);
268:         VecGetArrayRead(vecs[i],(const PetscScalar **)&petscvecarray);
269:         VecHYPRE_ParVectorReplacePointer(jac->hmnull[i],petscvecarray,jac->hmnull_hypre_data_array[i]);
270:         VecRestoreArrayRead(vecs[i],(const PetscScalar **)&petscvecarray);
271:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[i],(void**)&jac->phmnull[i]));
272:       }
273:       if (has_const) {
274:         MatCreateVecs(pc->pmat,&jac->hmnull_constant,NULL);
275:         PetscLogObjectParent((PetscObject)pc,(PetscObject)jac->hmnull_constant);
276:         VecSet(jac->hmnull_constant,1);
277:         VecNormalize(jac->hmnull_constant,NULL);
278:         VecHYPRE_IJVectorCreate(jac->hmnull_constant,&jac->hmnull[nvec]);
279:         VecGetArrayRead(jac->hmnull_constant,(const PetscScalar **)&petscvecarray);
280:         VecHYPRE_ParVectorReplacePointer(jac->hmnull[nvec],petscvecarray,jac->hmnull_hypre_data_array[nvec]);
281:         VecRestoreArrayRead(jac->hmnull_constant,(const PetscScalar **)&petscvecarray);
282:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[nvec],(void**)&jac->phmnull[nvec]));
283:         nvec++;
284:       }
285:       PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVectors,(jac->hsolver,nvec,jac->phmnull));
286:       jac->n_hmnull = nvec;
287:     }
288:   }

290:   /* special case for AMS */
291:   if (jac->setup == HYPRE_AMSSetup) {
292:     Mat_HYPRE          *hm;
293:     HYPRE_ParCSRMatrix parcsr;
294:     if (!jac->coords[0] && !jac->constants[0] && !(jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]))) {
295:       SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs either the coordinate vectors via PCSetCoordinates() or the edge constant vectors via PCHYPRESetEdgeConstantVectors() or the interpolation matrix via PCHYPRESetInterpolations");
296:     }
297:     if (jac->dim) {
298:       PetscStackCallStandard(HYPRE_AMSSetDimension,(jac->hsolver,jac->dim));
299:     }
300:     if (jac->constants[0]) {
301:       HYPRE_ParVector ozz,zoz,zzo = NULL;
302:       PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[0],(void**)(&ozz)));
303:       PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[1],(void**)(&zoz)));
304:       if (jac->constants[2]) {
305:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->constants[2],(void**)(&zzo)));
306:       }
307:       PetscStackCallStandard(HYPRE_AMSSetEdgeConstantVectors,(jac->hsolver,ozz,zoz,zzo));
308:     }
309:     if (jac->coords[0]) {
310:       HYPRE_ParVector coords[3];
311:       coords[0] = NULL;
312:       coords[1] = NULL;
313:       coords[2] = NULL;
314:       if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0],(void**)(&coords[0])));
315:       if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1],(void**)(&coords[1])));
316:       if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2],(void**)(&coords[2])));
317:       PetscStackCallStandard(HYPRE_AMSSetCoordinateVectors,(jac->hsolver,coords[0],coords[1],coords[2]));
318:     }
319:     if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE AMS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient");
320:     hm = (Mat_HYPRE*)(jac->G->data);
321:     PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
322:     PetscStackCallStandard(HYPRE_AMSSetDiscreteGradient,(jac->hsolver,parcsr));
323:     if (jac->alpha_Poisson) {
324:       hm = (Mat_HYPRE*)(jac->alpha_Poisson->data);
325:       PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
326:       PetscStackCallStandard(HYPRE_AMSSetAlphaPoissonMatrix,(jac->hsolver,parcsr));
327:     }
328:     if (jac->ams_beta_is_zero) {
329:       PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,NULL));
330:     } else if (jac->beta_Poisson) {
331:       hm = (Mat_HYPRE*)(jac->beta_Poisson->data);
332:       PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
333:       PetscStackCallStandard(HYPRE_AMSSetBetaPoissonMatrix,(jac->hsolver,parcsr));
334:     }
335:     if (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1])) {
336:       PetscInt           i;
337:       HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3];
338:       if (jac->ND_PiFull) {
339:         hm = (Mat_HYPRE*)(jac->ND_PiFull->data);
340:         PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&nd_parcsrfull)));
341:       } else {
342:         nd_parcsrfull = NULL;
343:       }
344:       for (i=0;i<3;++i) {
345:         if (jac->ND_Pi[i]) {
346:           hm = (Mat_HYPRE*)(jac->ND_Pi[i]->data);
347:           PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&nd_parcsr[i])));
348:         } else {
349:           nd_parcsr[i] = NULL;
350:         }
351:       }
352:       PetscStackCallStandard(HYPRE_AMSSetInterpolations,(jac->hsolver,nd_parcsrfull,nd_parcsr[0],nd_parcsr[1],nd_parcsr[2]));
353:     }
354:   }
355:   /* special case for ADS */
356:   if (jac->setup == HYPRE_ADSSetup) {
357:     Mat_HYPRE          *hm;
358:     HYPRE_ParCSRMatrix parcsr;
359:     if (!jac->coords[0] && !((jac->RT_PiFull || (jac->RT_Pi[0] && jac->RT_Pi[1])) && (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1])))) {
360:       SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs either the coordinate vectors via PCSetCoordinates() or the interpolation matrices via PCHYPRESetInterpolations");
361:     }
362:     else if (!jac->coords[1] || !jac->coords[2]) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner has been designed for three dimensional problems! For two dimensional problems, use HYPRE AMS instead");
363:     if (!jac->G) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs the discrete gradient operator via PCHYPRESetDiscreteGradient");
364:     if (!jac->C) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_USER,"HYPRE ADS preconditioner needs the discrete curl operator via PCHYPRESetDiscreteGradient");
365:     if (jac->coords[0]) {
366:       HYPRE_ParVector coords[3];
367:       coords[0] = NULL;
368:       coords[1] = NULL;
369:       coords[2] = NULL;
370:       if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[0],(void**)(&coords[0])));
371:       if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[1],(void**)(&coords[1])));
372:       if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->coords[2],(void**)(&coords[2])));
373:       PetscStackCallStandard(HYPRE_ADSSetCoordinateVectors,(jac->hsolver,coords[0],coords[1],coords[2]));
374:     }
375:     hm = (Mat_HYPRE*)(jac->G->data);
376:     PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
377:     PetscStackCallStandard(HYPRE_ADSSetDiscreteGradient,(jac->hsolver,parcsr));
378:     hm = (Mat_HYPRE*)(jac->C->data);
379:     PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&parcsr)));
380:     PetscStackCallStandard(HYPRE_ADSSetDiscreteCurl,(jac->hsolver,parcsr));
381:     if ((jac->RT_PiFull || (jac->RT_Pi[0] && jac->RT_Pi[1])) && (jac->ND_PiFull || (jac->ND_Pi[0] && jac->ND_Pi[1]))) {
382:       PetscInt           i;
383:       HYPRE_ParCSRMatrix rt_parcsrfull, rt_parcsr[3];
384:       HYPRE_ParCSRMatrix nd_parcsrfull, nd_parcsr[3];
385:       if (jac->RT_PiFull) {
386:         hm = (Mat_HYPRE*)(jac->RT_PiFull->data);
387:         PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&rt_parcsrfull)));
388:       } else {
389:         rt_parcsrfull = NULL;
390:       }
391:       for (i=0;i<3;++i) {
392:         if (jac->RT_Pi[i]) {
393:           hm = (Mat_HYPRE*)(jac->RT_Pi[i]->data);
394:           PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&rt_parcsr[i])));
395:         } else {
396:           rt_parcsr[i] = NULL;
397:         }
398:       }
399:       if (jac->ND_PiFull) {
400:         hm = (Mat_HYPRE*)(jac->ND_PiFull->data);
401:         PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&nd_parcsrfull)));
402:       } else {
403:         nd_parcsrfull = NULL;
404:       }
405:       for (i=0;i<3;++i) {
406:         if (jac->ND_Pi[i]) {
407:           hm = (Mat_HYPRE*)(jac->ND_Pi[i]->data);
408:           PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hm->ij,(void**)(&nd_parcsr[i])));
409:         } else {
410:           nd_parcsr[i] = NULL;
411:         }
412:       }
413:       PetscStackCallStandard(HYPRE_ADSSetInterpolations,(jac->hsolver,rt_parcsrfull,rt_parcsr[0],rt_parcsr[1],rt_parcsr[2],nd_parcsrfull,nd_parcsr[0],nd_parcsr[1],nd_parcsr[2]));
414:     }
415:   }
416:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
417:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b,(void**)&bv));
418:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x,(void**)&xv));
419:   PetscStackCallStandard(jac->setup,(jac->hsolver,hmat,bv,xv));
420:   return(0);
421: }

423: static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
424: {
425:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
426:   Mat_HYPRE          *hjac = (Mat_HYPRE*)(jac->hpmat->data);
427:   PetscErrorCode     ierr;
428:   HYPRE_ParCSRMatrix hmat;
429:   HYPRE_Complex      *xv,*sxv;
430:   HYPRE_Complex      *bv,*sbv;
431:   HYPRE_ParVector    jbv,jxv;
432:   PetscInt           hierr;

435:   PetscCitationsRegister(hypreCitation,&cite);
436:   if (!jac->applyrichardson) {VecSet(x,0.0);}
437:   VecGetArrayRead(b,(const PetscScalar **)&bv);
438:   VecGetArrayWrite(x,(PetscScalar **)&xv);
439:   VecHYPRE_ParVectorReplacePointer(hjac->b,bv,sbv);
440:   VecHYPRE_ParVectorReplacePointer(hjac->x,xv,sxv);

442:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
443:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b,(void**)&jbv));
444:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x,(void**)&jxv));
445:   PetscStackCall("Hypre solve",h(*jac->solve)(jac->hsolver,hmat,jbv,jxv);
446:   if (hierr && hierr != HYPRE_ERROR_CONV) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
447:   if (hierr) hypre__global_error = 0;);

449:   if (jac->setup == HYPRE_AMSSetup && jac->ams_beta_is_zero_part) {
450:     PetscStackCallStandard(HYPRE_AMSProjectOutGradients,(jac->hsolver,jxv));
451:   }
452:   VecHYPRE_ParVectorReplacePointer(hjac->b,sbv,bv);
453:   VecHYPRE_ParVectorReplacePointer(hjac->x,sxv,xv);
454:   VecRestoreArrayWrite(x,(PetscScalar **)&xv);
455:   VecRestoreArrayRead(b,(const PetscScalar **)&bv);
456:   return(0);
457: }

459: static PetscErrorCode PCReset_HYPRE(PC pc)
460: {
461:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;

465:   MatDestroy(&jac->hpmat);
466:   MatDestroy(&jac->G);
467:   MatDestroy(&jac->C);
468:   MatDestroy(&jac->alpha_Poisson);
469:   MatDestroy(&jac->beta_Poisson);
470:   MatDestroy(&jac->RT_PiFull);
471:   MatDestroy(&jac->RT_Pi[0]);
472:   MatDestroy(&jac->RT_Pi[1]);
473:   MatDestroy(&jac->RT_Pi[2]);
474:   MatDestroy(&jac->ND_PiFull);
475:   MatDestroy(&jac->ND_Pi[0]);
476:   MatDestroy(&jac->ND_Pi[1]);
477:   MatDestroy(&jac->ND_Pi[2]);
478:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0])); jac->coords[0] = NULL;
479:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1])); jac->coords[1] = NULL;
480:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2])); jac->coords[2] = NULL;
481:   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0])); jac->constants[0] = NULL;
482:   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1])); jac->constants[1] = NULL;
483:   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2])); jac->constants[2] = NULL;
484:   PCHYPREResetNearNullSpace_Private(pc);
485:   jac->ams_beta_is_zero = PETSC_FALSE;
486:   jac->dim = 0;
487:   return(0);
488: }

490: static PetscErrorCode PCDestroy_HYPRE(PC pc)
491: {
492:   PC_HYPRE                 *jac = (PC_HYPRE*)pc->data;
493:   PetscErrorCode           ierr;

496:   PCReset_HYPRE(pc);
497:   if (jac->destroy) PetscStackCallStandard(jac->destroy,(jac->hsolver));
498:   PetscFree(jac->hypre_type);
499:   if (jac->comm_hypre != MPI_COMM_NULL) {MPI_Comm_free(&(jac->comm_hypre));}
500:   PetscFree(pc->data);

502:   PetscObjectChangeTypeName((PetscObject)pc,0);
503:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);
504:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);
505:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetCoordinates_C",NULL);
506:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",NULL);
507:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",NULL);
508:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",NULL);
509:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetConstantEdgeVectors_C",NULL);
510:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",NULL);
511:   PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",NULL);
512:   PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",NULL);
513:   return(0);
514: }

516: /* --------------------------------------------------------------------------------------------*/
517: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PetscOptionItems *PetscOptionsObject,PC pc)
518: {
519:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
521:   PetscBool      flag;

524:   PetscOptionsHead(PetscOptionsObject,"HYPRE Pilut Options");
525:   PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
526:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter));
527:   PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
528:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol));
529:   PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
530:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize));
531:   PetscOptionsTail();
532:   return(0);
533: }

535: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
536: {
537:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
539:   PetscBool      iascii;

542:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
543:   if (iascii) {
544:     PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut preconditioning\n");
545:     if (jac->maxiter != PETSC_DEFAULT) {
546:       PetscViewerASCIIPrintf(viewer,"    maximum number of iterations %d\n",jac->maxiter);
547:     } else {
548:       PetscViewerASCIIPrintf(viewer,"    default maximum number of iterations \n");
549:     }
550:     if (jac->tol != PETSC_DEFAULT) {
551:       PetscViewerASCIIPrintf(viewer,"    drop tolerance %g\n",(double)jac->tol);
552:     } else {
553:       PetscViewerASCIIPrintf(viewer,"    default drop tolerance \n");
554:     }
555:     if (jac->factorrowsize != PETSC_DEFAULT) {
556:       PetscViewerASCIIPrintf(viewer,"    factor row size %d\n",jac->factorrowsize);
557:     } else {
558:       PetscViewerASCIIPrintf(viewer,"    default factor row size \n");
559:     }
560:   }
561:   return(0);
562: }

564: /* --------------------------------------------------------------------------------------------*/
565: static PetscErrorCode PCSetFromOptions_HYPRE_Euclid(PetscOptionItems *PetscOptionsObject,PC pc)
566: {
567:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
569:   PetscBool      flag,eu_bj = jac->eu_bj ? PETSC_TRUE : PETSC_FALSE;

572:   PetscOptionsHead(PetscOptionsObject,"HYPRE Euclid Options");
573:   PetscOptionsInt("-pc_hypre_euclid_level","Factorization levels","None",jac->eu_level,&jac->eu_level,&flag);
574:   if (flag) PetscStackCallStandard(HYPRE_EuclidSetLevel,(jac->hsolver,jac->eu_level));

576:   PetscOptionsReal("-pc_hypre_euclid_droptolerance","Drop tolerance for ILU(k) in Euclid","None",jac->eu_droptolerance,&jac->eu_droptolerance,&flag);
577:   if (flag){
578:     PetscMPIInt size;

580:     MPI_Comm_size(PetscObjectComm((PetscObject)pc),&size);
581:     if (size > 1) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"hypre's Euclid does not support a parallel drop tolerance");
582:     PetscStackCallStandard(HYPRE_EuclidSetILUT,(jac->hsolver,jac->eu_droptolerance));
583:   }

585:   PetscOptionsBool("-pc_hypre_euclid_bj", "Use Block Jacobi for ILU in Euclid", "None", eu_bj,&eu_bj,&flag);
586:   if (flag) {
587:     jac->eu_bj = eu_bj ? 1 : 0;
588:     PetscStackCallStandard(HYPRE_EuclidSetBJ,(jac->hsolver,jac->eu_bj));
589:   }
590:   PetscOptionsTail();
591:   return(0);
592: }

594: static PetscErrorCode PCView_HYPRE_Euclid(PC pc,PetscViewer viewer)
595: {
596:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
598:   PetscBool      iascii;

601:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
602:   if (iascii) {
603:     PetscViewerASCIIPrintf(viewer,"  HYPRE Euclid preconditioning\n");
604:     if (jac->eu_level != PETSC_DEFAULT) {
605:       PetscViewerASCIIPrintf(viewer,"    factorization levels %d\n",jac->eu_level);
606:     } else {
607:       PetscViewerASCIIPrintf(viewer,"    default factorization levels \n");
608:     }
609:     PetscViewerASCIIPrintf(viewer,"    drop tolerance %g\n",(double)jac->eu_droptolerance);
610:     PetscViewerASCIIPrintf(viewer,"    use Block-Jacobi? %D\n",jac->eu_bj);
611:   }
612:   return(0);
613: }

615: /* --------------------------------------------------------------------------------------------*/

617: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
618: {
619:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
620:   Mat_HYPRE          *hjac = (Mat_HYPRE*)(jac->hpmat->data);
621:   PetscErrorCode     ierr;
622:   HYPRE_ParCSRMatrix hmat;
623:   HYPRE_Complex      *xv,*bv;
624:   HYPRE_Complex      *sbv,*sxv;
625:   HYPRE_ParVector    jbv,jxv;
626:   PetscInt           hierr;

629:   PetscCitationsRegister(hypreCitation,&cite);
630:   VecSet(x,0.0);
631:   VecGetArrayRead(b,(const PetscScalar**)&bv);
632:   VecGetArray(x,(PetscScalar**)&xv);
633:   VecHYPRE_ParVectorReplacePointer(hjac->b,bv,sbv);
634:   VecHYPRE_ParVectorReplacePointer(hjac->x,xv,sxv);

636:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
637:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b,(void**)&jbv));
638:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x,(void**)&jxv));

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

645:   VecHYPRE_ParVectorReplacePointer(hjac->b,sbv,bv);
646:   VecHYPRE_ParVectorReplacePointer(hjac->x,sxv,xv);
647:   VecRestoreArray(x,(PetscScalar**)&xv);
648:   VecRestoreArrayRead(b,(const PetscScalar**)&bv);
649:   return(0);
650: }

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

655: static const char *HYPREBoomerAMGCycleType[]   = {"","V","W"};
656: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
657: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
658: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
659: static const char *HYPREBoomerAMGSmoothType[]  = {"Schwarz-smoothers","Pilut","ParaSails","Euclid"};
660: static const char *HYPREBoomerAMGRelaxType[]   = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi",
661:                                                   "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi",
662:                                                   "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination",
663:                                                   "" /* 10 */, "" /* 11 */, "" /* 12 */, "l1-Gauss-Seidel" /* nonsymmetric */, "backward-l1-Gauss-Seidel" /* nonsymmetric */,
664:                                                   "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"};
665: static const char *HYPREBoomerAMGInterpType[]  = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
666:                                                   "ext+i-cc", "standard", "standard-wts", "block", "block-wtd", "FF", "FF1",
667:                                                   "ext", "ad-wts", "ext-mm", "ext+i-mm", "ext+e-mm"};
668: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PetscOptionItems *PetscOptionsObject,PC pc)
669: {
670:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
672:   PetscInt       bs,n,indx,level;
673:   PetscBool      flg, tmp_truth;
674:   double         tmpdbl, twodbl[2];
675:   const char     *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};

678:   PetscOptionsHead(PetscOptionsObject,"HYPRE BoomerAMG Options");
679:   PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);
680:   if (flg) {
681:     jac->cycletype = indx+1;
682:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
683:   }
684:   PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
685:   if (flg) {
686:     if (jac->maxlevels < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
687:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
688:   }
689:   PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
690:   if (flg) {
691:     if (jac->maxiter < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
692:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
693:   }
694:   PetscOptionsReal("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);
695:   if (flg) {
696:     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);
697:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
698:   }
699:   bs = 1;
700:   if (pc->pmat) {
701:     MatGetBlockSize(pc->pmat,&bs);
702:   }
703:   PetscOptionsInt("-pc_hypre_boomeramg_numfunctions","Number of functions","HYPRE_BoomerAMGSetNumFunctions",bs,&bs,&flg);
704:   if (flg) {
705:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
706:   }

708:   PetscOptionsReal("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);
709:   if (flg) {
710:     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);
711:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
712:   }

714:   PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);
715:   if (flg) {
716:     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);
717:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
718:   }

720:   PetscOptionsInt("-pc_hypre_boomeramg_agg_nl","Number of levels of aggressive coarsening","None",jac->agg_nl,&jac->agg_nl,&flg);
721:   if (flg) {
722:     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);

724:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
725:   }

727:   PetscOptionsInt("-pc_hypre_boomeramg_agg_num_paths","Number of paths for aggressive coarsening","None",jac->agg_num_paths,&jac->agg_num_paths,&flg);
728:   if (flg) {
729:     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);

731:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
732:   }


735:   PetscOptionsReal("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
736:   if (flg) {
737:     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);
738:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
739:   }
740:   PetscOptionsReal("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
741:   if (flg) {
742:     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);
743:     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);
744:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
745:   }

747:   /* Grid sweeps */
748:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);
749:   if (flg) {
750:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx));
751:     /* modify the jac structure so we can view the updated options with PC_View */
752:     jac->gridsweeps[0] = indx;
753:     jac->gridsweeps[1] = indx;
754:     /*defaults coarse to 1 */
755:     jac->gridsweeps[2] = 1;
756:   }
757:   PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen","Use a nodal based coarsening 1-6","HYPRE_BoomerAMGSetNodal",jac->nodal_coarsening,&jac->nodal_coarsening,&flg);
758:   if (flg) {
759:     PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,jac->nodal_coarsening));
760:   }
761:   PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen_diag","Diagonal in strength matrix for nodal based coarsening 0-2","HYPRE_BoomerAMGSetNodalDiag",jac->nodal_coarsening_diag,&jac->nodal_coarsening_diag,&flg);
762:   if (flg) {
763:     PetscStackCallStandard(HYPRE_BoomerAMGSetNodalDiag,(jac->hsolver,jac->nodal_coarsening_diag));
764:   }
765:   PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_variant","Variant of algorithm 1-3","HYPRE_BoomerAMGSetInterpVecVariant",jac->vec_interp_variant, &jac->vec_interp_variant,&flg);
766:   if (flg) {
767:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecVariant,(jac->hsolver,jac->vec_interp_variant));
768:   }
769:   PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_qmax","Max elements per row for each Q","HYPRE_BoomerAMGSetInterpVecQMax",jac->vec_interp_qmax, &jac->vec_interp_qmax,&flg);
770:   if (flg) {
771:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecQMax,(jac->hsolver,jac->vec_interp_qmax));
772:   }
773:   PetscOptionsBool("-pc_hypre_boomeramg_vec_interp_smooth","Whether to smooth the interpolation vectors","HYPRE_BoomerAMGSetSmoothInterpVectors",jac->vec_interp_smooth, &jac->vec_interp_smooth,&flg);
774:   if (flg) {
775:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothInterpVectors,(jac->hsolver,jac->vec_interp_smooth));
776:   }
777:   PetscOptionsInt("-pc_hypre_boomeramg_interp_refine","Preprocess the interpolation matrix through iterative weight refinement","HYPRE_BoomerAMGSetInterpRefine",jac->interp_refine, &jac->interp_refine,&flg);
778:   if (flg) {
779:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpRefine,(jac->hsolver,jac->interp_refine));
780:   }
781:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);
782:   if (flg) {
783:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1));
784:     jac->gridsweeps[0] = indx;
785:   }
786:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);
787:   if (flg) {
788:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2));
789:     jac->gridsweeps[1] = indx;
790:   }
791:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);
792:   if (flg) {
793:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3));
794:     jac->gridsweeps[2] = indx;
795:   }

797:   /* Smooth type */
798:   PetscOptionsEList("-pc_hypre_boomeramg_smooth_type","Enable more complex smoothers","None",HYPREBoomerAMGSmoothType,ALEN(HYPREBoomerAMGSmoothType),HYPREBoomerAMGSmoothType[0],&indx,&flg);
799:   if (flg) {
800:     jac->smoothtype = indx;
801:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,indx+6));
802:     jac->smoothnumlevels = 25;
803:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,25));
804:   }

806:   /* Number of smoothing levels */
807:   PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_levels","Number of levels on which more complex smoothers are used","None",25,&indx,&flg);
808:   if (flg && (jac->smoothtype != -1)) {
809:     jac->smoothnumlevels = indx;
810:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,indx));
811:   }

813:   /* Number of levels for ILU(k) for Euclid */
814:   PetscOptionsInt("-pc_hypre_boomeramg_eu_level","Number of levels for ILU(k) in Euclid smoother","None",0,&indx,&flg);
815:   if (flg && (jac->smoothtype == 3)) {
816:     jac->eu_level = indx;
817:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,indx));
818:   }

820:   /* Filter for ILU(k) for Euclid */
821:   double droptolerance;
822:   PetscOptionsReal("-pc_hypre_boomeramg_eu_droptolerance","Drop tolerance for ILU(k) in Euclid smoother","None",0,&droptolerance,&flg);
823:   if (flg && (jac->smoothtype == 3)) {
824:     jac->eu_droptolerance = droptolerance;
825:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,droptolerance));
826:   }

828:   /* Use Block Jacobi ILUT for Euclid */
829:   PetscOptionsBool("-pc_hypre_boomeramg_eu_bj", "Use Block Jacobi for ILU in Euclid smoother?", "None", PETSC_FALSE, &tmp_truth, &flg);
830:   if (flg && (jac->smoothtype == 3)) {
831:     jac->eu_bj = tmp_truth;
832:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuBJ,(jac->hsolver,jac->eu_bj));
833:   }

835:   /* Relax type */
836:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
837:   if (flg) {
838:     jac->relaxtype[0] = jac->relaxtype[1]  = indx;
839:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx));
840:     /* by default, coarse type set to 9 */
841:     jac->relaxtype[2] = 9;
842:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, 9, 3));
843:   }
844:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
845:   if (flg) {
846:     jac->relaxtype[0] = indx;
847:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1));
848:   }
849:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
850:   if (flg) {
851:     jac->relaxtype[1] = indx;
852:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2));
853:   }
854:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);
855:   if (flg) {
856:     jac->relaxtype[2] = indx;
857:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3));
858:   }

860:   /* Relaxation Weight */
861:   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);
862:   if (flg) {
863:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl));
864:     jac->relaxweight = tmpdbl;
865:   }

867:   n         = 2;
868:   twodbl[0] = twodbl[1] = 1.0;
869:   PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
870:   if (flg) {
871:     if (n == 2) {
872:       indx =  (int)PetscAbsReal(twodbl[1]);
873:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx));
874:     } 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);
875:   }

877:   /* Outer relaxation Weight */
878:   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);
879:   if (flg) {
880:     PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl));
881:     jac->outerrelaxweight = tmpdbl;
882:   }

884:   n         = 2;
885:   twodbl[0] = twodbl[1] = 1.0;
886:   PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
887:   if (flg) {
888:     if (n == 2) {
889:       indx =  (int)PetscAbsReal(twodbl[1]);
890:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx));
891:     } 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);
892:   }

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

897:   if (flg && tmp_truth) {
898:     jac->relaxorder = 0;
899:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
900:   }
901:   PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);
902:   if (flg) {
903:     jac->measuretype = indx;
904:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
905:   }
906:   /* update list length 3/07 */
907:   PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);
908:   if (flg) {
909:     jac->coarsentype = indx;
910:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
911:   }

913:   PetscOptionsInt("-pc_hypre_boomeramg_max_coarse_size", "Maximum size of coarsest grid", "None", jac->maxc, &jac->maxc, &flg);
914:   if (flg) {
915:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxCoarseSize,(jac->hsolver, jac->maxc));
916:   }
917:   PetscOptionsInt("-pc_hypre_boomeramg_min_coarse_size", "Minimum size of coarsest grid", "None", jac->minc, &jac->minc, &flg);
918:   if (flg) {
919:     PetscStackCallStandard(HYPRE_BoomerAMGSetMinCoarseSize,(jac->hsolver, jac->minc));
920:   }

922:   /* AIR */
923: #if PETSC_PKG_HYPRE_VERSION_GE(2,18,0)
924:   PetscOptionsInt("-pc_hypre_boomeramg_restriction_type", "Type of AIR method (distance 1 or 2, 0 means no AIR)", "None", jac->Rtype, &jac->Rtype, NULL);
925:   PetscStackCallStandard(HYPRE_BoomerAMGSetRestriction,(jac->hsolver,jac->Rtype));
926:   if (jac->Rtype) {
927:     jac->interptype = 100; /* no way we can pass this with strings... Set it as default as in MFEM, then users can still customize it back to a different one */

929:     PetscOptionsReal("-pc_hypre_boomeramg_strongthresholdR","Threshold for R","None",jac->Rstrongthreshold,&jac->Rstrongthreshold,NULL);
930:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThresholdR,(jac->hsolver,jac->Rstrongthreshold));

932:     PetscOptionsReal("-pc_hypre_boomeramg_filterthresholdR","Filter threshold for R","None",jac->Rfilterthreshold,&jac->Rfilterthreshold,NULL);
933:     PetscStackCallStandard(HYPRE_BoomerAMGSetFilterThresholdR,(jac->hsolver,jac->Rfilterthreshold));

935:     PetscOptionsReal("-pc_hypre_boomeramg_Adroptol","Defines the drop tolerance for the A-matrices from the 2nd level of AMG","None",jac->Adroptol,&jac->Adroptol,NULL);
936:     PetscStackCallStandard(HYPRE_BoomerAMGSetADropTol,(jac->hsolver,jac->Adroptol));

938:     PetscOptionsInt("-pc_hypre_boomeramg_Adroptype","Drops the entries that are not on the diagonal and smaller than its row norm: type 1: 1-norm, 2: 2-norm, -1: infinity norm","None",jac->Adroptype,&jac->Adroptype,NULL);
939:     PetscStackCallStandard(HYPRE_BoomerAMGSetADropType,(jac->hsolver,jac->Adroptype));
940:   }
941: #endif

943:   /* new 3/07 */
944:   PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);
945:   if (flg || jac->Rtype) {
946:     if (flg) jac->interptype = indx;
947:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
948:   }

950:   PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
951:   if (flg) {
952:     level = 3;
953:     PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);

955:     jac->printstatistics = PETSC_TRUE;
956:     PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level));
957:   }

959:   PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);
960:   if (flg) {
961:     level = 3;
962:     PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);

964:     jac->printstatistics = PETSC_TRUE;
965:     PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level));
966:   }

968:   PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);
969:   if (flg && tmp_truth) {
970:     PetscInt tmp_int;
971:     PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);
972:     if (flg) jac->nodal_relax_levels = tmp_int;
973:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6));
974:     PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1));
975:     PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0));
976:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels));
977:   }

979:   PetscOptionsBool("-pc_hypre_boomeramg_keeptranspose", "Avoid transpose matvecs in preconditioner application", "None", jac->keeptranspose, &jac->keeptranspose, NULL);
980:   PetscStackCallStandard(HYPRE_BoomerAMGSetKeepTranspose,(jac->hsolver,jac->keeptranspose ? 1 : 0));

982:   /* options for ParaSails solvers */
983:   PetscOptionsEList("-pc_hypre_boomeramg_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flg);
984:   if (flg) {
985:     jac->symt = indx;
986:     PetscStackCallStandard(HYPRE_BoomerAMGSetSym,(jac->hsolver,jac->symt));
987:   }

989:   PetscOptionsTail();
990:   return(0);
991: }

993: 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)
994: {
995:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
997:   HYPRE_Int      oits;

1000:   PetscCitationsRegister(hypreCitation,&cite);
1001:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter));
1002:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol));
1003:   jac->applyrichardson = PETSC_TRUE;
1004:   PCApply_HYPRE(pc,b,y);
1005:   jac->applyrichardson = PETSC_FALSE;
1006:   PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,&oits));
1007:   *outits = oits;
1008:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
1009:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
1010:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
1011:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
1012:   return(0);
1013: }


1016: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
1017: {
1018:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1020:   PetscBool      iascii;

1023:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1024:   if (iascii) {
1025:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG preconditioning\n");
1026:     PetscViewerASCIIPrintf(viewer,"    Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);
1027:     PetscViewerASCIIPrintf(viewer,"    Maximum number of levels %D\n",jac->maxlevels);
1028:     PetscViewerASCIIPrintf(viewer,"    Maximum number of iterations PER hypre call %D\n",jac->maxiter);
1029:     PetscViewerASCIIPrintf(viewer,"    Convergence tolerance PER hypre call %g\n",(double)jac->tol);
1030:     PetscViewerASCIIPrintf(viewer,"    Threshold for strong coupling %g\n",(double)jac->strongthreshold);
1031:     PetscViewerASCIIPrintf(viewer,"    Interpolation truncation factor %g\n",(double)jac->truncfactor);
1032:     PetscViewerASCIIPrintf(viewer,"    Interpolation: max elements per row %D\n",jac->pmax);
1033:     if (jac->interp_refine) {
1034:       PetscViewerASCIIPrintf(viewer,"    Interpolation: number of steps of weighted refinement %D\n",jac->interp_refine);
1035:     }
1036:     PetscViewerASCIIPrintf(viewer,"    Number of levels of aggressive coarsening %D\n",jac->agg_nl);
1037:     PetscViewerASCIIPrintf(viewer,"    Number of paths for aggressive coarsening %D\n",jac->agg_num_paths);

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

1041:     PetscViewerASCIIPrintf(viewer,"    Sweeps down         %D\n",jac->gridsweeps[0]);
1042:     PetscViewerASCIIPrintf(viewer,"    Sweeps up           %D\n",jac->gridsweeps[1]);
1043:     PetscViewerASCIIPrintf(viewer,"    Sweeps on coarse    %D\n",jac->gridsweeps[2]);

1045:     PetscViewerASCIIPrintf(viewer,"    Relax down          %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[0]]);
1046:     PetscViewerASCIIPrintf(viewer,"    Relax up            %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[1]]);
1047:     PetscViewerASCIIPrintf(viewer,"    Relax on coarse     %s\n",HYPREBoomerAMGRelaxType[jac->relaxtype[2]]);

1049:     PetscViewerASCIIPrintf(viewer,"    Relax weight  (all)      %g\n",(double)jac->relaxweight);
1050:     PetscViewerASCIIPrintf(viewer,"    Outer relax weight (all) %g\n",(double)jac->outerrelaxweight);

1052:     if (jac->relaxorder) {
1053:       PetscViewerASCIIPrintf(viewer,"    Using CF-relaxation\n");
1054:     } else {
1055:       PetscViewerASCIIPrintf(viewer,"    Not using CF-relaxation\n");
1056:     }
1057:     if (jac->smoothtype!=-1) {
1058:       PetscViewerASCIIPrintf(viewer,"    Smooth type          %s\n",HYPREBoomerAMGSmoothType[jac->smoothtype]);
1059:       PetscViewerASCIIPrintf(viewer,"    Smooth num levels    %D\n",jac->smoothnumlevels);
1060:     } else {
1061:       PetscViewerASCIIPrintf(viewer,"    Not using more complex smoothers.\n");
1062:     }
1063:     if (jac->smoothtype==3) {
1064:       PetscViewerASCIIPrintf(viewer,"    Euclid ILU(k) levels %D\n",jac->eu_level);
1065:       PetscViewerASCIIPrintf(viewer,"    Euclid ILU(k) drop tolerance %g\n",(double)jac->eu_droptolerance);
1066:       PetscViewerASCIIPrintf(viewer,"    Euclid ILU use Block-Jacobi? %D\n",jac->eu_bj);
1067:     }
1068:     PetscViewerASCIIPrintf(viewer,"    Measure type        %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
1069:     PetscViewerASCIIPrintf(viewer,"    Coarsen type        %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
1070:     PetscViewerASCIIPrintf(viewer,"    Interpolation type  %s\n",jac->interptype != 100 ? HYPREBoomerAMGInterpType[jac->interptype] : "1pt");
1071:     if (jac->nodal_coarsening) {
1072:       PetscViewerASCIIPrintf(viewer,"    Using nodal coarsening (with HYPRE_BOOMERAMGSetNodal() %D\n",jac->nodal_coarsening);
1073:     }
1074:     if (jac->vec_interp_variant) {
1075:       PetscViewerASCIIPrintf(viewer,"    HYPRE_BoomerAMGSetInterpVecVariant() %D\n",jac->vec_interp_variant);
1076:       PetscViewerASCIIPrintf(viewer,"    HYPRE_BoomerAMGSetInterpVecQMax() %D\n",jac->vec_interp_qmax);
1077:       PetscViewerASCIIPrintf(viewer,"    HYPRE_BoomerAMGSetSmoothInterpVectors() %d\n",jac->vec_interp_smooth);
1078:     }
1079:     if (jac->nodal_relax) {
1080:       PetscViewerASCIIPrintf(viewer,"    Using nodal relaxation via Schwarz smoothing on levels %D\n",jac->nodal_relax_levels);
1081:     }

1083:     /* AIR */
1084:     if (jac->Rtype) {
1085:       PetscViewerASCIIPrintf(viewer,"    Using approximate ideal restriction type %D\n",jac->Rtype);
1086:       PetscViewerASCIIPrintf(viewer,"      Threshold for R %g\n",(double)jac->Rstrongthreshold);
1087:       PetscViewerASCIIPrintf(viewer,"      Filter for R %g\n",(double)jac->Rfilterthreshold);
1088:       PetscViewerASCIIPrintf(viewer,"      A drop tolerance %g\n",(double)jac->Adroptol);
1089:       PetscViewerASCIIPrintf(viewer,"      A drop type %D\n",jac->Adroptype);
1090:     }
1091:   }
1092:   return(0);
1093: }

1095: /* --------------------------------------------------------------------------------------------*/
1096: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PetscOptionItems *PetscOptionsObject,PC pc)
1097: {
1098:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1100:   PetscInt       indx;
1101:   PetscBool      flag;
1102:   const char     *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};

1105:   PetscOptionsHead(PetscOptionsObject,"HYPRE ParaSails Options");
1106:   PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);
1107:   PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshold,&jac->threshold,&flag);
1108:   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshold,jac->nlevels));

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

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

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

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

1122:   PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);
1123:   if (flag) {
1124:     jac->symt = indx;
1125:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
1126:   }

1128:   PetscOptionsTail();
1129:   return(0);
1130: }

1132: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
1133: {
1134:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1136:   PetscBool      iascii;
1137:   const char     *symt = 0;

1140:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1141:   if (iascii) {
1142:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails preconditioning\n");
1143:     PetscViewerASCIIPrintf(viewer,"    nlevels %d\n",jac->nlevels);
1144:     PetscViewerASCIIPrintf(viewer,"    threshold %g\n",(double)jac->threshold);
1145:     PetscViewerASCIIPrintf(viewer,"    filter %g\n",(double)jac->filter);
1146:     PetscViewerASCIIPrintf(viewer,"    load balance %g\n",(double)jac->loadbal);
1147:     PetscViewerASCIIPrintf(viewer,"    reuse nonzero structure %s\n",PetscBools[jac->ruse]);
1148:     PetscViewerASCIIPrintf(viewer,"    print info to screen %s\n",PetscBools[jac->logging]);
1149:     if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
1150:     else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
1151:     else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
1152:     else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
1153:     PetscViewerASCIIPrintf(viewer,"    %s\n",symt);
1154:   }
1155:   return(0);
1156: }
1157: /* --------------------------------------------------------------------------------------------*/
1158: static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PetscOptionItems *PetscOptionsObject,PC pc)
1159: {
1160:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1162:   PetscInt       n;
1163:   PetscBool      flag,flag2,flag3,flag4;

1166:   PetscOptionsHead(PetscOptionsObject,"HYPRE AMS Options");
1167:   PetscOptionsInt("-pc_hypre_ams_print_level","Debugging output level for AMS","None",jac->as_print,&jac->as_print,&flag);
1168:   if (flag) PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
1169:   PetscOptionsInt("-pc_hypre_ams_max_iter","Maximum number of AMS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
1170:   if (flag) PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1171:   PetscOptionsInt("-pc_hypre_ams_cycle_type","Cycle type for AMS multigrid","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag);
1172:   if (flag) PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1173:   PetscOptionsReal("-pc_hypre_ams_tol","Error tolerance for AMS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
1174:   if (flag) PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
1175:   PetscOptionsInt("-pc_hypre_ams_relax_type","Relaxation type for AMS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
1176:   PetscOptionsInt("-pc_hypre_ams_relax_times","Number of relaxation steps for AMS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
1177:   PetscOptionsReal("-pc_hypre_ams_relax_weight","Relaxation weight for AMS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
1178:   PetscOptionsReal("-pc_hypre_ams_omega","SSOR coefficient for AMS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
1179:   if (flag || flag2 || flag3 || flag4) {
1180:     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1181:                                                                       jac->as_relax_times,
1182:                                                                       jac->as_relax_weight,
1183:                                                                       jac->as_omega));
1184:   }
1185:   PetscOptionsReal("-pc_hypre_ams_amg_alpha_theta","Threshold for strong coupling of vector Poisson AMG solver","None",jac->as_amg_alpha_theta,&jac->as_amg_alpha_theta,&flag);
1186:   n = 5;
1187:   PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options","AMG options for vector Poisson","None",jac->as_amg_alpha_opts,&n,&flag2);
1188:   if (flag || flag2) {
1189:     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1190:                                                                      jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1191:                                                                      jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1192:                                                                      jac->as_amg_alpha_theta,
1193:                                                                      jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1194:                                                                      jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1195:   }
1196:   PetscOptionsReal("-pc_hypre_ams_amg_beta_theta","Threshold for strong coupling of scalar Poisson AMG solver","None",jac->as_amg_beta_theta,&jac->as_amg_beta_theta,&flag);
1197:   n = 5;
1198:   PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options","AMG options for scalar Poisson solver","None",jac->as_amg_beta_opts,&n,&flag2);
1199:   if (flag || flag2) {
1200:     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1201:                                                                     jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1202:                                                                     jac->as_amg_beta_opts[2],       /* AMG relax_type */
1203:                                                                     jac->as_amg_beta_theta,
1204:                                                                     jac->as_amg_beta_opts[3],       /* AMG interp_type */
1205:                                                                     jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1206:   }
1207:   PetscOptionsInt("-pc_hypre_ams_projection_frequency","Frequency at which a projection onto the compatible subspace for problems with zero conductivity regions is performed","None",jac->ams_proj_freq,&jac->ams_proj_freq,&flag);
1208:   if (flag) { /* override HYPRE's default only if the options is used */
1209:     PetscStackCallStandard(HYPRE_AMSSetProjectionFrequency,(jac->hsolver,jac->ams_proj_freq));
1210:   }
1211:   PetscOptionsTail();
1212:   return(0);
1213: }

1215: static PetscErrorCode PCView_HYPRE_AMS(PC pc,PetscViewer viewer)
1216: {
1217:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1219:   PetscBool      iascii;

1222:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1223:   if (iascii) {
1224:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS preconditioning\n");
1225:     PetscViewerASCIIPrintf(viewer,"    subspace iterations per application %d\n",jac->as_max_iter);
1226:     PetscViewerASCIIPrintf(viewer,"    subspace cycle type %d\n",jac->ams_cycle_type);
1227:     PetscViewerASCIIPrintf(viewer,"    subspace iteration tolerance %g\n",jac->as_tol);
1228:     PetscViewerASCIIPrintf(viewer,"    smoother type %d\n",jac->as_relax_type);
1229:     PetscViewerASCIIPrintf(viewer,"    number of smoothing steps %d\n",jac->as_relax_times);
1230:     PetscViewerASCIIPrintf(viewer,"    smoother weight %g\n",jac->as_relax_weight);
1231:     PetscViewerASCIIPrintf(viewer,"    smoother omega %g\n",jac->as_omega);
1232:     if (jac->alpha_Poisson) {
1233:       PetscViewerASCIIPrintf(viewer,"    vector Poisson solver (passed in by user)\n");
1234:     } else {
1235:       PetscViewerASCIIPrintf(viewer,"    vector Poisson solver (computed) \n");
1236:     }
1237:     PetscViewerASCIIPrintf(viewer,"        boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);
1238:     PetscViewerASCIIPrintf(viewer,"        boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
1239:     PetscViewerASCIIPrintf(viewer,"        boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);
1240:     PetscViewerASCIIPrintf(viewer,"        boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);
1241:     PetscViewerASCIIPrintf(viewer,"        boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
1242:     PetscViewerASCIIPrintf(viewer,"        boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);
1243:     if (!jac->ams_beta_is_zero) {
1244:       if (jac->beta_Poisson) {
1245:         PetscViewerASCIIPrintf(viewer,"    scalar Poisson solver (passed in by user)\n");
1246:       } else {
1247:         PetscViewerASCIIPrintf(viewer,"    scalar Poisson solver (computed) \n");
1248:       }
1249:       PetscViewerASCIIPrintf(viewer,"        boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);
1250:       PetscViewerASCIIPrintf(viewer,"        boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
1251:       PetscViewerASCIIPrintf(viewer,"        boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);
1252:       PetscViewerASCIIPrintf(viewer,"        boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);
1253:       PetscViewerASCIIPrintf(viewer,"        boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
1254:       PetscViewerASCIIPrintf(viewer,"        boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);
1255:       if (jac->ams_beta_is_zero_part) {
1256:         PetscViewerASCIIPrintf(viewer,"        compatible subspace projection frequency %d (-1 HYPRE uses default)\n",jac->ams_proj_freq);
1257:       }
1258:     } else {
1259:       PetscViewerASCIIPrintf(viewer,"    scalar Poisson solver not used (zero-conductivity everywhere) \n");
1260:     }
1261:   }
1262:   return(0);
1263: }

1265: static PetscErrorCode PCSetFromOptions_HYPRE_ADS(PetscOptionItems *PetscOptionsObject,PC pc)
1266: {
1267:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1269:   PetscInt       n;
1270:   PetscBool      flag,flag2,flag3,flag4;

1273:   PetscOptionsHead(PetscOptionsObject,"HYPRE ADS Options");
1274:   PetscOptionsInt("-pc_hypre_ads_print_level","Debugging output level for ADS","None",jac->as_print,&jac->as_print,&flag);
1275:   if (flag) PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
1276:   PetscOptionsInt("-pc_hypre_ads_max_iter","Maximum number of ADS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
1277:   if (flag) PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1278:   PetscOptionsInt("-pc_hypre_ads_cycle_type","Cycle type for ADS multigrid","None",jac->ads_cycle_type,&jac->ads_cycle_type,&flag);
1279:   if (flag) PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ads_cycle_type));
1280:   PetscOptionsReal("-pc_hypre_ads_tol","Error tolerance for ADS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
1281:   if (flag) PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
1282:   PetscOptionsInt("-pc_hypre_ads_relax_type","Relaxation type for ADS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
1283:   PetscOptionsInt("-pc_hypre_ads_relax_times","Number of relaxation steps for ADS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
1284:   PetscOptionsReal("-pc_hypre_ads_relax_weight","Relaxation weight for ADS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
1285:   PetscOptionsReal("-pc_hypre_ads_omega","SSOR coefficient for ADS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
1286:   if (flag || flag2 || flag3 || flag4) {
1287:     PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1288:                                                                       jac->as_relax_times,
1289:                                                                       jac->as_relax_weight,
1290:                                                                       jac->as_omega));
1291:   }
1292:   PetscOptionsReal("-pc_hypre_ads_ams_theta","Threshold for strong coupling of AMS solver inside ADS","None",jac->as_amg_alpha_theta,&jac->as_amg_alpha_theta,&flag);
1293:   n = 5;
1294:   PetscOptionsIntArray("-pc_hypre_ads_ams_options","AMG options for AMS solver inside ADS","None",jac->as_amg_alpha_opts,&n,&flag2);
1295:   PetscOptionsInt("-pc_hypre_ads_ams_cycle_type","Cycle type for AMS solver inside ADS","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag3);
1296:   if (flag || flag2 || flag3) {
1297:     PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type,             /* AMS cycle type */
1298:                                                                 jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1299:                                                                 jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1300:                                                                 jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1301:                                                                 jac->as_amg_alpha_theta,
1302:                                                                 jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1303:                                                                 jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1304:   }
1305:   PetscOptionsReal("-pc_hypre_ads_amg_theta","Threshold for strong coupling of vector AMG solver inside ADS","None",jac->as_amg_beta_theta,&jac->as_amg_beta_theta,&flag);
1306:   n = 5;
1307:   PetscOptionsIntArray("-pc_hypre_ads_amg_options","AMG options for vector AMG solver inside ADS","None",jac->as_amg_beta_opts,&n,&flag2);
1308:   if (flag || flag2) {
1309:     PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1310:                                                                 jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1311:                                                                 jac->as_amg_beta_opts[2],       /* AMG relax_type */
1312:                                                                 jac->as_amg_beta_theta,
1313:                                                                 jac->as_amg_beta_opts[3],       /* AMG interp_type */
1314:                                                                 jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1315:   }
1316:   PetscOptionsTail();
1317:   return(0);
1318: }

1320: static PetscErrorCode PCView_HYPRE_ADS(PC pc,PetscViewer viewer)
1321: {
1322:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1324:   PetscBool      iascii;

1327:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1328:   if (iascii) {
1329:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS preconditioning\n");
1330:     PetscViewerASCIIPrintf(viewer,"    subspace iterations per application %d\n",jac->as_max_iter);
1331:     PetscViewerASCIIPrintf(viewer,"    subspace cycle type %d\n",jac->ads_cycle_type);
1332:     PetscViewerASCIIPrintf(viewer,"    subspace iteration tolerance %g\n",jac->as_tol);
1333:     PetscViewerASCIIPrintf(viewer,"    smoother type %d\n",jac->as_relax_type);
1334:     PetscViewerASCIIPrintf(viewer,"    number of smoothing steps %d\n",jac->as_relax_times);
1335:     PetscViewerASCIIPrintf(viewer,"    smoother weight %g\n",jac->as_relax_weight);
1336:     PetscViewerASCIIPrintf(viewer,"    smoother omega %g\n",jac->as_omega);
1337:     PetscViewerASCIIPrintf(viewer,"    AMS solver using boomerAMG\n");
1338:     PetscViewerASCIIPrintf(viewer,"        subspace cycle type %d\n",jac->ams_cycle_type);
1339:     PetscViewerASCIIPrintf(viewer,"        coarsening type %d\n",jac->as_amg_alpha_opts[0]);
1340:     PetscViewerASCIIPrintf(viewer,"        levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
1341:     PetscViewerASCIIPrintf(viewer,"        relaxation type %d\n",jac->as_amg_alpha_opts[2]);
1342:     PetscViewerASCIIPrintf(viewer,"        interpolation type %d\n",jac->as_amg_alpha_opts[3]);
1343:     PetscViewerASCIIPrintf(viewer,"        max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
1344:     PetscViewerASCIIPrintf(viewer,"        strength threshold %g\n",jac->as_amg_alpha_theta);
1345:     PetscViewerASCIIPrintf(viewer,"    vector Poisson solver using boomerAMG\n");
1346:     PetscViewerASCIIPrintf(viewer,"        coarsening type %d\n",jac->as_amg_beta_opts[0]);
1347:     PetscViewerASCIIPrintf(viewer,"        levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
1348:     PetscViewerASCIIPrintf(viewer,"        relaxation type %d\n",jac->as_amg_beta_opts[2]);
1349:     PetscViewerASCIIPrintf(viewer,"        interpolation type %d\n",jac->as_amg_beta_opts[3]);
1350:     PetscViewerASCIIPrintf(viewer,"        max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
1351:     PetscViewerASCIIPrintf(viewer,"        strength threshold %g\n",jac->as_amg_beta_theta);
1352:   }
1353:   return(0);
1354: }

1356: static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE(PC pc, Mat G)
1357: {
1358:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1359:   PetscBool      ishypre;

1363:   PetscObjectTypeCompare((PetscObject)G,MATHYPRE,&ishypre);
1364:   if (ishypre) {
1365:     PetscObjectReference((PetscObject)G);
1366:     MatDestroy(&jac->G);
1367:     jac->G = G;
1368:   } else {
1369:     MatDestroy(&jac->G);
1370:     MatConvert(G,MATHYPRE,MAT_INITIAL_MATRIX,&jac->G);
1371:   }
1372:   return(0);
1373: }

1375: /*@
1376:  PCHYPRESetDiscreteGradient - Set discrete gradient matrix

1378:    Collective on PC

1380:    Input Parameters:
1381: +  pc - the preconditioning context
1382: -  G - the discrete gradient

1384:    Level: intermediate

1386:    Notes:
1387:     G should have as many rows as the number of edges and as many columns as the number of vertices in the mesh
1388:           Each row of G has 2 nonzeros, with column indexes being the global indexes of edge's endpoints: matrix entries are +1 and -1 depending on edge orientation

1390: .seealso:
1391: @*/
1392: PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G)
1393: {

1400:   PetscTryMethod(pc,"PCHYPRESetDiscreteGradient_C",(PC,Mat),(pc,G));
1401:   return(0);
1402: }

1404: static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C)
1405: {
1406:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1407:   PetscBool      ishypre;

1411:   PetscObjectTypeCompare((PetscObject)C,MATHYPRE,&ishypre);
1412:   if (ishypre) {
1413:     PetscObjectReference((PetscObject)C);
1414:     MatDestroy(&jac->C);
1415:     jac->C = C;
1416:   } else {
1417:     MatDestroy(&jac->C);
1418:     MatConvert(C,MATHYPRE,MAT_INITIAL_MATRIX,&jac->C);
1419:   }
1420:   return(0);
1421: }

1423: /*@
1424:  PCHYPRESetDiscreteCurl - Set discrete curl matrix

1426:    Collective on PC

1428:    Input Parameters:
1429: +  pc - the preconditioning context
1430: -  C - the discrete curl

1432:    Level: intermediate

1434:    Notes:
1435:     C should have as many rows as the number of faces and as many columns as the number of edges in the mesh
1436:           Each row of G has as many nonzeros as the number of edges of a face, with column indexes being the global indexes of the corresponding edge: matrix entries are +1 and -1 depending on edge orientation with respect to the face orientation

1438: .seealso:
1439: @*/
1440: PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C)
1441: {

1448:   PetscTryMethod(pc,"PCHYPRESetDiscreteCurl_C",(PC,Mat),(pc,C));
1449:   return(0);
1450: }

1452: static PetscErrorCode PCHYPRESetInterpolations_HYPRE(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[])
1453: {
1454:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1455:   PetscBool      ishypre;
1457:   PetscInt       i;

1460:   MatDestroy(&jac->RT_PiFull);
1461:   MatDestroy(&jac->ND_PiFull);
1462:   for (i=0;i<3;++i) {
1463:     MatDestroy(&jac->RT_Pi[i]);
1464:     MatDestroy(&jac->ND_Pi[i]);
1465:   }

1467:   jac->dim = dim;
1468:   if (RT_PiFull) {
1469:     PetscObjectTypeCompare((PetscObject)RT_PiFull,MATHYPRE,&ishypre);
1470:     if (ishypre) {
1471:       PetscObjectReference((PetscObject)RT_PiFull);
1472:       jac->RT_PiFull = RT_PiFull;
1473:     } else {
1474:       MatConvert(RT_PiFull,MATHYPRE,MAT_INITIAL_MATRIX,&jac->RT_PiFull);
1475:     }
1476:   }
1477:   if (RT_Pi) {
1478:     for (i=0;i<dim;++i) {
1479:       if (RT_Pi[i]) {
1480:         PetscObjectTypeCompare((PetscObject)RT_Pi[i],MATHYPRE,&ishypre);
1481:         if (ishypre) {
1482:           PetscObjectReference((PetscObject)RT_Pi[i]);
1483:           jac->RT_Pi[i] = RT_Pi[i];
1484:         } else {
1485:           MatConvert(RT_Pi[i],MATHYPRE,MAT_INITIAL_MATRIX,&jac->RT_Pi[i]);
1486:         }
1487:       }
1488:     }
1489:   }
1490:   if (ND_PiFull) {
1491:     PetscObjectTypeCompare((PetscObject)ND_PiFull,MATHYPRE,&ishypre);
1492:     if (ishypre) {
1493:       PetscObjectReference((PetscObject)ND_PiFull);
1494:       jac->ND_PiFull = ND_PiFull;
1495:     } else {
1496:       MatConvert(ND_PiFull,MATHYPRE,MAT_INITIAL_MATRIX,&jac->ND_PiFull);
1497:     }
1498:   }
1499:   if (ND_Pi) {
1500:     for (i=0;i<dim;++i) {
1501:       if (ND_Pi[i]) {
1502:         PetscObjectTypeCompare((PetscObject)ND_Pi[i],MATHYPRE,&ishypre);
1503:         if (ishypre) {
1504:           PetscObjectReference((PetscObject)ND_Pi[i]);
1505:           jac->ND_Pi[i] = ND_Pi[i];
1506:         } else {
1507:           MatConvert(ND_Pi[i],MATHYPRE,MAT_INITIAL_MATRIX,&jac->ND_Pi[i]);
1508:         }
1509:       }
1510:     }
1511:   }

1513:   return(0);
1514: }

1516: /*@
1517:  PCHYPRESetInterpolations - Set interpolation matrices for AMS/ADS preconditioner

1519:    Collective on PC

1521:    Input Parameters:
1522: +  pc - the preconditioning context
1523: -  dim - the dimension of the problem, only used in AMS
1524: -  RT_PiFull - Raviart-Thomas interpolation matrix
1525: -  RT_Pi - x/y/z component of Raviart-Thomas interpolation matrix
1526: -  ND_PiFull - Nedelec interpolation matrix
1527: -  ND_Pi - x/y/z component of Nedelec interpolation matrix

1529:    Notes:
1530:     For AMS, only Nedelec interpolation matrices are needed, the Raviart-Thomas interpolation matrices can be set to NULL.
1531:           For ADS, both type of interpolation matrices are needed.
1532:    Level: intermediate

1534: .seealso:
1535: @*/
1536: PetscErrorCode PCHYPRESetInterpolations(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[])
1537: {
1539:   PetscInt       i;

1543:   if (RT_PiFull) {
1546:   }
1547:   if (RT_Pi) {
1549:     for (i=0;i<dim;++i) {
1550:       if (RT_Pi[i]) {
1553:       }
1554:     }
1555:   }
1556:   if (ND_PiFull) {
1559:   }
1560:   if (ND_Pi) {
1562:     for (i=0;i<dim;++i) {
1563:       if (ND_Pi[i]) {
1566:       }
1567:     }
1568:   }
1569:   PetscTryMethod(pc,"PCHYPRESetInterpolations_C",(PC,PetscInt,Mat,Mat[],Mat,Mat[]),(pc,dim,RT_PiFull,RT_Pi,ND_PiFull,ND_Pi));
1570:   return(0);
1571: }

1573: static PetscErrorCode PCHYPRESetPoissonMatrix_HYPRE(PC pc, Mat A, PetscBool isalpha)
1574: {
1575:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1576:   PetscBool      ishypre;

1580:   PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
1581:   if (ishypre) {
1582:     if (isalpha) {
1583:       PetscObjectReference((PetscObject)A);
1584:       MatDestroy(&jac->alpha_Poisson);
1585:       jac->alpha_Poisson = A;
1586:     } else {
1587:       if (A) {
1588:         PetscObjectReference((PetscObject)A);
1589:       } else {
1590:         jac->ams_beta_is_zero = PETSC_TRUE;
1591:       }
1592:       MatDestroy(&jac->beta_Poisson);
1593:       jac->beta_Poisson = A;
1594:     }
1595:   } else {
1596:     if (isalpha) {
1597:       MatDestroy(&jac->alpha_Poisson);
1598:       MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&jac->alpha_Poisson);
1599:     } else {
1600:       if (A) {
1601:         MatDestroy(&jac->beta_Poisson);
1602:         MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&jac->beta_Poisson);
1603:       } else {
1604:         MatDestroy(&jac->beta_Poisson);
1605:         jac->ams_beta_is_zero = PETSC_TRUE;
1606:       }
1607:     }
1608:   }
1609:   return(0);
1610: }

1612: /*@
1613:  PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix

1615:    Collective on PC

1617:    Input Parameters:
1618: +  pc - the preconditioning context
1619: -  A - the matrix

1621:    Level: intermediate

1623:    Notes:
1624:     A should be obtained by discretizing the vector valued Poisson problem with linear finite elements

1626: .seealso:
1627: @*/
1628: PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A)
1629: {

1636:   PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_TRUE));
1637:   return(0);
1638: }

1640: /*@
1641:  PCHYPRESetBetaPoissonMatrix - Set Poisson matrix

1643:    Collective on PC

1645:    Input Parameters:
1646: +  pc - the preconditioning context
1647: -  A - the matrix

1649:    Level: intermediate

1651:    Notes:
1652:     A should be obtained by discretizing the Poisson problem with linear finite elements.
1653:           Following HYPRE convention, the scalar Poisson solver of AMS can be turned off by passing NULL.

1655: .seealso:
1656: @*/
1657: PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A)
1658: {

1663:   if (A) {
1666:   }
1667:   PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_FALSE));
1668:   return(0);
1669: }

1671: static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE(PC pc,Vec ozz, Vec zoz, Vec zzo)
1672: {
1673:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1674:   PetscErrorCode     ierr;

1677:   /* throw away any vector if already set */
1678:   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0]));
1679:   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1]));
1680:   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2]));
1681:   jac->constants[0] = NULL;
1682:   jac->constants[1] = NULL;
1683:   jac->constants[2] = NULL;
1684:   VecHYPRE_IJVectorCreate(ozz,&jac->constants[0]);
1685:   VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);
1686:   VecHYPRE_IJVectorCreate(zoz,&jac->constants[1]);
1687:   VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);
1688:   jac->dim = 2;
1689:   if (zzo) {
1690:     VecHYPRE_IJVectorCreate(zzo,&jac->constants[2]);
1691:     VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);
1692:     jac->dim++;
1693:   }
1694:   return(0);
1695: }

1697: /*@
1698:  PCHYPRESetEdgeConstantVectors - Set the representation of the constant vector fields in edge element basis

1700:    Collective on PC

1702:    Input Parameters:
1703: +  pc - the preconditioning context
1704: -  ozz - vector representing (1,0,0) (or (1,0) in 2D)
1705: -  zoz - vector representing (0,1,0) (or (0,1) in 2D)
1706: -  zzo - vector representing (0,0,1) (use NULL in 2D)

1708:    Level: intermediate

1710:    Notes:

1712: .seealso:
1713: @*/
1714: PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo)
1715: {

1726:   PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));
1727:   return(0);
1728: }

1730: static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1731: {
1732:   PC_HYPRE        *jac = (PC_HYPRE*)pc->data;
1733:   Vec             tv;
1734:   PetscInt        i;
1735:   PetscErrorCode  ierr;

1738:   /* throw away any coordinate vector if already set */
1739:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0]));
1740:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1]));
1741:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2]));
1742:   jac->dim = dim;

1744:   /* compute IJ vector for coordinates */
1745:   VecCreate(PetscObjectComm((PetscObject)pc),&tv);
1746:   VecSetType(tv,VECSTANDARD);
1747:   VecSetSizes(tv,nloc,PETSC_DECIDE);
1748:   for (i=0;i<dim;i++) {
1749:     PetscScalar *array;
1750:     PetscInt    j;

1752:     VecHYPRE_IJVectorCreate(tv,&jac->coords[i]);
1753:     VecGetArrayWrite(tv,&array);
1754:     for (j=0;j<nloc;j++) {
1755:       array[j] = coords[j*dim+i];
1756:     }
1757:     PetscStackCallStandard(HYPRE_IJVectorSetValues,(jac->coords[i],nloc,NULL,(HYPRE_Complex*)array));
1758:     PetscStackCallStandard(HYPRE_IJVectorAssemble,(jac->coords[i]));
1759:     VecRestoreArrayWrite(tv,&array);
1760:   }
1761:   VecDestroy(&tv);
1762:   return(0);
1763: }

1765: /* ---------------------------------------------------------------------------------*/

1767: static PetscErrorCode  PCHYPREGetType_HYPRE(PC pc,const char *name[])
1768: {
1769:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

1772:   *name = jac->hypre_type;
1773:   return(0);
1774: }

1776: static PetscErrorCode  PCHYPRESetType_HYPRE(PC pc,const char name[])
1777: {
1778:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1780:   PetscBool      flag;

1783:   if (jac->hypre_type) {
1784:     PetscStrcmp(jac->hypre_type,name,&flag);
1785:     if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set");
1786:     return(0);
1787:   } else {
1788:     PetscStrallocpy(name, &jac->hypre_type);
1789:   }

1791:   jac->maxiter         = PETSC_DEFAULT;
1792:   jac->tol             = PETSC_DEFAULT;
1793:   jac->printstatistics = PetscLogPrintInfo;

1795:   PetscStrcmp("pilut",jac->hypre_type,&flag);
1796:   if (flag) {
1797:     MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1798:     PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver));
1799:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
1800:     pc->ops->view           = PCView_HYPRE_Pilut;
1801:     jac->destroy            = HYPRE_ParCSRPilutDestroy;
1802:     jac->setup              = HYPRE_ParCSRPilutSetup;
1803:     jac->solve              = HYPRE_ParCSRPilutSolve;
1804:     jac->factorrowsize      = PETSC_DEFAULT;
1805:     return(0);
1806:   }
1807:   PetscStrcmp("euclid",jac->hypre_type,&flag);
1808:   if (flag) {
1809: #if defined(PETSC_HAVE_64BIT_INDICES)
1810:     SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_SUP,"Hypre Euclid not support with 64 bit indices");
1811: #endif
1812:     MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1813:     PetscStackCallStandard(HYPRE_EuclidCreate,(jac->comm_hypre,&jac->hsolver));
1814:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
1815:     pc->ops->view           = PCView_HYPRE_Euclid;
1816:     jac->destroy            = HYPRE_EuclidDestroy;
1817:     jac->setup              = HYPRE_EuclidSetup;
1818:     jac->solve              = HYPRE_EuclidSolve;
1819:     jac->factorrowsize      = PETSC_DEFAULT;
1820:     jac->eu_level           = PETSC_DEFAULT; /* default */
1821:     return(0);
1822:   }
1823:   PetscStrcmp("parasails",jac->hypre_type,&flag);
1824:   if (flag) {
1825:     MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1826:     PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver));
1827:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
1828:     pc->ops->view           = PCView_HYPRE_ParaSails;
1829:     jac->destroy            = HYPRE_ParaSailsDestroy;
1830:     jac->setup              = HYPRE_ParaSailsSetup;
1831:     jac->solve              = HYPRE_ParaSailsSolve;
1832:     /* initialize */
1833:     jac->nlevels   = 1;
1834:     jac->threshold = .1;
1835:     jac->filter    = .1;
1836:     jac->loadbal   = 0;
1837:     if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE;
1838:     else jac->logging = (int) PETSC_FALSE;

1840:     jac->ruse = (int) PETSC_FALSE;
1841:     jac->symt = 0;
1842:     PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshold,jac->nlevels));
1843:     PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
1844:     PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
1845:     PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
1846:     PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
1847:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
1848:     return(0);
1849:   }
1850:   PetscStrcmp("boomeramg",jac->hypre_type,&flag);
1851:   if (flag) {
1852:     HYPRE_BoomerAMGCreate(&jac->hsolver);
1853:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_BoomerAMG;
1854:     pc->ops->view            = PCView_HYPRE_BoomerAMG;
1855:     pc->ops->applytranspose  = PCApplyTranspose_HYPRE_BoomerAMG;
1856:     pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
1857:     PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",PCGetInterpolations_BoomerAMG);
1858:     PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",PCGetCoarseOperators_BoomerAMG);
1859:     jac->destroy             = HYPRE_BoomerAMGDestroy;
1860:     jac->setup               = HYPRE_BoomerAMGSetup;
1861:     jac->solve               = HYPRE_BoomerAMGSolve;
1862:     jac->applyrichardson     = PETSC_FALSE;
1863:     /* these defaults match the hypre defaults */
1864:     jac->cycletype        = 1;
1865:     jac->maxlevels        = 25;
1866:     jac->maxiter          = 1;
1867:     jac->tol              = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
1868:     jac->truncfactor      = 0.0;
1869:     jac->strongthreshold  = .25;
1870:     jac->maxrowsum        = .9;
1871:     jac->coarsentype      = 6;
1872:     jac->measuretype      = 0;
1873:     jac->gridsweeps[0]    = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
1874:     jac->smoothtype       = -1; /* Not set by default */
1875:     jac->smoothnumlevels  = 25;
1876:     jac->eu_level         = 0;
1877:     jac->eu_droptolerance = 0;
1878:     jac->eu_bj            = 0;
1879:     jac->relaxtype[0]     = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a PC - most likely with CG */
1880:     jac->relaxtype[2]     = 9; /*G.E. */
1881:     jac->relaxweight      = 1.0;
1882:     jac->outerrelaxweight = 1.0;
1883:     jac->relaxorder       = 1;
1884:     jac->interptype       = 0;
1885:     jac->Rtype            = 0;
1886:     jac->Rstrongthreshold = 0.25;
1887:     jac->Rfilterthreshold = 0.0;
1888:     jac->Adroptype        = -1;
1889:     jac->Adroptol         = 0.0;
1890:     jac->agg_nl           = 0;
1891:     jac->pmax             = 0;
1892:     jac->truncfactor      = 0.0;
1893:     jac->agg_num_paths    = 1;
1894:     jac->maxc             = 9;
1895:     jac->minc             = 1;

1897:     jac->nodal_coarsening      = 0;
1898:     jac->nodal_coarsening_diag = 0;
1899:     jac->vec_interp_variant    = 0;
1900:     jac->vec_interp_qmax       = 0;
1901:     jac->vec_interp_smooth     = PETSC_FALSE;
1902:     jac->interp_refine         = 0;
1903:     jac->nodal_relax           = PETSC_FALSE;
1904:     jac->nodal_relax_levels    = 1;
1905:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
1906:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
1907:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
1908:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
1909:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
1910:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
1911:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
1912:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
1913:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
1914:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
1915:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
1916:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
1917:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
1918:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
1919:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0]));  /*defaults coarse to 9*/
1920:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */
1921:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxCoarseSize,(jac->hsolver, jac->maxc));
1922:     PetscStackCallStandard(HYPRE_BoomerAMGSetMinCoarseSize,(jac->hsolver, jac->minc));

1924:     /* AIR */
1925:     PetscStackCallStandard(HYPRE_BoomerAMGSetRestriction,(jac->hsolver,jac->Rtype));
1926:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThresholdR,(jac->hsolver,jac->Rstrongthreshold));
1927:     PetscStackCallStandard(HYPRE_BoomerAMGSetFilterThresholdR,(jac->hsolver,jac->Rfilterthreshold));
1928:     PetscStackCallStandard(HYPRE_BoomerAMGSetADropTol,(jac->hsolver,jac->Adroptol));
1929:     PetscStackCallStandard(HYPRE_BoomerAMGSetADropType,(jac->hsolver,jac->Adroptype));
1930:     return(0);
1931:   }
1932:   PetscStrcmp("ams",jac->hypre_type,&flag);
1933:   if (flag) {
1934:     HYPRE_AMSCreate(&jac->hsolver);
1935:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_AMS;
1936:     pc->ops->view            = PCView_HYPRE_AMS;
1937:     jac->destroy             = HYPRE_AMSDestroy;
1938:     jac->setup               = HYPRE_AMSSetup;
1939:     jac->solve               = HYPRE_AMSSolve;
1940:     jac->coords[0]           = NULL;
1941:     jac->coords[1]           = NULL;
1942:     jac->coords[2]           = NULL;
1943:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */
1944:     jac->as_print           = 0;
1945:     jac->as_max_iter        = 1; /* used as a preconditioner */
1946:     jac->as_tol             = 0.; /* used as a preconditioner */
1947:     jac->ams_cycle_type     = 13;
1948:     /* Smoothing options */
1949:     jac->as_relax_type      = 2;
1950:     jac->as_relax_times     = 1;
1951:     jac->as_relax_weight    = 1.0;
1952:     jac->as_omega           = 1.0;
1953:     /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1954:     jac->as_amg_alpha_opts[0] = 10;
1955:     jac->as_amg_alpha_opts[1] = 1;
1956:     jac->as_amg_alpha_opts[2] = 6;
1957:     jac->as_amg_alpha_opts[3] = 6;
1958:     jac->as_amg_alpha_opts[4] = 4;
1959:     jac->as_amg_alpha_theta   = 0.25;
1960:     /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1961:     jac->as_amg_beta_opts[0] = 10;
1962:     jac->as_amg_beta_opts[1] = 1;
1963:     jac->as_amg_beta_opts[2] = 6;
1964:     jac->as_amg_beta_opts[3] = 6;
1965:     jac->as_amg_beta_opts[4] = 4;
1966:     jac->as_amg_beta_theta   = 0.25;
1967:     PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
1968:     PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1969:     PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1970:     PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
1971:     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1972:                                                                       jac->as_relax_times,
1973:                                                                       jac->as_relax_weight,
1974:                                                                       jac->as_omega));
1975:     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1976:                                                                      jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1977:                                                                      jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1978:                                                                      jac->as_amg_alpha_theta,
1979:                                                                      jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1980:                                                                      jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1981:     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1982:                                                                     jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1983:                                                                     jac->as_amg_beta_opts[2],       /* AMG relax_type */
1984:                                                                     jac->as_amg_beta_theta,
1985:                                                                     jac->as_amg_beta_opts[3],       /* AMG interp_type */
1986:                                                                     jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1987:     /* Zero conductivity */
1988:     jac->ams_beta_is_zero      = PETSC_FALSE;
1989:     jac->ams_beta_is_zero_part = PETSC_FALSE;
1990:     return(0);
1991:   }
1992:   PetscStrcmp("ads",jac->hypre_type,&flag);
1993:   if (flag) {
1994:     HYPRE_ADSCreate(&jac->hsolver);
1995:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_ADS;
1996:     pc->ops->view            = PCView_HYPRE_ADS;
1997:     jac->destroy             = HYPRE_ADSDestroy;
1998:     jac->setup               = HYPRE_ADSSetup;
1999:     jac->solve               = HYPRE_ADSSolve;
2000:     jac->coords[0]           = NULL;
2001:     jac->coords[1]           = NULL;
2002:     jac->coords[2]           = NULL;
2003:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */
2004:     jac->as_print           = 0;
2005:     jac->as_max_iter        = 1; /* used as a preconditioner */
2006:     jac->as_tol             = 0.; /* used as a preconditioner */
2007:     jac->ads_cycle_type     = 13;
2008:     /* Smoothing options */
2009:     jac->as_relax_type      = 2;
2010:     jac->as_relax_times     = 1;
2011:     jac->as_relax_weight    = 1.0;
2012:     jac->as_omega           = 1.0;
2013:     /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */
2014:     jac->ams_cycle_type       = 14;
2015:     jac->as_amg_alpha_opts[0] = 10;
2016:     jac->as_amg_alpha_opts[1] = 1;
2017:     jac->as_amg_alpha_opts[2] = 6;
2018:     jac->as_amg_alpha_opts[3] = 6;
2019:     jac->as_amg_alpha_opts[4] = 4;
2020:     jac->as_amg_alpha_theta   = 0.25;
2021:     /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
2022:     jac->as_amg_beta_opts[0] = 10;
2023:     jac->as_amg_beta_opts[1] = 1;
2024:     jac->as_amg_beta_opts[2] = 6;
2025:     jac->as_amg_beta_opts[3] = 6;
2026:     jac->as_amg_beta_opts[4] = 4;
2027:     jac->as_amg_beta_theta   = 0.25;
2028:     PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
2029:     PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
2030:     PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
2031:     PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
2032:     PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
2033:                                                                       jac->as_relax_times,
2034:                                                                       jac->as_relax_weight,
2035:                                                                       jac->as_omega));
2036:     PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type,             /* AMG coarsen type */
2037:                                                                 jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
2038:                                                                 jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
2039:                                                                 jac->as_amg_alpha_opts[2],       /* AMG relax_type */
2040:                                                                 jac->as_amg_alpha_theta,
2041:                                                                 jac->as_amg_alpha_opts[3],       /* AMG interp_type */
2042:                                                                 jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
2043:     PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
2044:                                                                 jac->as_amg_beta_opts[1],       /* AMG agg_levels */
2045:                                                                 jac->as_amg_beta_opts[2],       /* AMG relax_type */
2046:                                                                 jac->as_amg_beta_theta,
2047:                                                                 jac->as_amg_beta_opts[3],       /* AMG interp_type */
2048:                                                                 jac->as_amg_beta_opts[4]));     /* AMG Pmax */
2049:     return(0);
2050:   }
2051:   PetscFree(jac->hypre_type);

2053:   jac->hypre_type = NULL;
2054:   SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_UNKNOWN_TYPE,"Unknown HYPRE preconditioner %s; Choices are euclid, pilut, parasails, boomeramg, ams",name);
2055:   return(0);
2056: }

2058: /*
2059:     It only gets here if the HYPRE type has not been set before the call to
2060:    ...SetFromOptions() which actually is most of the time
2061: */
2062: PetscErrorCode PCSetFromOptions_HYPRE(PetscOptionItems *PetscOptionsObject,PC pc)
2063: {
2065:   PetscInt       indx;
2066:   const char     *type[] = {"euclid","pilut","parasails","boomeramg","ams","ads"};
2067:   PetscBool      flg;

2070:   PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");
2071:   PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,ALEN(type),"boomeramg",&indx,&flg);
2072:   if (flg) {
2073:     PCHYPRESetType_HYPRE(pc,type[indx]);
2074:   } else {
2075:     PCHYPRESetType_HYPRE(pc,"boomeramg");
2076:   }
2077:   if (pc->ops->setfromoptions) {
2078:     pc->ops->setfromoptions(PetscOptionsObject,pc);
2079:   }
2080:   PetscOptionsTail();
2081:   return(0);
2082: }

2084: /*@C
2085:      PCHYPRESetType - Sets which hypre preconditioner you wish to use

2087:    Input Parameters:
2088: +     pc - the preconditioner context
2089: -     name - either  euclid, pilut, parasails, boomeramg, ams, ads

2091:    Options Database Keys:
2092:    -pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads

2094:    Level: intermediate

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

2099: @*/
2100: PetscErrorCode  PCHYPRESetType(PC pc,const char name[])
2101: {

2107:   PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));
2108:   return(0);
2109: }

2111: /*@C
2112:      PCHYPREGetType - Gets which hypre preconditioner you are using

2114:    Input Parameter:
2115: .     pc - the preconditioner context

2117:    Output Parameter:
2118: .     name - either  euclid, pilut, parasails, boomeramg, ams, ads

2120:    Level: intermediate

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

2125: @*/
2126: PetscErrorCode  PCHYPREGetType(PC pc,const char *name[])
2127: {

2133:   PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));
2134:   return(0);
2135: }

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

2140:    Options Database Keys:
2141: +   -pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads
2142: -   Too many others to list, run with -pc_type hypre -pc_hypre_type XXX -help to see options for the XXX
2143:           preconditioner

2145:    Level: intermediate

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

2152:           The options -pc_hypre_boomeramg_max_iter and -pc_hypre_boomeramg_tol refer to the number of iterations
2153:           (V-cycles) and tolerance that boomeramg does EACH time it is called. So for example, if
2154:           -pc_hypre_boomeramg_max_iter is set to 2 then 2-V-cycles are being used to define the preconditioner
2155:           (-pc_hypre_boomeramg_tol should be set to 0.0 - the default - to strictly use a fixed number of
2156:           iterations per hypre call). -ksp_max_it and -ksp_rtol STILL determine the total number of iterations
2157:           and tolerance for the Krylov solver. For example, if -pc_hypre_boomeramg_max_iter is 2 and -ksp_max_it is 10
2158:           then AT MOST twenty V-cycles of boomeramg will be called.

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

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

2170:           MatSetNearNullSpace() - if you provide a near null space to your matrix it is ignored by hypre UNLESS you also use
2171:           the following two options:

2173:    Options Database Keys:
2174: +   -pc_hypre_boomeramg_nodal_coarsen <n> - where n is from 1 to 6 (see HYPRE_BOOMERAMGSetNodal())
2175: -   -pc_hypre_boomeramg_vec_interp_variant <v> - where v is from 1 to 3 (see HYPRE_BoomerAMGSetInterpVecVariant())

2177:           Depending on the linear system you may see the same or different convergence depending on the values you use.

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

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

2184: M*/

2186: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
2187: {
2188:   PC_HYPRE       *jac;

2192:   PetscNewLog(pc,&jac);

2194:   pc->data                = jac;
2195:   pc->ops->reset          = PCReset_HYPRE;
2196:   pc->ops->destroy        = PCDestroy_HYPRE;
2197:   pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
2198:   pc->ops->setup          = PCSetUp_HYPRE;
2199:   pc->ops->apply          = PCApply_HYPRE;
2200:   jac->comm_hypre         = MPI_COMM_NULL;
2201:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);
2202:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);
2203:   PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);
2204:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);
2205:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",PCHYPRESetDiscreteCurl_HYPRE);
2206:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",PCHYPRESetInterpolations_HYPRE);
2207:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE);
2208:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",PCHYPRESetPoissonMatrix_HYPRE);
2209:   return(0);
2210: }

2212: /* ---------------------------------------------------------------------------------------------------------------------------------*/

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

2218:   /* keep copy of PFMG options used so may view them */
2219:   PetscInt its;
2220:   double   tol;
2221:   PetscInt relax_type;
2222:   PetscInt rap_type;
2223:   PetscInt num_pre_relax,num_post_relax;
2224:   PetscInt max_levels;
2225: } PC_PFMG;

2227: PetscErrorCode PCDestroy_PFMG(PC pc)
2228: {
2230:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

2233:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
2234:   MPI_Comm_free(&ex->hcomm);
2235:   PetscFree(pc->data);
2236:   return(0);
2237: }

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

2242: PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
2243: {
2245:   PetscBool      iascii;
2246:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

2249:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2250:   if (iascii) {
2251:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG preconditioning\n");
2252:     PetscViewerASCIIPrintf(viewer,"    max iterations %d\n",ex->its);
2253:     PetscViewerASCIIPrintf(viewer,"    tolerance %g\n",ex->tol);
2254:     PetscViewerASCIIPrintf(viewer,"    relax type %s\n",PFMGRelaxType[ex->relax_type]);
2255:     PetscViewerASCIIPrintf(viewer,"    RAP type %s\n",PFMGRAPType[ex->rap_type]);
2256:     PetscViewerASCIIPrintf(viewer,"    number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
2257:     PetscViewerASCIIPrintf(viewer,"    max levels %d\n",ex->max_levels);
2258:   }
2259:   return(0);
2260: }

2262: PetscErrorCode PCSetFromOptions_PFMG(PetscOptionItems *PetscOptionsObject,PC pc)
2263: {
2265:   PC_PFMG        *ex = (PC_PFMG*) pc->data;
2266:   PetscBool      flg = PETSC_FALSE;

2269:   PetscOptionsHead(PetscOptionsObject,"PFMG options");
2270:   PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);
2271:   if (flg) {
2272:     PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,3));
2273:   }
2274:   PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);
2275:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its));
2276:   PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
2277:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax));
2278:   PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
2279:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax));

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

2284:   PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);
2285:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol));
2286:   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);
2287:   PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type));
2288:   PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);
2289:   PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type));
2290:   PetscOptionsTail();
2291:   return(0);
2292: }

2294: PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
2295: {
2296:   PetscErrorCode    ierr;
2297:   PC_PFMG           *ex = (PC_PFMG*) pc->data;
2298:   PetscScalar       *yy;
2299:   const PetscScalar *xx;
2300:   PetscInt          ilower[3],iupper[3];
2301:   HYPRE_Int         hlower[3],hupper[3];
2302:   Mat_HYPREStruct   *mx = (Mat_HYPREStruct*)(pc->pmat->data);

2305:   PetscCitationsRegister(hypreCitation,&cite);
2306:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
2307:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2308:   iupper[0] += ilower[0] - 1;
2309:   iupper[1] += ilower[1] - 1;
2310:   iupper[2] += ilower[2] - 1;
2311:   hlower[0]  = (HYPRE_Int)ilower[0];
2312:   hlower[1]  = (HYPRE_Int)ilower[1];
2313:   hlower[2]  = (HYPRE_Int)ilower[2];
2314:   hupper[0]  = (HYPRE_Int)iupper[0];
2315:   hupper[1]  = (HYPRE_Int)iupper[1];
2316:   hupper[2]  = (HYPRE_Int)iupper[2];

2318:   /* copy x values over to hypre */
2319:   PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
2320:   VecGetArrayRead(x,&xx);
2321:   PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,hlower,hupper,(HYPRE_Complex*)xx));
2322:   VecRestoreArrayRead(x,&xx);
2323:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
2324:   PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx));

2326:   /* copy solution values back to PETSc */
2327:   VecGetArray(y,&yy);
2328:   PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,hlower,hupper,(HYPRE_Complex*)yy));
2329:   VecRestoreArray(y,&yy);
2330:   return(0);
2331: }

2333: 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)
2334: {
2335:   PC_PFMG        *jac = (PC_PFMG*)pc->data;
2337:   HYPRE_Int      oits;

2340:   PetscCitationsRegister(hypreCitation,&cite);
2341:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its));
2342:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol));

2344:   PCApply_PFMG(pc,b,y);
2345:   PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,&oits));
2346:   *outits = oits;
2347:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2348:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
2349:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol));
2350:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its));
2351:   return(0);
2352: }


2355: PetscErrorCode PCSetUp_PFMG(PC pc)
2356: {
2357:   PetscErrorCode  ierr;
2358:   PC_PFMG         *ex = (PC_PFMG*) pc->data;
2359:   Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
2360:   PetscBool       flg;

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

2366:   /* create the hypre solver object and set its information */
2367:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
2368:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
2369:   PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
2370:   PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver));
2371:   return(0);
2372: }

2374: /*MC
2375:      PCPFMG - the hypre PFMG multigrid solver

2377:    Level: advanced

2379:    Options Database:
2380: + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner
2381: . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2382: . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2383: . -pc_pfmg_tol <tol> tolerance of PFMG
2384: . -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
2385: - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin

2387:    Notes:
2388:     This is for CELL-centered descretizations

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

2393: .seealso:  PCMG, MATHYPRESTRUCT
2394: M*/

2396: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
2397: {
2399:   PC_PFMG        *ex;

2402:   PetscNew(&ex); \
2403:   pc->data = ex;

2405:   ex->its            = 1;
2406:   ex->tol            = 1.e-8;
2407:   ex->relax_type     = 1;
2408:   ex->rap_type       = 0;
2409:   ex->num_pre_relax  = 1;
2410:   ex->num_post_relax = 1;
2411:   ex->max_levels     = 0;

2413:   pc->ops->setfromoptions  = PCSetFromOptions_PFMG;
2414:   pc->ops->view            = PCView_PFMG;
2415:   pc->ops->destroy         = PCDestroy_PFMG;
2416:   pc->ops->apply           = PCApply_PFMG;
2417:   pc->ops->applyrichardson = PCApplyRichardson_PFMG;
2418:   pc->ops->setup           = PCSetUp_PFMG;

2420:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
2421:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
2422:   return(0);
2423: }

2425: /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/

2427: /* we know we are working with a HYPRE_SStructMatrix */
2428: typedef struct {
2429:   MPI_Comm            hcomm;       /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
2430:   HYPRE_SStructSolver ss_solver;

2432:   /* keep copy of SYSPFMG options used so may view them */
2433:   PetscInt its;
2434:   double   tol;
2435:   PetscInt relax_type;
2436:   PetscInt num_pre_relax,num_post_relax;
2437: } PC_SysPFMG;

2439: PetscErrorCode PCDestroy_SysPFMG(PC pc)
2440: {
2442:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

2445:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2446:   MPI_Comm_free(&ex->hcomm);
2447:   PetscFree(pc->data);
2448:   return(0);
2449: }

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

2453: PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
2454: {
2456:   PetscBool      iascii;
2457:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

2460:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2461:   if (iascii) {
2462:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG preconditioning\n");
2463:     PetscViewerASCIIPrintf(viewer,"  max iterations %d\n",ex->its);
2464:     PetscViewerASCIIPrintf(viewer,"  tolerance %g\n",ex->tol);
2465:     PetscViewerASCIIPrintf(viewer,"  relax type %s\n",PFMGRelaxType[ex->relax_type]);
2466:     PetscViewerASCIIPrintf(viewer,"  number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
2467:   }
2468:   return(0);
2469: }

2471: PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptionItems *PetscOptionsObject,PC pc)
2472: {
2474:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;
2475:   PetscBool      flg = PETSC_FALSE;

2478:   PetscOptionsHead(PetscOptionsObject,"SysPFMG options");
2479:   PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);
2480:   if (flg) {
2481:     PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,3));
2482:   }
2483:   PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);
2484:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its));
2485:   PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
2486:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax));
2487:   PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
2488:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax));

2490:   PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);
2491:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol));
2492:   PetscOptionsEList("-pc_syspfmg_relax_type","Relax type for the up and down cycles","HYPRE_SStructSysPFMGSetRelaxType",SysPFMGRelaxType,ALEN(SysPFMGRelaxType),SysPFMGRelaxType[ex->relax_type],&ex->relax_type,NULL);
2493:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type));
2494:   PetscOptionsTail();
2495:   return(0);
2496: }

2498: PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
2499: {
2500:   PetscErrorCode    ierr;
2501:   PC_SysPFMG        *ex = (PC_SysPFMG*) pc->data;
2502:   PetscScalar       *yy;
2503:   const PetscScalar *xx;
2504:   PetscInt          ilower[3],iupper[3];
2505:   HYPRE_Int         hlower[3],hupper[3];
2506:   Mat_HYPRESStruct  *mx     = (Mat_HYPRESStruct*)(pc->pmat->data);
2507:   PetscInt          ordering= mx->dofs_order;
2508:   PetscInt          nvars   = mx->nvars;
2509:   PetscInt          part    = 0;
2510:   PetscInt          size;
2511:   PetscInt          i;

2514:   PetscCitationsRegister(hypreCitation,&cite);
2515:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
2516:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2517:   iupper[0] += ilower[0] - 1;
2518:   iupper[1] += ilower[1] - 1;
2519:   iupper[2] += ilower[2] - 1;
2520:   hlower[0]  = (HYPRE_Int)ilower[0];
2521:   hlower[1]  = (HYPRE_Int)ilower[1];
2522:   hlower[2]  = (HYPRE_Int)ilower[2];
2523:   hupper[0]  = (HYPRE_Int)iupper[0];
2524:   hupper[1]  = (HYPRE_Int)iupper[1];
2525:   hupper[2]  = (HYPRE_Int)iupper[2];

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

2530:   /* copy x values over to hypre for variable ordering */
2531:   if (ordering) {
2532:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2533:     VecGetArrayRead(x,&xx);
2534:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(xx+(size*i))));
2535:     VecRestoreArrayRead(x,&xx);
2536:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2537:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
2538:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

2540:     /* copy solution values back to PETSc */
2541:     VecGetArray(y,&yy);
2542:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(yy+(size*i))));
2543:     VecRestoreArray(y,&yy);
2544:   } else {      /* nodal ordering must be mapped to variable ordering for sys_pfmg */
2545:     PetscScalar *z;
2546:     PetscInt    j, k;

2548:     PetscMalloc1(nvars*size,&z);
2549:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2550:     VecGetArrayRead(x,&xx);

2552:     /* transform nodal to hypre's variable ordering for sys_pfmg */
2553:     for (i= 0; i< size; i++) {
2554:       k= i*nvars;
2555:       for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
2556:     }
2557:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i))));
2558:     VecRestoreArrayRead(x,&xx);
2559:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2560:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

2562:     /* copy solution values back to PETSc */
2563:     VecGetArray(y,&yy);
2564:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i))));
2565:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
2566:     for (i= 0; i< size; i++) {
2567:       k= i*nvars;
2568:       for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
2569:     }
2570:     VecRestoreArray(y,&yy);
2571:     PetscFree(z);
2572:   }
2573:   return(0);
2574: }

2576: 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)
2577: {
2578:   PC_SysPFMG     *jac = (PC_SysPFMG*)pc->data;
2580:   HYPRE_Int      oits;

2583:   PetscCitationsRegister(hypreCitation,&cite);
2584:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its));
2585:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol));
2586:   PCApply_SysPFMG(pc,b,y);
2587:   PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,&oits));
2588:   *outits = oits;
2589:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2590:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
2591:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol));
2592:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its));
2593:   return(0);
2594: }

2596: PetscErrorCode PCSetUp_SysPFMG(PC pc)
2597: {
2598:   PetscErrorCode   ierr;
2599:   PC_SysPFMG       *ex = (PC_SysPFMG*) pc->data;
2600:   Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
2601:   PetscBool        flg;

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

2607:   /* create the hypre sstruct solver object and set its information */
2608:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2609:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2610:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver));
2611:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2612:   return(0);
2613: }

2615: /*MC
2616:      PCSysPFMG - the hypre SysPFMG multigrid solver

2618:    Level: advanced

2620:    Options Database:
2621: + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner
2622: . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2623: . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2624: . -pc_syspfmg_tol <tol> tolerance of SysPFMG
2625: - -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel

2627:    Notes:
2628:     This is for CELL-centered descretizations

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

2634: .seealso:  PCMG, MATHYPRESSTRUCT
2635: M*/

2637: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
2638: {
2640:   PC_SysPFMG     *ex;

2643:   PetscNew(&ex); \
2644:   pc->data = ex;

2646:   ex->its            = 1;
2647:   ex->tol            = 1.e-8;
2648:   ex->relax_type     = 1;
2649:   ex->num_pre_relax  = 1;
2650:   ex->num_post_relax = 1;

2652:   pc->ops->setfromoptions  = PCSetFromOptions_SysPFMG;
2653:   pc->ops->view            = PCView_SysPFMG;
2654:   pc->ops->destroy         = PCDestroy_SysPFMG;
2655:   pc->ops->apply           = PCApply_SysPFMG;
2656:   pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
2657:   pc->ops->setup           = PCSetUp_SysPFMG;

2659:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
2660:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2661:   return(0);
2662: }