fix(math): preflight finite dense operands at dispatch boundary

include/fesa/math/matrix.h

@@ -41,11 +41,11 @@ class Matrix {
   const double& operator()(std::size_t row, std::size_t column) const;
 
   /// @brief Multiplies this row-major matrix by a dense vector.
-  /// @throws std::invalid_argument if the dimensions are incompatible.
+  /// @throws std::invalid_argument for incompatible shapes or nonfinite input.
   Vector Multiply(const Vector& rhs) const;
 
   /// @brief Multiplies this row-major matrix by another dense matrix.
-  /// @throws std::invalid_argument if the dimensions are incompatible.
+  /// @throws std::invalid_argument for incompatible shapes or nonfinite input.
   Matrix Multiply(const Matrix& rhs) const;
 
  private:

src/fesa/analysis/eigen_reduction.cpp

@@ -5,6 +5,7 @@
 #include <cmath>
 #include <limits>
 #include <map>
+#include <stdexcept>
 #include <string>
 #include <utility>
 #include <vector>
@@ -226,7 +227,7 @@ Result<EigenReduction> BuildEigenReduction(
     const SparseMatrix& stiffness, const std::vector<Matrix>& mass,
     const std::vector<ShellNodeInitialFrame>& frames,
     const std::vector<std::uint8_t>& mask, const std::vector<double>& drilling,
-    const double length) {
+    const double length) try {
   const std::size_t count = mass.size();
   const Status valid = stiffness.Validate();
   if (!valid.IsOk()) return Result<EigenReduction>::Failure(valid);
@@ -401,6 +402,10 @@ Result<EigenReduction> BuildEigenReduction(
   return Result<EigenReduction>::Success(
       {std::move(reduced.Value()), std::move(lifts), free_count, null_count,
        length, eta_mz, eta_kz, pz_norm});
+} catch (const std::invalid_argument&) {
+  return Result<EigenReduction>::Failure(
+      Failure("invalid-eigen-mass-space", 0U,
+              "Eigen product operands must be finite and shape compatible."));
 }
 
 Vector RestoreScaledMode(const EigenReduction& reduction, const Matrix& vectors,

src/fesa/analysis/linear_frequency_procedure.cpp

@@ -275,6 +275,10 @@ Result<EigenAnalysisState> ExecuteCandidate(const Domain& domain) {
     std::size_t pivot = 0U;
     for (std::size_t dof = 0; dof < y.Size(); ++dof) {
       physical[dof] = D(dof, length) * y[dof];
+      if (!std::isfinite(physical[dof]))
+        return Result<EigenAnalysisState>::Failure(
+            MetricFailure("initial_modal_mass", physical[dof],
+                          state.backend.column_ids[column]));
       if (std::abs(y[dof]) > std::abs(y[pivot])) pivot = dof;
     }
     Vector mass_times = eigen_internal::MultiplyNodalMass(nodal_mass, physical);
@@ -413,6 +417,10 @@ Result<EigenAnalysisState> LinearFrequencyProcedureExecutor::Execute(
           std::move(diagnostics)));
     }
     return result;
+  } catch (const std::invalid_argument&) {
+    return Result<EigenAnalysisState>::Failure(
+        Failure("invalid-eigenvalue-result", "modal-workspace",
+                "Modal product inputs must be finite and shape compatible."));
   } catch (const std::bad_alloc&) {
     return Result<EigenAnalysisState>::Failure(
         Failure("eigen-allocation-failed", "modal-workspace",

src/fesa/math/dense_blas_internal.h

@@ -31,7 +31,7 @@ class ScopedSingleEquationThreads {
 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.
+/// @throws std::invalid_argument for mismatched rows or nonfinite operands.
 Matrix MultiplyTransposeLeft(const Matrix& left, const Matrix& right);
 
 /// @brief Multiplies a compact diagonal matrix by a vector through MKL.

src/fesa/math/matrix.cpp

@@ -1,5 +1,6 @@
 #include "fesa/math/matrix.h"
 
+#include <cmath>
 #include <mkl.h>
 
 #include <limits>
@@ -107,6 +108,12 @@ Vector Matrix::Multiply(const Vector& rhs) const {
     throw std::invalid_argument{
         "Matrix-vector multiplication has incompatible dimensions."};
   }
+  for (const double value : values_)
+    if (!std::isfinite(value))
+      throw std::invalid_argument{"Matrix product inputs must be finite."};
+  for (std::size_t i = 0; i < rhs.Size(); ++i)
+    if (!std::isfinite(rhs.Data()[i]))
+      throw std::invalid_argument{"Matrix product inputs must be finite."};
 
   Vector result{rows_};
   if (rows_ == 0 || columns_ == 0) {
@@ -126,6 +133,12 @@ Matrix Matrix::Multiply(const Matrix& rhs) const {
     throw std::invalid_argument{
         "Matrix multiplication has incompatible dimensions."};
   }
+  for (const double value : values_)
+    if (!std::isfinite(value))
+      throw std::invalid_argument{"Matrix product inputs must be finite."};
+  for (const double value : rhs.values_)
+    if (!std::isfinite(value))
+      throw std::invalid_argument{"Matrix product inputs must be finite."};
 
   Matrix result{rows_, rhs.columns_};
   if (rows_ == 0 || columns_ == 0 || rhs.columns_ == 0) {

src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp

@@ -1194,9 +1194,14 @@ Result<FeastEigenResult> MklFeastEigenSolver::SolveImpl(
     for (std::size_t row = 0; row < n; ++row)
       mode[row] = ordered_vectors(row, right);
     const auto action = mass != nullptr ? mass->Multiply(mode) : std::move(mode);
-    for (std::size_t row = 0; row < n; ++row)
+    for (std::size_t row = 0; row < n; ++row) {
+      if (!std::isfinite(action[row]))
+        return Result<FeastEigenResult>::Failure(
+            TraceFailure("eigenpair-validation-failed", "gram",
+                         "The independent mass action is nonfinite."));
       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) {
@@ -1245,6 +1250,10 @@ Result<FeastEigenResult> MklFeastEigenSolver::SolveImpl(
                           solve_upper,
                           upper_search_termination};
   return Result<FeastEigenResult>::Success(std::move(result));
+} catch (const std::invalid_argument&) {
+  return Result<FeastEigenResult>::Failure(
+      Failure("eigenpair-validation-failed", "independent-action",
+              "Modal product inputs must be finite and shape compatible."));
 } catch (const std::length_error&) {
   return Result<FeastEigenResult>::Failure(
       Failure("feast-integer-overflow", "independent-action",

tests/unit/math/dense_blas_internal_test.cpp

@@ -56,6 +56,13 @@ TEST(DenseBlasInternal, TransposeProductChecksShapeAndKnownRectangularAction) {
   EXPECT_EQ(empty.Rows(), 2U);
   EXPECT_EQ(empty.Columns(), 3U);
   EXPECT_DOUBLE_EQ(empty(1U, 2U), 0.0);
+  left(0U, 0U) = std::numeric_limits<double>::quiet_NaN();
+  EXPECT_THROW(static_cast<void>(MultiplyTransposeLeft(left, right)),
+               std::invalid_argument);
+  left(0U, 0U) = 1.0;
+  right(0U, 0U) = std::numeric_limits<double>::infinity();
+  EXPECT_THROW(static_cast<void>(MultiplyTransposeLeft(left, right)),
+               std::invalid_argument);
 }
 
 TEST(DenseBlasInternal, DiagonalProductKeepsEntrywiseCancellationAndChecks) {

tests/unit/math/matrix_test.cpp

@@ -109,5 +109,23 @@ TEST(DenseMath, RowMajorMatrixMatchesKnownGemvGemm) {
                std::invalid_argument);
 }
 
+TEST(DenseMath, MatrixProductsRejectNonfiniteOperandsAndPreserveValidReuse) {
+  Matrix left{2U, 2U, 1.0};
+  Matrix right{2U, 1U, 2.0};
+  Vector rhs{2U, 2.0};
+  rhs[1U] = std::numeric_limits<double>::infinity();
+  EXPECT_THROW(static_cast<void>(left.Multiply(rhs)), std::invalid_argument);
+  rhs[1U] = 2.0;
+  left(0U, 1U) = std::numeric_limits<double>::quiet_NaN();
+  EXPECT_THROW(static_cast<void>(left.Multiply(rhs)), std::invalid_argument);
+  EXPECT_THROW(static_cast<void>(left.Multiply(right)), std::invalid_argument);
+  left(0U, 1U) = 1.0;
+  right(1U, 0U) = std::numeric_limits<double>::infinity();
+  EXPECT_THROW(static_cast<void>(left.Multiply(right)), std::invalid_argument);
+  right(1U, 0U) = 2.0;
+  EXPECT_DOUBLE_EQ(left.Multiply(rhs)[0U], 4.0);
+  EXPECT_DOUBLE_EQ(left.Multiply(right)(1U, 0U), 4.0);
+}
+
 }  // namespace
 }  // namespace fesa