devtools/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2026-04-10 11:34:33 +08:00
Code Issues Packages Projects Releases Wiki Activity

Fix ComplexEigenSolver NaN with flush-to-zero arithmetic

Browse Source 
libeigen/eigen!2196

Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
This commit is contained in:
Rasmus Munk Larsen
2026-02-23 11:15:31 -08:00
parent 667cabe3aa
commit d537b51ede
4 changed files with 8 additions and 6 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
Eigen/src/Eigenvalues/ComplexEigenSolver.h
View File

@@ -265,8 +265,6 @@ template <typename MatrixType>
void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm) {
const Index n = m_eivalues.size();
matrixnorm = numext::maxi(matrixnorm, (std::numeric_limits<RealScalar>::min)());
// Compute X such that T = X D X^(-1), where D is the diagonal of T.
// The matrix X is unit triangular.
m_matX = EigenvectorType::Zero(n, n);
@@ -282,7 +280,8 @@ void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm
if (z == ComplexScalar(0)) {
// If the i-th and k-th eigenvalue are equal, then z equals 0.
// Use a small value instead, to prevent division by zero.
numext::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm;
numext::real_ref(z) = numext::maxi(std::numeric_limits<RealScalar>::epsilon() * matrixnorm,
(std::numeric_limits<RealScalar>::min)());
}
m_matX.coeffRef(i, k) = m_matX.coeff(i, k) / z;
}

3
test/eigensolver_complex.cpp
View File

@@ -132,7 +132,8 @@ void eigensolver(const MatrixType& m) {
ComplexEigenSolver<MatrixType> ei3(a);
VERIFY_IS_EQUAL(ei3.info(), Success);
VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(), RealScalar(1));
VERIFY((ei3.eigenvectors().transpose() * ei3.eigenvectors().transpose()).eval().isIdentity());
RealScalar tol = 2 * a.cols() * NumTraits<RealScalar>::epsilon();
VERIFY((ei3.eigenvectors().adjoint() * ei3.eigenvectors()).eval().isIdentity(tol));
}
}

3
test/eigensolver_generic.cpp
View File

@@ -89,7 +89,8 @@ void eigensolver(const MatrixType& m) {
EigenSolver<MatrixType> ei3(a);
VERIFY_IS_EQUAL(ei3.info(), Success);
VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(), RealScalar(1));
VERIFY((ei3.eigenvectors().transpose() * ei3.eigenvectors().transpose()).eval().isIdentity());
RealScalar tol = 2 * a.cols() * NumTraits<RealScalar>::epsilon();
VERIFY((ei3.eigenvectors().adjoint() * ei3.eigenvectors()).eval().isIdentity(tol));
}
}

3
test/eigensolver_selfadjoint.cpp
View File

@@ -197,7 +197,8 @@ void selfadjointeigensolver(const MatrixType& m) {
SelfAdjointEigenSolver<MatrixType> ei3(a);
VERIFY_IS_EQUAL(ei3.info(), Success);
VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(), RealScalar(1));
VERIFY((ei3.eigenvectors().transpose() * ei3.eigenvectors().transpose()).eval().isIdentity());
RealScalar tol = 2 * a.cols() * NumTraits<RealScalar>::epsilon();
VERIFY((ei3.eigenvectors().adjoint() * ei3.eigenvectors()).eval().isIdentity(tol));
}
}