From 9c4e6854ae40536cbb882befd8274610afeed4d3 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Wed, 14 Oct 2026 16:02:17 +0900 Subject: [PATCH] fix(math): preflight finite dense operands at dispatch boundary --- include/fesa/math/matrix.h | 4 ++-- src/fesa/analysis/eigen_reduction.cpp | 7 ++++++- .../analysis/linear_frequency_procedure.cpp | 8 ++++++++ src/fesa/math/dense_blas_internal.h | 2 +- src/fesa/math/matrix.cpp | 13 +++++++++++++ .../solvers/eigen/mkl_feast_eigen_solver.cpp | 11 ++++++++++- tests/unit/math/dense_blas_internal_test.cpp | 7 +++++++ tests/unit/math/matrix_test.cpp | 18 ++++++++++++++++++ 8 files changed, 65 insertions(+), 5 deletions(-) diff --git a/include/fesa/math/matrix.h b/include/fesa/math/matrix.h index eb700a8..af51673 100644 --- a/include/fesa/math/matrix.h +++ b/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: diff --git a/src/fesa/analysis/eigen_reduction.cpp b/src/fesa/analysis/eigen_reduction.cpp index 35b2ccb..2b172e4 100644 --- a/src/fesa/analysis/eigen_reduction.cpp +++ b/src/fesa/analysis/eigen_reduction.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include #include @@ -226,7 +227,7 @@ Result BuildEigenReduction( const SparseMatrix& stiffness, const std::vector& mass, const std::vector& frames, const std::vector& mask, const std::vector& 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::Failure(valid); @@ -401,6 +402,10 @@ Result BuildEigenReduction( return Result::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::Failure( + Failure("invalid-eigen-mass-space", 0U, + "Eigen product operands must be finite and shape compatible.")); } Vector RestoreScaledMode(const EigenReduction& reduction, const Matrix& vectors, diff --git a/src/fesa/analysis/linear_frequency_procedure.cpp b/src/fesa/analysis/linear_frequency_procedure.cpp index c886871..fe3110c 100644 --- a/src/fesa/analysis/linear_frequency_procedure.cpp +++ b/src/fesa/analysis/linear_frequency_procedure.cpp @@ -275,6 +275,10 @@ Result 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::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 LinearFrequencyProcedureExecutor::Execute( std::move(diagnostics))); } return result; + } catch (const std::invalid_argument&) { + return Result::Failure( + Failure("invalid-eigenvalue-result", "modal-workspace", + "Modal product inputs must be finite and shape compatible.")); } catch (const std::bad_alloc&) { return Result::Failure( Failure("eigen-allocation-failed", "modal-workspace", diff --git a/src/fesa/math/dense_blas_internal.h b/src/fesa/math/dense_blas_internal.h index 5934506..78f9e77 100644 --- a/src/fesa/math/dense_blas_internal.h +++ b/src/fesa/math/dense_blas_internal.h @@ -31,7 +31,7 @@ class ScopedSingleEquationThreads { Result 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. diff --git a/src/fesa/math/matrix.cpp b/src/fesa/math/matrix.cpp index fa9ce89..a28e167 100644 --- a/src/fesa/math/matrix.cpp +++ b/src/fesa/math/matrix.cpp @@ -1,5 +1,6 @@ #include "fesa/math/matrix.h" +#include #include #include @@ -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) { diff --git a/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp b/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp index 5104a6a..7e3a237 100644 --- a/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp +++ b/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp @@ -1194,8 +1194,13 @@ Result 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::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); @@ -1245,6 +1250,10 @@ Result MklFeastEigenSolver::SolveImpl( solve_upper, upper_search_termination}; return Result::Success(std::move(result)); +} catch (const std::invalid_argument&) { + return Result::Failure( + Failure("eigenpair-validation-failed", "independent-action", + "Modal product inputs must be finite and shape compatible.")); } catch (const std::length_error&) { return Result::Failure( Failure("feast-integer-overflow", "independent-action", diff --git a/tests/unit/math/dense_blas_internal_test.cpp b/tests/unit/math/dense_blas_internal_test.cpp index c365287..eaebe46 100644 --- a/tests/unit/math/dense_blas_internal_test.cpp +++ b/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::quiet_NaN(); + EXPECT_THROW(static_cast(MultiplyTransposeLeft(left, right)), + std::invalid_argument); + left(0U, 0U) = 1.0; + right(0U, 0U) = std::numeric_limits::infinity(); + EXPECT_THROW(static_cast(MultiplyTransposeLeft(left, right)), + std::invalid_argument); } TEST(DenseBlasInternal, DiagonalProductKeepsEntrywiseCancellationAndChecks) { diff --git a/tests/unit/math/matrix_test.cpp b/tests/unit/math/matrix_test.cpp index 575b848..b315b44 100644 --- a/tests/unit/math/matrix_test.cpp +++ b/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::infinity(); + EXPECT_THROW(static_cast(left.Multiply(rhs)), std::invalid_argument); + rhs[1U] = 2.0; + left(0U, 1U) = std::numeric_limits::quiet_NaN(); + EXPECT_THROW(static_cast(left.Multiply(rhs)), std::invalid_argument); + EXPECT_THROW(static_cast(left.Multiply(right)), std::invalid_argument); + left(0U, 1U) = 1.0; + right(1U, 0U) = std::numeric_limits::infinity(); + EXPECT_THROW(static_cast(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