diff options
author | Gael Guennebaud <g.gael@free.fr> | 2015-06-08 16:16:42 +0200 |
---|---|---|
committer | Gael Guennebaud <g.gael@free.fr> | 2015-06-08 16:16:42 +0200 |
commit | cd8b996f99de67035a0504cbaf0a627fb68f0f1d (patch) | |
tree | 7321a99c33157a2a696c7b4f5c2efc7d4d7dd6b7 | |
parent | 913a61870da36dab308b38c48cb3553487dca8ff (diff) |
Extend unit test and documentation of SelfAdjointEigenSolver::computeDirect
-rw-r--r-- | Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 10 | ||||
-rw-r--r-- | test/eigensolver_selfadjoint.cpp | 32 |
2 files changed, 29 insertions, 13 deletions
diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 1830b99c1..27a014a96 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -198,17 +198,21 @@ template<typename _MatrixType> class SelfAdjointEigenSolver EIGEN_DEVICE_FUNC SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors); - /** \brief Computes eigendecomposition of given matrix using a direct algorithm + /** \brief Computes eigendecomposition of given matrix using a closed-form algorithm * * This is a variant of compute(const MatrixType&, int options) which * directly solves the underlying polynomial equation. * - * Currently only 3x3 matrices for which the sizes are known at compile time are supported (e.g., Matrix3d). + * Currently only 2x2 and 3x3 matrices for which the sizes are known at compile time are supported (e.g., Matrix3d). * - * This method is usually significantly faster than the QR algorithm + * This method is usually significantly faster than the QR iterative algorithm * but it might also be less accurate. It is also worth noting that * for 3x3 matrices it involves trigonometric operations which are * not necessarily available for all scalar types. + * + * For the 3x3 case, we observed the following worst case relative error regarding the eigenvalues: + * - double: 1e-8 + * - float: 1e-3 * * \sa compute(const MatrixType&, int options) */ diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp index 6750c7609..41b6d99ab 100644 --- a/test/eigensolver_selfadjoint.cpp +++ b/test/eigensolver_selfadjoint.cpp @@ -18,22 +18,34 @@ template<typename MatrixType> void selfadjointeigensolver_essential_check(const { typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits<Scalar>::Real RealScalar; - RealScalar largerEps = 10*test_precision<RealScalar>(); + RealScalar eival_eps = (std::min)(test_precision<RealScalar>(), NumTraits<Scalar>::dummy_precision()*20000); SelfAdjointEigenSolver<MatrixType> eiSymm(m); VERIFY_IS_EQUAL(eiSymm.info(), Success); - VERIFY((m.template selfadjointView<Lower>() * eiSymm.eigenvectors()).isApprox( - eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps)); + VERIFY_IS_APPROX(m.template selfadjointView<Lower>() * eiSymm.eigenvectors(), + eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal()); VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues(), eiSymm.eigenvalues()); VERIFY_IS_UNITARY(eiSymm.eigenvectors()); - SelfAdjointEigenSolver<MatrixType> eiDirect; - eiDirect.computeDirect(m); - VERIFY_IS_EQUAL(eiDirect.info(), Success); - VERIFY((m.template selfadjointView<Lower>() * eiDirect.eigenvectors()).isApprox( - eiDirect.eigenvectors() * eiDirect.eigenvalues().asDiagonal(), largerEps)); - VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues(), eiDirect.eigenvalues()); - VERIFY_IS_UNITARY(eiDirect.eigenvectors()); + if(m.cols()<=4) + { + SelfAdjointEigenSolver<MatrixType> eiDirect; + eiDirect.computeDirect(m); + VERIFY_IS_EQUAL(eiDirect.info(), Success); + VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiDirect.eigenvalues()); + if(! eiSymm.eigenvalues().isApprox(eiDirect.eigenvalues(), eival_eps) ) + { + std::cerr << "reference eigenvalues: " << eiSymm.eigenvalues().transpose() << "\n" + << "obtained eigenvalues: " << eiDirect.eigenvalues().transpose() << "\n" + << "diff: " << (eiSymm.eigenvalues()-eiDirect.eigenvalues()).transpose() << "\n" + << "error (eps): " << (eiSymm.eigenvalues()-eiDirect.eigenvalues()).norm() / eiSymm.eigenvalues().norm() << " (" << eival_eps << ")\n"; + } + VERIFY(eiSymm.eigenvalues().isApprox(eiDirect.eigenvalues(), eival_eps)); + VERIFY_IS_APPROX(m.template selfadjointView<Lower>() * eiDirect.eigenvectors(), + eiDirect.eigenvectors() * eiDirect.eigenvalues().asDiagonal()); + VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues(), eiDirect.eigenvalues()); + VERIFY_IS_UNITARY(eiDirect.eigenvectors()); + } } template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m) |