KrylovRunner.hpp

#pragma once
 
#include <string>
 
// MTL4 headers
#include <boost/numeric/itl/itl.hpp>
#include <boost/numeric/mtl/mtl.hpp>
 
// AMDiS headers
#include <amdis/Environment.hpp>
#include <amdis/linearalgebra/LinearSolverInterface.hpp>
#include <amdis/linearalgebra/mtl/Preconditioners.hpp>
 
namespace AMDiS
{
  template <class M, class X, class Y, class ITLSolver>
  class KrylovRunner
      : public LinearSolverInterface<M,X,Y>
  {
    using Self = KrylovRunner;
 
    using PreconBase = PreconditionerInterface<M,X,Y>;
    using FloatType = typename M::value_type;
 
  public:
    struct Creator final : CreatorInterfaceName<LinearSolverInterface<M,X,Y>>
    {
      std::unique_ptr<LinearSolverInterface<M,X,Y>>
      createWithString(std::string prefix) final
      {
        return std::make_unique<Self>(prefix);
      }
    };
 
  public:
    explicit KrylovRunner(std::string const& prefix)
      : itlSolver_(prefix)
    {
      Parameters::get(prefix + "->absolute tolerance", aTol_);
      Parameters::get(prefix + "->relative tolerance", rTol_);
      Parameters::get(prefix + "->max iteration", maxIter_);
      Parameters::get(prefix + "->print cycle", printCycle_);
 
      createPrecon(prefix);
    }
 
    void init(M const& A) override
    {
      P_->init(A);
      A_ = &A;
    }
 
    void finish() override
    {
      P_->finish();
    }
 
    void apply(X& x, Y const& b, Dune::InverseOperatorResult& stat) override
    {
      Dune::Timer t;
 
      // print information about the solution process
      FloatType r0 = two_norm(b - (*A_)*x);
      itl::cyclic_iteration<FloatType> iter(r0, maxIter_, rTol_, aTol_, printCycle_);
      iter.set_quite(Environment::infoLevel() == 0);
      iter.suppress_resume(Environment::infoLevel() == 0);
 
      [[maybe_unused]] int error = itlSolver_(*A_, x, b, *P_, iter);
 
      stat.iterations = iter.iterations();
      stat.converged = iter.is_converged();
      stat.reduction = iter.relresid();
      stat.elapsed = t.elapsed();
    }
 
  protected:
    void createPrecon(std::string const& prefix)
    {
      // Creator for the left preconditioner
      std::string preconName = "default";
      Parameters::get(prefix + "->precon", preconName);
 
      auto creator
        = named(CreatorMap<PreconBase>::get(preconName, prefix + "->precon"));
      P_ = creator->createWithString(prefix + "->precon");
    }
 
  private:
    ITLSolver itlSolver_;
 
    std::string prefix_;
 
    std::shared_ptr<PreconBase> P_;
 
    M const* A_ = nullptr;
 
    FloatType aTol_ = 0;
 
    FloatType rTol_ = 1.e-6;
 
    std::size_t maxIter_ = 1000;
 
    std::size_t printCycle_ = 100;
  };
 
} // end namespace AMDiS