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

#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 "sparsemat.hpp"

namespace mfem
{

class ParFiniteElementSpace;
class HypreParMatrix;

class HypreParVector : public Vector
{
private:
   int own_ParVector;

   hypre_ParVector *x;

   friend class HypreParMatrix;

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

   operator hypre_ParVector*() const;
#ifndef HYPRE_PAR_VECTOR_STRUCT
   operator HYPRE_ParVector() const;
#endif
   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)
   { height = GetNumRows(); width = GetNumCols(); X = Y = 0; CommPkg = 0; }
   HypreParMatrix(MPI_Comm comm, int glob_size, int *row_starts,
                  SparseMatrix *diag);
   HypreParMatrix(MPI_Comm comm, int global_num_rows, int global_num_cols,
                  int *row_starts, int *col_starts, SparseMatrix *diag);
   HypreParMatrix(MPI_Comm comm, int global_num_rows, int global_num_cols,
                  int *row_starts, int *col_starts,
                  SparseMatrix *diag, SparseMatrix *offd, int *cmap);

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

   HypreParMatrix(MPI_Comm comm, int global_num_rows, int global_num_cols,
                  int *row_starts, int *col_starts, 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();

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

   operator hypre_ParCSRMatrix*();
#ifndef HYPRE_PAR_CSR_MATRIX_STRUCT
   operator HYPRE_ParCSRMatrix();
#endif
   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 GetNumCols() const
   { return hypre_CSRMatrixNumCols(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);

   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 ScaleRows(const Vector & s);
   void InvScaleRows(const Vector & s);
   void operator*=(double s);

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

   virtual ~HypreParMatrix();
};

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,
                 Array<int> &ess_dof_list,
                 HypreParVector &x, HypreParVector &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:
   int print_level;
   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_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;
   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;

public:
   HypreBoomerAMG(HypreParMatrix &A);

   void SetSystemsOptions(int dim);

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

   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 // MFEM_USE_MPI

#endif