Actual source code: hypre.c

petsc-master 2019-06-19
Report Typos and Errors

  2: /*
  3:    Provides an interface to the LLNL package hypre
  4: */

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

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

 17: static PetscBool cite = PETSC_FALSE;
 18: 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";

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

 27:   HYPRE_Int (*destroy)(HYPRE_Solver);
 28:   HYPRE_Int (*solve)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);
 29:   HYPRE_Int (*setup)(HYPRE_Solver,HYPRE_ParCSRMatrix,HYPRE_ParVector,HYPRE_ParVector);

 31:   MPI_Comm comm_hypre;
 32:   char     *hypre_type;

 34:   /* options for Pilut and BoomerAMG*/
 35:   PetscInt  maxiter;
 36:   PetscReal tol;

 38:   /* options for Pilut */
 39:   PetscInt factorrowsize;

 41:   /* options for ParaSails */
 42:   PetscInt  nlevels;
 43:   PetscReal threshhold;
 44:   PetscReal filter;
 45:   PetscInt  sym;
 46:   PetscReal loadbal;
 47:   PetscInt  logging;
 48:   PetscInt  ruse;
 49:   PetscInt  symt;

 51:   /* options for BoomerAMG */
 52:   PetscBool printstatistics;

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

 80:   PetscInt  nodal_coarsening;
 81:   PetscInt  nodal_coarsening_diag;
 82:   PetscInt  vec_interp_variant;
 83:   PetscInt  vec_interp_qmax;
 84:   PetscBool vec_interp_smooth;
 85:   PetscInt  interp_refine;

 87:   HYPRE_IJVector  *hmnull;
 88:   HYPRE_ParVector *phmnull;  /* near null space passed to hypre */
 89:   PetscInt        n_hmnull;
 90:   Vec             hmnull_constant;
 91:   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 */

 93:   /* options for AS (Auxiliary Space preconditioners) */
 94:   PetscInt  as_print;
 95:   PetscInt  as_max_iter;
 96:   PetscReal as_tol;
 97:   PetscInt  as_relax_type;
 98:   PetscInt  as_relax_times;
 99:   PetscReal as_relax_weight;
100:   PetscReal as_omega;
101:   PetscInt  as_amg_alpha_opts[5]; /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for vector Poisson (AMS) or Curl problem (ADS) */
102:   PetscReal as_amg_alpha_theta;   /* AMG strength for vector Poisson (AMS) or Curl problem (ADS) */
103:   PetscInt  as_amg_beta_opts[5];  /* AMG coarsen type, agg_levels, relax_type, interp_type, Pmax for scalar Poisson (AMS) or vector Poisson (ADS) */
104:   PetscReal as_amg_beta_theta;    /* AMG strength for scalar Poisson (AMS) or vector Poisson (ADS)  */
105:   PetscInt  ams_cycle_type;
106:   PetscInt  ads_cycle_type;

108:   /* additional data */
109:   Mat G;             /* MatHYPRE */
110:   Mat C;             /* MatHYPRE */
111:   Mat alpha_Poisson; /* MatHYPRE */
112:   Mat beta_Poisson;  /* MatHYPRE */

114:   /* extra information for AMS */
115:   PetscInt       dim; /* geometrical dimension */
116:   HYPRE_IJVector coords[3];
117:   HYPRE_IJVector constants[3];
118:   Mat            RT_PiFull, RT_Pi[3];
119:   Mat            ND_PiFull, ND_Pi[3];
120:   PetscBool      ams_beta_is_zero;
121:   PetscBool      ams_beta_is_zero_part;
122:   PetscInt       ams_proj_freq;
123: } PC_HYPRE;

125: PetscErrorCode PCHYPREGetSolver(PC pc,HYPRE_Solver *hsolver)
126: {
127:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

130:   *hsolver = jac->hsolver;
131:   return(0);
132: }

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

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

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

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

194: /* Resets (frees) Hypre's representation of the near null space */
195: static PetscErrorCode PCHYPREResetNearNullSpace_Private(PC pc)
196: {
197:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
198:   PetscInt       i;

202:   for (i=0; i<jac->n_hmnull; i++) {
203:     PETSC_UNUSED HYPRE_Complex *harray;
204:     VecHYPRE_ParVectorReplacePointer(jac->hmnull[i],jac->hmnull_hypre_data_array[i],harray);
205:     PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->hmnull[i]));
206:   }
207:   PetscFree(jac->hmnull);
208:   PetscFree(jac->hmnull_hypre_data_array);
209:   PetscFree(jac->phmnull);
210:   VecDestroy(&jac->hmnull_constant);
211:   jac->n_hmnull = 0;
212:   return(0);
213: }

215: static PetscErrorCode PCSetUp_HYPRE(PC pc)
216: {
217:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
218:   Mat_HYPRE          *hjac;
219:   HYPRE_ParCSRMatrix hmat;
220:   HYPRE_ParVector    bv,xv;
221:   PetscBool          ishypre;
222:   PetscErrorCode     ierr;

225:   if (!jac->hypre_type) {
226:     PCHYPRESetType(pc,"boomeramg");
227:   }

229:   PetscObjectTypeCompare((PetscObject)pc->pmat,MATHYPRE,&ishypre);
230:   if (!ishypre) {
231:     MatDestroy(&jac->hpmat);
232:     MatConvert(pc->pmat,MATHYPRE,MAT_INITIAL_MATRIX,&jac->hpmat);
233:   } else {
234:     PetscObjectReference((PetscObject)pc->pmat);
235:     MatDestroy(&jac->hpmat);
236:     jac->hpmat = pc->pmat;
237:   }
238:   hjac = (Mat_HYPRE*)(jac->hpmat->data);

240:   /* special case for BoomerAMG */
241:   if (jac->setup == HYPRE_BoomerAMGSetup) {
242:     MatNullSpace    mnull;
243:     PetscBool       has_const;
244:     PetscInt        bs,nvec,i;
245:     const Vec       *vecs;
246:     HYPRE_Complex   *petscvecarray;

248:     MatGetBlockSize(pc->pmat,&bs);
249:     if (bs > 1) PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
250:     MatGetNearNullSpace(pc->mat, &mnull);
251:     if (mnull) {
252:       PCHYPREResetNearNullSpace_Private(pc);
253:       MatNullSpaceGetVecs(mnull, &has_const, &nvec, &vecs);
254:       PetscMalloc1(nvec+1,&jac->hmnull);
255:       PetscMalloc1(nvec+1,&jac->hmnull_hypre_data_array);
256:       PetscMalloc1(nvec+1,&jac->phmnull);
257:       for (i=0; i<nvec; i++) {
258:         VecHYPRE_IJVectorCreate(vecs[i],&jac->hmnull[i]);
259:         VecGetArrayRead(vecs[i],(const PetscScalar **)&petscvecarray);
260:         VecHYPRE_ParVectorReplacePointer(jac->hmnull[i],petscvecarray,jac->hmnull_hypre_data_array[i]);
261:         VecRestoreArrayRead(vecs[i],(const PetscScalar **)&petscvecarray);
262:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[i],(void**)&jac->phmnull[i]));
263:       }
264:       if (has_const) {
265:         MatCreateVecs(pc->pmat,&jac->hmnull_constant,NULL);
266:         VecSet(jac->hmnull_constant,1);
267:         VecNormalize(jac->hmnull_constant,NULL);
268:         VecHYPRE_IJVectorCreate(jac->hmnull_constant,&jac->hmnull[nvec]);
269:         VecGetArrayRead(jac->hmnull_constant,(const PetscScalar **)&petscvecarray);
270:         VecHYPRE_ParVectorReplacePointer(jac->hmnull[nvec],petscvecarray,jac->hmnull_hypre_data_array[nvec]);
271:         VecRestoreArrayRead(jac->hmnull_constant,(const PetscScalar **)&petscvecarray);
272:         PetscStackCallStandard(HYPRE_IJVectorGetObject,(jac->hmnull[nvec],(void**)&jac->phmnull[nvec]));
273:         nvec++;
274:       }
275:       PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVectors,(jac->hsolver,nvec,jac->phmnull));
276:       jac->n_hmnull = nvec;
277:     }
278:   }

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

413: static PetscErrorCode PCApply_HYPRE(PC pc,Vec b,Vec x)
414: {
415:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
416:   Mat_HYPRE          *hjac = (Mat_HYPRE*)(jac->hpmat->data);
417:   PetscErrorCode     ierr;
418:   HYPRE_ParCSRMatrix hmat;
419:   HYPRE_Complex      *xv,*sxv;
420:   HYPRE_Complex      *bv,*sbv;
421:   HYPRE_ParVector    jbv,jxv;
422:   PetscInt           hierr;

425:   PetscCitationsRegister(hypreCitation,&cite);
426:   if (!jac->applyrichardson) {VecSet(x,0.0);}
427:   VecGetArrayRead(b,(const PetscScalar **)&bv);
428:   VecGetArray(x,(PetscScalar **)&xv);
429:   VecHYPRE_ParVectorReplacePointer(hjac->b,bv,sbv);
430:   VecHYPRE_ParVectorReplacePointer(hjac->x,xv,sxv);

432:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
433:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b,(void**)&jbv));
434:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x,(void**)&jxv));
435:   PetscStackCall("Hypre solve",h(*jac->solve)(jac->hsolver,hmat,jbv,jxv);
436:   if (hierr && hierr != HYPRE_ERROR_CONV) SETERRQ1(PETSC_COMM_SELF,PETSC_ERR_LIB,"Error in HYPRE solver, error code %d",hierr);
437:   if (hierr) hypre__global_error = 0;);

439:   if (jac->setup == HYPRE_AMSSetup && jac->ams_beta_is_zero_part) {
440:     PetscStackCallStandard(HYPRE_AMSProjectOutGradients,(jac->hsolver,jxv));
441:   }
442:   VecHYPRE_ParVectorReplacePointer(hjac->b,sbv,bv);
443:   VecHYPRE_ParVectorReplacePointer(hjac->x,sxv,xv);
444:   VecRestoreArray(x,(PetscScalar **)&xv);
445:   VecRestoreArrayRead(b,(const PetscScalar **)&bv);
446:   return(0);
447: }

449: static PetscErrorCode PCReset_HYPRE(PC pc)
450: {
451:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;

455:   MatDestroy(&jac->hpmat);
456:   MatDestroy(&jac->G);
457:   MatDestroy(&jac->C);
458:   MatDestroy(&jac->alpha_Poisson);
459:   MatDestroy(&jac->beta_Poisson);
460:   MatDestroy(&jac->RT_PiFull);
461:   MatDestroy(&jac->RT_Pi[0]);
462:   MatDestroy(&jac->RT_Pi[1]);
463:   MatDestroy(&jac->RT_Pi[2]);
464:   MatDestroy(&jac->ND_PiFull);
465:   MatDestroy(&jac->ND_Pi[0]);
466:   MatDestroy(&jac->ND_Pi[1]);
467:   MatDestroy(&jac->ND_Pi[2]);
468:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0])); jac->coords[0] = NULL;
469:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1])); jac->coords[1] = NULL;
470:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2])); jac->coords[2] = NULL;
471:   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0])); jac->constants[0] = NULL;
472:   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1])); jac->constants[1] = NULL;
473:   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2])); jac->constants[2] = NULL;
474:   PCHYPREResetNearNullSpace_Private(pc);
475:   jac->ams_beta_is_zero = PETSC_FALSE;
476:   jac->dim = 0;
477:   return(0);
478: }

480: static PetscErrorCode PCDestroy_HYPRE(PC pc)
481: {
482:   PC_HYPRE                 *jac = (PC_HYPRE*)pc->data;
483:   PetscErrorCode           ierr;

486:   PCReset_HYPRE(pc);
487:   if (jac->destroy) PetscStackCallStandard(jac->destroy,(jac->hsolver));
488:   PetscFree(jac->hypre_type);
489:   if (jac->comm_hypre != MPI_COMM_NULL) { MPI_Comm_free(&(jac->comm_hypre));}
490:   PetscFree(pc->data);

492:   PetscObjectChangeTypeName((PetscObject)pc,0);
493:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",NULL);
494:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",NULL);
495:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetCoordinates_C",NULL);
496:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",NULL);
497:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",NULL);
498:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",NULL);
499:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetConstantEdgeVectors_C",NULL);
500:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",NULL);
501:   PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",NULL);
502:   PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",NULL);
503:   return(0);
504: }

506: /* --------------------------------------------------------------------------------------------*/
507: static PetscErrorCode PCSetFromOptions_HYPRE_Pilut(PetscOptionItems *PetscOptionsObject,PC pc)
508: {
509:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
511:   PetscBool      flag;

514:   PetscOptionsHead(PetscOptionsObject,"HYPRE Pilut Options");
515:   PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
516:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetMaxIter,(jac->hsolver,jac->maxiter));
517:   PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
518:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetDropTolerance,(jac->hsolver,jac->tol));
519:   PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
520:   if (flag) PetscStackCallStandard(HYPRE_ParCSRPilutSetFactorRowSize,(jac->hsolver,jac->factorrowsize));
521:   PetscOptionsTail();
522:   return(0);
523: }

