diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h index 19002e0eb..63fd2d48b 100644 --- a/Eigen/src/QR/QR.h +++ b/Eigen/src/QR/QR.h @@ -49,11 +49,20 @@ template class QR typedef Matrix MatrixTypeR; typedef Matrix VectorType; + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via QR::compute(const MatrixType&). + */ + QR() : m_qr(), m_hCoeffs(), m_isInitialized(false) {} + QR(const MatrixType& matrix) : m_qr(matrix.rows(), matrix.cols()), - m_hCoeffs(matrix.cols()) + m_hCoeffs(matrix.cols()), + m_isInitialized(false) { - _compute(matrix); + compute(matrix); } /** \deprecated use isInjective() @@ -62,7 +71,11 @@ template class QR * \note Since the rank is computed only once, i.e. the first time it is needed, this * method almost does not perform any further computation. */ - EIGEN_DEPRECATED bool isFullRank() const { return rank() == m_qr.cols(); } + EIGEN_DEPRECATED bool isFullRank() const + { + ei_assert(m_isInitialized && "QR is not initialized."); + return rank() == m_qr.cols(); + } /** \returns the rank of the matrix of which *this is the QR decomposition. * @@ -78,6 +91,7 @@ template class QR */ inline int dimensionOfKernel() const { + ei_assert(m_isInitialized && "QR is not initialized."); return m_qr.cols() - rank(); } @@ -89,6 +103,7 @@ template class QR */ inline bool isInjective() const { + ei_assert(m_isInitialized && "QR is not initialized."); return rank() == m_qr.cols(); } @@ -100,6 +115,7 @@ template class QR */ inline bool isSurjective() const { + ei_assert(m_isInitialized && "QR is not initialized."); return rank() == m_qr.rows(); } @@ -110,6 +126,7 @@ template class QR */ inline bool isInvertible() const { + ei_assert(m_isInitialized && "QR is not initialized."); return isInjective() && isSurjective(); } @@ -117,6 +134,7 @@ template class QR const Part, UpperTriangular> matrixR(void) const { + ei_assert(m_isInitialized && "QR is not initialized."); int cols = m_qr.cols(); return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template part(); } @@ -149,21 +167,21 @@ template class QR MatrixType matrixQ(void) const; - private: - - void _compute(const MatrixType& matrix); + void compute(const MatrixType& matrix); protected: MatrixType m_qr; VectorType m_hCoeffs; mutable int m_rank; mutable bool m_rankIsUptodate; + bool m_isInitialized; }; /** \returns the rank of the matrix of which *this is the QR decomposition. */ template int QR::rank() const { + ei_assert(m_isInitialized && "QR is not initialized."); if (!m_rankIsUptodate) { RealScalar maxCoeff = m_qr.diagonal().cwise().abs().maxCoeff(); @@ -179,10 +197,12 @@ int QR::rank() const #ifndef EIGEN_HIDE_HEAVY_CODE template -void QR::_compute(const MatrixType& matrix) -{ +void QR::compute(const MatrixType& matrix) +{ m_rankIsUptodate = false; m_qr = matrix; + m_hCoeffs.resize(matrix.cols()); + int rows = matrix.rows(); int cols = matrix.cols(); RealScalar eps2 = precision()*precision(); @@ -237,6 +257,7 @@ void QR::_compute(const MatrixType& matrix) m_hCoeffs.coeffRef(k) = 0; } } + m_isInitialized = true; } template @@ -246,6 +267,7 @@ bool QR::solve( ResultType *result ) const { + ei_assert(m_isInitialized && "QR is not initialized."); const int rows = m_qr.rows(); ei_assert(b.rows() == rows); result->resize(rows, b.cols()); @@ -274,6 +296,7 @@ bool QR::solve( template MatrixType QR::matrixQ() const { + ei_assert(m_isInitialized && "QR is not initialized."); // compute the product Q_0 Q_1 ... Q_n-1, // where Q_k is the k-th Householder transformation I - h_k v_k v_k' // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...] diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h index 0a52acf3d..0073a0ccb 100644 --- a/Eigen/src/SVD/SVD.h +++ b/Eigen/src/SVD/SVD.h @@ -61,10 +61,19 @@ template class SVD public: + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via QR::compute(const MatrixType&). + */ + SVD() : m_matU(), m_matV(), m_sigma(), m_isInitialized(false) {} + SVD(const MatrixType& matrix) : m_matU(matrix.rows(), std::min(matrix.rows(), matrix.cols())), m_matV(matrix.cols(),matrix.cols()), - m_sigma(std::min(matrix.rows(),matrix.cols())) + m_sigma(std::min(matrix.rows(),matrix.cols())), + m_isInitialized(false) { compute(matrix); } @@ -72,9 +81,23 @@ template class SVD template bool solve(const MatrixBase &b, ResultType* result) const; - const MatrixUType& matrixU() const { return m_matU; } - const SingularValuesType& singularValues() const { return m_sigma; } - const MatrixVType& matrixV() const { return m_matV; } + const MatrixUType& matrixU() const + { + ei_assert(m_isInitialized && "SVD is not initialized."); + return m_matU; + } + + const SingularValuesType& singularValues() const + { + ei_assert(m_isInitialized && "SVD is not initialized."); + return m_sigma; + } + + const MatrixVType& matrixV() const + { + ei_assert(m_isInitialized && "SVD is not initialized."); + return m_matV; + } void compute(const MatrixType& matrix); SVD& sort(); @@ -95,6 +118,7 @@ template class SVD MatrixVType m_matV; /** \internal */ SingularValuesType m_sigma; + bool m_isInitialized; }; /** Computes / recomputes the SVD decomposition A = U S V^* of \a matrix @@ -473,11 +497,15 @@ void SVD::compute(const MatrixType& matrix) break; } // end big switch } // end iterations + + m_isInitialized = true; } template SVD& SVD::sort() { + ei_assert(m_isInitialized && "SVD is not initialized."); + int mu = m_matU.rows(); int mv = m_matV.rows(); int n = m_matU.cols(); @@ -521,6 +549,8 @@ template template bool SVD::solve(const MatrixBase &b, ResultType* result) const { + ei_assert(m_isInitialized && "SVD is not initialized."); + const int rows = m_matU.rows(); ei_assert(b.rows() == rows); @@ -556,6 +586,7 @@ template void SVD::computeUnitaryPositive(UnitaryType *unitary, PositiveType *positive) const { + ei_assert(m_isInitialized && "SVD is not initialized."); ei_assert(m_matU.cols() == m_matV.cols() && "Polar decomposition is only for square matrices"); if(unitary) *unitary = m_matU * m_matV.adjoint(); if(positive) *positive = m_matV * m_sigma.asDiagonal() * m_matV.adjoint(); @@ -574,6 +605,7 @@ template void SVD::computePositiveUnitary(UnitaryType *positive, PositiveType *unitary) const { + ei_assert(m_isInitialized && "SVD is not initialized."); ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices"); if(unitary) *unitary = m_matU * m_matV.adjoint(); if(positive) *positive = m_matU * m_sigma.asDiagonal() * m_matU.adjoint(); @@ -592,6 +624,7 @@ template template void SVD::computeRotationScaling(RotationType *rotation, ScalingType *scaling) const { + ei_assert(m_isInitialized && "SVD is not initialized."); ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices"); Scalar x = (m_matU * m_matV.adjoint()).determinant(); // so x has absolute value 1 Matrix sv(m_sigma); @@ -618,6 +651,7 @@ template template void SVD::computeScalingRotation(ScalingType *scaling, RotationType *rotation) const { + ei_assert(m_isInitialized && "SVD is not initialized."); ei_assert(m_matU.rows() == m_matV.rows() && "Polar decomposition is only for square matrices"); Scalar x = (m_matU * m_matV.adjoint()).determinant(); // so x has absolute value 1 Matrix sv(m_sigma); diff --git a/test/qr.cpp b/test/qr.cpp index 8fe356444..b71abb782 100644 --- a/test/qr.cpp +++ b/test/qr.cpp @@ -121,6 +121,22 @@ template void qr_invertible() VERIFY(lu.solve(m3, &m2)); } +template void qr_verify_assert() +{ + MatrixType tmp; + + QR qr; + VERIFY_RAISES_ASSERT(qr.isFullRank()) + VERIFY_RAISES_ASSERT(qr.rank()) + VERIFY_RAISES_ASSERT(qr.dimensionOfKernel()) + VERIFY_RAISES_ASSERT(qr.isInjective()) + VERIFY_RAISES_ASSERT(qr.isSurjective()) + VERIFY_RAISES_ASSERT(qr.isInvertible()) + VERIFY_RAISES_ASSERT(qr.matrixR()) + VERIFY_RAISES_ASSERT(qr.solve(tmp,&tmp)) + VERIFY_RAISES_ASSERT(qr.matrixQ()) +} + void test_qr() { for(int i = 0; i < 1; i++) { @@ -144,4 +160,11 @@ void test_qr() // CALL_SUBTEST( qr_invertible() ); // CALL_SUBTEST( qr_invertible() ); } + + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); + CALL_SUBTEST(qr_verify_assert()); } diff --git a/test/svd.cpp b/test/svd.cpp index 688c3f402..4a703369e 100644 --- a/test/svd.cpp +++ b/test/svd.cpp @@ -24,6 +24,7 @@ #include "main.h" #include +#include template void svd(const MatrixType& m) { @@ -85,6 +86,22 @@ template void svd(const MatrixType& m) } } +template void svd_verify_assert() +{ + MatrixType tmp; + + SVD svd; + VERIFY_RAISES_ASSERT(svd.solve(tmp, &tmp)) + VERIFY_RAISES_ASSERT(svd.matrixU()) + VERIFY_RAISES_ASSERT(svd.singularValues()) + VERIFY_RAISES_ASSERT(svd.matrixV()) + VERIFY_RAISES_ASSERT(svd.sort()) + VERIFY_RAISES_ASSERT(svd.computeUnitaryPositive(&tmp,&tmp)) + VERIFY_RAISES_ASSERT(svd.computePositiveUnitary(&tmp,&tmp)) + VERIFY_RAISES_ASSERT(svd.computeRotationScaling(&tmp,&tmp)) + VERIFY_RAISES_ASSERT(svd.computeScalingRotation(&tmp,&tmp)) +} + void test_svd() { for(int i = 0; i < g_repeat; i++) { @@ -96,4 +113,9 @@ void test_svd() // CALL_SUBTEST( svd(MatrixXcd(6,6)) ); // CALL_SUBTEST( svd(MatrixXcf(3,3)) ); } + + CALL_SUBTEST( svd_verify_assert() ); + CALL_SUBTEST( svd_verify_assert() ); + CALL_SUBTEST( svd_verify_assert() ); + CALL_SUBTEST( svd_verify_assert() ); }