diff options
| -rw-r--r-- | Eigen/src/Eigenvalues/EigenSolver.h | 15 | ||||
| -rw-r--r-- | test/eigensolver_generic.cpp | 13 |
2 files changed, 27 insertions, 1 deletions
diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h index bf20e03ef..739466949 100644 --- a/Eigen/src/Eigenvalues/EigenSolver.h +++ b/Eigen/src/Eigenvalues/EigenSolver.h @@ -366,6 +366,7 @@ EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvect { using std::sqrt; using std::abs; + using std::max; eigen_assert(matrix.cols() == matrix.rows()); // Reduce to real Schur form. @@ -390,7 +391,19 @@ EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvect else { Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1)); - Scalar z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1))); + Scalar z; + // Compute z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1))); + // without overflow + { + Scalar t0 = m_matT.coeff(i+1, i); + Scalar t1 = m_matT.coeff(i, i+1); + Scalar maxval = (max)(abs(p),(max)(abs(t0),abs(t1))); + t0 /= maxval; + t1 /= maxval; + Scalar p0 = p/maxval; + z = maxval * sqrt(abs(p0 * p0 + t0 * t1)); + } + m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z); m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z); i += 2; diff --git a/test/eigensolver_generic.cpp b/test/eigensolver_generic.cpp index 005af81eb..91383b5cf 100644 --- a/test/eigensolver_generic.cpp +++ b/test/eigensolver_generic.cpp @@ -121,5 +121,18 @@ void test_eigensolver_generic() } ); + // regression test for bug 793 +#ifdef EIGEN_TEST_PART_2 + { + MatrixXd a(3,3); + a << 0, 0, 1, + 1, 1, 1, + 1, 1e+200, 1; + Eigen::EigenSolver<MatrixXd> eig(a); + VERIFY_IS_APPROX(a * eig.pseudoEigenvectors(), eig.pseudoEigenvectors() * eig.pseudoEigenvalueMatrix()); + VERIFY_IS_APPROX(a * eig.eigenvectors(), eig.eigenvectors() * eig.eigenvalues().asDiagonal()); + } +#endif + TEST_SET_BUT_UNUSED_VARIABLE(s) } |