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.org.
//
// 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

#include "../config/config.hpp"

#ifdef MFEM_USE_MPI

#include <mpi.h>

// 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"
#include "temp_multivector.h"

#ifdef HYPRE_COMPLEX
#error "MFEM does not work with HYPRE's complex numbers support"
#endif

#include "sparsemat.hpp"
#include "hypre_parcsr.hpp"

namespace mfem
{

class ParFiniteElementSpace;
class HypreParMatrix;

namespace internal
{

// Convert a HYPRE_Int to int
inline int to_int(HYPRE_Int i)
{
#ifdef HYPRE_BIGINT
   MFEM_ASSERT(HYPRE_Int(int(i)) == i, "overflow converting HYPRE_Int to int");
#endif
   return int(i);
}

}

class HypreParVector : public Vector
{
private:
   int own_ParVector;

   hypre_ParVector *x;

   friend class HypreParMatrix;

   // Set Vector::data and Vector::size from *x
   inline void _SetDataAndSize_();

public:
   HypreParVector(MPI_Comm comm, HYPRE_Int glob_size, HYPRE_Int *col);
   HypreParVector(MPI_Comm comm, HYPRE_Int glob_size, double *_data,
                  HYPRE_Int *col);
   HypreParVector(const HypreParVector &y);
   explicit HypreParVector(const HypreParMatrix &A, int transpose = 0);
   explicit HypreParVector(HYPRE_ParVector y);
   explicit HypreParVector(ParFiniteElementSpace *pfes);

   MPI_Comm GetComm() { return x->comm; }

   inline HYPRE_Int *Partitioning() { return x->partitioning; }

   inline HYPRE_Int GlobalSize() { return x->global_size; }

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

   void SetOwnership(int own) { own_ParVector = own; }

   int GetOwnership() const { return own_ParVector; }

   Vector* GlobalVector() const;

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

   void SetData(double *_data);

   HYPRE_Int Randomize(HYPRE_Int seed);

   void Print(const char *fname) const;

   ~HypreParVector();
};

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


double ParNormlp(const Vector &vec, double p, MPI_Comm comm);


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

   mutable HypreParVector *X, *Y;

   // Flags indicating ownership of A->diag->{i,j,data}, A->offd->{i,j,data},
   // and A->col_map_offd.
   // The possible values for diagOwner are:
   //  -1: no special treatment of A->diag (default)
   //   0: prevent hypre from destroying A->diag->{i,j,data}
   //   1: same as 0, plus take ownership of A->diag->{i,j}
   //   2: same as 0, plus take ownership of A->diag->data
   //   3: same as 0, plus take ownership of A->diag->{i,j,data}
   // The same values and rules apply to offdOwner and A->offd.
   // The possible values for colMapOwner are:
   //  -1: no special treatment of A->col_map_offd (default)
   //   0: prevent hypre from destroying A->col_map_offd
   //   1: same as 0, plus take ownership of A->col_map_offd
   // All owned arrays are destroyed with 'delete []'.
   char diagOwner, offdOwner, colMapOwner;

   // Does the object own the pointer A?
   char ParCSROwner;

   // Initialize with defaults. Does not initialize inherited members.
   void Init();

   // Delete all owned data. Does not perform re-initialization with defaults.
   void Destroy();

   // Copy (shallow or deep, based on HYPRE_BIGINT) the I and J arrays from csr
   // to hypre_csr. Shallow copy the data. Return the appropriate ownership
   // flag.
   static char CopyCSR(SparseMatrix *csr, hypre_CSRMatrix *hypre_csr);
   // Copy (shallow or deep, based on HYPRE_BIGINT) the I and J arrays from
   // bool_csr to hypre_csr. Allocate the data array and set it to all ones.
   // Return the appropriate ownership flag.
   static char CopyBoolCSR(Table *bool_csr, hypre_CSRMatrix *hypre_csr);

   // Copy the j array of a hypre_CSRMatrix to the given J array, converting
   // the indices from HYPRE_Int to int.
   static void CopyCSR_J(hypre_CSRMatrix *hypre_csr, int *J);

public:
   HypreParMatrix();

