From 8346cc3410ca63f85eff516a82bce42b849ad9e0 Mon Sep 17 00:00:00 2001 From: Rasmus Munk Larsen Date: Sat, 4 Apr 2026 18:02:37 -0700 Subject: [PATCH] Fix three bugs in SelfAdjointEigenSolver and improve test coverage Bug fixes: 1. computeDirect 3x3: eigenvalues from computeRoots() are theoretically sorted via the trigonometric formula, but floating-point rounding (especially in float) can break the ordering. Add a 3-element sorting network after computeRoots() to guarantee sorted output. 2. compute(): Inf input silently returns Success with garbage results, unlike NaN which correctly returns NoConvergence. Add an isfinite() check on the scaling factor (which is maxCoeff of the matrix) to detect Inf/NaN early and return NoConvergence. 3. computeFromTridiagonal(): lacks the equilibration scaling that compute() applies, causing overflow and NoConvergence for tridiagonal matrices with large entries. Add the same scale-to-[-1,1] pattern. New tests: - Eigenvalue sorting verification (both iterative and direct solvers) - Repeated/degenerate eigenvalues (all equal, multiplicity n-1, two clusters, nearly repeated separated by O(epsilon)) - Extreme eigenvalue ranges (high condition number spanning many orders of magnitude, near-underflow, near-overflow, mixed positive/negative, rank-deficient with zero eigenvalue) - computeFromTridiagonal with large and tiny values - Diagonal matrices (eigenvalues must match sorted diagonal) - operatorInverseSqrt accuracy (sqrtA*invSqrtA=I, invSqrtA*A*invSqrtA=I, symmetry) - RowMajor storage for computeDirect (2x2, 3x3, float, double) and iterative solver (dynamic RowMajor) - Inf input detection - Tridiagonal structure verification (off-tridiagonal entries are zero) - Direct solver stress tests: 3x3 (near-planar covariance, triple eigenvalue, double eigenvalue, large off-diagonal, nearly singular) and 2x2 (equal eigenvalues, tiny off-diagonal, huge diagonal ratio, anti-diagonal dominant, negative entries) - Tightened unitary tolerance from fixed 32*test_precision to 4*n*epsilon (scales with matrix size) - Fixed typo: "expponential" -> "exponential" Co-Authored-By: Claude Opus 4.6 (1M context) --- .../src/Eigenvalues/SelfAdjointEigenSolver.h | 51 +- test/eigensolver_selfadjoint.cpp | 473 +++++++++++++++++- 2 files changed, 516 insertions(+), 8 deletions(-) diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index b7b57ec49..6c46e840e 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -451,6 +451,13 @@ EIGEN_DEVICE_FUNC SelfAdjointEigenSolver& SelfAdjointEigenSolver(); RealScalar scale = mat.cwiseAbs().maxCoeff(); + if (!(numext::isfinite)(scale)) { + // Input contains Inf or NaN. + m_info = NoConvergence; + m_isInitialized = true; + m_eigenvectorsOk = false; + return *this; + } if (numext::is_exactly_zero(scale)) scale = RealScalar(1); mat.template triangularView() /= scale; m_subdiag.resize(n - 1); @@ -470,16 +477,36 @@ EIGEN_DEVICE_FUNC SelfAdjointEigenSolver& SelfAdjointEigenSolver SelfAdjointEigenSolver& SelfAdjointEigenSolver::computeFromTridiagonal( const RealVectorType& diag, const SubDiagonalType& subdiag, int options) { - // TODO : Add an option to scale the values beforehand bool computeEigenvectors = (options & ComputeEigenvectors) == ComputeEigenvectors; m_eivalues = diag; m_subdiag = subdiag; + + // Scale the tridiagonal matrix to [-1:1] to avoid over- and underflow, + // just like compute() does for the full matrix. + RealScalar scale = m_eivalues.cwiseAbs().maxCoeff(); + if (m_subdiag.size() > 0) scale = numext::maxi(scale, m_subdiag.cwiseAbs().maxCoeff()); + if (!(numext::isfinite)(scale)) { + m_info = NoConvergence; + m_isInitialized = true; + m_eigenvectorsOk = false; + return *this; + } + if (numext::is_exactly_zero(scale)) scale = RealScalar(1); + const bool needsScaling = scale != RealScalar(1); + if (needsScaling) { + m_eivalues /= scale; + m_subdiag /= scale; + } + if (computeEigenvectors) { m_eivec.setIdentity(diag.size(), diag.size()); } m_info = internal::computeFromTridiagonal_impl(m_eivalues, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec); + // Scale back the eigenvalues. + if (needsScaling) m_eivalues *= scale; + m_isInitialized = true; m_eigenvectorsOk = computeEigenvectors; return *this; @@ -662,6 +689,28 @@ struct direct_selfadjoint_eigenvalues { // compute the eigenvalues computeRoots(scaledMat, eivals); + // computeRoots produces theoretically sorted roots, but floating-point + // rounding in the trigonometric formulas can break the ordering. + // Enforce sorting with a branchless min/max network (3 elements). + { + Scalar tmp; + if (eivals(0) > eivals(1)) { + tmp = eivals(0); + eivals(0) = eivals(1); + eivals(1) = tmp; + } + if (eivals(1) > eivals(2)) { + tmp = eivals(1); + eivals(1) = eivals(2); + eivals(2) = tmp; + } + if (eivals(0) > eivals(1)) { + tmp = eivals(0); + eivals(0) = eivals(1); + eivals(1) = tmp; + } + } + // compute the eigenvectors if (computeEigenvectors) { if ((eivals(2) - eivals(0)) <= Eigen::NumTraits::epsilon()) { diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp index 4126d5204..11eded8e1 100644 --- a/test/eigensolver_selfadjoint.cpp +++ b/test/eigensolver_selfadjoint.cpp @@ -25,8 +25,8 @@ void selfadjointeigensolver_essential_check(const MatrixType& m) { SelfAdjointEigenSolver eiSymm(m); VERIFY_IS_EQUAL(eiSymm.info(), Success); + Index n = m.cols(); RealScalar scaling = m.cwiseAbs().maxCoeff(); - RealScalar unitary_error_factor = RealScalar(32); if (scaling < (std::numeric_limits::min)()) { VERIFY(eiSymm.eigenvalues().cwiseAbs().maxCoeff() <= (std::numeric_limits::min)()); @@ -35,7 +35,18 @@ void selfadjointeigensolver_essential_check(const MatrixType& m) { (eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal()) / scaling); } VERIFY_IS_APPROX(m.template selfadjointView().eigenvalues(), eiSymm.eigenvalues()); - VERIFY(eiSymm.eigenvectors().isUnitary(test_precision() * unitary_error_factor)); + + // Eigenvectors must be unitary. Use a tolerance proportional to n*epsilon, + // which is the expected rounding error for Householder-based orthogonal transformations. + RealScalar unitary_tol = RealScalar(4) * RealScalar(numext::maxi(Index(1), n)) * NumTraits::epsilon(); + // But don't go below the test_precision floor (matters for float). + unitary_tol = numext::maxi(unitary_tol, test_precision()); + VERIFY(eiSymm.eigenvectors().isUnitary(unitary_tol)); + + // Verify eigenvalues are sorted in non-decreasing order. + for (Index i = 1; i < n; ++i) { + VERIFY(eiSymm.eigenvalues()(i) >= eiSymm.eigenvalues()(i - 1)); + } if (m.cols() <= 4) { SelfAdjointEigenSolver eiDirect; @@ -58,7 +69,13 @@ void selfadjointeigensolver_essential_check(const MatrixType& m) { VERIFY_IS_APPROX(m.template selfadjointView().eigenvalues() / scaling, eiDirect.eigenvalues() / scaling); } - VERIFY(eiDirect.eigenvectors().isUnitary(test_precision() * unitary_error_factor)); + // Direct solver eigenvectors must also be unitary. + VERIFY(eiDirect.eigenvectors().isUnitary(unitary_tol)); + + // Direct solver eigenvalues must also be sorted. + for (Index i = 1; i < n; ++i) { + VERIFY(eiDirect.eigenvalues()(i) >= eiDirect.eigenvalues()(i - 1)); + } } } @@ -149,9 +166,15 @@ void selfadjointeigensolver(const MatrixType& m) { VERIFY_IS_APPROX(tridiag.diagonal(), tridiag.matrixT().diagonal()); VERIFY_IS_APPROX(tridiag.subDiagonal(), tridiag.matrixT().template diagonal<-1>()); Matrix T = tridiag.matrixT(); - if (rows > 1 && cols > 1) { - // FIXME check that upper and lower part are 0: - // VERIFY(T.topRightCorner(rows-2, cols-2).template triangularView().isZero()); + if (rows > 2) { + // Verify that the tridiagonal matrix is actually tridiagonal (zero outside the three central diagonals). + for (Index i = 0; i < rows; ++i) { + for (Index j = 0; j < cols; ++j) { + if (numext::abs(i - j) > 1) { + VERIFY(numext::is_exactly_zero(T(i, j))); + } + } + } } VERIFY_IS_APPROX(tridiag.diagonal(), T.diagonal()); VERIFY_IS_APPROX(tridiag.subDiagonal(), T.template diagonal<1>()); @@ -170,7 +193,7 @@ void selfadjointeigensolver(const MatrixType& m) { eiSymmTridiag.eigenvectors().real().transpose()); } - // Test matrix expponential from eigendecomposition. + // Test matrix exponential from eigendecomposition. // First scale to avoid overflow. symmB = symmB / symmB.norm(); eiSymm.compute(symmB); @@ -202,6 +225,320 @@ void selfadjointeigensolver(const MatrixType& m) { } } +// Test matrices with exact eigenvalue multiplicities. +template +void selfadjointeigensolver_repeated_eigenvalues(const MatrixType& m) { + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + Index n = m.rows(); + if (n < 2) return; + + // Create a random unitary matrix via QR. + MatrixType q = MatrixType::Random(n, n); + HouseholderQR qr(q); + q = qr.householderQ(); + + // All eigenvalues equal (scalar multiple of identity). + { + RealScalar lambda = internal::random(-10, 10); + MatrixType A = lambda * MatrixType::Identity(n, n); + selfadjointeigensolver_essential_check(A); + } + + // Eigenvalue of multiplicity n-1 (one distinct, rest equal). + { + Matrix d = Matrix::Constant(n, RealScalar(3)); + d(0) = RealScalar(-2); + MatrixType A = (q * d.template cast().asDiagonal() * q.adjoint()).eval(); + A.template triangularView().setZero(); + selfadjointeigensolver_essential_check(A); + } + + // Two clusters: first half one value, second half another. + if (n >= 4) { + Matrix d(n); + for (Index i = 0; i < n / 2; ++i) d(i) = RealScalar(1); + for (Index i = n / 2; i < n; ++i) d(i) = RealScalar(5); + MatrixType A = (q * d.template cast().asDiagonal() * q.adjoint()).eval(); + A.template triangularView().setZero(); + selfadjointeigensolver_essential_check(A); + } + + // Nearly repeated eigenvalues: separated by O(epsilon). + { + Matrix d(n); + for (Index i = 0; i < n; ++i) { + d(i) = RealScalar(1) + RealScalar(i) * NumTraits::epsilon() * RealScalar(10); + } + MatrixType A = (q * d.template cast().asDiagonal() * q.adjoint()).eval(); + A.template triangularView().setZero(); + selfadjointeigensolver_essential_check(A); + } +} + +// Test matrices with extreme condition numbers and eigenvalue ranges. +template +void selfadjointeigensolver_extreme_eigenvalues(const MatrixType& m) { + using std::pow; + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + Index n = m.rows(); + if (n < 2) return; + + // Create a random unitary matrix. + MatrixType q = MatrixType::Random(n, n); + HouseholderQR qr(q); + q = qr.householderQ(); + + // Eigenvalues spanning many orders of magnitude (high condition number). + { + RealScalar maxExp = RealScalar(std::numeric_limits::max_exponent10) / RealScalar(4); + Matrix d(n); + for (Index i = 0; i < n; ++i) { + RealScalar exponent = -maxExp + RealScalar(2) * maxExp * RealScalar(i) / RealScalar(n - 1); + d(i) = pow(RealScalar(10), exponent); + } + MatrixType A = (q * d.template cast().asDiagonal() * q.adjoint()).eval(); + A.template triangularView().setZero(); + + SelfAdjointEigenSolver eig(A); + VERIFY_IS_EQUAL(eig.info(), Success); + // For ill-conditioned matrices we can only check the relative residual. + // ||A*V - V*D|| / ||A|| should be O(n * epsilon). + RealScalar Anorm = A.template selfadjointView().operatorNorm(); + if (Anorm > (std::numeric_limits::min)()) { + MatrixType residual = A.template selfadjointView() * eig.eigenvectors() - + eig.eigenvectors() * eig.eigenvalues().asDiagonal(); + RealScalar rel_err = residual.norm() / Anorm; + RealScalar tol = RealScalar(4) * RealScalar(n) * NumTraits::epsilon(); + VERIFY(rel_err <= tol); + } + // Eigenvalues must still be sorted. + for (Index i = 1; i < n; ++i) { + VERIFY(eig.eigenvalues()(i) >= eig.eigenvalues()(i - 1)); + } + } + + // Very tiny eigenvalues (near underflow). + { + RealScalar tiny = (std::numeric_limits::min)() * RealScalar(100); + Matrix d(n); + for (Index i = 0; i < n; ++i) { + d(i) = tiny * (RealScalar(1) + RealScalar(i)); + } + MatrixType A = (q * d.template cast().asDiagonal() * q.adjoint()).eval(); + A.template triangularView().setZero(); + selfadjointeigensolver_essential_check(A); + } + + // Very large eigenvalues (near overflow). + { + RealScalar huge = (std::numeric_limits::max)() / (RealScalar(n) * RealScalar(100)); + Matrix d(n); + for (Index i = 0; i < n; ++i) { + d(i) = huge * (RealScalar(1) + RealScalar(i) * RealScalar(0.01)); + } + MatrixType A = (q * d.template cast().asDiagonal() * q.adjoint()).eval(); + A.template triangularView().setZero(); + selfadjointeigensolver_essential_check(A); + } + + // Mix of positive and negative eigenvalues. + { + Matrix d(n); + for (Index i = 0; i < n; ++i) { + d(i) = (i % 2 == 0) ? RealScalar(i + 1) : RealScalar(-(i + 1)); + } + MatrixType A = (q * d.template cast().asDiagonal() * q.adjoint()).eval(); + A.template triangularView().setZero(); + selfadjointeigensolver_essential_check(A); + } + + // One zero eigenvalue among non-zero ones (rank-deficient). + { + Matrix d = Matrix::LinSpaced(n, RealScalar(0), RealScalar(n - 1)); + MatrixType A = (q * d.template cast().asDiagonal() * q.adjoint()).eval(); + A.template triangularView().setZero(); + selfadjointeigensolver_essential_check(A); + } +} + +// Test computeFromTridiagonal with scaled inputs (regression for missing scaling). +template +void selfadjointeigensolver_tridiagonal_scaled(const MatrixType& m) { + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + Index n = m.rows(); + if (n < 2) return; + + // Create a tridiagonal matrix with large entries. + typedef Matrix RealVectorType; + RealVectorType diag(n), subdiag(n - 1); + + // Case 1: Large values. + RealScalar scale = (std::numeric_limits::max)() / (RealScalar(n) * RealScalar(100)); + for (Index i = 0; i < n; ++i) diag(i) = scale * RealScalar(i + 1); + for (Index i = 0; i < n - 1; ++i) subdiag(i) = scale * RealScalar(0.5); + + SelfAdjointEigenSolver eig1; + eig1.computeFromTridiagonal(diag, subdiag, ComputeEigenvectors); + VERIFY_IS_EQUAL(eig1.info(), Success); + + // Reconstruct tridiagonal and check residual. + Matrix T = Matrix::Zero(n, n); + T.diagonal() = diag; + T.template diagonal<1>() = subdiag; + T.template diagonal<-1>() = subdiag; + VERIFY_IS_APPROX( + T, eig1.eigenvectors().real() * eig1.eigenvalues().asDiagonal() * eig1.eigenvectors().real().transpose()); + + // Case 2: Tiny values. + scale = (std::numeric_limits::min)() * RealScalar(100); + for (Index i = 0; i < n; ++i) diag(i) = scale * RealScalar(i + 1); + for (Index i = 0; i < n - 1; ++i) subdiag(i) = scale * RealScalar(0.5); + + SelfAdjointEigenSolver eig2; + eig2.computeFromTridiagonal(diag, subdiag, ComputeEigenvectors); + VERIFY_IS_EQUAL(eig2.info(), Success); + + // Eigenvalues-only mode should produce the same eigenvalues. + SelfAdjointEigenSolver eig2v; + eig2v.computeFromTridiagonal(diag, subdiag, EigenvaluesOnly); + VERIFY_IS_EQUAL(eig2v.info(), Success); + VERIFY_IS_APPROX(eig2.eigenvalues(), eig2v.eigenvalues()); +} + +// Test with diagonal matrices (tridiagonalization is trivial). +template +void selfadjointeigensolver_diagonal(const MatrixType& m) { + typedef typename MatrixType::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + Index n = m.rows(); + + // Random diagonal matrix. + MatrixType diag = MatrixType::Zero(n, n); + for (Index i = 0; i < n; ++i) { + diag(i, i) = internal::random(-100, 100); + } + selfadjointeigensolver_essential_check(diag); + + // The eigenvalues should be the diagonal entries, sorted. + SelfAdjointEigenSolver eig(diag); + VERIFY_IS_EQUAL(eig.info(), Success); + + Matrix expected_evals(n); + for (Index i = 0; i < n; ++i) expected_evals(i) = numext::real(diag(i, i)); + std::sort(expected_evals.data(), expected_evals.data() + n); + VERIFY_IS_APPROX(eig.eigenvalues(), expected_evals); +} + +// Test operatorInverseSqrt more thoroughly. +template +void selfadjointeigensolver_inverse_sqrt(const MatrixType& m) { + Index n = m.rows(); + if (n < 1) return; + + // Create a positive-definite matrix. + MatrixType a = MatrixType::Random(n, n); + MatrixType spd = a.adjoint() * a + MatrixType::Identity(n, n); + spd.template triangularView().setZero(); + + SelfAdjointEigenSolver eig(spd); + VERIFY_IS_EQUAL(eig.info(), Success); + + MatrixType sqrtA = eig.operatorSqrt(); + MatrixType invSqrtA = eig.operatorInverseSqrt(); + + // sqrtA * invSqrtA should be identity. + VERIFY_IS_APPROX(sqrtA * invSqrtA, MatrixType::Identity(n, n)); + + // invSqrtA * A * invSqrtA should be identity. + VERIFY_IS_APPROX(invSqrtA * spd.template selfadjointView() * invSqrtA, MatrixType::Identity(n, n)); + + // invSqrtA should be symmetric/selfadjoint. + VERIFY_IS_APPROX(invSqrtA, invSqrtA.adjoint()); +} + +// Test that RowMajor matrices work correctly with computeDirect. +template +void selfadjointeigensolver_rowmajor() { + typedef Matrix RowMajorMatrix3d; + typedef Matrix RowMajorMatrix2d; + typedef Matrix RowMajorMatrix3f; + typedef Matrix RowMajorMatrix2f; + + // 3x3 RowMajor double + { + RowMajorMatrix3d a = RowMajorMatrix3d::Random(); + RowMajorMatrix3d symmA = a.transpose() * a; + SelfAdjointEigenSolver eig; + eig.computeDirect(symmA); + VERIFY_IS_EQUAL(eig.info(), Success); + // Compare with iterative solver. + SelfAdjointEigenSolver eigRef(symmA); + VERIFY_IS_APPROX(eigRef.eigenvalues(), eig.eigenvalues()); + } + + // 2x2 RowMajor double + { + RowMajorMatrix2d a = RowMajorMatrix2d::Random(); + RowMajorMatrix2d symmA = a.transpose() * a; + SelfAdjointEigenSolver eig; + eig.computeDirect(symmA); + VERIFY_IS_EQUAL(eig.info(), Success); + SelfAdjointEigenSolver eigRef(symmA); + VERIFY_IS_APPROX(eigRef.eigenvalues(), eig.eigenvalues()); + } + + // 3x3 RowMajor float + { + RowMajorMatrix3f a = RowMajorMatrix3f::Random(); + RowMajorMatrix3f symmA = a.transpose() * a; + SelfAdjointEigenSolver eig; + eig.computeDirect(symmA); + VERIFY_IS_EQUAL(eig.info(), Success); + SelfAdjointEigenSolver eigRef(symmA); + VERIFY_IS_APPROX(eigRef.eigenvalues(), eig.eigenvalues()); + } + + // 2x2 RowMajor float + { + RowMajorMatrix2f a = RowMajorMatrix2f::Random(); + RowMajorMatrix2f symmA = a.transpose() * a; + SelfAdjointEigenSolver eig; + eig.computeDirect(symmA); + VERIFY_IS_EQUAL(eig.info(), Success); + SelfAdjointEigenSolver eigRef(symmA); + VERIFY_IS_APPROX(eigRef.eigenvalues(), eig.eigenvalues()); + } + + // Dynamic RowMajor with iterative solver + { + typedef Matrix RowMajorMatrixXd; + int s = internal::random(2, 20); + RowMajorMatrixXd a = RowMajorMatrixXd::Random(s, s); + RowMajorMatrixXd symmA = a.transpose() * a; + SelfAdjointEigenSolver eig(symmA); + VERIFY_IS_EQUAL(eig.info(), Success); + double scaling = symmA.cwiseAbs().maxCoeff(); + if (scaling > (std::numeric_limits::min)()) { + VERIFY_IS_APPROX((symmA.template selfadjointView() * eig.eigenvectors()) / scaling, + (eig.eigenvectors() * eig.eigenvalues().asDiagonal()) / scaling); + } + } +} + +// Test matrix with Inf entries returns NoConvergence (similar to NaN test). +template +void selfadjointeigensolver_inf() { + Matrix3d m; + m.setRandom(); + m = m * m.transpose(); + m(1, 1) = std::numeric_limits::infinity(); + SelfAdjointEigenSolver eig(m); + VERIFY_IS_EQUAL(eig.info(), NoConvergence); +} + template void bug_854() { Matrix3d m; @@ -234,6 +571,88 @@ void bug_1204() { SelfAdjointEigenSolver > eig(A); } +// Specific 3x3 test cases that stress the direct solver. +template +void direct_3x3_stress() { + // Near-planar point cloud covariance: two large eigenvalues, one near-zero. + { + Matrix3d m; + m << 100, 50, 0.001, 50, 100, 0.002, 0.001, 0.002, 1e-10; + selfadjointeigensolver_essential_check(m); + } + + // All equal diagonal entries (triple eigenvalue). + { + Matrix3d m = Matrix3d::Identity() * 7.0; + selfadjointeigensolver_essential_check(m); + } + + // Two exactly equal eigenvalues (from explicit construction). + { + Matrix3d q; + q << 1, 0, 0, 0, 1.0 / std::sqrt(2.0), 1.0 / std::sqrt(2.0), 0, -1.0 / std::sqrt(2.0), 1.0 / std::sqrt(2.0); + Vector3d d(1.0, 5.0, 5.0); + Matrix3d m = q * d.asDiagonal() * q.transpose(); + selfadjointeigensolver_essential_check(m); + } + + // Large off-diagonal relative to diagonal. + { + Matrix3d m; + m << 1, 1000, 1000, 1000, 1, 1000, 1000, 1000, 1; + selfadjointeigensolver_essential_check(m); + } + + // Nearly singular: one eigenvalue much smaller than others. + { + Matrix3d m; + m << 1, 0.5, 0.3, 0.5, 1, 0.4, 0.3, 0.4, 1; + m *= 1e15; + Matrix3d perturbation = Matrix3d::Zero(); + perturbation(0, 0) = 1e-15; + m += perturbation; + selfadjointeigensolver_essential_check(m); + } +} + +// Specific 2x2 test cases that stress the direct solver. +template +void direct_2x2_stress() { + // Equal eigenvalues. + { + Matrix2d m = Matrix2d::Identity() * 42.0; + selfadjointeigensolver_essential_check(m); + } + + // Very small off-diagonal. + { + Matrix2d m; + m << 1.0, 1e-15, 1e-15, 1.0; + selfadjointeigensolver_essential_check(m); + } + + // Huge ratio between diagonal entries. + { + Matrix2d m; + m << 1e100, 0, 0, 1e-100; + selfadjointeigensolver_essential_check(m); + } + + // Anti-diagonal dominant. + { + Matrix2d m; + m << 0, 1e10, 1e10, 0; + selfadjointeigensolver_essential_check(m); + } + + // Negative entries. + { + Matrix2d m; + m << -5.0, 3.0, 3.0, -5.0; + selfadjointeigensolver_essential_check(m); + } +} + EIGEN_DECLARE_TEST(eigensolver_selfadjoint) { int s = 0; for (int i = 0; i < g_repeat; i++) { @@ -266,6 +685,39 @@ EIGEN_DECLARE_TEST(eigensolver_selfadjoint) { CALL_SUBTEST_5(selfadjointeigensolver(MatrixXcd(2, 2))); CALL_SUBTEST_6(selfadjointeigensolver(Matrix())); CALL_SUBTEST_7(selfadjointeigensolver(Matrix())); + + // repeated eigenvalues + CALL_SUBTEST_17(selfadjointeigensolver_repeated_eigenvalues(Matrix3d())); + CALL_SUBTEST_15(selfadjointeigensolver_repeated_eigenvalues(Matrix2d())); + CALL_SUBTEST_2(selfadjointeigensolver_repeated_eigenvalues(Matrix4d())); + CALL_SUBTEST_4(selfadjointeigensolver_repeated_eigenvalues(MatrixXd(s, s))); + CALL_SUBTEST_13(selfadjointeigensolver_repeated_eigenvalues(Matrix3f())); + CALL_SUBTEST_12(selfadjointeigensolver_repeated_eigenvalues(Matrix2f())); + CALL_SUBTEST_18(selfadjointeigensolver_repeated_eigenvalues(Matrix3cd())); + + // extreme eigenvalues (near overflow/underflow, high condition number) + CALL_SUBTEST_17(selfadjointeigensolver_extreme_eigenvalues(Matrix3d())); + CALL_SUBTEST_2(selfadjointeigensolver_extreme_eigenvalues(Matrix4d())); + CALL_SUBTEST_4(selfadjointeigensolver_extreme_eigenvalues(MatrixXd(s, s))); + CALL_SUBTEST_13(selfadjointeigensolver_extreme_eigenvalues(Matrix3f())); + CALL_SUBTEST_3(selfadjointeigensolver_extreme_eigenvalues(MatrixXf(s, s))); + + // computeFromTridiagonal with scaled inputs + CALL_SUBTEST_4(selfadjointeigensolver_tridiagonal_scaled(MatrixXd(s, s))); + CALL_SUBTEST_3(selfadjointeigensolver_tridiagonal_scaled(MatrixXf(s, s))); + + // diagonal matrices + CALL_SUBTEST_17(selfadjointeigensolver_diagonal(Matrix3d())); + CALL_SUBTEST_4(selfadjointeigensolver_diagonal(MatrixXd(s, s))); + + // operatorInverseSqrt + CALL_SUBTEST_17(selfadjointeigensolver_inverse_sqrt(Matrix3d())); + CALL_SUBTEST_2(selfadjointeigensolver_inverse_sqrt(Matrix4d())); + CALL_SUBTEST_4(selfadjointeigensolver_inverse_sqrt(MatrixXd(s, s))); + CALL_SUBTEST_13(selfadjointeigensolver_inverse_sqrt(Matrix3f())); + + // RowMajor + CALL_SUBTEST_19(selfadjointeigensolver_rowmajor<0>()); } CALL_SUBTEST_17(bug_854<0>()); @@ -273,6 +725,13 @@ EIGEN_DECLARE_TEST(eigensolver_selfadjoint) { CALL_SUBTEST_17(bug_1204<0>()); CALL_SUBTEST_17(bug_1225<0>()); + // Stress tests for direct 3x3 and 2x2 solvers. + CALL_SUBTEST_17(direct_3x3_stress<0>()); + CALL_SUBTEST_15(direct_2x2_stress<0>()); + + // Test Inf input handling. + CALL_SUBTEST_17(selfadjointeigensolver_inf<0>()); + // Test problem size constructors s = internal::random(1, EIGEN_TEST_MAX_SIZE / 4); CALL_SUBTEST_8(SelfAdjointEigenSolver tmp1(s));