525: static PetscErrorCode PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
526: {
527:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
529:   PetscBool      iascii;

532:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
533:   if (iascii) {
534:     PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut preconditioning\n");
535:     if (jac->maxiter != PETSC_DEFAULT) {
536:       PetscViewerASCIIPrintf(viewer,"    maximum number of iterations %d\n",jac->maxiter);
537:     } else {
538:       PetscViewerASCIIPrintf(viewer,"    default maximum number of iterations \n");
539:     }
540:     if (jac->tol != PETSC_DEFAULT) {
541:       PetscViewerASCIIPrintf(viewer,"    drop tolerance %g\n",(double)jac->tol);
542:     } else {
543:       PetscViewerASCIIPrintf(viewer,"    default drop tolerance \n");
544:     }
545:     if (jac->factorrowsize != PETSC_DEFAULT) {
546:       PetscViewerASCIIPrintf(viewer,"    factor row size %d\n",jac->factorrowsize);
547:     } else {
548:       PetscViewerASCIIPrintf(viewer,"    default factor row size \n");
549:     }
550:   }
551:   return(0);
552: }

554: /* --------------------------------------------------------------------------------------------*/
555: static PetscErrorCode PCSetFromOptions_HYPRE_Euclid(PetscOptionItems *PetscOptionsObject,PC pc)
556: {
557:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
559:   PetscBool      flag;

562:   PetscOptionsHead(PetscOptionsObject,"HYPRE Euclid Options");
563:   PetscOptionsInt("-pc_hypre_euclid_level","Factorization levels","None",jac->eu_level,&jac->eu_level,&flag);
564:   if (flag) PetscStackCallStandard(HYPRE_EuclidSetLevel,(jac->hsolver,jac->eu_level));
565:   PetscOptionsTail();
566:   return(0);
567: }

569: static PetscErrorCode PCView_HYPRE_Euclid(PC pc,PetscViewer viewer)
570: {
571:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
573:   PetscBool      iascii;

576:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
577:   if (iascii) {
578:     PetscViewerASCIIPrintf(viewer,"  HYPRE Euclid preconditioning\n");
579:     if (jac->eu_level != PETSC_DEFAULT) {
580:       PetscViewerASCIIPrintf(viewer,"    factorization levels %d\n",jac->eu_level);
581:     } else {
582:       PetscViewerASCIIPrintf(viewer,"    default factorization levels \n");
583:     }
584:   }
585:   return(0);
586: }

588: /* --------------------------------------------------------------------------------------------*/

590: static PetscErrorCode PCApplyTranspose_HYPRE_BoomerAMG(PC pc,Vec b,Vec x)
591: {
592:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
593:   Mat_HYPRE          *hjac = (Mat_HYPRE*)(jac->hpmat->data);
594:   PetscErrorCode     ierr;
595:   HYPRE_ParCSRMatrix hmat;
596:   HYPRE_Complex      *xv,*bv;
597:   HYPRE_Complex      *sbv,*sxv;
598:   HYPRE_ParVector    jbv,jxv;
599:   PetscInt           hierr;

602:   PetscCitationsRegister(hypreCitation,&cite);
603:   VecSet(x,0.0);
604:   VecGetArrayRead(b,(const PetscScalar**)&bv);
605:   VecGetArray(x,(PetscScalar**)&xv);
606:   VecHYPRE_ParVectorReplacePointer(hjac->b,bv,sbv);
607:   VecHYPRE_ParVectorReplacePointer(hjac->x,xv,sxv);

609:   PetscStackCallStandard(HYPRE_IJMatrixGetObject,(hjac->ij,(void**)&hmat));
610:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->b,(void**)&jbv));
611:   PetscStackCallStandard(HYPRE_IJVectorGetObject,(hjac->x,(void**)&jxv));

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

618:   VecHYPRE_ParVectorReplacePointer(hjac->b,sbv,bv);
619:   VecHYPRE_ParVectorReplacePointer(hjac->x,sxv,xv);
620:   VecRestoreArray(x,(PetscScalar**)&xv);
621:   VecRestoreArrayRead(b,(const PetscScalar**)&bv);
622:   return(0);
623: }

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

628: static const char *HYPREBoomerAMGCycleType[]   = {"","V","W"};
629: static const char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout", "", "PMIS", "", "HMIS"};
630: static const char *HYPREBoomerAMGMeasureType[] = {"local","global"};
631: /* The following corresponds to HYPRE_BoomerAMGSetRelaxType which has many missing numbers in the enum */
632: static const char *HYPREBoomerAMGSmoothType[]   = {"Schwarz-smoothers","Pilut","ParaSails","Euclid"};
633: static const char *HYPREBoomerAMGRelaxType[]   = {"Jacobi","sequential-Gauss-Seidel","seqboundary-Gauss-Seidel","SOR/Jacobi","backward-SOR/Jacobi",
634:                                                   "" /* [5] hybrid chaotic Gauss-Seidel (works only with OpenMP) */,"symmetric-SOR/Jacobi",
635:                                                   "" /* 7 */,"l1scaled-SOR/Jacobi","Gaussian-elimination",
636:                                                   "" /* 10 */, "" /* 11 */, "" /* 12 */, "l1-Gauss-Seidel" /* nonsymmetric */, "backward-l1-Gauss-Seidel" /* nonsymmetric */,
637:                                                   "CG" /* non-stationary */,"Chebyshev","FCF-Jacobi","l1scaled-Jacobi"};
638: static const char *HYPREBoomerAMGInterpType[]  = {"classical", "", "", "direct", "multipass", "multipass-wts", "ext+i",
639:                                                   "ext+i-cc", "standard", "standard-wts", "block", "block-wtd", "FF", "FF1"};
640: static PetscErrorCode PCSetFromOptions_HYPRE_BoomerAMG(PetscOptionItems *PetscOptionsObject,PC pc)
641: {
642:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
644:   PetscInt       bs,n,indx,level;
645:   PetscBool      flg, tmp_truth;
646:   double         tmpdbl, twodbl[2];

649:   PetscOptionsHead(PetscOptionsObject,"HYPRE BoomerAMG Options");
650:   PetscOptionsEList("-pc_hypre_boomeramg_cycle_type","Cycle type","None",HYPREBoomerAMGCycleType+1,2,HYPREBoomerAMGCycleType[jac->cycletype],&indx,&flg);
651:   if (flg) {
652:     jac->cycletype = indx+1;
653:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
654:   }
655:   PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
656:   if (flg) {
657:     if (jac->maxlevels < 2) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of levels %d must be at least two",jac->maxlevels);
658:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
659:   }
660:   PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used PER hypre call","None",jac->maxiter,&jac->maxiter,&flg);
661:   if (flg) {
662:     if (jac->maxiter < 1) SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_OUTOFRANGE,"Number of iterations %d must be at least one",jac->maxiter);
663:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
664:   }
665:   PetscOptionsReal("-pc_hypre_boomeramg_tol","Convergence tolerance PER hypre call (0.0 = use a fixed number of iterations)","None",jac->tol,&jac->tol,&flg);
666:   if (flg) {
667:     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);
668:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
669:   }
670:   bs = 1;
671:   if (pc->pmat) {
672:     MatGetBlockSize(pc->pmat,&bs);
673:   }
674:   PetscOptionsInt("-pc_hypre_boomeramg_numfunctions","Number of functions","HYPRE_BoomerAMGSetNumFunctions",bs,&bs,&flg);
675:   if (flg) {
676:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumFunctions,(jac->hsolver,bs));
677:   }

679:   PetscOptionsReal("-pc_hypre_boomeramg_truncfactor","Truncation factor for interpolation (0=no truncation)","None",jac->truncfactor,&jac->truncfactor,&flg);
680:   if (flg) {
681:     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);
682:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
683:   }

685:   PetscOptionsInt("-pc_hypre_boomeramg_P_max","Max elements per row for interpolation operator (0=unlimited)","None",jac->pmax,&jac->pmax,&flg);
686:   if (flg) {
687:     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);
688:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
689:   }

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

695:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
696:   }

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

702:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
703:   }


706:   PetscOptionsReal("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
707:   if (flg) {
708:     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);
709:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
710:   }
711:   PetscOptionsReal("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
712:   if (flg) {
713:     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);
714:     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);
715:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
716:   }

718:   /* Grid sweeps */
719:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_all","Number of sweeps for the up and down grid levels","None",jac->gridsweeps[0],&indx,&flg);
720:   if (flg) {
721:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver,indx));
722:     /* modify the jac structure so we can view the updated options with PC_View */
723:     jac->gridsweeps[0] = indx;
724:     jac->gridsweeps[1] = indx;
725:     /*defaults coarse to 1 */
726:     jac->gridsweeps[2] = 1;
727:   }
728:   PetscOptionsInt("-pc_hypre_boomeramg_nodal_coarsen","Use a nodal based coarsening 1-6","HYPRE_BoomerAMGSetNodal",jac->nodal_coarsening,&jac->nodal_coarsening,&flg);
729:   if (flg) {
730:     PetscStackCallStandard(HYPRE_BoomerAMGSetNodal,(jac->hsolver,jac->nodal_coarsening));
731:   }
732:   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);
733:   if (flg) {
734:     PetscStackCallStandard(HYPRE_BoomerAMGSetNodalDiag,(jac->hsolver,jac->nodal_coarsening_diag));
735:   }
736:   PetscOptionsInt("-pc_hypre_boomeramg_vec_interp_variant","Variant of algorithm 1-3","HYPRE_BoomerAMGSetInterpVecVariant",jac->vec_interp_variant, &jac->vec_interp_variant,&flg);
737:   if (flg) {
738:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecVariant,(jac->hsolver,jac->vec_interp_variant));
739:   }
740:   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);
741:   if (flg) {
742:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpVecQMax,(jac->hsolver,jac->vec_interp_qmax));
743:   }
744:   PetscOptionsBool("-pc_hypre_boomeramg_vec_interp_smooth","Whether to smooth the interpolation vectors","HYPRE_BoomerAMGSetSmoothInterpVectors",jac->vec_interp_smooth, &jac->vec_interp_smooth,&flg);
745:   if (flg) {
746:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothInterpVectors,(jac->hsolver,jac->vec_interp_smooth));
747:   }
748:   PetscOptionsInt("-pc_hypre_boomeramg_interp_refine","Preprocess the interpolation matrix through iterative weight refinement","HYPRE_BoomerAMGSetInterpRefine",jac->interp_refine, &jac->interp_refine,&flg);
749:   if (flg) {
750:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpRefine,(jac->hsolver,jac->interp_refine));
751:   }
752:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_down","Number of sweeps for the down cycles","None",jac->gridsweeps[0], &indx,&flg);
753:   if (flg) {
754:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 1));
755:     jac->gridsweeps[0] = indx;
756:   }
757:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_up","Number of sweeps for the up cycles","None",jac->gridsweeps[1],&indx,&flg);
758:   if (flg) {
759:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 2));
760:     jac->gridsweeps[1] = indx;
761:   }
762:   PetscOptionsInt("-pc_hypre_boomeramg_grid_sweeps_coarse","Number of sweeps for the coarse level","None",jac->gridsweeps[2],&indx,&flg);
763:   if (flg) {
764:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleNumSweeps,(jac->hsolver,indx, 3));
765:     jac->gridsweeps[2] = indx;
766:   }

768:   /* Smooth type */
769:   PetscOptionsEList("-pc_hypre_boomeramg_smooth_type","Enable more complex smoothers","None",HYPREBoomerAMGSmoothType,ALEN(HYPREBoomerAMGSmoothType),HYPREBoomerAMGSmoothType[0],&indx,&flg);
770:   if (flg) {
771:     jac->smoothtype = indx;
772:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,indx+6));
773:     jac->smoothnumlevels = 25;
774:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,25));
775:   }

777:   /* Number of smoothing levels */
778:   PetscOptionsInt("-pc_hypre_boomeramg_smooth_num_levels","Number of levels on which more complex smoothers are used","None",25,&indx,&flg);
779:   if (flg && (jac->smoothtype != -1)) {
780:     jac->smoothnumlevels = indx;
781:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,indx));
782:   }

784:   /* Number of levels for ILU(k) for Euclid */
785:   PetscOptionsInt("-pc_hypre_boomeramg_eu_level","Number of levels for ILU(k) in Euclid smoother","None",0,&indx,&flg);
786:   if (flg && (jac->smoothtype == 3)) {
787:     jac->eu_level = indx;
788:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,indx));
789:   }

791:   /* Filter for ILU(k) for Euclid */
792:   double droptolerance;
793:   PetscOptionsReal("-pc_hypre_boomeramg_eu_droptolerance","Drop tolerance for ILU(k) in Euclid smoother","None",0,&droptolerance,&flg);
794:   if (flg && (jac->smoothtype == 3)) {
795:     jac->eu_droptolerance = droptolerance;
796:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuLevel,(jac->hsolver,droptolerance));
797:   }

799:   /* Use Block Jacobi ILUT for Euclid */
800:   PetscOptionsBool("-pc_hypre_boomeramg_eu_bj", "Use Block Jacobi for ILU in Euclid smoother?", "None", PETSC_FALSE, &tmp_truth, &flg);
801:   if (flg && (jac->smoothtype == 3)) {
802:     jac->eu_bj = tmp_truth;
803:     PetscStackCallStandard(HYPRE_BoomerAMGSetEuBJ,(jac->hsolver,jac->eu_bj));
804:   }

806:   /* Relax type */
807:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_all","Relax type for the up and down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
808:   if (flg) {
809:     jac->relaxtype[0] = jac->relaxtype[1]  = indx;
810:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, indx));
811:     /* by default, coarse type set to 9 */
812:     jac->relaxtype[2] = 9;
813:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, 9, 3));
814:   }
815:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_down","Relax type for the down cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
816:   if (flg) {
817:     jac->relaxtype[0] = indx;
818:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 1));
819:   }
820:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_up","Relax type for the up cycles","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[6],&indx,&flg);
821:   if (flg) {
822:     jac->relaxtype[1] = indx;
823:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 2));
824:   }
825:   PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,ALEN(HYPREBoomerAMGRelaxType),HYPREBoomerAMGRelaxType[9],&indx,&flg);
826:   if (flg) {
827:     jac->relaxtype[2] = indx;
828:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleRelaxType,(jac->hsolver, indx, 3));
829:   }

