hypre.hpp

Go to the documentation of this file.

// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
// reserved. See file COPYRIGHT for details.
//
// This file is part of the MFEM library. For more information and source code
// availability see http://mfem.googlecode.com.
//
// MFEM is free software; you can redistribute it and/or modify it under the
// terms of the GNU Lesser General Public License (as published by the Free
// Software Foundation) version 2.1 dated February 1999.

#ifndef MFEM_HYPRE
#define MFEM_HYPRE

// Enable internal hypre timing routines
#define HYPRE_TIMING

// hypre header files
#include "seq_mv.h"
#include "_hypre_parcsr_mv.h"
#include "_hypre_parcsr_ls.h"

class HypreParVector : public Vector
{
private:
   int own_ParVector;

   hypre_ParVector *x;

   friend class HypreParMatrix;

public:
   HypreParVector(int glob_size, int *col);
   HypreParVector(int glob_size, double *_data, int *col);
   HypreParVector(const HypreParVector &y);
   HypreParVector(HYPRE_ParVector y);

   operator hypre_ParVector*() const;
   operator HYPRE_ParVector() const;
   hypre_ParVector *StealParVector() { own_ParVector = 0; return x; }

   Vector *GlobalVector();

   HypreParVector& operator= (double d);
   HypreParVector& operator= (const HypreParVector &y);

   void SetData(double *_data);

   int Randomize(int seed);

   void Print(const char *fname);

   ~HypreParVector();
};

double InnerProduct(HypreParVector &x, HypreParVector &y);
double InnerProduct(HypreParVector *x, HypreParVector *y);


class HypreParMatrix : public Operator
{
private:
   hypre_ParCSRMatrix *A;

   hypre_ParCSRCommPkg *CommPkg;

   mutable HypreParVector *X, *Y;

public:
   HypreParMatrix(hypre_ParCSRMatrix *a) : A(a)
   { X = Y = 0; CommPkg = 0; }
   HypreParMatrix(int size, int *row, SparseMatrix *diag);
   HypreParMatrix(int M, int N, int *row, int *col, SparseMatrix *diag);
   HypreParMatrix(int M, int N, int *row, int *col, SparseMatrix *diag,
                  SparseMatrix *offd, int *cmap);

   HypreParMatrix(int *row, int *col, SparseMatrix *a);

   HypreParMatrix(int M, int N, int *row, int *col, Table *diag);
   HypreParMatrix(MPI_Comm comm, int id, int np, int *row, int *col,
                  int *i_diag, int *j_diag, int *i_offd, int *j_offd,
                  int *cmap, int cmap_size);

   HypreParMatrix(MPI_Comm comm, int nrows, int glob_nrows, int glob_ncols,
                  int *I, int *J, double *data, int *rows, int *cols);

   // hypre's communication package object
   void SetCommPkg(hypre_ParCSRCommPkg *comm_pkg);
   void CheckCommPkg();
   void DestroyCommPkg();

   operator hypre_ParCSRMatrix*();
   operator HYPRE_ParCSRMatrix();
   hypre_ParCSRMatrix* StealData();

   inline int NNZ() { return A->num_nonzeros; }
   inline int * RowPart() { return A->row_starts; }
   inline int * ColPart() { return A->col_starts; }
   inline int M() { return A -> global_num_rows; }
   inline int N() { return A -> global_num_cols; }

   void GetDiag(Vector &diag);
   HypreParMatrix * Transpose();

   int GetNumRows() const
   { return hypre_CSRMatrixNumRows(hypre_ParCSRMatrixDiag(A)); }

   int GetGlobalNumRows() const { return hypre_ParCSRMatrixGlobalNumRows(A); }

   int GetGlobalNumCols() const { return hypre_ParCSRMatrixGlobalNumCols(A); }

   int *GetRowStarts() const { return hypre_ParCSRMatrixRowStarts(A); }

   int *GetColStarts() const { return hypre_ParCSRMatrixColStarts(A); }

