fix(math): preserve certified sums and complete MKL algebra

src/fesa/analysis/eigen_reduction.cpp

@@ -122,12 +122,11 @@ Result<NodalLift> MakeLift(const Matrix& mass,
                   .Multiply(basis);
   double scale = 0.0;
   double basis_error = 0.0;
+  const auto basis_gram =
+      dense_blas_internal::MultiplyTransposeLeft(basis, basis);
   for (std::size_t i = 0; i < positive; ++i) {
     for (std::size_t j = 0; j < positive; ++j) {
-      double gram = 0.0;
-      for (std::size_t a = 0; a < 6U; ++a) {
-        gram += basis(a, i) * basis(a, j);
-      }
+      const double gram = basis_gram(i, j);
       basis_error = std::hypot(basis_error, gram - (i == j ? 1.0 : 0.0));
       scale = (std::max)(scale, std::abs(mp(i, j)));
     }
@@ -137,32 +136,34 @@ Result<NodalLift> MakeLift(const Matrix& mass,
     return Result<NodalLift>::Failure(Failure("invalid-eigen-mass-space", node,
                                               "Invalid positive mass basis."));
   }
+  Matrix normalized_mp(positive, positive);
+  for (std::size_t i = 0; i < positive; ++i)
+    for (std::size_t j = 0; j < positive; ++j) {
+      normalized_mp(i, j) = mp(i, j) / scale;
+      if (!std::isfinite(normalized_mp(i, j)))
+        return Result<NodalLift>::Failure(
+            Failure("eigen-mass-conditioning", node,
+                    "A positive-space mass pivot is not resolvable."));
+    }
   Matrix lower(positive, positive);
-  for (std::size_t i = 0; i < positive; ++i) {
-    for (std::size_t j = 0; j <= i; ++j) {
-      double value = mp(i, j) / scale;
-      for (std::size_t k = 0; k < j; ++k) value -= lower(i, k) * lower(j, k);
-      if (i == j) {
-        if (!(value > 0.0) || !std::isfinite(value)) {
-          return Result<NodalLift>::Failure(
-              Failure("eigen-mass-conditioning", node,
-                      "A positive-space mass pivot is not resolvable."));
-        }
-        lower(i, j) = std::sqrt(value);
-      } else {
-        lower(i, j) = value / lower(j, j);
-      }
-    }
-  }
   Matrix inverse_lower(positive, positive);
-  for (std::size_t j = 0; j < positive; ++j) {
-    for (std::size_t i = 0; i < positive; ++i) {
-      double value = i == j ? 1.0 : 0.0;
-      for (std::size_t k = 0; k < i; ++k)
-        value -= lower(i, k) * inverse_lower(k, j);
-      inverse_lower(i, j) = value / lower(i, i);
-    }
+  if (!dense_blas_internal::CholeskyLowerAndInverse(normalized_mp, lower,
+                                                  inverse_lower)) {
+    return Result<NodalLift>::Failure(
+        Failure("eigen-mass-conditioning", node,
+                "A positive-space mass pivot is not resolvable."));
   }
+  Matrix transposed_lower(positive, positive);
+  Matrix scaled_inverse_transpose(positive, positive);
+  const double root_scale = std::sqrt(scale);
+  for (std::size_t i = 0; i < positive; ++i)
+    for (std::size_t j = 0; j < positive; ++j) {
+      transposed_lower(i, j) = lower(j, i);
+      scaled_inverse_transpose(i, j) = inverse_lower(j, i) / root_scale;
+    }
+  const auto inverse_mass = dense_blas_internal::MultiplyTransposeLeft(
+      inverse_lower, inverse_lower);
+  const auto rebuilt_mass = lower.Multiply(transposed_lower);
   double norm_mp = 0.0;
   double norm_inverse = 0.0;
   double reconstruction = 0.0;
@@ -172,12 +173,8 @@ Result<NodalLift> MakeLift(const Matrix& mass,
     double row_mp = 0.0;
     double row_inverse = 0.0;
     for (std::size_t j = 0; j < positive; ++j) {
-      double inverse = 0.0;
-      double rebuilt = 0.0;
-      for (std::size_t k = 0; k < positive; ++k) {
-        inverse += inverse_lower(k, i) * inverse_lower(k, j);
-        rebuilt += lower(i, k) * lower(j, k);
-      }
+      const double inverse = inverse_mass(i, j);
+      const double rebuilt = rebuilt_mass(i, j);
       row_inverse += std::abs(inverse);
       row_mp += std::abs(mp(i, j) / scale);
       reconstruction = std::hypot(reconstruction, rebuilt - mp(i, j) / scale);
@@ -195,13 +192,9 @@ Result<NodalLift> MakeLift(const Matrix& mass,
         "eigen-mass-conditioning", node,
         "Nodal mass fails rcond_inf >= 2^-32 or Cholesky reconstruction."));
   }
-  const double root_scale = std::sqrt(scale);
+  result.scaled_lift = basis.Multiply(scaled_inverse_transpose);
   for (std::size_t c = 0; c < 6U; ++c) {
     for (std::size_t j = 0; j < positive; ++j) {
-      for (std::size_t k = 0; k < positive; ++k) {
-        result.scaled_lift(c, j) +=
-            basis(c, k) * inverse_lower(j, k) / root_scale;
-      }
       if (!std::isfinite(result.scaled_lift(c, j))) {
         return Result<NodalLift>::Failure(Failure(
             "eigen-mass-conditioning", node, "Nodal whitening is nonfinite."));
@@ -280,19 +273,23 @@ Result<EigenReduction> BuildEigenReduction(
     free_count += lift.Value().record.free_count;
     null_count += lift.Value().record.null_count;
     const auto& z = lift.Value().null_vector;
+    Vector null_coordinate{6U};
+    Matrix null_column{6U, 1U};
+    for (std::size_t a = 0; a < 6U; ++a) {
+      null_coordinate[a] = z[a];
+      null_column(a, 0U) = z[a];
+    }
+    const auto mass_null_action = scaled.Multiply(null_coordinate);
+    const auto positive_null = dense_blas_internal::MultiplyTransposeLeft(
+        lift.Value().positive_basis, null_column);
     for (std::size_t a = 0; a < 6U; ++a) {
       null_norm = std::hypot(null_norm, z[a]);
       expected_kz_norm = std::hypot(expected_kz_norm, drilling[node] * z[a]);
       if (mask[6U * node + a] != 0U) continue;
-      double action = 0.0;
-      for (std::size_t b = 0; b < 6U; ++b) action += scaled(a, b) * z[b];
-      mass_null_error = std::hypot(mass_null_error, action);
+      mass_null_error = std::hypot(mass_null_error, mass_null_action[a]);
     }
     for (std::size_t j = 0; j < lift.Value().positive_basis.Columns(); ++j) {
-      double dot = 0.0;
-      for (std::size_t a = 0; a < 6U; ++a)
-        dot += lift.Value().positive_basis(a, j) * z[a];
-      pz_norm = std::hypot(pz_norm, dot);
+      pz_norm = std::hypot(pz_norm, positive_null(j, 0U));
     }
     lifts.push_back(std::move(lift.Value()));
   }
@@ -318,10 +315,6 @@ Result<EigenReduction> BuildEigenReduction(
     for (std::size_t a = 0; a < 6U; ++a) {
       const auto row = node * 6U + a;
       if (mask[row] != 0U) continue;
-      std::map<std::size_t, double> null_action;
-      if (lifts[node].record.null_count != 0U) {
-        null_action[node] = -drilling[node] * lifts[node].null_vector[a];
-      }
       for (std::size_t pos = offsets[row]; pos < offsets[row + 1U]; ++pos) {
         const auto col = columns[pos];
         if (mask[col] != 0U) continue;
@@ -334,14 +327,28 @@ Result<EigenReduction> BuildEigenReduction(
                       "Scaled stiffness is nonfinite."));
         }
         stiffness_norm = std::hypot(stiffness_norm, value);
-        null_action[other] += value * lifts[other].null_vector[b];
         auto block = blocks.try_emplace(other, 6U, 6U).first;
         block->second(a, b) = value;
       }
-      for (const auto& item : null_action)
-        kz_error = std::hypot(kz_error, item.second);
     }
     const auto& left = lifts[node];
+    std::map<std::size_t, Vector> null_actions;
+    for (const auto& entry : blocks) {
+      Vector null_coordinate{6U};
+      for (std::size_t b = 0; b < 6U; ++b)
+        null_coordinate[b] = lifts[entry.first].null_vector[b];
+      null_actions.emplace(entry.first, entry.second.Multiply(null_coordinate));
+    }
+    if (left.record.null_count != 0U) {
+      auto& local = null_actions.try_emplace(node, 6U).first->second;
+      for (std::size_t a = 0; a < 6U; ++a)
+        local[a] -= drilling[node] * left.null_vector[a];
+    }
+    for (std::size_t a = 0; a < 6U; ++a) {
+      if (mask[node * 6U + a] != 0U) continue;
+      for (const auto& action : null_actions)
+        kz_error = std::hypot(kz_error, action.second[a]);
+    }
     std::map<std::size_t, Matrix> projected_blocks;
     for (const auto& entry : blocks) {
       projected_blocks.emplace(

src/fesa/analysis/thermal_transient_algebra.cpp

@@ -8,6 +8,7 @@
 #include <utility>
 
 #include "thermal_transient_bounds.h"
+#include "math/dense_blas_internal.h"
 
 namespace fesa::thermal_transient_internal {
 namespace {
@@ -117,14 +118,9 @@ 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];
-  }
+  const auto entries = s.c.Values().size();
+  Vector capacity_diagonal{entries}, conductivity_diagonal{entries};
+  Vector capacity_input{entries}, conductivity_input{entries};
   Vector b{n};
   std::vector<I> exact_a, exact_b(n);
   for (std::size_t f = 0; f < n; ++f) {
@@ -135,43 +131,39 @@ 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});
+      capacity_diagonal[j] = c / delta;
       if (free[col] != n) {
+        capacity_input[j] = old[col];
         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}));
         exact_b[f] = Add(exact_b[f], Mul(cd, I{old[col]}));
       } else {
         const double g = p.prescribed[constrained[col]];
+        capacity_input[j] = old[col] - g;
+        conductivity_diagonal[j] = k;
+        conductivity_input[j] = g;
         // Preserve raw old constrained capacity; never project frame zero.
         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];
+  const auto capacity_action = dense_blas_internal::MultiplyDiagonal(
+      capacity_diagonal, capacity_input);
+  const auto conductivity_action = dense_blas_internal::MultiplyDiagonal(
+      conductivity_diagonal, conductivity_input);
+  for (std::size_t f = 0; f < n; ++f) {
+    const auto row = p.free[f];
+    for (auto j = s.c.RowOffsets()[row]; j < s.c.RowOffsets()[row + 1]; ++j) {
+      const auto col = s.c.ColumnIndices()[j];
+      // Preserve the original per-entry subtraction before the row sum.
+      // Separately summing C*delta and K*g loses small stable RHS terms.
+      b[f] += free[col] != n ? capacity_action[j]
+                             : capacity_action[j] - conductivity_action[j];
+    }
   }
-  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))
@@ -493,6 +485,9 @@ 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::invalid_argument&) {
+  return Result<ThermalStepCandidate>::Failure(
+      Failure("Nonfinite backward Euler RHS."));
 } catch (const std::length_error&) {
   return Result<ThermalStepCandidate>::Failure(
       Failure("Thermal sparse action exceeds the backend integer range."));

src/fesa/math/dense_blas_internal.cpp

@@ -1,10 +1,72 @@
 #include "math/dense_blas_internal.h"
 
+#include <cmath>
 #include <limits>
 #include <stdexcept>
+#include <utility>
+#include <vector>
 
 namespace fesa::dense_blas_internal {
 
+bool CholeskyLowerAndInverse(const Matrix& matrix, Matrix& lower,
+                            Matrix& inverse_lower) {
+  if (matrix.Rows() != matrix.Columns())
+    throw std::invalid_argument{"Cholesky requires a square matrix."};
+  const auto size = ToMklSize(matrix.Rows());
+  if (!size.HasValue())
+    throw std::length_error{"Cholesky size exceeds the MKL integer range."};
+  Matrix candidate_lower{matrix.Rows(), matrix.Columns()};
+  Matrix candidate_inverse{matrix.Rows(), matrix.Columns()};
+  if (matrix.Rows() == 0U) {
+    lower = std::move(candidate_lower);
+    inverse_lower = std::move(candidate_inverse);
+    return true;
+  }
+  std::vector<double> values(matrix.Rows() * matrix.Columns());
+  for (std::size_t row = 0; row < matrix.Rows(); ++row)
+    for (std::size_t column = 0; column < matrix.Columns(); ++column) {
+      if (!std::isfinite(matrix(row, column)))
+        throw std::invalid_argument{"Cholesky inputs must be finite."};
+      values[row * matrix.Columns() + column] = matrix(row, column);
+    }
+  if (LAPACKE_dpotrf(LAPACK_ROW_MAJOR, 'L', size.Value(), values.data(),
+                     size.Value()) != 0) return false;
+  for (std::size_t row = 0; row < matrix.Rows(); ++row)
+    for (std::size_t column = 0; column <= row; ++column)
+      candidate_lower(row, column) = values[row * matrix.Columns() + column];
+  if (LAPACKE_dtrtri(LAPACK_ROW_MAJOR, 'L', 'N', size.Value(), values.data(),
+                     size.Value()) != 0) return false;
+  for (std::size_t row = 0; row < matrix.Rows(); ++row)
+    for (std::size_t column = 0; column <= row; ++column) {
+      const auto value = values[row * matrix.Columns() + column];
+      if (!std::isfinite(candidate_lower(row, column)) || !std::isfinite(value))
+        return false;
+      candidate_inverse(row, column) = value;
+    }
+  lower = std::move(candidate_lower);
+  inverse_lower = std::move(candidate_inverse);
+  return true;
+}
+
+Vector MultiplyDiagonal(const Vector& diagonal, const Vector& rhs) {
+  if (diagonal.Size() != rhs.Size())
+    throw std::invalid_argument{"Diagonal product dimensions differ."};
+  const auto size = ToMklSize(diagonal.Size());
+  if (!size.HasValue())
+    throw std::length_error{"Diagonal size exceeds the MKL integer range."};
+  for (std::size_t i = 0; i < rhs.Size(); ++i) {
+    if (!std::isfinite(diagonal[i]) || !std::isfinite(rhs[i]))
+      throw std::invalid_argument{"Diagonal product inputs must be finite."};
+  }
+  Vector candidate{rhs.Size()};
+  if (rhs.Size() == 0U) return candidate;
+  // k=0 stores exactly one coefficient per independent row. MKL performs
+  // each product without regrouping a caller's subsequent certified sums.
+  cblas_dsbmv(CblasColMajor, CblasUpper, size.Value(), 0, 1.0, diagonal.Data(),
+              1, rhs.Data(), 1, 0.0, candidate.Data(), 1);
+  return candidate;
+}
+
 Matrix MultiplyTransposeLeft(const Matrix& left, const Matrix& right) {
   if (left.Rows() != right.Rows()) {
     throw std::invalid_argument{"Transpose product dimensions differ."};

src/fesa/math/dense_blas_internal.h

@@ -7,6 +7,7 @@
 
 #include "fesa/core/status.h"
 #include "fesa/math/matrix.h"
+#include "fesa/math/vector.h"
 
 namespace fesa::dense_blas_internal {
 
@@ -18,6 +19,19 @@ Result<MKL_INT> ToMklSize(std::size_t size);
 /// @throws std::invalid_argument if the row counts differ.
 Matrix MultiplyTransposeLeft(const Matrix& left, const Matrix& right);
 
+/// @brief Multiplies a compact diagonal matrix by a vector through MKL.
+/// @throws std::invalid_argument if shapes differ or inputs are nonfinite.
+/// @throws std::length_error if the size exceeds the backend integer range.
+Vector MultiplyDiagonal(const Vector& diagonal, const Vector& rhs);
+
+/// @brief Computes a lower Cholesky factor and its inverse with MKL LAPACK.
+/// @return False for a nonpositive pivot or nonfinite factor/inverse.
+/// @note Outputs are replaced only after both finite candidates succeed.
+/// @throws std::invalid_argument for a nonsquare or nonfinite input.
+/// @throws std::length_error for a backend size overflow.
+bool CholeskyLowerAndInverse(const Matrix& matrix, Matrix& lower,
+                            Matrix& inverse_lower);
+
 /// @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.

src/fesa/math/sparse_matrix.cpp

@@ -211,6 +211,10 @@ Vector SparseMatrix::Multiply(const Vector& rhs) const {
   if (!Validate().IsOk()) {
     throw std::invalid_argument{"Sparse multiplication requires valid CSR."};
   }
+  for (std::size_t i = 0; i < rhs.Size(); ++i) {
+    if (!std::isfinite(rhs[i]))
+      throw std::invalid_argument{"Sparse multiplication RHS must be finite."};
+  }
   const MKL_INT rows = SparseSize(rows_);
   const MKL_INT columns = SparseSize(columns_);
   static_cast<void>(SparseSize(values_.size()));

src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp

@@ -22,6 +22,7 @@
 #include "fesa/math/vector.h"
 #include "fesa/solvers/linear/mkl_pardiso_solver.h"
 #include "solvers/eigen/mkl_feast_eigen_solver_test_seam.h"
+#include "math/dense_blas_internal.h"
 
 namespace fesa {
 
@@ -1095,6 +1096,10 @@ Result<FeastEigenResult> MklFeastEigenSolver::SolveImpl(
       const double difference = action[row] - mu * mass_action[row];
       residual_norm += difference * difference;
       ordered_vectors(row, output) = vectors[row + input * n] / norm;
+      if (!std::isfinite(ordered_vectors(row, output)))
+        return Result<FeastEigenResult>::Failure(
+            TraceFailure("eigenpair-validation-failed", std::to_string(input),
+                         "FEAST returned a non-normalizable eigenvector."));
     }
     double residual_denominator = norm;
     if (mass != nullptr) {
@@ -1183,20 +1188,20 @@ Result<FeastEigenResult> MklFeastEigenSolver::SolveImpl(
   }
   double gram_error = 0.0;
   double gram_squared = 0.0;
-  std::vector<Vector> mass_modes;
-  mass_modes.reserve(count);
+  Matrix mass_modes{n, 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));
+    const auto action = mass != nullptr ? mass->Multiply(mode) : std::move(mode);
+    for (std::size_t row = 0; row < n; ++row)
+      mass_modes(row, right) = action[row];
   }
+  const auto gram = dense_blas_internal::MultiplyTransposeLeft(ordered_vectors,
+                                                              mass_modes);
   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) * mass_modes[right][row];
+      const double dot = gram(left, right);
       const double expected = left == right ? 1.0 : 0.0;
       const double difference = std::abs(dot - expected);
       gram_squared += difference * difference;