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

Tridiagonalization::diagonal() and ::subDiagonal() did not work. Added unit-test

Browse Source
This commit is contained in:
Christoph Hertzberg 2014-09-24 14:37:13 +02:00
parent 446001ef51
commit 27d6b4daf9
2 changed files with 14 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
Eigen/src/Eigenvalues/Tridiagonalization.h
View File

@ -91,10 +91,8 @@ template<typename _MatrixType> class Tridiagonalization
>::type DiagonalReturnType; >::type DiagonalReturnType;
typedef typename internal::conditional<NumTraits<Scalar>::IsComplex, typedef typename internal::conditional<NumTraits<Scalar>::IsComplex,
typename internal::add_const_on_value_type<typename Diagonal< typename internal::add_const_on_value_type<typename Diagonal<const MatrixType, -1>::RealReturnType>::type,
Block<const MatrixType,SizeMinusOne,SizeMinusOne> >::RealReturnType>::type, const Diagonal<const MatrixType, -1>
const Diagonal<
Block<const MatrixType,SizeMinusOne,SizeMinusOne> >
>::type SubDiagonalReturnType; >::type SubDiagonalReturnType;
/** \brief Return type of matrixQ() */ /** \brief Return type of matrixQ() */
@ -307,7 +305,7 @@ typename Tridiagonalization<MatrixType>::DiagonalReturnType
Tridiagonalization<MatrixType>::diagonal() const Tridiagonalization<MatrixType>::diagonal() const
{ {
eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
return m_matrix.diagonal(); return m_matrix.diagonal().real();
} }
template<typename MatrixType> template<typename MatrixType>
@ -315,8 +313,7 @@ typename Tridiagonalization<MatrixType>::SubDiagonalReturnType
Tridiagonalization<MatrixType>::subDiagonal() const Tridiagonalization<MatrixType>::subDiagonal() const
{ {
eigen_assert(m_isInitialized && "Tridiagonalization is not initialized."); eigen_assert(m_isInitialized && "Tridiagonalization is not initialized.");
Index n = m_matrix.rows(); return m_matrix.template diagonal<-1>().real();
return Block<const MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1).diagonal();
} }
namespace internal { namespace internal {

11
test/eigensolver_selfadjoint.cpp
View File

@ -111,8 +111,17 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
// test Tridiagonalization's methods // test Tridiagonalization's methods
Tridiagonalization<MatrixType> tridiag(symmC); Tridiagonalization<MatrixType> tridiag(symmC);
// FIXME tridiag.matrixQ().adjoint() does not work VERIFY_IS_APPROX(tridiag.diagonal(), tridiag.matrixT().template diagonal());
VERIFY_IS_APPROX(tridiag.subDiagonal(), tridiag.matrixT().template diagonal<-1>());
MatrixType T = tridiag.matrixT();
if(rows>1 && cols>1) {
// FIXME check that upper and lower part are 0:
//VERIFY(T.topRightCorner(rows-2, cols-2).template triangularView<Upper>().isZero());
}
VERIFY_IS_APPROX(tridiag.diagonal(), T.diagonal().real());
VERIFY_IS_APPROX(tridiag.subDiagonal(), T.template diagonal<1>().real());
VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint()); VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT().eval() * MatrixType(tridiag.matrixQ()).adjoint());
VERIFY_IS_APPROX(MatrixType(symmC.template selfadjointView<Lower>()), tridiag.matrixQ() * tridiag.matrixT() * tridiag.matrixQ().adjoint());
// Test computation of eigenvalues from tridiagonal matrix // Test computation of eigenvalues from tridiagonal matrix
if(rows > 1) if(rows > 1)