From 0e07b4df080db9ee266ce9bf4e7d084c643eb959 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Fri, 9 Oct 2026 02:24:55 +0900 Subject: [PATCH] fix(eigen): reject nonfinite modal validation scales --- .../analysis/solid_frequency_procedure.cpp | 24 +++++++++++++++---- .../solvers/eigen/mkl_feast_eigen_solver.cpp | 9 +++++-- .../solid_frequency_procedure_test.cpp | 16 +++++++++++++ 3 files changed, 43 insertions(+), 6 deletions(-) diff --git a/src/fesa/analysis/solid_frequency_procedure.cpp b/src/fesa/analysis/solid_frequency_procedure.cpp index 556f5e7..0c92757 100644 --- a/src/fesa/analysis/solid_frequency_procedure.cpp +++ b/src/fesa/analysis/solid_frequency_procedure.cpp @@ -265,11 +265,27 @@ Result ExecuteCandidate(const Domain& domain) { Vector residual = kv; residual.Axpy(-lambda, mv); const double norm_phi = vector.Norm(); - const double denominator = (norm_k + std::abs(lambda) * norm_m) * norm_phi; - const double eta_r = Ratio(residual.Norm(), denominator); + const double residual_norm = residual.Norm(); + const double lambda_mass = lambda * modal_mass; + const double energy_error = std::abs(energy - lambda_mass); + const double matrix_scale = norm_k + std::abs(lambda) * norm_m; + const double denominator = matrix_scale * norm_phi; + const double energy_denominator = denominator * norm_phi; + // Finite / infinity is zero in IEEE arithmetic but cannot certify a pair. + // Validate both original-pencil scales before either division, preserving + // the exact-zero branch only for representable finite numerator/scale. + if (!std::isfinite(norm_k) || !std::isfinite(norm_m) || + !std::isfinite(norm_phi) || !std::isfinite(residual_norm) || + !std::isfinite(lambda_mass) || !std::isfinite(energy_error) || + !std::isfinite(matrix_scale) || !std::isfinite(denominator) || + !std::isfinite(energy_denominator)) + return Result::Failure( + Failure("solid-eigenpair-validation-failed", + "A nonfinite original K/M validation numerator or scale " + "cannot certify a pair.")); + const double eta_r = Ratio(residual_norm, denominator); const double eta_norm = std::abs(modal_mass - 1.0); - const double eta_energy = - Ratio(std::abs(energy - lambda * modal_mass), denominator * norm_phi); + const double eta_energy = Ratio(energy_error, energy_denominator); if (!std::isfinite(lambda) || !std::isfinite(eta_r) || eta_r > 1e-10 || !std::isfinite(eta_norm) || eta_norm > 1e-10 || !std::isfinite(eta_energy) || eta_energy > 1e-10 || diff --git a/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp b/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp index f72bb2e..af2179f 100644 --- a/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp +++ b/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp @@ -1091,12 +1091,17 @@ Result MklFeastEigenSolver::SolveImpl( euclidean = std::hypot(euclidean, vectors[row + input * n]); residual_denominator = (norm_k + std::abs(mu) * norm_m) * euclidean; } + if (!std::isfinite(residual_norm) || !std::isfinite(residual_denominator)) + return Result::Failure( + TraceFailure("eigenpair-validation-failed", std::to_string(input), + "A nonfinite independent residual numerator or scale " + "cannot certify a pair.")); const double rho_b = residual_denominator == 0.0 ? (residual_norm == 0.0 ? 0.0 : (std::numeric_limits::infinity)()) : std::sqrt(residual_norm) / residual_denominator; - if (rho_b > 1.0e-10) { + if (!std::isfinite(rho_b) || rho_b > 1.0e-10) { if (mass != nullptr) return Result::Failure(TraceFailure( "eigenpair-validation-failed", std::to_string(input), @@ -1198,7 +1203,7 @@ Result MklFeastEigenSolver::SolveImpl( } gram_error = std::sqrt(gram_squared / (mass != nullptr ? 1.0 : static_cast(count))); - if (gram_error > 1.0e-10) { + if (!std::isfinite(gram_error) || gram_error > 1.0e-10) { return Result::Failure(TraceFailure( "eigenpair-validation-failed", "gram", "The normalized all-c Frobenius Gram error exceeds 1e-10: " + diff --git a/tests/unit/analysis/solid_frequency_procedure_test.cpp b/tests/unit/analysis/solid_frequency_procedure_test.cpp index a5a2bba..62479b5 100644 --- a/tests/unit/analysis/solid_frequency_procedure_test.cpp +++ b/tests/unit/analysis/solid_frequency_procedure_test.cpp @@ -276,3 +276,19 @@ TEST(SolidFrequencyProcedure, // Original modal Value output is pending. This proves only our independent // checks on the supplied input; it is not an external frequency comparison. } + +TEST(SolidFrequencyProcedure, RejectsOverflowingOriginalPencilValidationScale) { + // Every local K entry, eigenvalue, and spectral enclosure is representable. + // The all-odd x Walsh mode still overflows the original residual/energy + // denominators: finite / infinity must never become passing zero evidence. + auto domain = Brick(24U, "", 1.0, 2.0, 3e306); + ASSERT_TRUE(domain.HasValue()) << StatusText(domain.GetStatus()); + const auto result = + fesa::SolidFrequencyProcedureExecutor{}.Execute(domain.Value()); + ASSERT_FALSE(result.HasValue()); + ASSERT_FALSE(result.GetStatus().Diagnostics().empty()); + EXPECT_EQ(result.GetStatus().Diagnostics()[0].code, + "solid-eigenpair-validation-failed"); + EXPECT_NE(result.GetStatus().Diagnostics()[0].message.find("nonfinite"), + std::string::npos); +}