Actual source code: hyppilut.c

  1: /*$Id: bvec2.c,v 1.202 2001/09/12 03:26:24 bsmith Exp $*/
  2: /*

  4: */

 6:  #include src/sles/pc/pcimpl.h
  7: EXTERN_C_BEGIN
  8: #include "HYPRE.h"
  9: #include "IJ_mv.h"
 10: #include "HYPRE_parcsr_ls.h"
 11: EXTERN_C_END

 13: extern int MatHYPRE_IJMatrixCreate(Mat,HYPRE_IJMatrix*);
 14: extern int MatHYPRE_IJMatrixCopy(Mat,HYPRE_IJMatrix);
 15: extern int VecHYPRE_IJVectorCreate(Vec,HYPRE_IJVector*);

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

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

 29:   /* options for pilut and BoomerAMG*/
 30:   int                maxiter;
 31:   double             tol;

 33:   /* options for pilut */
 34:   int                factorrowsize;

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

 46:   /* options for euclid */
 47:   PetscTruth         bjilu;
 48:   int                levels;

 50:   /* options for euclid and BoomerAMG */
 51:   PetscTruth         printstatistics;

 53:   /* options for BoomerAMG */
 54:   int                maxlevels;
 55:   double             strongthreshold;
 56:   double             maxrowsum;
 57:   int                *gridsweeps;
 58:   int                coarsentype;
 59:   int                measuretype;
 60:   int                *relaxtype;
 61:   int                **gridrelaxpoints;
 62: } PC_HYPRE;


 65: #undef __FUNCT__  
 67: static int PCSetUp_HYPRE(PC pc)
 68: {
 69:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
 70:   int                ierr;
 71:   HYPRE_ParCSRMatrix hmat;
 72:   HYPRE_ParVector    bv,xv;

 75:   if (!jac->ij) { /* create the matrix the first time through */
 76:     MatHYPRE_IJMatrixCreate(pc->pmat,&jac->ij);
 77:   }
 78:   if (!jac->b) {
 79:     VecHYPRE_IJVectorCreate(pc->vec,&jac->b);
 80:     VecHYPRE_IJVectorCreate(pc->vec,&jac->x);
 81:   }
 82:   MatHYPRE_IJMatrixCopy(pc->pmat,jac->ij);
 83:   HYPRE_IJMatrixGetObject(jac->ij,(void**)&hmat);
 84:   HYPRE_IJVectorGetObject(jac->b,(void**)&bv);
 85:   HYPRE_IJVectorGetObject(jac->x,(void**)&xv);
 86:   (*jac->setup)(jac->hsolver,hmat,bv,xv);
 87:   if (ierr) SETERRQ1(1,"Error in HYPRE setup, error code %d",ierr);
 88:   return(0);
 89: }

 91: /*
 92:     Replaces the address where the HYPRE vector points to its data with the address of
 93:   PETSc's data. Saves the old address so it can be reset when we are finished with it.
 94:   Allows use to get the data into a HYPRE vector without the cost of memcopies 
 95: */
 96: #define HYPREReplacePointer(b,newvalue,savedvalue) {
 97:    hypre_ParVector *par_vector   = (hypre_ParVector *)hypre_IJVectorObject(((hypre_IJVector*)b));
 98:    hypre_Vector    *local_vector = hypre_ParVectorLocalVector(par_vector);
 99:    savedvalue         = local_vector->data;
100:    local_vector->data = newvalue;}

102: #undef __FUNCT__  
104: static int PCApply_HYPRE(PC pc,Vec b,Vec x)
105: {
106:   PC_HYPRE           *jac = (PC_HYPRE*)pc->data;
107:   int                ierr;
108:   HYPRE_ParCSRMatrix hmat;
109:   PetscScalar        *bv,*xv;
110:   HYPRE_ParVector    jbv,jxv;
111:   PetscScalar        *sbv,*sxv;

114:   VecGetArray(b,&bv);
115:   VecGetArray(x,&xv);
116:   HYPREReplacePointer(jac->b,bv,sbv);
117:   HYPREReplacePointer(jac->x,xv,sxv);

119:   HYPRE_IJMatrixGetObject(jac->ij,(void**)&hmat);
120:   HYPRE_IJVectorGetObject(jac->b,(void**)&jbv);
121:   HYPRE_IJVectorGetObject(jac->x,(void**)&jxv);
122:   (*jac->solve)(jac->hsolver,hmat,jbv,jxv);
123:   /* error code of 1 in boomerAMG merely means convergence not achieved */
124:   if (ierr && (ierr != 1 || jac->solve != HYPRE_BoomerAMGSolve)) SETERRQ1(1,"Error in HYPRE solver, error code %d",ierr);
125: 

127:   HYPREReplacePointer(jac->b,sbv,bv);
128:   HYPREReplacePointer(jac->x,sxv,xv);
129:   VecRestoreArray(x,&xv);
130:   VecRestoreArray(b,&bv);

132:   return(0);
133: }