   explicit HypreParMatrix(hypre_ParCSRMatrix *a)
   {
      Init();
      A = a;
      height = GetNumRows();
      width = GetNumCols();
   }

   HypreParMatrix(MPI_Comm comm, HYPRE_Int glob_size, HYPRE_Int *row_starts,
                  SparseMatrix *diag);

   HypreParMatrix(MPI_Comm comm, HYPRE_Int global_num_rows,
                  HYPRE_Int global_num_cols, HYPRE_Int *row_starts,
                  HYPRE_Int *col_starts, SparseMatrix *diag);

   HypreParMatrix(MPI_Comm comm, HYPRE_Int global_num_rows,
                  HYPRE_Int global_num_cols, HYPRE_Int *row_starts,
                  HYPRE_Int *col_starts, SparseMatrix *diag, SparseMatrix *offd,
                  HYPRE_Int *cmap);

   HypreParMatrix(MPI_Comm comm,
                  HYPRE_Int global_num_rows, HYPRE_Int global_num_cols,
                  HYPRE_Int *row_starts, HYPRE_Int *col_starts,
                  HYPRE_Int *diag_i, HYPRE_Int *diag_j, double *diag_data,
                  HYPRE_Int *offd_i, HYPRE_Int *offd_j, double *offd_data,
                  HYPRE_Int offd_num_cols, HYPRE_Int *offd_col_map);

   HypreParMatrix(MPI_Comm comm, HYPRE_Int *row_starts, HYPRE_Int *col_starts,
                  SparseMatrix *a);

   HypreParMatrix(MPI_Comm comm, HYPRE_Int global_num_rows,
                  HYPRE_Int global_num_cols, HYPRE_Int *row_starts,
                  HYPRE_Int *col_starts, Table *diag);

   HypreParMatrix(MPI_Comm comm, int id, int np, HYPRE_Int *row, HYPRE_Int *col,
                  HYPRE_Int *i_diag, HYPRE_Int *j_diag, HYPRE_Int *i_offd,
                  HYPRE_Int *j_offd, HYPRE_Int *cmap, HYPRE_Int cmap_size);

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

   void MakeRef(const HypreParMatrix &master);

   MPI_Comm GetComm() const { return A->comm; }

   operator hypre_ParCSRMatrix*() { return A; }
#ifndef HYPRE_PAR_CSR_MATRIX_STRUCT
   operator HYPRE_ParCSRMatrix() { return (HYPRE_ParCSRMatrix) A; }
#endif
   hypre_ParCSRMatrix* StealData();

   void SetOwnerFlags(char diag, char offd, char colmap)
   { diagOwner = diag, offdOwner = offd, colMapOwner = colmap; }

   char OwnsDiag() const { return diagOwner; }
   char OwnsOffd() const { return offdOwner; }
   char OwnsColMap() const { return colMapOwner; }

   void CopyRowStarts();
   void CopyColStarts();

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

   void GetDiag(Vector &diag) const;
   void GetDiag(SparseMatrix &diag) const;
   void GetOffd(SparseMatrix &offd, HYPRE_Int* &cmap) const;

   void GetBlocks(Array2D<HypreParMatrix*> &blocks,
                  bool interleaved_rows = false,
                  bool interleaved_cols = false) const;

   HypreParMatrix * Transpose();

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

   int GetNumCols() const
   {
      return internal::to_int(
                hypre_CSRMatrixNumCols(hypre_ParCSRMatrixDiag(A)));
   }

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

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

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

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

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

   void Mult(double a, const Vector &x, double b, Vector &y) const;
   void MultTranspose(double a, const Vector &x, double b, Vector &y) const;

   virtual void Mult(const Vector &x, Vector &y) const
   { Mult(1.0, x, 0.0, y); }
   virtual void MultTranspose(const Vector &x, Vector &y) const
   { MultTranspose(1.0, x, 0.0, y); }

   void BooleanMult(int alpha, int *x, int beta, int *y)
   {
      internal::hypre_ParCSRMatrixBooleanMatvec(A, alpha, x, beta, y);
   }

   HypreParMatrix* LeftDiagMult(const SparseMatrix &D,
                                HYPRE_Int* row_starts = NULL) const;

