OGS
EigenLinearSolver.cpp
Go to the documentation of this file.
1
11#include "EigenLinearSolver.h"
12
13#include <Eigen/Sparse>
14
15#include "BaseLib/Error.h"
16#include "BaseLib/Logging.h"
17
18#ifdef USE_MKL
19#include <Eigen/PardisoSupport>
20#endif
21
22#ifdef USE_EIGEN_UNSUPPORTED
23// TODO(naumov); Don't include header files directly from Eigen's source dir
24// clang-format off
25#include <unsupported/Eigen/IterativeSolvers>
26#include <unsupported/Eigen/src/IterativeSolvers/Scaling.h>
27// clang-format on
28#endif
29
30#include "EigenMatrix.h"
31#include "EigenVector.h"
32
33namespace MathLib
34{
35// TODO change to LinearSolver
37{
38public:
41
42 virtual ~EigenLinearSolverBase() = default;
43
44 bool solve(Vector& b, Vector& x, EigenOption& opt)
45 {
46#ifdef USE_EIGEN_UNSUPPORTED
47 if (scaling_)
48 {
49 b = scaling_->LeftScaling().cwiseProduct(b);
50 }
51#endif
52
53 auto const success = solveImpl(b, x, opt);
54
55#ifdef USE_EIGEN_UNSUPPORTED
56 if (scaling_)
57 {
58 x = scaling_->RightScaling().cwiseProduct(x);
59 }
60#endif
61
62 return success;
63 }
64
65 bool compute(Matrix& A, EigenOption& opt,
66 MathLib::LinearSolverBehaviour linear_solver_behaviour)
67 {
68 if (linear_solver_behaviour == MathLib::LinearSolverBehaviour::REUSE)
69 {
70 // Checking if this is the first compute() call should work both for
71 // direct solvers (that factorize the matrix A) and for iterative
72 // solvers (that store a reference to A in the preconditioner).
74 {
75 // There is nothing there, yet, to be reused. Hence, we compute
76 // and store the result.
77 linear_solver_behaviour =
79 }
80 else
81 {
82 return true;
83 }
84 }
85
86 // TODO (CL) That might not work under all circumstances, e.g., if the
87 // linear solver fails subsequently. Time will tell.
89
90#ifdef USE_EIGEN_UNSUPPORTED
91 if (opt.scaling)
92 {
93 INFO("-> scale");
94 scaling_ = std::make_unique<
95 Eigen::IterScaling<EigenMatrix::RawMatrixType>>();
96 scaling_->computeRef(A);
97 }
98#endif
99
100 return computeImpl(A, opt, linear_solver_behaviour);
101 }
102
103protected:
104 virtual bool solveImpl(Vector const& b, Vector& x, EigenOption& opt) = 0;
105
106 virtual bool computeImpl(
107 Matrix& A, EigenOption& opt,
108 MathLib::LinearSolverBehaviour const linear_solver_behaviour) = 0;
109
110private:
111#ifdef USE_EIGEN_UNSUPPORTED
112 std::unique_ptr<Eigen::IterScaling<EigenMatrix::RawMatrixType>> scaling_;
113#endif
114
116};
117
118namespace details
119{
121template <class T_SOLVER>
123{
124public:
125 bool solveImpl(Vector const& b, Vector& x, EigenOption& opt) override
126 {
127 INFO("-> solve with Eigen direct linear solver {:s}",
129
130 x = solver_.solve(b);
131 if (solver_.info() != Eigen::Success || x.hasNaN())
132 {
133 ERR("Failed during Eigen linear solve");
134 return false;
135 }
136
137 return true;
138 }
139
141 [[maybe_unused]] MathLib::LinearSolverBehaviour const
142 linear_solver_behaviour) override
143 {
144 INFO("-> compute with Eigen direct linear solver {:s}",
146 if (!A.isCompressed())
147 {
148 A.makeCompressed();
149 }
150
151 solver_.compute(A);
152 if (solver_.info() != Eigen::Success)
153 {
154 ERR("Failed during Eigen linear solver initialization");
155 return false;
156 }
157
158 return true;
159 }
160
161private:
162 T_SOLVER solver_;
163};
164
165#ifdef USE_EIGEN_UNSUPPORTED
166// implementations for some iterative linear solver methods --------------------
167
168// restart
169template <typename Solver>
170void setRestartImpl(Solver&, int const)
171{
172 DBUG("-> restart is not implemented for this linear solver.");
173}
174
175template <typename Matrix, typename Precon>
176void setRestartImpl(Eigen::GMRES<Matrix, Precon>& solver, int const restart)
177{
178 solver.set_restart(restart);
179 INFO("-> set restart value: {:d}", solver.get_restart());
180}
181
182// L
183template <typename Solver>
184void setLImpl(Solver&, int const)
185{
186 DBUG("-> setL() is not implemented for this linear solver.");
187}
188
189template <typename Matrix, typename Precon>
190void setLImpl(Eigen::BiCGSTABL<Matrix, Precon>& solver, int const l)
191{
192 solver.setL(l);
193}
194
195template <typename Matrix, typename Precon>
196void setLImpl(Eigen::IDRSTABL<Matrix, Precon>& solver, int const l)
197{
198 solver.setL(l);
199}
200
201// S
202template <typename Solver>
203void setSImpl(Solver&, int const)
204{
205 DBUG("-> setS() is not implemented for this linear solver.");
206}
207
208template <typename Matrix, typename Precon>
209void setSImpl(Eigen::IDRS<Matrix, Precon>& solver, int const s)
210{
211 solver.setS(s);
212}
213
214template <typename Matrix, typename Precon>
215void setSImpl(Eigen::IDRSTABL<Matrix, Precon>& solver, int const s)
216{
217 solver.setS(s);
218}
219
220// angle
221template <typename Solver>
222void setAngleImpl(Solver&, double const)
223{
224 DBUG("-> setAngle() is not implemented for this linear solver.");
225}
226
227template <typename Matrix, typename Precon>
228void setAngleImpl(Eigen::IDRS<Matrix, Precon>& solver, double const angle)
229{
230 solver.setAngle(angle);
231}
232
233// smoothing
234template <typename Solver>
235void setSmoothingImpl(Solver&, bool const)
236{
237 DBUG("-> setSmoothing() is not implemented for this linear solver.");
238}
239
240template <typename Matrix, typename Precon>
241void setSmoothingImpl(Eigen::IDRS<Matrix, Precon>& solver, bool const smoothing)
242{
243 solver.setSmoothing(smoothing);
244}
245
246// residual update
247template <typename Solver>
248void setResidualUpdateImpl(Solver&, bool const)
249{
250 DBUG("-> setResidualUpdate() is not implemented for this linear solver.");
251}
252
253template <typename Matrix, typename Precon>
254void setResidualUpdateImpl(Eigen::IDRS<Matrix, Precon>& solver,
255 bool const residual_update)
256{
257 solver.setResidualUpdate(residual_update);
258}
259
260// -----------------------------------------------------------------------------
261#endif
262
264template <class T_SOLVER>
266{
267public:
269 Matrix& A, EigenOption& opt,
270 MathLib::LinearSolverBehaviour const linear_solver_behaviour) override
271 {
272 INFO("-> compute with Eigen iterative linear solver {:s} (precon {:s})",
275 solver_.setTolerance(opt.error_tolerance);
276 solver_.setMaxIterations(opt.max_iterations);
277
278#ifdef USE_EIGEN_UNSUPPORTED
280 opt.restart);
284 opt.smoothing);
286 opt.angle);
288 T_SOLVER>::setResidualUpdate(opt.residualupdate);
289#endif
290
291 auto compute = [this](Matrix& m)
292 {
293 if (!m.isCompressed())
294 {
295 m.makeCompressed();
296 }
297
298 solver_.compute(m);
299 };
300
301 switch (linear_solver_behaviour)
302 {
304 {
305 // matrix must be copied, because Eigen's linear solver stores a
306 // reference to it cf.
307 // https://libeigen.gitlab.io/docs/classEigen_1_1IterativeSolverBase.html#a7dfa55c55e82d697bde227696a630914
308 A_ = A;
309 compute(A_);
310 break;
311 }
313 {
314 compute(A);
315 break;
316 }
318 OGS_FATAL(
319 "If NumLib::LinearSolverBehaviour::REUSE is set then "
320 "EigenLinearSolver::compute() should never be executed");
321 };
322
323 if (solver_.info() != Eigen::Success)
324 {
325 ERR("Failed during Eigen linear solver initialization");
326 return false;
327 }
328
329 return true;
330 }
331
332 bool solveImpl(Vector const& b, Vector& x, EigenOption& opt) override
333 {
334 INFO("-> solve with Eigen iterative linear solver {:s} (precon {:s})",
337
338 x = solver_.solveWithGuess(b, x);
339 INFO("\t iteration: {:d}/{:d}", solver_.iterations(),
340 opt.max_iterations);
341 INFO("\t residual: {:e}\n", solver_.error());
342
343 if (solver_.info() != Eigen::Success)
344 {
345 ERR("Failed during Eigen linear solve");
346 return false;
347 }
348
349 return true;
350 }
351
352private:
353 T_SOLVER solver_;
355#ifdef USE_EIGEN_UNSUPPORTED
356 void setRestart(int const restart) { setRestartImpl(solver_, restart); }
357 void setL(int const l) { setLImpl(solver_, l); }
358 void setS(int const s) { setSImpl(solver_, s); }
359 void setAngle(double const angle) { setAngleImpl(solver_, angle); }
360 void setSmoothing(bool const smoothing)
361 {
362 setSmoothingImpl(solver_, smoothing);
363 }
364 void setResidualUpdate(bool const residual_update)
365 {
366 setResidualUpdateImpl(solver_, residual_update);
367 }
368#endif
369};
370
371template <template <typename, typename> class Solver, typename Precon>
372std::unique_ptr<EigenLinearSolverBase> createIterativeSolver()
373{
374 using Slv =
376 return std::make_unique<Slv>();
377}
378
379template <template <typename, typename> class Solver>
380std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
381 EigenOption::PreconType precon_type)
382{
383 switch (precon_type)
384 {
386 return createIterativeSolver<Solver,
387 Eigen::IdentityPreconditioner>();
390 Solver, Eigen::DiagonalPreconditioner<double>>();
392 // TODO for this preconditioner further options can be passed.
393 // see
394 // https://libeigen.gitlab.io/docs/classEigen_1_1IncompleteLUT.html
395 return createIterativeSolver<Solver,
396 Eigen::IncompleteLUT<double>>();
397 default:
398 OGS_FATAL("Invalid Eigen preconditioner type.");
399 }
400}
401
402template <typename Mat, typename Precon>
403using EigenCGSolver = Eigen::ConjugateGradient<Mat, Eigen::Lower, Precon>;
404
405std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
406 EigenOption::SolverType solver_type, EigenOption::PreconType precon_type)
407{
408 switch (solver_type)
409 {
411 {
412 return createIterativeSolver<Eigen::BiCGSTAB>(precon_type);
413 }
415 {
416#ifdef USE_EIGEN_UNSUPPORTED
417#if EIGEN_VERSION_AT_LEAST(3, 4, 90)
418 return createIterativeSolver<Eigen::BiCGSTABL>(precon_type);
419#else
420 OGS_FATAL("BiCGSTABL requires at least Eigen version 3.4.90");
421#endif
422#else
423 OGS_FATAL(
424 "The code is not compiled with the Eigen unsupported modules. "
425 "Linear solver type BiCGSTABL is not available.");
426#endif
427 }
429 {
430 return createIterativeSolver<EigenCGSolver>(precon_type);
431 }
433 {
434#ifdef USE_EIGEN_UNSUPPORTED
435 return createIterativeSolver<Eigen::GMRES>(precon_type);
436#else
437 OGS_FATAL(
438 "The code is not compiled with the Eigen unsupported modules. "
439 "Linear solver type GMRES is not available.");
440#endif
441 }
443 {
444#ifdef USE_EIGEN_UNSUPPORTED
445#if EIGEN_VERSION_AT_LEAST(3, 4, 90)
446 return createIterativeSolver<Eigen::IDRS>(precon_type);
447#else
448 OGS_FATAL("IDRS requires at least Eigen version 3.4.90");
449#endif
450#else
451 OGS_FATAL(
452 "The code is not compiled with the Eigen unsupported modules. "
453 "Linear solver type IDRS is not available.");
454#endif
455 }
457 {
458#ifdef USE_EIGEN_UNSUPPORTED
459#if EIGEN_VERSION_AT_LEAST(3, 4, 90)
460 return createIterativeSolver<Eigen::IDRSTABL>(precon_type);
461#else
462 OGS_FATAL("IDRSTABL requires at least Eigen version 3.4.90");
463#endif
464#else
465 OGS_FATAL(
466 "The code is not compiled with the Eigen unsupported modules. "
467 "Linear solver type IDRSTABL is not available.");
468#endif
469 }
470 default:
471 OGS_FATAL("Invalid Eigen iterative linear solver type. Aborting.");
472 }
473}
474
475} // namespace details
476
477EigenLinearSolver::EigenLinearSolver(std::string const& /*solver_name*/,
478 EigenOption const& option)
479 : option_(option)
480{
481 using Matrix = EigenMatrix::RawMatrixType;
482
483 switch (option_.solver_type)
484 {
486 {
487 using SolverType =
488 Eigen::SparseLU<Matrix, Eigen::COLAMDOrdering<int>>;
489 solver_ = std::make_unique<
491 return;
492 }
501 return;
503 {
504#ifdef USE_MKL
505 using SolverType = Eigen::PardisoLU<EigenMatrix::RawMatrixType>;
507 return;
508#else
509 OGS_FATAL(
510 "The code is not compiled with Intel MKL. Linear solver type "
511 "PardisoLU is not available.");
512#endif
513 }
514 }
515
516 OGS_FATAL("Invalid Eigen linear solver type. Aborting.");
517}
518
520
522 EigenMatrix& A,
523 MathLib::LinearSolverBehaviour const linear_solver_behaviour)
524{
525 INFO("------------------------------------------------------------------");
526 INFO("*** Eigen solver compute()");
527
528 return solver_->compute(A.getRawMatrix(), option_, linear_solver_behaviour);
529}
530
532{
533 INFO("------------------------------------------------------------------");
534 INFO("*** Eigen solver solve()");
535
536 return solver_->solve(b.getRawVector(), x.getRawVector(), option_);
537}
538
541 MathLib::LinearSolverBehaviour const linear_solver_behaviour)
542{
543 return solver_->compute(A.getRawMatrix(), option_,
544 linear_solver_behaviour) &&
545 solver_->solve(b.getRawVector(), x.getRawVector(), option_);
546}
547
548} // namespace MathLib
#define OGS_FATAL(...)
Definition Error.h:26
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:30
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
bool compute(Matrix &A, EigenOption &opt, MathLib::LinearSolverBehaviour linear_solver_behaviour)
EigenMatrix::RawMatrixType Matrix
EigenVector::RawVectorType Vector
virtual bool computeImpl(Matrix &A, EigenOption &opt, MathLib::LinearSolverBehaviour const linear_solver_behaviour)=0
virtual bool solveImpl(Vector const &b, Vector &x, EigenOption &opt)=0
virtual ~EigenLinearSolverBase()=default
bool solve(Vector &b, Vector &x, EigenOption &opt)
bool compute(EigenMatrix &A, MathLib::LinearSolverBehaviour const linear_solver_behaviour)
bool solve(EigenVector &b, EigenVector &x)
EigenLinearSolver(std::string const &solver_name, EigenOption const &option)
std::unique_ptr< EigenLinearSolverBase > solver_
Eigen::SparseMatrix< double, Eigen::RowMajor > RawMatrixType
Definition EigenMatrix.h:31
RawMatrixType & getRawMatrix()
Global vector based on Eigen vector.
Definition EigenVector.h:25
Eigen::VectorXd RawVectorType
Definition EigenVector.h:27
RawVectorType & getRawVector()
return a raw Eigen vector object
Template class for Eigen direct linear solvers.
bool solveImpl(Vector const &b, Vector &x, EigenOption &opt) override
bool computeImpl(Matrix &A, EigenOption &opt, MathLib::LinearSolverBehaviour const linear_solver_behaviour) override
Template class for Eigen iterative linear solvers.
bool solveImpl(Vector const &b, Vector &x, EigenOption &opt) override
bool computeImpl(Matrix &A, EigenOption &opt, MathLib::LinearSolverBehaviour const linear_solver_behaviour) override
std::unique_ptr< EigenLinearSolverBase > createIterativeSolver()
Eigen::ConjugateGradient< Mat, Eigen::Lower, Precon > EigenCGSolver
static const double s
Option for Eigen sparse solver.
Definition EigenOption.h:19
PreconType precon_type
Preconditioner type.
Definition EigenOption.h:44
static std::string getSolverName(SolverType const solver_type)
return a linear solver name from the solver type
PreconType
Preconditioner type.
Definition EigenOption.h:35
SolverType solver_type
Linear solver type.
Definition EigenOption.h:42
double error_tolerance
Error tolerance.
Definition EigenOption.h:48
SolverType
Solver type.
Definition EigenOption.h:22
static std::string getPreconName(PreconType const precon_type)
return a preconditioner name from the preconditioner type
int max_iterations
Maximum iteration count.
Definition EigenOption.h:46