From 590a93793f0e547f2d6c331341f4892247f0ff81 Mon Sep 17 00:00:00 2001 From: "KOKO\\Mimi" Date: Wed, 14 Oct 2026 15:40:52 +0900 Subject: [PATCH] fix(math): preserve certified sums and complete MKL algebra --- include/fesa/math/sparse_matrix.h | 4 +- src/fesa/analysis/eigen_reduction.cpp | 111 ++++++++++-------- .../analysis/thermal_transient_algebra.cpp | 55 ++++----- src/fesa/math/dense_blas_internal.cpp | 62 ++++++++++ src/fesa/math/dense_blas_internal.h | 14 +++ src/fesa/math/sparse_matrix.cpp | 4 + .../solvers/eigen/mkl_feast_eigen_solver.cpp | 19 +-- tests/unit/math/dense_blas_internal_test.cpp | 51 ++++++++ tests/unit/math/sparse_matrix_test.cpp | 11 +- 9 files changed, 236 insertions(+), 95 deletions(-) diff --git a/include/fesa/math/sparse_matrix.h b/include/fesa/math/sparse_matrix.h index bbebfc1..1aabf21 100644 --- a/include/fesa/math/sparse_matrix.h +++ b/include/fesa/math/sparse_matrix.h @@ -47,11 +47,11 @@ class SparseMatrix { const std::vector& Values() const noexcept; /// @brief Multiplies canonical CSR by a dense vector through private MKL. - /// @throws std::invalid_argument if dimensions or CSR are invalid. + /// @throws std::invalid_argument if dimensions, CSR, or RHS 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. + /// Finite-input overflow is retained for consumer output 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 7edb030..b7e98fc 100644 --- a/src/fesa/analysis/eigen_reduction.cpp +++ b/src/fesa/analysis/eigen_reduction.cpp @@ -122,12 +122,11 @@ Result 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 MakeLift(const Matrix& mass, return Result::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::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::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::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 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 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::Failure(Failure( "eigen-mass-conditioning", node, "Nodal whitening is nonfinite.")); @@ -280,19 +273,23 @@ Result 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 BuildEigenReduction( for (std::size_t a = 0; a < 6U; ++a) { const auto row = node * 6U + a; if (mask[row] != 0U) continue; - std::map 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 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 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 projected_blocks; for (const auto& entry : blocks) { projected_blocks.emplace( diff --git a/src/fesa/analysis/thermal_transient_algebra.cpp b/src/fesa/analysis/thermal_transient_algebra.cpp index be029fa..2784d43 100644 --- a/src/fesa/analysis/thermal_transient_algebra.cpp +++ b/src/fesa/analysis/thermal_transient_algebra.cpp @@ -8,6 +8,7 @@ #include #include "thermal_transient_bounds.h" +#include "math/dense_blas_internal.h" namespace fesa::thermal_transient_internal { namespace { @@ -117,14 +118,9 @@ Result Build(const ThermalSystem& s, const ThermalPartition& p, constrained[p.constrained[i]] = i; SparsePattern pattern{{0}, {}}; std::vector coo; - SparsePattern rhs_pattern{{0U}, {}}; - std::vector capacity_rhs, conductivity_rhs; - Vector capacity_input = old; - Vector prescribed_full{m}; - std::vector 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 exact_a, exact_b(n); for (std::size_t f = 0; f < n; ++f) { @@ -135,43 +131,39 @@ Result 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::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::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::Failure(a.GetStatus()); if (!Finite(b)) @@ -493,6 +485,9 @@ Result SolveThermalIncrement( status = Check(s, p, old, delta, equations.Value(), candidate); if (!status.IsOk()) return Result::Failure(status); return Result::Success(std::move(candidate)); +} catch (const std::invalid_argument&) { + return Result::Failure( + Failure("Nonfinite backward Euler RHS.")); } catch (const std::length_error&) { return Result::Failure( Failure("Thermal sparse action exceeds the backend integer range.")); diff --git a/src/fesa/math/dense_blas_internal.cpp b/src/fesa/math/dense_blas_internal.cpp index 6c13333..40d53d4 100644 --- a/src/fesa/math/dense_blas_internal.cpp +++ b/src/fesa/math/dense_blas_internal.cpp @@ -1,10 +1,72 @@ #include "math/dense_blas_internal.h" +#include #include #include +#include +#include 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 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."}; diff --git a/src/fesa/math/dense_blas_internal.h b/src/fesa/math/dense_blas_internal.h index e96a05e..21e7749 100644 --- a/src/fesa/math/dense_blas_internal.h +++ b/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 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. diff --git a/src/fesa/math/sparse_matrix.cpp b/src/fesa/math/sparse_matrix.cpp index 5c8bb75..988964e 100644 --- a/src/fesa/math/sparse_matrix.cpp +++ b/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(SparseSize(values_.size())); diff --git a/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp b/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp index 1b42810..5104a6a 100644 --- a/src/fesa/solvers/eigen/mkl_feast_eigen_solver.cpp +++ b/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 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::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 MklFeastEigenSolver::SolveImpl( } double gram_error = 0.0; double gram_squared = 0.0; - std::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; diff --git a/tests/unit/math/dense_blas_internal_test.cpp b/tests/unit/math/dense_blas_internal_test.cpp index c6cf25b..6e32e17 100644 --- a/tests/unit/math/dense_blas_internal_test.cpp +++ b/tests/unit/math/dense_blas_internal_test.cpp @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -57,5 +58,55 @@ TEST(DenseBlasInternal, TransposeProductChecksShapeAndKnownRectangularAction) { EXPECT_DOUBLE_EQ(empty(1U, 2U), 0.0); } +TEST(DenseBlasInternal, DiagonalProductKeepsEntrywiseCancellationAndChecks) { + Vector diagonal{3U}; + diagonal[0U] = 4; + diagonal[1U] = 2; + diagonal[2U] = -1; + Vector input{3U}; + const double large = std::ldexp(1.0, 56); + input[0U] = 1; + input[1U] = -large; + input[2U] = 2 * large; + const auto product = MultiplyDiagonal(diagonal, input); + EXPECT_DOUBLE_EQ(product[0U] + (product[1U] - product[2U]), 4); + EXPECT_EQ(MultiplyDiagonal(Vector{0U}, Vector{0U}).Size(), 0U); + EXPECT_THROW(static_cast(MultiplyDiagonal(diagonal, Vector{2U})), + std::invalid_argument); + input[0U] = std::numeric_limits::infinity(); + EXPECT_THROW(static_cast(MultiplyDiagonal(diagonal, input)), + std::invalid_argument); + input[0U] = 1; + diagonal[0U] = std::numeric_limits::quiet_NaN(); + EXPECT_THROW(static_cast(MultiplyDiagonal(diagonal, input)), + std::invalid_argument); +} + +TEST(DenseBlasInternal, CholeskyAndInverseMatchKnownFactorsAndKeepFailedOutput) { + Matrix matrix{2U, 2U}; + matrix(0U, 0U) = 4; + matrix(0U, 1U) = 2; + matrix(1U, 0U) = 2; + matrix(1U, 1U) = 5; + Matrix lower{0U, 0U}, inverse{0U, 0U}; + ASSERT_TRUE(CholeskyLowerAndInverse(matrix, lower, inverse)); + EXPECT_DOUBLE_EQ(lower(0U, 0U), 2); + EXPECT_DOUBLE_EQ(lower(1U, 0U), 1); + EXPECT_DOUBLE_EQ(lower(1U, 1U), 2); + EXPECT_DOUBLE_EQ(lower(0U, 1U), 0); + EXPECT_DOUBLE_EQ(inverse(0U, 0U), .5); + EXPECT_DOUBLE_EQ(inverse(1U, 0U), -.25); + EXPECT_DOUBLE_EQ(inverse(1U, 1U), .5); + matrix(1U, 1U) = 0; + EXPECT_FALSE(CholeskyLowerAndInverse(matrix, lower, inverse)); + EXPECT_DOUBLE_EQ(lower(1U, 1U), 2); + EXPECT_DOUBLE_EQ(inverse(1U, 1U), .5); + EXPECT_THROW(static_cast(CholeskyLowerAndInverse( + Matrix{2U, 1U}, lower, inverse)), std::invalid_argument); + matrix(0U, 0U) = std::numeric_limits::quiet_NaN(); + EXPECT_THROW(static_cast(CholeskyLowerAndInverse(matrix, lower, inverse)), + std::invalid_argument); +} + } // namespace } // namespace fesa::dense_blas_internal diff --git a/tests/unit/math/sparse_matrix_test.cpp b/tests/unit/math/sparse_matrix_test.cpp index 525453e..7185f91 100644 --- a/tests/unit/math/sparse_matrix_test.cpp +++ b/tests/unit/math/sparse_matrix_test.cpp @@ -181,16 +181,19 @@ TEST(SparseAssembly, PreservesZeroShapesAndEmptyRows) { } } -TEST(SparseAssembly, KeepsIeeeNonfiniteResultsAndRecoversAfterBadRhs) { +TEST(SparseAssembly, RejectsNonfiniteRhsAndRecoversAfterFiniteOverflow) { 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::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. + EXPECT_THROW(static_cast(matrix.Value().Multiply(rhs)), + std::invalid_argument); + rhs[1U] = 4.0; rhs[0U] = (std::numeric_limits::max)(); EXPECT_TRUE(std::isinf(matrix.Value().Multiply(rhs)[0U])); + rhs[0U] = (std::numeric_limits::quiet_NaN)(); + EXPECT_THROW(static_cast(matrix.Value().Multiply(rhs)), + std::invalid_argument); EXPECT_THROW(static_cast(matrix.Value().Multiply(Vector{1U})), std::invalid_argument); rhs[0U] = 3.0;