135: #undef __FUNCT__  
137: static int PCDestroy_HYPRE(PC pc)
138: {
139:   PC_HYPRE *jac = (PC_HYPRE*)pc->data;
140:   int      ierr;

143:   HYPRE_IJMatrixDestroy(jac->ij);
144:   HYPRE_IJVectorDestroy(jac->b);
145:   HYPRE_IJVectorDestroy(jac->x);
146:   (*jac->destroy)(jac->hsolver);
147:   PetscFree(jac);
148:   return(0);
149: }

151: /* --------------------------------------------------------------------------------------------*/
152: #undef __FUNCT__  
154: static int PCSetFromOptions_HYPRE_Pilut(PC pc)
155: {
156:   PC_HYPRE  *jac = (PC_HYPRE*)pc->data;
157:   int        ierr;
158:   PetscTruth flag;

161:   PetscOptionsHead("HYPRE Pilut Options");
162:     PetscOptionsInt("-pc_hypre_pilut_maxiter","Number of iterations","None",jac->maxiter,&jac->maxiter,&flag);
163:     if (flag) {
164:       HYPRE_ParCSRPilutSetMaxIter(jac->hsolver,jac->maxiter);
165:     }
166:     PetscOptionsReal("-pc_hypre_pilut_tol","Drop tolerance","None",jac->tol,&jac->tol,&flag);
167:     if (flag) {
168:       HYPRE_ParCSRPilutSetDropTolerance(jac->hsolver,jac->tol);
169:     }
170:     PetscOptionsInt("-pc_hypre_pilut_factorrowsize","FactorRowSize","None",jac->factorrowsize,&jac->factorrowsize,&flag);
171:     if (flag) {
172:       HYPRE_ParCSRPilutSetFactorRowSize(jac->hsolver,jac->factorrowsize);
173:     }
174:   PetscOptionsTail();
175:   return(0);
176: }

178: #undef __FUNCT__  
180: static int PCView_HYPRE_Pilut(PC pc,PetscViewer viewer)
181: {
182:   PC_HYPRE    *jac = (PC_HYPRE*)pc->data;
183:   int         ierr;
184:   PetscTruth  isascii;

187:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
188:   if (isascii) {
189:     PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut preconditioningn");
190:     if (jac->maxiter != PETSC_DEFAULT) {
191:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: maximum number of iterations %dn",jac->maxiter);
192:     } else {
193:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default maximum number of iterations n");
194:     }
195:     if (jac->tol != PETSC_DEFAULT) {
196:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: drop tolerance %gn",jac->tol);
197:     } else {
198:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default drop tolerance n");
199:     }
200:     if (jac->factorrowsize != PETSC_DEFAULT) {
201:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: factor row size %dn",jac->factorrowsize);
202:     } else {
203:       PetscViewerASCIIPrintf(viewer,"  HYPRE Pilut: default factor row size n");
204:     }
205:   }
206:   return(0);
207: }

