From 0344fff12ed5d39d672b20fb246318ebb33f03c3 Mon Sep 17 00:00:00 2001
From: Norihiro Watanabe <norihiro.watanabe@ufz.de>
Date: Wed, 2 Nov 2016 11:25:38 +0100
Subject: [PATCH] [Math/Eigen] support scaling
---
MathLib/LinAlg/Eigen/EigenLinearSolver.cpp | 25 ++++++++++++++++++++++
MathLib/LinAlg/Eigen/EigenOption.cpp | 3 +++
MathLib/LinAlg/Eigen/EigenOption.h | 4 ++++
3 files changed, 32 insertions(+)
diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
index 858987d7f8e..f308a482925 100644
--- a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
+++ b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
@@ -15,6 +15,10 @@
#include <Eigen/PardisoSupport>
#endif
+#ifdef USE_EIGEN_UNSUPPORTED
+#include <unsupported/Eigen/src/IterativeSolvers/Scaling.h>
+#endif
+
#include "BaseLib/ConfigTree.h"
#include "EigenVector.h"
#include "EigenMatrix.h"
@@ -224,6 +228,13 @@ void EigenLinearSolver::setOption(BaseLib::ConfigTree const& option)
ptSolver->getConfigParameterOptional<int>("max_iteration_step")) {
_option.max_iterations = *max_iteration_step;
}
+#ifdef USE_EIGEN_UNSUPPORTED
+ if (auto scaling =
+ //! \ogs_file_param{linear_solver__eigen__scaling}
+ ptSolver->getConfigParameterOptional<int>("scaling")) {
+ _option.scaling = static_cast<bool>(*scaling);
+ }
+#endif
}
bool EigenLinearSolver::solve(EigenMatrix &A, EigenVector& b, EigenVector &x)
@@ -231,8 +242,22 @@ bool EigenLinearSolver::solve(EigenMatrix &A, EigenVector& b, EigenVector &x)
INFO("------------------------------------------------------------------");
INFO("*** Eigen solver computation");
+#ifdef USE_EIGEN_UNSUPPORTED
+ std::unique_ptr<Eigen::IterScaling<EigenMatrix::RawMatrixType>> scal;
+ if (_option.scaling)
+ {
+ INFO("-> scale");
+ scal.reset(new Eigen::IterScaling<EigenMatrix::RawMatrixType>());
+ scal->computeRef(A.getRawMatrix());
+ b.getRawVector() = scal->LeftScaling().cwiseProduct(b.getRawVector());
+ }
+#endif
auto const success = _solver->solve(A.getRawMatrix(), b.getRawVector(),
x.getRawVector(), _option);
+#ifdef USE_EIGEN_UNSUPPORTED
+ if (scal)
+ x.getRawVector() = scal->RightScaling().cwiseProduct(x.getRawVector());
+#endif
INFO("------------------------------------------------------------------");
diff --git a/MathLib/LinAlg/Eigen/EigenOption.cpp b/MathLib/LinAlg/Eigen/EigenOption.cpp
index 8bd700d93b6..4fa228a245a 100644
--- a/MathLib/LinAlg/Eigen/EigenOption.cpp
+++ b/MathLib/LinAlg/Eigen/EigenOption.cpp
@@ -19,6 +19,9 @@ EigenOption::EigenOption()
precon_type = PreconType::NONE;
max_iterations = static_cast<int>(1e6);
error_tolerance = 1.e-16;
+#ifdef USE_EIGEN_UNSUPPORTED
+ scaling = false;
+#endif
}
EigenOption::SolverType EigenOption::getSolverType(const std::string &solver_name)
diff --git a/MathLib/LinAlg/Eigen/EigenOption.h b/MathLib/LinAlg/Eigen/EigenOption.h
index 70b4329908f..46ccf4d4d09 100644
--- a/MathLib/LinAlg/Eigen/EigenOption.h
+++ b/MathLib/LinAlg/Eigen/EigenOption.h
@@ -43,6 +43,10 @@ struct EigenOption final
int max_iterations;
/// Error tolerance
double error_tolerance;
+#ifdef USE_EIGEN_UNSUPPORTED
+ /// Scaling the coefficient matrix and the RHS bector
+ bool scaling;
+#endif
/// Constructor
///
--
GitLab