diff options
author | Gael Guennebaud <g.gael@free.fr> | 2015-05-07 15:55:12 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2015-05-07 15:55:12 +0200 |
commit | 4a936974a5b895f61d0f80a61b27ea642037cb11 (patch) | |
tree | 311f4b521eecb64b95d54494b64d07ae2561d943 /Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | |
parent | c2107d30ce9b9f30ff1f2d436667f3d09a4d9bd5 (diff) |
bug #1013: fix 2x2 direct eigensolver for identical eiegnvalues
Diffstat (limited to 'Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h')
-rw-r--r-- | Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 30 |
1 files changed, 19 insertions, 11 deletions
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 1dcfacf0b..8240b9465 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -739,6 +739,7 @@ struct direct_selfadjoint_eigenvalues<SolverType,2,false> static inline void run(SolverType& solver, const MatrixType& mat, int options) { EIGEN_USING_STD_MATH(sqrt); + EIGEN_USING_STD_MATH(abs); eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows()); eigen_assert((options&~(EigVecMask|GenEigMask))==0 @@ -760,22 +761,29 @@ struct direct_selfadjoint_eigenvalues<SolverType,2,false> // compute the eigen vectors if(computeEigenvectors) { - scaledMat.diagonal().array () -= eivals(1); - Scalar a2 = numext::abs2(scaledMat(0,0)); - Scalar c2 = numext::abs2(scaledMat(1,1)); - Scalar b2 = numext::abs2(scaledMat(1,0)); - if(a2>c2) + if((eivals(1)-eivals(0))<=abs(eivals(1))*Eigen::NumTraits<Scalar>::epsilon()) { - eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0); - eivecs.col(1) /= sqrt(a2+b2); + eivecs.setIdentity(); } else { - eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0); - eivecs.col(1) /= sqrt(c2+b2); - } + scaledMat.diagonal().array () -= eivals(1); + Scalar a2 = numext::abs2(scaledMat(0,0)); + Scalar c2 = numext::abs2(scaledMat(1,1)); + Scalar b2 = numext::abs2(scaledMat(1,0)); + if(a2>c2) + { + eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0); + eivecs.col(1) /= sqrt(a2+b2); + } + else + { + eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0); + eivecs.col(1) /= sqrt(c2+b2); + } - eivecs.col(0) << eivecs.col(1).unitOrthogonal(); + eivecs.col(0) << eivecs.col(1).unitOrthogonal(); + } } // Rescale back to the original size. |