209: /* --------------------------------------------------------------------------------------------*/
210: #undef __FUNCT__  
212: static int PCSetFromOptions_HYPRE_Euclid(PC pc)
213: {
214:   PC_HYPRE  *jac = (PC_HYPRE*)pc->data;
215:   int        ierr;
216:   PetscTruth flag;
217:   char       *args[2];

220:   jac->bjilu              = PETSC_FALSE;
221:   jac->levels             = 1;

223:   PetscOptionsHead("HYPRE Euclid Options");
224:     PetscOptionsInt("-pc_hypre_euclid_levels","Number of levels of fill ILU(k)","None",jac->levels,&jac->levels,&flag);
225:     if (flag) {
226:       char levels[16];
227:       if (jac->levels < 0) SETERRQ1(1,"Number of levels %d must be nonegative",jac->levels);
228:       sprintf(levels,"%d",jac->levels);
229:       args[0] = "-level"; args[1] = levels;
230:       HYPRE_EuclidSetParams(jac->hsolver,2,args);
231:     }
232:     PetscOptionsLogical("-pc_hypre_euclid_bj","Use block Jacobi ILU(k)","None",jac->bjilu,&jac->bjilu,PETSC_NULL);
233:     if (jac->bjilu) {
234:       args[0] = "-bj"; args[1] = "1";
235:       HYPRE_EuclidSetParams(jac->hsolver,2,args);
236:     }
237: 
238:     PetscOptionsLogical("-pc_hypre_euclid_print_statistics","Print statistics","None",jac->printstatistics,&jac->printstatistics,PETSC_NULL);
239:     if (jac->printstatistics) {
240:       args[0] = "-eu_stats"; args[1] = "1";
241:       HYPRE_EuclidSetParams(jac->hsolver,2,args);
242:       args[0] = "-eu_mem"; args[1] = "1";
243:       HYPRE_EuclidSetParams(jac->hsolver,2,args);
244:     }
245:   PetscOptionsTail();
246:   return(0);
247: }

249: #undef __FUNCT__  
251: static int PCView_HYPRE_Euclid(PC pc,PetscViewer viewer)
252: {
253:   PC_HYPRE    *jac = (PC_HYPRE*)pc->data;
254:   int         ierr;
255:   PetscTruth  isascii;

258:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
259:   if (isascii) {
260:     PetscViewerASCIIPrintf(viewer,"  HYPRE Euclid preconditioningn");
261:     PetscViewerASCIIPrintf(viewer,"  HYPRE Euclid: number of levels %dn",jac->levels);
262:     if (jac->bjilu) {
263:       PetscViewerASCIIPrintf(viewer,"  HYPRE Euclid: Using block Jacobi ILU instead of parallel ILUn");
264:     }
265:   }
266:   return(0);
267: }

269: /* --------------------------------------------------------------------------------------------*/