831:   /* Relaxation Weight */
832:   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);
833:   if (flg) {
834:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxWt,(jac->hsolver,tmpdbl));
835:     jac->relaxweight = tmpdbl;
836:   }

838:   n         = 2;
839:   twodbl[0] = twodbl[1] = 1.0;
840:   PetscOptionsRealArray("-pc_hypre_boomeramg_relax_weight_level","Set the relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
841:   if (flg) {
842:     if (n == 2) {
843:       indx =  (int)PetscAbsReal(twodbl[1]);
844:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelRelaxWt,(jac->hsolver,twodbl[0],indx));
845:     } 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);
846:   }

848:   /* Outer relaxation Weight */
849:   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);
850:   if (flg) {
851:     PetscStackCallStandard(HYPRE_BoomerAMGSetOuterWt,(jac->hsolver, tmpdbl));
852:     jac->outerrelaxweight = tmpdbl;
853:   }

855:   n         = 2;
856:   twodbl[0] = twodbl[1] = 1.0;
857:   PetscOptionsRealArray("-pc_hypre_boomeramg_outer_relax_weight_level","Set the outer relaxation weight for a particular level (weight,level)","None",twodbl, &n, &flg);
858:   if (flg) {
859:     if (n == 2) {
860:       indx =  (int)PetscAbsReal(twodbl[1]);
861:       PetscStackCallStandard(HYPRE_BoomerAMGSetLevelOuterWt,(jac->hsolver, twodbl[0], indx));
862:     } 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);
863:   }

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

868:   if (flg && tmp_truth) {
869:     jac->relaxorder = 0;
870:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
871:   }
872:   PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,ALEN(HYPREBoomerAMGMeasureType),HYPREBoomerAMGMeasureType[0],&indx,&flg);
873:   if (flg) {
874:     jac->measuretype = indx;
875:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
876:   }
877:   /* update list length 3/07 */
878:   PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,ALEN(HYPREBoomerAMGCoarsenType),HYPREBoomerAMGCoarsenType[6],&indx,&flg);
879:   if (flg) {
880:     jac->coarsentype = indx;
881:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
882:   }

884:   /* new 3/07 */
885:   PetscOptionsEList("-pc_hypre_boomeramg_interp_type","Interpolation type","None",HYPREBoomerAMGInterpType,ALEN(HYPREBoomerAMGInterpType),HYPREBoomerAMGInterpType[0],&indx,&flg);
886:   if (flg) {
887:     jac->interptype = indx;
888:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
889:   }

891:   PetscOptionsName("-pc_hypre_boomeramg_print_statistics","Print statistics","None",&flg);
892:   if (flg) {
893:     level = 3;
894:     PetscOptionsInt("-pc_hypre_boomeramg_print_statistics","Print statistics","None",level,&level,NULL);

896:     jac->printstatistics = PETSC_TRUE;
897:     PetscStackCallStandard(HYPRE_BoomerAMGSetPrintLevel,(jac->hsolver,level));
898:   }

900:   PetscOptionsName("-pc_hypre_boomeramg_print_debug","Print debug information","None",&flg);
901:   if (flg) {
902:     level = 3;
903:     PetscOptionsInt("-pc_hypre_boomeramg_print_debug","Print debug information","None",level,&level,NULL);

905:     jac->printstatistics = PETSC_TRUE;
906:     PetscStackCallStandard(HYPRE_BoomerAMGSetDebugFlag,(jac->hsolver,level));
907:   }

909:   PetscOptionsBool("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None", PETSC_FALSE, &tmp_truth, &flg);
910:   if (flg && tmp_truth) {
911:     PetscInt tmp_int;
912:     PetscOptionsInt("-pc_hypre_boomeramg_nodal_relaxation", "Nodal relaxation via Schwarz", "None",jac->nodal_relax_levels,&tmp_int,&flg);
913:     if (flg) jac->nodal_relax_levels = tmp_int;
914:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothType,(jac->hsolver,6));
915:     PetscStackCallStandard(HYPRE_BoomerAMGSetDomainType,(jac->hsolver,1));
916:     PetscStackCallStandard(HYPRE_BoomerAMGSetOverlap,(jac->hsolver,0));
917:     PetscStackCallStandard(HYPRE_BoomerAMGSetSmoothNumLevels,(jac->hsolver,jac->nodal_relax_levels));
918:   }

920:   PetscOptionsTail();
921:   return(0);
922: }

924: 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)
925: {
926:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
928:   HYPRE_Int      oits;

931:   PetscCitationsRegister(hypreCitation,&cite);
932:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,its*jac->maxiter));
933:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,rtol));
934:   jac->applyrichardson = PETSC_TRUE;
935:   PCApply_HYPRE(pc,b,y);
936:   jac->applyrichardson = PETSC_FALSE;
937:   PetscStackCallStandard(HYPRE_BoomerAMGGetNumIterations,(jac->hsolver,&oits));
938:   *outits = oits;
939:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
940:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
941:   PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
942:   PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
943:   return(0);
944: }


947: static PetscErrorCode PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
948: {
949:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
951:   PetscBool      iascii;

954:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
955:   if (iascii) {
956:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG preconditioning\n");
957:     PetscViewerASCIIPrintf(viewer,"    Cycle type %s\n",HYPREBoomerAMGCycleType[jac->cycletype]);
958:     PetscViewerASCIIPrintf(viewer,"    Maximum number of levels %D\n",jac->maxlevels);
959:     PetscViewerASCIIPrintf(viewer,"    Maximum number of iterations PER hypre call %D\n",jac->maxiter);
960:     PetscViewerASCIIPrintf(viewer,"    Convergence tolerance PER hypre call %g\n",(double)jac->tol);
961:     PetscViewerASCIIPrintf(viewer,"    Threshold for strong coupling %g\n",(double)jac->strongthreshold);
962:     PetscViewerASCIIPrintf(viewer,"    Interpolation truncation factor %g\n",(double)jac->truncfactor);
963:     PetscViewerASCIIPrintf(viewer,"    Interpolation: max elements per row %D\n",jac->pmax);
964:     if (jac->interp_refine) {
965:       PetscViewerASCIIPrintf(viewer,"    Interpolation: number of steps of weighted refinement %D\n",jac->interp_refine);
966:     }
967:     PetscViewerASCIIPrintf(viewer,"    Number of levels of aggressive coarsening %D\n",jac->agg_nl);
968:     PetscViewerASCIIPrintf(viewer,"    Number of paths for aggressive coarsening %D\n",jac->agg_num_paths);

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

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

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

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

983:     if (jac->relaxorder) {
984:       PetscViewerASCIIPrintf(viewer,"    Using CF-relaxation\n");
985:     } else {
986:       PetscViewerASCIIPrintf(viewer,"    Not using CF-relaxation\n");
987:     }
988:     if (jac->smoothtype!=-1) {
989:       PetscViewerASCIIPrintf(viewer,"    Smooth type          %s\n",HYPREBoomerAMGSmoothType[jac->smoothtype]);
990:       PetscViewerASCIIPrintf(viewer,"    Smooth num levels    %D\n",jac->smoothnumlevels);
991:     } else {
992:       PetscViewerASCIIPrintf(viewer,"    Not using more complex smoothers.\n");
993:     }
994:     if (jac->smoothtype==3) {
995:       PetscViewerASCIIPrintf(viewer,"    Euclid ILU(k) levels %D\n",jac->eu_level);
996:       PetscViewerASCIIPrintf(viewer,"    Euclid ILU(k) drop tolerance %g\n",(double)jac->eu_droptolerance);
997:       PetscViewerASCIIPrintf(viewer,"    Euclid ILU use Block-Jacobi? %D\n",jac->eu_bj);
998:     }
999:     PetscViewerASCIIPrintf(viewer,"    Measure type        %s\n",HYPREBoomerAMGMeasureType[jac->measuretype]);
1000:     PetscViewerASCIIPrintf(viewer,"    Coarsen type        %s\n",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
1001:     PetscViewerASCIIPrintf(viewer,"    Interpolation type  %s\n",HYPREBoomerAMGInterpType[jac->interptype]);
1002:     if (jac->nodal_coarsening) {
1003:       PetscViewerASCIIPrintf(viewer,"    Using nodal coarsening (with HYPRE_BOOMERAMGSetNodal() %D\n",jac->nodal_coarsening);
1004:     }
1005:     if (jac->vec_interp_variant) {
1006:       PetscViewerASCIIPrintf(viewer,"    HYPRE_BoomerAMGSetInterpVecVariant() %D\n",jac->vec_interp_variant);
1007:       PetscViewerASCIIPrintf(viewer,"    HYPRE_BoomerAMGSetInterpVecQMax() %D\n",jac->vec_interp_qmax);
1008:       PetscViewerASCIIPrintf(viewer,"    HYPRE_BoomerAMGSetSmoothInterpVectors() %d\n",jac->vec_interp_smooth);
1009:     }
1010:     if (jac->nodal_relax) {
1011:       PetscViewerASCIIPrintf(viewer,"    Using nodal relaxation via Schwarz smoothing on levels %D\n",jac->nodal_relax_levels);
1012:     }
1013:   }
1014:   return(0);
1015: }

1017: /* --------------------------------------------------------------------------------------------*/
1018: static PetscErrorCode PCSetFromOptions_HYPRE_ParaSails(PetscOptionItems *PetscOptionsObject,PC pc)
1019: {
1020:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1022:   PetscInt       indx;
1023:   PetscBool      flag;
1024:   const char     *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"};

1027:   PetscOptionsHead(PetscOptionsObject,"HYPRE ParaSails Options");
1028:   PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);
1029:   PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshhold,&jac->threshhold,&flag);
1030:   if (flag) PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));

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

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

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

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

1044:   PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,ALEN(symtlist),symtlist[0],&indx,&flag);
1045:   if (flag) {
1046:     jac->symt = indx;
1047:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
1048:   }

1050:   PetscOptionsTail();
1051:   return(0);
1052: }

