From 770fa2a0d1e162bc710646934a4ee79af030f4e4 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Wed, 14 Oct 2026 15:56:03 +0900 Subject: [PATCH] perf(solvers): scope default one-equation calls to one thread --- src/fesa/math/dense_blas_internal.cpp | 18 +++++++++ src/fesa/math/dense_blas_internal.h | 15 ++++++++ .../solvers/linear/mkl_pardiso_solver.cpp | 3 ++ .../linear/mkl_pardiso_tangent_solver.cpp | 3 ++ tests/unit/math/dense_blas_internal_test.cpp | 31 +++++++++++++++ .../linear/mkl_pardiso_solver_test.cpp | 38 +++++++++++++++++++ .../solvers/linear/tangent_solver_test.cpp | 36 ++++++++++++++++++ 7 files changed, 144 insertions(+) diff --git a/src/fesa/math/dense_blas_internal.cpp b/src/fesa/math/dense_blas_internal.cpp index 40d53d4..413a71c 100644 --- a/src/fesa/math/dense_blas_internal.cpp +++ b/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(mkl_set_num_threads_local(previous)); + } +} + +ScopedSingleEquationThreads::~ScopedSingleEquationThreads() { + if (restore_default_) static_cast(mkl_set_num_threads_local(0)); +} + bool CholeskyLowerAndInverse(const Matrix& matrix, Matrix& lower, Matrix& inverse_lower) { if (matrix.Rows() != matrix.Columns()) diff --git a/src/fesa/math/dense_blas_internal.h b/src/fesa/math/dense_blas_internal.h index 21e7749..5934506 100644 --- a/src/fesa/math/dense_blas_internal.h +++ b/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 ToMklSize(std::size_t size); diff --git a/src/fesa/solvers/linear/mkl_pardiso_solver.cpp b/src/fesa/solvers/linear/mkl_pardiso_solver.cpp index add8672..0cc9adc 100644 --- a/src/fesa/solvers/linear/mkl_pardiso_solver.cpp +++ b/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(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_, diff --git a/src/fesa/solvers/linear/mkl_pardiso_tangent_solver.cpp b/src/fesa/solvers/linear/mkl_pardiso_tangent_solver.cpp index 1554b67..b544312 100644 --- a/src/fesa/solvers/linear/mkl_pardiso_tangent_solver.cpp +++ b/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(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_, diff --git a/tests/unit/math/dense_blas_internal_test.cpp b/tests/unit/math/dense_blas_internal_test.cpp index 6e32e17..c365287 100644 --- a/tests/unit/math/dense_blas_internal_test.cpp +++ b/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(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(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 diff --git a/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp b/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp index 1a6953e..b77e150 100644 --- a/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp +++ b/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -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(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::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; diff --git a/tests/unit/solvers/linear/tangent_solver_test.cpp b/tests/unit/solvers/linear/tangent_solver_test.cpp index 955d368..16795e7 100644 --- a/tests/unit/solvers/linear/tangent_solver_test.cpp +++ b/tests/unit/solvers/linear/tangent_solver_test.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -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(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::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