Merge branch 'FlexibilizationMatrixStorageInEigenLinearSolver' into '…

…master' Flexibilization of storage of intermediate results in EigenLinearSolver solves memory See merge request ogs/ogs!4977
ufz · Apr 12, 2024 · 06c4053 · 06c4053
2 parents 23dd666 + 0fbeeb7
commit 06c4053
Show file tree

Hide file tree

Showing 6 changed files with 107 additions and 30 deletions.
diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
@@ -62,7 +62,8 @@ class EigenLinearSolverBase
         return success;
     }
 
-    bool compute(Matrix& A, EigenOption& opt)
+    bool compute(Matrix& A, EigenOption& opt,
+                 MathLib::LinearSolverBehaviour const linear_solver_behaviour)
     {
 #ifdef USE_EIGEN_UNSUPPORTED
         if (opt.scaling)
@@ -74,13 +75,15 @@ class EigenLinearSolverBase
         }
 #endif
 
-        return computeImpl(A, opt);
+        return computeImpl(A, opt, linear_solver_behaviour);
     }
 
 protected:
     virtual bool solveImpl(Vector const& b, Vector& x, EigenOption& opt) = 0;
 
-    virtual bool computeImpl(Matrix& A, EigenOption& opt) = 0;
+    virtual bool computeImpl(
+        Matrix& A, EigenOption& opt,
+        MathLib::LinearSolverBehaviour const linear_solver_behaviour) = 0;
 
 private:
 #ifdef USE_EIGEN_UNSUPPORTED
@@ -110,7 +113,9 @@ class EigenDirectLinearSolver final : public EigenLinearSolverBase
         return true;
     }
 
