From e777bc01112240f0ceb352ae1d9f9f6903a5797e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B6rg=20Buchwald?= <joerg.buchwald@ufz.de>
Date: Fri, 7 Jan 2022 20:33:37 +0100
Subject: [PATCH] add IDRS, BiCGSTABL and IDRSTABL iterative solvers
---
MathLib/LinAlg/Eigen/EigenLinearSolver.cpp | 221 +++++++++++++++++-
MathLib/LinAlg/Eigen/EigenLinearSolver.h | 5 +
MathLib/LinAlg/Eigen/EigenOption.cpp | 23 ++
MathLib/LinAlg/Eigen/EigenOption.h | 8 +
.../LinAlg/Eigen/LinearSolverOptionsParser.h | 55 +++++
5 files changed, 311 insertions(+), 1 deletion(-)
diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
index 6154ee1393c..9e316e1e4ef 100644
--- a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
+++ b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
@@ -19,7 +19,7 @@
#endif
#ifdef USE_EIGEN_UNSUPPORTED
-#include <unsupported/Eigen/src/IterativeSolvers/GMRES.h>
+#include <unsupported/Eigen/IterativeSolvers>
#include <unsupported/Eigen/src/IterativeSolvers/Scaling.h>
#endif
@@ -94,6 +94,16 @@ public:
solver_.setMaxIterations(opt.max_iterations);
MathLib::details::EigenIterativeLinearSolver<T_SOLVER>::setRestart(
opt.restart);
+ MathLib::details::EigenIterativeLinearSolver<T_SOLVER>::setL(
+ opt.l);
+ MathLib::details::EigenIterativeLinearSolver<T_SOLVER>::setS(
+ opt.s);
+ MathLib::details::EigenIterativeLinearSolver<T_SOLVER>::setSmoothing(
+ opt.smoothing);
+ MathLib::details::EigenIterativeLinearSolver<T_SOLVER>::setAngle(
+ opt.angle);
+ MathLib::details::EigenIterativeLinearSolver<T_SOLVER>::setResidualUpdate(
+ opt.residualupdate);
if (!A.isCompressed())
{
@@ -124,6 +134,11 @@ public:
private:
T_SOLVER solver_;
void setRestart(int const /*restart*/) {}
+ void setL(int const /*l*/) {}
+ void setS(int const /*s*/) {}
+ void setAngle(double const /*angle*/) {}
+ void setSmoothing(bool const /*smoothing*/) {}
+ void setResidualUpdate(bool const /*residual update*/) {}
};
/// Specialization for (all) three preconditioners separately
@@ -154,6 +169,177 @@ void EigenIterativeLinearSolver<
INFO("-> set restart value: {:d}", solver_.get_restart());
}
+/// BiCGSTABL
+template <>
+void EigenIterativeLinearSolver<
+ Eigen::BiCGSTABL<EigenMatrix::RawMatrixType,
+ Eigen::IdentityPreconditioner>>::setL(int const l)
+{
+ solver_.setL(l);
+}
+
+template <>
+void EigenIterativeLinearSolver<Eigen::BiCGSTABL<
+ EigenMatrix::RawMatrixType,
+ Eigen::DiagonalPreconditioner<double>>>::setL(int const l)
+{
+ solver_.setL(l);
+}
+
+template <>
+void EigenIterativeLinearSolver<
+ Eigen::BiCGSTABL<EigenMatrix::RawMatrixType,
+ Eigen::IncompleteLUT<double>>>::setL(int const l)
+{
+ solver_.setL(l);
+}
+
+/// IDRS
+template <>
+void EigenIterativeLinearSolver<
+ Eigen::IDRS<EigenMatrix::RawMatrixType,
+ Eigen::IdentityPreconditioner>>::setS(int const s)
+{
+ solver_.setS(s);
+}
+
+template <>
+void EigenIterativeLinearSolver<Eigen::IDRS<
+ EigenMatrix::RawMatrixType,
+ Eigen::DiagonalPreconditioner<double>>>::setS(int const s)
+{
+ solver_.setS(s);
+}
+
+template <>
+void EigenIterativeLinearSolver<
+ Eigen::IDRS<EigenMatrix::RawMatrixType,
+ Eigen::IncompleteLUT<double>>>::setS(int const s)
+{
+ solver_.setS(s);
+}
+
+template <>
+void EigenIterativeLinearSolver<
+ Eigen::IDRS<EigenMatrix::RawMatrixType,
+ Eigen::IdentityPreconditioner>>::setAngle(double const angle)
+{
+ solver_.setAngle(angle);
+}
+
+template <>
+void EigenIterativeLinearSolver<Eigen::IDRS<
+ EigenMatrix::RawMatrixType,
+ Eigen::DiagonalPreconditioner<double>>>::setAngle(double const angle)
+{
+ solver_.setAngle(angle);
+}
+
+template <>
+void EigenIterativeLinearSolver<
+ Eigen::IDRS<EigenMatrix::RawMatrixType,
+ Eigen::IncompleteLUT<double>>>::setAngle(double const angle)
+{
+ solver_.setAngle(angle);
+}
+
+template <>
+void EigenIterativeLinearSolver<
+ Eigen::IDRS<EigenMatrix::RawMatrixType,
+ Eigen::IdentityPreconditioner>>::setSmoothing(bool const smoothing)
+{
+ solver_.setSmoothing(smoothing);
+}
+
+template <>
+void EigenIterativeLinearSolver<Eigen::IDRS<
+ EigenMatrix::RawMatrixType,
+ Eigen::DiagonalPreconditioner<double>>>::setSmoothing(bool const smoothing)
+{
+ solver_.setSmoothing(smoothing);
+}
+
+template <>
+void EigenIterativeLinearSolver<
+ Eigen::IDRS<EigenMatrix::RawMatrixType,
+ Eigen::IncompleteLUT<double>>>::setSmoothing(bool const smoothing)
+{
+ solver_.setSmoothing(smoothing);
+}
+
+template <>
+void EigenIterativeLinearSolver<
+ Eigen::IDRS<EigenMatrix::RawMatrixType,
+ Eigen::IdentityPreconditioner>>::setResidualUpdate(bool const residualupdate)
+{
+ solver_.setResidualUpdate(residualupdate);
+}
+
+template <>
+void EigenIterativeLinearSolver<Eigen::IDRS<
+ EigenMatrix::RawMatrixType,
+ Eigen::DiagonalPreconditioner<double>>>::setResidualUpdate(bool const residualupdate)
+{
+ solver_.setResidualUpdate(residualupdate);
+}
+
+template <>
+void EigenIterativeLinearSolver<
+ Eigen::IDRS<EigenMatrix::RawMatrixType,
+ Eigen::IncompleteLUT<double>>>::setResidualUpdate(bool const residualupdate)
+{
+ solver_.setResidualUpdate(residualupdate);
+}
+
+/// IDRSTABL
+template <>
+void EigenIterativeLinearSolver<
+ Eigen::IDRSTABL<EigenMatrix::RawMatrixType,
+ Eigen::IdentityPreconditioner>>::setS(int const s)
+{
+ solver_.setS(s);
+}
+
+template <>
+void EigenIterativeLinearSolver<Eigen::IDRSTABL<
+ EigenMatrix::RawMatrixType,
+ Eigen::DiagonalPreconditioner<double>>>::setS(int const s)
+{
+ solver_.setS(s);
+}
+
+template <>
+void EigenIterativeLinearSolver<
+ Eigen::IDRSTABL<EigenMatrix::RawMatrixType,
+ Eigen::IncompleteLUT<double>>>::setS(int const s)
+{
+ solver_.setS(s);
+}
+
+template <>
+void EigenIterativeLinearSolver<
+ Eigen::IDRSTABL<EigenMatrix::RawMatrixType,
+ Eigen::IdentityPreconditioner>>::setL(int const l)
+{
+ solver_.setL(l);
+}
+
+template <>
+void EigenIterativeLinearSolver<Eigen::IDRSTABL<
+ EigenMatrix::RawMatrixType,
+ Eigen::DiagonalPreconditioner<double>>>::setL(int const l)
+{
+ solver_.setL(l);
+}
+
+template <>
+void EigenIterativeLinearSolver<
+ Eigen::IDRSTABL<EigenMatrix::RawMatrixType,
+ Eigen::IncompleteLUT<double>>>::setL(int const l)
+{
+ solver_.setL(l);
+}
+
template <template <typename, typename> class Solver, typename Precon>
std::unique_ptr<EigenLinearSolverBase> createIterativeSolver()
{
@@ -197,6 +383,16 @@ std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
{
return createIterativeSolver<Eigen::BiCGSTAB>(precon_type);
}
+ case EigenOption::SolverType::BiCGSTABL:
+ {
+#ifdef USE_EIGEN_UNSUPPORTED
+ return createIterativeSolver<Eigen::BiCGSTABL>(precon_type);
+#else
+ OGS_FATAL(
+ "The code is not compiled with the Eigen unsupported modules. "
+ "Linear solver type BiCGSTABL is not available.");
+#endif
+ }
case EigenOption::SolverType::CG:
{
return createIterativeSolver<EigenCGSolver>(precon_type);
@@ -209,6 +405,26 @@ std::unique_ptr<EigenLinearSolverBase> createIterativeSolver(
OGS_FATAL(
"The code is not compiled with the Eigen unsupported modules. "
"Linear solver type GMRES is not available.");
+#endif
+ }
+ case EigenOption::SolverType::IDRS:
+ {
+#ifdef USE_EIGEN_UNSUPPORTED
+ return createIterativeSolver<Eigen::IDRS>(precon_type);
+#else
+ OGS_FATAL(
+ "The code is not compiled with the Eigen unsupported modules. "
+ "Linear solver type IDRS is not available.");
+#endif
+ }
+ case EigenOption::SolverType::IDRSTABL:
+ {
+#ifdef USE_EIGEN_UNSUPPORTED
+ return createIterativeSolver<Eigen::IDRSTABL>(precon_type);
+#else
+ OGS_FATAL(
+ "The code is not compiled with the Eigen unsupported modules. "
+ "Linear solver type IDRSTABL is not available.");
#endif
}
default:
@@ -235,8 +451,11 @@ EigenLinearSolver::EigenLinearSolver(std::string const& /*solver_name*/,
return;
}
case EigenOption::SolverType::BiCGSTAB:
+ case EigenOption::SolverType::BiCGSTABL:
case EigenOption::SolverType::CG:
case EigenOption::SolverType::GMRES:
+ case EigenOption::SolverType::IDRS:
+ case EigenOption::SolverType::IDRSTABL:
solver_ = details::createIterativeSolver(option_.solver_type,
option_.precon_type);
return;
diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.h b/MathLib/LinAlg/Eigen/EigenLinearSolver.h
index 7585f9528a1..772818d9216 100644
--- a/MathLib/LinAlg/Eigen/EigenLinearSolver.h
+++ b/MathLib/LinAlg/Eigen/EigenLinearSolver.h
@@ -52,6 +52,11 @@ protected:
EigenOption option_;
std::unique_ptr<EigenLinearSolverBase> solver_;
void setRestart();
+ void setL();
+ void setS();
+ void setSmoothing();
+ void setAngle();
+ void setResidualUpdate();
};
} // namespace MathLib
diff --git a/MathLib/LinAlg/Eigen/EigenOption.cpp b/MathLib/LinAlg/Eigen/EigenOption.cpp
index b0587329f1d..15f89cd55bc 100644
--- a/MathLib/LinAlg/Eigen/EigenOption.cpp
+++ b/MathLib/LinAlg/Eigen/EigenOption.cpp
@@ -23,6 +23,11 @@ EigenOption::EigenOption()
#ifdef USE_EIGEN_UNSUPPORTED
scaling = false;
restart = 30;
+ l = 2;
+ s = 4;
+ angle = 0.7;
+ smoothing = false;
+ residualupdate = false;
#endif
}
@@ -37,6 +42,18 @@ EigenOption::SolverType EigenOption::getSolverType(
{
return SolverType::BiCGSTAB;
}
+ if (solver_name == "BiCGSTABL")
+ {
+ return SolverType::BiCGSTABL;
+ }
+ if (solver_name == "IDRS")
+ {
+ return SolverType::IDRS;
+ }
+ if (solver_name == "IDRSTABL")
+ {
+ return SolverType::IDRSTABL;
+ }
if (solver_name == "SparseLU")
{
return SolverType::SparseLU;
@@ -80,6 +97,12 @@ std::string EigenOption::getSolverName(SolverType const solver_type)
return "CG";
case SolverType::BiCGSTAB:
return "BiCGSTAB";
+ case SolverType::BiCGSTABL:
+ return "BiCGSTABL";
+ case SolverType::IDRS:
+ return "IDRS";
+ case SolverType::IDRSTABL:
+ return "IDRSTABL";
case SolverType::SparseLU:
return "SparseLU";
case SolverType::PardisoLU:
diff --git a/MathLib/LinAlg/Eigen/EigenOption.h b/MathLib/LinAlg/Eigen/EigenOption.h
index 24ced185fad..61a88797cf9 100644
--- a/MathLib/LinAlg/Eigen/EigenOption.h
+++ b/MathLib/LinAlg/Eigen/EigenOption.h
@@ -22,6 +22,9 @@ struct EigenOption final
{
CG,
BiCGSTAB,
+ BiCGSTABL,
+ IDRS,
+ IDRSTABL,
SparseLU,
PardisoLU,
GMRES
@@ -48,6 +51,11 @@ struct EigenOption final
bool scaling;
/// Restart value for the GMRES solver
int restart;
+ int l;
+ int s;
+ double angle;
+ bool smoothing;
+ bool residualupdate;
#endif
/// Constructor
diff --git a/MathLib/LinAlg/Eigen/LinearSolverOptionsParser.h b/MathLib/LinAlg/Eigen/LinearSolverOptionsParser.h
index 775ca11ffd6..87c2ed0df17 100644
--- a/MathLib/LinAlg/Eigen/LinearSolverOptionsParser.h
+++ b/MathLib/LinAlg/Eigen/LinearSolverOptionsParser.h
@@ -101,6 +101,61 @@ struct LinearSolverOptionsParser<EigenLinearSolver> final
OGS_FATAL(
"The code is not compiled with the Eigen unsupported modules. "
"GMRES/GMRES option restart is not available.");
+#endif
+ }
+ if (auto l =
+ //! \ogs_file_param{prj__linear_solvers__linear_solver__eigen__l}
+ config->getConfigParameterOptional<int>("l"))
+ {
+#ifdef USE_EIGEN_UNSUPPORTED
+ options.l = *l;
+#else
+ OGS_FATAL(
+ "The code is not compiled with the Eigen unsupported modules.");
+#endif
+ }
+ if (auto s =
+ //! \ogs_file_param{prj__linear_solvers__linear_solver__eigen__s}
+ config->getConfigParameterOptional<int>("s"))
+ {
+#ifdef USE_EIGEN_UNSUPPORTED
+ options.s = *s;
+#else
+ OGS_FATAL(
+ "The code is not compiled with the Eigen unsupported modules.");
+#endif
+ }
+ if (auto smoothing =
+ //! \ogs_file_param{prj__linear_solvers__linear_solver__eigen__smoothing}
+ config->getConfigParameterOptional<int>("smoothing"))
+ {
+#ifdef USE_EIGEN_UNSUPPORTED
+ options.smoothing = *smoothing;
+#else
+ OGS_FATAL(
+ "The code is not compiled with the Eigen unsupported modules.");
+#endif
+ }
+ if (auto angle =
+ //! \ogs_file_param{prj__linear_solvers__linear_solver__eigen__angle}
+ config->getConfigParameterOptional<int>("angle"))
+ {
+#ifdef USE_EIGEN_UNSUPPORTED
+ options.angle = *angle;
+#else
+ OGS_FATAL(
+ "The code is not compiled with the Eigen unsupported modules.");
+#endif
+ }
+ if (auto residualupdate =
+ //! \ogs_file_param{prj__linear_solvers__linear_solver__eigen__residual_update}
+ config->getConfigParameterOptional<int>("residual_update"))
+ {
+#ifdef USE_EIGEN_UNSUPPORTED
+ options.residualupdate = *residualupdate;
+#else
+ OGS_FATAL(
+ "The code is not compiled with the Eigen unsupported modules.");
#endif
}
return {prefix, options};
--
GitLab