1054: static PetscErrorCode PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
1055: {
1056:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1058:   PetscBool      iascii;
1059:   const char     *symt = 0;;

1062:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1063:   if (iascii) {
1064:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails preconditioning\n");
1065:     PetscViewerASCIIPrintf(viewer,"    nlevels %d\n",jac->nlevels);
1066:     PetscViewerASCIIPrintf(viewer,"    threshold %g\n",(double)jac->threshhold);
1067:     PetscViewerASCIIPrintf(viewer,"    filter %g\n",(double)jac->filter);
1068:     PetscViewerASCIIPrintf(viewer,"    load balance %g\n",(double)jac->loadbal);
1069:     PetscViewerASCIIPrintf(viewer,"    reuse nonzero structure %s\n",PetscBools[jac->ruse]);
1070:     PetscViewerASCIIPrintf(viewer,"    print info to screen %s\n",PetscBools[jac->logging]);
1071:     if (!jac->symt) symt = "nonsymmetric matrix and preconditioner";
1072:     else if (jac->symt == 1) symt = "SPD matrix and preconditioner";
1073:     else if (jac->symt == 2) symt = "nonsymmetric matrix but SPD preconditioner";
1074:     else SETERRQ1(PetscObjectComm((PetscObject)pc),PETSC_ERR_ARG_WRONG,"Unknown HYPRE ParaSails symmetric option %d",jac->symt);
1075:     PetscViewerASCIIPrintf(viewer,"    %s\n",symt);
1076:   }
1077:   return(0);
1078: }
1079: /* --------------------------------------------------------------------------------------------*/
1080: static PetscErrorCode PCSetFromOptions_HYPRE_AMS(PetscOptionItems *PetscOptionsObject,PC pc)
1081: {
1082:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1084:   PetscInt       n;
1085:   PetscBool      flag,flag2,flag3,flag4;

1088:   PetscOptionsHead(PetscOptionsObject,"HYPRE AMS Options");
1089:   PetscOptionsInt("-pc_hypre_ams_print_level","Debugging output level for AMS","None",jac->as_print,&jac->as_print,&flag);
1090:   if (flag) PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
1091:   PetscOptionsInt("-pc_hypre_ams_max_iter","Maximum number of AMS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
1092:   if (flag) PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1093:   PetscOptionsInt("-pc_hypre_ams_cycle_type","Cycle type for AMS multigrid","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag);
1094:   if (flag) PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1095:   PetscOptionsReal("-pc_hypre_ams_tol","Error tolerance for AMS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
1096:   if (flag) PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
1097:   PetscOptionsInt("-pc_hypre_ams_relax_type","Relaxation type for AMS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
1098:   PetscOptionsInt("-pc_hypre_ams_relax_times","Number of relaxation steps for AMS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
1099:   PetscOptionsReal("-pc_hypre_ams_relax_weight","Relaxation weight for AMS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
1100:   PetscOptionsReal("-pc_hypre_ams_omega","SSOR coefficient for AMS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
1101:   if (flag || flag2 || flag3 || flag4) {
1102:     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1103:                                                                       jac->as_relax_times,
1104:                                                                       jac->as_relax_weight,
1105:                                                                       jac->as_omega));
1106:   }
1107:   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);
1108:   n = 5;
1109:   PetscOptionsIntArray("-pc_hypre_ams_amg_alpha_options","AMG options for vector Poisson","None",jac->as_amg_alpha_opts,&n,&flag2);
1110:   if (flag || flag2) {
1111:     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1112:                                                                      jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1113:                                                                      jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1114:                                                                      jac->as_amg_alpha_theta,
1115:                                                                      jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1116:                                                                      jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1117:   }
1118:   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);
1119:   n = 5;
1120:   PetscOptionsIntArray("-pc_hypre_ams_amg_beta_options","AMG options for scalar Poisson solver","None",jac->as_amg_beta_opts,&n,&flag2);
1121:   if (flag || flag2) {
1122:     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1123:                                                                     jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1124:                                                                     jac->as_amg_beta_opts[2],       /* AMG relax_type */
1125:                                                                     jac->as_amg_beta_theta,
1126:                                                                     jac->as_amg_beta_opts[3],       /* AMG interp_type */
1127:                                                                     jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1128:   }
1129:   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);
1130:   if (flag) { /* override HYPRE's default only if the options is used */
1131:     PetscStackCallStandard(HYPRE_AMSSetProjectionFrequency,(jac->hsolver,jac->ams_proj_freq));
1132:   }
1133:   PetscOptionsTail();
1134:   return(0);
1135: }

1137: static PetscErrorCode PCView_HYPRE_AMS(PC pc,PetscViewer viewer)
1138: {
1139:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1141:   PetscBool      iascii;

1144:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1145:   if (iascii) {
1146:     PetscViewerASCIIPrintf(viewer,"  HYPRE AMS preconditioning\n");
1147:     PetscViewerASCIIPrintf(viewer,"    subspace iterations per application %d\n",jac->as_max_iter);
1148:     PetscViewerASCIIPrintf(viewer,"    subspace cycle type %d\n",jac->ams_cycle_type);
1149:     PetscViewerASCIIPrintf(viewer,"    subspace iteration tolerance %g\n",jac->as_tol);
1150:     PetscViewerASCIIPrintf(viewer,"    smoother type %d\n",jac->as_relax_type);
1151:     PetscViewerASCIIPrintf(viewer,"    number of smoothing steps %d\n",jac->as_relax_times);
1152:     PetscViewerASCIIPrintf(viewer,"    smoother weight %g\n",jac->as_relax_weight);
1153:     PetscViewerASCIIPrintf(viewer,"    smoother omega %g\n",jac->as_omega);
1154:     if (jac->alpha_Poisson) {
1155:       PetscViewerASCIIPrintf(viewer,"    vector Poisson solver (passed in by user)\n");
1156:     } else {
1157:       PetscViewerASCIIPrintf(viewer,"    vector Poisson solver (computed) \n");
1158:     }
1159:     PetscViewerASCIIPrintf(viewer,"        boomerAMG coarsening type %d\n",jac->as_amg_alpha_opts[0]);
1160:     PetscViewerASCIIPrintf(viewer,"        boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
1161:     PetscViewerASCIIPrintf(viewer,"        boomerAMG relaxation type %d\n",jac->as_amg_alpha_opts[2]);
1162:     PetscViewerASCIIPrintf(viewer,"        boomerAMG interpolation type %d\n",jac->as_amg_alpha_opts[3]);
1163:     PetscViewerASCIIPrintf(viewer,"        boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
1164:     PetscViewerASCIIPrintf(viewer,"        boomerAMG strength threshold %g\n",jac->as_amg_alpha_theta);
1165:     if (!jac->ams_beta_is_zero) {
1166:       if (jac->beta_Poisson) {
1167:         PetscViewerASCIIPrintf(viewer,"    scalar Poisson solver (passed in by user)\n");
1168:       } else {
1169:         PetscViewerASCIIPrintf(viewer,"    scalar Poisson solver (computed) \n");
1170:       }
1171:       PetscViewerASCIIPrintf(viewer,"        boomerAMG coarsening type %d\n",jac->as_amg_beta_opts[0]);
1172:       PetscViewerASCIIPrintf(viewer,"        boomerAMG levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
1173:       PetscViewerASCIIPrintf(viewer,"        boomerAMG relaxation type %d\n",jac->as_amg_beta_opts[2]);
1174:       PetscViewerASCIIPrintf(viewer,"        boomerAMG interpolation type %d\n",jac->as_amg_beta_opts[3]);
1175:       PetscViewerASCIIPrintf(viewer,"        boomerAMG max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
1176:       PetscViewerASCIIPrintf(viewer,"        boomerAMG strength threshold %g\n",jac->as_amg_beta_theta);
1177:       if (jac->ams_beta_is_zero_part) {
1178:         PetscViewerASCIIPrintf(viewer,"        compatible subspace projection frequency %d (-1 HYPRE uses default)\n",jac->ams_proj_freq);
1179:       }
1180:     } else {
1181:       PetscViewerASCIIPrintf(viewer,"    scalar Poisson solver not used (zero-conductivity everywhere) \n");
1182:     }
1183:   }
1184:   return(0);
1185: }

1187: static PetscErrorCode PCSetFromOptions_HYPRE_ADS(PetscOptionItems *PetscOptionsObject,PC pc)
1188: {
1189:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1191:   PetscInt       n;
1192:   PetscBool      flag,flag2,flag3,flag4;

1195:   PetscOptionsHead(PetscOptionsObject,"HYPRE ADS Options");
1196:   PetscOptionsInt("-pc_hypre_ads_print_level","Debugging output level for ADS","None",jac->as_print,&jac->as_print,&flag);
1197:   if (flag) PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
1198:   PetscOptionsInt("-pc_hypre_ads_max_iter","Maximum number of ADS multigrid iterations within PCApply","None",jac->as_max_iter,&jac->as_max_iter,&flag);
1199:   if (flag) PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1200:   PetscOptionsInt("-pc_hypre_ads_cycle_type","Cycle type for ADS multigrid","None",jac->ads_cycle_type,&jac->ads_cycle_type,&flag);
1201:   if (flag) PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ads_cycle_type));
1202:   PetscOptionsReal("-pc_hypre_ads_tol","Error tolerance for ADS multigrid","None",jac->as_tol,&jac->as_tol,&flag);
1203:   if (flag) PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
1204:   PetscOptionsInt("-pc_hypre_ads_relax_type","Relaxation type for ADS smoother","None",jac->as_relax_type,&jac->as_relax_type,&flag);
1205:   PetscOptionsInt("-pc_hypre_ads_relax_times","Number of relaxation steps for ADS smoother","None",jac->as_relax_times,&jac->as_relax_times,&flag2);
1206:   PetscOptionsReal("-pc_hypre_ads_relax_weight","Relaxation weight for ADS smoother","None",jac->as_relax_weight,&jac->as_relax_weight,&flag3);
1207:   PetscOptionsReal("-pc_hypre_ads_omega","SSOR coefficient for ADS smoother","None",jac->as_omega,&jac->as_omega,&flag4);
1208:   if (flag || flag2 || flag3 || flag4) {
1209:     PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1210:                                                                       jac->as_relax_times,
1211:                                                                       jac->as_relax_weight,
1212:                                                                       jac->as_omega));
1213:   }
1214:   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);
1215:   n = 5;
1216:   PetscOptionsIntArray("-pc_hypre_ads_ams_options","AMG options for AMS solver inside ADS","None",jac->as_amg_alpha_opts,&n,&flag2);
1217:   PetscOptionsInt("-pc_hypre_ads_ams_cycle_type","Cycle type for AMS solver inside ADS","None",jac->ams_cycle_type,&jac->ams_cycle_type,&flag3);
1218:   if (flag || flag2 || flag3) {
1219:     PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type,             /* AMS cycle type */
1220:                                                                 jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1221:                                                                 jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1222:                                                                 jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1223:                                                                 jac->as_amg_alpha_theta,
1224:                                                                 jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1225:                                                                 jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1226:   }
1227:   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);
1228:   n = 5;
1229:   PetscOptionsIntArray("-pc_hypre_ads_amg_options","AMG options for vector AMG solver inside ADS","None",jac->as_amg_beta_opts,&n,&flag2);
1230:   if (flag || flag2) {
1231:     PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1232:                                                                 jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1233:                                                                 jac->as_amg_beta_opts[2],       /* AMG relax_type */
1234:                                                                 jac->as_amg_beta_theta,
1235:                                                                 jac->as_amg_beta_opts[3],       /* AMG interp_type */
1236:                                                                 jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1237:   }
1238:   PetscOptionsTail();
1239:   return(0);
1240: }

1242: static PetscErrorCode PCView_HYPRE_ADS(PC pc,PetscViewer viewer)
1243: {
1244:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1246:   PetscBool      iascii;

1249:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
1250:   if (iascii) {
1251:     PetscViewerASCIIPrintf(viewer,"  HYPRE ADS preconditioning\n");
1252:     PetscViewerASCIIPrintf(viewer,"    subspace iterations per application %d\n",jac->as_max_iter);
1253:     PetscViewerASCIIPrintf(viewer,"    subspace cycle type %d\n",jac->ads_cycle_type);
1254:     PetscViewerASCIIPrintf(viewer,"    subspace iteration tolerance %g\n",jac->as_tol);
1255:     PetscViewerASCIIPrintf(viewer,"    smoother type %d\n",jac->as_relax_type);
1256:     PetscViewerASCIIPrintf(viewer,"    number of smoothing steps %d\n",jac->as_relax_times);
1257:     PetscViewerASCIIPrintf(viewer,"    smoother weight %g\n",jac->as_relax_weight);
1258:     PetscViewerASCIIPrintf(viewer,"    smoother omega %g\n",jac->as_omega);
1259:     PetscViewerASCIIPrintf(viewer,"    AMS solver using boomerAMG\n");
1260:     PetscViewerASCIIPrintf(viewer,"        subspace cycle type %d\n",jac->ams_cycle_type);
1261:     PetscViewerASCIIPrintf(viewer,"        coarsening type %d\n",jac->as_amg_alpha_opts[0]);
1262:     PetscViewerASCIIPrintf(viewer,"        levels of aggressive coarsening %d\n",jac->as_amg_alpha_opts[1]);
1263:     PetscViewerASCIIPrintf(viewer,"        relaxation type %d\n",jac->as_amg_alpha_opts[2]);
1264:     PetscViewerASCIIPrintf(viewer,"        interpolation type %d\n",jac->as_amg_alpha_opts[3]);
1265:     PetscViewerASCIIPrintf(viewer,"        max nonzero elements in interpolation rows %d\n",jac->as_amg_alpha_opts[4]);
1266:     PetscViewerASCIIPrintf(viewer,"        strength threshold %g\n",jac->as_amg_alpha_theta);
1267:     PetscViewerASCIIPrintf(viewer,"    vector Poisson solver using boomerAMG\n");
1268:     PetscViewerASCIIPrintf(viewer,"        coarsening type %d\n",jac->as_amg_beta_opts[0]);
1269:     PetscViewerASCIIPrintf(viewer,"        levels of aggressive coarsening %d\n",jac->as_amg_beta_opts[1]);
1270:     PetscViewerASCIIPrintf(viewer,"        relaxation type %d\n",jac->as_amg_beta_opts[2]);
1271:     PetscViewerASCIIPrintf(viewer,"        interpolation type %d\n",jac->as_amg_beta_opts[3]);
1272:     PetscViewerASCIIPrintf(viewer,"        max nonzero elements in interpolation rows %d\n",jac->as_amg_beta_opts[4]);
1273:     PetscViewerASCIIPrintf(viewer,"        strength threshold %g\n",jac->as_amg_beta_theta);
1274:   }
1275:   return(0);
1276: }

1278: static PetscErrorCode PCHYPRESetDiscreteGradient_HYPRE(PC pc, Mat G)
1279: {
1280:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1281:   PetscBool      ishypre;

1285:   PetscObjectTypeCompare((PetscObject)G,MATHYPRE,&ishypre);
1286:   if (ishypre) {
1287:     PetscObjectReference((PetscObject)G);
1288:     MatDestroy(&jac->G);
1289:     jac->G = G;
1290:   } else {
1291:     MatDestroy(&jac->G);
1292:     MatConvert(G,MATHYPRE,MAT_INITIAL_MATRIX,&jac->G);
1293:   }
1294:   return(0);
1295: }

1297: /*@
1298:  PCHYPRESetDiscreteGradient - Set discrete gradient matrix

1300:    Collective on PC

1302:    Input Parameters:
1303: +  pc - the preconditioning context
1304: -  G - the discrete gradient

1306:    Level: intermediate

1308:    Notes:
1309:     G should have as many rows as the number of edges and as many columns as the number of vertices in the mesh
1310:           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

1312: .seealso:
1313: @*/
1314: PetscErrorCode PCHYPRESetDiscreteGradient(PC pc, Mat G)
1315: {

1322:   PetscTryMethod(pc,"PCHYPRESetDiscreteGradient_C",(PC,Mat),(pc,G));
1323:   return(0);
1324: }

1326: static PetscErrorCode PCHYPRESetDiscreteCurl_HYPRE(PC pc, Mat C)
1327: {
1328:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1329:   PetscBool      ishypre;

1333:   PetscObjectTypeCompare((PetscObject)C,MATHYPRE,&ishypre);
1334:   if (ishypre) {
1335:     PetscObjectReference((PetscObject)C);
1336:     MatDestroy(&jac->C);
1337:     jac->C = C;
1338:   } else {
1339:     MatDestroy(&jac->C);
1340:     MatConvert(C,MATHYPRE,MAT_INITIAL_MATRIX,&jac->C);
1341:   }
1342:   return(0);
1343: }