   int Mult(HypreParVector &x, HypreParVector &y,
            double alpha = 1.0, double beta = 0.0);
   int Mult(HYPRE_ParVector x, HYPRE_ParVector y,
            double alpha = 1.0, double beta = 0.0);
   int MultTranspose(HypreParVector &x, HypreParVector &y,
                     double alpha = 1.0, double beta = 0.0);

   virtual void Mult(const Vector &x, Vector &y) const;
   virtual void MultTranspose(const Vector &x, Vector &y) const;

   void Print(const char *fname, int offi = 0, int offj = 0);
   void Read(const char *fname);

   virtual ~HypreParMatrix();
};

HypreParMatrix * ParMult(HypreParMatrix *A, HypreParMatrix *B);

HypreParMatrix * RAP(HypreParMatrix *A, HypreParMatrix *P);

void EliminateBC(HypreParMatrix &A, HypreParMatrix &Ae,
                 Array<int> &ess_dof_list,
                 HypreParVector &x, HypreParVector &b);

class HypreSolver : public Operator
{
protected:
   HypreParMatrix *A;

   mutable HypreParVector *B, *X;

   mutable int setup_called;

public:
   HypreSolver();

   HypreSolver(HypreParMatrix *_A);

   virtual operator HYPRE_Solver() const = 0;

   virtual HYPRE_PtrToParSolverFcn SetupFcn() const = 0;
   virtual HYPRE_PtrToParSolverFcn SolveFcn() const = 0;

   virtual void Mult(const HypreParVector &b, HypreParVector &x) const;
   virtual void Mult(const Vector &b, Vector &x) const;

   virtual ~HypreSolver();
};

class HyprePCG : public HypreSolver
{
private:
   int print_level, use_zero_initial_iterate;
   HYPRE_Solver pcg_solver;

public:
   HyprePCG(HypreParMatrix &_A);

   void SetTol(double tol);
   void SetMaxIter(int max_iter);
   void SetLogging(int logging);
   void SetPrintLevel(int print_lvl);

   void SetPreconditioner(HypreSolver &precond);

   void SetResidualConvergenceOptions(int res_frequency=-1, double rtol=0.0);

   void SetZeroInintialIterate() { use_zero_initial_iterate = 1; }

   void GetNumIterations(int &num_iterations)
   { HYPRE_ParCSRPCGGetNumIterations(pcg_solver, &num_iterations); }

   virtual operator HYPRE_Solver() const { return pcg_solver; }

   virtual HYPRE_PtrToParSolverFcn SetupFcn() const
   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRPCGSetup; }
   virtual HYPRE_PtrToParSolverFcn SolveFcn() const
   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRPCGSolve; }

   virtual void Mult(const HypreParVector &b, HypreParVector &x) const;
   using HypreSolver::Mult;

   virtual ~HyprePCG();
};

class HypreGMRES : public HypreSolver
{
private:
   int print_level, use_zero_initial_iterate;
   HYPRE_Solver gmres_solver;

public:
   HypreGMRES(HypreParMatrix &_A);

   void SetTol(double tol);
   void SetMaxIter(int max_iter);
   void SetKDim(int dim);
   void SetLogging(int logging);
   void SetPrintLevel(int print_lvl);

   void SetPreconditioner(HypreSolver &precond);

   void SetZeroInintialIterate() { use_zero_initial_iterate = 1; }

   virtual operator HYPRE_Solver() const  { return gmres_solver; }

   virtual HYPRE_PtrToParSolverFcn SetupFcn() const
   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRGMRESSetup; }
   virtual HYPRE_PtrToParSolverFcn SolveFcn() const
   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRGMRESSolve; }

   virtual void Mult (const HypreParVector &b, HypreParVector &x) const;
   using HypreSolver::Mult;

   virtual ~HypreGMRES();
};

class HypreIdentity : public HypreSolver
{
public:
   virtual operator HYPRE_Solver() const { return NULL; }

   virtual HYPRE_PtrToParSolverFcn SetupFcn() const
   { return (HYPRE_PtrToParSolverFcn) hypre_ParKrylovIdentitySetup; }
   virtual HYPRE_PtrToParSolverFcn SolveFcn() const
   { return (HYPRE_PtrToParSolverFcn) hypre_ParKrylovIdentity; }

