diff options
| -rw-r--r-- | Eigen/src/LU/LU.h | 68 | ||||
| -rw-r--r-- | Eigen/src/LU/PartialLU.h | 46 | ||||
| -rw-r--r-- | test/lu.cpp | 38 |
3 files changed, 138 insertions, 14 deletions
diff --git a/Eigen/src/LU/LU.h b/Eigen/src/LU/LU.h index e6333b990..560f8d06a 100644 --- a/Eigen/src/LU/LU.h +++ b/Eigen/src/LU/LU.h @@ -92,12 +92,22 @@ template<typename MatrixType> class LU MatrixType::MaxColsAtCompileTime // so it has the same number of rows and at most as many columns. > ImageResultType; + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via LU::compute(const MatrixType&). + */ + LU(); + /** Constructor. * * \param matrix the matrix of which to compute the LU decomposition. */ LU(const MatrixType& matrix); + void compute(const MatrixType& matrix); + /** \returns the LU decomposition matrix: the upper-triangular part is U, the * unit-lower-triangular part is L (at least for square matrices; in the non-square * case, special care is needed, see the documentation of class LU). @@ -106,6 +116,7 @@ template<typename MatrixType> class LU */ inline const MatrixType& matrixLU() const { + ei_assert(m_originalMatrix != 0 && "LU is not initialized."); return m_lu; } @@ -117,6 +128,7 @@ template<typename MatrixType> class LU */ inline const IntColVectorType& permutationP() const { + ei_assert(m_originalMatrix != 0 && "LU is not initialized."); return m_p; } @@ -128,6 +140,7 @@ template<typename MatrixType> class LU */ inline const IntRowVectorType& permutationQ() const { + ei_assert(m_originalMatrix != 0 && "LU is not initialized."); return m_q; } @@ -243,6 +256,7 @@ template<typename MatrixType> class LU */ inline int rank() const { + ei_assert(m_originalMatrix != 0 && "LU is not initialized."); return m_rank; } @@ -253,6 +267,7 @@ template<typename MatrixType> class LU */ inline int dimensionOfKernel() const { + ei_assert(m_originalMatrix != 0 && "LU is not initialized."); return m_lu.cols() - m_rank; } @@ -264,6 +279,7 @@ template<typename MatrixType> class LU */ inline bool isInjective() const { + ei_assert(m_originalMatrix != 0 && "LU is not initialized."); return m_rank == m_lu.cols(); } @@ -275,6 +291,7 @@ template<typename MatrixType> class LU */ inline bool isSurjective() const { + ei_assert(m_originalMatrix != 0 && "LU is not initialized."); return m_rank == m_lu.rows(); } @@ -285,6 +302,7 @@ template<typename MatrixType> class LU */ inline bool isInvertible() const { + ei_assert(m_originalMatrix != 0 && "LU is not initialized."); return isInjective() && isSurjective(); } @@ -317,7 +335,7 @@ template<typename MatrixType> class LU } protected: - const MatrixType& m_originalMatrix; + const MatrixType* m_originalMatrix; MatrixType m_lu; IntColVectorType m_p; IntRowVectorType m_q; @@ -327,12 +345,38 @@ template<typename MatrixType> class LU }; template<typename MatrixType> +LU<MatrixType>::LU() + : m_originalMatrix(0), + m_lu(), + m_p(), + m_q(), + m_det_pq(0), + m_rank(-1), + m_precision(precision<RealScalar>()) +{ +} + +template<typename MatrixType> LU<MatrixType>::LU(const MatrixType& matrix) - : m_originalMatrix(matrix), - m_lu(matrix), - m_p(matrix.rows()), - m_q(matrix.cols()) + : m_originalMatrix(0), + m_lu(), + m_p(), + m_q(), + m_det_pq(0), + m_rank(-1), + m_precision(precision<RealScalar>()) { + compute(matrix); +} + +template<typename MatrixType> +void LU<MatrixType>::compute(const MatrixType& matrix) +{ + m_originalMatrix = &matrix; + m_lu = matrix; + m_p.resize(matrix.rows()); + m_q.resize(matrix.cols()); + const int size = matrix.diagonalSize(); const int rows = matrix.rows(); const int cols = matrix.cols(); @@ -409,6 +453,7 @@ LU<MatrixType>::LU(const MatrixType& matrix) template<typename MatrixType> typename ei_traits<MatrixType>::Scalar LU<MatrixType>::determinant() const { + ei_assert(m_originalMatrix != 0 && "LU is not initialized."); return Scalar(m_det_pq) * m_lu.diagonal().prod(); } @@ -462,17 +507,20 @@ template<typename MatrixType> template<typename ImageMatrixType> void LU<MatrixType>::computeImage(ImageMatrixType *result) const { - ei_assert(m_rank > 0); - result->resize(m_originalMatrix.rows(), m_rank); + // can be caused by a rank deficient matrix or by calling computeImage on an + // unitialized LU object + ei_assert(m_rank > 0 && "LU is not initialized or Matrix has rank zero."); + result->resize(m_originalMatrix->rows(), m_rank); for(int i = 0; i < m_rank; ++i) - result->col(i) = m_originalMatrix.col(m_q.coeff(i)); + result->col(i) = m_originalMatrix->col(m_q.coeff(i)); } template<typename MatrixType> const typename LU<MatrixType>::ImageResultType LU<MatrixType>::image() const { - ImageResultType result(m_originalMatrix.rows(), m_rank); + ei_assert(m_originalMatrix != 0 && "LU is not initialized."); + ImageResultType result(m_originalMatrix->rows(), m_rank); computeImage(&result); return result; } @@ -484,6 +532,8 @@ bool LU<MatrixType>::solve( ResultType *result ) const { + ei_assert(m_originalMatrix != 0 && "LU is not initialized."); + /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}. * So we proceed as follows: * Step 1: compute c = Pb. diff --git a/Eigen/src/LU/PartialLU.h b/Eigen/src/LU/PartialLU.h index 7ca4f4f51..c78e40b8f 100644 --- a/Eigen/src/LU/PartialLU.h +++ b/Eigen/src/LU/PartialLU.h @@ -70,6 +70,14 @@ template<typename MatrixType> class PartialLU MatrixType::MaxRowsAtCompileTime) }; + /** + * \brief Default Constructor. + * + * The default constructor is useful in cases in which the user intends to + * perform decompositions via PartialLU::compute(const MatrixType&). + */ + PartialLU(); + /** Constructor. * * \param matrix the matrix of which to compute the LU decomposition. @@ -79,6 +87,8 @@ template<typename MatrixType> class PartialLU */ PartialLU(const MatrixType& matrix); + void compute(const MatrixType& matrix); + /** \returns the LU decomposition matrix: the upper-triangular part is U, the * unit-lower-triangular part is L (at least for square matrices; in the non-square * case, special care is needed, see the documentation of class LU). @@ -87,8 +97,9 @@ template<typename MatrixType> class PartialLU */ inline const MatrixType& matrixLU() const { + ei_assert(m_isInitialized && "PartialLU is not initialized."); return m_lu; - } + } /** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed, * representing the P permutation i.e. the permutation of the rows. For its precise meaning, @@ -96,6 +107,7 @@ template<typename MatrixType> class PartialLU */ inline const IntColVectorType& permutationP() const { + ei_assert(m_isInitialized && "PartialLU is not initialized."); return m_p; } @@ -164,18 +176,37 @@ template<typename MatrixType> class PartialLU } protected: - const MatrixType& m_originalMatrix; MatrixType m_lu; IntColVectorType m_p; int m_det_p; + bool m_isInitialized; }; template<typename MatrixType> +PartialLU<MatrixType>::PartialLU() + : m_lu(), + m_p(), + m_det_p(0), + m_isInitialized(false) +{ +} + +template<typename MatrixType> PartialLU<MatrixType>::PartialLU(const MatrixType& matrix) - : m_originalMatrix(matrix), - m_lu(matrix), - m_p(matrix.rows()) + : m_lu(), + m_p(), + m_det_p(0), + m_isInitialized(false) { + compute(matrix); +} + +template<typename MatrixType> +void PartialLU<MatrixType>::compute(const MatrixType& matrix) +{ + m_lu = matrix; + m_p.resize(matrix.rows()); + ei_assert(matrix.rows() == matrix.cols() && "PartialLU is only for square (and moreover invertible) matrices"); const int size = matrix.rows(); @@ -213,11 +244,14 @@ PartialLU<MatrixType>::PartialLU(const MatrixType& matrix) std::swap(m_p.coeffRef(k), m_p.coeffRef(rows_transpositions.coeff(k))); m_det_p = (number_of_transpositions%2) ? -1 : 1; + + m_isInitialized = true; } template<typename MatrixType> typename ei_traits<MatrixType>::Scalar PartialLU<MatrixType>::determinant() const { + ei_assert(m_isInitialized && "PartialLU is not initialized."); return Scalar(m_det_p) * m_lu.diagonal().prod(); } @@ -228,6 +262,8 @@ void PartialLU<MatrixType>::solve( ResultType *result ) const { + ei_assert(m_isInitialized && "PartialLU is not initialized."); + /* The decomposition PA = LU can be rewritten as A = P^{-1} L U. * So we proceed as follows: * Step 1: compute c = Pb. diff --git a/test/lu.cpp b/test/lu.cpp index 625241330..923acee3e 100644 --- a/test/lu.cpp +++ b/test/lu.cpp @@ -92,6 +92,37 @@ template<typename MatrixType> void lu_invertible() VERIFY(lu.solve(m3, &m2)); } +template<typename MatrixType> void lu_verify_assert() +{ + MatrixType tmp; + + LU<MatrixType> lu; + VERIFY_RAISES_ASSERT(lu.matrixLU()) + VERIFY_RAISES_ASSERT(lu.permutationP()) + VERIFY_RAISES_ASSERT(lu.permutationQ()) + VERIFY_RAISES_ASSERT(lu.computeKernel(&tmp)) + VERIFY_RAISES_ASSERT(lu.computeImage(&tmp)) + VERIFY_RAISES_ASSERT(lu.kernel()) + VERIFY_RAISES_ASSERT(lu.image()) + VERIFY_RAISES_ASSERT(lu.solve(tmp,&tmp)) + VERIFY_RAISES_ASSERT(lu.determinant()) + VERIFY_RAISES_ASSERT(lu.rank()) + VERIFY_RAISES_ASSERT(lu.dimensionOfKernel()) + VERIFY_RAISES_ASSERT(lu.isInjective()) + VERIFY_RAISES_ASSERT(lu.isSurjective()) + VERIFY_RAISES_ASSERT(lu.isInvertible()) + VERIFY_RAISES_ASSERT(lu.computeInverse(&tmp)) + VERIFY_RAISES_ASSERT(lu.inverse()) + + PartialLU<MatrixType> plu; + VERIFY_RAISES_ASSERT(plu.matrixLU()) + VERIFY_RAISES_ASSERT(plu.permutationP()) + VERIFY_RAISES_ASSERT(plu.solve(tmp,&tmp)) + VERIFY_RAISES_ASSERT(plu.determinant()) + VERIFY_RAISES_ASSERT(plu.computeInverse(&tmp)) + VERIFY_RAISES_ASSERT(plu.inverse()) +} + void test_lu() { for(int i = 0; i < g_repeat; i++) { @@ -104,4 +135,11 @@ void test_lu() CALL_SUBTEST( lu_invertible<MatrixXcf>() ); CALL_SUBTEST( lu_invertible<MatrixXcd>() ); } + + CALL_SUBTEST( lu_verify_assert<Matrix3f>() ); + CALL_SUBTEST( lu_verify_assert<Matrix3d>() ); + CALL_SUBTEST( lu_verify_assert<MatrixXf>() ); + CALL_SUBTEST( lu_verify_assert<MatrixXd>() ); + CALL_SUBTEST( lu_verify_assert<MatrixXcf>() ); + CALL_SUBTEST( lu_verify_assert<MatrixXcd>() ); } |