SolverCreator.hpp 5.64 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
14

#include <amdis/CreatorMap.hpp>
#include <amdis/Initfile.hpp>

15
16
17
#include <amdis/linearalgebra/LinearSolver.hpp>
#include <amdis/linearalgebra/eigen/DirectRunner.hpp>
#include <amdis/linearalgebra/eigen/IterativeRunner.hpp>
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32

namespace AMDiS
{
  template <class Matrix, class VectorX, class VectorB, template <class> class IterativeSolver>
  struct IterativeSolverCreator
      : public CreatorInterfaceName<LinearSolverInterface<Matrix, VectorX, VectorB>>
  {
    template <class Precon>
    using SolverCreator
      = typename LinearSolver<Matrix, VectorX, VectorB,
          IterativeRunner<IterativeSolver<Precon>, Matrix, VectorX, VectorB>>::Creator;

    using SolverBase = LinearSolverInterface<Matrix, VectorX, VectorB>;
    using Scalar = typename Matrix::Scalar;

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

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


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

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

      if (ordering == "amd") {
        using AMD = Eigen::AMDOrdering<typename Matrix::StorageIndex>;
        return IncompleteCholesky<AMD>{}.create(prefix);
      }
      else {
        using NATURAL = Eigen::NaturalOrdering<typename Matrix::StorageIndex>;
        return IncompleteCholesky<NATURAL>{}.create(prefix);
      }
    }
79
80
81
82
83
84
85
86
  };


  /// 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
87
88
89
   * - *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
90
91
92
93
94
95
96
97
98
99
100
101
102
   * - *umfpack*: external UMFPACK solver, \see Eigen::UmfPackLU
   * - *superlu*: external SparseLU solver, \see Eigen::SparseLU
   **/
  template <class T, int O, class VectorX, class VectorB>
  class DefaultCreators<LinearSolverInterface<Eigen::SparseMatrix<T,O>, VectorX, VectorB>>
  {
    using Matrix = Eigen::SparseMatrix<T,O>;
    using SolverBase = LinearSolverInterface<Matrix, VectorX, VectorB>;

    template <template <class> class IterativeSolver>
    using SolverCreator
      = IterativeSolverCreator<Matrix, VectorX, VectorB, IterativeSolver>;

Praetorius, Simon's avatar
Praetorius, Simon committed
103
104
105
106
107
108
    template <class DirectSolver>
    using DirectSolverCreator
      = typename LinearSolver<Matrix, VectorX, VectorB,
          DirectRunner<DirectSolver, Matrix, VectorX, VectorB>>::Creator;

    //@{
109
110
111
112
113
114
    template <class Precon>
    using CG = Eigen::ConjugateGradient<Matrix, Eigen::Lower|Eigen::Upper, Precon>;

    template <class Precon>
    using BCGS = Eigen::BiCGSTAB<Matrix, Precon>;

Praetorius, Simon's avatar
Praetorius, Simon committed
115
116
    template <class Precon>
    using MINRES = Eigen::MINRES<Matrix, Eigen::Lower|Eigen::Upper, Precon>;
117

Praetorius, Simon's avatar
Praetorius, Simon committed
118
119
120
121
122
123
    template <class Precon>
    using GMRES = Eigen::GMRES<Matrix, Precon>;

    template <class Precon>
    using DGMRES = Eigen::DGMRES<Matrix, Precon>;
    // @}
124
125
126
127
128
129

    using Map = CreatorMap<SolverBase>;

  public:
    static void init()
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
130
      // conjugate gradient solver
131
132
133
      auto cg = new SolverCreator<CG>;
      Map::addCreator("cg", cg);

Praetorius, Simon's avatar
Praetorius, Simon committed
134
      // bi-conjugate gradient stabilized solver
135
136
137
138
      auto bicgstab = new SolverCreator<BCGS>;
      Map::addCreator("bicgstab", bicgstab);
      Map::addCreator("bcgs", bicgstab);

Praetorius, Simon's avatar
Praetorius, Simon committed
139
140
141
142
143
144
145
146
147
148
149
150
      // 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);

151
#if HAVE_SUITESPARSE_UMFPACK
Praetorius, Simon's avatar
Praetorius, Simon committed
152
      // sparse LU factorization and solver based on UmfPack
153
154
155
156
      auto umfpack = new DirectSolverCreator<Eigen::UmfPackLU<Matrix>>;
      Map::addCreator("umfpack", umfpack);
#endif

Praetorius, Simon's avatar
Praetorius, Simon committed
157
158
159
160
#if HAVE_SUPERLU
      // sparse direct LU factorization and solver based on the SuperLU library
      auto superlu = new DirectSolverCreator<Eigen::SuperLU<Matrix>>;
      Map::addCreator("superlu", superlu);
161
162
163
#endif

      // default iterative solver
Praetorius, Simon's avatar
Praetorius, Simon committed
164
      Map::addCreator("default", gmres);
165
166
167
168

      // default direct solvers
#if HAVE_SUITESPARSE_UMFPACK
      Map::addCreator("direct", umfpack);
Praetorius, Simon's avatar
Praetorius, Simon committed
169
170
#elif HAVE_SUPERLU
      Map::addCreator("direct", superlu);
171
172
173
174
175
#endif
    }
  };

} // end namespace AMDiS