IterativeRunner.hpp 1.75 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#pragma once

#include <dune/istl/preconditioner.hh>

#include <amdis/linear_algebra/RunnerInterface.hpp>
#include <amdis/linear_algebra/SolverInfo.hpp>
#include <amdis/linear_algebra/eigen/PreconConfig.hpp>
#include <amdis/linear_algebra/eigen/SolverConfig.hpp>

namespace AMDiS
{
  template <class IterativeSolver, class Matrix, class VectorX, class VectorB>
  class IterativeRunner
      : public RunnerInterface<Matrix, VectorX, VectorB>
  {
    using SolverCfg = SolverConfig<IterativeSolver>;
    using PreconCfg = PreconConfig<typename IterativeSolver::Preconditioner>;

  public:
Praetorius, Simon's avatar
Praetorius, Simon committed
20
    IterativeRunner(std::string const& prefix)
21
22
23
24
25
26
27
28
29
30
31
      : solver_{}
    {
      SolverCfg::init(prefix, solver_);
      PreconCfg::init(prefix + "->precon", solver_.preconditioner());
    }


    /// Implementes \ref RunnerInterface::init()
    virtual void init(Matrix const& A) override
    {
      solver_.compute(A);
Praetorius, Simon's avatar
Praetorius, Simon committed
32
33
34

      test_exit(solver_.info() == Eigen::Success,
        "Error in solver.compute(matrix)");
35
36
37
38
39
40
41
42
43
    }

    /// Implementes \ref RunnerInterface::exit()
    virtual void exit() override {}

    /// Implementes \ref RunnerInterface::solve()
    virtual int solve(Matrix const& A, VectorX& x, VectorB const& b,
                      SolverInfo& solverInfo) override
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
44
      solver_.setTolerance(solverInfo.relTolerance());
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
      x = solver_.solveWithGuess(b, x);

      auto r = VectorB(b);
      if (x.norm() != 0)
        r -= A * x;

      solverInfo.setAbsResidual(r.norm());
      solverInfo.setRelResidual(solver_.error());
      solverInfo.setError(solver_.info());
      msg("number of iteration: {}", solver_.iterations());

      return solver_.info() == Eigen::Success ? 0 : int(solver_.info());
    }

  private:
    IterativeSolver solver_;
  };

} // end namespace AMDiS