-    bool computeImpl(Matrix& A, EigenOption& opt) override
+    bool computeImpl(Matrix& A, EigenOption& opt,
+                     [[maybe_unused]] MathLib::LinearSolverBehaviour const
+                         linear_solver_behaviour) override
     {
         INFO("-> compute with Eigen direct linear solver {:s}",
              EigenOption::getSolverName(opt.solver_type));
@@ -236,7 +241,9 @@ template <class T_SOLVER>
 class EigenIterativeLinearSolver final : public EigenLinearSolverBase
 {
 public:
-    bool computeImpl(Matrix& A, EigenOption& opt) override
+    bool computeImpl(
+        Matrix& A, EigenOption& opt,
+        MathLib::LinearSolverBehaviour const linear_solver_behaviour) override
     {
         INFO("-> compute with Eigen iterative linear solver {:s} (precon {:s})",
              EigenOption::getSolverName(opt.solver_type),
@@ -257,17 +264,38 @@ class EigenIterativeLinearSolver final : public EigenLinearSolverBase
             T_SOLVER>::setResidualUpdate(opt.residualupdate);
 #endif
 
-        // matrix must be copied, because Eigen's linear solver stores a
-        // reference to it cf.
-        // https://libeigen.gitlab.io/docs/classEigen_1_1IterativeSolverBase.html#a7dfa55c55e82d697bde227696a630914
-        A_ = A;
+        auto compute = [this](Matrix& m)
+        {
+            if (!m.isCompressed())
+            {
+                m.makeCompressed();
+            }
+
+            solver_.compute(m);
+        };
 
-        if (!A_.isCompressed())
+        switch (linear_solver_behaviour)
         {
-            A_.makeCompressed();
-        }
+            case MathLib::LinearSolverBehaviour::RECOMPUTE_AND_STORE:
+            {
+                // matrix must be copied, because Eigen's linear solver stores a
+                // reference to it cf.
+                // https://libeigen.gitlab.io/docs/classEigen_1_1IterativeSolverBase.html#a7dfa55c55e82d697bde227696a630914
+                A_ = A;
+                compute(A_);
+                break;
+            }
+            case MathLib::LinearSolverBehaviour::RECOMPUTE:
+            {
+                compute(A);
+                break;
+            }
+            case MathLib::LinearSolverBehaviour::REUSE:
+                OGS_FATAL(
+                    "If NumLib::LinearSolverBehaviour::REUSE is set then "
+                    "EigenLinearSolver::compute() should never be executed");
+        };
 
-        solver_.compute(A_);
         if (solver_.info() != Eigen::Success)
         {
             ERR("Failed during Eigen linear solver initialization");
@@ -466,12 +494,14 @@ EigenLinearSolver::EigenLinearSolver(std::string const& /*solver_name*/,
 
 EigenLinearSolver::~EigenLinearSolver() = default;
 
-bool EigenLinearSolver::compute(EigenMatrix& A)
+bool EigenLinearSolver::compute(
+    EigenMatrix& A,
+    MathLib::LinearSolverBehaviour const linear_solver_behaviour)
 {
     INFO("------------------------------------------------------------------");
     INFO("*** Eigen solver compute()");
 
-    return solver_->compute(A.getRawMatrix(), option_);
+    return solver_->compute(A.getRawMatrix(), option_, linear_solver_behaviour);
 }
 
 bool EigenLinearSolver::solve(EigenVector& b, EigenVector& x)
@@ -482,9 +512,12 @@ bool EigenLinearSolver::solve(EigenVector& b, EigenVector& x)
     return solver_->solve(b.getRawVector(), x.getRawVector(), option_);
 }
 
-bool EigenLinearSolver::solve(EigenMatrix& A, EigenVector& b, EigenVector& x)
+bool EigenLinearSolver::solve(
+    EigenMatrix& A, EigenVector& b, EigenVector& x,
+    MathLib::LinearSolverBehaviour const linear_solver_behaviour)
 {
-    return solver_->compute(A.getRawMatrix(), option_) &&
+    return solver_->compute(A.getRawMatrix(), option_,
+                            linear_solver_behaviour) &&
            solver_->solve(b.getRawVector(), x.getRawVector(), option_);
 }
 

diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.h b/MathLib/LinAlg/Eigen/EigenLinearSolver.h
@@ -14,6 +14,7 @@
 #include <vector>
 
 #include "EigenOption.h"
+#include "MathLib/LinAlg/LinearSolverBehaviour.h"
 
 namespace MathLib
 {
@@ -52,8 +53,8 @@ class EigenLinearSolver final
      * I.e., computes the (LU) decomposition in case of a direct solver, or
      * computes the preconditioner of an iterative solver.
      */
-    bool compute(EigenMatrix& A);
-
+    bool compute(EigenMatrix& A,
+                 MathLib::LinearSolverBehaviour const linear_solver_behaviour);
     /**
      * Solves the linear system for the given right-hand side \c b and initial
      * guess \c x.
@@ -64,7 +65,11 @@ class EigenLinearSolver final
     bool solve(EigenVector& b, EigenVector& x);
 
     /// Computes and solves in a single call.
-    bool solve(EigenMatrix& A, EigenVector& b, EigenVector& x);
+    bool solve(EigenMatrix& A,
+               EigenVector& b,
+               EigenVector& x,
+               MathLib::LinearSolverBehaviour const linear_solver_behaviour =
+                   MathLib::LinearSolverBehaviour::RECOMPUTE);
 
 protected:
     EigenOption option_;

diff --git a/MathLib/LinAlg/LinearSolverBehaviour.h b/MathLib/LinAlg/LinearSolverBehaviour.h
@@ -0,0 +1,22 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2024, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#pragma once
+
+namespace MathLib
+{
+enum class LinearSolverBehaviour : int
+{
+    RECOMPUTE,
+    RECOMPUTE_AND_STORE,
+    REUSE
+};
+
+}  // namespace MathLib
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -29,20 +29,23 @@ namespace detail
 #if !defined(USE_PETSC) && !defined(USE_LIS)
 bool solvePicard(GlobalLinearSolver& linear_solver, GlobalMatrix& A,
                  GlobalVector& rhs, GlobalVector& x,
-                 bool const compute_necessary)
+                 MathLib::LinearSolverBehaviour const linear_solver_behaviour)
 {
     BaseLib::RunTime time_linear_solver;
     time_linear_solver.start();
 
-    if (compute_necessary)
+    if (linear_solver_behaviour == MathLib::LinearSolverBehaviour::RECOMPUTE ||
+        linear_solver_behaviour ==
+            MathLib::LinearSolverBehaviour::RECOMPUTE_AND_STORE)
     {
-        if (!linear_solver.compute(A))
+        if (!linear_solver.compute(A, linear_solver_behaviour))
         {
             ERR("Picard: The linear solver failed in the compute() step.");
             return false;
         }
     }
 
+    // REUSE the previously computed preconditioner or LU decomposition
     bool const iteration_succeeded = linear_solver.solve(rhs, x);
 
     INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
@@ -58,9 +61,11 @@ bool solvePicard(GlobalLinearSolver& linear_solver, GlobalMatrix& A,
 #else
 bool solvePicard(GlobalLinearSolver& linear_solver, GlobalMatrix& A,
                  GlobalVector& rhs, GlobalVector& x,
-                 bool const compute_necessary)
+                 MathLib::LinearSolverBehaviour const linear_solver_behaviour)
 {
-    if (!compute_necessary)
+    if (linear_solver_behaviour ==
+            MathLib::LinearSolverBehaviour::RECOMPUTE_AND_STORE ||
+        linear_solver_behaviour == MathLib::LinearSolverBehaviour::REUSE)
     {
         WARN(
             "The performance optimization to skip the linear solver compute() "

diff --git a/NumLib/ODESolver/NonlinearSystem.h b/NumLib/ODESolver/NonlinearSystem.h
@@ -11,6 +11,7 @@
 #pragma once
 
 #include "EquationSystem.h"
+#include "MathLib/LinAlg/LinearSolverBehaviour.h"
 #include "Types.h"
 
 namespace NumLib
@@ -129,7 +130,8 @@ class NonlinearSystem<NonlinearSolverTag::Picard> : public EquationSystem
 
     //! Returns whether the assembled matrix \f$A\f$ has changed and the linear
     //! solver must perform the MathLib::EigenLinearSolver::compute() step.
-    virtual bool linearSolverNeedsToCompute() const = 0;
+    virtual MathLib::LinearSolverBehaviour linearSolverNeedsToCompute()
+        const = 0;
 };
 
 //! @}

diff --git a/NumLib/ODESolver/TimeDiscretizedODESystem.h b/NumLib/ODESolver/TimeDiscretizedODESystem.h
@@ -234,11 +234,21 @@ class TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
         return _ode.getMatrixSpecifications(process_id);
     }
 
-    bool linearSolverNeedsToCompute() const override
+    MathLib::LinearSolverBehaviour linearSolverNeedsToCompute() const override
     {
-        return !_ode.shouldLinearSolverComputeOnlyUponTimestepChange() ||
-               _time_disc.getCurrentTimeIncrement() !=
-                   _time_disc.getPreviousTimeIncrement();
+        if (_ode.shouldLinearSolverComputeOnlyUponTimestepChange() &&
+            _time_disc.getCurrentTimeIncrement() !=
+                _time_disc.getPreviousTimeIncrement())
+        {
+            return MathLib::LinearSolverBehaviour::RECOMPUTE_AND_STORE;
+        }
+        if (_ode.shouldLinearSolverComputeOnlyUponTimestepChange() &&
+            _time_disc.getCurrentTimeIncrement() ==
+                _time_disc.getPreviousTimeIncrement())
+        {
+            return MathLib::LinearSolverBehaviour::REUSE;
+        }
+        return MathLib::LinearSolverBehaviour::RECOMPUTE;
     }
 
 private: