SolverCreator.hpp 5.89 KB
Newer Older
1
2
3
4
5
6
7
#pragma once

#include <memory>

#include <Eigen/IterativeLinearSolvers>
#include <Eigen/MetisSupport>
#include <Eigen/SparseLU>
Praetorius, Simon's avatar
Praetorius, Simon committed
8
#include <Eigen/SuperLUSupport>
9
#include <Eigen/UmfPackSupport>
Praetorius, Simon's avatar
Praetorius, Simon committed
10
#include <unsupported/Eigen/IterativeSolvers>
11
12
13

#include <amdis/CreatorMap.hpp>
#include <amdis/Initfile.hpp>
14
#include <amdis/Output.hpp>
15

16
17
18
#include <amdis/linearalgebra/LinearSolver.hpp>
#include <amdis/linearalgebra/eigen/DirectRunner.hpp>
#include <amdis/linearalgebra/eigen/IterativeRunner.hpp>
19
#include <amdis/linearalgebra/eigen/Traits.hpp>
20
21
22

namespace AMDiS
{
Praetorius, Simon's avatar
Praetorius, Simon committed
23
  template <class Mat, class Vec, template <class> class IterativeSolver>
24
  struct IterativeSolverCreator
Praetorius, Simon's avatar
Praetorius, Simon committed
25
      : public CreatorInterfaceName<LinearSolverInterface<Mat,Vec>>
26
27
28
  {
    template <class Precon>
    using SolverCreator
Praetorius, Simon's avatar
Praetorius, Simon committed
29
      = typename LinearSolver<Mat,Vec, IterativeRunner<Mat, Vec, IterativeSolver<Precon>>>::Creator;
30

Praetorius, Simon's avatar
Praetorius, Simon committed
31
32
33
    using SolverBase = LinearSolverInterface<Mat,Vec>;
    using M = typename Mat::BaseMatrix;
    using Scalar = typename M::Scalar;
34

35
    std::unique_ptr<SolverBase> createWithString(std::string prefix) override
36
37
    {
      // get creator string for preconditioner
Praetorius, Simon's avatar
Praetorius, Simon committed
38
39
      std::string precon = "no";
      Parameters::get(prefix + "->precon->name", precon);
40

Praetorius, Simon's avatar
Praetorius, Simon committed
41
42
      if (precon == "diag" ||
          precon == "jacobi")
43
44
      {
        auto creator = SolverCreator<Eigen::DiagonalPreconditioner<Scalar>>{};
45
        return creator.createWithString(prefix);
46
      }
Praetorius, Simon's avatar
Praetorius, Simon committed
47
      else if (precon == "ilu")
48
49
      {
        auto creator = SolverCreator<Eigen::IncompleteLUT<Scalar>>{};
50
        return creator.createWithString(prefix);
51
      }
Praetorius, Simon's avatar
Praetorius, Simon committed
52
53
54
55
      else if (precon == "ic")
      {
        return createIncompleteCholesky(prefix);
      }
56
57
      else {
        auto creator = SolverCreator<Eigen::IdentityPreconditioner>{};
58
        return creator.createWithString(prefix);
59
60
      }
    }
Praetorius, Simon's avatar
Praetorius, Simon committed
61
62
63
64
65
66


    template <class Ordering>
    using IncompleteCholesky =
      SolverCreator<Eigen::IncompleteCholesky<Scalar, Eigen::Lower|Eigen::Upper, Ordering>>;

67
    std::unique_ptr<SolverBase> createIncompleteCholesky(std::string prefix) const
Praetorius, Simon's avatar
Praetorius, Simon committed
68
69
70
71
72
    {
      std::string ordering = "amd";
      Parameters::get(prefix + "->precon->ordering", ordering);

      if (ordering == "amd") {
Praetorius, Simon's avatar
Praetorius, Simon committed
73
        using AMD = Eigen::AMDOrdering<typename M::StorageIndex>;
74
        return IncompleteCholesky<AMD>{}.createWithString(prefix);
Praetorius, Simon's avatar
Praetorius, Simon committed
75
76
      }
      else {
Praetorius, Simon's avatar
Praetorius, Simon committed
77
        using NATURAL = Eigen::NaturalOrdering<typename M::StorageIndex>;
78
        return IncompleteCholesky<NATURAL>{}.createWithString(prefix);
Praetorius, Simon's avatar
Praetorius, Simon committed
79
80
      }
    }
81
82
83
84
85
86
87
  };

