perf(solvers): scope default one-equation calls to one thread

src/fesa/math/dense_blas_internal.cpp

@@ -8,6 +8,24 @@
 
 namespace fesa::dense_blas_internal {
 
+ScopedSingleEquationThreads::ScopedSingleEquationThreads(
+    const std::size_t equation_count) noexcept {
+  // The measured 1x1 PARDISO case has backend/thread overhead and an
+  // intermittently failing default-thread lifecycle. Broader sizes are not
+  // inferred from this auxiliary evidence.
+  if (equation_count != 1U) return;
+  const int previous = mkl_set_num_threads_local(1);
+  if (previous == 0) {
+    restore_default_ = true;
+  } else {
+    static_cast<void>(mkl_set_num_threads_local(previous));
+  }
+}
+
+ScopedSingleEquationThreads::~ScopedSingleEquationThreads() {
+  if (restore_default_) static_cast<void>(mkl_set_num_threads_local(0));
+}
+
 bool CholeskyLowerAndInverse(const Matrix& matrix, Matrix& lower,
                             Matrix& inverse_lower) {
   if (matrix.Rows() != matrix.Columns())

src/fesa/math/dense_blas_internal.h

@@ -11,6 +11,21 @@
 
 namespace fesa::dense_blas_internal {
 
+/// @brief Limits only default-policy one-equation backend calls to one thread.
+/// @note Existing caller TLS overrides are retained. Destruction restores the
+/// default request; process-wide MKL/OpenMP settings are never changed.
+class ScopedSingleEquationThreads {
+ public:
+  explicit ScopedSingleEquationThreads(std::size_t equation_count) noexcept;
+  ScopedSingleEquationThreads(const ScopedSingleEquationThreads&) = delete;
+  ScopedSingleEquationThreads& operator=(const ScopedSingleEquationThreads&) =
+      delete;
+  ~ScopedSingleEquationThreads();
+
+ private:
+  bool restore_default_{false};
+};
+
 /// @brief Converts a dense storage length to the private MKL integer type.
 /// @return The converted length or an uncategorized overflow failure.
 Result<MKL_INT> ToMklSize(std::size_t size);

src/fesa/solvers/linear/mkl_pardiso_solver.cpp

@@ -14,6 +14,7 @@
 
 #include "fesa/math/sparse_matrix.h"
 #include "fesa/math/vector.h"
+#include "math/dense_blas_internal.h"
 
 namespace fesa {
 namespace {
@@ -286,6 +287,8 @@ class MklPardisoSolver::Impl {
   /// @brief Calls PARDISO with retained private arrays and phase state.
   void CallPardiso(const MKL_INT phase, double* rhs, double* solution,
                    MKL_INT& error) {
+    const dense_blas_internal::ScopedSingleEquationThreads thread_scope{
+        static_cast<std::size_t>(equation_count_)};
     pardiso(pt_.data(), &max_factorizations_, &matrix_number_, &mtype_, &phase,
             &equation_count_, values_.data(), row_offsets_.data(),
             column_indices_.data(), permutation_.data(), &rhs_count_,

src/fesa/solvers/linear/mkl_pardiso_tangent_solver.cpp

@@ -13,6 +13,7 @@
 
 #include "fesa/math/sparse_matrix.h"
 #include "fesa/math/vector.h"
+#include "math/dense_blas_internal.h"
 
 namespace fesa {
 namespace {
@@ -230,6 +231,8 @@ class MklPardisoTangentSolver::Impl {
   /// @brief Calls PARDISO with retained private arrays and phase state.
   void CallPardiso(const MKL_INT phase, double* rhs, double* correction,
                    MKL_INT& error) {
+    const dense_blas_internal::ScopedSingleEquationThreads thread_scope{
+        static_cast<std::size_t>(equation_count_)};
     pardiso(pt_.data(), &max_factorizations_, &matrix_number_, &mtype_, &phase,
             &equation_count_, values_.data(), row_offsets_.data(),
             column_indices_.data(), permutation_.data(), &rhs_count_,

tests/unit/math/dense_blas_internal_test.cpp

@@ -108,5 +108,36 @@ TEST(DenseBlasInternal, CholeskyAndInverseMatchKnownFactorsAndKeepFailedOutput)
                std::invalid_argument);
 }
 
+TEST(DenseBlasInternal, SingleEquationScopeRestoresTlsAndRespectsOverrides) {
+  const int previous = mkl_set_num_threads_local(0);
+  struct Restore {
+    int previous;
+    ~Restore() { static_cast<void>(mkl_set_num_threads_local(previous)); }
+  } restore{previous};
+  const int default_threads = mkl_get_max_threads();
+  {
+    const ScopedSingleEquationThreads scope{1U};
+    EXPECT_EQ(mkl_get_max_threads(), 1);
+  }
+  EXPECT_EQ(mkl_get_max_threads(), default_threads);
+  EXPECT_EQ(mkl_set_num_threads_local(0), 0);
+  {
+    const ScopedSingleEquationThreads scope{2U};
+    EXPECT_EQ(mkl_get_max_threads(), default_threads);
+  }
+  try {
+    const ScopedSingleEquationThreads scope{1U};
+    throw std::runtime_error{"scope failure"};
+  } catch (const std::runtime_error&) {
+  }
+  EXPECT_EQ(mkl_set_num_threads_local(0), 0);
+  static_cast<void>(mkl_set_num_threads_local(3));
+  {
+    const ScopedSingleEquationThreads scope{1U};
+    EXPECT_EQ(mkl_get_max_threads(), 3);
+  }
+  EXPECT_EQ(mkl_get_max_threads(), 3);
+}
+
 }  // namespace
 }  // namespace fesa::dense_blas_internal

tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp

@@ -3,6 +3,7 @@
 #include <gtest/gtest.h>
 
 #include <cmath>
+#include <future>
 #include <initializer_list>
 #include <limits>
 #include <string>
@@ -11,6 +12,7 @@
 
 #include "fesa/fem/dof_manager.h"
 #include "fesa/math/sparse_matrix.h"
+#include "math/dense_blas_internal.h"
 
 namespace {
 
@@ -123,6 +125,42 @@ TEST(MklPardisoSolver, ReusesOneFactorizationForRepeatedRhs) {
   EXPECT_LE(NormalizedResidual(matrix, second, rhs_second), 1.0e-10);
 }
 
+TEST(MklPardisoSolver, SingleEquationPolicyPreservesLifecycleAndConcurrentTls) {
+  const auto worker = [] {
+    const int previous = mkl_set_num_threads_local(0);
+    struct Restore {
+      int previous;
+      ~Restore() { static_cast<void>(mkl_set_num_threads_local(previous)); }
+    } restore{previous};
+    const int default_threads = mkl_get_max_threads();
+    fesa::MklPardisoSolver solver;
+    for (int repeat = 0; repeat < 10; ++repeat) {
+      if (!solver.Factorize(MakeDenseCsr(1U, {2.0})).IsOk()) return false;
+      auto solution = MakeVector({9.0});
+      if (!solver.Solve(MakeVector({4.0}), solution).IsOk() ||
+          solution[0U] != 2.0) return false;
+      if (!solver.Factorize(MakeDenseCsr(1U, {4.0})).IsOk()) return false;
+      if (!solver.Solve(MakeVector({4.0}), solution).IsOk() ||
+          solution[0U] != 1.0) return false;
+      auto invalid_rhs = MakeVector({std::numeric_limits<double>::quiet_NaN()});
+      if (solver.Solve(invalid_rhs, solution).IsOk() || solution[0U] != 1.0)
+        return false;
+      if (!solver.Solve(MakeVector({8.0}), solution).IsOk() ||
+          solution[0U] != 2.0) return false;
+      if (solver.Factorize(MakeDenseCsr(1U, {-1.0})).IsOk()) return false;
+      if (solver.Solve(MakeVector({8.0}), solution).IsOk() ||
+          solution[0U] != 2.0) return false;
+      if (mkl_get_max_threads() != default_threads ||
+          mkl_set_num_threads_local(0) != 0) return false;
+    }
+    return true;
+  };
+  auto first = std::async(std::launch::async, worker);
+  auto second = std::async(std::launch::async, worker);
+  EXPECT_TRUE(first.get());
+  EXPECT_TRUE(second.get());
+}
+
 TEST(MklPardisoSolver, StillRejectsNonsymmetricPublicCsr) {
   const auto nonsymmetric = MakeDenseCsr(2U, {2.0, 1.0, 0.0, 3.0});
   fesa::MklPardisoSolver solver;

tests/unit/solvers/linear/tangent_solver_test.cpp

@@ -3,6 +3,7 @@
 #include <gtest/gtest.h>
 
 #include <cstddef>
+#include <future>
 #include <initializer_list>
 #include <limits>
 #include <type_traits>
@@ -13,6 +14,7 @@
 #include "fesa/math/sparse_matrix.h"
 #include "fesa/math/vector.h"
 #include "fesa/solvers/linear/tangent_solver.h"
+#include "math/dense_blas_internal.h"
 
 namespace {
 
@@ -169,4 +171,38 @@ TEST(TangentSolver, ReusesRhsButFreshValuesAndFailureCannotReuseStaleState) {
   EXPECT_TRUE(solver.SolveCorrection(MakeVector({10, -7}), correction).IsOk());
 }
 
+TEST(TangentSolver, SingleEquationScopeRestoresTlsAcrossFreshFailureAndReuse) {
+  const auto worker = [] {
+    const int previous = mkl_set_num_threads_local(0);
+    struct Restore {
+      int previous;
+      ~Restore() { static_cast<void>(mkl_set_num_threads_local(previous)); }
+    } restore{previous};
+    const int default_threads = mkl_get_max_threads();
+    fesa::MklPardisoTangentSolver solver;
+    for (int repeat = 0; repeat < 10; ++repeat) {
+      auto correction = MakeVector({17.0});
+      if (!solver.FactorizeFresh(MakeDenseCsr(1U, {2.0})).IsOk()) return false;
+      if (!solver.SolveCorrection(MakeVector({4.0}), correction).IsOk() ||
+          correction[0U] != 2.0) return false;
+      if (!solver.FactorizeFresh(MakeDenseCsr(1U, {4.0})).IsOk()) return false;
+      if (!solver.SolveCorrection(MakeVector({4.0}), correction).IsOk() ||
+          correction[0U] != 1.0) return false;
+      if (solver.SolveCorrection(MakeVector({
+              std::numeric_limits<double>::infinity()}), correction).IsOk() ||
+          correction[0U] != 1.0) return false;
+      if (solver.FactorizeFresh(MakeDenseCsr(1U, {0.0})).IsOk()) return false;
+      if (solver.SolveCorrection(MakeVector({4.0}), correction).IsOk() ||
+          correction[0U] != 1.0) return false;
+      if (mkl_get_max_threads() != default_threads ||
+          mkl_set_num_threads_local(0) != 0) return false;
+    }
+    return true;
+  };
+  auto first = std::async(std::launch::async, worker);
+  auto second = std::async(std::launch::async, worker);
+  EXPECT_TRUE(first.get());
+  EXPECT_TRUE(second.get());
+}
+
 }  // namespace