diff --git a/include/fesa/math/sparse_matrix.h b/include/fesa/math/sparse_matrix.h
index dedf73e..bbebfc1 100644
--- a/include/fesa/math/sparse_matrix.h
+++ b/include/fesa/math/sparse_matrix.h
@@ -46,8 +46,12 @@ class SparseMatrix {
   /// @brief Returns CSR values including preserved structural zeros.
   const std::vector<double>& Values() const noexcept;
 
-  /// @brief Multiplies this matrix by a dense vector in stable CSR order.
-  /// @throws std::invalid_argument if the dimensions are incompatible.
+  /// @brief Multiplies canonical CSR by a dense vector through private MKL.
+  /// @throws std::invalid_argument if dimensions or CSR are invalid.
+  /// @throws std::length_error if a size exceeds the backend integer range.
+  /// @throws std::runtime_error if the sparse backend fails.
+  /// @note Calls retain no handle; concurrent const calls own their resources.
+  /// Nonfinite products retain IEEE values for consumer finite validation.
   Vector Multiply(const Vector& rhs) const;
 
   /// @brief Validates shape, indices, ordering, and finite CSR values.
diff --git a/src/fesa/analysis/eigen_reduction.cpp b/src/fesa/analysis/eigen_reduction.cpp
index c02514e..7edb030 100644
--- a/src/fesa/analysis/eigen_reduction.cpp
+++ b/src/fesa/analysis/eigen_reduction.cpp
@@ -10,6 +10,7 @@
 #include <vector>
 
 #include "fesa/fem/dof_manager.h"
+#include "math/dense_blas_internal.h"
 
 namespace fesa::eigen_internal {
 namespace {
@@ -117,7 +118,8 @@ Result<NodalLift> MakeLift(const Matrix& mass,
   if (positive == 0U) return Result<NodalLift>::Success(std::move(result));
 
   const Matrix scaled = ScaledMass(mass, length);
-  Matrix mp(positive, positive);
+  Matrix mp = dense_blas_internal::MultiplyTransposeLeft(basis, scaled)
+                  .Multiply(basis);
   double scale = 0.0;
   double basis_error = 0.0;
   for (std::size_t i = 0; i < positive; ++i) {
@@ -125,9 +127,6 @@ Result<NodalLift> MakeLift(const Matrix& mass,
       double gram = 0.0;
       for (std::size_t a = 0; a < 6U; ++a) {
         gram += basis(a, i) * basis(a, j);
-        for (std::size_t b = 0; b < 6U; ++b) {
-          mp(i, j) += basis(a, i) * scaled(a, b) * basis(b, j);
-        }
       }
       basis_error = std::hypot(basis_error, gram - (i == j ? 1.0 : 0.0));
       scale = (std::max)(scale, std::abs(mp(i, j)));
@@ -343,17 +342,18 @@ Result<EigenReduction> BuildEigenReduction(
         kz_error = std::hypot(kz_error, item.second);
     }
     const auto& left = lifts[node];
+    std::map<std::size_t, Matrix> projected_blocks;
+    for (const auto& entry : blocks) {
+      projected_blocks.emplace(
+          entry.first, dense_blas_internal::MultiplyTransposeLeft(
+                           left.scaled_lift, entry.second)
+                           .Multiply(lifts[entry.first].scaled_lift));
+    }
     for (std::size_t i = 0; i < left.scaled_lift.Columns(); ++i) {
-      for (const auto& entry : blocks) {
+      for (const auto& entry : projected_blocks) {
         const auto& right = lifts[entry.first];
         for (std::size_t j = 0; j < right.scaled_lift.Columns(); ++j) {
-          double value = 0.0;
-          for (std::size_t a = 0; a < 6U; ++a) {
-            for (std::size_t b = 0; b < 6U; ++b) {
-              value += left.scaled_lift(a, i) * entry.second(a, b) *
-                       right.scaled_lift(b, j);
-            }
-          }
+          const double value = entry.second(i, j);
           pattern.column_indices.push_back(right.offset + j);
           contributions.push_back({left.offset + i, right.offset + j, value,
                                    node, pattern.column_indices.size() - 1U});
@@ -385,11 +385,12 @@ Vector RestoreScaledMode(const EigenReduction& reduction, const Matrix& vectors,
   Vector result(reduction.nodes.size() * 6U);
   for (std::size_t node = 0; node < reduction.nodes.size(); ++node) {
     const auto& lift = reduction.nodes[node];
+    Vector reduced{lift.scaled_lift.Columns()};
+    for (std::size_t j = 0; j < reduced.Size(); ++j)
+      reduced[j] = vectors(lift.offset + j, column);
+    const auto restored = lift.scaled_lift.Multiply(reduced);
     for (std::size_t a = 0; a < 6U; ++a) {
-      for (std::size_t j = 0; j < lift.scaled_lift.Columns(); ++j) {
-        result[6U * node + a] +=
-            lift.scaled_lift(a, j) * vectors(lift.offset + j, column);
-      }
+      result[6U * node + a] = restored[a];
     }
   }
   return result;
@@ -399,10 +400,11 @@ Vector MultiplyNodalMass(const std::vector<Matrix>& mass,
                          const Vector& vector) {
   Vector result(vector.Size());
   for (std::size_t node = 0; node < mass.size(); ++node) {
+    Vector local{6U};
+    for (std::size_t b = 0; b < 6U; ++b) local[b] = vector[6U * node + b];
+    const auto product = mass[node].Multiply(local);
     for (std::size_t a = 0; a < 6U; ++a) {
-      for (std::size_t b = 0; b < 6U; ++b) {
-        result[6U * node + a] += mass[node](a, b) * vector[6U * node + b];
-      }
+      result[6U * node + a] = product[a];
     }
   }
   return result;
diff --git a/src/fesa/analysis/linear_frequency_procedure.cpp b/src/fesa/analysis/linear_frequency_procedure.cpp
index 1064c4e..c886871 100644
--- a/src/fesa/analysis/linear_frequency_procedure.cpp
+++ b/src/fesa/analysis/linear_frequency_procedure.cpp
@@ -32,11 +32,7 @@ Status Failure(const std::string& code, const std::string& identity,
 }
 
 double Dot(const Vector& left, const Vector& right) {
-  double result = 0.0;
-  for (std::size_t index = 0U; index < left.Size(); ++index) {
-    result += left[index] * right[index];
-  }
-  return result;
+  return left.Dot(right);
 }
 
 double Norm(const Vector& vector) {
@@ -318,8 +314,13 @@ Result<EigenAnalysisState> ExecuteCandidate(const Domain& domain) {
                                     (pencil_scale * y_norm) * y_norm);
     const double rayleigh = energy / modal_mass;
     Vector x(a.Rows());
-    for (std::size_t row = 0; row < x.Size(); ++row)
+    for (std::size_t row = 0; row < x.Size(); ++row) {
       x[row] = state.backend.eigenvectors(row, column) * normalization;
+      if (!std::isfinite(x[row]))
+        return Result<EigenAnalysisState>::Failure(
+            Failure("invalid-eigenvalue-result", std::to_string(column),
+                    "Recovered mode is nonfinite or violates constraints."));
+    }
     Vector ax = a.Multiply(x);
     Vector whitened_residual = ax;
     whitened_residual.Axpy(-lambda, x);
@@ -421,6 +422,10 @@ Result<EigenAnalysisState> LinearFrequencyProcedureExecutor::Execute(
     return Result<EigenAnalysisState>::Failure(
         Failure("eigen-index-overflow", "modal-workspace",
                 "Requested workspace exceeds container capacity."));
+  } catch (const std::runtime_error&) {
+    return Result<EigenAnalysisState>::Failure(
+        Failure("eigen-multiplication-failed", "modal-workspace",
+                "Modal sparse multiplication failed."));
   }
 }
 
diff --git a/src/fesa/analysis/solid_frequency_procedure.cpp b/src/fesa/analysis/solid_frequency_procedure.cpp
index 0c92757..f3a33dc 100644
--- a/src/fesa/analysis/solid_frequency_procedure.cpp
+++ b/src/fesa/analysis/solid_frequency_procedure.cpp
@@ -248,8 +248,13 @@ Result<SolidFrequencyState> ExecuteCandidate(const Domain& domain) {
   mass_phi.reserve(count);
   for (std::size_t c = 0; c < count; ++c) {
     Vector vector(k.Rows());
-    for (std::size_t row = 0; row < k.Rows(); ++row)
+    for (std::size_t row = 0; row < k.Rows(); ++row) {
       vector[row] = state.backend.eigenvectors(row, c);
+      if (!std::isfinite(vector[row]))
+        return Result<SolidFrequencyState>::Failure(
+            Failure("solid-eigenpair-validation-failed",
+                    "Modal mass is not finite positive."));
+    }
     Vector mv = m.Multiply(vector);
     const double raw_mass = vector.Dot(mv);
     if (!(raw_mass > 0.0) || !std::isfinite(raw_mass))
@@ -257,6 +262,12 @@ Result<SolidFrequencyState> ExecuteCandidate(const Domain& domain) {
           Failure("solid-eigenpair-validation-failed",
                   "Modal mass is not finite positive."));
     vector.Scale(1.0 / std::sqrt(raw_mass));
+    for (std::size_t row = 0; row < vector.Size(); ++row) {
+      if (!std::isfinite(vector[row]))
+        return Result<SolidFrequencyState>::Failure(
+            Failure("solid-eigenpair-validation-failed",
+                    "Normalized modal vector is not finite."));
+    }
     mv = m.Multiply(vector);
     const Vector kv = k.Multiply(vector);
     const double lambda = state.backend.eigenvalues[c];
@@ -362,6 +373,10 @@ Result<SolidFrequencyState> SolidFrequencyProcedureExecutor::Execute(
     return Result<SolidFrequencyState>::Failure(
         Failure("eigen-allocation-failed",
                 "Solid modal workspace exceeds addressable storage."));
+  } catch (const std::runtime_error&) {
+    return Result<SolidFrequencyState>::Failure(
+        Failure("eigen-multiplication-failed",
+                "Solid modal sparse multiplication failed."));
   }
 }
 }  // namespace fesa
diff --git a/src/fesa/analysis/thermal_transient_algebra.cpp b/src/fesa/analysis/thermal_transient_algebra.cpp
index 7552e83..be029fa 100644
--- a/src/fesa/analysis/thermal_transient_algebra.cpp
+++ b/src/fesa/analysis/thermal_transient_algebra.cpp
@@ -3,6 +3,7 @@
 #include <algorithm>
 #include <cmath>
 #include <limits>
+#include <stdexcept>
 #include <string>
 #include <utility>
 
@@ -116,6 +117,14 @@ Result<Equations> Build(const ThermalSystem& s, const ThermalPartition& p,
     constrained[p.constrained[i]] = i;
   SparsePattern pattern{{0}, {}};
   std::vector<CooContribution> coo;
+  SparsePattern rhs_pattern{{0U}, {}};
+  std::vector<CooContribution> capacity_rhs, conductivity_rhs;
+  Vector capacity_input = old;
+  Vector prescribed_full{m};
+  std::vector<bool> capacity_column_used(m, false);
+  for (std::size_t i = 0; i < p.constrained.size(); ++i) {
+    prescribed_full[p.constrained[i]] = p.prescribed[i];
+  }
   Vector b{n};
   std::vector<I> exact_a, exact_b(n);
   for (std::size_t f = 0; f < n; ++f) {
@@ -126,22 +135,43 @@ Result<Equations> Build(const ThermalSystem& s, const ThermalPartition& p,
       const auto col = s.c.ColumnIndices()[j];
       const double c = s.c.Values()[j], k = s.k.Values()[j];
       const I cd = Div(I{c}, I{delta});
+      rhs_pattern.column_indices.push_back(col);
+      capacity_column_used[col] = true;
+      capacity_rhs.push_back({f, col, c / delta, 0U, j});
+      conductivity_rhs.push_back({f, col, k, 0U, j});
       if (free[col] != n) {
         coo.push_back({f, free[col], c / delta + k, 0, coo.size()});
         pattern.column_indices.push_back(free[col]);
         exact_a.push_back(Add(cd, I{k}));
-        b[f] += (c / delta) * old[col];
         exact_b[f] = Add(exact_b[f], Mul(cd, I{old[col]}));
       } else {
         const double g = p.prescribed[constrained[col]];
         // Preserve raw old constrained capacity; never project frame zero.
-        b[f] += (c / delta) * (old[col] - g) - k * g;
         exact_b[f] = Add(exact_b[f],
                          Sub(Mul(cd, Sub(I{old[col]}, I{g})), Mul(I{k}, I{g})));
       }
     }
     pattern.row_offsets.push_back(coo.size());
+    rhs_pattern.row_offsets.push_back(capacity_rhs.size());
   }
+  // Do not create overflow in an absent constrained column: the original
+  // free-row RHS never evaluated that subtraction. Referenced structural
+  // zeros still evaluate it and reject nonfinite intermediate values.
+  for (std::size_t i = 0; i < p.constrained.size(); ++i) {
+    const auto column = p.constrained[i];
+    if (capacity_column_used[column]) capacity_input[column] -= p.prescribed[i];
+  }
+  if (!Finite(capacity_input))
+    return Result<Equations>::Failure(Failure("Nonfinite backward Euler RHS."));
+  auto c_rhs = SparseMatrix::FromCoo(n, m, std::move(capacity_rhs), rhs_pattern);
+  auto k_rhs = SparseMatrix::FromCoo(n, m, std::move(conductivity_rhs),
+                                     rhs_pattern);
+  if (!c_rhs.HasValue() || !k_rhs.HasValue())
+    return Result<Equations>::Failure(Failure("Nonfinite backward Euler RHS."));
+  const auto capacity_action = c_rhs.Value().Multiply(capacity_input);
+  const auto conductivity_action = k_rhs.Value().Multiply(prescribed_full);
+  for (std::size_t f = 0; f < n; ++f)
+    b[f] += capacity_action[f] - conductivity_action[f];
   auto a = SparseMatrix::FromCoo(n, n, std::move(coo), pattern);
   if (!a.HasValue()) return Result<Equations>::Failure(a.GetStatus());
   if (!Finite(b))
@@ -178,13 +208,15 @@ Status Check(const ThermalSystem& s, const ThermalPartition& p,
   I a_norm, t_norm, b_norm = Norm(eq.b);
   Vector red{n}, red_error{n};
   std::vector<I> construction(n);
+  Vector free_temperature{n};
+  for (std::size_t f = 0; f < n; ++f)
+    free_temperature[f] = candidate.temperature[p.free[f]];
+  const auto reduced_action = eq.a.Multiply(free_temperature);
   for (std::size_t f = 0; f < n; ++f) {
     I row_norm, scale, ea, identity;
-    double dot = 0;
     for (auto j = eq.a.RowOffsets()[f]; j < eq.a.RowOffsets()[f + 1]; ++j) {
       const double a = eq.a.Values()[j],
                    t = candidate.temperature[p.free[eq.a.ColumnIndices()[j]]];
-      dot += a * t;
       row_norm = Add(row_norm, Magnitude(a));
       scale = Add(scale, Mul(Magnitude(a), Magnitude(t)));
       const I error{Error(eq.exact_a[j], a)};
@@ -193,7 +225,7 @@ Status Check(const ThermalSystem& s, const ThermalPartition& p,
     }
     a_norm = Max(a_norm, row_norm);
     t_norm = Max(t_norm, Magnitude(candidate.temperature[p.free[f]]));
-    red[f] = dot - eq.b[f];
+    red[f] = reduced_action[f] - eq.b[f];
     scale = Add(scale, Magnitude(eq.b[f]));
     const double nu =
         static_cast<double>(eq.a.RowOffsets()[f + 1] - eq.a.RowOffsets()[f]);
@@ -216,18 +248,19 @@ Status Check(const ThermalSystem& s, const ThermalPartition& p,
   for (std::size_t i = 0; i < m; ++i)
     change[i] = candidate.temperature[i] - old[i];
   if (!Finite(change)) return Failure("Nonfinite temperature change.");
+  const auto capacity_change = s.c.Multiply(change);
+  const auto conductivity_new = s.k.Multiply(candidate.temperature);
   I c_change, k_new, q_abs, c_temperatures, underflow_sum;
   std::vector<I> k_columns(m);
   double energy = 0;
   for (std::size_t i = 0; i < m; ++i) {
-    double cd = 0, kt = 0, ksum = 0;
+    const double cd = capacity_change[i], kt = conductivity_new[i];
+    double ksum = 0;
     I cs, ks, ts, cabs, nullscale, nullerror, coefficient;
     for (auto j = s.c.RowOffsets()[i]; j < s.c.RowOffsets()[i + 1]; ++j) {
       const auto col = s.c.ColumnIndices()[j];
       const double c = s.c.Values()[j], k = s.k.Values()[j],
                    t = candidate.temperature[col];
-      cd += c * change[col];
-      kt += k * t;
       ksum += k;
       cs = Add(cs, Mul(Magnitude(c), Magnitude(change[col])));
       ks = Add(ks, Mul(Magnitude(k), Magnitude(t)));
@@ -334,10 +367,13 @@ Status Check(const ThermalSystem& s, const ThermalPartition& p,
   if (p.constrained.empty() && Zero(q_abs)) {
     ev.dissipation_checked = true;
     double cn = 0, co = 0, cd = 0, kn = 0;
+    const auto capacity_new = s.c.Multiply(candidate.temperature);
+    const auto capacity_old = s.c.Multiply(old);
     I sc, subscale, t_one, d_square, o_square, skew_square, quadratic_underflow;
     for (std::size_t i = 0; i < m; ++i) {
       const double t = candidate.temperature[i];
-      double row_cn = 0, row_co = 0, row_cd = 0, row_kn = 0;
+      const double row_cn = capacity_new[i], row_co = capacity_old[i],
+                   row_cd = capacity_change[i], row_kn = conductivity_new[i];
       I scale_cn, scale_co, scale_cd, scale_kn;
       t_one = Add(t_one, Magnitude(t));
       d_square = Add(d_square, Mul(Magnitude(change[i]), Magnitude(change[i])));
@@ -345,10 +381,6 @@ Status Check(const ThermalSystem& s, const ThermalPartition& p,
       for (auto j = s.c.RowOffsets()[i]; j < s.c.RowOffsets()[i + 1]; ++j) {
         const auto col = s.c.ColumnIndices()[j];
         const double c = s.c.Values()[j], k = s.k.Values()[j];
-        row_cn += c * candidate.temperature[col];
-        row_co += c * old[col];
-        row_cd += c * change[col];
-        row_kn += k * candidate.temperature[col];
         scale_cn = Add(
             scale_cn, Mul(Magnitude(c), Magnitude(candidate.temperature[col])));
         scale_co = Add(scale_co, Mul(Magnitude(c), Magnitude(old[col])));
@@ -421,7 +453,7 @@ Status Check(const ThermalSystem& s, const ThermalPartition& p,
 
 Result<ThermalStepCandidate> SolveThermalIncrement(
     const ThermalSystem& s, const ThermalPartition& p, const Vector& old,
-    double delta, LinearSolver& solver, ThermalFactorization& cache) {
+    double delta, LinearSolver& solver, ThermalFactorization& cache) try {
   auto status = Validate(s, p, old, delta);
   if (!status.IsOk()) return Result<ThermalStepCandidate>::Failure(status);
   auto equations = Build(s, p, old, delta);
@@ -461,6 +493,12 @@ Result<ThermalStepCandidate> SolveThermalIncrement(
   status = Check(s, p, old, delta, equations.Value(), candidate);
   if (!status.IsOk()) return Result<ThermalStepCandidate>::Failure(status);
   return Result<ThermalStepCandidate>::Success(std::move(candidate));
+} catch (const std::length_error&) {
+  return Result<ThermalStepCandidate>::Failure(
+      Failure("Thermal sparse action exceeds the backend integer range."));
+} catch (const std::runtime_error&) {
+  return Result<ThermalStepCandidate>::Failure(
+      Failure("Thermal sparse action failed."));
 }
 
 Status RunThermalHistory(const ThermalAssembly& assembly,
diff --git a/src/fesa/elements/euler_beam_3d.cpp b/src/fesa/elements/euler_beam_3d.cpp
index ea43ced..82ca250 100644
--- a/src/fesa/elements/euler_beam_3d.cpp
+++ b/src/fesa/elements/euler_beam_3d.cpp
@@ -12,6 +12,7 @@
 
 #include "fesa/math/vector3.h"
 #include "fesa/model/model_types.h"
+#include "math/dense_blas_internal.h"
 
 namespace fesa {
 
@@ -228,10 +229,9 @@ double NormalizedMatrixError(const Matrix& lhs, const Matrix& rhs) {
 std::array<double, kGeneralizedComponentCount> GeneralizedStrain(
     const Matrix& b, const Vector& local_displacement) {
   std::array<double, kGeneralizedComponentCount> strain{};
+  const auto product = b.Multiply(local_displacement);
   for (std::size_t component = 0; component < strain.size(); ++component) {
-    for (std::size_t dof = 0; dof < local_displacement.Size(); ++dof) {
-      strain[component] += b(component, dof) * local_displacement[dof];
-    }
+    strain[component] = product[component];
   }
   return strain;
 }
@@ -395,13 +395,15 @@ Matrix EulerBeam3D::LocalStiffness() const {
   // Both Gauss points are required: a one-point bending rule loses two ranks.
   for (const double xi : gauss_points) {
     const Matrix b = StrainDisplacement(xi, length_);
+    Matrix weighted_b = b;
+    for (std::size_t component = 0; component < diagonal.size(); ++component)
+      for (std::size_t column = 0; column < kElementDofCount; ++column)
+        weighted_b(component, column) *= diagonal[component];
+    const auto contribution =
+        dense_blas_internal::MultiplyTransposeLeft(b, weighted_b);
     for (std::size_t row = 0; row < kElementDofCount; ++row) {
       for (std::size_t column = 0; column < kElementDofCount; ++column) {
-        for (std::size_t component = 0; component < diagonal.size();
-             ++component) {
-          stiffness(row, column) += b(component, row) * diagonal[component] *
-                                    b(component, column) * jacobian;
-        }
+        stiffness(row, column) += contribution(row, column) * jacobian;
       }
     }
   }
@@ -418,22 +420,17 @@ Matrix EulerBeam3D::GlobalStiffness() const {
   const Matrix local = LocalStiffness();
   const Matrix transform = Transformation(rotation_);
   const Matrix local_times_transform = local.Multiply(transform);
-  Matrix global{kElementDofCount, kElementDofCount};
   // Kg=T^T*Kl*T while dl=T*dg.
-  for (std::size_t row = 0; row < kElementDofCount; ++row) {
-    for (std::size_t column = 0; column < kElementDofCount; ++column) {
-      for (std::size_t inner = 0; inner < kElementDofCount; ++inner) {
-        global(row, column) +=
-            transform(inner, row) * local_times_transform(inner, column);
-      }
-    }
-  }
-  return global;
+  return dense_blas_internal::MultiplyTransposeLeft(transform,
+                                                    local_times_transform);
 }
 
 Vector EulerBeam3D::LocalEquivalentLoad(
     const ConstantLocalLineLoad& load) const {
   const std::array<double, 4> components = {load.px, load.py, load.pz, load.mx};
+  Matrix component_column{components.size(), 1U};
+  for (std::size_t component = 0; component < components.size(); ++component)
+    component_column(component, 0U) = components[component];
   Vector equivalent{kElementDofCount};
   const double inverse_sqrt_three = 1.0 / std::sqrt(3.0);
   const std::array<double, 2> gauss_points = {-inverse_sqrt_three,
@@ -441,12 +438,10 @@ Vector EulerBeam3D::LocalEquivalentLoad(
   const double jacobian = 0.5 * length_;
   for (const double xi : gauss_points) {
     const Matrix interpolation = KinematicInterpolation(xi, length_);
+    const auto contribution = dense_blas_internal::MultiplyTransposeLeft(
+        interpolation, component_column);
     for (std::size_t dof = 0; dof < equivalent.Size(); ++dof) {
-      for (std::size_t component = 0; component < components.size();
-           ++component) {
-        equivalent[dof] +=
-            interpolation(component, dof) * components[component] * jacobian;
-      }
+      equivalent[dof] += contribution(dof, 0U) * jacobian;
     }
   }
   return equivalent;
diff --git a/src/fesa/math/dense_blas_internal.cpp b/src/fesa/math/dense_blas_internal.cpp
index e19948c..6c13333 100644
--- a/src/fesa/math/dense_blas_internal.cpp
+++ b/src/fesa/math/dense_blas_internal.cpp
@@ -5,6 +5,17 @@
 
 namespace fesa::dense_blas_internal {
 
+Matrix MultiplyTransposeLeft(const Matrix& left, const Matrix& right) {
+  if (left.Rows() != right.Rows()) {
+    throw std::invalid_argument{"Transpose product dimensions differ."};
+  }
+  Matrix transposed{left.Columns(), left.Rows()};
+  for (std::size_t row = 0; row < left.Rows(); ++row)
+    for (std::size_t column = 0; column < left.Columns(); ++column)
+      transposed(column, row) = left(row, column);
+  return transposed.Multiply(right);
+}
+
 Result<MKL_INT> ToMklSize(const std::size_t size) {
   if (size > static_cast<std::size_t>((std::numeric_limits<MKL_INT>::max)())) {
     return Result<MKL_INT>::Failure(Status::Failure(
diff --git a/src/fesa/math/dense_blas_internal.h b/src/fesa/math/dense_blas_internal.h
index 056a20a..e96a05e 100644
--- a/src/fesa/math/dense_blas_internal.h
+++ b/src/fesa/math/dense_blas_internal.h
@@ -6,6 +6,7 @@
 #include <cstddef>
 
 #include "fesa/core/status.h"
+#include "fesa/math/matrix.h"
 
 namespace fesa::dense_blas_internal {
 
@@ -13,6 +14,10 @@ namespace fesa::dense_blas_internal {
 /// @return The converted length or an uncategorized overflow failure.
 Result<MKL_INT> ToMklSize(std::size_t size);
 
+/// @brief Computes left-transpose times right through ordinary dense MKL.
+/// @throws std::invalid_argument if the row counts differ.
+Matrix MultiplyTransposeLeft(const Matrix& left, const Matrix& right);
+
 /// @brief Copies contiguous double values through the private BLAS backend.
 /// @param source Readable storage for size values when size is nonzero.
 /// @param size Number of contiguous values to copy.
diff --git a/src/fesa/math/high_low_arithmetic.h b/src/fesa/math/high_low_arithmetic.h
index c724d1e..05ccad6 100644
--- a/src/fesa/math/high_low_arithmetic.h
+++ b/src/fesa/math/high_low_arithmetic.h
@@ -4,6 +4,8 @@
 #include <cmath>
 #include <cfenv>
 #include <cstdint>
+#include <cstring>
+#include <initializer_list>
 #include <limits>
 #if defined(_M_X64) || defined(__SSE2__)
 #include <xmmintrin.h>
@@ -56,9 +58,9 @@ inline bool IsNormalized(HighLow value) {
   if (!IsFinite(value)) return false;
   if (value.high == 0 && value.low == 0)
     return !std::signbit(value.high) && !std::signbit(value.low);
-  const auto normalized = TwoSum(value.high, value.low);
-  return IsFinite(normalized) && normalized.high == value.high &&
-         normalized.low == value.low;
+  // If RN(high + low) == high, TwoSum's bridge is zero and its error is
+  // exactly low. The finite and canonical exact-zero guards remain above.
+  return value.high + value.low == value.high;
 }
 /// @brief Rounds the full physical value to one finite output value.
 inline double Rounded(HighLow value) {
@@ -92,15 +94,27 @@ inline HighLow Subtract(HighLow left, HighLow right) {
 /// @brief Returns the exponent of the least nonzero binary64 significand bit.
 /// @pre value is finite and nonzero; gradual underflow is enabled.
 inline int LowestBitExponent(double value) {
-  int exponent = 0;
-  const double fraction = std::frexp(std::abs(value), &exponent);
-  auto significand = static_cast<std::uint64_t>(std::ldexp(fraction, 53));
+  static_assert(sizeof(double) == sizeof(std::uint64_t) &&
+                std::numeric_limits<double>::is_iec559 &&
+                std::numeric_limits<double>::digits == 53);
+  std::uint64_t bits = 0;
+  std::memcpy(&bits, &value, sizeof(bits));
+  const auto encoded_exponent = (bits >> 52U) & 0x7ffU;
+  auto significand = bits & 0xfffffffffffffULL;
+  const int exponent = encoded_exponent == 0U
+                           ? -1074
+                           : static_cast<int>(encoded_exponent) - 1075;
+  if (encoded_exponent != 0U) significand |= 1ULL << 52U;
   int trailing = 0;
-  while ((significand & 1U) == 0U) {
-    significand >>= 1U;
-    ++trailing;
+  // Strip large zero groups first without vendor/compiler intrinsics.
+  for (const int shift : {32, 16, 8, 4, 2, 1}) {
+    const auto mask = (1ULL << shift) - 1U;
+    if ((significand & mask) == 0U) {
+      significand >>= shift;
+      trailing += shift;
+    }
   }
-  return exponent - 53 + trailing;
+  return exponent + trailing;
 }
 /// @brief Gives a product and its FMA residual when the dyadic tail fits.
 /// @note Rejects underflow-risk exact tails, even if the rounded high is finite.
diff --git a/src/fesa/math/sparse_matrix.cpp b/src/fesa/math/sparse_matrix.cpp
index f0ff8c8..5c8bb75 100644
--- a/src/fesa/math/sparse_matrix.cpp
+++ b/src/fesa/math/sparse_matrix.cpp
@@ -8,11 +8,49 @@
 #include <tuple>
 #include <utility>
 
+#include <mkl.h>
+
 #include "fesa/fem/dof_manager.h"
+#include "math/dense_blas_internal.h"
 
 namespace fesa {
 namespace {
 
+/// @brief Checks the actual private backend integer width before conversion.
+MKL_INT SparseSize(const std::size_t size) {
+  const auto converted = dense_blas_internal::ToMklSize(size);
+  if (!converted.HasValue()) {
+    throw std::length_error{"Sparse size exceeds the MKL integer range."};
+  }
+  return converted.Value();
+}
+
+/// @brief Releases a call-local sparse handle before its borrowed arrays die.
+class SparseHandle {
+ public:
+  SparseHandle() = default;
+  SparseHandle(const SparseHandle&) = delete;
+  SparseHandle& operator=(const SparseHandle&) = delete;
+  ~SparseHandle() {
+    if (handle != nullptr) static_cast<void>(mkl_sparse_destroy(handle));
+  }
+  /// @brief Reports normal release failure before publishing the candidate.
+  sparse_status_t Destroy() {
+    const auto owned = handle;
+    handle = nullptr;
+    return mkl_sparse_destroy(owned);
+  }
+  sparse_matrix_t handle{nullptr};
+};
+
+/// @brief Fails explicitly instead of substituting a scalar backend.
+void CheckSparseStatus(const sparse_status_t status, const char* operation) {
+  if (status != SPARSE_STATUS_SUCCESS) {
+    throw std::runtime_error{std::string{"Sparse backend "} + operation +
+                             " failed."};
+  }
+}
+
 /// @brief Builds a structured sparse-matrix model failure.
 Status SparseFailure(const std::string& code, const std::string& identity,
                      const std::string& message) {
@@ -170,16 +208,43 @@ Vector SparseMatrix::Multiply(const Vector& rhs) const {
     throw std::invalid_argument{
         "Sparse matrix-vector multiplication has incompatible dimensions."};
   }
-
-  Vector result{rows_};
-  for (std::size_t row = 0U; row < rows_; ++row) {
-    double value = 0.0;
-    for (std::size_t position = row_offsets_[row];
-         position < row_offsets_[row + 1U]; ++position) {
-      value += values_[position] * rhs[column_indices_[position]];
-    }
-    result[row] = value;
+  if (!Validate().IsOk()) {
+    throw std::invalid_argument{"Sparse multiplication requires valid CSR."};
   }
+  const MKL_INT rows = SparseSize(rows_);
+  const MKL_INT columns = SparseSize(columns_);
+  static_cast<void>(SparseSize(values_.size()));
+  Vector result{rows_};
+  if (rows == 0 || columns == 0 || values_.empty()) return result;
+
+  // Each const call owns its converted arrays and handle. Neither copy/move
+  // nor concurrent calls retain pointers into another call's storage.
+  std::vector<MKL_INT> row_offsets;
+  std::vector<MKL_INT> column_indices;
+  row_offsets.reserve(row_offsets_.size());
+  column_indices.reserve(column_indices_.size());
+  for (const auto offset : row_offsets_) {
+    row_offsets.push_back(SparseSize(offset));
+  }
+  for (const auto column : column_indices_) {
+    column_indices.push_back(SparseSize(column));
+  }
+  auto values = values_;
+  SparseHandle sparse;
+  CheckSparseStatus(
+      mkl_sparse_d_create_csr(&sparse.handle, SPARSE_INDEX_BASE_ZERO, rows,
+                              columns, row_offsets.data(),
+                              row_offsets.data() + 1, column_indices.data(),
+                              values.data()),
+      "create CSR");
+  const matrix_descr descriptor{SPARSE_MATRIX_TYPE_GENERAL,
+                                SPARSE_FILL_MODE_FULL,
+                                SPARSE_DIAG_NON_UNIT};
+  CheckSparseStatus(
+      mkl_sparse_d_mv(SPARSE_OPERATION_NON_TRANSPOSE, 1.0, sparse.handle,
+                     descriptor, rhs.Data(), 0.0, result.Data()),
+      "multiply");
+  CheckSparseStatus(sparse.Destroy(), "destroy CSR");
   return result;
 }
 
diff --git a/src/fesa/results/global_equilibrium_recovery.cpp b/src/fesa/results/global_equilibrium_recovery.cpp
index 1fa096b..3f36bca 100644
--- a/src/fesa/results/global_equilibrium_recovery.cpp
+++ b/src/fesa/results/global_equilibrium_recovery.cpp
@@ -5,6 +5,7 @@
 #include <cmath>
 #include <cstddef>
 #include <limits>
+#include <stdexcept>
 #include <string>
 #include <vector>
 
@@ -26,33 +27,25 @@ double IndexedNorm(const Vector& values,
   return result;
 }
 
-/// @brief Computes separate Kff*df and Kfc*dc norms in stable CSR order.
+/// @brief Computes separate Kff*df and Kfc*dc norms through sparse MKL.
 std::array<double, 2> FreeEquationInternalTermNorms(
     const SparseMatrix& stiffness, const Vector& displacement,
     const DofManager& dofs) {
-  std::vector<unsigned char> free_column(dofs.FullDofCount(), 0U);
+  Vector free_displacement{displacement.Size()};
+  Vector constrained_displacement = displacement;
   for (const std::size_t full_dof : dofs.FreeDofs()) {
-    free_column[full_dof] = 1U;
+    free_displacement[full_dof] = displacement[full_dof];
+    constrained_displacement[full_dof] = 0.0;
   }
-
+  const auto free_term = stiffness.Multiply(free_displacement);
+  const auto constrained_term = stiffness.Multiply(constrained_displacement);
   double free_term_norm = 0.0;
   double constrained_term_norm = 0.0;
   for (const std::size_t row : dofs.FreeDofs()) {
-    double free_term = 0.0;
-    double constrained_term = 0.0;
-    for (std::size_t position = stiffness.RowOffsets()[row];
-         position < stiffness.RowOffsets()[row + 1U]; ++position) {
-      const std::size_t column = stiffness.ColumnIndices()[position];
-      const double contribution =
-          stiffness.Values()[position] * displacement[column];
-      if (free_column[column] != 0U) {
-        free_term += contribution;
-      } else {
-        constrained_term += contribution;
-      }
-    }
-    free_term_norm = std::hypot(free_term_norm, free_term);
-    constrained_term_norm = std::hypot(constrained_term_norm, constrained_term);
+    // Only the original free-row terms enter these acceptance scales.
+    free_term_norm = std::hypot(free_term_norm, free_term[row]);
+    constrained_term_norm =
+        std::hypot(constrained_term_norm, constrained_term[row]);
   }
   return {free_term_norm, constrained_term_norm};
 }
@@ -96,7 +89,7 @@ double NormalizedBalance(const std::array<double, 3>& balance,
 Result<GlobalEquilibriumResult> GlobalEquilibriumRecovery(
     const Domain& domain, const DofManager& dofs,
     const SparseMatrix& full_stiffness, const Vector& full_displacement,
-    const Vector& full_external_force) {
+    const Vector& full_external_force) try {
   const std::size_t full_count = dofs.FullDofCount();
   if (full_stiffness.Rows() != full_count ||
       full_stiffness.Columns() != full_count ||
@@ -171,6 +164,15 @@ Result<GlobalEquilibriumResult> GlobalEquilibriumRecovery(
   // shell global equilibrium below separately filters constrained reactions.
   result.reaction = result.residual;
   return Result<GlobalEquilibriumResult>::Success(std::move(result));
+} catch (const std::length_error&) {
+  return RecoveryResultFailure<GlobalEquilibriumResult>(
+      "recovery-index-overflow", {domain.SourcePath(), 0U},
+      domain.SourceContentIdentity(),
+      "Recovery multiplication exceeds the backend integer range.");
+} catch (const std::runtime_error&) {
+  return RecoveryResultFailure<GlobalEquilibriumResult>(
+      "recovery-multiplication-failed", {domain.SourcePath(), 0U},
+      domain.SourceContentIdentity(), "Recovery sparse multiplication failed.");
 }
 
 Status RecoverShellGlobalEvidence(const Domain& domain, const DofManager& dofs,
diff --git a/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp b/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp
index af2179f..1b42810 100644
--- a/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp
+++ b/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp
@@ -12,11 +12,13 @@
 #include <new>
 #include <numeric>
 #include <sstream>
+#include <stdexcept>
 #include <string>
 #include <utility>
 #include <vector>
 
 #include "fesa/math/sparse_matrix.h"
+#include "fesa/fem/dof_manager.h"
 #include "fesa/math/vector.h"
 #include "fesa/solvers/linear/mkl_pardiso_solver.h"
 #include "solvers/eigen/mkl_feast_eigen_solver_test_seam.h"
@@ -186,6 +188,37 @@ Result<Csr> CopyUpperOneBased(const SparseMatrix& matrix, const double scale) {
   return Result<Csr>::Success(std::move(copy));
 }
 
+/// @brief Keeps the independently scaled full or upper-mirrored operator.
+Result<SparseMatrix> ScaledActionMatrix(const SparseMatrix& matrix,
+                                       const double scale,
+                                       const bool mirror_upper = false) {
+  std::vector<CooContribution> contributions;
+  std::vector<std::vector<std::size_t>> columns(matrix.Rows());
+  for (std::size_t row = 0; row < matrix.Rows(); ++row) {
+    for (auto position = matrix.RowOffsets()[row];
+         position < matrix.RowOffsets()[row + 1U]; ++position) {
+      const auto column = matrix.ColumnIndices()[position];
+      if (mirror_upper && column < row) continue;
+      const auto value = matrix.Values()[position] / scale;
+      columns[row].push_back(column);
+      contributions.push_back({row, column, value, row, position});
+      if (mirror_upper && column != row) {
+        columns[column].push_back(row);
+        contributions.push_back({column, row, value, row, position});
+      }
+    }
+  }
+  SparsePattern pattern{{0U}, {}};
+  for (auto& row : columns) {
+    std::sort(row.begin(), row.end());
+    pattern.column_indices.insert(pattern.column_indices.end(), row.begin(),
+                                  row.end());
+    pattern.row_offsets.push_back(pattern.column_indices.size());
+  }
+  return SparseMatrix::FromCoo(matrix.Rows(), matrix.Columns(),
+                                std::move(contributions), pattern);
+}
+
 Result<Csr> IdentityOneBased(const std::size_t size) {
   if (!FitsMkl(size) || !FitsMkl(size + 1U)) {
     return Result<Csr>::Failure(
@@ -321,13 +354,11 @@ bool FillInitialBasis(std::vector<double>& vectors, const std::size_t n,
     }
     for (double& value : reflector) value /= reflector_norm;
     for (std::size_t target = column; target < columns; ++target) {
-      double dot = 0.0;
-      for (std::size_t row = column; row < n; ++row) {
-        dot += reflector[row - column] * vectors[row + target * n];
-      }
-      for (std::size_t row = column; row < n; ++row) {
-        vectors[row + target * n] -= 2.0 * reflector[row - column] * dot;
-      }
+      const auto length = static_cast<MKL_INT>(n - column);
+      auto* target_column = vectors.data() + column + target * n;
+      const double dot = cblas_ddot(length, reflector.data(), 1,
+                                    target_column, 1);
+      cblas_daxpy(length, -2.0 * dot, reflector.data(), 1, target_column, 1);
     }
     reflectors.push_back(std::move(reflector));
   }
@@ -338,13 +369,11 @@ bool FillInitialBasis(std::vector<double>& vectors, const std::size_t n,
     const std::size_t column = index - 1U;
     const auto& reflector = reflectors[column];
     for (std::size_t target = 0U; target < columns; ++target) {
-      double dot = 0.0;
-      for (std::size_t row = column; row < n; ++row) {
-        dot += reflector[row - column] * vectors[row + target * n];
-      }
-      for (std::size_t row = column; row < n; ++row) {
-        vectors[row + target * n] -= 2.0 * reflector[row - column] * dot;
-      }
+      const auto length = static_cast<MKL_INT>(n - column);
+      auto* target_column = vectors.data() + column + target * n;
+      const double dot = cblas_ddot(length, reflector.data(), 1,
+                                    target_column, 1);
+      cblas_daxpy(length, -2.0 * dot, reflector.data(), 1, target_column, 1);
     }
   }
   for (std::size_t column = 0U; column < columns; ++column)
@@ -353,10 +382,9 @@ bool FillInitialBasis(std::vector<double>& vectors, const std::size_t n,
   double gram_error = 0.0;
   for (std::size_t left = 0U; left < columns; ++left) {
     for (std::size_t right = 0U; right < columns; ++right) {
-      double dot = 0.0;
-      for (std::size_t row = 0U; row < n; ++row) {
-        dot += vectors[row + left * n] * vectors[row + right * n];
-      }
+      const double dot = cblas_ddot(static_cast<MKL_INT>(n),
+                                    vectors.data() + left * n, 1,
+                                    vectors.data() + right * n, 1);
       gram_error = std::hypot(gram_error, dot - (left == right ? 1.0 : 0.0));
     }
   }
@@ -390,7 +418,7 @@ Result<FeastEigenResult> MklFeastEigenSolver::Solve(
 Result<FeastEigenResult> MklFeastEigenSolver::SolveImpl(
     const SparseMatrix& whitened_operator, const SparseMatrix* mass,
     const std::size_t requested_modes, const double zero_operator_scale,
-    const double mass_lower_bound) const {
+    const double mass_lower_bound) const try {
   static_assert(sizeof(MKL_INT) == 8U,
                 "FEAST requires the installed ILP64 ABI.");
   const Status csr_status = whitened_operator.Validate();
@@ -472,6 +500,9 @@ Result<FeastEigenResult> MklFeastEigenSolver::SolveImpl(
   auto b =
       mass != nullptr ? CopyUpperOneBased(*mass, 1.0) : IdentityOneBased(n);
   if (!b.HasValue()) return Result<FeastEigenResult>::Failure(b.GetStatus());
+  auto action_matrix = ScaledActionMatrix(whitened_operator, operator_scale);
+  if (!action_matrix.HasValue())
+    return Result<FeastEigenResult>::Failure(action_matrix.GetStatus());
 
   double lower = -1.0;
   double upper = 1.0;
@@ -1045,13 +1076,11 @@ Result<FeastEigenResult> MklFeastEigenSolver::SolveImpl(
           "eigen-endpoint-separation", std::to_string(input),
           "A returned FEAST eigenvalue is too close to a solve endpoint."));
     }
-    double norm = 0.0;
+    Vector raw(n);
     for (std::size_t row = 0U; row < n; ++row)
-      norm += vectors[row + input * n] * vectors[row + input * n];
+      raw[row] = vectors[row + input * n];
+    double norm = raw.Dot(raw);
     if (mass != nullptr) {
-      Vector raw(n);
-      for (std::size_t row = 0U; row < n; ++row)
-        raw[row] = vectors[row + input * n];
       norm = raw.Dot(mass->Multiply(raw));
     }
     norm = std::sqrt(norm);
@@ -1060,21 +1089,10 @@ Result<FeastEigenResult> MklFeastEigenSolver::SolveImpl(
           TraceFailure("eigenpair-validation-failed", std::to_string(input),
                        "FEAST returned a non-normalizable eigenvector."));
     double residual_norm = 0.0;
+    const auto action = action_matrix.Value().Multiply(raw);
+    const auto mass_action = mass != nullptr ? mass->Multiply(raw) : raw;
     for (std::size_t row = 0U; row < n; ++row) {
-      double ax = 0.0;
-      for (std::size_t p = offsets[row]; p < offsets[row + 1U]; ++p) {
-        ax += values[p] / operator_scale *
-              vectors[whitened_operator.ColumnIndices()[p] + input * n];
-      }
-      double mx = vectors[row + input * n];
-      if (mass != nullptr) {
-        mx = 0.0;
-        for (std::size_t p = mass->RowOffsets()[row];
-             p < mass->RowOffsets()[row + 1U]; ++p)
-          mx +=
-              mass->Values()[p] * vectors[mass->ColumnIndices()[p] + input * n];
-      }
-      const double difference = ax - mu * mx;
+      const double difference = action[row] - mu * mass_action[row];
       residual_norm += difference * difference;
       ordered_vectors(row, output) = vectors[row + input * n] / norm;
     }
@@ -1109,20 +1127,11 @@ Result<FeastEigenResult> MklFeastEigenSolver::SolveImpl(
                 FormatDouble(rho_b) + ", mu=" + FormatDouble(mu) +
                 "; calls=" + CallHistory(calls)));
       const FeastBackendCall& last_call = calls.back();
-      std::vector<double> mirrored_action(n, 0.0);
-      std::vector<double> full_action(n, 0.0);
-      for (std::size_t row = 0U; row < n; ++row) {
-        for (std::size_t p = offsets[row]; p < offsets[row + 1U]; ++p) {
-          const std::size_t column = whitened_operator.ColumnIndices()[p];
-          const double entry = values[p] / operator_scale;
-          full_action[row] += entry * vectors[column + input * n];
-          if (column < row) continue;
-          mirrored_action[row] += entry * vectors[column + input * n];
-          if (column != row) {
-            mirrored_action[column] += entry * vectors[row + input * n];
-          }
-        }
-      }
+      auto mirrored_matrix =
+          ScaledActionMatrix(whitened_operator, operator_scale, true);
+      if (!mirrored_matrix.HasValue())
+        return Result<FeastEigenResult>::Failure(mirrored_matrix.GetStatus());
+      const auto mirrored_action = mirrored_matrix.Value().Multiply(raw);
       double mirrored_residual_squared = 0.0;
       double action_difference_squared = 0.0;
       for (std::size_t row = 0U; row < n; ++row) {
@@ -1130,7 +1139,7 @@ Result<FeastEigenResult> MklFeastEigenSolver::SolveImpl(
             mirrored_action[row] - mu * vectors[row + input * n];
         mirrored_residual_squared += mirrored_residual * mirrored_residual;
         const double action_difference =
-            full_action[row] - mirrored_action[row];
+            action[row] - mirrored_action[row];
         action_difference_squared += action_difference * action_difference;
       }
       const double rho_b_upper = std::sqrt(mirrored_residual_squared) / norm;
@@ -1174,22 +1183,20 @@ Result<FeastEigenResult> MklFeastEigenSolver::SolveImpl(
   }
   double gram_error = 0.0;
   double gram_squared = 0.0;
+  std::vector<Vector> mass_modes;
+  mass_modes.reserve(count);
+  for (std::size_t right = 0; right < count; ++right) {
+    Vector mode(n);
+    for (std::size_t row = 0; row < n; ++row)
+      mode[row] = ordered_vectors(row, right);
+    mass_modes.push_back(mass != nullptr ? mass->Multiply(mode)
+                                         : std::move(mode));
+  }
   for (std::size_t left = 0U; left < count; ++left) {
     for (std::size_t right = 0U; right < count; ++right) {
       double dot = 0.0;
       for (std::size_t row = 0U; row < n; ++row)
-        dot += ordered_vectors(row, left) * ordered_vectors(row, right);
-      if (mass != nullptr) {
-        dot = 0.0;
-        for (std::size_t row = 0U; row < n; ++row) {
-          double action = 0.0;
-          for (std::size_t p = mass->RowOffsets()[row];
-               p < mass->RowOffsets()[row + 1U]; ++p)
-            action += mass->Values()[p] *
-                      ordered_vectors(mass->ColumnIndices()[p], right);
-          dot += ordered_vectors(row, left) * action;
-        }
-      }
+        dot += ordered_vectors(row, left) * mass_modes[right][row];
       const double expected = left == right ? 1.0 : 0.0;
       const double difference = std::abs(dot - expected);
       gram_squared += difference * difference;
@@ -1233,6 +1240,14 @@ Result<FeastEigenResult> MklFeastEigenSolver::SolveImpl(
                           solve_upper,
                           upper_search_termination};
   return Result<FeastEigenResult>::Success(std::move(result));
+} catch (const std::length_error&) {
+  return Result<FeastEigenResult>::Failure(
+      Failure("feast-integer-overflow", "independent-action",
+              "Independent sparse action exceeds the backend range."));
+} catch (const std::runtime_error&) {
+  return Result<FeastEigenResult>::Failure(
+      Failure("eigen-multiplication-failed", "independent-action",
+              "Independent sparse action failed."));
 }
 
 }  // namespace fesa
diff --git a/tests/unit/math/dense_blas_internal_test.cpp b/tests/unit/math/dense_blas_internal_test.cpp
index d036800..c6cf25b 100644
--- a/tests/unit/math/dense_blas_internal_test.cpp
+++ b/tests/unit/math/dense_blas_internal_test.cpp
@@ -5,6 +5,7 @@
 #include <array>
 #include <cstddef>
 #include <limits>
+#include <stdexcept>
 
 namespace fesa::dense_blas_internal {
 namespace {
@@ -35,5 +36,26 @@ TEST(DenseBlasInternal, CDup003CopiesZeroAndNonzeroContiguousValues) {
   EXPECT_EQ(destination, source);
 }
 
+TEST(DenseBlasInternal, TransposeProductChecksShapeAndKnownRectangularAction) {
+  Matrix left{3U, 2U};
+  Matrix right{3U, 1U};
+  for (std::size_t row = 0; row < 3U; ++row) {
+    left(row, 0U) = static_cast<double>(row + 1U);
+    left(row, 1U) = 2.0;
+    right(row, 0U) = static_cast<double>(row);
+  }
+  const auto product = MultiplyTransposeLeft(left, right);
+  ASSERT_EQ(product.Rows(), 2U);
+  ASSERT_EQ(product.Columns(), 1U);
+  EXPECT_DOUBLE_EQ(product(0U, 0U), 8.0);
+  EXPECT_DOUBLE_EQ(product(1U, 0U), 6.0);
+  EXPECT_THROW(static_cast<void>(MultiplyTransposeLeft(left, Matrix{2U, 1U})),
+               std::invalid_argument);
+  const auto empty = MultiplyTransposeLeft(Matrix{0U, 2U}, Matrix{0U, 3U});
+  EXPECT_EQ(empty.Rows(), 2U);
+  EXPECT_EQ(empty.Columns(), 3U);
+  EXPECT_DOUBLE_EQ(empty(1U, 2U), 0.0);
+}
+
 }  // namespace
 }  // namespace fesa::dense_blas_internal
diff --git a/tests/unit/math/high_low_arithmetic_original.h b/tests/unit/math/high_low_arithmetic_original.h
new file mode 100644
index 0000000..072fe6b
--- /dev/null
+++ b/tests/unit/math/high_low_arithmetic_original.h
@@ -0,0 +1,162 @@
+#ifndef FESA_TEST_HIGH_LOW_ARITHMETIC_ORIGINAL_H_
+#define FESA_TEST_HIGH_LOW_ARITHMETIC_ORIGINAL_H_
+
+#include <cmath>
+#include <cfenv>
+#include <cstdint>
+#include <limits>
+#if defined(_M_X64) || defined(__SSE2__)
+#include <xmmintrin.h>
+#endif
+
+namespace fesa::high_low_original {
+/// @brief Checks the rounding and gradual-underflow premises without mutation.
+inline bool ArithmeticEnvironmentSupported() {
+  if (std::fegetround() != FE_TONEAREST) return false;
+#if defined(_M_X64) || defined(__SSE2__)
+  // MXCSR FTZ, DAZ and rounding-control bits must all be clear.
+  if ((_mm_getcsr() & 0xE040U) != 0U) return false;
+#endif
+  return true;
+}
+/// @brief Two binary64 components of one physical value; raw parts are kept.
+/// @note Arithmetic requires normalized finite parts and nearest rounding.
+/// This is an expansion with bounded rounding, not an IEEE quad type.
+struct HighLow {
+  double high{0};
+  double low{0};
+  HighLow() = default;
+  // Scalar lifting canonicalizes exact signed zero; raw pair inputs do not.
+  HighLow(double value) : high(value == 0 ? 0.0 : value) {}
+  HighLow(double high_part, double low_part)
+      : high(high_part), low(low_part) {}
+};
+
+/// @brief Gives an invalid arithmetic result for a checked caller to reject.
+inline HighLow InvalidHighLow() {
+  return {std::numeric_limits<double>::quiet_NaN(), 0};
+}
+/// @brief Checks both physical components, before comparisons or publication.
+inline bool IsFinite(HighLow value) {
+  return std::isfinite(value.high) && std::isfinite(value.low);
+}
+/// @brief Recovers the exact error of supported finite binary64 addition.
+inline HighLow TwoSum(double left, double right) {
+  const double sum = left + right;
+  const double bridge = sum - left;
+  const double restored_left = sum - bridge;
+  const double error = (left - restored_left) + (right - bridge);
+  if (!std::isfinite(sum) || !std::isfinite(bridge) ||
+      !std::isfinite(restored_left) || !std::isfinite(error))
+    return InvalidHighLow();
+  return sum == 0 && error == 0 ? HighLow{} : HighLow{sum, error};
+}
+/// @brief Requires the raw pair to equal its nearest-rounding normalization.
+inline bool IsNormalized(HighLow value) {
+  if (!IsFinite(value)) return false;
+  if (value.high == 0 && value.low == 0)
+    return !std::signbit(value.high) && !std::signbit(value.low);
+  const auto normalized = TwoSum(value.high, value.low);
+  return IsFinite(normalized) && normalized.high == value.high &&
+         normalized.low == value.low;
+}
+/// @brief Rounds the full physical value to one finite output value.
+inline double Rounded(HighLow value) {
+  return IsNormalized(value) ? value.high + value.low
+                             : std::numeric_limits<double>::quiet_NaN();
+}
+/// @brief Adds normalized expansions; rounding of discarded tail is bounded.
+inline HighLow Add(HighLow left, HighLow right) {
+  if (!IsNormalized(left) || !IsNormalized(right)) return InvalidHighLow();
+  const auto head = TwoSum(left.high, right.high);
+  const auto tail = TwoSum(left.low, right.low);
+  const auto middle = TwoSum(head.low, tail.high);
+  const auto joined = TwoSum(head.high, middle.high);
+  const double remainder = (joined.low + middle.low) + tail.low;
+  if (!IsFinite(head) || !IsFinite(tail) || !IsFinite(middle) ||
+      !IsFinite(joined) || !std::isfinite(remainder))
+    return InvalidHighLow();
+  return TwoSum(joined.high, remainder);
+}
+/// @brief Negates both components without collapsing the physical value.
+inline HighLow Negate(HighLow value) {
+  if (!IsNormalized(value)) return InvalidHighLow();
+  return value.high == 0 && value.low == 0
+             ? HighLow{}
+             : HighLow{-value.high, -value.low};
+}
+/// @brief Subtracts full expansions, including cancellations between highs.
+inline HighLow Subtract(HighLow left, HighLow right) {
+  return Add(left, Negate(right));
+}
+/// @brief Returns the exponent of the least nonzero binary64 significand bit.
+/// @pre value is finite and nonzero; gradual underflow is enabled.
+inline int LowestBitExponent(double value) {
+  int exponent = 0;
+  const double fraction = std::frexp(std::abs(value), &exponent);
+  auto significand = static_cast<std::uint64_t>(std::ldexp(fraction, 53));
+  int trailing = 0;
+  while ((significand & 1U) == 0U) {
+    significand >>= 1U;
+    ++trailing;
+  }
+  return exponent - 53 + trailing;
+}
+/// @brief Gives a product and its FMA residual when the dyadic tail fits.
+/// @note Rejects underflow-risk exact tails, even if the rounded high is finite.
+inline HighLow TwoProduct(double left, double right) {
+  if (!std::isfinite(left) || !std::isfinite(right)) return InvalidHighLow();
+  if (left == 0 || right == 0) return {};
+  if (LowestBitExponent(left) + LowestBitExponent(right) < -1074)
+    return InvalidHighLow();
+  const double product = left * right;
+  if (!std::isfinite(product)) return InvalidHighLow();
+  const double error = std::fma(left, right, -product);
+  if (!std::isfinite(error)) return InvalidHighLow();
+  return TwoSum(product, error);
+}
+/// @brief Retains all four component products in deterministic order.
+inline HighLow Multiply(HighLow left, HighLow right) {
+  if (!IsNormalized(left) || !IsNormalized(right)) return InvalidHighLow();
+  auto product = TwoProduct(left.high, right.high);
+  product = Add(product, TwoProduct(left.high, right.low));
+  product = Add(product, TwoProduct(left.low, right.high));
+  return Add(product, TwoProduct(left.low, right.low));
+}
+/// @brief Orders canonical values; 2 means unordered/unsupported input.
+inline int Compare(HighLow left, HighLow right) {
+  if (!IsNormalized(left) || !IsNormalized(right)) return 2;
+  if (left.high < right.high) return -1;
+  if (left.high > right.high) return 1;
+  if (left.low < right.low) return -1;
+  if (left.low > right.low) return 1;
+  return 0;
+}
+/// @brief Takes the absolute value using the full physical sign.
+inline HighLow Abs(HighLow value) {
+  const int sign = Compare(value, HighLow{});
+  return sign == 2 ? InvalidHighLow() : (sign < 0 ? Negate(value) : value);
+}
+/// @brief Refines a binary64 quotient twice using the full residual.
+/// @note Unsupported products/quotient underflow propagate an invalid result.
+inline HighLow Divide(HighLow numerator, HighLow denominator) {
+  if (!IsNormalized(numerator) || !IsNormalized(denominator) ||
+      Compare(denominator, HighLow{}) == 0)
+    return InvalidHighLow();
+  if (Compare(numerator, HighLow{}) == 0) return {};
+  const double first = numerator.high / denominator.high;
+  if (!std::isfinite(first) || first == 0) return InvalidHighLow();
+  auto quotient = HighLow{first};
+  for (int iteration = 0; iteration < 2; ++iteration) {
+    const auto residual = Subtract(numerator, Multiply(denominator, quotient));
+    if (!IsNormalized(residual)) return InvalidHighLow();
+    const double correction = residual.high / denominator.high;
+    if (!std::isfinite(correction) || (residual.high != 0 && correction == 0))
+      return InvalidHighLow();
+    quotient = Add(quotient, HighLow{correction});
+  }
+  return quotient;
+}
+}  // namespace fesa::high_low_original
+
+#endif  // FESA_TEST_HIGH_LOW_ARITHMETIC_ORIGINAL_H_
diff --git a/tests/unit/math/high_low_arithmetic_test.cpp b/tests/unit/math/high_low_arithmetic_test.cpp
index 092410d..03a94ec 100644
--- a/tests/unit/math/high_low_arithmetic_test.cpp
+++ b/tests/unit/math/high_low_arithmetic_test.cpp
@@ -3,7 +3,13 @@
 #include <gtest/gtest.h>
 
 #include <cmath>
+#include <cstdint>
+#include <cstring>
 #include <limits>
+#include <random>
+#include <vector>
+
+#include "high_low_arithmetic_original.h"
 
 namespace {
 using namespace fesa::math_internal;
@@ -125,4 +131,92 @@ TEST(HighLowArithmetic, RejectsChangedRoundingAndFlushToZeroEnvironment) {
   EXPECT_FALSE(ArithmeticEnvironmentSupported());
 #endif
 }
+
+std::uint64_t Bits(const double value) {
+  std::uint64_t result = 0;
+  std::memcpy(&result, &value, sizeof(result));
+  return result;
+}
+
+double FromBits(const std::uint64_t bits) {
+  double result = 0;
+  std::memcpy(&result, &bits, sizeof(result));
+  return result;
+}
+
+void ExpectSamePair(const HighLow actual,
+                    const fesa::high_low_original::HighLow expected) {
+  EXPECT_EQ(Bits(actual.high), Bits(expected.high));
+  EXPECT_EQ(Bits(actual.low), Bits(expected.low));
+}
+
+TEST(HighLowArithmetic, BitExponentMatchesFrozenOracleAcrossEveryExponent) {
+  ASSERT_TRUE(ArithmeticEnvironmentSupported());
+  // Every normal exponent, every possible trailing significand-bit count,
+  // both signs; then every subnormal single-bit and mixed-bit boundary.
+  for (std::uint64_t exponent = 1; exponent < 2047; ++exponent) {
+    for (unsigned trailing = 0; trailing <= 52U; ++trailing) {
+      const auto fraction = trailing == 52U ? 0ULL : 1ULL << trailing;
+      for (const auto sign : {0ULL, 1ULL << 63U}) {
+        const double value = FromBits(sign | (exponent << 52U) | fraction);
+        ASSERT_EQ(LowestBitExponent(value),
+                  fesa::high_low_original::LowestBitExponent(value));
+      }
+    }
+  }
+  for (unsigned bit = 0; bit < 52U; ++bit) {
+    for (const auto significand : {1ULL << bit, (1ULL << 52U) - (1ULL << bit)}) {
+      const double value = FromBits(significand);
+      EXPECT_EQ(LowestBitExponent(value),
+                fesa::high_low_original::LowestBitExponent(value));
+    }
+  }
+}
+
+TEST(HighLowArithmetic, OptimizedArithmeticMatchesFrozenOriginalBits) {
+  ASSERT_TRUE(ArithmeticEnvironmentSupported());
+  std::vector<double> values{0.0, -0.0,
+      std::numeric_limits<double>::denorm_min(),
+      -std::numeric_limits<double>::denorm_min(),
+      std::numeric_limits<double>::min(),
+      std::nextafter(std::numeric_limits<double>::min(), 0.0),
+      std::numeric_limits<double>::max(),
+      -std::numeric_limits<double>::max(), 1.0, -1.0,
+      std::nextafter(1.0, 2.0), std::ldexp(1.0, -53),
+      std::ldexp(1.0, -54), std::numeric_limits<double>::infinity(),
+      -std::numeric_limits<double>::infinity(),
+      std::numeric_limits<double>::quiet_NaN()};
+  std::mt19937_64 random{0x58a0fU};
+  for (int sample = 0; sample < 4000; ++sample) {
+    values.push_back(FromBits(random()));
+  }
+  for (std::size_t index = 0; index < values.size(); ++index) {
+    const double high = values[index];
+    const double low = values[(index * 13U + 7U) % values.size()];
+    const fesa::high_low_original::HighLow original{high, low};
+    EXPECT_EQ(IsNormalized({high, low}),
+              fesa::high_low_original::IsNormalized(original));
+    // Normalize a raw pair with the frozen operation, never the new predicate.
+    const auto a = fesa::high_low_original::TwoSum(high, low);
+    const auto b = fesa::high_low_original::TwoSum(
+        values[(index + 1U) % values.size()],
+        values[(index * 3U + 2U) % values.size()]);
+    ExpectSamePair(Add({a.high, a.low}, {b.high, b.low}),
+                   fesa::high_low_original::Add(a, b));
+    ExpectSamePair(Multiply({a.high, a.low}, {b.high, b.low}),
+                   fesa::high_low_original::Multiply(a, b));
+    ExpectSamePair(Divide({a.high, a.low}, {b.high, b.low}),
+                   fesa::high_low_original::Divide(a, b));
+    ExpectSamePair(TwoProduct(high, low),
+                   fesa::high_low_original::TwoProduct(high, low));
+  }
+  for (double high : {1.0, std::nextafter(1.0, 2.0), -1.0,
+                      std::numeric_limits<double>::max()}) {
+    for (double low : {0.0, -0.0, std::ldexp(1.0, -53),
+                       -std::ldexp(1.0, -53)}) {
+      EXPECT_EQ(IsNormalized({high, low}),
+                fesa::high_low_original::IsNormalized({high, low}));
+    }
+  }
+}
 }  // namespace
diff --git a/tests/unit/math/sparse_matrix_test.cpp b/tests/unit/math/sparse_matrix_test.cpp
index 7754a2f..525453e 100644
--- a/tests/unit/math/sparse_matrix_test.cpp
+++ b/tests/unit/math/sparse_matrix_test.cpp
@@ -3,6 +3,8 @@
 #include <gtest/gtest.h>
 
 #include <algorithm>
+#include <cmath>
+#include <future>
 #include <limits>
 #include <type_traits>
 #include <utility>
@@ -131,4 +133,101 @@ TEST(SparseAssembly, PreservesExpectedStructuralZeros) {
             4);
 }
 
+TEST(SparseAssembly, MultipliesRectangularCopyMoveAndConcurrentConstCalls) {
+  auto assembled = SparseMatrix::FromCoo(
+      3U, 5U, {{0U, 4U, 3.0, 0U, 0U}, {2U, 1U, -2.0, 1U, 0U}},
+      {{0U, 1U, 1U, 2U}, {4U, 1U}});
+  ASSERT_TRUE(assembled.HasValue());
+  const auto original = assembled.Value();
+  auto copied = original;
+  auto moved = std::move(copied);
+  Vector rhs{5U, 2.0};
+  const auto evaluate = [&rhs](const SparseMatrix& matrix) {
+    for (int repeat = 0; repeat < 20; ++repeat) {
+      const auto product = matrix.Multiply(rhs);
+      if (product.Size() != 3U || product[0U] != 6.0 ||
+          product[1U] != 0.0 || product[2U] != -4.0) return false;
+    }
+    return true;
+  };
+  auto first = std::async(std::launch::async, [&] { return evaluate(moved); });
+  auto second = std::async(std::launch::async, [&] { return evaluate(moved); });
+  EXPECT_TRUE(first.get());
+  EXPECT_TRUE(second.get());
+  EXPECT_TRUE(evaluate(original));
+  auto replacement = SparseMatrix::FromCoo(
+      1U, 1U, {{0U, 0U, 7.0, 0U, 0U}}, {{0U, 1U}, {0U}});
+  ASSERT_TRUE(replacement.HasValue());
+  moved = replacement.Value();
+  EXPECT_DOUBLE_EQ(moved.Multiply(Vector{1U, 3.0})[0U], 21.0);
+  EXPECT_TRUE(evaluate(original));
+}
+
+TEST(SparseAssembly, PreservesZeroShapesAndEmptyRows) {
+  auto no_rows = SparseMatrix::FromCoo(0U, 3U, {}, {{0U}, {}});
+  auto no_columns = SparseMatrix::FromCoo(2U, 0U, {}, {{0U, 0U, 0U}, {}});
+  auto no_entries = SparseMatrix::FromCoo(2U, 3U, {}, {{0U, 0U, 0U}, {}});
+  ASSERT_TRUE(no_rows.HasValue());
+  ASSERT_TRUE(no_columns.HasValue());
+  ASSERT_TRUE(no_entries.HasValue());
+  EXPECT_EQ(no_rows.Value().Multiply(Vector{3U}).Size(), 0U);
+  const auto empty_columns = no_columns.Value().Multiply(Vector{0U});
+  const auto empty_entries = no_entries.Value().Multiply(Vector{3U});
+  ASSERT_EQ(empty_columns.Size(), 2U);
+  ASSERT_EQ(empty_entries.Size(), 2U);
+  for (std::size_t row = 0; row < 2U; ++row) {
+    EXPECT_DOUBLE_EQ(empty_columns[row], 0.0);
+    EXPECT_DOUBLE_EQ(empty_entries[row], 0.0);
+  }
+}
+
+TEST(SparseAssembly, KeepsIeeeNonfiniteResultsAndRecoversAfterBadRhs) {
+  auto matrix = SparseMatrix::FromCoo(
+      2U, 2U, {{0U, 0U, 2.0, 0U, 0U}}, {{0U, 1U, 2U}, {0U, 1U}});
+  ASSERT_TRUE(matrix.HasValue());
+  Vector rhs{2U, (std::numeric_limits<double>::infinity)()};
+  const auto nonfinite = matrix.Value().Multiply(rhs);
+  EXPECT_TRUE(std::isinf(nonfinite[0U]));
+  EXPECT_TRUE(std::isnan(nonfinite[1U]));  // Stored structural zero times Inf.
+  rhs[0U] = (std::numeric_limits<double>::max)();
+  EXPECT_TRUE(std::isinf(matrix.Value().Multiply(rhs)[0U]));
+  EXPECT_THROW(static_cast<void>(matrix.Value().Multiply(Vector{1U})),
+               std::invalid_argument);
+  rhs[0U] = 3.0;
+  rhs[1U] = 4.0;
+  const auto finite = matrix.Value().Multiply(rhs);
+  EXPECT_DOUBLE_EQ(finite[0U], 6.0);
+  EXPECT_DOUBLE_EQ(finite[1U], 0.0);
+}
+
+TEST(SparseAssembly, MatchesIndependentLongDoubleResidualOnUnevenRows) {
+  const std::vector<double> entries{1.25, -0.375, 4.5, -3.0, 0.125, 2.0};
+  const SparsePattern pattern{{0U, 3U, 3U, 5U, 6U}, {0U, 2U, 4U, 1U, 3U, 2U}};
+  std::vector<CooContribution> contributions;
+  for (std::size_t row = 0; row < 4U; ++row) {
+    for (auto position = pattern.row_offsets[row];
+         position < pattern.row_offsets[row + 1U]; ++position) {
+      contributions.push_back(
+          {row, pattern.column_indices[position], entries[position], row,
+           position});
+    }
+  }
+  auto matrix = SparseMatrix::FromCoo(4U, 5U, contributions, pattern);
+  ASSERT_TRUE(matrix.HasValue());
+  Vector rhs{5U};
+  for (std::size_t column = 0; column < 5U; ++column) {
+    rhs[column] = 0.3 + static_cast<double>(column) * 0.7;
+  }
+  const auto product = matrix.Value().Multiply(rhs);
+  for (std::size_t row = 0; row < 4U; ++row) {
+    long double expected = 0;
+    for (auto position = pattern.row_offsets[row];
+         position < pattern.row_offsets[row + 1U]; ++position) {
+      expected += static_cast<long double>(entries[position]) *
+                  rhs[pattern.column_indices[position]];
+    }
+    EXPECT_NEAR(product[row], static_cast<double>(expected), 2.0e-14);
+  }
+}
+
 }  // namespace
diff --git a/tests/unit/results/result_recovery_components_test.cpp b/tests/unit/results/result_recovery_components_test.cpp
index 91bcc79..1f9a057 100644
--- a/tests/unit/results/result_recovery_components_test.cpp
+++ b/tests/unit/results/result_recovery_components_test.cpp
@@ -6,6 +6,7 @@
 #include <filesystem>
 #include <limits>
 #include <memory>
+#include <optional>
 #include <stdexcept>
 #include <string>
 #include <utility>
@@ -353,6 +354,41 @@ TEST(ResultRecoveryComponents,
   }
 }
 
+TEST(ResultRecoveryComponents, PartitionedScaleKeepsCancelledInternalTerms) {
+  auto definition = MakeShellDefinition();
+  definition.node_sets = {{"Fixed", std::optional<std::string>{"Shell-1"},
+                           {0U}, {definition.source_path, 60U}}};
+  definition.steps[0].boundaries = {
+      {"Fixed", 1, 1, 1.0, {definition.source_path, 61U}}};
+  const auto fixture = MakeShellFixture(std::move(definition));
+  const auto count = fixture.dofs->FullDofCount();
+  fesa::SparsePattern pattern;
+  pattern.row_offsets = {0U, 1U, 3U};
+  pattern.row_offsets.resize(count + 1U, 3U);
+  pattern.column_indices = {0U, 0U, 1U};
+  auto matrix = fesa::SparseMatrix::FromCoo(
+      count, count,
+      {{0U, 0U, 1.0, 0U, 0U}, {1U, 0U, -1.0, 0U, 1U},
+       {1U, 1U, 1.0, 0U, 2U}}, pattern);
+  ASSERT_TRUE(matrix.HasValue());
+  fesa::Vector displacement{count};
+  displacement[0U] = 1.0;
+  displacement[1U] = 1.0;
+  fesa::Vector external{count};
+  external[1U] = 1.0e-12;
+  const auto result = fesa::results_internal::GlobalEquilibriumRecovery(
+      *fixture.domain, *fixture.dofs, matrix.Value(), displacement, external);
+  ASSERT_TRUE(result.HasValue());
+  EXPECT_DOUBLE_EQ(result.Value().internal_force[1U], 0.0);
+  EXPECT_DOUBLE_EQ(result.Value().normalized_free_residual, 1.0e-12);
+  displacement[1U] = std::numeric_limits<double>::quiet_NaN();
+  const auto rejected = fesa::results_internal::GlobalEquilibriumRecovery(
+      *fixture.domain, *fixture.dofs, matrix.Value(), displacement, external);
+  ASSERT_FALSE(rejected.HasValue());
+  ExpectStatusCode(rejected.GetStatus(), "nonfinite-recovery-value");
+  EXPECT_DOUBLE_EQ(result.Value().reaction[1U], -1.0e-12);
+}
+
 // C-MODULE-002
 TEST(ResultRecoveryComponents,
      BeamResultRecoveryKeepsEndpointGaussAndStressIdentities) {
diff --git a/tests/unit/solvers/linear/tangent_solver_test.cpp b/tests/unit/solvers/linear/tangent_solver_test.cpp
index 0dab893..955d368 100644
--- a/tests/unit/solvers/linear/tangent_solver_test.cpp
+++ b/tests/unit/solvers/linear/tangent_solver_test.cpp
@@ -129,4 +129,44 @@ TEST(TangentSolver, RejectsSolveBeforeSuccessfulFactorization) {
   EXPECT_DOUBLE_EQ(correction[1U], -9.0);
 }
 
+TEST(TangentSolver, ReusesRhsButFreshValuesAndFailureCannotReuseStaleState) {
+  fesa::MklPardisoTangentSolver solver;
+  ASSERT_TRUE(solver.FactorizeFresh(MakeDenseCsr(2U, {2.0, 1.0, 0.0, 3.0}))
+                  .IsOk());
+  fesa::Vector correction{2U};
+  for (double scale : {1.0, -2.0, 0.5}) {
+    ASSERT_TRUE(solver.SolveCorrection(MakeVector({4.0 * scale, 6.0 * scale}),
+                                      correction).IsOk());
+    EXPECT_NEAR(correction[0U], scale, 1.0e-14);
+    EXPECT_NEAR(correction[1U], 2.0 * scale, 1.0e-14);
+  }
+  const auto bad_rhs = MakeVector({std::numeric_limits<double>::quiet_NaN(), 0});
+  correction = MakeVector({17, -4});
+  EXPECT_FALSE(solver.SolveCorrection(bad_rhs, correction).IsOk());
+  EXPECT_DOUBLE_EQ(correction[0U], 17);
+  EXPECT_DOUBLE_EQ(correction[1U], -4);
+  ASSERT_TRUE(solver.SolveCorrection(MakeVector({4, 6}), correction).IsOk());
+  ASSERT_TRUE(solver.FactorizeFresh(MakeDenseCsr(2U, {4, 0, 1, 2})).IsOk());
+  ASSERT_TRUE(solver.SolveCorrection(MakeVector({4, 5}), correction).IsOk());
+  EXPECT_NEAR(correction[0U], 1, 1.0e-14);
+  EXPECT_NEAR(correction[1U], 2, 1.0e-14);
+  auto diagonal = fesa::SparseMatrix::FromCoo(
+      2U, 2U, {{0U, 0U, 5, 0U, 0U}, {1U, 1U, 7, 1U, 0U}},
+      {{0U, 1U, 2U}, {0U, 1U}});
+  ASSERT_TRUE(diagonal.HasValue());
+  ASSERT_TRUE(solver.FactorizeFresh(diagonal.Value()).IsOk());
+  ASSERT_TRUE(solver.SolveCorrection(MakeVector({10, -7}), correction).IsOk());
+  EXPECT_NEAR(correction[0U], 2, 1.0e-14);
+  EXPECT_NEAR(correction[1U], -1, 1.0e-14);
+  EXPECT_FALSE(solver.FactorizeFresh(MakeDenseCsr(2U, {0, 0, 0, 0})).IsOk());
+  correction = MakeVector({19, -8});
+  const auto stale = solver.SolveCorrection(MakeVector({10, -7}), correction);
+  ASSERT_FALSE(stale.IsOk());
+  EXPECT_EQ(stale.Diagnostics().front().code, "solver-not-factorized");
+  EXPECT_DOUBLE_EQ(correction[0U], 19);
+  EXPECT_DOUBLE_EQ(correction[1U], -8);
+  ASSERT_TRUE(solver.FactorizeFresh(diagonal.Value()).IsOk());
+  EXPECT_TRUE(solver.SolveCorrection(MakeVector({10, -7}), correction).IsOk());
+}
+
 }  // namespace