BlockJacobi.hpp 2.05 KB
Newer Older
Praetorius, Simon's avatar
Praetorius, Simon committed
1
2
3
4
#pragma once

#include <dune/common/tuplevector.hh>

5
6
#include <dune/dec/common/Compatible.hpp>
#include <dune/dec/linear_algebra/block.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
7
8
9

namespace Dec
{
10
11
12
13
14
15
16
17
18
19
20
21
22
  /**
    * \addtogroup linear_solver
    * @{
    **/

  /// \brief Block-Jacobi method, can be used as solver and preconditioner
  /**
   * The (variadic) template parameters `Solver...` must model `LinearSolver`, i.e.
   * must provide the methods
   * - `compute()` (initialization of the solver)
   * - `solve()` (solve with initial solution u=0)
   * - `solveWithGuess()` (solve with provided initial solution)
   **/
23
  template <class... Solver>
Praetorius, Simon's avatar
Praetorius, Simon committed
24
25
  struct BlockJacobi
  {
26
27
    static constexpr int N = sizeof...(Solver);

Praetorius, Simon's avatar
Praetorius, Simon committed
28
29
    using Matrix = BlockMatrix<float_type, N, N>;
    using Vector = BlockVector<float_type, N>;
30

Praetorius, Simon's avatar
Praetorius, Simon committed
31
    /// Constructor, stores references to all solvers
32
    BlockJacobi(Solver&... solver)
Praetorius, Simon's avatar
Praetorius, Simon committed
33
34
      : solver_(&solver...)
    {}
35

36
    /// Calls `compute` on all solvers, with the diagonal matrix \f$ A_{ii}\f$ and stores
Praetorius, Simon's avatar
Praetorius, Simon committed
37
    /// pointer to the matrix.
Praetorius, Simon's avatar
Praetorius, Simon committed
38
39
    BlockJacobi& compute(Matrix const& mat)
    {
40
41
      forEach(range_<0,N>, [this,&mat](auto const _i)
      {
Praetorius, Simon's avatar
Praetorius, Simon committed
42
43
44
        this->solver_[_i]->compute(mat(_i,_i));
      });
      mat_ = &mat;
45

Praetorius, Simon's avatar
Praetorius, Simon committed
46
47
48
      return *this;
    }

49
    /// Solve \f$ A_{ii} u_i = b_i\f$
Praetorius, Simon's avatar
Praetorius, Simon committed
50
51
    void solve(Vector const& b, Vector& u) const
    {
52
53
54
      forEach(range_<0,N>, [this,&b,&u](auto const _i)
      {
        this->solver_[_i]->solve(b[_i], u[_i]);
Praetorius, Simon's avatar
Praetorius, Simon committed
55
56
      });
    }
57

58
    /// Solve \f$ A_{ii} d = b_i - (A\cdot u)_i,\; u += d\f$
Praetorius, Simon's avatar
Praetorius, Simon committed
59
60
    void solveWithGuess(Vector const& b, Vector& u) const
    {
61
62
      assert( mat_ != nullptr && "Call BlockJacobi::compute() first" );

Praetorius, Simon's avatar
Praetorius, Simon committed
63
64
      Vector r(b);
      r -= (*mat_)*u;
65

Praetorius, Simon's avatar
Praetorius, Simon committed
66
      Vector d(u, tag::ressource{});
Praetorius, Simon's avatar
Praetorius, Simon committed
67
      solve(r, d);
68

Praetorius, Simon's avatar
Praetorius, Simon committed
69
70
      u += d;
    }
71

Praetorius, Simon's avatar
Praetorius, Simon committed
72
  private:
73
74
75
    Dune::TupleVector<Solver*...> solver_;
    Matrix const* mat_ = nullptr;

Praetorius, Simon's avatar
Praetorius, Simon committed
76
77
  };

78
  /// Generator for Block-Jacobi solver, \relates BlockJacobi
79
80
81
82
83
84
  template <class... Solver>
  BlockJacobi<Solver...> makeBlockJacobi(Solver&... solver)
  {
    return {solver...};
  }

85
86
  /** @} */

Praetorius, Simon's avatar
Praetorius, Simon committed
87
} // end namespace Dec