Fix flaky product and eigensolver_selfadjoint tests

libeigen/eigen!2326

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
This commit is contained in:
Rasmus Munk Larsen
2026-03-20 13:44:03 -07:00
parent a72172e563
commit a0b16a7e1b
2 changed files with 39 additions and 16 deletions

View File

@@ -26,7 +26,7 @@ void selfadjointeigensolver_essential_check(const MatrixType& m) {
VERIFY_IS_EQUAL(eiSymm.info(), Success);
RealScalar scaling = m.cwiseAbs().maxCoeff();
RealScalar unitary_error_factor = RealScalar(16);
RealScalar unitary_error_factor = RealScalar(32);
if (scaling < (std::numeric_limits<RealScalar>::min)()) {
VERIFY(eiSymm.eigenvalues().cwiseAbs().maxCoeff() <= (std::numeric_limits<RealScalar>::min)());
@@ -80,7 +80,7 @@ void selfadjointeigensolver(const MatrixType& m) {
MatrixType symmA = a.adjoint() * a + a1.adjoint() * a1;
MatrixType symmC = symmA;
svd_fill_random(symmA, Symmetric);
svd_fill_random(symmA, SelfAdjoint);
symmA.template triangularView<StrictlyUpper>().setZero();
symmC.template triangularView<StrictlyUpper>().setZero();

View File

@@ -111,8 +111,12 @@ void product(const MatrixType& m) {
VERIFY_IS_APPROX(m3, m1 * (m1.transpose() * m2));
// continue testing Product.h: distributivity
VERIFY_IS_APPROX(square * (m1 + m2), square * m1 + square * m2);
VERIFY_IS_APPROX(square * (m1 - m2), square * m1 - square * m2);
{
// Increase tolerance, since coefficients here can get relatively large.
RealScalar tol = RealScalar(2) * get_test_precision(m1);
VERIFY(verifyIsApprox(square * (m1 + m2), square * m1 + square * m2, tol));
VERIFY(verifyIsApprox(square * (m1 - m2), square * m1 - square * m2, tol));
}
// continue testing Product.h: compatibility with ScalarMultiple.h
VERIFY_IS_APPROX(s1 * (square * m1), (s1 * square) * m1);
@@ -130,7 +134,10 @@ void product(const MatrixType& m) {
// test the previous tests were not screwed up because operator* returns 0
// (we use the more accurate default epsilon)
if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1) {
// Skip non-commutativity check for very low-precision types (e.g. bfloat16) where
// random small matrices can produce products that are approximately equal.
if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1 &&
NumTraits<RealScalar>::dummy_precision() < RealScalar(0.04)) {
VERIFY(areNotApprox(m1.transpose() * m2, m2.transpose() * m1, not_approx_epsilon));
}
@@ -138,7 +145,8 @@ void product(const MatrixType& m) {
res = square;
res.noalias() += m1 * m2.transpose();
VERIFY_IS_APPROX(res, square + m1 * m2.transpose());
if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1) {
if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1 &&
NumTraits<RealScalar>::dummy_precision() < RealScalar(0.04)) {
VERIFY(areNotApprox(res, square + m2 * m1.transpose(), not_approx_epsilon));
}
vcres = vc2;
@@ -149,7 +157,8 @@ void product(const MatrixType& m) {
res = square;
res.noalias() -= m1 * m2.transpose();
VERIFY_IS_APPROX(res, square - (m1 * m2.transpose()));
if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1) {
if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1 &&
NumTraits<RealScalar>::dummy_precision() < RealScalar(0.04)) {
VERIFY(areNotApprox(res, square - m2 * m1.transpose(), not_approx_epsilon));
}
vcres = vc2;
@@ -183,21 +192,35 @@ void product(const MatrixType& m) {
res.noalias() -= square - m1 * m2.transpose();
VERIFY_IS_APPROX(res, square - m1 * m2.transpose());
tm1 = m1;
VERIFY_IS_APPROX(tm1.transpose() * v1, m1.transpose() * v1);
VERIFY_IS_APPROX(v1.transpose() * tm1, v1.transpose() * m1);
// Row-major vs col-major GEMV: different accumulation orders produce
// rounding differences bounded by O(sqrt(n)) * epsilon per inner product.
// tm1.transpose() * v1 has inner dimension rows; v1^T * tm1 has inner dim cols.
{
RealScalar gemv_tol = (std::max)(get_test_precision(m1),
numext::sqrt(RealScalar((std::max)(rows, cols))) * NumTraits<Scalar>::epsilon());
tm1 = m1;
VERIFY(verifyIsApprox(tm1.transpose() * v1, m1.transpose() * v1, gemv_tol));
VERIFY(verifyIsApprox(v1.transpose() * tm1, v1.transpose() * m1, gemv_tol));
}
// test submatrix and matrix/vector product
for (int i = 0; i < rows; ++i) res.row(i) = m1.row(i) * m2.transpose();
VERIFY_IS_APPROX(res, m1 * m2.transpose());
// the other way round:
for (int i = 0; i < rows; ++i) res.col(i) = m1 * m2.transpose().col(i);
VERIFY_IS_APPROX(res, m1 * m2.transpose());
// Row-by-row vs full GEMM: different evaluation strategies can differ
// by O(sqrt(n)) * epsilon for low-precision types.
{
RealScalar prod_tol =
(std::max)(get_test_precision(m1), numext::sqrt(RealScalar(cols)) * NumTraits<Scalar>::epsilon());
for (int i = 0; i < rows; ++i) res.row(i) = m1.row(i) * m2.transpose();
VERIFY(verifyIsApprox(res, m1 * m2.transpose(), prod_tol));
// the other way round:
for (int i = 0; i < rows; ++i) res.col(i) = m1 * m2.transpose().col(i);
VERIFY(verifyIsApprox(res, m1 * m2.transpose(), prod_tol));
}
res2 = square2;
res2.noalias() += m1.transpose() * m2;
VERIFY_IS_APPROX(res2, square2 + m1.transpose() * m2);
if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1) {
if (!NumTraits<Scalar>::IsInteger && (std::min)(rows, cols) > 1 &&
NumTraits<RealScalar>::dummy_precision() < RealScalar(0.04)) {
VERIFY(areNotApprox(res2, square2 + m2.transpose() * m1, not_approx_epsilon));
}