diff options
author | Jitse Niesen <jitse@maths.leeds.ac.uk> | 2010-05-30 21:49:35 +0100 |
---|---|---|
committer | Jitse Niesen <jitse@maths.leeds.ac.uk> | 2010-05-30 21:49:35 +0100 |
commit | db8631b66a4f5cb9957a58cc629be5ef438e3059 (patch) | |
tree | 50a3ef1425a13e183f42ae36dd34087899e06460 | |
parent | 6ce22a61b3a2560a5464b368bee64555645da59f (diff) |
Guard with assert against using decomposition objects uninitialized.
-rw-r--r-- | Eigen/src/Eigenvalues/HessenbergDecomposition.h | 29 | ||||
-rw-r--r-- | Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h | 51 | ||||
-rw-r--r-- | Eigen/src/Eigenvalues/Tridiagonalization.h | 27 | ||||
-rw-r--r-- | test/eigensolver_complex.cpp | 12 | ||||
-rw-r--r-- | test/eigensolver_selfadjoint.cpp | 14 | ||||
-rw-r--r-- | 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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> class HessenbergDecomposition */ HessenbergDecompositionMatrixHReturnType<MatrixType> matrixH() const { + ei_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); return HessenbergDecompositionMatrixHReturnType<MatrixType>(*this); } @@ -257,6 +277,7 @@ template<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> 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<typename MatrixType> SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>::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<MatrixType>& SelfAdjointEigenSolver<MatrixType>::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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> 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<typename _MatrixType> class Tridiagonalization MatrixType m_matrix; CoeffVectorType m_hCoeffs; + bool m_isInitialized; }; template<typename MatrixType> const typename Tridiagonalization<MatrixType>::DiagonalReturnType Tridiagonalization<MatrixType>::diagonal() const { + ei_assert(m_isInitialized && "Tridiagonalization is not initialized."); return m_matrix.diagonal(); } @@ -325,6 +342,7 @@ template<typename MatrixType> const typename Tridiagonalization<MatrixType>::SubDiagonalReturnType Tridiagonalization<MatrixType>::subDiagonal() const { + ei_assert(m_isInitialized && "Tridiagonalization is not initialized."); Index n = m_matrix.rows(); return Block<MatrixType,SizeMinusOne,SizeMinusOne>(m_matrix, 1, 0, n-1,n-1).diagonal(); } @@ -335,6 +353,7 @@ Tridiagonalization<MatrixType>::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<Scalar>().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<typename MatrixType> void eigensolver(const MatrixType& m) VERIFY_IS_APPROX(id.operatorNorm(), RealScalar(1)); } +template<typename MatrixType> void eigensolver_verify_assert() +{ + ComplexEigenSolver<MatrixType> 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<std::complex<float>, 1, 1>()) ); + CALL_SUBTEST_4( eigensolver_verify_assert(Matrix3f()) ); + // Test problem size constructors CALL_SUBTEST_5(ComplexEigenSolver<MatrixXf>(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<typename MatrixType> void selfadjointeigensolver(const MatrixType& m) eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps)); VERIFY_IS_APPROX(symmA.template selfadjointView<Lower>().eigenvalues(), eiSymm.eigenvalues()); + SelfAdjointEigenSolver<MatrixType> 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<typename MatrixType> void selfadjointeigensolver(const MatrixType& m) MatrixType id = MatrixType::Identity(rows, cols); VERIFY_IS_APPROX(id.template selfadjointView<Lower>().operatorNorm(), RealScalar(1)); + + SelfAdjointEigenSolver<MatrixType> 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<typename Scalar,int Size> void hessenberg(int size = Size) MatrixType cs2Q = cs2.matrixQ(); VERIFY_IS_EQUAL(cs1Q, cs2Q); + // Test assertions for when used uninitialized + HessenbergDecomposition<MatrixType> 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() } |