diff options
Diffstat (limited to 'Eigen/src/Cholesky/LLT.h')
| -rw-r--r-- | Eigen/src/Cholesky/LLT.h | 88 |
1 files changed, 56 insertions, 32 deletions
diff --git a/Eigen/src/Cholesky/LLT.h b/Eigen/src/Cholesky/LLT.h index 87ca8d423..8c9b2b398 100644 --- a/Eigen/src/Cholesky/LLT.h +++ b/Eigen/src/Cholesky/LLT.h @@ -13,6 +13,16 @@ namespace Eigen { namespace internal{ + +template<typename _MatrixType, int _UpLo> struct traits<LLT<_MatrixType, _UpLo> > + : traits<_MatrixType> +{ + typedef MatrixXpr XprKind; + typedef SolverStorage StorageKind; + typedef int StorageIndex; + enum { Flags = 0 }; +}; + template<typename MatrixType, int UpLo> struct LLT_Traits; } @@ -24,7 +34,7 @@ template<typename MatrixType, int UpLo> struct LLT_Traits; * * \tparam _MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition * \tparam _UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper. - * The other triangular part won't be read. + * The other triangular part won't be read. * * This class performs a LL^T Cholesky decomposition of a symmetric, positive definite * matrix A such that A = LL^* = U^*U, where L is lower triangular. @@ -41,27 +51,30 @@ template<typename MatrixType, int UpLo> struct LLT_Traits; * Example: \include LLT_example.cpp * Output: \verbinclude LLT_example.out * + * \b Performance: for best performance, it is recommended to use a column-major storage format + * with the Lower triangular part (the default), or, equivalently, a row-major storage format + * with the Upper triangular part. Otherwise, you might get a 20% slowdown for the full factorization + * step, and rank-updates can be up to 3 times slower. + * * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. * + * Note that during the decomposition, only the lower (or upper, as defined by _UpLo) triangular part of A is considered. + * Therefore, the strict lower part does not have to store correct values. + * * \sa MatrixBase::llt(), SelfAdjointView::llt(), class LDLT */ - /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH) - * 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<typename _MatrixType, int _UpLo> class LLT + : public SolverBase<LLT<_MatrixType, _UpLo> > { public: typedef _MatrixType MatrixType; + typedef SolverBase<LLT> Base; + friend class SolverBase<LLT>; + + EIGEN_GENERIC_PUBLIC_INTERFACE(LLT) enum { - RowsAtCompileTime = MatrixType::RowsAtCompileTime, - ColsAtCompileTime = MatrixType::ColsAtCompileTime, MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime }; - typedef typename MatrixType::Scalar Scalar; - typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; - typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 - typedef typename MatrixType::StorageIndex StorageIndex; enum { PacketSize = internal::packet_traits<Scalar>::size, @@ -96,7 +109,7 @@ template<typename _MatrixType, int _UpLo> class LLT compute(matrix.derived()); } - /** \brief Constructs a LDLT factorization from a given matrix + /** \brief Constructs a LLT factorization from a given matrix * * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when * \c MatrixType is a Eigen::Ref. @@ -125,6 +138,7 @@ template<typename _MatrixType, int _UpLo> class LLT return Traits::getL(m_matrix); } + #ifdef EIGEN_PARSED_BY_DOXYGEN /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. * * Since this LLT class assumes anyway that the matrix A is invertible, the solution @@ -137,16 +151,11 @@ template<typename _MatrixType, int _UpLo> class LLT */ template<typename Rhs> inline const Solve<LLT, Rhs> - solve(const MatrixBase<Rhs>& b) const - { - eigen_assert(m_isInitialized && "LLT is not initialized."); - eigen_assert(m_matrix.rows()==b.rows() - && "LLT::solve(): invalid number of rows of the right hand side matrix b"); - return Solve<LLT, Rhs>(*this, b.derived()); - } + solve(const MatrixBase<Rhs>& b) const; + #endif template<typename Derived> - void solveInPlace(MatrixBase<Derived> &bAndX) const; + void solveInPlace(const MatrixBase<Derived> &bAndX) const; template<typename InputType> LLT& compute(const EigenBase<InputType>& matrix); @@ -176,8 +185,8 @@ template<typename _MatrixType, int _UpLo> class LLT /** \brief Reports whether previous computation was successful. * - * \returns \c Success if computation was succesful, - * \c NumericalIssue if the matrix.appears to be negative. + * \returns \c Success if computation was successful, + * \c NumericalIssue if the matrix.appears not to be positive definite. */ ComputationInfo info() const { @@ -190,18 +199,20 @@ template<typename _MatrixType, int _UpLo> class LLT * This method is provided for compatibility with other matrix decompositions, thus enabling generic code such as: * \code x = decomposition.adjoint().solve(b) \endcode */ - const LLT& adjoint() const { return *this; }; + const LLT& adjoint() const EIGEN_NOEXCEPT { return *this; }; - inline Index rows() const { return m_matrix.rows(); } - inline Index cols() const { return m_matrix.cols(); } + inline EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_matrix.rows(); } + inline EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_matrix.cols(); } template<typename VectorType> - LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1); + LLT & rankUpdate(const VectorType& vec, const RealScalar& sigma = 1); #ifndef EIGEN_PARSED_BY_DOXYGEN template<typename RhsType, typename DstType> - EIGEN_DEVICE_FUNC void _solve_impl(const RhsType &rhs, DstType &dst) const; + + template<bool Conjugate, typename RhsType, typename DstType> + void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const; #endif protected: @@ -425,7 +436,8 @@ LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const EigenBase<InputType> eigen_assert(a.rows()==a.cols()); const Index size = a.rows(); m_matrix.resize(size, size); - m_matrix = a.derived(); + if (!internal::is_same_dense(m_matrix, a.derived())) + m_matrix = a.derived(); // Compute matrix L1 norm = max abs column sum. m_l1_norm = RealScalar(0); @@ -454,7 +466,7 @@ LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const EigenBase<InputType> */ template<typename _MatrixType, int _UpLo> template<typename VectorType> -LLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma) +LLT<_MatrixType,_UpLo> & LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma) { EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType); eigen_assert(v.size()==m_matrix.cols()); @@ -472,8 +484,17 @@ template<typename _MatrixType,int _UpLo> template<typename RhsType, typename DstType> void LLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const { - dst = rhs; - solveInPlace(dst); + _solve_impl_transposed<true>(rhs, dst); +} + +template<typename _MatrixType,int _UpLo> +template<bool Conjugate, typename RhsType, typename DstType> +void LLT<_MatrixType,_UpLo>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const +{ + dst = rhs; + + matrixL().template conjugateIf<!Conjugate>().solveInPlace(dst); + matrixU().template conjugateIf<!Conjugate>().solveInPlace(dst); } #endif @@ -485,11 +506,14 @@ void LLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const * * This version avoids a copy when the right hand side matrix b is not needed anymore. * + * \warning The parameter is only marked 'const' to make the C++ compiler accept a temporary expression here. + * This function will const_cast it, so constness isn't honored here. + * * \sa LLT::solve(), MatrixBase::llt() */ template<typename MatrixType, int _UpLo> template<typename Derived> -void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const +void LLT<MatrixType,_UpLo>::solveInPlace(const MatrixBase<Derived> &bAndX) const { eigen_assert(m_isInitialized && "LLT is not initialized."); eigen_assert(m_matrix.rows()==bAndX.rows()); |