diff options
| -rw-r--r-- | Eigen/src/Eigenvalues/ComplexEigenSolver.h | 6 | ||||
| -rw-r--r-- | Eigen/src/Eigenvalues/RealSchur.h | 12 | ||||
| -rw-r--r-- | test/eigensolver_complex.cpp | 9 | ||||
| -rw-r--r-- | test/eigensolver_generic.cpp | 9 | ||||
| -rw-r--r-- | test/eigensolver_selfadjoint.cpp | 9 |
5 files changed, 43 insertions, 2 deletions
diff --git a/Eigen/src/Eigenvalues/ComplexEigenSolver.h b/Eigen/src/Eigenvalues/ComplexEigenSolver.h index ec3b1633e..dc5fae06a 100644 --- a/Eigen/src/Eigenvalues/ComplexEigenSolver.h +++ b/Eigen/src/Eigenvalues/ComplexEigenSolver.h @@ -250,7 +250,7 @@ template<typename _MatrixType> class ComplexEigenSolver EigenvectorType m_matX; private: - void doComputeEigenvectors(const RealScalar& matrixnorm); + void doComputeEigenvectors(RealScalar matrixnorm); void sortEigenvalues(bool computeEigenvectors); }; @@ -284,10 +284,12 @@ ComplexEigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool template<typename MatrixType> -void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(const RealScalar& matrixnorm) +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); diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h index d6a339f07..f5c86041d 100644 --- a/Eigen/src/Eigenvalues/RealSchur.h +++ b/Eigen/src/Eigenvalues/RealSchur.h @@ -248,12 +248,24 @@ template<typename MatrixType> template<typename InputType> RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU) { + const Scalar considerAsZero = (std::numeric_limits<Scalar>::min)(); + eigen_assert(matrix.cols() == matrix.rows()); Index maxIters = m_maxIters; if (maxIters == -1) maxIters = m_maxIterationsPerRow * matrix.rows(); Scalar scale = matrix.derived().cwiseAbs().maxCoeff(); + if(scale<considerAsZero) + { + m_matT.setZero(matrix.rows(),matrix.cols()); + if(computeU) + m_matU.setIdentity(matrix.rows(),matrix.cols()); + m_info = Success; + m_isInitialized = true; + m_matUisUptodate = computeU; + return *this; + } // Step 1. Reduce to Hessenberg form m_hess.compute(matrix.derived()/scale); diff --git a/test/eigensolver_complex.cpp b/test/eigensolver_complex.cpp index 8e2bb9ef0..293b1b265 100644 --- a/test/eigensolver_complex.cpp +++ b/test/eigensolver_complex.cpp @@ -131,6 +131,15 @@ template<typename MatrixType> void eigensolver(const MatrixType& m) ComplexEigenSolver<MatrixType> eig(a.adjoint() * a); eig.compute(a.adjoint() * a); } + + // regression test for bug 478 + { + a.setZero(); + 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()); + } } template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m) diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp index e18fbf687..d0e644d4b 100644 --- a/test/eigensolver_generic.cpp +++ b/test/eigensolver_generic.cpp @@ -76,6 +76,15 @@ template<typename MatrixType> void eigensolver(const MatrixType& m) EigenSolver<MatrixType> eig(a.adjoint() * a); eig.compute(a.adjoint() * a); } + + // regression test for bug 478 + { + a.setZero(); + 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()); + } } template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m) diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp index 4ed126116..39ad4130e 100644 --- a/test/eigensolver_selfadjoint.cpp +++ b/test/eigensolver_selfadjoint.cpp @@ -180,6 +180,15 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m) SelfAdjointEigenSolver<MatrixType> eig(a.adjoint() * a); eig.compute(a.adjoint() * a); } + + // regression test for bug 478 + { + a.setZero(); + 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()); + } } template<int> |