From db8631b66a4f5cb9957a58cc629be5ef438e3059 Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Sun, 30 May 2010 21:49:35 +0100 Subject: [PATCH] Guard with assert against using decomposition objects uninitialized. --- .../src/Eigenvalues/HessenbergDecomposition.h | 29 +++++++++-- .../src/Eigenvalues/SelfAdjointEigenSolver.h | 51 ++++++++++++------- Eigen/src/Eigenvalues/Tridiagonalization.h | 27 ++++++++-- test/eigensolver_complex.cpp | 12 +++++ test/eigensolver_selfadjoint.cpp | 14 +++++ test/hessenberg.cpp | 7 +++ 6 files changed, 113 insertions(+), 27 deletions(-) diff --git a/Eigen/src/Eigenvalues/HessenbergDecomposition.h b/Eigen/src/Eigenvalues/HessenbergDecomposition.h index 1111ffb12..220531bf5 100644 --- a/Eigen/src/Eigenvalues/HessenbergDecomposition.h +++ b/Eigen/src/Eigenvalues/HessenbergDecomposition.h @@ -107,7 +107,8 @@ template class HessenbergDecomposition */ HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size) : m_matrix(size,size), - m_temp(size) + m_temp(size), + m_isInitialized(false) { if(size>1) m_hCoeffs.resize(size-1); @@ -124,12 +125,17 @@ template class HessenbergDecomposition */ HessenbergDecomposition(const MatrixType& matrix) : m_matrix(matrix), - m_temp(matrix.rows()) + m_temp(matrix.rows()), + m_isInitialized(false) { if(matrix.rows()<2) + { + m_isInitialized = true; return; + } m_hCoeffs.resize(matrix.rows()-1,1); _compute(m_matrix, m_hCoeffs, m_temp); + m_isInitialized = true; } /** \brief Computes Hessenberg decomposition of given matrix. @@ -152,9 +158,13 @@ template class HessenbergDecomposition { m_matrix = matrix; if(matrix.rows()<2) + { + m_isInitialized = true; return; + } m_hCoeffs.resize(matrix.rows()-1,1); _compute(m_matrix, m_hCoeffs, m_temp); + m_isInitialized = true; } /** \brief Returns the Householder coefficients. @@ -170,7 +180,11 @@ template class HessenbergDecomposition * * \sa packedMatrix(), \ref Householder_Module "Householder module" */ - const CoeffVectorType& householderCoefficients() const { return m_hCoeffs; } + const CoeffVectorType& householderCoefficients() const + { + ei_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); + return m_hCoeffs; + } /** \brief Returns the internal representation of the decomposition * @@ -201,7 +215,11 @@ template class HessenbergDecomposition * * \sa householderCoefficients() */ - const MatrixType& packedMatrix() const { return m_matrix; } + const MatrixType& packedMatrix() const + { + ei_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); + return m_matrix; + } /** \brief Reconstructs the orthogonal matrix Q in the decomposition * @@ -219,6 +237,7 @@ template class HessenbergDecomposition */ HouseholderSequenceType matrixQ() const { + ei_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate(), false, m_matrix.rows() - 1, 1); } @@ -244,6 +263,7 @@ template class HessenbergDecomposition */ HessenbergDecompositionMatrixHReturnType matrixH() const { + ei_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); return HessenbergDecompositionMatrixHReturnType(*this); } @@ -257,6 +277,7 @@ template class HessenbergDecomposition MatrixType m_matrix; CoeffVectorType m_hCoeffs; VectorType m_temp; + bool m_isInitialized; }; #ifndef EIGEN_HIDE_HEAVY_CODE diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h index 2c53655d1..20773a549 100644 --- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h +++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h @@ -115,10 +115,9 @@ template class SelfAdjointEigenSolver : m_eivec(), m_eivalues(), m_tridiag(), - m_subdiag() - { - ei_assert(Size!=Dynamic); - } + m_subdiag(), + m_isInitialized(false) + { } /** \brief Constructor, pre-allocates memory for dynamic-size matrices. * @@ -137,7 +136,8 @@ template class SelfAdjointEigenSolver : m_eivec(size, size), m_eivalues(size), m_tridiag(size), - m_subdiag(size > 1 ? size - 1 : 1) + m_subdiag(size > 1 ? size - 1 : 1), + m_isInitialized(false) {} /** \brief Constructor; computes eigendecomposition of given matrix. @@ -162,7 +162,8 @@ template class SelfAdjointEigenSolver : m_eivec(matrix.rows(), matrix.cols()), m_eivalues(matrix.cols()), m_tridiag(matrix.rows()), - m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1) + m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1), + m_isInitialized(false) { compute(matrix, computeEigenvectors); } @@ -192,7 +193,8 @@ template class SelfAdjointEigenSolver : m_eivec(matA.rows(), matA.cols()), m_eivalues(matA.cols()), m_tridiag(matA.rows()), - m_subdiag(matA.rows() > 1 ? matA.rows() - 1 : 1) + m_subdiag(matA.rows() > 1 ? matA.rows() - 1 : 1), + m_isInitialized(false) { compute(matA, matB, computeEigenvectors); } @@ -283,9 +285,8 @@ template class SelfAdjointEigenSolver */ const MatrixType& eigenvectors() const { - #ifndef NDEBUG - ei_assert(m_eigenvectorsOk); - #endif + ei_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); + ei_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); return m_eivec; } @@ -303,7 +304,11 @@ template class SelfAdjointEigenSolver * * \sa eigenvectors(), MatrixBase::eigenvalues() */ - const RealVectorType& eigenvalues() const { return m_eivalues; } + const RealVectorType& eigenvalues() const + { + ei_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); + return m_eivalues; + } /** \brief Computes the positive-definite square root of the matrix. * @@ -325,6 +330,8 @@ template class SelfAdjointEigenSolver */ MatrixType operatorSqrt() const { + ei_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); + ei_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint(); } @@ -348,6 +355,8 @@ template class SelfAdjointEigenSolver */ MatrixType operatorInverseSqrt() const { + ei_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); + ei_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint(); } @@ -357,9 +366,8 @@ template class SelfAdjointEigenSolver RealVectorType m_eivalues; TridiagonalizationType m_tridiag; typename TridiagonalizationType::SubDiagonalType m_subdiag; - #ifndef NDEBUG + bool m_isInitialized; bool m_eigenvectorsOk; - #endif }; #ifndef EIGEN_HIDE_HEAVY_CODE @@ -386,18 +394,19 @@ static void ei_tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index template SelfAdjointEigenSolver& SelfAdjointEigenSolver::compute(const MatrixType& matrix, bool computeEigenvectors) { - #ifndef NDEBUG - m_eigenvectorsOk = computeEigenvectors; - #endif assert(matrix.cols() == matrix.rows()); Index n = matrix.cols(); m_eivalues.resize(n,1); - m_eivec.resize(n,n); + if(computeEigenvectors) + m_eivec.resize(n,n); if(n==1) { m_eivalues.coeffRef(0,0) = ei_real(matrix.coeff(0,0)); - m_eivec.setOnes(); + if(computeEigenvectors) + m_eivec.setOnes(); + m_isInitialized = true; + m_eigenvectorsOk = computeEigenvectors; return *this; } @@ -438,9 +447,13 @@ SelfAdjointEigenSolver& SelfAdjointEigenSolver::compute( if (k > 0) { std::swap(m_eivalues[i], m_eivalues[k+i]); - m_eivec.col(i).swap(m_eivec.col(k+i)); + if(computeEigenvectors) + m_eivec.col(i).swap(m_eivec.col(k+i)); } } + + m_isInitialized = true; + m_eigenvectorsOk = computeEigenvectors; return *this; } diff --git a/Eigen/src/Eigenvalues/Tridiagonalization.h b/Eigen/src/Eigenvalues/Tridiagonalization.h index 02917f2e6..e204c30ea 100644 --- a/Eigen/src/Eigenvalues/Tridiagonalization.h +++ b/Eigen/src/Eigenvalues/Tridiagonalization.h @@ -109,7 +109,9 @@ template class Tridiagonalization * \sa compute() for an example. */ Tridiagonalization(Index size = Size==Dynamic ? 2 : Size) - : m_matrix(size,size), m_hCoeffs(size > 1 ? size-1 : 1) + : m_matrix(size,size), + m_hCoeffs(size > 1 ? size-1 : 1), + m_isInitialized(false) {} /** \brief Constructor; computes tridiagonal decomposition of given matrix. @@ -123,9 +125,12 @@ template class Tridiagonalization * Output: \verbinclude Tridiagonalization_Tridiagonalization_MatrixType.out */ Tridiagonalization(const MatrixType& matrix) - : m_matrix(matrix), m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1) + : m_matrix(matrix), + m_hCoeffs(matrix.cols() > 1 ? matrix.cols()-1 : 1), + m_isInitialized(false) { _compute(m_matrix, m_hCoeffs); + m_isInitialized = true; } /** \brief Computes tridiagonal decomposition of given matrix. @@ -149,6 +154,7 @@ template class Tridiagonalization m_matrix = matrix; m_hCoeffs.resize(matrix.rows()-1, 1); _compute(m_matrix, m_hCoeffs); + m_isInitialized = true; } /** \brief Returns the Householder coefficients. @@ -167,7 +173,11 @@ template class Tridiagonalization * * \sa packedMatrix(), \ref Householder_Module "Householder module" */ - inline CoeffVectorType householderCoefficients() const { return m_hCoeffs; } + inline CoeffVectorType householderCoefficients() const + { + ei_assert(m_isInitialized && "Tridiagonalization is not initialized."); + return m_hCoeffs; + } /** \brief Returns the internal representation of the decomposition * @@ -200,7 +210,11 @@ template class Tridiagonalization * * \sa householderCoefficients() */ - inline const MatrixType& packedMatrix() const { return m_matrix; } + inline const MatrixType& packedMatrix() const + { + ei_assert(m_isInitialized && "Tridiagonalization is not initialized."); + return m_matrix; + } /** \brief Returns the unitary matrix Q in the decomposition * @@ -219,6 +233,7 @@ template class Tridiagonalization */ HouseholderSequenceType matrixQ() const { + ei_assert(m_isInitialized && "Tridiagonalization is not initialized."); return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate(), false, m_matrix.rows() - 1, 1); } @@ -312,12 +327,14 @@ template class Tridiagonalization MatrixType m_matrix; CoeffVectorType m_hCoeffs; + bool m_isInitialized; }; template const typename Tridiagonalization::DiagonalReturnType Tridiagonalization::diagonal() const { + ei_assert(m_isInitialized && "Tridiagonalization is not initialized."); return m_matrix.diagonal(); } @@ -325,6 +342,7 @@ template const typename Tridiagonalization::SubDiagonalReturnType Tridiagonalization::subDiagonal() const { + ei_assert(m_isInitialized && "Tridiagonalization is not initialized."); Index n = m_matrix.rows(); return Block(m_matrix, 1, 0, n-1,n-1).diagonal(); } @@ -335,6 +353,7 @@ Tridiagonalization::matrixT() const { // FIXME should this function (and other similar ones) rather take a matrix as argument // and fill it ? (to avoid temporaries) + ei_assert(m_isInitialized && "Tridiagonalization is not initialized."); Index n = m_matrix.rows(); MatrixType matT = m_matrix; matT.topRightCorner(n-1, n-1).diagonal() = subDiagonal().template cast().conjugate(); diff --git a/test/eigensolver_complex.cpp b/test/eigensolver_complex.cpp index 5c5d7b38f..4437d9811 100644 --- a/test/eigensolver_complex.cpp +++ b/test/eigensolver_complex.cpp @@ -76,6 +76,13 @@ template void eigensolver(const MatrixType& m) VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1)); } +template void eigensolver_verify_assert() +{ + ComplexEigenSolver eig; + VERIFY_RAISES_ASSERT(eig.eigenvectors()) + VERIFY_RAISES_ASSERT(eig.eigenvalues()) +} + void test_eigensolver_complex() { for(int i = 0; i < g_repeat; i++) { @@ -85,6 +92,11 @@ void test_eigensolver_complex() CALL_SUBTEST_4( eigensolver(Matrix3f()) ); } + CALL_SUBTEST_1( eigensolver_verify_assert(Matrix4cf()) ); + CALL_SUBTEST_2( eigensolver_verify_assert(MatrixXcd(14,14)) ); + CALL_SUBTEST_3( eigensolver_verify_assert(Matrix, 1, 1>()) ); + CALL_SUBTEST_4( eigensolver_verify_assert(Matrix3f()) ); + // Test problem size constructors CALL_SUBTEST_5(ComplexEigenSolver(10)); } diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp index 25ef280a1..3ff84c4e0 100644 --- a/test/eigensolver_selfadjoint.cpp +++ b/test/eigensolver_selfadjoint.cpp @@ -105,6 +105,9 @@ template void selfadjointeigensolver(const MatrixType& m) eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps)); VERIFY_IS_APPROX(symmA.template selfadjointView().eigenvalues(), eiSymm.eigenvalues()); + SelfAdjointEigenSolver eiSymmNoEivecs(symmA, false); + VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiSymmNoEivecs.eigenvalues()); + // generalized eigen problem Ax = lBx VERIFY((symmA * eiSymmGen.eigenvectors()).isApprox( symmB * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps)); @@ -115,6 +118,17 @@ template void selfadjointeigensolver(const MatrixType& m) MatrixType id = MatrixType::Identity(rows, cols); VERIFY_IS_APPROX(id.template selfadjointView().operatorNorm(), RealScalar(1)); + + SelfAdjointEigenSolver eiSymmUninitialized; + VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvalues()); + VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors()); + VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt()); + VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt()); + + eiSymmUninitialized.compute(symmA, false); + VERIFY_RAISES_ASSERT(eiSymmUninitialized.eigenvectors()); + VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorSqrt()); + VERIFY_RAISES_ASSERT(eiSymmUninitialized.operatorInverseSqrt()); } void test_eigensolver_selfadjoint() diff --git a/test/hessenberg.cpp b/test/hessenberg.cpp index 16aa564ae..b4cfe288d 100644 --- a/test/hessenberg.cpp +++ b/test/hessenberg.cpp @@ -54,6 +54,13 @@ template void hessenberg(int size = Size) MatrixType cs2Q = cs2.matrixQ(); VERIFY_IS_EQUAL(cs1Q, cs2Q); + // Test assertions for when used uninitialized + HessenbergDecomposition hessUninitialized; + VERIFY_RAISES_ASSERT( hessUninitialized.matrixH() ); + VERIFY_RAISES_ASSERT( hessUninitialized.matrixQ() ); + VERIFY_RAISES_ASSERT( hessUninitialized.householderCoefficients() ); + VERIFY_RAISES_ASSERT( hessUninitialized.packedMatrix() ); + // TODO: Add tests for packedMatrix() and householderCoefficients() }