HyprePrecon.hpp

#pragma once
 
#if AMDIS_HAS_HYPRE && HAVE_MPI
 
#include <HYPRE.h>
 
#include <amdis/Initfile.hpp>
#include <amdis/linearalgebra/mtl/PreconditionerInterface.hpp>
#include <amdis/linearalgebra/mtl/itl/hypre.hpp>
 
namespace AMDiS
{
 
  template <class M, class X, class Y>
  class HyprePrecon
      : public PreconditionerInterface<M,X,Y>
  {
    using Self = HyprePrecon;
    using Super = PreconditionerInterface<M,X,Y>;
 
  public:
    struct Creator : CreatorInterfaceName<Super>
    {
      std::unique_ptr<Super> createWithString(std::string prefix) override
      {
        return std::make_unique<Self>(prefix);
      }
    };
 
  public:
    HyprePrecon(std::string const& prefix)
    {
      // read HYPRE parameters
      int maxIter = 1, cycleMode = -1, interpolation = -1, relaxation = -1, info = 1;
      Parameters::get(prefix + "->max iteration", maxIter);
      Parameters::get(prefix + "->cycle mode", cycleMode);
      Parameters::get(prefix + "->interpolation type", interpolation);
      Parameters::get(prefix + "->relax type", relaxation);
      Parameters::get(prefix + "->info", info);
 
      config_.maxIter(maxIter);
      config_.interpolationType(interpolation);
      config_.relaxType(relaxation);
      config_.cycle(cycleMode);
      config_.printLevel(info);
    }
 
    void init(M const& mat) override
    {
      hypreMatrix_.init(mat);
      HYPRE_IJMatrixGetObject(hypreMatrix_, (void**)(&matrix_));
      HYPRE_BoomerAMGCreate(&solver_);
 
      // apply solver parameters
      config_(solver_);
 
      mtl::dense_vector<double> swap(1, 0.0);
      mtl::HypreParVector x(swap);
      HYPRE_BoomerAMGSetup(solver_, matrix_, x, x);
 
      solverCreated_ = true;
    }
 
    void finish() override
    {
      if (solverCreated_)
        HYPRE_BoomerAMGDestroy(solver_);
      solverCreated_ = false;
    }
 
    void solve(X const& in, Y& out) const override
    {
      assert(solverCreated_);
      mtl::HypreParVector x(in); // use b as initial guess
      mtl::HypreParVector b(in);
 
      [[maybe_unused]] int error = HYPRE_BoomerAMGSolve(solver_, matrix_, b, x);
      assert(error != 0);
 
      // write output back to MTL vector
      mtl::convert(x.getHypreVector(), out);
    }
 
    void adjoint_solve(X const& in, Y& out) const override
    {
      assert(solverCreated_);
      mtl::HypreParVector x(in); // use b as initial guess
      mtl::HypreParVector b(in);
 
      [[maybe_unused]] int error = HYPRE_BoomerAMGSolveT(solver_, matrix_, b, x);
      assert(error != 0);
 
      // write output back to MTL vector
      mtl::convert(x.getHypreVector(), out);
    }
 
  private:
    mtl::HypreMatrix hypreMatrix_;
    HYPRE_ParCSRMatrix matrix_;
    HYPRE_Solver solver_;
 
    mtl::AMGConfigurator config_;
    bool solverCreated_ = false;
  };
 
} // end namespace AMDiS
 
#endif // AMDIS_HAS_HYPRE && HAVE_MPI