diff options
author | Gael Guennebaud <g.gael@free.fr> | 2014-08-21 10:49:09 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2014-08-21 10:49:09 +0200 |
commit | 9c0aa81fbfc4abd0dc297d75305b9d579dad2754 (patch) | |
tree | cc132ae596a596a652c53d189364d79a77eb9abf /Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | |
parent | eeadc06e838a737492625cc3e4d5d5555bb40ff7 (diff) |
bug #854: fix numerical issue in SelfAdjointEigenSolver::computeDirect for 3x3 matrices. The tolerance to detect stable cross products was too optimistic.
Add respective unit tests.
Diffstat (limited to 'Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h')
-rw-r--r-- | Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 18 |
1 files changed, 14 insertions, 4 deletions
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index fc8ecaa6f..a6bbdac6b 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -605,7 +605,6 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 if(computeEigenvectors) { Scalar safeNorm2 = Eigen::NumTraits<Scalar>::epsilon(); - safeNorm2 *= safeNorm2; if((eivals(2)-eivals(0))<=Eigen::NumTraits<Scalar>::epsilon()) { eivecs.setIdentity(); @@ -619,7 +618,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 Scalar d0 = eivals(2) - eivals(1); Scalar d1 = eivals(1) - eivals(0); int k = d0 > d1 ? 2 : 0; - d0 = d0 > d1 ? d1 : d0; + d0 = d0 > d1 ? d0 : d1; tmp.diagonal().array () -= eivals(k); VectorType cross; @@ -627,19 +626,25 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 n = (cross = tmp.row(0).cross(tmp.row(1))).squaredNorm(); if(n>safeNorm2) + { eivecs.col(k) = cross / sqrt(n); + } else { n = (cross = tmp.row(0).cross(tmp.row(2))).squaredNorm(); if(n>safeNorm2) + { eivecs.col(k) = cross / sqrt(n); + } else { n = (cross = tmp.row(1).cross(tmp.row(2))).squaredNorm(); if(n>safeNorm2) + { eivecs.col(k) = cross / sqrt(n); + } else { // the input matrix and/or the eigenvaues probably contains some inf/NaN, @@ -659,12 +664,16 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 tmp.diagonal().array() -= eivals(1); if(d0<=Eigen::NumTraits<Scalar>::epsilon()) + { eivecs.col(1) = eivecs.col(k).unitOrthogonal(); + } else { - n = (cross = eivecs.col(k).cross(tmp.row(0).normalized())).squaredNorm(); + n = (cross = eivecs.col(k).cross(tmp.row(0))).squaredNorm(); if(n>safeNorm2) + { eivecs.col(1) = cross / sqrt(n); + } else { n = (cross = eivecs.col(k).cross(tmp.row(1))).squaredNorm(); @@ -678,13 +687,14 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3 else { // we should never reach this point, - // if so the last two eigenvalues are likely to ve very closed to each other + // if so the last two eigenvalues are likely to be very close to each other eivecs.col(1) = eivecs.col(k).unitOrthogonal(); } } } // make sure that eivecs[1] is orthogonal to eivecs[2] + // FIXME: this step should not be needed Scalar d = eivecs.col(1).dot(eivecs.col(k)); eivecs.col(1) = (eivecs.col(1) - d * eivecs.col(k)).normalized(); } |