diff --git a/Eigen/src/Cholesky/LDLT.h b/Eigen/src/Cholesky/LDLT.h index 81f3aaa32..60cb98307 100644 --- a/Eigen/src/Cholesky/LDLT.h +++ b/Eigen/src/Cholesky/LDLT.h @@ -27,6 +27,8 @@ #ifndef EIGEN_LDLT_H #define EIGEN_LDLT_H +template struct LDLT_Traits; + /** \ingroup cholesky_Module * * \class LDLT @@ -52,7 +54,7 @@ * Note that during the decomposition, only the upper triangular part of A is considered. Therefore, * the strict lower part does not have to store correct values. */ -template class LDLT +template class LDLT { public: typedef _MatrixType MatrixType; @@ -61,7 +63,8 @@ template class LDLT ColsAtCompileTime = MatrixType::ColsAtCompileTime, Options = MatrixType::Options, MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, - MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime + MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, + UpLo = _UpLo }; typedef typename MatrixType::Scalar Scalar; typedef typename NumTraits::Real RealScalar; @@ -69,6 +72,8 @@ template class LDLT typedef typename ei_plain_col_type::type IntColVectorType; typedef Matrix TmpMatrixType; + typedef LDLT_Traits Traits; + /** \brief Default Constructor. * * The default constructor is useful in cases in which the user intends to @@ -100,11 +105,18 @@ template class LDLT compute(matrix); } - /** \returns an expression of the lower triangular matrix L */ - inline TriangularView matrixL(void) const + /** \returns a view of the upper triangular matrix U */ + inline typename Traits::MatrixU matrixU() const { ei_assert(m_isInitialized && "LDLT is not initialized."); - return m_matrix; + return Traits::getU(m_matrix); + } + + /** \returns a view of the lower triangular matrix L */ + inline typename Traits::MatrixL matrixL() const + { + ei_assert(m_isInitialized && "LDLT is not initialized."); + return Traits::getL(m_matrix); } /** \returns a vector of integers, whose size is the number of rows of the matrix being decomposed, @@ -186,14 +198,115 @@ template class LDLT IntColVectorType m_p; IntColVectorType m_transpositions; // FIXME do we really need to store permanently the transpositions? TmpMatrixType m_temporary; - Index m_sign; + int m_sign; bool m_isInitialized; }; +template struct ei_ldlt_inplace; + +template<> struct ei_ldlt_inplace +{ + template + static bool unblocked(MatrixType& mat, Transpositions& transpositions, Workspace& temp, int* sign=0) + { + typedef typename MatrixType::Scalar Scalar; + typedef typename MatrixType::RealScalar RealScalar; + typedef typename MatrixType::Index Index; + ei_assert(mat.rows()==mat.cols()); + const Index size = mat.rows(); + + if (size <= 1) + { + transpositions.setZero(); + if(sign) + *sign = ei_real(mat.coeff(0,0))>0 ? 1:-1; + return true; + } + + RealScalar cutoff = 0, biggest_in_corner; + + for (Index k = 0; k < size; ++k) + { + // Find largest diagonal element + Index index_of_biggest_in_corner; + biggest_in_corner = mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner); + index_of_biggest_in_corner += k; + + if(k == 0) + { + // The biggest overall is the point of reference to which further diagonals + // are compared; if any diagonal is negligible compared + // to the largest overall, the algorithm bails. + cutoff = ei_abs(NumTraits::epsilon() * biggest_in_corner); + + if(sign) + *sign = ei_real(mat.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1; + } + + // Finish early if the matrix is not full rank. + if(biggest_in_corner < cutoff) + { + for(Index i = k; i < size; i++) transpositions.coeffRef(i) = i; + break; + } + + transpositions.coeffRef(k) = index_of_biggest_in_corner; + if(k != index_of_biggest_in_corner) + { + mat.row(k).swap(mat.row(index_of_biggest_in_corner)); + mat.col(k).swap(mat.col(index_of_biggest_in_corner)); + } + + Index rs = size - k - 1; + Block A21(mat,k+1,k,rs,1); + Block A10(mat,k,0,1,k); + Block A20(mat,k+1,0,rs,k); + + if(k>0) + { + temp.head(k) = mat.diagonal().head(k).asDiagonal() * A10.adjoint(); + mat.coeffRef(k,k) -= (A10 * temp.head(k)).value(); + if(rs>0) + A21.noalias() -= A20 * temp.head(k); + } + if((rs>0) && (ei_abs(mat.coeffRef(k,k)) > cutoff)) + A21 /= mat.coeffRef(k,k); + } + + return true; + } +}; + +template<> struct ei_ldlt_inplace +{ + template + static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, Transpositions& transpositions, Workspace& temp, int* sign=0) + { + Transpose matt(mat); + return ei_ldlt_inplace::unblocked(matt, transpositions, temp, sign); + } +}; + +template struct LDLT_Traits +{ + typedef TriangularView MatrixL; + typedef TriangularView MatrixU; + inline static MatrixL getL(const MatrixType& m) { return m; } + inline static MatrixU getU(const MatrixType& m) { return m.adjoint(); } +}; + +template struct LDLT_Traits +{ + typedef TriangularView MatrixL; + typedef TriangularView MatrixU; + inline static MatrixL getL(const MatrixType& m) { return m.adjoint(); } + inline static MatrixU getU(const MatrixType& m) { return m; } +}; + /** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix */ -template -LDLT& LDLT::compute(const MatrixType& a) +template +LDLT& LDLT::compute(const MatrixType& a) { ei_assert(a.rows()==a.cols()); const Index size = a.rows(); @@ -203,85 +316,29 @@ LDLT& LDLT::compute(const MatrixType& a) m_p.resize(size); m_transpositions.resize(size); m_isInitialized = false; - - if (size <= 1) - { - m_p.setZero(); - m_transpositions.setZero(); - m_sign = ei_real(a.coeff(0,0))>0 ? 1:-1; - m_isInitialized = true; - return *this; - } - - RealScalar cutoff = 0, biggest_in_corner; - - // By using a temporary, packet-aligned products are guarenteed. In the LLT - // case this is unnecessary because the diagonal is included and will always - // have optimal alignment. m_temporary.resize(size); - for (Index j = 0; j < size; ++j) - { - // Find largest diagonal element - Index index_of_biggest_in_corner; - biggest_in_corner = m_matrix.diagonal().tail(size-j).cwiseAbs().maxCoeff(&index_of_biggest_in_corner); - index_of_biggest_in_corner += j; + ei_ldlt_inplace::unblocked(m_matrix, m_transpositions, m_temporary, &m_sign); - if(j == 0) - { - // The biggest overall is the point of reference to which further diagonals - // are compared; if any diagonal is negligible compared - // to the largest overall, the algorithm bails. - cutoff = ei_abs(NumTraits::epsilon() * biggest_in_corner); - - m_sign = ei_real(m_matrix.diagonal().coeff(index_of_biggest_in_corner)) > 0 ? 1 : -1; - } - - // Finish early if the matrix is not full rank. - if(biggest_in_corner < cutoff) - { - for(Index i = j; i < size; i++) m_transpositions.coeffRef(i) = i; - break; - } - - m_transpositions.coeffRef(j) = index_of_biggest_in_corner; - if(j != index_of_biggest_in_corner) - { - m_matrix.row(j).swap(m_matrix.row(index_of_biggest_in_corner)); - m_matrix.col(j).swap(m_matrix.col(index_of_biggest_in_corner)); - } - - Index rs = size - j - 1; - Block A21(m_matrix,j+1,j,rs,1); - Block A10(m_matrix,j,0,1,j); - Block A20(m_matrix,j+1,0,rs,j); - - if(j>0) - { - m_temporary.head(j) = m_matrix.diagonal().head(j).asDiagonal() * A10.adjoint(); - m_matrix.coeffRef(j,j) -= (A10 * m_temporary.head(j)).value(); - if(rs>0) - A21.noalias() -= A20 * m_temporary.head(j); - } - if((rs>0) && (ei_abs(m_matrix.coeffRef(j,j)) > cutoff)) - A21 /= m_matrix.coeffRef(j,j); - } + if(size==0) + m_p.setZero(); // Reverse applied swaps to get P matrix. for(Index k = 0; k < size; ++k) m_p.coeffRef(k) = k; for(Index k = size-1; k >= 0; --k) { std::swap(m_p.coeffRef(k), m_p.coeffRef(m_transpositions.coeff(k))); } - + m_isInitialized = true; return *this; } -template -struct ei_solve_retval, Rhs> - : ei_solve_retval_base, Rhs> +template +struct ei_solve_retval, Rhs> + : ei_solve_retval_base, Rhs> { - EIGEN_MAKE_SOLVE_HELPERS(LDLT<_MatrixType>,Rhs) + typedef LDLT<_MatrixType,_UpLo> LDLTType; + EIGEN_MAKE_SOLVE_HELPERS(LDLTType,Rhs) template void evalTo(Dest& dst) const { @@ -301,9 +358,9 @@ struct ei_solve_retval, Rhs> * * \sa LDLT::solve(), MatrixBase::ldlt() */ -template +template template -bool LDLT::solveInPlace(MatrixBase &bAndX) const +bool LDLT::solveInPlace(MatrixBase &bAndX) const { ei_assert(m_isInitialized && "LDLT is not initialized."); const Index size = m_matrix.rows(); @@ -314,13 +371,13 @@ bool LDLT::solveInPlace(MatrixBase &bAndX) const // y = L^-1 z //matrixL().solveInPlace(bAndX); - m_matrix.template triangularView().solveInPlace(bAndX); + matrixL().solveInPlace(bAndX); // w = D^-1 y bAndX = m_matrix.diagonal().asDiagonal().inverse() * bAndX; // u = L^-T w - m_matrix.adjoint().template triangularView().solveInPlace(bAndX); + matrixU().solveInPlace(bAndX); // x = P^T u for (Index i = size-1; i >= 0; --i) bAndX.row(m_transpositions.coeff(i)).swap(bAndX.row(i)); @@ -331,8 +388,8 @@ bool LDLT::solveInPlace(MatrixBase &bAndX) const /** \returns the matrix represented by the decomposition, * i.e., it returns the product: P^T L D L^* P. * This function is provided for debug purpose. */ -template -MatrixType LDLT::reconstructedMatrix() const +template +MatrixType LDLT::reconstructedMatrix() const { ei_assert(m_isInitialized && "LDLT is not initialized."); const Index size = m_matrix.rows(); @@ -342,7 +399,7 @@ MatrixType LDLT::reconstructedMatrix() const // PI for(Index i = 0; i < size; ++i) res.row(m_transpositions.coeff(i)).swap(res.row(i)); // L^* P - res = matrixL().adjoint() * res; + res = matrixU() * res; // D(L^*P) res = vectorD().asDiagonal() * res; // L(DL^*P) @@ -353,6 +410,16 @@ MatrixType LDLT::reconstructedMatrix() const return res; } +/** \cholesky_module + * \returns the Cholesky decomposition with full pivoting without square root of \c *this + */ +template +inline const LDLT::PlainObject, UpLo> +SelfAdjointView::ldlt() const +{ + return LDLT(m_matrix); +} + /** \cholesky_module * \returns the Cholesky decomposition with full pivoting without square root of \c *this */ diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 29fa465e1..6e853436c 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -162,7 +162,6 @@ template class LLT bool m_isInitialized; }; -// forward declaration (defined at the end of this file) template struct ei_llt_inplace; template<> struct ei_llt_inplace diff --git a/Eigen/src/Core/SelfAdjointView.h b/Eigen/src/Core/SelfAdjointView.h index eed3f9336..84c4dc521 100644 --- a/Eigen/src/Core/SelfAdjointView.h +++ b/Eigen/src/Core/SelfAdjointView.h @@ -153,7 +153,7 @@ template class SelfAdjointView /////////// Cholesky module /////////// const LLT llt() const; - const LDLT ldlt() const; + const LDLT ldlt() const; /////////// Eigenvalue module /////////// diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h index b3bc9c161..5cf62f4c6 100644 --- a/Eigen/src/Core/util/ForwardDeclarations.h +++ b/Eigen/src/Core/util/ForwardDeclarations.h @@ -158,7 +158,7 @@ template class FullPivHouseholderQR; template class SVD; template class JacobiSVD; template class LLT; -template class LDLT; +template class LDLT; template class HouseholderSequence; template class PlanarRotation; diff --git a/test/cholesky.cpp b/test/cholesky.cpp index 0ae26c7d5..d403af7ba 100644 --- a/test/cholesky.cpp +++ b/test/cholesky.cpp @@ -118,11 +118,18 @@ template void cholesky(const MatrixType& m) } { - LDLT ldlt(symm); - VERIFY_IS_APPROX(symm, ldlt.reconstructedMatrix()); - vecX = ldlt.solve(vecB); + LDLT ldltlo(symm); + VERIFY_IS_APPROX(symm, ldltlo.reconstructedMatrix()); + vecX = ldltlo.solve(vecB); VERIFY_IS_APPROX(symm * vecX, vecB); - matX = ldlt.solve(matB); + matX = ldltlo.solve(matB); + VERIFY_IS_APPROX(symm * matX, matB); + + LDLT ldltup(symm); + VERIFY_IS_APPROX(symm, ldltup.reconstructedMatrix()); + vecX = ldltup.solve(vecB); + VERIFY_IS_APPROX(symm * vecX, vecB); + matX = ldltup.solve(matB); VERIFY_IS_APPROX(symm * matX, matB); }