OGS 6.2.0-97-g4a610c866
EigenLinearSolver.cpp
Go to the documentation of this file.
1 
10 #include "EigenLinearSolver.h"
11 
12 #include <logog/include/logog.hpp>
13 
14 #ifdef USE_MKL
15 #include <Eigen/PardisoSupport>
16 #endif
17 
18 #ifdef USE_EIGEN_UNSUPPORTED
19 #include <Eigen/Sparse>
20 #include <unsupported/Eigen/src/IterativeSolvers/GMRES.h>
21 #include <unsupported/Eigen/src/IterativeSolvers/Scaling.h>
22 #endif
23 
24 #include "BaseLib/ConfigTree.h"
25 #include "EigenVector.h"
26 #include "EigenMatrix.h"
27 #include "EigenTools.h"
28 
30 
31 namespace MathLib
32 {
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, EigenOption &opt) = 0;
45 };
46 
47 namespace details
48 {
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",
59  if (!A.isCompressed())
60  {
61  A.makeCompressed();
62  }
63 
64  _solver.compute(A);
65  if(_solver.info()!=Eigen::Success) {
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  ERR("Failed during Eigen linear solve");
73  return false;
74  }
75 
76  return true;
77  }
78 
79 private:
80  T_SOLVER _solver;
81 };
82 
84 template <class T_SOLVER>
86 {
87 public:
88  bool solve(Matrix& A, Vector const& b, Vector& x, EigenOption& opt) override
89  {
90  INFO("-> solve with %s (precon %s)",
93  _solver.setTolerance(opt.error_tolerance);
94  _solver.setMaxIterations(opt.max_iterations);
95 
96  if (!A.isCompressed())
97  {
98  A.makeCompressed();
99  }
100 
101  _solver.compute(A);
102  if(_solver.info()!=Eigen::Success) {
103  ERR("Failed during Eigen linear solver initialization");
104  return false;
105  }
106 
107  x = _solver.solveWithGuess(b, x);
108  INFO("\t iteration: %d/%ld", _solver.iterations(), opt.max_iterations);
109  INFO("\t residual: %e\n", _solver.error());
110 
111  if(_solver.info()!=Eigen::Success) {
112  ERR("Failed during Eigen linear solve");
113  return false;
114  }
115 
116  return true;
117  }
118 
119 private:
120  T_SOLVER _solver;
121 };
122 
123 template <template <typename, typename> class Solver, typename Precon>
124 std::unique_ptr<EigenLinearSolverBase> createIterativeSolver()
125 {
126  using Slv = EigenIterativeLinearSolver<
127  Solver<EigenMatrix::RawMatrixType, Precon>>;
128  return std::make_unique<Slv>();
129 }
130 
131 template <template <typename, typename> class Solver>
132 std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
133  EigenOption::PreconType precon_type)
134 {
135  switch (precon_type) {
137  return createIterativeSolver<Solver,
138  Eigen::IdentityPreconditioner>();
140  return createIterativeSolver<
141  Solver, Eigen::DiagonalPreconditioner<double>>();
143  // TODO for this preconditioner further options can be passed.
144  // see https://eigen.tuxfamily.org/dox/classEigen_1_1IncompleteLUT.html
145  return createIterativeSolver<
146  Solver, Eigen::IncompleteLUT<double>>();
147  default:
148  OGS_FATAL("Invalid Eigen preconditioner type.");
149  }
150 }
151 
152 template <typename Mat, typename Precon>
153 using EigenCGSolver = Eigen::ConjugateGradient<Mat, Eigen::Lower, Precon>;
154 
155 std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
156  EigenOption::SolverType solver_type, EigenOption::PreconType precon_type)
157 {
158  switch (solver_type) {
160  return createIterativeSolver<Eigen::BiCGSTAB>(precon_type);
161  }
163  return createIterativeSolver<EigenCGSolver>(precon_type);
164  }
166 #ifdef USE_EIGEN_UNSUPPORTED
167  return createIterativeSolver<Eigen::GMRES>(precon_type);
168 #else
169  OGS_FATAL(
170  "The code is not compiled with the Eigen unsupported modules. "
171  "Linear solver type GMRES is not available.");
172 #endif
173  }
174  default:
175  OGS_FATAL("Invalid Eigen iterative linear solver type. Aborting.");
176  }
177 }
178 
179 } // namespace details
180 
182  const std::string& /*solver_name*/,
183  const BaseLib::ConfigTree* const option)
184 {
186 
187  if (option)
188  {
189  setOption(*option);
190  }
191 
192  // TODO for my taste it is much too unobvious that the default solver type
193  // currently is SparseLU.
194  switch (_option.solver_type) {
196  using SolverType =
197  Eigen::SparseLU<Matrix, Eigen::COLAMDOrdering<int>>;
198  _solver = std::make_unique<
200  return;
201  }
205  _solver = details::createIterativeSolver(_option.solver_type,
206  _option.precon_type);
207  return;
209 #ifdef USE_MKL
210  using SolverType = Eigen::PardisoLU<EigenMatrix::RawMatrixType>;
212  return;
213 #else
214  OGS_FATAL(
215  "The code is not compiled with Intel MKL. Linear solver type "
216  "PardisoLU is not available.");
217 #endif
218  }
219  }
220 
221  OGS_FATAL("Invalid Eigen linear solver type. Aborting.");
222 }
223 
225 
227 {
228  ignoreOtherLinearSolvers(option, "eigen");
230  auto const ptSolver = option.getConfigSubtreeOptional("eigen");
231  if (!ptSolver)
232  {
233  return;
234  }
235 
236  if (auto solver_type =
238  ptSolver->getConfigParameterOptional<std::string>("solver_type")) {
239  _option.solver_type = _option.getSolverType(*solver_type);
240  }
241  if (auto precon_type =
243  ptSolver->getConfigParameterOptional<std::string>("precon_type")) {
244  _option.precon_type = _option.getPreconType(*precon_type);
245  }
246  if (auto error_tolerance =
248  ptSolver->getConfigParameterOptional<double>("error_tolerance")) {
249  _option.error_tolerance = *error_tolerance;
250  }
251  if (auto max_iteration_step =
253  ptSolver->getConfigParameterOptional<int>("max_iteration_step")) {
254  _option.max_iterations = *max_iteration_step;
255  }
256  if (auto scaling =
258  ptSolver->getConfigParameterOptional<bool>("scaling")) {
259 #ifdef USE_EIGEN_UNSUPPORTED
260  _option.scaling = *scaling;
261 #else
262  OGS_FATAL(
263  "The code is not compiled with the Eigen unsupported modules. "
264  "scaling is not available.");
265 #endif
266  }
267 }
268 
270 {
271  INFO("------------------------------------------------------------------");
272  INFO("*** Eigen solver computation");
273 
274 #ifdef USE_EIGEN_UNSUPPORTED
275  std::unique_ptr<Eigen::IterScaling<EigenMatrix::RawMatrixType>> scal;
276  if (_option.scaling)
277  {
278  INFO("-> scale");
279  scal =
280  std::make_unique<Eigen::IterScaling<EigenMatrix::RawMatrixType>>();
281  scal->computeRef(A.getRawMatrix());
282  b.getRawVector() = scal->LeftScaling().cwiseProduct(b.getRawVector());
283  }
284 #endif
285  auto const success = _solver->solve(A.getRawMatrix(), b.getRawVector(),
286  x.getRawVector(), _option);
287 #ifdef USE_EIGEN_UNSUPPORTED
288  if (scal)
289  {
290  x.getRawVector() = scal->RightScaling().cwiseProduct(x.getRawVector());
291  }
292 #endif
293 
294  INFO("------------------------------------------------------------------");
295 
296  return success;
297 }
298 
299 } // namespace MathLib
Eigen::ConjugateGradient< Mat, Eigen::Lower, Precon > EigenCGSolver
static std::string getPreconName(PreconType const precon_type)
return a preconditioner name from the preconditioner type
Definition: EigenOption.cpp:88
Template class for Eigen direct linear solvers.
PreconType
Preconditioner type.
Definition: EigenOption.h:31
Eigen::VectorXd RawVectorType
Definition: EigenVector.h:28
RawVectorType & getRawVector()
return a raw Eigen vector object
Definition: EigenVector.h:115
double error_tolerance
Error tolerance.
Definition: EigenOption.h:45
SolverType
Solver type.
Definition: EigenOption.h:21
EigenMatrix::RawMatrixType Matrix
virtual bool solve(Matrix &A, Vector const &b, Vector &x, EigenOption &opt)=0
Solves the linear equation system for .
int max_iterations
Maximum iteration count.
Definition: EigenOption.h:43
static std::string getSolverName(SolverType const solver_type)
return a linear solver name from the solver type
Definition: EigenOption.cpp:71
EigenLinearSolver(const std::string &solver_name, BaseLib::ConfigTree const *const option)
std::unique_ptr< EigenLinearSolverBase > createIterativeSolver(EigenOption::SolverType solver_type, EigenOption::PreconType precon_type)
Option for Eigen sparse solver.
Definition: EigenOption.h:18
PreconType precon_type
Preconditioner type.
Definition: EigenOption.h:41
Global vector based on Eigen vector.
Definition: EigenVector.h:25
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.
void ignoreOtherLinearSolvers(const BaseLib::ConfigTree &config, const std::string &solver_name)
void setOption(const BaseLib::ConfigTree &option)
bool solve(Matrix &A, Vector const &b, Vector &x, EigenOption &opt) override
Solves the linear equation system for .
EigenVector::RawVectorType Vector
RawMatrixType & getRawMatrix()
Definition: EigenMatrix.h:174
#define OGS_FATAL(fmt,...)
Definition: Error.h:63
boost::optional< ConfigTree > getConfigSubtreeOptional(std::string const &root) const
Definition: ConfigTree.cpp:160
SolverType solver_type
Linear solver type.
Definition: EigenOption.h:39
std::unique_ptr< EigenLinearSolverBase > createIterativeSolver()
virtual ~EigenLinearSolverBase()=default
bool solve(EigenMatrix &A, EigenVector &b, EigenVector &x)
Eigen::SparseMatrix< double, Eigen::RowMajor > RawMatrixType
Definition: EigenMatrix.h:32