hypre.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/base_header.hpp>
#include <kernel/solver/adp_solver_base.hpp>
 
#if defined(FEAT_HAVE_HYPRE) || defined(DOXYGEN)
 
// includes, system
#include <vector>
 
namespace FEAT
{
  namespace Solver
  {
    namespace Hypre
    {
      typedef double HypreDataType;
 
#ifdef FEAT_TPL_HYPRE_INT64
      typedef long long HypreIndexType;
#else
      typedef int HypreIndexType;
#endif
 
      void* create_core(const void* comm, Index dof_offset, Index num_owned_dofs, Index num_nonzeros);
 
      void destroy_core(void* core);
 
      HypreIndexType* get_row_ptr(void* core);
      HypreIndexType* get_col_idx(void* core);
      HypreDataType* get_mat_val(void* core);
 
      HypreDataType* get_vec_def(void* core);
      HypreDataType* get_vec_cor(void* core);
 
      void upload_symbolic(void* core);
      void upload_mat_val(void* core);
      void upload_vec_def(void* core);
      void download_vec_cor(void* core);
      void format_vec_cor(void* core);
 
      // ParaSails Wrappers
      void* create_parasails(void* core, int* iparam, double* dparam);
      void destroy_parasails(void* solver);
      void solve_parasails(void* core, void* solver);
 
      // Euclid Wrappers
      void* create_euclid(void* core, int* iparam, double* dparam);
      void destroy_euclid(void* solver);
      void solve_euclid(void* core, void* solver);
 
      // BoomerAMG Wrappers
      void* create_boomeramg(void* core, int* iparam, double* dparam);
      void destroy_boomeramg(void* solver);
      void solve_boomeramg(void* core, void* solver);
    } // namespace Hypre
 
    template<typename Matrix_, typename Filter_, typename SolverBase_ = Solver::SolverBase<typename Matrix_::VectorTypeL>>
    class HypreSolverBase :
      public ADPSolverBase<Matrix_, Filter_, SolverBase_>
    {
    public:
      typedef ADPSolverBase<Matrix_, Filter_, SolverBase_> BaseClass;
      // the vector type
      typedef typename BaseClass::VectorType VectorType;
 
    protected:
      void* _core;
 
      explicit HypreSolverBase(const Matrix_& matrix, const Filter_& filter) :
        BaseClass(matrix, filter),
        _core(nullptr)
      {
      }
 
      void _upload_def(const VectorType& vec_def)
      {
        this->_upload_vector(Hypre::get_vec_def(this->_core), vec_def.local());
        Hypre::upload_vec_def(this->_core);
      }
 
      void _download_cor(VectorType& vec_cor)
      {
        Hypre::download_vec_cor(this->_core);
        this->_download_vector(vec_cor.local(), Hypre::get_vec_cor(this->_core));
 
        // apply correction filter
        this->_system_filter.filter_cor(vec_cor);
      }
 
      void _format_cor()
      {
        Hypre::format_vec_cor(this->_core);
      }
 
    public:
      virtual void init_symbolic() override
      {
        BaseClass::init_symbolic();
 
        XASSERT(this->_core == nullptr);
 
        // create our HYPRE core wrapper object
        this->_core = Hypre::create_core(
#ifdef FEAT_HAVE_MPI
          &this->_get_comm()->mpi_comm(),
#else
          nullptr, // no communicator for non-MPI builds
#endif
          this->_get_global_dof_offset(), this->_get_num_owned_dofs(), this->_get_adp_matrix_num_nzes());
 
 
        // upload matrix structure
        BaseClass::_upload_symbolic(Hypre::get_row_ptr(this->_core), Hypre::get_col_idx(this->_core));
 
        // assemble system
        Hypre::upload_symbolic(this->_core);
      }
 
      virtual void init_numeric() override
      {
        BaseClass::init_numeric();
 
        XASSERT(this->_core != nullptr);
 
        // upload matrix values
        this->_upload_numeric(Hypre::get_mat_val(this->_core),
          Hypre::get_row_ptr(this->_core), Hypre::get_col_idx(this->_core));
 
        Hypre::upload_mat_val(this->_core);
      }
 
      virtual void done_symbolic() override
      {
        XASSERT(this->_core != nullptr);
 
        Hypre::destroy_core(this->_core);
        this->_core = nullptr;
 
        BaseClass::done_symbolic();
      }
    }; // class HypreSolverBase<...>
 
