@@ -24,7 +24,7 @@ template<typename MatrixType, int UpLo> struct LLT_Traits;
2424 *
2525 * \tparam _MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition
2626 * \tparam _UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper.
27- * The other triangular part won't be read.
27+ * The other triangular part won't be read.
2828 *
2929 * This class performs a LL^T Cholesky decomposition of a symmetric, positive definite
3030 * matrix A such that A = LL^* = U^*U, where L is lower triangular.
@@ -41,14 +41,18 @@ template<typename MatrixType, int UpLo> struct LLT_Traits;
4141 * Example: \include LLT_example.cpp
4242 * Output: \verbinclude LLT_example.out
4343 *
44+ * \b Performance: for best performance, it is recommended to use a column-major storage format
45+ * with the Lower triangular part (the default), or, equivalently, a row-major storage format
46+ * with the Upper triangular part. Otherwise, you might get a 20% slowdown for the full factorization
47+ * step, and rank-updates can be up to 3 times slower.
48+ *
4449 * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
4550 *
51+ * Note that during the decomposition, only the lower (or upper, as defined by _UpLo) triangular part of A is considered.
52+ * Therefore, the strict lower part does not have to store correct values.
53+ *
4654 * \sa MatrixBase::llt(), SelfAdjointView::llt(), class LDLT
4755 */
48- /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH)
49- * Note that during the decomposition, only the upper triangular part of A is considered. Therefore,
50- * the strict lower part does not have to store correct values.
51- */
5256template <typename _MatrixType, int _UpLo> class LLT
5357{
5458 public:
@@ -146,7 +150,7 @@ template<typename _MatrixType, int _UpLo> class LLT
146150 }
147151
148152 template <typename Derived>
149- void solveInPlace (MatrixBase<Derived> &bAndX) const ;
153+ void solveInPlace (const MatrixBase<Derived> &bAndX) const ;
150154
151155 template <typename InputType>
152156 LLT& compute (const EigenBase<InputType>& matrix);
@@ -177,7 +181,7 @@ template<typename _MatrixType, int _UpLo> class LLT
177181 /* * \brief Reports whether previous computation was successful.
178182 *
179183 * \returns \c Success if computation was succesful,
180- * \c NumericalIssue if the matrix.appears to be negative .
184+ * \c NumericalIssue if the matrix.appears not to be positive definite .
181185 */
182186 ComputationInfo info () const
183187 {
@@ -425,7 +429,8 @@ LLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const EigenBase<InputType>
425429 eigen_assert (a.rows ()==a.cols ());
426430 const Index size = a.rows ();
427431 m_matrix.resize (size, size);
428- m_matrix = a.derived ();
432+ if (!internal::is_same_dense (m_matrix, a.derived ()))
433+ m_matrix = a.derived ();
429434
430435 // Compute matrix L1 norm = max abs column sum.
431436 m_l1_norm = RealScalar (0 );
@@ -485,11 +490,14 @@ void LLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const
485490 *
486491 * This version avoids a copy when the right hand side matrix b is not needed anymore.
487492 *
493+ * \warning The parameter is only marked 'const' to make the C++ compiler accept a temporary expression here.
494+ * This function will const_cast it, so constness isn't honored here.
495+ *
488496 * \sa LLT::solve(), MatrixBase::llt()
489497 */
490498template <typename MatrixType, int _UpLo>
491499template <typename Derived>
492- void LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const
500+ void LLT<MatrixType,_UpLo>::solveInPlace(const MatrixBase<Derived> &bAndX) const
493501{
494502 eigen_assert (m_isInitialized && " LLT is not initialized."
495503 eigen_assert (m_matrix.rows ()==bAndX.rows ());
0 commit comments