   void ScaleRows(const Vector & s);
   void InvScaleRows(const Vector & s);
   void operator*=(double s);

   void Threshold(double threshold = 0.0);

   void EliminateZeroRows() { hypre_ParCSRMatrixFixZeroRows(A); }

   void EliminateRowsCols(const Array<int> &rows_cols, const HypreParVector &X,
                          HypreParVector &B);

   HypreParMatrix* EliminateRowsCols(const Array<int> &rows_cols);

   void Print(const char *fname, HYPRE_Int offi = 0, HYPRE_Int offj = 0);
   void Read(MPI_Comm comm, const char *fname);

   virtual ~HypreParMatrix() { Destroy(); }
};

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

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

void EliminateBC(HypreParMatrix &A, HypreParMatrix &Ae,
                 const Array<int> &ess_dof_list, const Vector &X, Vector &B);


class HypreSmoother : public Solver
{
protected:
   HypreParMatrix *A;
   mutable HypreParVector *B, *X;
   mutable HypreParVector *V, *Z;
   mutable HypreParVector *X0, *X1;

   int type;
   int relax_times;
   double relax_weight;
   double omega;
   int poly_order;
   double poly_fraction;
   int poly_scale;

   double lambda;
   double mu;
   int taubin_iter;

   double *l1_norms;
   double max_eig_est;
   double min_eig_est;
   double window_params[3];

   double* fir_coeffs;

public:
   enum Type { Jacobi = 0, l1Jacobi = 1, l1GS = 2, l1GStr = 4, lumpedJacobi = 5,
               GS = 6, Chebyshev = 16, Taubin = 1001, FIR = 1002
             };

   HypreSmoother();

   HypreSmoother(HypreParMatrix &_A, int type = l1GS,
                 int relax_times = 1, double relax_weight = 1.0,
                 double omega = 1.0, int poly_order = 2,
                 double poly_fraction = .3);

   void SetType(HypreSmoother::Type type, int relax_times = 1);
   void SetSOROptions(double relax_weight, double omega);
   void SetPolyOptions(int poly_order, double poly_fraction);
   void SetTaubinOptions(double lambda, double mu, int iter);

   void SetWindowByName(const char* window_name);
   void SetWindowParameters(double a, double b, double c);
   void SetFIRCoefficients(double max_eig);

   virtual void SetOperator(const Operator &op);

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

   virtual ~HypreSmoother();
};


class HypreSolver : public Solver
{
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 SetOperator(const Operator &op)
   { mfem_error("HypreSolvers do not support SetOperator!"); }

   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:
   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() { iterative_mode = false; }

   void GetNumIterations(int &num_iterations)
   {
      HYPRE_Int num_it;
      HYPRE_ParCSRPCGGetNumIterations(pcg_solver, &num_it);
      num_iterations = internal::to_int(num_it);
   }

   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:
   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() { iterative_mode = false; }

   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:
   HypreDiagScale() : HypreSolver() { }
   explicit HypreDiagScale(HypreParMatrix &A) : HypreSolver(&A) { }
   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;

   Array<HYPRE_ParVector> rbms;

   ParFiniteElementSpace *fespace;

   void RecomputeRBMs();

   void SetDefaultOptions();

   // If amg_precond is NULL, allocates it and sets default options.
   // Otherwise saves the options from amg_precond, destroys it, allocates a new
   // one, and sets its options to the saved values.
   void ResetAMGPrecond();

public:
   HypreBoomerAMG();

   HypreBoomerAMG(HypreParMatrix &A);

   virtual void SetOperator(const Operator &op);

   void SetSystemsOptions(int dim);

   void SetElasticityOptions(ParFiniteElementSpace *fespace);

   void SetPrintLevel(int print_level)
   { HYPRE_BoomerAMGSetPrintLevel(amg_precond, print_level); }

   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();
};

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);

   void SetPrintLevel(int print_lvl);

   void SetSingularProblem() { HYPRE_AMSSetBetaPoissonMatrix(ams, NULL); }

   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);

   void SetPrintLevel(int print_lvl);

   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();
};

class HypreLOBPCG
{
private:
   MPI_Comm comm;
   int myid;
   int numProcs;
   int nev;   // Number of desired eigenmodes
   int seed;  // Random seed used for initial vectors