    template<typename Matrix_, typename Filter_>
    class ParaSailsPrecond :
      public HypreSolverBase<Matrix_, Filter_>
    {
    public:
      typedef HypreSolverBase<Matrix_, Filter_> BaseClass;
      typedef typename BaseClass::VectorType VectorType;
 
    protected:
      void* _solver;
 
      // 0: num levels, third argument for HYPRE_ParaSailsSetParams, default: 1
      // 1: symmetry, second argument for HYPRE_ParaSailsSetSym, default: 2
      int _iparam[2];
 
      // 0: threshold, second argument for HYPRE_ParaSailsSetParams, default: 0.1
      // 1: filter, second argument for HYPRE_ParaSailsSetFilter, default: 0.05
      double _dparam[2];
 
    public:
      explicit ParaSailsPrecond(const Matrix_& matrix, const Filter_& filter) :
        BaseClass(matrix, filter),
        _solver(nullptr)
      {
        _iparam[0] = 1;
        _iparam[1] = 2;
 
        _dparam[0] = 0.1;
        _dparam[1] = 0.05;
      }
 
      explicit ParaSailsPrecond(const String& section_name, const PropertyMap* section,
        const Matrix_& matrix, const Filter_& filter)
         :
        ParaSailsPrecond(matrix, filter)
      {
        auto level_p = section->query("level");
        if(level_p.second && !level_p.first.parse(_iparam[0]))
          throw ParseError(section_name + ".level", level_p.first, "a non-negative integer");
 
        auto sym_p = section->query("sym");
        if(sym_p.second && !sym_p.first.parse(_iparam[1]))
          throw ParseError(section_name + ".sym", sym_p.first, "a non-negative integer");
 
        auto thresh_p = section->query("thresh");
        if(thresh_p.second && !thresh_p.first.parse(_dparam[0]))
          throw ParseError(section_name + ".thresh", thresh_p.first, "a positive float");
 
        auto filter_p = section->query("filter");
        if(filter_p.second && !filter_p.first.parse(_dparam[1]))
          throw ParseError(section_name + ".filter", filter_p.first, "a positive float");
      }
 
      virtual String name() const override
      {
        return "ParaSailsPrecond";
      }
 
      virtual void init_numeric() override
      {
        BaseClass::init_numeric();
 
        // create ParaSails preconditioner
        this->_solver = Hypre::create_parasails(this->_core, this->_iparam, this->_dparam);
      }
 
      virtual void done_numeric() override
      {
        if(_solver != nullptr)
        {
          Hypre::destroy_parasails(_solver);
          _solver = nullptr;
        }
 
        BaseClass::done_numeric();
      }
 
      virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) override
      {
        // upload defect vector and format correction
        this->_upload_def(vec_def);
        this->_format_cor();
 
        // apply ParaSails preconditioner
        Hypre::solve_parasails(this->_core, this->_solver);
 
        // download correction
        this->_download_cor(vec_cor);
 
        // okay
        return Status::success;
      }
    }; // class ParaSailsPrecond
 
    template<typename Matrix_, typename Filter_>
    inline std::shared_ptr<ParaSailsPrecond<Matrix_, Filter_>> new_parasails_precond(
      const Matrix_& matrix, const Filter_& filter)
    {
      return std::make_shared<ParaSailsPrecond<Matrix_, Filter_>>(matrix, filter);
    }
 
    template<typename Matrix_, typename Filter_>
    inline std::shared_ptr<ParaSailsPrecond<Matrix_, Filter_>> new_parasails_precond(
      const String& section_name, const PropertyMap* section,
      const Matrix_& matrix, const Filter_& filter)
    {
      return std::make_shared<ParaSailsPrecond<Matrix_, Filter_>>(section_name, section, matrix, filter);
    }
 
