diff --git a/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp b/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp
index b77e150..4758e70 100644
--- a/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp
+++ b/tests/unit/solvers/linear/mkl_pardiso_solver_test.cpp
@@ -193,6 +193,32 @@ TEST(MklPardisoSolver, RefactorizesWithoutLeakingState) {
             1.0e-10);
 }
 
+TEST(MklPardisoSolver, SingleEquationScopeTransitionsBetweenRetainedSizes) {
+  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;
+  ASSERT_TRUE(solver.Factorize(MakeDenseCsr(1U, {2.0})).IsOk());
+  auto single = MakeVector({0.0});
+  ASSERT_TRUE(solver.Solve(MakeVector({4.0}), single).IsOk());
+  EXPECT_DOUBLE_EQ(single[0U], 2.0);
+  ASSERT_TRUE(solver.Factorize(MakeDenseCsr(2U, {2, 0, 0, 3})).IsOk());
+  fesa::Vector pair{2U};
+  ASSERT_TRUE(solver.Solve(MakeVector({4, 6}), pair).IsOk());
+  EXPECT_NEAR(pair[0U], 2, 1e-14);
+  EXPECT_NEAR(pair[1U], 2, 1e-14);
+  EXPECT_EQ(mkl_get_max_threads(), default_threads);
+  EXPECT_EQ(mkl_set_num_threads_local(0), 0);
+  ASSERT_TRUE(solver.Factorize(MakeDenseCsr(1U, {4.0})).IsOk());
+  ASSERT_TRUE(solver.Solve(MakeVector({8.0}), single).IsOk());
+  EXPECT_DOUBLE_EQ(single[0U], 2.0);
+  EXPECT_EQ(mkl_get_max_threads(), default_threads);
+  EXPECT_EQ(mkl_set_num_threads_local(0), 0);
+}
+
 TEST(MklPardisoSolver, ClassifiesSingularIndefiniteAndNonfiniteFailures) {
   fesa::MklPardisoSolver solver;
 
diff --git a/tests/unit/solvers/linear/tangent_solver_test.cpp b/tests/unit/solvers/linear/tangent_solver_test.cpp
index 16795e7..bb7096f 100644
--- a/tests/unit/solvers/linear/tangent_solver_test.cpp
+++ b/tests/unit/solvers/linear/tangent_solver_test.cpp
@@ -205,4 +205,30 @@ TEST(TangentSolver, SingleEquationScopeRestoresTlsAcrossFreshFailureAndReuse) {
   EXPECT_TRUE(second.get());
 }
 
+TEST(TangentSolver, SingleEquationScopeTransitionsBetweenRetainedSizes) {
+  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;
+  ASSERT_TRUE(solver.FactorizeFresh(MakeDenseCsr(1U, {2.0})).IsOk());
+  auto single = MakeVector({0.0});
+  ASSERT_TRUE(solver.SolveCorrection(MakeVector({4.0}), single).IsOk());
+  EXPECT_DOUBLE_EQ(single[0U], 2.0);
+  ASSERT_TRUE(solver.FactorizeFresh(MakeDenseCsr(2U, {2, 1, 0, 3})).IsOk());
+  fesa::Vector pair{2U};
+  ASSERT_TRUE(solver.SolveCorrection(MakeVector({6, 6}), pair).IsOk());
+  EXPECT_NEAR(pair[0U], 2, 1e-14);
+  EXPECT_NEAR(pair[1U], 2, 1e-14);
+  EXPECT_EQ(mkl_get_max_threads(), default_threads);
+  EXPECT_EQ(mkl_set_num_threads_local(0), 0);
+  ASSERT_TRUE(solver.FactorizeFresh(MakeDenseCsr(1U, {4.0})).IsOk());
+  ASSERT_TRUE(solver.SolveCorrection(MakeVector({8.0}), single).IsOk());
+  EXPECT_DOUBLE_EQ(single[0U], 2.0);
+  EXPECT_EQ(mkl_get_max_threads(), default_threads);
+  EXPECT_EQ(mkl_set_num_threads_local(0), 0);
+}
+
 }  // namespace