   HYPRE_Int glbSize; // Global number of DoFs in the linear system
   HYPRE_Int * part;  // Row partitioning of the linear system

   // Pointer to HYPRE's solver struct
   HYPRE_Solver lobpcg_solver;

   // Interface for matrix storage type
   mv_InterfaceInterpreter interpreter;

   // Interface for setting up and performing matrix-vector products
   HYPRE_MatvecFunctions matvec_fn;

   // Eigenvalues
   Array<double> eigenvalues;

   // Forward declaration
   class HypreMultiVector;

   // MultiVector to store eigenvectors
   HypreMultiVector * multi_vec;

   // Empty vectors used to setup the matrices and preconditioner
   HypreParVector * x;

   class HypreMultiVector
   {
   private:
      // Pointer to hypre's multi-vector object
      mv_MultiVectorPtr mv_ptr;

      // Wrappers for each member of the multivector
      HypreParVector ** hpv;

      // Number of vectors in the multivector
      int nv;

   public:
      HypreMultiVector(int n, HypreParVector & v,
                       mv_InterfaceInterpreter & interpreter);
      ~HypreMultiVector();

      void Randomize(HYPRE_Int seed);

      HypreParVector & GetVector(unsigned int i);

      HypreParVector ** StealVectors();

      operator mv_MultiVectorPtr() const { return mv_ptr; }

      mv_MultiVectorPtr & GetMultiVector() { return mv_ptr; }
   };

   static void    * OperatorMatvecCreate( void *A, void *x );
   static HYPRE_Int OperatorMatvec( void *matvec_data,
                                    HYPRE_Complex alpha,
                                    void *A,
                                    void *x,
                                    HYPRE_Complex beta,
                                    void *y );
   static HYPRE_Int OperatorMatvecDestroy( void *matvec_data );

   static HYPRE_Int PrecondSolve(void *solver,
                                 void *A,
                                 void *b,
                                 void *x);
   static HYPRE_Int PrecondSetup(void *solver,
                                 void *A,
                                 void *b,
                                 void *x);

public:
   HypreLOBPCG(MPI_Comm comm);
   ~HypreLOBPCG();

   void SetTol(double tol);
   void SetMaxIter(int max_iter);
   void SetPrintLevel(int logging);
   void SetNumModes(int num_eigs) { nev = num_eigs; }
   void SetPrecondUsageMode(int pcg_mode);
   void SetRandomSeed(int s) { seed = s; }

   // The following four methods support general operators
   void SetPreconditioner(Solver & precond);
   void SetOperator(Operator & A);
   void SetMassMatrix(Operator & M);

   void Solve();

   void GetEigenvalues(Array<double> & eigenvalues);

   HypreParVector & GetEigenvector(unsigned int i);

   HypreParVector ** StealEigenvectors() { return multi_vec->StealVectors(); }
};

class HypreAME
{
private:
   int myid;
   int numProcs;
   int nev;   // Number of desired eigenmodes
   bool setT;

   // Pointer to HYPRE's AME solver struct
   HYPRE_Solver ame_solver;

   // Pointer to HYPRE's AMS solver struct
   HypreSolver * ams_precond;

   // Eigenvalues
   HYPRE_Real * eigenvalues;

   // MultiVector to store eigenvectors
   HYPRE_ParVector * multi_vec;

   HypreParVector ** eigenvectors;

   void createDummyVectors();

public:
   HypreAME(MPI_Comm comm);
   ~HypreAME();

   void SetTol(double tol);
   void SetMaxIter(int max_iter);
   void SetPrintLevel(int logging);
   void SetNumModes(int num_eigs);

   // The following four methods support operators of type HypreParMatrix.
   void SetPreconditioner(HypreSolver & precond);
   void SetOperator(HypreParMatrix & A);
   void SetMassMatrix(HypreParMatrix & M);

   void Solve();

   void GetEigenvalues(Array<double> & eigenvalues);

   HypreParVector & GetEigenvector(unsigned int i);

   HypreParVector ** StealEigenvectors();
};

}

#endif // MFEM_USE_MPI

#endif