  /// Adds default creators for linear solvers based on `Eigen::SparseMatrix`.
  /**
   * Adds creators for full-matrix aware solvers.
   * - *cg*: conjugate gradient method, \see Eigen::ConjugateGradient
   * - *bcgs*: stabilized bi-conjugate gradient method, \see Eigen::BiCGSTAB
Praetorius, Simon's avatar
Praetorius, Simon committed
88
89
90
   * - *minres*: Minimal Residual Algorithm (for symmetric matrices), \see Eigen::MINRES
   * - *gmres*: Generalized Minimal Residual Algorithm, \see Eigen::GMRES
   * - *dgmres*: stabilized bi-conjugate gradient method, \see Eigen::DGMRES
91
92
93
   * - *umfpack*: external UMFPACK solver, \see Eigen::UmfPackLU
   * - *superlu*: external SparseLU solver, \see Eigen::SparseLU
   **/
Praetorius, Simon's avatar
Praetorius, Simon committed
94
95
  template <class Mat, class Vec>
  class DefaultCreators< LinearSolverInterface<Mat,Vec> >
96
  {
Praetorius, Simon's avatar
Praetorius, Simon committed
97
98
    using SolverBase = LinearSolverInterface<Mat,Vec>;
    using M = typename Mat::BaseMatrix;
99
100
101

    template <template <class> class IterativeSolver>
    using SolverCreator
Praetorius, Simon's avatar
Praetorius, Simon committed
102
103
104
105
106
107
108
109
      = IterativeSolverCreator<Mat, Vec, IterativeSolver>;

    template <template <class> class DirectSolver>
    struct RunnerCreator
    {
      template <class M, class V>
      using Impl = DirectRunner<M, V, DirectSolver>;
    };
110

111
    template <template <class> class DirectSolver>
Praetorius, Simon's avatar
Praetorius, Simon committed
112
    using DirectSolverCreator
Praetorius, Simon's avatar
Praetorius, Simon committed
113
      = typename LinearSolver<Mat, Vec, DirectRunner<Mat, Vec, DirectSolver>>::Creator;
Praetorius, Simon's avatar
Praetorius, Simon committed
114
115

    //@{
116
    template <class Precon>
Praetorius, Simon's avatar
Praetorius, Simon committed
117
    using CG = Eigen::ConjugateGradient<M, Eigen::Lower|Eigen::Upper, Precon>;
118
119

    template <class Precon>
Praetorius, Simon's avatar
Praetorius, Simon committed
120
    using BCGS = Eigen::BiCGSTAB<M, Precon>;
121

Praetorius, Simon's avatar
Praetorius, Simon committed
122
    template <class Precon>
Praetorius, Simon's avatar
Praetorius, Simon committed
123
    using MINRES = Eigen::MINRES<M, Eigen::Lower|Eigen::Upper, Precon>;
124

Praetorius, Simon's avatar
Praetorius, Simon committed
125
    template <class Precon>
Praetorius, Simon's avatar
Praetorius, Simon committed
126
    using GMRES = Eigen::GMRES<M, Precon>;
Praetorius, Simon's avatar
Praetorius, Simon committed
127
128

    template <class Precon>
Praetorius, Simon's avatar
Praetorius, Simon committed
129
    using DGMRES = Eigen::DGMRES<M, Precon>;
Praetorius, Simon's avatar
Praetorius, Simon committed
130
    // @}
131
132
133
134
135
136

    using Map = CreatorMap<SolverBase>;

  public:
    static void init()
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
137
      // conjugate gradient solver
138
139
140
      auto cg = new SolverCreator<CG>;
      Map::addCreator("cg", cg);

Praetorius, Simon's avatar
Praetorius, Simon committed
141
      // bi-conjugate gradient stabilized solver
142
143
144
145
      auto bicgstab = new SolverCreator<BCGS>;
      Map::addCreator("bicgstab", bicgstab);
      Map::addCreator("bcgs", bicgstab);

Praetorius, Simon's avatar
Praetorius, Simon committed
146
147
148
149
150
151
152
153
154
155
156
157
      // Minimal Residual Algorithm (for symmetric matrices)
      auto minres = new SolverCreator<MINRES>;
      Map::addCreator("minres", minres);

      // Generalized Minimal Residual Algorithm
      auto gmres = new SolverCreator<GMRES>;
      Map::addCreator("gmres", gmres);

      // Restarted GMRES with deflation.
      auto dgmres = new SolverCreator<DGMRES>;
      Map::addCreator("dgmres", dgmres);

158
159
160
      // default iterative solver
      Map::addCreator("default", gmres);

Praetorius, Simon's avatar
Praetorius, Simon committed
161
      init_direct(std::is_same<typename Dune::FieldTraits<typename M::Scalar>::real_type, double>{});
162
163
164
165
166
    }

    static void init_direct(std::false_type) {}
    static void init_direct(std::true_type)
    {
167
#if HAVE_SUITESPARSE_UMFPACK
Praetorius, Simon's avatar
Praetorius, Simon committed
168
      // sparse LU factorization and solver based on UmfPack
169
      auto umfpack = new DirectSolverCreator<Eigen::UmfPackLU>;
170
171
172
      Map::addCreator("umfpack", umfpack);
#endif

Praetorius, Simon's avatar
Praetorius, Simon committed
173
174
#if HAVE_SUPERLU
      // sparse direct LU factorization and solver based on the SuperLU library
175
      auto superlu = new DirectSolverCreator<Eigen::SuperLU>;
Praetorius, Simon's avatar
Praetorius, Simon committed
176
      Map::addCreator("superlu", superlu);
177
178
179
180
181
#endif

      // default direct solvers
#if HAVE_SUITESPARSE_UMFPACK
      Map::addCreator("direct", umfpack);
Praetorius, Simon's avatar
Praetorius, Simon committed
182
183
#elif HAVE_SUPERLU
      Map::addCreator("direct", superlu);
184
185
186
187
188
#endif
    }
  };

} // end namespace AMDiS