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

Extend unit test and documentation of SelfAdjointEigenSolver::computeDirect

Browse Source
This commit is contained in:
Gael Guennebaud 2015-06-08 16:16:42 +02:00
parent 913a61870d
commit cd8b996f99
2 changed files with 29 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -198,18 +198,22 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
EIGEN_DEVICE_FUNC EIGEN_DEVICE_FUNC
SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors); 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 * This is a variant of compute(const MatrixType&, int options) which
* directly solves the underlying polynomial equation. * 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 * but it might also be less accurate. It is also worth noting that
* for 3x3 matrices it involves trigonometric operations which are * for 3x3 matrices it involves trigonometric operations which are
* not necessarily available for all scalar types. * 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) * \sa compute(const MatrixType&, int options)
*/ */
EIGEN_DEVICE_FUNC EIGEN_DEVICE_FUNC

22
test/eigensolver_selfadjoint.cpp
View File

@ -18,23 +18,35 @@ template<typename MatrixType> void selfadjointeigensolver_essential_check(const
{ {
typedef typename MatrixType::Scalar Scalar; typedef typename MatrixType::Scalar Scalar;
typedef typename NumTraits<Scalar>::Real RealScalar; 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); SelfAdjointEigenSolver<MatrixType> eiSymm(m);
VERIFY_IS_EQUAL(eiSymm.info(), Success); VERIFY_IS_EQUAL(eiSymm.info(), Success);
VERIFY((m.template selfadjointView<Lower>() * eiSymm.eigenvectors()).isApprox( VERIFY_IS_APPROX(m.template selfadjointView<Lower>() * eiSymm.eigenvectors(),
eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps)); eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal());
VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues(), eiSymm.eigenvalues()); VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues(), eiSymm.eigenvalues());
VERIFY_IS_UNITARY(eiSymm.eigenvectors()); VERIFY_IS_UNITARY(eiSymm.eigenvectors());
if(m.cols()<=4)
{
SelfAdjointEigenSolver<MatrixType> eiDirect; SelfAdjointEigenSolver<MatrixType> eiDirect;
eiDirect.computeDirect(m); eiDirect.computeDirect(m);
VERIFY_IS_EQUAL(eiDirect.info(), Success); VERIFY_IS_EQUAL(eiDirect.info(), Success);
VERIFY((m.template selfadjointView<Lower>() * eiDirect.eigenvectors()).isApprox( VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiDirect.eigenvalues());
eiDirect.eigenvectors() * eiDirect.eigenvalues().asDiagonal(), largerEps)); 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_APPROX(m.template selfadjointView<Lower>().eigenvalues(), eiDirect.eigenvalues());
VERIFY_IS_UNITARY(eiDirect.eigenvectors()); VERIFY_IS_UNITARY(eiDirect.eigenvectors());
} }
}
template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m) template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
{ {