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

Guard with assert against using decomposition objects uninitialized.

Browse Source
This commit is contained in:
Jitse Niesen 2010-05-30 21:49:35 +01:00
parent 6ce22a61b3
commit db8631b66a
6 changed files with 113 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

29
Eigen/src/Eigenvalues/HessenbergDecomposition.h
View File

@ -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

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

@ -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);
if(computeEigenvectors)
m_eivec.resize(n,n);
if(n==1)
{
m_eivalues.coeffRef(0,0) = ei_real(matrix.coeff(0,0));
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]);
if(computeEigenvectors)
m_eivec.col(i).swap(m_eivec.col(k+i));
}
}
m_isInitialized = true;
m_eigenvectorsOk = computeEigenvectors;
return *this;
}

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

@ -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();

12
test/eigensolver_complex.cpp
View File

@ -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));
}

14
test/eigensolver_selfadjoint.cpp
View File

@ -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()

7
test/hessenberg.cpp
View File

@ -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()
}