GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-06-04 18:54:00 +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 c2107d30ce
commit 4a936974a5
1 changed files with 19 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -739,6 +739,7 @@ struct direct_selfadjoint_eigenvalues<SolverType,2,false>
static inline void run(SolverType& solver, const MatrixType& mat, int options) static inline void run(SolverType& solver, const MatrixType& mat, int options)
{ {
EIGEN_USING_STD_MATH(sqrt); EIGEN_USING_STD_MATH(sqrt);
EIGEN_USING_STD_MATH(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
@ -760,22 +761,29 @@ struct direct_selfadjoint_eigenvalues<SolverType,2,false>
// 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.