KrylovRunner.hpp

#pragma once

#include <string>

// MTL4 headers
#include <boost/numeric/itl/itl.hpp>
#include <boost/numeric/mtl/mtl.hpp>

// AMDiS headers
#include <amdis/linearalgebra/RunnerInterface.hpp>
#include <amdis/linearalgebra/SolverInfo.hpp>

#include <amdis/linearalgebra/mtl/ITL_Preconditioner.hpp>

namespace AMDiS
{
  template <class Mat, class Vec, class ITLSolver>
  class KrylovRunner
      : public RunnerInterface<Mat,Vec>
  {
    using M = typename Mat::BaseMatrix;
    using X = typename Vec::BaseVector;
    using Y = typename Vec::BaseVector;

    using PreconBase = PreconditionerInterface<Mat,Vec>;
    using FloatType = typename Mat::value_type;

  public:
    explicit KrylovRunner(std::string const& prefix)
      : itlSolver_(prefix)
      , prefix_(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_);
    }

    void init(M const& A) override
    {
      createPrecon(prefix_);
      P_->init(A);
    }

    void exit() override
    {
      P_->exit();
    }

    int solve(M const& A, X& x, Y const& b, SolverInfo& solverInfo) override
    {
      test_exit(bool(P_), "There is no preconditioner. Call init() before solve().");

      // print information about the solution process
      itl::cyclic_iteration<FloatType> iter(two_norm(b - A*x), maxIter_, rTol_, aTol_, printCycle_);
      iter.set_quite(solverInfo.info() == 0);
      iter.suppress_resume(solverInfo.info() == 0);

      int error = itlSolver_(A, x, b, *P_, iter);
      solverInfo.setAbsResidual(iter.resid());
      solverInfo.setRelResidual(iter.relresid());

      return error;
    }

    int adjointSolve(M const& A, X& x, Y const& b, SolverInfo& solverInfo) override
    {
      test_exit(bool(P_), "There is no preconditioner. Call init() before adjointSolve().");

      auto At = trans(A);

      // print information about the solution process
      itl::cyclic_iteration<FloatType> iter(two_norm(b - At*x), maxIter_, rTol_, aTol_, printCycle_);
      iter.set_quite(solverInfo.info() == 0);
      iter.suppress_resume(solverInfo.info() == 0);

      int error = itlSolver_(At, x, b, *P_, iter);
      solverInfo.setAbsResidual(iter.resid());
      solverInfo.setRelResidual(iter.relresid());

      return error;
    }


  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>::getCreator(preconName, prefix + "->precon"));
      P_ = creator->createWithString(prefix + "->precon");
    }

  private:
    ITLSolver itlSolver_;

    std::string prefix_;

    std::shared_ptr<PreconBase> P_;

    FloatType aTol_ = 0;

    FloatType rTol_ = 1.e-6;

    std::size_t maxIter_ = 1000;

    std::size_t printCycle_ = 100;
  };

} // end namespace AMDiS