1345: /*@
1346:  PCHYPRESetDiscreteCurl - Set discrete curl matrix

1348:    Collective on PC

1350:    Input Parameters:
1351: +  pc - the preconditioning context
1352: -  C - the discrete curl

1354:    Level: intermediate

1356:    Notes:
1357:     C should have as many rows as the number of faces and as many columns as the number of edges in the mesh
1358:           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

1360: .seealso:
1361: @*/
1362: PetscErrorCode PCHYPRESetDiscreteCurl(PC pc, Mat C)
1363: {

1370:   PetscTryMethod(pc,"PCHYPRESetDiscreteCurl_C",(PC,Mat),(pc,C));
1371:   return(0);
1372: }

1374: static PetscErrorCode PCHYPRESetInterpolations_HYPRE(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[])
1375: {
1376:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1377:   PetscBool      ishypre;
1379:   PetscInt       i;

1382:   MatDestroy(&jac->RT_PiFull);
1383:   MatDestroy(&jac->ND_PiFull);
1384:   for (i=0;i<3;++i) {
1385:     MatDestroy(&jac->RT_Pi[i]);
1386:     MatDestroy(&jac->ND_Pi[i]);
1387:   }

1389:   jac->dim = dim;
1390:   if (RT_PiFull) {
1391:     PetscObjectTypeCompare((PetscObject)RT_PiFull,MATHYPRE,&ishypre);
1392:     if (ishypre) {
1393:       PetscObjectReference((PetscObject)RT_PiFull);
1394:       jac->RT_PiFull = RT_PiFull;
1395:     } else {
1396:       MatConvert(RT_PiFull,MATHYPRE,MAT_INITIAL_MATRIX,&jac->RT_PiFull);
1397:     }
1398:   }
1399:   if (RT_Pi) {
1400:     for (i=0;i<dim;++i) {
1401:       if (RT_Pi[i]) {
1402:         PetscObjectTypeCompare((PetscObject)RT_Pi[i],MATHYPRE,&ishypre);
1403:         if (ishypre) {
1404:           PetscObjectReference((PetscObject)RT_Pi[i]);
1405:           jac->RT_Pi[i] = RT_Pi[i];
1406:         } else {
1407:           MatConvert(RT_Pi[i],MATHYPRE,MAT_INITIAL_MATRIX,&jac->RT_Pi[i]);
1408:         }
1409:       }
1410:     }
1411:   }
1412:   if (ND_PiFull) {
1413:     PetscObjectTypeCompare((PetscObject)ND_PiFull,MATHYPRE,&ishypre);
1414:     if (ishypre) {
1415:       PetscObjectReference((PetscObject)ND_PiFull);
1416:       jac->ND_PiFull = ND_PiFull;
1417:     } else {
1418:       MatConvert(ND_PiFull,MATHYPRE,MAT_INITIAL_MATRIX,&jac->ND_PiFull);
1419:     }
1420:   }
1421:   if (ND_Pi) {
1422:     for (i=0;i<dim;++i) {
1423:       if (ND_Pi[i]) {
1424:         PetscObjectTypeCompare((PetscObject)ND_Pi[i],MATHYPRE,&ishypre);
1425:         if (ishypre) {
1426:           PetscObjectReference((PetscObject)ND_Pi[i]);
1427:           jac->ND_Pi[i] = ND_Pi[i];
1428:         } else {
1429:           MatConvert(ND_Pi[i],MATHYPRE,MAT_INITIAL_MATRIX,&jac->ND_Pi[i]);
1430:         }
1431:       }
1432:     }
1433:   }

1435:   return(0);
1436: }

1438: /*@
1439:  PCHYPRESetInterpolations - Set interpolation matrices for AMS/ADS preconditioner

1441:    Collective on PC

1443:    Input Parameters:
1444: +  pc - the preconditioning context
1445: -  dim - the dimension of the problem, only used in AMS
1446: -  RT_PiFull - Raviart-Thomas interpolation matrix
1447: -  RT_Pi - x/y/z component of Raviart-Thomas interpolation matrix
1448: -  ND_PiFull - Nedelec interpolation matrix
1449: -  ND_Pi - x/y/z component of Nedelec interpolation matrix

1451:    Notes:
1452:     For AMS, only Nedelec interpolation matrices are needed, the Raviart-Thomas interpolation matrices can be set to NULL.
1453:           For ADS, both type of interpolation matrices are needed.
1454:    Level: intermediate

1456: .seealso:
1457: @*/
1458: PetscErrorCode PCHYPRESetInterpolations(PC pc, PetscInt dim, Mat RT_PiFull, Mat RT_Pi[], Mat ND_PiFull, Mat ND_Pi[])
1459: {
1461:   PetscInt       i;

1465:   if (RT_PiFull) {
1468:   }
1469:   if (RT_Pi) {
1471:     for (i=0;i<dim;++i) {
1472:       if (RT_Pi[i]) {
1475:       }
1476:     }
1477:   }
1478:   if (ND_PiFull) {
1481:   }
1482:   if (ND_Pi) {
1484:     for (i=0;i<dim;++i) {
1485:       if (ND_Pi[i]) {
1488:       }
1489:     }
1490:   }
1491:   PetscTryMethod(pc,"PCHYPRESetInterpolations_C",(PC,PetscInt,Mat,Mat[],Mat,Mat[]),(pc,dim,RT_PiFull,RT_Pi,ND_PiFull,ND_Pi));
1492:   return(0);
1493: }

1495: static PetscErrorCode PCHYPRESetPoissonMatrix_HYPRE(PC pc, Mat A, PetscBool isalpha)
1496: {
1497:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1498:   PetscBool      ishypre;

1502:   PetscObjectTypeCompare((PetscObject)A,MATHYPRE,&ishypre);
1503:   if (ishypre) {
1504:     if (isalpha) {
1505:       PetscObjectReference((PetscObject)A);
1506:       MatDestroy(&jac->alpha_Poisson);
1507:       jac->alpha_Poisson = A;
1508:     } else {
1509:       if (A) {
1510:         PetscObjectReference((PetscObject)A);
1511:       } else {
1512:         jac->ams_beta_is_zero = PETSC_TRUE;
1513:       }
1514:       MatDestroy(&jac->beta_Poisson);
1515:       jac->beta_Poisson = A;
1516:     }
1517:   } else {
1518:     if (isalpha) {
1519:       MatDestroy(&jac->alpha_Poisson);
1520:       MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&jac->alpha_Poisson);
1521:     } else {
1522:       if (A) {
1523:         MatDestroy(&jac->beta_Poisson);
1524:         MatConvert(A,MATHYPRE,MAT_INITIAL_MATRIX,&jac->beta_Poisson);
1525:       } else {
1526:         MatDestroy(&jac->beta_Poisson);
1527:         jac->ams_beta_is_zero = PETSC_TRUE;
1528:       }
1529:     }
1530:   }
1531:   return(0);
1532: }

1534: /*@
1535:  PCHYPRESetAlphaPoissonMatrix - Set vector Poisson matrix

1537:    Collective on PC

1539:    Input Parameters:
1540: +  pc - the preconditioning context
1541: -  A - the matrix

1543:    Level: intermediate

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

1548: .seealso:
1549: @*/
1550: PetscErrorCode PCHYPRESetAlphaPoissonMatrix(PC pc, Mat A)
1551: {

1558:   PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_TRUE));
1559:   return(0);
1560: }

1562: /*@
1563:  PCHYPRESetBetaPoissonMatrix - Set Poisson matrix

1565:    Collective on PC

1567:    Input Parameters:
1568: +  pc - the preconditioning context
1569: -  A - the matrix

1571:    Level: intermediate

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

1577: .seealso:
1578: @*/
1579: PetscErrorCode PCHYPRESetBetaPoissonMatrix(PC pc, Mat A)
1580: {

1585:   if (A) {
1588:   }
1589:   PetscTryMethod(pc,"PCHYPRESetPoissonMatrix_C",(PC,Mat,PetscBool),(pc,A,PETSC_FALSE));
1590:   return(0);
1591: }

1593: static PetscErrorCode PCHYPRESetEdgeConstantVectors_HYPRE(PC pc,Vec ozz, Vec zoz, Vec zzo)
1594: {
1595:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
1596:   PetscErrorCode     ierr;

1599:   /* throw away any vector if already set */
1600:   if (jac->constants[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[0]));
1601:   if (jac->constants[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[1]));
1602:   if (jac->constants[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->constants[2]));
1603:   jac->constants[0] = NULL;
1604:   jac->constants[1] = NULL;
1605:   jac->constants[2] = NULL;
1606:   VecHYPRE_IJVectorCreate(ozz,&jac->constants[0]);
1607:   VecHYPRE_IJVectorCopy(ozz,jac->constants[0]);
1608:   VecHYPRE_IJVectorCreate(zoz,&jac->constants[1]);
1609:   VecHYPRE_IJVectorCopy(zoz,jac->constants[1]);
1610:   jac->dim = 2;
1611:   if (zzo) {
1612:     VecHYPRE_IJVectorCreate(zzo,&jac->constants[2]);
1613:     VecHYPRE_IJVectorCopy(zzo,jac->constants[2]);
1614:     jac->dim++;
1615:   }
1616:   return(0);
1617: }

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

1622:    Collective on PC

1624:    Input Parameters:
1625: +  pc - the preconditioning context
1626: -  ozz - vector representing (1,0,0) (or (1,0) in 2D)
1627: -  zoz - vector representing (0,1,0) (or (0,1) in 2D)
1628: -  zzo - vector representing (0,0,1) (use NULL in 2D)

1630:    Level: intermediate

1632:    Notes:

1634: .seealso:
1635: @*/
1636: PetscErrorCode PCHYPRESetEdgeConstantVectors(PC pc, Vec ozz, Vec zoz, Vec zzo)
1637: {

1648:   PetscTryMethod(pc,"PCHYPRESetEdgeConstantVectors_C",(PC,Vec,Vec,Vec),(pc,ozz,zoz,zzo));
1649:   return(0);
1650: }

1652: static PetscErrorCode PCSetCoordinates_HYPRE(PC pc, PetscInt dim, PetscInt nloc, PetscReal *coords)
1653: {
1654:   PC_HYPRE        *jac = (PC_HYPRE*)pc->data;
1655:   Vec             tv;
1656:   PetscInt        i;
1657:   PetscErrorCode  ierr;

1660:   /* throw away any coordinate vector if already set */
1661:   if (jac->coords[0]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[0]));
1662:   if (jac->coords[1]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[1]));
1663:   if (jac->coords[2]) PetscStackCallStandard(HYPRE_IJVectorDestroy,(jac->coords[2]));
1664:   jac->dim = dim;

1666:   /* compute IJ vector for coordinates */
1667:   VecCreate(PetscObjectComm((PetscObject)pc),&tv);
1668:   VecSetType(tv,VECSTANDARD);
1669:   VecSetSizes(tv,nloc,PETSC_DECIDE);
1670:   for (i=0;i<dim;i++) {
1671:     PetscScalar *array;
1672:     PetscInt    j;

1674:     VecHYPRE_IJVectorCreate(tv,&jac->coords[i]);
1675:     VecGetArray(tv,&array);
1676:     for (j=0;j<nloc;j++) {
1677:       array[j] = coords[j*dim+i];
1678:     }
1679:     PetscStackCallStandard(HYPRE_IJVectorSetValues,(jac->coords[i],nloc,NULL,(HYPRE_Complex*)array));
1680:     PetscStackCallStandard(HYPRE_IJVectorAssemble,(jac->coords[i]));
1681:     VecRestoreArray(tv,&array);
1682:   }
1683:   VecDestroy(&tv);
1684:   return(0);
1685: }

1687: /* ---------------------------------------------------------------------------------*/

1689: static PetscErrorCode  PCHYPREGetType_HYPRE(PC pc,const char *name[])
1690: {
1691:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;

1694:   *name = jac->hypre_type;
1695:   return(0);
1696: }

