SolverPrecon.hpp

#pragma once

#include <memory>

#include <dune/istl/solver.hh>

#include <amdis/CreatorMap.hpp>
#include <amdis/Initfile.hpp>
#include <amdis/Output.hpp>
#include <amdis/linearalgebra/LinearSolverInterface.hpp>
#include <amdis/linearalgebra/mtl/PreconditionerInterface.hpp>

namespace AMDiS
{
  template <class M, class X, class Y>
  class SolverPrecon
      : public PreconditionerInterface<M,X,Y>
  {
    using Self = SolverPrecon;
    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>(std::move(prefix));
      }
    };

  public:
    SolverPrecon(std::string const& prefix)
    {
      std::string solverName = "default";
      Parameters::get(prefix + "->solver", solverName);

      auto solverCreator = named(CreatorMap<LinearSolverInterface<M,X,Y>>::get(solverName, prefix + "->solver"));
      solver_ = solverCreator->createWithString(prefix + "->solver");
    }

    void init(M const& matrix) override
    {
      matrix_ = &matrix;
      solver_->init(matrix);
      stat_.clear();
    }

    void finish() override
    {
      matrix_ = nullptr;
    }

    void solve(X const& x, Y& y) const override
    {
      test_exit_dbg(bool(solver_), "No solver initialized!");
      test_exit_dbg(bool(matrix_), "No matrix initialized!");

      y.checked_change_resource(x);
      test_exit(size(y) == num_cols(*matrix_), "incompatible size");
      solver_->apply(y, x, stat_);
    }

    void adjoint_solve(X const& x, Y& y) const override
    {
      error_exit("Not Implemented.");
    }

  private:
    mutable Dune::InverseOperatorResult stat_;

    std::shared_ptr<LinearSolverInterface<M,X,Y>> solver_;
    M const* matrix_ = nullptr;
  };

} // namespace AMDiS