OGS
EigenLinearSolver.cpp
Go to the documentation of this file.
1 
11 #include "EigenLinearSolver.h"
12 
13 #include <Eigen/Sparse>
14 
15 #include "BaseLib/Logging.h"
16 
17 #ifdef USE_MKL
18 #include <Eigen/PardisoSupport>
19 #endif
20 
21 #ifdef USE_EIGEN_UNSUPPORTED
22 #include <unsupported/Eigen/src/IterativeSolvers/GMRES.h>
23 #include <unsupported/Eigen/src/IterativeSolvers/Scaling.h>
24 #endif
25 
26 #include "BaseLib/ConfigTree.h"
27 #include "EigenMatrix.h"
28 #include "EigenTools.h"
29 #include "EigenVector.h"
31 
32 namespace MathLib
33 {
34 // TODO change to LinearSolver
36 {
37 public:
40 
41  virtual ~EigenLinearSolverBase() = default;
42 
44  virtual bool solve(Matrix& A, Vector const& b, Vector& x,
45  EigenOption& opt) = 0;
46 };
47 
48 namespace details
49 {
51 template <class T_SOLVER>
53 {
54 public:
55  bool solve(Matrix& A, Vector const& b, Vector& x, EigenOption& opt) override
56  {
57  INFO("-> solve with {:s}", EigenOption::getSolverName(opt.solver_type));
58  if (!A.isCompressed())
59  {
60  A.makeCompressed();
61  }
62 
63  solver_.compute(A);
64  if (solver_.info() != Eigen::Success)
65  {
66  ERR("Failed during Eigen linear solver initialization");
67  return false;
68  }
69 
70  x = solver_.solve(b);
71  if (solver_.info() != Eigen::Success)
72  {
73  ERR("Failed during Eigen linear solve");
74  return false;
75  }
76 
77  return true;
78  }
79 
80 private:
81  T_SOLVER solver_;
82 };
83 
85 template <class T_SOLVER>
87 {
88 public:
89  bool solve(Matrix& A, Vector const& b, Vector& x, EigenOption& opt) override
90  {
91  INFO("-> solve with {:s} (precon {:s})",
94  solver_.setTolerance(opt.error_tolerance);
95  solver_.setMaxIterations(opt.max_iterations);
97  opt.restart);
98 
99  if (!A.isCompressed())
100  {
101  A.makeCompressed();
102  }
103 
104  solver_.compute(A);
105  if (solver_.info() != Eigen::Success)
106  {
107  ERR("Failed during Eigen linear solver initialization");
108  return false;
109  }
110 
111  x = solver_.solveWithGuess(b, x);
112  INFO("\t iteration: {:d}/{:d}", solver_.iterations(),
113  opt.max_iterations);
114  INFO("\t residual: {:e}\n", solver_.error());
115 
116  if (solver_.info() != Eigen::Success)
117  {
118  ERR("Failed during Eigen linear solve");
119  return false;
120  }
121 
122  return true;
123  }
124 
125 private:
126  T_SOLVER solver_;
127  void setRestart(int const /*restart*/) {}
128 };
129 
131 template <>
132 void EigenIterativeLinearSolver<
133  Eigen::GMRES<EigenMatrix::RawMatrixType,
134  Eigen::IdentityPreconditioner>>::setRestart(int const restart)
135 {
136  solver_.set_restart(restart);
137  INFO("-> set restart value: {:d}", solver_.get_restart());
138 }
139 
140 template <>
141 void EigenIterativeLinearSolver<Eigen::GMRES<
143  Eigen::DiagonalPreconditioner<double>>>::setRestart(int const restart)
144 {
145  solver_.set_restart(restart);
146  INFO("-> set restart value: {:d}", solver_.get_restart());
147 }
148 
149 template <>
151  Eigen::GMRES<EigenMatrix::RawMatrixType,
152  Eigen::IncompleteLUT<double>>>::setRestart(int const restart)
153 {
154  solver_.set_restart(restart);
155  INFO("-> set restart value: {:d}", solver_.get_restart());
156 }
157 
158 template <template <typename, typename> class Solver, typename Precon>
159 std::unique_ptr<EigenLinearSolverBase> createIterativeSolver()
160 {
161  using Slv =
163  return std::make_unique<Slv>();
164 }
165 
166 template <template <typename, typename> class Solver>
167 std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
168  EigenOption::PreconType precon_type)
169 {
170  switch (precon_type)
171  {
173  return createIterativeSolver<Solver,
174  Eigen::IdentityPreconditioner>();
176  return createIterativeSolver<
177  Solver, Eigen::DiagonalPreconditioner<double>>();
179  // TODO for this preconditioner further options can be passed.
180  // see
181  // https://eigen.tuxfamily.org/dox/classEigen_1_1IncompleteLUT.html
182  return createIterativeSolver<Solver,
183  Eigen::IncompleteLUT<double>>();
184  default:
185  OGS_FATAL("Invalid Eigen preconditioner type.");
186  }
187 }
188 
189 template <typename Mat, typename Precon>
190 using EigenCGSolver = Eigen::ConjugateGradient<Mat, Eigen::Lower, Precon>;
191 
192 std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
193  EigenOption::SolverType solver_type, EigenOption::PreconType precon_type)
194 {
195  switch (solver_type)
196  {
198  {
199  return createIterativeSolver<Eigen::BiCGSTAB>(precon_type);
200  }
202  {
203  return createIterativeSolver<EigenCGSolver>(precon_type);
204  }
206  {
207 #ifdef USE_EIGEN_UNSUPPORTED
208  return createIterativeSolver<Eigen::GMRES>(precon_type);
209 #else
210  OGS_FATAL(
211  "The code is not compiled with the Eigen unsupported modules. "
212  "Linear solver type GMRES is not available.");
213 #endif
214  }
215  default:
216  OGS_FATAL("Invalid Eigen iterative linear solver type. Aborting.");
217  }
218 }
219 
220 } // namespace details
221 
222 EigenLinearSolver::EigenLinearSolver(const std::string& /*solver_name*/,
223  const BaseLib::ConfigTree* const option)
224 {
225  using Matrix = EigenMatrix::RawMatrixType;
226 
227  if (option)
228  {
229  setOption(*option);
230  }
231 
232  // TODO for my taste it is much too unobvious that the default solver type
233  // currently is SparseLU.
234  switch (option_.solver_type)
235  {
237  {
238  using SolverType =
239  Eigen::SparseLU<Matrix, Eigen::COLAMDOrdering<int>>;
240  solver_ = std::make_unique<
242  return;
243  }
249  return;
251  {
252 #ifdef USE_MKL
253  using SolverType = Eigen::PardisoLU<EigenMatrix::RawMatrixType>;
255  return;
256 #else
257  OGS_FATAL(
258  "The code is not compiled with Intel MKL. Linear solver type "
259  "PardisoLU is not available.");
260 #endif
261  }
262  }
263 
264  OGS_FATAL("Invalid Eigen linear solver type. Aborting.");
265 }
266 
268 
270 {
271  ignoreOtherLinearSolvers(option, "eigen");
273  auto const ptSolver = option.getConfigSubtreeOptional("eigen");
274  if (!ptSolver)
275  {
276  return;
277  }
278 
279  if (auto solver_type =
281  ptSolver->getConfigParameterOptional<std::string>("solver_type"))
282  {
284  }
285  if (auto precon_type =
287  ptSolver->getConfigParameterOptional<std::string>("precon_type"))
288  {
290  }
291  if (auto error_tolerance =
293  ptSolver->getConfigParameterOptional<double>("error_tolerance"))
294  {
295  option_.error_tolerance = *error_tolerance;
296  }
297  if (auto max_iteration_step =
299  ptSolver->getConfigParameterOptional<int>("max_iteration_step"))
300  {
301  option_.max_iterations = *max_iteration_step;
302  }
303  if (auto scaling =
305  ptSolver->getConfigParameterOptional<bool>("scaling"))
306  {
307 #ifdef USE_EIGEN_UNSUPPORTED
308  option_.scaling = *scaling;
309 #else
310  OGS_FATAL(
311  "The code is not compiled with the Eigen unsupported modules. "
312  "scaling is not available.");
313 #endif
314  }
315  if (auto restart =
317  ptSolver->getConfigParameterOptional<int>("restart"))
318  {
319 #ifdef USE_EIGEN_UNSUPPORTED
320  option_.restart = *restart;
321 #else
322  OGS_FATAL(
323  "The code is not compiled with the Eigen unsupported modules. "
324  "GMRES/GMRES option restart is not available.");
325 #endif
326  }
327 }
328 
330 {
331  INFO("------------------------------------------------------------------");
332  INFO("*** Eigen solver computation");
333 
334 #ifdef USE_EIGEN_UNSUPPORTED
335  std::unique_ptr<Eigen::IterScaling<EigenMatrix::RawMatrixType>> scal;
336  if (option_.scaling)
337  {
338  INFO("-> scale");
339  scal =
340  std::make_unique<Eigen::IterScaling<EigenMatrix::RawMatrixType>>();
341  scal->computeRef(A.getRawMatrix());
342  b.getRawVector() = scal->LeftScaling().cwiseProduct(b.getRawVector());
343  }
344 #endif
345  auto const success = solver_->solve(A.getRawMatrix(), b.getRawVector(),
346  x.getRawVector(), option_);
347 #ifdef USE_EIGEN_UNSUPPORTED
348  if (scal)
349  {
350  x.getRawVector() = scal->RightScaling().cwiseProduct(x.getRawVector());
351  }
352 #endif
353 
354  INFO("------------------------------------------------------------------");
355 
356  return success;
357 }
358 
359 } // namespace MathLib
#define OGS_FATAL(...)
Definition: Error.h:26
void INFO(char const *fmt, Args const &... args)
Definition: Logging.h:32
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:42
std::optional< ConfigTree > getConfigSubtreeOptional(std::string const &root) const
Definition: ConfigTree.cpp:155
EigenMatrix::RawMatrixType Matrix
EigenVector::RawVectorType Vector
virtual bool solve(Matrix &A, Vector const &b, Vector &x, EigenOption &opt)=0
Solves the linear equation system for .
virtual ~EigenLinearSolverBase()=default
void setOption(const BaseLib::ConfigTree &option)
EigenLinearSolver(const std::string &solver_name, BaseLib::ConfigTree const *const option)
bool solve(EigenMatrix &A, EigenVector &b, EigenVector &x)
std::unique_ptr< EigenLinearSolverBase > solver_
Eigen::SparseMatrix< double, Eigen::RowMajor > RawMatrixType
Definition: EigenMatrix.h:31
RawMatrixType & getRawMatrix()
Definition: EigenMatrix.h:162
Global vector based on Eigen vector.
Definition: EigenVector.h:26
RawVectorType & getRawVector()
return a raw Eigen vector object
Definition: EigenVector.h:110
Eigen::VectorXd RawVectorType
Definition: EigenVector.h:28
Template class for Eigen direct linear solvers.
bool solve(Matrix &A, Vector const &b, Vector &x, EigenOption &opt) override
Solves the linear equation system for .
Template class for Eigen iterative linear solvers.
bool solve(Matrix &A, Vector const &b, Vector &x, EigenOption &opt) override
Solves the linear equation system for .
std::unique_ptr< EigenLinearSolverBase > createIterativeSolver()
Eigen::ConjugateGradient< Mat, Eigen::Lower, Precon > EigenCGSolver
void ignoreOtherLinearSolvers(const BaseLib::ConfigTree &config, const std::string &solver_name)
Option for Eigen sparse solver.
Definition: EigenOption.h:19
PreconType precon_type
Preconditioner type.
Definition: EigenOption.h:41
static std::string getSolverName(SolverType const solver_type)
return a linear solver name from the solver type
Definition: EigenOption.cpp:75
PreconType
Preconditioner type.
Definition: EigenOption.h:32
static PreconType getPreconType(const std::string &precon_name)
Definition: EigenOption.cpp:56
SolverType solver_type
Linear solver type.
Definition: EigenOption.h:39
double error_tolerance
Error tolerance.
Definition: EigenOption.h:45
SolverType
Solver type.
Definition: EigenOption.h:22
static SolverType getSolverType(const std::string &solver_name)
Definition: EigenOption.cpp:29
static std::string getPreconName(PreconType const precon_type)
return a preconditioner name from the preconditioner type
Definition: EigenOption.cpp:93
int max_iterations
Maximum iteration count.
Definition: EigenOption.h:43