GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-07-14 17:11:50 +08:00
Code Issues Packages Projects Releases Wiki Activity

bug #1013: fix 2x2 direct eigensolver for identical eiegnvalues

Browse Source
This commit is contained in:
Gael Guennebaud 2015-05-07 15:55:12 +02:00
parent 62d334c7d3
commit 847bb317cd
1 changed files with 20 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

31
Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
View File

@ -695,6 +695,8 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2
static inline void run(SolverType& solver, const MatrixType& mat, int options) static inline void run(SolverType& solver, const MatrixType& mat, int options)
{ {
using std::sqrt; using std::sqrt;
using std::abs;
eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows()); eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
eigen_assert((options&~(EigVecMask|GenEigMask))==0 eigen_assert((options&~(EigVecMask|GenEigMask))==0
&& (options&EigVecMask)!=EigVecMask && (options&EigVecMask)!=EigVecMask
@ -715,22 +717,29 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2
// compute the eigen vectors // compute the eigen vectors
if(computeEigenvectors) if(computeEigenvectors)
{ {
scaledMat.diagonal().array () -= eivals(1); if((eivals(1)-eivals(0))<=abs(eivals(1))*Eigen::NumTraits<Scalar>::epsilon())
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.setIdentity();
eivecs.col(1) /= sqrt(a2+b2);
} }
else else
{ {
eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0); scaledMat.diagonal().array () -= eivals(1);
eivecs.col(1) /= sqrt(c2+b2); 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. // Rescale back to the original size.