   virtual ~HypreIdentity() { }
};

class HypreDiagScale : public HypreSolver
{
public:
   virtual operator HYPRE_Solver() const { return NULL; }

   virtual HYPRE_PtrToParSolverFcn SetupFcn() const
   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRDiagScaleSetup; }
   virtual HYPRE_PtrToParSolverFcn SolveFcn() const
   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParCSRDiagScale; }

   HypreParMatrix* GetData() { return A; }
   virtual ~HypreDiagScale() { }
};

class HypreParaSails : public HypreSolver
{
private:
   HYPRE_Solver sai_precond;

public:
   HypreParaSails(HypreParMatrix &A);

   void SetSymmetry(int sym);

   virtual operator HYPRE_Solver() const { return sai_precond; }

   virtual HYPRE_PtrToParSolverFcn SetupFcn() const
   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParaSailsSetup; }
   virtual HYPRE_PtrToParSolverFcn SolveFcn() const
   { return (HYPRE_PtrToParSolverFcn) HYPRE_ParaSailsSolve; }

   virtual ~HypreParaSails();
};

class HypreBoomerAMG : public HypreSolver
{
private:
   HYPRE_Solver amg_precond;

public:
   HypreBoomerAMG(HypreParMatrix &A);

   void SetSystemsOptions(int dim);

   virtual operator HYPRE_Solver() const { return amg_precond; }

   virtual HYPRE_PtrToParSolverFcn SetupFcn() const
   { return (HYPRE_PtrToParSolverFcn) HYPRE_BoomerAMGSetup; }
   virtual HYPRE_PtrToParSolverFcn SolveFcn() const
   { return (HYPRE_PtrToParSolverFcn) HYPRE_BoomerAMGSolve; }

   virtual ~HypreBoomerAMG();
};

class ParFiniteElementSpace;

HypreParMatrix* DiscreteGrad(ParFiniteElementSpace *edge_fespace,
                             ParFiniteElementSpace *vert_fespace);
HypreParMatrix* DiscreteCurl(ParFiniteElementSpace *face_fespace,
                             ParFiniteElementSpace *edge_fespace);

class HypreAMS : public HypreSolver
{
private:
   HYPRE_Solver ams;

   HypreParVector *x, *y, *z;
   HypreParMatrix *G;
   HypreParMatrix *Pi, *Pix, *Piy, *Piz;

public:
   HypreAMS(HypreParMatrix &A, ParFiniteElementSpace *edge_fespace);

   virtual operator HYPRE_Solver() const { return ams; }

   virtual HYPRE_PtrToParSolverFcn SetupFcn() const
   { return (HYPRE_PtrToParSolverFcn) HYPRE_AMSSetup; }
   virtual HYPRE_PtrToParSolverFcn SolveFcn() const
   { return (HYPRE_PtrToParSolverFcn) HYPRE_AMSSolve; }

   virtual ~HypreAMS();
};

class HypreADS : public HypreSolver
{
private:
   HYPRE_Solver ads;

   HypreParVector *x, *y, *z;
   HypreParMatrix *G;
   HypreParMatrix *C;
   HypreParMatrix *ND_Pi, *ND_Pix, *ND_Piy, *ND_Piz;
   HypreParMatrix *RT_Pi, *RT_Pix, *RT_Piy, *RT_Piz;

public:
   HypreADS(HypreParMatrix &A, ParFiniteElementSpace *face_fespace);

   virtual operator HYPRE_Solver() const { return ads; }

   virtual HYPRE_PtrToParSolverFcn SetupFcn() const
   { return (HYPRE_PtrToParSolverFcn) HYPRE_ADSSetup; }
   virtual HYPRE_PtrToParSolverFcn SolveFcn() const
   { return (HYPRE_PtrToParSolverFcn) HYPRE_ADSSolve; }

   virtual ~HypreADS();
};

#endif