OGS 6.1.0-1721-g6382411ad
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()) A.makeCompressed();
60 
61  _solver.compute(A);
62  if(_solver.info()!=Eigen::Success) {
63  ERR("Failed during Eigen linear solver initialization");
64  return false;
65  }
66 
67  x = _solver.solve(b);
68  if(_solver.info()!=Eigen::Success) {
69  ERR("Failed during Eigen linear solve");
70  return false;
71  }
72 
73  return true;
74  }
75 
76 private:
77  T_SOLVER _solver;
78 };
79 
81 template <class T_SOLVER>
83 {
84 public:
85  bool solve(Matrix& A, Vector const& b, Vector& x, EigenOption& opt) override
86  {
87  INFO("-> solve with %s (precon %s)",
90  _solver.setTolerance(opt.error_tolerance);
91  _solver.setMaxIterations(opt.max_iterations);
92 
93  if (!A.isCompressed())
94  A.makeCompressed();
95 
96  _solver.compute(A);
97  if(_solver.info()!=Eigen::Success) {
98  ERR("Failed during Eigen linear solver initialization");
99  return false;
100  }
101 
102  x = _solver.solveWithGuess(b, x);
103  INFO("\t iteration: %d/%ld", _solver.iterations(), opt.max_iterations);
104  INFO("\t residual: %e\n", _solver.error());
105 
106  if(_solver.info()!=Eigen::Success) {
107  ERR("Failed during Eigen linear solve");
108  return false;
109  }
110 
111  return true;
112  }
113 
114 private:
115  T_SOLVER _solver;
116 };
117 
118 template <template <typename, typename> class Solver, typename Precon>
119 std::unique_ptr<EigenLinearSolverBase> createIterativeSolver()
120 {
121  using Slv = EigenIterativeLinearSolver<
122  Solver<EigenMatrix::RawMatrixType, Precon>>;
123  return std::make_unique<Slv>();
124 }
125 
126 template <template <typename, typename> class Solver>
127 std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
128  EigenOption::PreconType precon_type)
129 {
130  switch (precon_type) {
132  return createIterativeSolver<Solver,
133  Eigen::IdentityPreconditioner>();
135  return createIterativeSolver<
136  Solver, Eigen::DiagonalPreconditioner<double>>();
138  // TODO for this preconditioner further options can be passed.
139  // see https://eigen.tuxfamily.org/dox/classEigen_1_1IncompleteLUT.html
140  return createIterativeSolver<
141  Solver, Eigen::IncompleteLUT<double>>();
142  default:
143  OGS_FATAL("Invalid Eigen preconditioner type.");
144  }
145 }
146 
147 template <typename Mat, typename Precon>
148 using EigenCGSolver = Eigen::ConjugateGradient<Mat, Eigen::Lower, Precon>;
149 
150 std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
151  EigenOption::SolverType solver_type, EigenOption::PreconType precon_type)
152 {
153  switch (solver_type) {
155  return createIterativeSolver<Eigen::BiCGSTAB>(precon_type);
156  }
158  return createIterativeSolver<EigenCGSolver>(precon_type);
159  }
161 #ifdef USE_EIGEN_UNSUPPORTED
162  return createIterativeSolver<Eigen::GMRES>(precon_type);
163 #else
164  OGS_FATAL(
165  "The code is not compiled with the Eigen unsupported modules. "
166  "Linear solver type GMRES is not available.");
167 #endif
168  }
169  default:
170  OGS_FATAL("Invalid Eigen iterative linear solver type. Aborting.");
171  }
172 }
173 
174 } // details
175 
177  const std::string& /*solver_name*/,
178  const BaseLib::ConfigTree* const option)
179 {
181 
182  if (option)
183  setOption(*option);
184 
185  // TODO for my taste it is much too unobvious that the default solver type
186  // currently is SparseLU.
187  switch (_option.solver_type) {
189  using SolverType =
190  Eigen::SparseLU<Matrix, Eigen::COLAMDOrdering<int>>;
191  _solver = std::make_unique<
193  return;
194  }
198  _solver = details::createIterativeSolver(_option.solver_type,
199  _option.precon_type);
200  return;
202 #ifdef USE_MKL
203  using SolverType = Eigen::PardisoLU<EigenMatrix::RawMatrixType>;
205  return;
206 #else
207  OGS_FATAL(
208  "The code is not compiled with Intel MKL. Linear solver type "
209  "PardisoLU is not available.");
210 #endif
211  }
212  }
213 
214  OGS_FATAL("Invalid Eigen linear solver type. Aborting.");
215 }
216 
218 
220 {
221  ignoreOtherLinearSolvers(option, "eigen");
223  auto const ptSolver = option.getConfigSubtreeOptional("eigen");
224  if (!ptSolver)
225  return;
226 
227  if (auto solver_type =
229  ptSolver->getConfigParameterOptional<std::string>("solver_type")) {
230  _option.solver_type = _option.getSolverType(*solver_type);
231  }
232  if (auto precon_type =
234  ptSolver->getConfigParameterOptional<std::string>("precon_type")) {
235  _option.precon_type = _option.getPreconType(*precon_type);
236  }
237  if (auto error_tolerance =
239  ptSolver->getConfigParameterOptional<double>("error_tolerance")) {
240  _option.error_tolerance = *error_tolerance;
241  }
242  if (auto max_iteration_step =
244  ptSolver->getConfigParameterOptional<int>("max_iteration_step")) {
245  _option.max_iterations = *max_iteration_step;
246  }
247  if (auto scaling =
249  ptSolver->getConfigParameterOptional<bool>("scaling")) {
250 #ifdef USE_EIGEN_UNSUPPORTED
251  _option.scaling = *scaling;
252 #else
253  OGS_FATAL(
254  "The code is not compiled with the Eigen unsupported modules. "
255  "scaling is not available.");
256 #endif
257  }
258 }
259 
261 {
262  INFO("------------------------------------------------------------------");
263  INFO("*** Eigen solver computation");
264 
265 #ifdef USE_EIGEN_UNSUPPORTED
266  std::unique_ptr<Eigen::IterScaling<EigenMatrix::RawMatrixType>> scal;
267  if (_option.scaling)
268  {
269  INFO("-> scale");
270  scal =
271  std::make_unique<Eigen::IterScaling<EigenMatrix::RawMatrixType>>();
272  scal->computeRef(A.getRawMatrix());
273  b.getRawVector() = scal->LeftScaling().cwiseProduct(b.getRawVector());
274  }
275 #endif
276  auto const success = _solver->solve(A.getRawMatrix(), b.getRawVector(),
277  x.getRawVector(), _option);
278 #ifdef USE_EIGEN_UNSUPPORTED
279  if (scal)
280  x.getRawVector() = scal->RightScaling().cwiseProduct(x.getRawVector());
281 #endif
282 
283  INFO("------------------------------------------------------------------");
284 
285  return success;
286 }
287 
288 } //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:72
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:55
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:166
#define OGS_FATAL(fmt,...)
Definition: Error.h:71
boost::optional< ConfigTree > getConfigSubtreeOptional(std::string const &root) const
Definition: ConfigTree.cpp:156
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