271: static char *HYPREBoomerAMGCoarsenType[] = {"CLJP","Ruge-Stueben","","modifiedRuge-Stueben","","","Falgout"};
272: static char *HYPREBoomerAMGMeasureType[] = {"local","global"};
273: static char *HYPREBoomerAMGRelaxType[]   = {"Jacobi","sequential-Gauss-Seidel","","Gauss-Seidel/Jacobi","","","symmetric-Gauss-Seidel/Jacobi",
274:                                             "","","Gaussian-elimination"};
275: #undef __FUNCT__  
277: static int PCSetFromOptions_HYPRE_BoomerAMG(PC pc)
278: {
279:   PC_HYPRE  *jac = (PC_HYPRE*)pc->data;
280:   int        ierr,n = 4,i;
281:   PetscTruth flg;
282:   char       result[32];

285:   jac->maxlevels       = 25;
286:   jac->maxiter         = 20;
287:   jac->tol             = 1.e-7;
288:   jac->strongthreshold = .25;
289:   jac->maxrowsum       = .9;
290:   jac->coarsentype     = 6;
291:   jac->measuretype     = 0;
292: 

294:   /* this is terrible; HYPRE frees this array so we have to malloc it */
295:   jac->gridsweeps    = (int*)malloc(4*sizeof(int));
296:   jac->gridsweeps[0] = jac->gridsweeps[1] = jac->gridsweeps[2] = 2;
297:   jac->gridsweeps[3] = 1;

299:   jac->relaxtype     = (int*)malloc(4*sizeof(int));
300:   jac->relaxtype[0]  = jac->relaxtype[1] = jac->relaxtype[2] = 3;
301:   jac->relaxtype[3]  = 9;

303:   PetscOptionsHead("HYPRE BoomerAMG Options");
304:     PetscOptionsInt("-pc_hypre_boomeramg_max_levels","Number of levels (of grids) allowed","None",jac->maxlevels,&jac->maxlevels,&flg);
305:     if (flg) {
306:       if (jac->maxlevels < 2) SETERRQ1(1,"Number of levels %d must be at least two",jac->maxlevels);
307:       HYPRE_BoomerAMGSetMaxLevels(jac->hsolver,jac->maxlevels);
308:     }
309:     PetscOptionsInt("-pc_hypre_boomeramg_max_iter","Maximum iterations used","None",jac->maxiter,&jac->maxiter,&flg);
310:     if (flg) {
311:       if (jac->maxiter < 1) SETERRQ1(1,"Number of iterations %d must be at least one",jac->maxiter);
312:       HYPRE_BoomerAMGSetMaxIter(jac->hsolver,jac->maxiter);
313:     }
314:     PetscOptionsScalar("-pc_hypre_boomeramg_tol","Convergence tolerance","None",jac->tol,&jac->tol,&flg);
315:     if (flg) {
316:       if (jac->tol < 0.0) SETERRQ1(1,"Tolerance %g must be great than or equal zero",jac->tol);
317:       HYPRE_BoomerAMGSetTol(jac->hsolver,jac->tol);
318:     }
319:     PetscOptionsScalar("-pc_hypre_boomeramg_strong_threshold","Threshold for being strongly connected","None",jac->strongthreshold,&jac->strongthreshold,&flg);
320:     if (flg) {
321:       if (jac->strongthreshold < 0.0) SETERRQ1(1,"Strong threshold %g must be great than or equal zero",jac->strongthreshold);
322:       HYPRE_BoomerAMGSetStrongThreshold(jac->hsolver,jac->strongthreshold);
323:     }
324:     PetscOptionsScalar("-pc_hypre_boomeramg_max_row_sum","Maximum row sum","None",jac->maxrowsum,&jac->maxrowsum,&flg);
325:     if (flg) {
326:       if (jac->maxrowsum < 0.0) SETERRQ1(1,"Maximum row sum %g must be greater than zero",jac->maxrowsum);
327:       if (jac->maxrowsum > 1.0) SETERRQ1(1,"Maximum row sum %g must be less than or equal one",jac->maxrowsum);
328:       HYPRE_BoomerAMGSetMaxRowSum(jac->hsolver,jac->maxrowsum);
329:     }
330: 
331:     n = 4;
332:     PetscOptionsIntArray("-pc_hypre_boomeramg_grid_sweeps","Grid sweeps for fine,down,up,coarse","None",jac->gridsweeps,&n,&flg);
333:     if (flg) {
334:       if (n == 1) {
335:         jac->gridsweeps[1] = jac->gridsweeps[2] =  jac->gridsweeps[3] = jac->gridsweeps[0];
336:         n = 4;
337:       }
338:       if (n != 4) SETERRQ1(1,"You must provide either 1 or 4 values seperated by commas, you provided %d",n);
339:       HYPRE_BoomerAMGSetNumGridSweeps(jac->hsolver,jac->gridsweeps);
340:       CHKMEMQ;
341:     }

343:     /*
344:          Suggested by QUANDALLE Philippe <Philippe.QUANDALLE@ifp.fr>

346:         gridrelaxpoints[i][j] are for i=0,1,2,3 (fine,down,up,coarse) and j=sweep number
347:         0 indicates smooth all points
348:         1 indicates smooth coarse points
349:        -1 indicates smooth fine points

351:          Here when j=1 it first smooths all the coarse points, then all the fine points.
352:     */
353:     jac->gridrelaxpoints    = (int**)malloc(4*sizeof(int*));
354:     if(jac->gridsweeps[0]>0) jac->gridrelaxpoints[0] = (int*)malloc(jac->gridsweeps[0]*sizeof(int));
355:     if(jac->gridsweeps[1]>0) jac->gridrelaxpoints[1] = (int*)malloc(jac->gridsweeps[1]*sizeof(int));
356:     if(jac->gridsweeps[2]>0) jac->gridrelaxpoints[2] = (int*)malloc(jac->gridsweeps[2]*sizeof(int));
357:     if(jac->gridsweeps[3]>0) jac->gridrelaxpoints[3] = (int*)malloc(jac->gridsweeps[3]*sizeof(int));

359:     PetscOptionsLogical("-pc_hypre_boomeramg_sweep_all","Sweep all points","None",PETSC_FALSE,&flg,0);
360:     if(jac->gridsweeps[0] == 1) jac->gridrelaxpoints[0][0] = 0;
361:     else if(jac->gridsweeps[0] == 2) {
362:       if (flg) {
363:         jac->gridrelaxpoints[0][0] = 0; jac->gridrelaxpoints[0][1] = 0;
364:       } else {
365:         jac->gridrelaxpoints[0][0] = 1; jac->gridrelaxpoints[0][1] = -1;
366:       }
367:     } else if (jac->gridsweeps[0] > 2) {
368:       SETERRQ(1,"Grid sweeps can only be 0, 1, or 2");
369:     }
370: 
371:     if(jac->gridsweeps[1] == 1) jac->gridrelaxpoints[1][0] = 0;
372:     else if(jac->gridsweeps[1] == 2) {
373:       if (flg) {
374:         jac->gridrelaxpoints[1][0] = 0;  jac->gridrelaxpoints[1][1] = 0;
375:       } else {
376:         jac->gridrelaxpoints[1][0] = 1; jac->gridrelaxpoints[1][1] = -1;
377:       }
378:     } else if (jac->gridsweeps[1] > 2) {
379:       SETERRQ(1,"Grid sweeps can only be 0, 1, or 2");
380:     }

382:     if(jac->gridsweeps[2] == 1) jac->gridrelaxpoints[2][0] = 0;
383:     else if(jac->gridsweeps[2] == 2) {
384:       if (flg) {
385:         jac->gridrelaxpoints[2][0] = 0; jac->gridrelaxpoints[2][1] = 0;
386:       } else {
387:         jac->gridrelaxpoints[2][0] = -1; jac->gridrelaxpoints[2][1] = 1;
388:       }
389:     } else if (jac->gridsweeps[2] > 2) {
390:       SETERRQ(1,"Grid sweeps can only be 0, 1, or 2");
391:     }

393:     for (i=0; i<jac->gridsweeps[3]; i++) {
394:       jac->gridrelaxpoints[3][i] = 0;
395:     }
396:     HYPRE_BoomerAMGSetGridRelaxPoints(jac->hsolver,jac->gridrelaxpoints);


399:     PetscOptionsEList("-pc_hypre_boomeramg_measure_type","Measure type","None",HYPREBoomerAMGMeasureType,2,HYPREBoomerAMGMeasureType[0],result,16,&flg);
400:     if (flg) {
401:       int type = -1;
402:       for (i=0; i<2; i++) {
403:         PetscStrcmp(result,HYPREBoomerAMGMeasureType[i],&flg);
404:         if (flg) {
405:           type = i;
406:           break;
407:         }
408:       }
409:       if (type == -1) SETERRQ1(1,"Unknown measure type %s",result);
410:       jac->measuretype = type;
411:       HYPRE_BoomerAMGSetMeasureType(jac->hsolver,jac->measuretype);
412:     }
413:     PetscOptionsEList("-pc_hypre_boomeramg_coarsen_type","Coarsen type","None",HYPREBoomerAMGCoarsenType,7,HYPREBoomerAMGCoarsenType[6],result,25,&flg);
414:     if (flg) {
415:       int type = -1;
416:       for (i=0; i<7; i++) {
417:         PetscStrcmp(result,HYPREBoomerAMGCoarsenType[i],&flg);
418:         if (flg) {
419:           type = i;
420:           break;
421:         }
422:       }
423:       if (type == -1) SETERRQ1(1,"Unknown coarsen type %s",result);
424:       jac->coarsentype = type;
425:       HYPRE_BoomerAMGSetCoarsenType(jac->hsolver,jac->coarsentype);
426:     }
427:     PetscOptionsEList("-pc_hypre_boomeramg_relax_type","Relax type","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[3],result,32,&flg);
428:     if (flg) {
429:       int type = -1;
430:       for (i=0; i<10; i++) {
431:         PetscStrcmp(result,HYPREBoomerAMGRelaxType[i],&flg);
432:         if (flg) {
433:           type = i;
434:           break;
435:         }
436:       }
437:       if (type == -1) SETERRQ1(1,"Unknown relax type %s",result);
438:       jac->relaxtype[0] = jac->relaxtype[1] = jac->relaxtype[2] = type;
439:     }
440:     PetscOptionsEList("-pc_hypre_boomeramg_relax_type_coarse","Relax type on coarse grid","None",HYPREBoomerAMGRelaxType,10,HYPREBoomerAMGRelaxType[3],result,32,&flg);
441:     if (flg) {
442:       int type = -1;
443:       for (i=0; i<10; i++) {
444:         PetscStrcmp(result,HYPREBoomerAMGRelaxType[i],&flg);
445:         if (flg) {
446:           type = i;
447:           break;
448:         }
449:       }
450:       if (type == -1) SETERRQ1(1,"Unknown relax type %s",result);
451:       jac->relaxtype[3] = type;
452:     }
453:     HYPRE_BoomerAMGSetGridRelaxType(jac->hsolver,jac->relaxtype);
454: 
455:     PetscOptionsLogical("-pc_hypre_boomeramg_print_statistics","Print statistics","None",jac->printstatistics,&jac->printstatistics,PETSC_NULL);
456:     if (jac->printstatistics) {
457:       HYPRE_BoomerAMGSetPrintLevel(jac->hsolver,3);
458:     }
459:   PetscOptionsTail();
460:   return(0);
461: }

463: #undef __FUNCT__  
465: static int PCView_HYPRE_BoomerAMG(PC pc,PetscViewer viewer)
466: {
467:   PC_HYPRE    *jac = (PC_HYPRE*)pc->data;
468:   int         ierr;
469:   PetscTruth  isascii;

472:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
473:   if (isascii) {
474:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG preconditioningn");
475:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of levels %dn",jac->maxlevels);
476:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum number of iterations %dn",jac->maxiter);
477:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Convergence tolerance %gn",jac->tol);
478:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Threshold for strong coupling %gn",jac->strongthreshold);
479:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Maximum row sums %gn",jac->maxrowsum);
480:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps on fine grid %dn",jac->gridsweeps[0]);
481:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps down         %dn",jac->gridsweeps[1]);
482:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps up           %dn",jac->gridsweeps[2]);
483:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Sweeps on coarse    %dn",jac->gridsweeps[3]);

485:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax on fine grid %sn",HYPREBoomerAMGRelaxType[jac->relaxtype[0]]);
486:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax down         %sn",HYPREBoomerAMGRelaxType[jac->relaxtype[1]]);
487:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax up           %sn",HYPREBoomerAMGRelaxType[jac->relaxtype[2]]);
488:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Relax on coarse    %sn",HYPREBoomerAMGRelaxType[jac->relaxtype[3]]);

490:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Measure type    %sn",HYPREBoomerAMGMeasureType[jac->measuretype]);
491:     PetscViewerASCIIPrintf(viewer,"  HYPRE BoomerAMG: Coarsen type    %sn",HYPREBoomerAMGCoarsenType[jac->coarsentype]);
492:   }
493:   return(0);
494: }

496: /* --------------------------------------------------------------------------------------------*/
497: #undef __FUNCT__  
499: static int PCSetFromOptions_HYPRE_ParaSails(PC pc)
500: {
501:   PC_HYPRE  *jac = (PC_HYPRE*)pc->data;
502:   int        ierr;
503:   PetscTruth flag;
504:   char       *symtlist[] = {"nonsymmetric","SPD","nonsymmetric,SPD"},buff[32];

507:   jac->nlevels     = 1;
508:   jac->threshhold  = .1;
509:   jac->filter      = .1;
510:   jac->loadbal     = 0;
511:   if (PetscLogPrintInfo) {
512:     jac->logging     = (int) PETSC_TRUE;
513:   } else {
514:     jac->logging     = (int) PETSC_FALSE;
515:   }
516:   jac->ruse = (int) PETSC_FALSE;
517:   jac->symt   = 0;

519:   PetscOptionsHead("HYPRE ParaSails Options");
520:     PetscOptionsInt("-pc_hypre_parasails_nlevels","Number of number of levels","None",jac->nlevels,&jac->nlevels,0);
521:     PetscOptionsReal("-pc_hypre_parasails_thresh","Threshold","None",jac->threshhold,&jac->threshhold,0);
522:     HYPRE_ParaSailsSetParams(jac->hsolver,jac->threshhold,jac->nlevels);

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

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

530:     PetscOptionsLogical("-pc_hypre_parasails_logging","Print info to screen","None",(PetscTruth)jac->logging,(PetscTruth*)&jac->logging,0);
531:     HYPRE_ParaSailsSetLogging(jac->hsolver,jac->logging);

533:     PetscOptionsLogical("-pc_hypre_parasails_reuse","Reuse nonzero pattern in preconditioner","None",(PetscTruth)jac->ruse,(PetscTruth*)&jac->ruse,0);
534:     HYPRE_ParaSailsSetReuse(jac->hsolver,jac->ruse);

536:     PetscOptionsEList("-pc_hypre_parasails_sym","Symmetry of matrix and preconditioner","None",symtlist,3,symtlist[0],buff,32,&flag);
537:     while (flag) {
538:       PetscStrcmp(symtlist[0],buff,&flag);
539:       if (flag) {
540:         jac->symt = 0;
541:         break;
542:       }
543:       PetscStrcmp(symtlist[1],buff,&flag);
544:       if (flag) {
545:         jac->symt = 1;
546:         break;
547:       }
548:       PetscStrcmp(symtlist[2],buff,&flag);
549:       if (flag) {
550:         jac->symt = 2;
551:         break;
552:       }
553:       SETERRQ1(1,"Unknown HYPRE ParaSails Sym option %s",buff);
554:     }
555:     HYPRE_ParaSailsSetSym(jac->hsolver,jac->symt);

557:   PetscOptionsTail();
558:   return(0);
559: }

561: #undef __FUNCT__  
563: static int PCView_HYPRE_ParaSails(PC pc,PetscViewer viewer)
564: {
565:   PC_HYPRE    *jac = (PC_HYPRE*)pc->data;
566:   int         ierr;
567:   PetscTruth  isascii;
568:   char        *symt;

571:   PetscTypeCompare((PetscObject)viewer,PETSC_VIEWER_ASCII,&isascii);
572:   if (isascii) {
573:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails preconditioningn");
574:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: nlevels %dn",jac->nlevels);
575:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: threshold %gn",jac->threshhold);
576:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: filter %gn",jac->filter);
577:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: load balance %gn",jac->loadbal);
578:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: reuse nonzero structure %sn",jac->ruse ? "true" : "false");
579:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: print info to screen %sn",jac->logging ? "true" : "false");
580:     if (jac->symt == 0) {
581:       symt = "nonsymmetric matrix and preconditioner";
582:     } else if (jac->symt == 1) {
583:       symt = "SPD matrix and preconditioner";
584:     } else if (jac->symt == 2) {
585:       symt = "nonsymmetric matrix but SPD preconditioner";
586:     } else {
587:       SETERRQ1(1,"Unknown HYPRE ParaSails sym option %d",jac->symt);
588:     }
589:     PetscViewerASCIIPrintf(viewer,"  HYPRE ParaSails: %sn",symt);
590:   }
591:   return(0);
592: }
593: /* ---------------------------------------------------------------------------------*/

595: #undef __FUNCT__  
597: static int PCHYPRESetType_HYPRE(PC pc,char *name)
598: {
599:   PC_HYPRE   *jac = (PC_HYPRE*)pc->data;
600:   int        ierr;
601:   PetscTruth flag;

604:   if (pc->ops->setup) {
605:     SETERRQ(1,"Cannot set the HYPRE preconditioner type once it has been set");
606:   }

608:   pc->ops->setup          = PCSetUp_HYPRE;
609:   pc->ops->apply          = PCApply_HYPRE;
610:   pc->ops->destroy        = PCDestroy_HYPRE;

612:   jac->maxiter            = PETSC_DEFAULT;
613:   jac->tol                = PETSC_DEFAULT;
614:   jac->printstatistics    = PetscLogPrintInfo;

616:   PetscStrcmp("pilut",name,&flag);
617:   if (flag) {
618:     ierr                    = HYPRE_ParCSRPilutCreate(pc->comm,&jac->hsolver);
619:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Pilut;
620:     pc->ops->view           = PCView_HYPRE_Pilut;
621:     jac->destroy            = HYPRE_ParCSRPilutDestroy;
622:     jac->setup              = HYPRE_ParCSRPilutSetup;
623:     jac->solve              = HYPRE_ParCSRPilutSolve;
624:     jac->factorrowsize      = PETSC_DEFAULT;
625:     return(0);
626:   }
627:   PetscStrcmp("parasails",name,&flag);
628:   if (flag) {
629:     ierr                    = HYPRE_ParaSailsCreate(pc->comm,&jac->hsolver);
630:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_ParaSails;
631:     pc->ops->view           = PCView_HYPRE_ParaSails;
632:     jac->destroy            = HYPRE_ParaSailsDestroy;
633:     jac->setup              = HYPRE_ParaSailsSetup;
634:     jac->solve              = HYPRE_ParaSailsSolve;
635:     return(0);
636:   }
637:   PetscStrcmp("euclid",name,&flag);
638:   if (flag) {
639:     ierr                    = HYPRE_EuclidCreate(pc->comm,&jac->hsolver);
640:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_Euclid;
641:     pc->ops->view           = PCView_HYPRE_Euclid;
642:     jac->destroy            = HYPRE_EuclidDestroy;
643:     jac->setup              = HYPRE_EuclidSetup;
644:     jac->solve              = HYPRE_EuclidSolve;
645:     return(0);
646:   }
647:   PetscStrcmp("boomeramg",name,&flag);
648:   if (flag) {
649:     ierr                    = HYPRE_BoomerAMGCreate(&jac->hsolver);
650:     pc->ops->setfromoptions = PCSetFromOptions_HYPRE_BoomerAMG;
651:     pc->ops->view           = PCView_HYPRE_BoomerAMG;
652:     jac->destroy            = HYPRE_BoomerAMGDestroy;
653:     jac->setup              = HYPRE_BoomerAMGSetup;
654:     jac->solve              = HYPRE_BoomerAMGSolve;
655:     return(0);
656:   }
657:   SETERRQ1(1,"Unknown HYPRE preconditioner %s; Choices are pilut, parasails, euclid, boomeramg",name);
658:   return(0);
659: }

661: /*
662:     It only gets here if the HYPRE type has not been set before the call to 
663:    ...SetFromOptions() which actually is most of the time
664: */
665: #undef __FUNCT__  
667: static int PCSetFromOptions_HYPRE(PC pc)
668: {
669:   int        ierr;
670:   char       buff[32],*type[] = {"pilut","parasails","boomerAMG","euclid"};
671:   PetscTruth flg;

674:   PetscOptionsHead("HYPRE preconditioner options");
675:     PetscOptionsEList("-pc_hypre_type","HYPRE preconditioner type","PCHYPRESetType",type,4,"pilut",buff,32,&flg);

677: 
678:     if (PetscOptionsPublishCount) {   /* force the default if it was not yet set and user did not set with option */
679:       if (!flg && !pc->ops->apply) {
680:         flg  = PETSC_TRUE;
681:         PetscStrcpy(buff,"pilut");
682:       }
683:     }

685:     if (flg) {
686:       PCHYPRESetType_HYPRE(pc,buff);
687:     }
688:     if (pc->ops->setfromoptions) {
689:       pc->ops->setfromoptions(pc);
690:     }
691:   PetscOptionsTail();
692:   return(0);
693: }

695: #undef __FUNCT__  
697: /*@
698:      PCHYPRESetType - Sets which hypre preconditioner you wish to use

700:    Input Parameters:
701: +     pc - the preconditioner context
702: -     name - either  pilut, parasails, boomerAMG, euclid

704:    Options Database Keys:
705:    -pc_hypre_type - One of pilut, parasails, boomerAMG, euclid
706:  
707:    Level: intermediate

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

712: @*/
713: int PCHYPRESetType(PC pc,char *name)
714: {
715:   int ierr,(*f)(PC,char*);

719:   PetscObjectQueryFunction((PetscObject)pc,"PCHYPRESetType_C",(void (**)(void))&f);
720:   if (f) {
721:     (*f)(pc,name);
722:   }
723:   return(0);
724: }

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

729:    Options Database Keys:
730: +   -pc_hypre_type - One of pilut, parasails, boomerAMG, euclid
731: -   Too many others to list, run with -pc_type hypre -pc_hypre_type XXX to see options for the XXX
732:           preconditioner
733:  
734:    Level: intermediate

736:    Notes: The many hypre options can ONLY be set via the options database (e.g. the command line
737:           or with PetscOptionsSetValue(), there are no functions to set them)

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

742: M*/

744: EXTERN_C_BEGIN
745: #undef __FUNCT__  
747: int PCCreate_HYPRE(PC pc)
748: {
749:   PC_HYPRE *jac;
750:   int       ierr;

753:   ierr                     = PetscNew(PC_HYPRE,&jac);
754:   ierr                     = PetscMemzero(jac,sizeof(PC_HYPRE));
755:   pc->data                 = jac;
756:   pc->ops->setfromoptions  = PCSetFromOptions_HYPRE;
757:   PetscObjectComposeFunctionDynamic((PetscObject)pc,"PCHYPRESetType_C",
758:                                     "PCHYPRESetType_HYPRE",PCHYPRESetType_HYPRE);
759:   return(0);
760: }
761: EXTERN_C_END