KrylovRunner.hpp 3.83 KB
Newer Older
1
2
#pragma once

3
4
#include <string>

5
6
7
8
9
// MTL4 headers
#include <boost/numeric/itl/itl.hpp>
#include <boost/numeric/mtl/mtl.hpp>

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

13
#include <amdis/linearalgebra/mtl/ITL_Preconditioner.hpp>
14
15
16
17
18
19
20

namespace AMDiS
{
  /** \ingroup Solver
   *
   * \brief
   * Wrapper class for different MTL4 itl-solvers. These solvers
21
   * are parametrized by Matrix and Vector.
22
   **/
Praetorius, Simon's avatar
Praetorius, Simon committed
23
  template <class Mat, class Vec, class ITLSolver>
24
  class KrylovRunner
Praetorius, Simon's avatar
Praetorius, Simon committed
25
      : public RunnerInterface<Mat,Vec>
26
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
27
28
29
30
31
32
    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;
33

34
35
  public:
    /// Constructor.
Praetorius, Simon's avatar
Praetorius, Simon committed
36
    explicit KrylovRunner(std::string const& prefix)
37
      : itlSolver_(prefix)
38
      , prefix_(prefix)
39
    {
40
41
42
43
      Parameters::get(prefix_ + "->absolute tolerance", aTol_);
      Parameters::get(prefix_ + "->relative tolerance", rTol_);
      Parameters::get(prefix_ + "->max iteration", maxIter_);
      Parameters::get(prefix_ + "->print cycle", printCycle_);
44
45
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
46
47
    /// Implementation of \ref RunnerInterface::initImpl()
    void init(M const& A) override
48
    {
49
      createPrecon(prefix_);
50
      P_->init(A);
51
    }
52

53
    /// Implementation of \ref RunnerInterface::exit()
54
    void exit() override
55
    {
56
      P_->exit();
57
58
59
    }

    /// Implementation of \ref RunnerInterface::solve()
Praetorius, Simon's avatar
Praetorius, Simon committed
60
    int solve(M const& A, X& x, Y const& b, SolverInfo& solverInfo) override
61
    {
62
      test_exit(bool(P_), "There is no preconditioner. Call init() before solve().");
63

64
      // print information about the solution process
Praetorius, Simon's avatar
Praetorius, Simon committed
65
      itl::cyclic_iteration<FloatType> iter(two_norm(b - A*x), maxIter_, rTol_, aTol_, printCycle_);
Praetorius, Simon's avatar
Praetorius, Simon committed
66
67
      iter.set_quite(solverInfo.info() == 0);
      iter.suppress_resume(solverInfo.info() == 0);
68

69
      int error = itlSolver_(A, x, b, *P_, iter);
70
71
      solverInfo.setAbsResidual(iter.resid());
      solverInfo.setRelResidual(iter.relresid());
72
73
74
75

      return error;
    }

Praetorius, Simon's avatar
Praetorius, Simon committed
76
77
78
    /// Implementation of \ref RunnerInterface::solve()
    int adjointSolve(M const& A, X& x, Y const& b, SolverInfo& solverInfo) override
    {
79
      test_exit(bool(P_), "There is no preconditioner. Call init() before adjointSolve().");
Praetorius, Simon's avatar
Praetorius, Simon committed
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95

      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;
    }


96
  protected:
97
    /// Create left/right preconditioners from parameters given in the init-file
98
    void createPrecon(std::string const& prefix)
99
    {
100
      // Creator for the left preconditioner
101
102
      std::string preconName = "default";
      Parameters::get(prefix + "->precon", preconName);
103

104
105
106
      auto creator
        = named(CreatorMap<PreconBase>::getCreator(preconName, prefix + "->precon"));
      P_ = creator->createWithString(prefix + "->precon");
107
108
    }

109
  private:
110
    /// The itl solver object that performes the actual solve
111
    ITLSolver itlSolver_;
112

113
114
115
116
117
118
119
120
121
122
123
    /// The parameter prefix
    std::string prefix_;

    /// Pointer to the left/right preconditioner
    std::shared_ptr<PreconBase> P_;

    /// Absolute solver tolerance ||b - A*x|| < aTol
    FloatType aTol_ = 0;

    /// Relative solver tolerance ||b - A*x||/||b - A*x0|| < rTol
    FloatType rTol_ = 1.e-6;
124

125
    /// The maximal number of iterations
126
    std::size_t maxIter_ = 1000;
127

128
    /// Print solver residuals every \ref printCycle 's iteration
129
    std::size_t printCycle_ = 100;
130
131
132
  };

} // end namespace AMDiS