GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-04-19 16:19:37 +08:00
Code Issues Packages Projects Releases Wiki Activity

bug #478: fix regression in the eigen decomposition of zero matrices.

Browse Source
This commit is contained in:
Gael Guennebaud 2017-01-31 14:22:42 +01:00
parent 63de19c000
commit 53026d29d4
5 changed files with 43 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
Eigen/src/Eigenvalues/ComplexEigenSolver.h
View File

@ -250,7 +250,7 @@ template<typename _MatrixType> class ComplexEigenSolver
EigenvectorType m_matX; EigenvectorType m_matX;
private: private:
void doComputeEigenvectors(const RealScalar& matrixnorm); void doComputeEigenvectors(RealScalar matrixnorm);
void sortEigenvalues(bool computeEigenvectors); void sortEigenvalues(bool computeEigenvectors);
}; };
@ -284,10 +284,12 @@ ComplexEigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool
template<typename MatrixType> template<typename MatrixType>
void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(const RealScalar& matrixnorm) void ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm)
{ {
const Index n = m_eivalues.size(); const Index n = m_eivalues.size();
matrixnorm = numext::maxi(matrixnorm,(std::numeric_limits<RealScalar>::min)());
// Compute X such that T = X D X^(-1), where D is the diagonal of T. // Compute X such that T = X D X^(-1), where D is the diagonal of T.
// The matrix X is unit triangular. // The matrix X is unit triangular.
m_matX = EigenvectorType::Zero(n, n); m_matX = EigenvectorType::Zero(n, n);

12
Eigen/src/Eigenvalues/RealSchur.h
View File

@ -248,12 +248,24 @@ template<typename MatrixType>
template<typename InputType> template<typename InputType>
RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU) RealSchur<MatrixType>& RealSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU)
{ {
const Scalar considerAsZero = (std::numeric_limits<Scalar>::min)();
eigen_assert(matrix.cols() == matrix.rows()); eigen_assert(matrix.cols() == matrix.rows());
Index maxIters = m_maxIters; Index maxIters = m_maxIters;
if (maxIters == -1) if (maxIters == -1)
maxIters = m_maxIterationsPerRow * matrix.rows(); maxIters = m_maxIterationsPerRow * matrix.rows();
Scalar scale = matrix.derived().cwiseAbs().maxCoeff(); Scalar scale = matrix.derived().cwiseAbs().maxCoeff();
if(scale<considerAsZero)
{
m_matT.setZero(matrix.rows(),matrix.cols());
if(computeU)
m_matU.setIdentity(matrix.rows(),matrix.cols());
m_info = Success;
m_isInitialized = true;
m_matUisUptodate = computeU;
return *this;
}
// Step 1. Reduce to Hessenberg form // Step 1. Reduce to Hessenberg form
m_hess.compute(matrix.derived()/scale); m_hess.compute(matrix.derived()/scale);

9
test/eigensolver_complex.cpp
View File

@ -131,6 +131,15 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
ComplexEigenSolver<MatrixType> eig(a.adjoint() * a); ComplexEigenSolver<MatrixType> eig(a.adjoint() * a);
eig.compute(a.adjoint() * a); eig.compute(a.adjoint() * a);
} }
// regression test for bug 478
{
a.setZero();
ComplexEigenSolver<MatrixType> ei3(a);
VERIFY_IS_EQUAL(ei3.info(), Success);
VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(),RealScalar(1));
VERIFY((ei3.eigenvectors().transpose()*ei3.eigenvectors().transpose()).eval().isIdentity());
}
} }
template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m) template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m)

9
test/eigensolver_generic.cpp
View File

@ -76,6 +76,15 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
EigenSolver<MatrixType> eig(a.adjoint() * a); EigenSolver<MatrixType> eig(a.adjoint() * a);
eig.compute(a.adjoint() * a); eig.compute(a.adjoint() * a);
} }
// regression test for bug 478
{
a.setZero();
EigenSolver<MatrixType> ei3(a);
VERIFY_IS_EQUAL(ei3.info(), Success);
VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(),RealScalar(1));
VERIFY((ei3.eigenvectors().transpose()*ei3.eigenvectors().transpose()).eval().isIdentity());
}
} }
template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m) template<typename MatrixType> void eigensolver_verify_assert(const MatrixType& m)

9
test/eigensolver_selfadjoint.cpp
View File

@ -180,6 +180,15 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
SelfAdjointEigenSolver<MatrixType> eig(a.adjoint() * a); SelfAdjointEigenSolver<MatrixType> eig(a.adjoint() * a);
eig.compute(a.adjoint() * a); eig.compute(a.adjoint() * a);
} }
// regression test for bug 478
{
a.setZero();
SelfAdjointEigenSolver<MatrixType> ei3(a);
VERIFY_IS_EQUAL(ei3.info(), Success);
VERIFY_IS_MUCH_SMALLER_THAN(ei3.eigenvalues().norm(),RealScalar(1));
VERIFY((ei3.eigenvectors().transpose()*ei3.eigenvectors().transpose()).eval().isIdentity());
}
} }
template<int> template<int>