diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
index 858987d7f8e9c750b84e7f334fdbe355e1db0595..f308a482925b8e06151026be529e855ec0c4d448 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 8bd700d93b6298ad3b1628476725a6ba1881bb4a..4fa228a245ac09319de712084eae0887b2bb73e6 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 70b4329908f1dabcf75e189fb2e567c147f68e5a..46ccf4d4d092147e5fecd0490b1db4a9e7faab50 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
     ///