ITL_Preconditioner.hpp 4.41 KB
Newer Older
1
2
3
4
5
6
7
8
#pragma once

// MTL4 headers
#include <boost/numeric/itl/itl.hpp>
#include <boost/numeric/itl/pc/ilu_0.hpp>
#include <boost/numeric/itl/pc/ic_0.hpp>
#include <boost/numeric/mtl/vector/assigner.hpp>

9
#include <amdis/CreatorMap.hpp>
Praetorius, Simon's avatar
Praetorius, Simon committed
10
#include <amdis/linearalgebra/mtl/HyprePrecon.hpp>
11
#include <amdis/linearalgebra/mtl/Preconditioner.hpp>
12
#include <amdis/linearalgebra/mtl/SolverPrecon.hpp>
13
#include <amdis/linearalgebra/mtl/Traits.hpp>
14
15
#include <amdis/linearalgebra/mtl/itl/masslumping.hpp>

16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
namespace AMDiS
{
  /**
   * \ingroup Solver
   * \class AMDiS::DiagonalPreconditioner
   * \brief ITL_Preconditioner implementation of diagonal (jacobi) preconditioner,
   * \implements ITL_Preconditioner
   *
   * Diagonal preconditioner \f$ M^{-1} \f$ for the system \f$ Ax=b \f$ is defined as: \f$ M=diag(A) \f$.
   */
  template <class Matrix>
  using DiagonalPreconditioner = itl::pc::diagonal<Matrix>;

  /**
   * \ingroup Solver
   * \class AMDiS::DiagonalPreconditioner
   * \brief ITL_Preconditioner implementation of diagonal (jacobi) preconditioner,
   * \implements ITL_Preconditioner
   *
   * Diagonal preconditioner \f$ M^{-1} \f$ for the system \f$ Ax=b \f$ is defined as: \f$ M_ii=sum_j(A_ij) \f$.
   */
  template <class Matrix>
  using MassLumpingPreconditioner = itl::pc::masslumping<Matrix>;


  /**
   * \ingroup Solver
   * \class AMDiS::IdentityPreconditioner
   * \brief ITL_Preconditioner implementation of identity preconditioner,
   * \implements ITL_Preconditioner
   *
   * Identity preconditioner. Behaves like no preconditioning.
   */
  template <class Matrix>
  using IdentityPreconditioner = itl::pc::identity<Matrix>;


  /**
   * \ingroup Solver
   * \class AMDiS::ILUPreconditioner
   * \brief ITL_Preconditioner implementation of ILU (Incomplete LU factorization)
   * preconditioner. \implements ITL_Preconditioner
   *
   * The preconditioner is used from ITL. It corresponds for instance to
   * "Iterative Methods for Sparce Linear Systems", second edition, Yousef Saad.
   *  The preconditioner is described in chapter 10.3 (algorithm 10.4).
   */
  template <class Matrix>
  using ILUPreconditioner = itl::pc::ilu_0<Matrix>;


  /**
   * \ingroup Solver
   * \class AMDiS::ICPreconditioner
   * \brief ITL_Preconditioner implementation of IC (Incomplete Cholesky factorization)
   * preconditioner. \implements ITL_Preconditioner
   *
   * IC (Incomplete Cholesky factorization) preconditioner.
   */
  template <class Matrix>
  using ICPreconditioner = itl::pc::ic_0<Matrix>;
77

78

Praetorius, Simon's avatar
Praetorius, Simon committed
79
80
  template <class Mat, class Vec>
  class DefaultCreators< PreconditionerInterface<Mat, Vec> >
81
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
82
    using PreconBase = PreconditionerInterface<Mat, Vec>;
83

Praetorius, Simon's avatar
Praetorius, Simon committed
84
    template <template <class> class ITLPrecon>
85
    using PreconCreator
Praetorius, Simon's avatar
Praetorius, Simon committed
86
      = typename Preconditioner<Mat, Vec, ITLPrecon>::Creator;
87

88
    using Map = CreatorMap<PreconBase>;
89
90

  public:
91
92
    static void init()
    {
93
      auto pc_diag = new PreconCreator<DiagonalPreconditioner>;
94
95
      Map::addCreator("diag", pc_diag);
      Map::addCreator("jacobi", pc_diag);
96

97
      auto pc_mass = new PreconCreator<MassLumpingPreconditioner>;
98
      Map::addCreator("masslumping", pc_mass);
99

100
      auto pc_ilu = new PreconCreator<ILUPreconditioner>;
101
102
      Map::addCreator("ilu", pc_ilu);
      Map::addCreator("ilu0", pc_ilu);
103

104
      auto pc_ic = new PreconCreator<ICPreconditioner>;
105
106
      Map::addCreator("ic", pc_ic);
      Map::addCreator("ic0", pc_ic);
107

108
      auto pc_id = new PreconCreator<IdentityPreconditioner>;
109
110
      Map::addCreator("identity", pc_id);
      Map::addCreator("no", pc_id);
111

Praetorius, Simon's avatar
Praetorius, Simon committed
112
      auto pc_solver = new typename SolverPrecon<Mat,Vec>::Creator;
113
114
      Map::addCreator("solver", pc_solver);

Praetorius, Simon's avatar
Praetorius, Simon committed
115
      Map::addCreator("default", pc_id);
Praetorius, Simon's avatar
Praetorius, Simon committed
116
117

      init_hypre(std::is_same<typename Dune::FieldTraits<typename Mat::value_type>::real_type, double>{});
Praetorius, Simon's avatar
Praetorius, Simon committed
118
119
120
121
122
123
    }

    static void init_hypre(std::false_type) {}
    static void init_hypre(std::true_type)
    {
#if AMDIS_HAS_HYPRE && HAVE_MPI
Praetorius, Simon's avatar
Praetorius, Simon committed
124
      auto pc_hypre = new typename HyprePrecon<Mat,Vec>::Creator;
125
126
      Map::addCreator("hypre", pc_hypre);
#endif
127
128
    }
  };
129

Praetorius, Simon's avatar
Praetorius, Simon committed
130
131
132
  template <class Mat, class Vec, class X>
  itl::pc::solver<PreconditionerInterface<Mat,Vec>, X, false>
  solve(PreconditionerInterface<Mat,Vec> const& P, X const& vin)
133
134
135
136
  {
    return {P, vin};
  }

Praetorius, Simon's avatar
Praetorius, Simon committed
137
138
139
  template <class Mat, class Vec, class X>
  itl::pc::solver<PreconditionerInterface<Mat,Vec>, X, true>
  adjoint_solve(PreconditionerInterface<Mat,Vec> const& P, X const& vin)
140
141
142
  {
    return {P, vin};
  }
143

144
} // namespace AMDiS