    template<typename Matrix_, typename Filter_>
    class EuclidPrecond :
      public HypreSolverBase<Matrix_, Filter_>
    {
    public:
      typedef HypreSolverBase<Matrix_, Filter_> BaseClass;
      typedef typename BaseClass::VectorType VectorType;
 
    protected:
      void* _solver;
 
      // 0: factorization level for ILU(k), default: 1
      int _iparam[1];
 
      // 0: drop tolerance for ILU(k), default: 0.0
      double _dparam[1];
 
    public:
      explicit EuclidPrecond(const Matrix_& matrix, const Filter_& filter) :
        BaseClass(matrix, filter),
        _solver(nullptr)
      {
        _iparam[0] = 1;
        _dparam[0] = 0.0;
      }
 
      explicit EuclidPrecond(const String& section_name, const PropertyMap* section,
        const Matrix_& matrix, const Filter_& filter)
         :
        EuclidPrecond(matrix, filter)
      {
        auto level_p = section->query("level");
        if(level_p.second && !level_p.first.parse(_iparam[0]))
          throw ParseError(section_name + ".level", level_p.first, "a non-negative integer");
 
        auto drop_p = section->query("drop");
        if(drop_p.second && !drop_p.first.parse(_dparam[0]))
          throw ParseError(section_name + ".drop", drop_p.first, "a non-negative float");
      }
 
      virtual String name() const override
      {
        return "EuclidPrecond";
      }
 
      virtual void init_numeric() override
      {
        BaseClass::init_numeric();
 
        // create Euclid preconditioner
        this->_solver = Hypre::create_euclid(this->_core, this->_iparam, this->_dparam);
      }
 
      virtual void done_numeric() override
      {
        if(_solver != nullptr)
        {
          Hypre::destroy_euclid(_solver);
          _solver = nullptr;
        }
      }
 
      virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) override
      {
        // upload defect vector and format correction
        this->_upload_def(vec_def);
        this->_format_cor();
 
        // apply Euclid preconditioner
        Hypre::solve_euclid(this->_core, this->_solver);
 
        // download correction
        this->_download_cor(vec_cor);
 
        // okay
        return Status::success;
      }
    }; // class EuclidPrecond
 
    template<typename Matrix_, typename Filter_>
    inline std::shared_ptr<EuclidPrecond<Matrix_, Filter_>> new_euclid_precond(
      const Matrix_& matrix, const Filter_& filter)
    {
      return std::make_shared<EuclidPrecond<Matrix_, Filter_>>(matrix, filter);
    }
 
    template<typename Matrix_, typename Filter_>
    inline std::shared_ptr<EuclidPrecond<Matrix_, Filter_>> new_euclid_precond(
      const String& section_name, const PropertyMap* section,
      const Matrix_& matrix, const Filter_& filter)
    {
      return std::make_shared<EuclidPrecond<Matrix_, Filter_>>(section_name, section, matrix, filter);
    }
 
    template<typename Matrix_, typename Filter_>
    class BoomerAMG :
      public HypreSolverBase<Matrix_, Filter_>
    {
    public:
      typedef HypreSolverBase<Matrix_, Filter_> BaseClass;
      typedef typename BaseClass::VectorType VectorType;
 
    protected:
      void* _solver;
 
      int _iparam[1];
 
      double _dparam[1];
 
    public:
      explicit BoomerAMG(const Matrix_& matrix, const Filter_& filter) :
        BaseClass(matrix, filter),
        _solver(nullptr)
      {
        // set tolerance to 0 and maximum iterations to 1, so that BoomerAMG
        // acts as a preconditioner rather than a real iterative solver
        _iparam[0] = 1;
        _dparam[0] = 0.0;
      }
 
      explicit BoomerAMG(const String& DOXY(section_name), const PropertyMap* DOXY(section),
        const Matrix_& matrix, const Filter_& filter)
         :
        BoomerAMG(matrix, filter)
      {
      }
 
      virtual String name() const override
      {
        return "BoomerAMG";
      }
 
      virtual void init_numeric() override
      {
        BaseClass::init_numeric();
 
        // create BoomerAMG preconditioner
        this->_solver = Hypre::create_boomeramg(this->_core, this->_iparam, this->_dparam);
      }
 
      virtual void done_numeric() override
      {
        if(_solver != nullptr)
        {
          Hypre::destroy_boomeramg(_solver);
          _solver = nullptr;
        }
      }
 
      virtual Status apply(VectorType& vec_cor, const VectorType& vec_def) override
      {
        // upload defect vector and format correction
        this->_upload_def(vec_def);
        this->_format_cor();
 
        // apply BoomerAMG preconditioner
        Hypre::solve_boomeramg(this->_core, this->_solver);
 
        // download correction
        this->_download_cor(vec_cor);
 
        // okay
        return Status::success;
      }
    }; // class BoomerAMG
 
    template<typename Matrix_, typename Filter_>
    inline std::shared_ptr<BoomerAMG<Matrix_, Filter_>> new_boomeramg(
      const Matrix_& matrix, const Filter_& filter)
    {
      return std::make_shared<BoomerAMG<Matrix_, Filter_>>(matrix, filter);
    }
 
    template<typename Matrix_, typename Filter_>
    inline std::shared_ptr<BoomerAMG<Matrix_, Filter_>> new_boomeramg(
      const String& section_name, const PropertyMap* section,
      const Matrix_& matrix, const Filter_& filter)
    {
      return std::make_shared<BoomerAMG<Matrix_, Filter_>>(section_name, section, matrix, filter);
    }
  } // namespace Solver
} // namespace FEAT
 
#endif // defined(FEAT_HAVE_HYPRE) || defined(DOXYGEN)