From b97a21ddee2e72010a1780cf44131cce665e487f Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Thu, 2 May 2024 20:55:25 +0200
Subject: [PATCH] Try to fix accelerated linear PDE assembly
---
MathLib/LinAlg/Eigen/EigenLinearSolver.cpp | 26 +++++++++++++++++++++-
NumLib/ODESolver/NonlinearSolver.cpp | 12 +++-------
2 files changed, 28 insertions(+), 10 deletions(-)
diff --git a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
index db3ffc563bd..1f0c394355f 100644
--- a/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
+++ b/MathLib/LinAlg/Eigen/EigenLinearSolver.cpp
@@ -63,8 +63,30 @@ public:
}
bool compute(Matrix& A, EigenOption& opt,
- MathLib::LinearSolverBehaviour const linear_solver_behaviour)
+ MathLib::LinearSolverBehaviour linear_solver_behaviour)
{
+ if (linear_solver_behaviour == MathLib::LinearSolverBehaviour::REUSE)
+ {
+ // Checking if this is the first compute() call should work both for
+ // direct solvers (that factorize the matrix A) and for iterative
+ // solvers (that store a reference to A in the preconditioner).
+ if (is_first_compute_call_)
+ {
+ // There is nothing there, yet, to be reused. Hence, we compute
+ // and store the result.
+ linear_solver_behaviour =
+ MathLib::LinearSolverBehaviour::RECOMPUTE_AND_STORE;
+ }
+ else
+ {
+ return true;
+ }
+ }
+
+ // TODO (CL) That might not work under all circumstances, e.g., if the
+ // linear solver fails subsequently. Time will tell.
+ is_first_compute_call_ = false;
+
#ifdef USE_EIGEN_UNSUPPORTED
if (opt.scaling)
{
@@ -89,6 +111,8 @@ private:
#ifdef USE_EIGEN_UNSUPPORTED
std::unique_ptr<Eigen::IterScaling<EigenMatrix::RawMatrixType>> scaling_;
#endif
+
+ bool is_first_compute_call_ = true;
};
namespace details
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index 539aa469f85..a84376c2394 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -34,18 +34,12 @@ bool solvePicard(GlobalLinearSolver& linear_solver, GlobalMatrix& A,
BaseLib::RunTime time_linear_solver;
time_linear_solver.start();
- if (linear_solver_behaviour == MathLib::LinearSolverBehaviour::RECOMPUTE ||
- linear_solver_behaviour ==
- MathLib::LinearSolverBehaviour::RECOMPUTE_AND_STORE)
+ if (!linear_solver.compute(A, linear_solver_behaviour))
{
- if (!linear_solver.compute(A, linear_solver_behaviour))
- {
- ERR("Picard: The linear solver failed in the compute() step.");
- return false;
- }
+ 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());
--
GitLab