1698: static PetscErrorCode  PCHYPRESetType_HYPRE(PC pc,const char name[])
1699: {
1700:   PC_HYPRE       *jac = (PC_HYPRE*)pc->data;
1702:   PetscBool      flag;

1705:   if (jac->hypre_type) {
1706:     PetscStrcmp(jac->hypre_type,name,&flag);
1707:     if (!flag) SETERRQ(PetscObjectComm((PetscObject)pc),PETSC_ERR_ORDER,"Cannot reset the HYPRE preconditioner type once it has been set");
1708:     return(0);
1709:   } else {
1710:     PetscStrallocpy(name, &jac->hypre_type);
1711:   }

1713:   jac->maxiter         = PETSC_DEFAULT;
1714:   jac->tol             = PETSC_DEFAULT;
1715:   jac->printstatistics = PetscLogPrintInfo;

1717:   PetscStrcmp("pilut",jac->hypre_type,&flag);
1718:   if (flag) {
1719:     MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1720:     PetscStackCallStandard(HYPRE_ParCSRPilutCreate,(jac->comm_hypre,&jac->hsolver));
1721:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
1722:     pc->ops->view           = PCView_HYPRE_Pilut;
1723:     jac->destroy            = HYPRE_ParCSRPilutDestroy;
1724:     jac->setup              = HYPRE_ParCSRPilutSetup;
1725:     jac->solve              = HYPRE_ParCSRPilutSolve;
1726:     jac->factorrowsize      = PETSC_DEFAULT;
1727:     return(0);
1728:   }
1729:   PetscStrcmp("euclid",jac->hypre_type,&flag);
1730:   if (flag) {
1731:     MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1732:     PetscStackCallStandard(HYPRE_EuclidCreate,(jac->comm_hypre,&jac->hsolver));
1733:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
1734:     pc->ops->view           = PCView_HYPRE_Euclid;
1735:     jac->destroy            = HYPRE_EuclidDestroy;
1736:     jac->setup              = HYPRE_EuclidSetup;
1737:     jac->solve              = HYPRE_EuclidSolve;
1738:     jac->factorrowsize      = PETSC_DEFAULT;
1739:     jac->eu_level           = PETSC_DEFAULT; /* default */
1740:     return(0);
1741:   }
1742:   PetscStrcmp("parasails",jac->hypre_type,&flag);
1743:   if (flag) {
1744:     MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(jac->comm_hypre));
1745:     PetscStackCallStandard(HYPRE_ParaSailsCreate,(jac->comm_hypre,&jac->hsolver));
1746:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
1747:     pc->ops->view           = PCView_HYPRE_ParaSails;
1748:     jac->destroy            = HYPRE_ParaSailsDestroy;
1749:     jac->setup              = HYPRE_ParaSailsSetup;
1750:     jac->solve              = HYPRE_ParaSailsSolve;
1751:     /* initialize */
1752:     jac->nlevels    = 1;
1753:     jac->threshhold = .1;
1754:     jac->filter     = .1;
1755:     jac->loadbal    = 0;
1756:     if (PetscLogPrintInfo) jac->logging = (int) PETSC_TRUE;
1757:     else jac->logging = (int) PETSC_FALSE;

1759:     jac->ruse = (int) PETSC_FALSE;
1760:     jac->symt = 0;
1761:     PetscStackCallStandard(HYPRE_ParaSailsSetParams,(jac->hsolver,jac->threshhold,jac->nlevels));
1762:     PetscStackCallStandard(HYPRE_ParaSailsSetFilter,(jac->hsolver,jac->filter));
1763:     PetscStackCallStandard(HYPRE_ParaSailsSetLoadbal,(jac->hsolver,jac->loadbal));
1764:     PetscStackCallStandard(HYPRE_ParaSailsSetLogging,(jac->hsolver,jac->logging));
1765:     PetscStackCallStandard(HYPRE_ParaSailsSetReuse,(jac->hsolver,jac->ruse));
1766:     PetscStackCallStandard(HYPRE_ParaSailsSetSym,(jac->hsolver,jac->symt));
1767:     return(0);
1768:   }
1769:   PetscStrcmp("boomeramg",jac->hypre_type,&flag);
1770:   if (flag) {
1771:     HYPRE_BoomerAMGCreate(&jac->hsolver);
1772:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_BoomerAMG;
1773:     pc->ops->view            = PCView_HYPRE_BoomerAMG;
1774:     pc->ops->applytranspose  = PCApplyTranspose_HYPRE_BoomerAMG;
1775:     pc->ops->applyrichardson = PCApplyRichardson_HYPRE_BoomerAMG;
1776:     PetscObjectComposeFunction((PetscObject)pc,"PCGetInterpolations_C",PCGetInterpolations_BoomerAMG);
1777:     PetscObjectComposeFunction((PetscObject)pc,"PCGetCoarseOperators_C",PCGetCoarseOperators_BoomerAMG);
1778:     jac->destroy             = HYPRE_BoomerAMGDestroy;
1779:     jac->setup               = HYPRE_BoomerAMGSetup;
1780:     jac->solve               = HYPRE_BoomerAMGSolve;
1781:     jac->applyrichardson     = PETSC_FALSE;
1782:     /* these defaults match the hypre defaults */
1783:     jac->cycletype        = 1;
1784:     jac->maxlevels        = 25;
1785:     jac->maxiter          = 1;
1786:     jac->tol              = 0.0; /* tolerance of zero indicates use as preconditioner (suppresses convergence errors) */
1787:     jac->truncfactor      = 0.0;
1788:     jac->strongthreshold  = .25;
1789:     jac->maxrowsum        = .9;
1790:     jac->coarsentype      = 6;
1791:     jac->measuretype      = 0;
1792:     jac->gridsweeps[0]    = jac->gridsweeps[1] = jac->gridsweeps[2] = 1;
1793:     jac->smoothtype       = -1; /* Not set by default */
1794:     jac->smoothnumlevels  = 25;
1795:     jac->eu_level         = 0;
1796:     jac->eu_droptolerance = 0;
1797:     jac->eu_bj            = 0;
1798:     jac->relaxtype[0]     = jac->relaxtype[1] = 6; /* Defaults to SYMMETRIC since in PETSc we are using a a PC - most likely with CG */
1799:     jac->relaxtype[2]     = 9; /*G.E. */
1800:     jac->relaxweight      = 1.0;
1801:     jac->outerrelaxweight = 1.0;
1802:     jac->relaxorder       = 1;
1803:     jac->interptype       = 0;
1804:     jac->agg_nl           = 0;
1805:     jac->pmax             = 0;
1806:     jac->truncfactor      = 0.0;
1807:     jac->agg_num_paths    = 1;

1809:     jac->nodal_coarsening      = 0;
1810:     jac->nodal_coarsening_diag = 0;
1811:     jac->vec_interp_variant    = 0;
1812:     jac->vec_interp_qmax       = 0;
1813:     jac->vec_interp_smooth     = PETSC_FALSE;
1814:     jac->interp_refine         = 0;
1815:     jac->nodal_relax           = PETSC_FALSE;
1816:     jac->nodal_relax_levels    = 1;
1817:     PetscStackCallStandard(HYPRE_BoomerAMGSetCycleType,(jac->hsolver,jac->cycletype));
1818:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxLevels,(jac->hsolver,jac->maxlevels));
1819:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxIter,(jac->hsolver,jac->maxiter));
1820:     PetscStackCallStandard(HYPRE_BoomerAMGSetTol,(jac->hsolver,jac->tol));
1821:     PetscStackCallStandard(HYPRE_BoomerAMGSetTruncFactor,(jac->hsolver,jac->truncfactor));
1822:     PetscStackCallStandard(HYPRE_BoomerAMGSetStrongThreshold,(jac->hsolver,jac->strongthreshold));
1823:     PetscStackCallStandard(HYPRE_BoomerAMGSetMaxRowSum,(jac->hsolver,jac->maxrowsum));
1824:     PetscStackCallStandard(HYPRE_BoomerAMGSetCoarsenType,(jac->hsolver,jac->coarsentype));
1825:     PetscStackCallStandard(HYPRE_BoomerAMGSetMeasureType,(jac->hsolver,jac->measuretype));
1826:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxOrder,(jac->hsolver, jac->relaxorder));
1827:     PetscStackCallStandard(HYPRE_BoomerAMGSetInterpType,(jac->hsolver,jac->interptype));
1828:     PetscStackCallStandard(HYPRE_BoomerAMGSetAggNumLevels,(jac->hsolver,jac->agg_nl));
1829:     PetscStackCallStandard(HYPRE_BoomerAMGSetPMaxElmts,(jac->hsolver,jac->pmax));
1830:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumPaths,(jac->hsolver,jac->agg_num_paths));
1831:     PetscStackCallStandard(HYPRE_BoomerAMGSetRelaxType,(jac->hsolver, jac->relaxtype[0]));  /*defaults coarse to 9*/
1832:     PetscStackCallStandard(HYPRE_BoomerAMGSetNumSweeps,(jac->hsolver, jac->gridsweeps[0])); /*defaults coarse to 1 */
1833:     return(0);
1834:   }
1835:   PetscStrcmp("ams",jac->hypre_type,&flag);
1836:   if (flag) {
1837:     HYPRE_AMSCreate(&jac->hsolver);
1838:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_AMS;
1839:     pc->ops->view            = PCView_HYPRE_AMS;
1840:     jac->destroy             = HYPRE_AMSDestroy;
1841:     jac->setup               = HYPRE_AMSSetup;
1842:     jac->solve               = HYPRE_AMSSolve;
1843:     jac->coords[0]           = NULL;
1844:     jac->coords[1]           = NULL;
1845:     jac->coords[2]           = NULL;
1846:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE AMS */
1847:     jac->as_print           = 0;
1848:     jac->as_max_iter        = 1; /* used as a preconditioner */
1849:     jac->as_tol             = 0.; /* used as a preconditioner */
1850:     jac->ams_cycle_type     = 13;
1851:     /* Smoothing options */
1852:     jac->as_relax_type      = 2;
1853:     jac->as_relax_times     = 1;
1854:     jac->as_relax_weight    = 1.0;
1855:     jac->as_omega           = 1.0;
1856:     /* Vector valued Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1857:     jac->as_amg_alpha_opts[0] = 10;
1858:     jac->as_amg_alpha_opts[1] = 1;
1859:     jac->as_amg_alpha_opts[2] = 6;
1860:     jac->as_amg_alpha_opts[3] = 6;
1861:     jac->as_amg_alpha_opts[4] = 4;
1862:     jac->as_amg_alpha_theta   = 0.25;
1863:     /* Scalar Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1864:     jac->as_amg_beta_opts[0] = 10;
1865:     jac->as_amg_beta_opts[1] = 1;
1866:     jac->as_amg_beta_opts[2] = 6;
1867:     jac->as_amg_beta_opts[3] = 6;
1868:     jac->as_amg_beta_opts[4] = 4;
1869:     jac->as_amg_beta_theta   = 0.25;
1870:     PetscStackCallStandard(HYPRE_AMSSetPrintLevel,(jac->hsolver,jac->as_print));
1871:     PetscStackCallStandard(HYPRE_AMSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1872:     PetscStackCallStandard(HYPRE_AMSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1873:     PetscStackCallStandard(HYPRE_AMSSetTol,(jac->hsolver,jac->as_tol));
1874:     PetscStackCallStandard(HYPRE_AMSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1875:                                                                       jac->as_relax_times,
1876:                                                                       jac->as_relax_weight,
1877:                                                                       jac->as_omega));
1878:     PetscStackCallStandard(HYPRE_AMSSetAlphaAMGOptions,(jac->hsolver,jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1879:                                                                      jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1880:                                                                      jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1881:                                                                      jac->as_amg_alpha_theta,
1882:                                                                      jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1883:                                                                      jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1884:     PetscStackCallStandard(HYPRE_AMSSetBetaAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1885:                                                                     jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1886:                                                                     jac->as_amg_beta_opts[2],       /* AMG relax_type */
1887:                                                                     jac->as_amg_beta_theta,
1888:                                                                     jac->as_amg_beta_opts[3],       /* AMG interp_type */
1889:                                                                     jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1890:     /* Zero conductivity */
1891:     jac->ams_beta_is_zero      = PETSC_FALSE;
1892:     jac->ams_beta_is_zero_part = PETSC_FALSE;
1893:     return(0);
1894:   }
1895:   PetscStrcmp("ads",jac->hypre_type,&flag);
1896:   if (flag) {
1897:     HYPRE_ADSCreate(&jac->hsolver);
1898:     pc->ops->setfromoptions  = PCSetFromOptions_HYPRE_ADS;
1899:     pc->ops->view            = PCView_HYPRE_ADS;
1900:     jac->destroy             = HYPRE_ADSDestroy;
1901:     jac->setup               = HYPRE_ADSSetup;
1902:     jac->solve               = HYPRE_ADSSolve;
1903:     jac->coords[0]           = NULL;
1904:     jac->coords[1]           = NULL;
1905:     jac->coords[2]           = NULL;
1906:     /* solver parameters: these are borrowed from mfem package, and they are not the default values from HYPRE ADS */
1907:     jac->as_print           = 0;
1908:     jac->as_max_iter        = 1; /* used as a preconditioner */
1909:     jac->as_tol             = 0.; /* used as a preconditioner */
1910:     jac->ads_cycle_type     = 13;
1911:     /* Smoothing options */
1912:     jac->as_relax_type      = 2;
1913:     jac->as_relax_times     = 1;
1914:     jac->as_relax_weight    = 1.0;
1915:     jac->as_omega           = 1.0;
1916:     /* AMS solver parameters: cycle_type, coarsen type, agg_levels, relax_type, interp_type, Pmax */
1917:     jac->ams_cycle_type       = 14;
1918:     jac->as_amg_alpha_opts[0] = 10;
1919:     jac->as_amg_alpha_opts[1] = 1;
1920:     jac->as_amg_alpha_opts[2] = 6;
1921:     jac->as_amg_alpha_opts[3] = 6;
1922:     jac->as_amg_alpha_opts[4] = 4;
1923:     jac->as_amg_alpha_theta   = 0.25;
1924:     /* Vector Poisson AMG solver parameters: coarsen type, agg_levels, relax_type, interp_type, Pmax */
1925:     jac->as_amg_beta_opts[0] = 10;
1926:     jac->as_amg_beta_opts[1] = 1;
1927:     jac->as_amg_beta_opts[2] = 6;
1928:     jac->as_amg_beta_opts[3] = 6;
1929:     jac->as_amg_beta_opts[4] = 4;
1930:     jac->as_amg_beta_theta   = 0.25;
1931:     PetscStackCallStandard(HYPRE_ADSSetPrintLevel,(jac->hsolver,jac->as_print));
1932:     PetscStackCallStandard(HYPRE_ADSSetMaxIter,(jac->hsolver,jac->as_max_iter));
1933:     PetscStackCallStandard(HYPRE_ADSSetCycleType,(jac->hsolver,jac->ams_cycle_type));
1934:     PetscStackCallStandard(HYPRE_ADSSetTol,(jac->hsolver,jac->as_tol));
1935:     PetscStackCallStandard(HYPRE_ADSSetSmoothingOptions,(jac->hsolver,jac->as_relax_type,
1936:                                                                       jac->as_relax_times,
1937:                                                                       jac->as_relax_weight,
1938:                                                                       jac->as_omega));
1939:     PetscStackCallStandard(HYPRE_ADSSetAMSOptions,(jac->hsolver,jac->ams_cycle_type,             /* AMG coarsen type */
1940:                                                                 jac->as_amg_alpha_opts[0],       /* AMG coarsen type */
1941:                                                                 jac->as_amg_alpha_opts[1],       /* AMG agg_levels */
1942:                                                                 jac->as_amg_alpha_opts[2],       /* AMG relax_type */
1943:                                                                 jac->as_amg_alpha_theta,
1944:                                                                 jac->as_amg_alpha_opts[3],       /* AMG interp_type */
1945:                                                                 jac->as_amg_alpha_opts[4]));     /* AMG Pmax */
1946:     PetscStackCallStandard(HYPRE_ADSSetAMGOptions,(jac->hsolver,jac->as_amg_beta_opts[0],       /* AMG coarsen type */
1947:                                                                 jac->as_amg_beta_opts[1],       /* AMG agg_levels */
1948:                                                                 jac->as_amg_beta_opts[2],       /* AMG relax_type */
1949:                                                                 jac->as_amg_beta_theta,
1950:                                                                 jac->as_amg_beta_opts[3],       /* AMG interp_type */
1951:                                                                 jac->as_amg_beta_opts[4]));     /* AMG Pmax */
1952:     return(0);
1953:   }
1954:   PetscFree(jac->hypre_type);

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

1961: /*
1962:     It only gets here if the HYPRE type has not been set before the call to
1963:    ...SetFromOptions() which actually is most of the time
1964: */
1965: PetscErrorCode PCSetFromOptions_HYPRE(PetscOptionItems *PetscOptionsObject,PC pc)
1966: {
1968:   PetscInt       indx;
1969:   const char     *type[] = {"euclid","pilut","parasails","boomeramg","ams","ads"};
1970:   PetscBool      flg;

1973:   PetscOptionsHead(PetscOptionsObject,"HYPRE preconditioner options");
1974:   PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"boomeramg",&indx,&flg);
1975:   if (flg) {
1976:     PCHYPRESetType_HYPRE(pc,type[indx]);
1977:   } else {
1978:     PCHYPRESetType_HYPRE(pc,"boomeramg");
1979:   }
1980:   if (pc->ops->setfromoptions) {
1981:     pc->ops->setfromoptions(PetscOptionsObject,pc);
1982:   }
1983:   PetscOptionsTail();
1984:   return(0);
1985: }

1987: /*@C
1988:      PCHYPRESetType - Sets which hypre preconditioner you wish to use

1990:    Input Parameters:
1991: +     pc - the preconditioner context
1992: -     name - either  euclid, pilut, parasails, boomeramg, ams, ads

1994:    Options Database Keys:
1995:    -pc_hypre_type - One of euclid, pilut, parasails, boomeramg, ams, ads

1997:    Level: intermediate

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

2002: @*/
2003: PetscErrorCode  PCHYPRESetType(PC pc,const char name[])
2004: {

2010:   PetscTryMethod(pc,"PCHYPRESetType_C",(PC,const char[]),(pc,name));
2011:   return(0);
2012: }

2014: /*@C
2015:      PCHYPREGetType - Gets which hypre preconditioner you are using

2017:    Input Parameter:
2018: .     pc - the preconditioner context

2020:    Output Parameter:
2021: .     name - either  euclid, pilut, parasails, boomeramg, ams, ads

2023:    Level: intermediate

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

2028: @*/
2029: PetscErrorCode  PCHYPREGetType(PC pc,const char *name[])
2030: {

2036:   PetscTryMethod(pc,"PCHYPREGetType_C",(PC,const char*[]),(pc,name));
2037:   return(0);
2038: }

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

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

2048:    Level: intermediate

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

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

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

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

2073:           MatSetNearNullSpace() - if you provide a near null space to your matrix it is ignored by hypre UNLESS you also use
2074:           the two options:
2075: +   -pc_hypre_boomeramg_nodal_coarsen <n> - where n is from 1 to 6 (see HYPRE_BOOMERAMGSetNodal())
2076: -   -pc_hypre_boomeramg_vec_interp_variant <v> where v is from 1 to 3 (see HYPRE_BoomerAMGSetInterpVecVariant())

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

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

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

2085: M*/

2087: PETSC_EXTERN PetscErrorCode PCCreate_HYPRE(PC pc)
2088: {
2089:   PC_HYPRE       *jac;

2093:   PetscNewLog(pc,&jac);

2095:   pc->data                = jac;
2096:   pc->ops->reset          = PCReset_HYPRE;
2097:   pc->ops->destroy        = PCDestroy_HYPRE;
2098:   pc->ops->setfromoptions = PCSetFromOptions_HYPRE;
2099:   pc->ops->setup          = PCSetUp_HYPRE;
2100:   pc->ops->apply          = PCApply_HYPRE;
2101:   jac->comm_hypre         = MPI_COMM_NULL;
2102:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetType_C",PCHYPRESetType_HYPRE);
2103:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPREGetType_C",PCHYPREGetType_HYPRE);
2104:   PetscObjectComposeFunction((PetscObject)pc,"PCSetCoordinates_C",PCSetCoordinates_HYPRE);
2105:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteGradient_C",PCHYPRESetDiscreteGradient_HYPRE);
2106:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetDiscreteCurl_C",PCHYPRESetDiscreteCurl_HYPRE);
2107:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetInterpolations_C",PCHYPRESetInterpolations_HYPRE);
2108:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetEdgeConstantVectors_C",PCHYPRESetEdgeConstantVectors_HYPRE);
2109:   PetscObjectComposeFunction((PetscObject)pc,"PCHYPRESetPoissonMatrix_C",PCHYPRESetPoissonMatrix_HYPRE);
2110:   return(0);
2111: }

2113: /* ---------------------------------------------------------------------------------------------------------------------------------*/

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

2119:   /* keep copy of PFMG options used so may view them */
2120:   PetscInt its;
2121:   double   tol;
2122:   PetscInt relax_type;
2123:   PetscInt rap_type;
2124:   PetscInt num_pre_relax,num_post_relax;
2125:   PetscInt max_levels;
2126: } PC_PFMG;

2128: PetscErrorCode PCDestroy_PFMG(PC pc)
2129: {
2131:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

2134:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
2135:   MPI_Comm_free(&ex->hcomm);
2136:   PetscFree(pc->data);
2137:   return(0);
2138: }

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

2143: PetscErrorCode PCView_PFMG(PC pc,PetscViewer viewer)
2144: {
2146:   PetscBool      iascii;
2147:   PC_PFMG        *ex = (PC_PFMG*) pc->data;

2150:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2151:   if (iascii) {
2152:     PetscViewerASCIIPrintf(viewer,"  HYPRE PFMG preconditioning\n");
2153:     PetscViewerASCIIPrintf(viewer,"    max iterations %d\n",ex->its);
2154:     PetscViewerASCIIPrintf(viewer,"    tolerance %g\n",ex->tol);
2155:     PetscViewerASCIIPrintf(viewer,"    relax type %s\n",PFMGRelaxType[ex->relax_type]);
2156:     PetscViewerASCIIPrintf(viewer,"    RAP type %s\n",PFMGRAPType[ex->rap_type]);
2157:     PetscViewerASCIIPrintf(viewer,"    number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
2158:     PetscViewerASCIIPrintf(viewer,"    max levels %d\n",ex->max_levels);
2159:   }
2160:   return(0);
2161: }

2163: PetscErrorCode PCSetFromOptions_PFMG(PetscOptionItems *PetscOptionsObject,PC pc)
2164: {
2166:   PC_PFMG        *ex = (PC_PFMG*) pc->data;
2167:   PetscBool      flg = PETSC_FALSE;

2170:   PetscOptionsHead(PetscOptionsObject,"PFMG options");
2171:   PetscOptionsBool("-pc_pfmg_print_statistics","Print statistics","HYPRE_StructPFMGSetPrintLevel",flg,&flg,NULL);
2172:   if (flg) {
2173:     int level=3;
2174:     PetscStackCallStandard(HYPRE_StructPFMGSetPrintLevel,(ex->hsolver,level));
2175:   }
2176:   PetscOptionsInt("-pc_pfmg_its","Number of iterations of PFMG to use as preconditioner","HYPRE_StructPFMGSetMaxIter",ex->its,&ex->its,NULL);
2177:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(ex->hsolver,ex->its));
2178:   PetscOptionsInt("-pc_pfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_StructPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
2179:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPreRelax,(ex->hsolver,ex->num_pre_relax));
2180:   PetscOptionsInt("-pc_pfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_StructPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
2181:   PetscStackCallStandard(HYPRE_StructPFMGSetNumPostRelax,(ex->hsolver,ex->num_post_relax));

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

2186:   PetscOptionsReal("-pc_pfmg_tol","Tolerance of PFMG","HYPRE_StructPFMGSetTol",ex->tol,&ex->tol,NULL);
2187:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(ex->hsolver,ex->tol));
2188:   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);
2189:   PetscStackCallStandard(HYPRE_StructPFMGSetRelaxType,(ex->hsolver, ex->relax_type));
2190:   PetscOptionsEList("-pc_pfmg_rap_type","RAP type","HYPRE_StructPFMGSetRAPType",PFMGRAPType,ALEN(PFMGRAPType),PFMGRAPType[ex->rap_type],&ex->rap_type,NULL);
2191:   PetscStackCallStandard(HYPRE_StructPFMGSetRAPType,(ex->hsolver, ex->rap_type));
2192:   PetscOptionsTail();
2193:   return(0);
2194: }

2196: PetscErrorCode PCApply_PFMG(PC pc,Vec x,Vec y)
2197: {
2198:   PetscErrorCode    ierr;
2199:   PC_PFMG           *ex = (PC_PFMG*) pc->data;
2200:   PetscScalar       *yy;
2201:   const PetscScalar *xx;
2202:   PetscInt          ilower[3],iupper[3];
2203:   HYPRE_Int         hlower[3],hupper[3];
2204:   Mat_HYPREStruct   *mx = (Mat_HYPREStruct*)(pc->pmat->data);

2207:   PetscCitationsRegister(hypreCitation,&cite);
2208:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
2209:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2210:   iupper[0] += ilower[0] - 1;
2211:   iupper[1] += ilower[1] - 1;
2212:   iupper[2] += ilower[2] - 1;
2213:   hlower[0]  = (HYPRE_Int)ilower[0];
2214:   hlower[1]  = (HYPRE_Int)ilower[1];
2215:   hlower[2]  = (HYPRE_Int)ilower[2];
2216:   hupper[0]  = (HYPRE_Int)iupper[0];
2217:   hupper[1]  = (HYPRE_Int)iupper[1];
2218:   hupper[2]  = (HYPRE_Int)iupper[2];

2220:   /* copy x values over to hypre */
2221:   PetscStackCallStandard(HYPRE_StructVectorSetConstantValues,(mx->hb,0.0));
2222:   VecGetArrayRead(x,&xx);
2223:   PetscStackCallStandard(HYPRE_StructVectorSetBoxValues,(mx->hb,hlower,hupper,(HYPRE_Complex*)xx));
2224:   VecRestoreArrayRead(x,&xx);
2225:   PetscStackCallStandard(HYPRE_StructVectorAssemble,(mx->hb));
2226:   PetscStackCallStandard(HYPRE_StructPFMGSolve,(ex->hsolver,mx->hmat,mx->hb,mx->hx));

2228:   /* copy solution values back to PETSc */
2229:   VecGetArray(y,&yy);
2230:   PetscStackCallStandard(HYPRE_StructVectorGetBoxValues,(mx->hx,hlower,hupper,(HYPRE_Complex*)yy));
2231:   VecRestoreArray(y,&yy);
2232:   return(0);
2233: }

2235: 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)
2236: {
2237:   PC_PFMG        *jac = (PC_PFMG*)pc->data;
2239:   HYPRE_Int      oits;

2242:   PetscCitationsRegister(hypreCitation,&cite);
2243:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,its*jac->its));
2244:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,rtol));

2246:   PCApply_PFMG(pc,b,y);
2247:   PetscStackCallStandard(HYPRE_StructPFMGGetNumIterations,(jac->hsolver,&oits));
2248:   *outits = oits;
2249:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2250:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
2251:   PetscStackCallStandard(HYPRE_StructPFMGSetTol,(jac->hsolver,jac->tol));
2252:   PetscStackCallStandard(HYPRE_StructPFMGSetMaxIter,(jac->hsolver,jac->its));
2253:   return(0);
2254: }


2257: PetscErrorCode PCSetUp_PFMG(PC pc)
2258: {
2259:   PetscErrorCode  ierr;
2260:   PC_PFMG         *ex = (PC_PFMG*) pc->data;
2261:   Mat_HYPREStruct *mx = (Mat_HYPREStruct*)(pc->pmat->data);
2262:   PetscBool       flg;

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

2268:   /* create the hypre solver object and set its information */
2269:   if (ex->hsolver) PetscStackCallStandard(HYPRE_StructPFMGDestroy,(ex->hsolver));
2270:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
2271:   PetscStackCallStandard(HYPRE_StructPFMGSetup,(ex->hsolver,mx->hmat,mx->hb,mx->hx));
2272:   PetscStackCallStandard(HYPRE_StructPFMGSetZeroGuess,(ex->hsolver));
2273:   return(0);
2274: }

2276: /*MC
2277:      PCPFMG - the hypre PFMG multigrid solver

2279:    Level: advanced

2281:    Options Database:
2282: + -pc_pfmg_its <its> number of iterations of PFMG to use as preconditioner
2283: . -pc_pfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2284: . -pc_pfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2285: . -pc_pfmg_tol <tol> tolerance of PFMG
2286: . -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
2287: - -pc_pfmg_rap_type - type of coarse matrix generation, one of Galerkin,non-Galerkin

2289:    Notes:
2290:     This is for CELL-centered descretizations

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

2295: .seealso:  PCMG, MATHYPRESTRUCT
2296: M*/

2298: PETSC_EXTERN PetscErrorCode PCCreate_PFMG(PC pc)
2299: {
2301:   PC_PFMG        *ex;

2304:   PetscNew(&ex); \
2305:   pc->data = ex;

2307:   ex->its            = 1;
2308:   ex->tol            = 1.e-8;
2309:   ex->relax_type     = 1;
2310:   ex->rap_type       = 0;
2311:   ex->num_pre_relax  = 1;
2312:   ex->num_post_relax = 1;
2313:   ex->max_levels     = 0;

2315:   pc->ops->setfromoptions  = PCSetFromOptions_PFMG;
2316:   pc->ops->view            = PCView_PFMG;
2317:   pc->ops->destroy         = PCDestroy_PFMG;
2318:   pc->ops->apply           = PCApply_PFMG;
2319:   pc->ops->applyrichardson = PCApplyRichardson_PFMG;
2320:   pc->ops->setup           = PCSetUp_PFMG;

2322:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
2323:   PetscStackCallStandard(HYPRE_StructPFMGCreate,(ex->hcomm,&ex->hsolver));
2324:   return(0);
2325: }

2327: /* ---------------------------------------------------------------------------------------------------------------------------------------------------*/

2329: /* we know we are working with a HYPRE_SStructMatrix */
2330: typedef struct {
2331:   MPI_Comm            hcomm;       /* does not share comm with HYPRE_SStructMatrix because need to create solver before getting matrix */
2332:   HYPRE_SStructSolver ss_solver;

2334:   /* keep copy of SYSPFMG options used so may view them */
2335:   PetscInt its;
2336:   double   tol;
2337:   PetscInt relax_type;
2338:   PetscInt num_pre_relax,num_post_relax;
2339: } PC_SysPFMG;

2341: PetscErrorCode PCDestroy_SysPFMG(PC pc)
2342: {
2344:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

2347:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2348:   MPI_Comm_free(&ex->hcomm);
2349:   PetscFree(pc->data);
2350:   return(0);
2351: }

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

2355: PetscErrorCode PCView_SysPFMG(PC pc,PetscViewer viewer)
2356: {
2358:   PetscBool      iascii;
2359:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;

2362:   PetscObjectTypeCompare((PetscObject)viewer,PETSCVIEWERASCII,&iascii);
2363:   if (iascii) {
2364:     PetscViewerASCIIPrintf(viewer,"  HYPRE SysPFMG preconditioning\n");
2365:     PetscViewerASCIIPrintf(viewer,"  max iterations %d\n",ex->its);
2366:     PetscViewerASCIIPrintf(viewer,"  tolerance %g\n",ex->tol);
2367:     PetscViewerASCIIPrintf(viewer,"  relax type %s\n",PFMGRelaxType[ex->relax_type]);
2368:     PetscViewerASCIIPrintf(viewer,"  number pre-relax %d post-relax %d\n",ex->num_pre_relax,ex->num_post_relax);
2369:   }
2370:   return(0);
2371: }

2373: PetscErrorCode PCSetFromOptions_SysPFMG(PetscOptionItems *PetscOptionsObject,PC pc)
2374: {
2376:   PC_SysPFMG     *ex = (PC_SysPFMG*) pc->data;
2377:   PetscBool      flg = PETSC_FALSE;

2380:   PetscOptionsHead(PetscOptionsObject,"SysPFMG options");
2381:   PetscOptionsBool("-pc_syspfmg_print_statistics","Print statistics","HYPRE_SStructSysPFMGSetPrintLevel",flg,&flg,NULL);
2382:   if (flg) {
2383:     int level=3;
2384:     PetscStackCallStandard(HYPRE_SStructSysPFMGSetPrintLevel,(ex->ss_solver,level));
2385:   }
2386:   PetscOptionsInt("-pc_syspfmg_its","Number of iterations of SysPFMG to use as preconditioner","HYPRE_SStructSysPFMGSetMaxIter",ex->its,&ex->its,NULL);
2387:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(ex->ss_solver,ex->its));
2388:   PetscOptionsInt("-pc_syspfmg_num_pre_relax","Number of smoothing steps before coarse grid","HYPRE_SStructSysPFMGSetNumPreRelax",ex->num_pre_relax,&ex->num_pre_relax,NULL);
2389:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPreRelax,(ex->ss_solver,ex->num_pre_relax));
2390:   PetscOptionsInt("-pc_syspfmg_num_post_relax","Number of smoothing steps after coarse grid","HYPRE_SStructSysPFMGSetNumPostRelax",ex->num_post_relax,&ex->num_post_relax,NULL);
2391:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetNumPostRelax,(ex->ss_solver,ex->num_post_relax));

2393:   PetscOptionsReal("-pc_syspfmg_tol","Tolerance of SysPFMG","HYPRE_SStructSysPFMGSetTol",ex->tol,&ex->tol,NULL);
2394:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(ex->ss_solver,ex->tol));
2395:   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);
2396:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetRelaxType,(ex->ss_solver, ex->relax_type));
2397:   PetscOptionsTail();
2398:   return(0);
2399: }

2401: PetscErrorCode PCApply_SysPFMG(PC pc,Vec x,Vec y)
2402: {
2403:   PetscErrorCode    ierr;
2404:   PC_SysPFMG        *ex = (PC_SysPFMG*) pc->data;
2405:   PetscScalar       *yy;
2406:   const PetscScalar *xx;
2407:   PetscInt          ilower[3],iupper[3];
2408:   HYPRE_Int         hlower[3],hupper[3];
2409:   Mat_HYPRESStruct  *mx     = (Mat_HYPRESStruct*)(pc->pmat->data);
2410:   PetscInt          ordering= mx->dofs_order;
2411:   PetscInt          nvars   = mx->nvars;
2412:   PetscInt          part    = 0;
2413:   PetscInt          size;
2414:   PetscInt          i;

2417:   PetscCitationsRegister(hypreCitation,&cite);
2418:   DMDAGetCorners(mx->da,&ilower[0],&ilower[1],&ilower[2],&iupper[0],&iupper[1],&iupper[2]);
2419:   /* when HYPRE_MIXEDINT is defined, sizeof(HYPRE_Int) == 32 */
2420:   iupper[0] += ilower[0] - 1;
2421:   iupper[1] += ilower[1] - 1;
2422:   iupper[2] += ilower[2] - 1;
2423:   hlower[0]  = (HYPRE_Int)ilower[0];
2424:   hlower[1]  = (HYPRE_Int)ilower[1];
2425:   hlower[2]  = (HYPRE_Int)ilower[2];
2426:   hupper[0]  = (HYPRE_Int)iupper[0];
2427:   hupper[1]  = (HYPRE_Int)iupper[1];
2428:   hupper[2]  = (HYPRE_Int)iupper[2];

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

2433:   /* copy x values over to hypre for variable ordering */
2434:   if (ordering) {
2435:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2436:     VecGetArrayRead(x,&xx);
2437:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(xx+(size*i))));
2438:     VecRestoreArrayRead(x,&xx);
2439:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2440:     PetscStackCallStandard(HYPRE_SStructMatrixMatvec,(1.0,mx->ss_mat,mx->ss_b,0.0,mx->ss_x));
2441:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

2443:     /* copy solution values back to PETSc */
2444:     VecGetArray(y,&yy);
2445:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(yy+(size*i))));
2446:     VecRestoreArray(y,&yy);
2447:   } else {      /* nodal ordering must be mapped to variable ordering for sys_pfmg */
2448:     PetscScalar *z;
2449:     PetscInt    j, k;

2451:     PetscMalloc1(nvars*size,&z);
2452:     PetscStackCallStandard(HYPRE_SStructVectorSetConstantValues,(mx->ss_b,0.0));
2453:     VecGetArrayRead(x,&xx);

2455:     /* transform nodal to hypre's variable ordering for sys_pfmg */
2456:     for (i= 0; i< size; i++) {
2457:       k= i*nvars;
2458:       for (j= 0; j< nvars; j++) z[j*size+i]= xx[k+j];
2459:     }
2460:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorSetBoxValues,(mx->ss_b,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i))));
2461:     VecRestoreArrayRead(x,&xx);
2462:     PetscStackCallStandard(HYPRE_SStructVectorAssemble,(mx->ss_b));
2463:     PetscStackCallStandard(HYPRE_SStructSysPFMGSolve,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));

2465:     /* copy solution values back to PETSc */
2466:     VecGetArray(y,&yy);
2467:     for (i= 0; i< nvars; i++) PetscStackCallStandard(HYPRE_SStructVectorGetBoxValues,(mx->ss_x,part,hlower,hupper,i,(HYPRE_Complex*)(z+(size*i))));
2468:     /* transform hypre's variable ordering for sys_pfmg to nodal ordering */
2469:     for (i= 0; i< size; i++) {
2470:       k= i*nvars;
2471:       for (j= 0; j< nvars; j++) yy[k+j]= z[j*size+i];
2472:     }
2473:     VecRestoreArray(y,&yy);
2474:     PetscFree(z);
2475:   }
2476:   return(0);
2477: }

2479: 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)
2480: {
2481:   PC_SysPFMG     *jac = (PC_SysPFMG*)pc->data;
2483:   HYPRE_Int      oits;

2486:   PetscCitationsRegister(hypreCitation,&cite);
2487:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,its*jac->its));
2488:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,rtol));
2489:   PCApply_SysPFMG(pc,b,y);
2490:   PetscStackCallStandard(HYPRE_SStructSysPFMGGetNumIterations,(jac->ss_solver,&oits));
2491:   *outits = oits;
2492:   if (oits == its) *reason = PCRICHARDSON_CONVERGED_ITS;
2493:   else             *reason = PCRICHARDSON_CONVERGED_RTOL;
2494:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetTol,(jac->ss_solver,jac->tol));
2495:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetMaxIter,(jac->ss_solver,jac->its));
2496:   return(0);
2497: }

2499: PetscErrorCode PCSetUp_SysPFMG(PC pc)
2500: {
2501:   PetscErrorCode   ierr;
2502:   PC_SysPFMG       *ex = (PC_SysPFMG*) pc->data;
2503:   Mat_HYPRESStruct *mx = (Mat_HYPRESStruct*)(pc->pmat->data);
2504:   PetscBool        flg;

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

2510:   /* create the hypre sstruct solver object and set its information */
2511:   if (ex->ss_solver) PetscStackCallStandard(HYPRE_SStructSysPFMGDestroy,(ex->ss_solver));
2512:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2513:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetZeroGuess,(ex->ss_solver));
2514:   PetscStackCallStandard(HYPRE_SStructSysPFMGSetup,(ex->ss_solver,mx->ss_mat,mx->ss_b,mx->ss_x));
2515:   return(0);
2516: }

2518: /*MC
2519:      PCSysPFMG - the hypre SysPFMG multigrid solver

2521:    Level: advanced

2523:    Options Database:
2524: + -pc_syspfmg_its <its> number of iterations of SysPFMG to use as preconditioner
2525: . -pc_syspfmg_num_pre_relax <steps> number of smoothing steps before coarse grid
2526: . -pc_syspfmg_num_post_relax <steps> number of smoothing steps after coarse grid
2527: . -pc_syspfmg_tol <tol> tolerance of SysPFMG
2528: . -pc_syspfmg_relax_type -relaxation type for the up and down cycles, one of Weighted-Jacobi,Red/Black-Gauss-Seidel

2530:    Notes:
2531:     This is for CELL-centered descretizations

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

2537: .seealso:  PCMG, MATHYPRESSTRUCT
2538: M*/

2540: PETSC_EXTERN PetscErrorCode PCCreate_SysPFMG(PC pc)
2541: {
2543:   PC_SysPFMG     *ex;

2546:   PetscNew(&ex); \
2547:   pc->data = ex;

2549:   ex->its            = 1;
2550:   ex->tol            = 1.e-8;
2551:   ex->relax_type     = 1;
2552:   ex->num_pre_relax  = 1;
2553:   ex->num_post_relax = 1;

2555:   pc->ops->setfromoptions  = PCSetFromOptions_SysPFMG;
2556:   pc->ops->view            = PCView_SysPFMG;
2557:   pc->ops->destroy         = PCDestroy_SysPFMG;
2558:   pc->ops->apply           = PCApply_SysPFMG;
2559:   pc->ops->applyrichardson = PCApplyRichardson_SysPFMG;
2560:   pc->ops->setup           = PCSetUp_SysPFMG;

2562:   MPI_Comm_dup(PetscObjectComm((PetscObject)pc),&(ex->hcomm));
2563:   PetscStackCallStandard(HYPRE_SStructSysPFMGCreate,(ex->hcomm,&ex->ss_solver));
2564:   return(0);
2565: }