13#include <Eigen/Sparse>
19#include <Eigen/PardisoSupport>
22#ifdef USE_EIGEN_UNSUPPORTED
25#include <unsupported/Eigen/IterativeSolvers>
26#include <unsupported/Eigen/src/IterativeSolvers/Scaling.h>
46#ifdef USE_EIGEN_UNSUPPORTED
49 b = scaling_->LeftScaling().cwiseProduct(b);
53 auto const success =
solveImpl(b, x, opt);
55#ifdef USE_EIGEN_UNSUPPORTED
58 x = scaling_->RightScaling().cwiseProduct(x);
68#ifdef USE_EIGEN_UNSUPPORTED
72 scaling_ = std::make_unique<
73 Eigen::IterScaling<EigenMatrix::RawMatrixType>>();
74 scaling_->computeRef(A);
78 return computeImpl(A, opt, linear_solver_behaviour);
89#ifdef USE_EIGEN_UNSUPPORTED
90 std::unique_ptr<Eigen::IterScaling<EigenMatrix::RawMatrixType>> scaling_;
97template <
class T_SOLVER>
103 INFO(
"-> solve with Eigen direct linear solver {:s}",
107 if (
solver_.info() != Eigen::Success || x.hasNaN())
109 ERR(
"Failed during Eigen linear solve");
118 linear_solver_behaviour)
override
120 INFO(
"-> compute with Eigen direct linear solver {:s}",
122 if (!A.isCompressed())
128 if (
solver_.info() != Eigen::Success)
130 ERR(
"Failed during Eigen linear solver initialization");
141#ifdef USE_EIGEN_UNSUPPORTED
145template <
typename Solver>
146void setRestartImpl(Solver&,
int const)
148 DBUG(
"-> restart is not implemented for this linear solver.");
151template <
typename Matrix,
typename Precon>
152void setRestartImpl(Eigen::GMRES<Matrix, Precon>& solver,
int const restart)
154 solver.set_restart(restart);
155 INFO(
"-> set restart value: {:d}", solver.get_restart());
159template <
typename Solver>
160void setLImpl(Solver&,
int const)
162 DBUG(
"-> setL() is not implemented for this linear solver.");
165template <
typename Matrix,
typename Precon>
166void setLImpl(Eigen::BiCGSTABL<Matrix, Precon>& solver,
int const l)
171template <
typename Matrix,
typename Precon>
172void setLImpl(Eigen::IDRSTABL<Matrix, Precon>& solver,
int const l)
178template <
typename Solver>
179void setSImpl(Solver&,
int const)
181 DBUG(
"-> setS() is not implemented for this linear solver.");
184template <
typename Matrix,
typename Precon>
185void setSImpl(Eigen::IDRS<Matrix, Precon>& solver,
int const s)
190template <
typename Matrix,
typename Precon>
191void setSImpl(Eigen::IDRSTABL<Matrix, Precon>& solver,
int const s)
197template <
typename Solver>
198void setAngleImpl(Solver&,
double const)
200 DBUG(
"-> setAngle() is not implemented for this linear solver.");
203template <
typename Matrix,
typename Precon>
204void setAngleImpl(Eigen::IDRS<Matrix, Precon>& solver,
double const angle)
206 solver.setAngle(angle);
210template <
typename Solver>
211void setSmoothingImpl(Solver&,
bool const)
213 DBUG(
"-> setSmoothing() is not implemented for this linear solver.");
216template <
typename Matrix,
typename Precon>
217void setSmoothingImpl(Eigen::IDRS<Matrix, Precon>& solver,
bool const smoothing)
219 solver.setSmoothing(smoothing);
223template <
typename Solver>
224void setResidualUpdateImpl(Solver&,
bool const)
226 DBUG(
"-> setResidualUpdate() is not implemented for this linear solver.");
229template <
typename Matrix,
typename Precon>
230void setResidualUpdateImpl(Eigen::IDRS<Matrix, Precon>& solver,
231 bool const residual_update)
233 solver.setResidualUpdate(residual_update);
240template <
class T_SOLVER>
248 INFO(
"-> compute with Eigen iterative linear solver {:s} (precon {:s})",
254#ifdef USE_EIGEN_UNSUPPORTED
264 T_SOLVER>::setResidualUpdate(opt.residualupdate);
269 if (!m.isCompressed())
277 switch (linear_solver_behaviour)
295 "If NumLib::LinearSolverBehaviour::REUSE is set then "
296 "EigenLinearSolver::compute() should never be executed");
299 if (
solver_.info() != Eigen::Success)
301 ERR(
"Failed during Eigen linear solver initialization");
310 INFO(
"-> solve with Eigen iterative linear solver {:s} (precon {:s})",
314 x =
solver_.solveWithGuess(b, x);
315 INFO(
"\t iteration: {:d}/{:d}",
solver_.iterations(),
319 if (
solver_.info() != Eigen::Success)
321 ERR(
"Failed during Eigen linear solve");
331#ifdef USE_EIGEN_UNSUPPORTED
332 void setRestart(
int const restart) { setRestartImpl(
solver_, restart); }
333 void setL(
int const l) { setLImpl(
solver_, l); }
334 void setS(
int const s) { setSImpl(
solver_,
s); }
335 void setAngle(
double const angle) { setAngleImpl(
solver_, angle); }
336 void setSmoothing(
bool const smoothing)
338 setSmoothingImpl(
solver_, smoothing);
340 void setResidualUpdate(
bool const residual_update)
342 setResidualUpdateImpl(
solver_, residual_update);
347template <
template <
typename,
typename>
class Solver,
typename Precon>
352 return std::make_unique<Slv>();
355template <
template <
typename,
typename>
class Solver>
363 Eigen::IdentityPreconditioner>();
366 Solver, Eigen::DiagonalPreconditioner<double>>();
372 Eigen::IncompleteLUT<double>>();
374 OGS_FATAL(
"Invalid Eigen preconditioner type.");
378template <
typename Mat,
typename Precon>
388 return createIterativeSolver<Eigen::BiCGSTAB>(precon_type);
392#ifdef USE_EIGEN_UNSUPPORTED
393#if EIGEN_VERSION_AT_LEAST(3, 4, 90)
394 return createIterativeSolver<Eigen::BiCGSTABL>(precon_type);
396 OGS_FATAL(
"BiCGSTABL requires at least Eigen version 3.4.90");
400 "The code is not compiled with the Eigen unsupported modules. "
401 "Linear solver type BiCGSTABL is not available.");
406 return createIterativeSolver<EigenCGSolver>(precon_type);
410#ifdef USE_EIGEN_UNSUPPORTED
411 return createIterativeSolver<Eigen::GMRES>(precon_type);
414 "The code is not compiled with the Eigen unsupported modules. "
415 "Linear solver type GMRES is not available.");
420#ifdef USE_EIGEN_UNSUPPORTED
421#if EIGEN_VERSION_AT_LEAST(3, 4, 90)
422 return createIterativeSolver<Eigen::IDRS>(precon_type);
424 OGS_FATAL(
"IDRS requires at least Eigen version 3.4.90");
428 "The code is not compiled with the Eigen unsupported modules. "
429 "Linear solver type IDRS is not available.");
434#ifdef USE_EIGEN_UNSUPPORTED
435#if EIGEN_VERSION_AT_LEAST(3, 4, 90)
436 return createIterativeSolver<Eigen::IDRSTABL>(precon_type);
438 OGS_FATAL(
"IDRSTABL requires at least Eigen version 3.4.90");
442 "The code is not compiled with the Eigen unsupported modules. "
443 "Linear solver type IDRSTABL is not available.");
447 OGS_FATAL(
"Invalid Eigen iterative linear solver type. Aborting.");
464 Eigen::SparseLU<Matrix, Eigen::COLAMDOrdering<int>>;
481 using SolverType = Eigen::PardisoLU<EigenMatrix::RawMatrixType>;
486 "The code is not compiled with Intel MKL. Linear solver type "
487 "PardisoLU is not available.");
492 OGS_FATAL(
"Invalid Eigen linear solver type. Aborting.");
501 INFO(
"------------------------------------------------------------------");
502 INFO(
"*** Eigen solver compute()");
509 INFO(
"------------------------------------------------------------------");
510 INFO(
"*** Eigen solver solve()");
520 linear_solver_behaviour) &&
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
void DBUG(fmt::format_string< Args... > fmt, Args &&... args)
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
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(Matrix &A, EigenOption &opt, MathLib::LinearSolverBehaviour const linear_solver_behaviour)
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
RawMatrixType & getRawMatrix()
Global vector based on Eigen vector.
Eigen::VectorXd RawVectorType
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
Option for Eigen sparse solver.
PreconType precon_type
Preconditioner type.
static std::string getSolverName(SolverType const solver_type)
return a linear solver name from the solver type
PreconType
Preconditioner type.
SolverType solver_type
Linear solver type.
double error_tolerance
Error tolerance.
static std::string getPreconName(PreconType const precon_type)
return a preconditioner name from the preconditioner type
int max_iterations
Maximum iteration count.