1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_LLT_H 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_LLT_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal{ 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int UpLo> struct LLT_Traits; 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup Cholesky_Module 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class LLT 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Standard Cholesky decomposition (LL^T) of a matrix and associated features 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param MatrixType the type of the matrix of which we are computing the LL^T Cholesky decomposition 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper. 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The other triangular part won't be read. 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class performs a LL^T Cholesky decomposition of a symmetric, positive definite 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix A such that A = LL^* = U^*U, where L is lower triangular. 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * While the Cholesky decomposition is particularly useful to solve selfadjoint problems like D^*D x = b, 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * for that purpose, we recommend the Cholesky decomposition without square root which is more stable 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and even faster. Nevertheless, this standard Cholesky decomposition remains useful in many other 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * situations like generalised eigen problems with hermitian matrices. 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Remember that Cholesky decompositions are not rank-revealing. This LLT decomposition is only stable on positive definite matrices, 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * use LDLT instead for the semidefinite case. Also, do not use a Cholesky decomposition to determine whether a system of equations 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * has a solution. 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include LLT_example.cpp 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude LLT_example.out 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::llt(), class LDLT 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /* HEY THIS DOX IS DISABLED BECAUSE THERE's A BUG EITHER HERE OR IN LDLT ABOUT THAT (OR BOTH) 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Note that during the decomposition, only the upper triangular part of A is considered. Therefore, 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the strict lower part does not have to store correct values. 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int _UpLo> class LLT 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowsAtCompileTime = MatrixType::RowsAtCompileTime, 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColsAtCompileTime = MatrixType::ColsAtCompileTime, 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Options = MatrixType::Options, 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath PacketSize = internal::packet_traits<Scalar>::size, 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath AlignmentMask = int(PacketSize)-1, 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath UpLo = _UpLo 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef internal::LLT_Traits<MatrixType,UpLo> Traits; 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Default Constructor. 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The default constructor is useful in cases in which the user intends to 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * perform decompositions via LLT::compute(const MatrixType&). 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LLT() : m_matrix(), m_isInitialized(false) {} 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Default Constructor with memory preallocation 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Like the default constructor but with preallocation of the internal data 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * according to the specified problem \a size. 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa LLT() 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LLT(Index size) : m_matrix(size, size), 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false) {} 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LLT(const MatrixType& matrix) 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_matrix(matrix.rows(), matrix.cols()), 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false) 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix); 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns a view of the upper triangular matrix U */ 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline typename Traits::MatrixU matrixU() const 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LLT is not initialized."); 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Traits::getU(m_matrix); 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns a view of the lower triangular matrix L */ 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline typename Traits::MatrixL matrixL() const 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LLT is not initialized."); 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Traits::getL(m_matrix); 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Since this LLT class assumes anyway that the matrix A is invertible, the solution 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * theoretically exists and is unique regardless of b. 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include LLT_solve.cpp 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude LLT_solve.out 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa solveInPlace(), MatrixBase::llt() 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs> 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const internal::solve_retval<LLT, Rhs> 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solve(const MatrixBase<Rhs>& b) const 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LLT is not initialized."); 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_matrix.rows()==b.rows() 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath && "LLT::solve(): invalid number of rows of the right hand side matrix b"); 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return internal::solve_retval<LLT, Rhs>(*this, b.derived()); 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath #ifdef EIGEN2_SUPPORT 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename OtherDerived, typename ResultType> 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *result = this->solve(b); 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool isPositiveDefinite() const { return true; } 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath #endif 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Derived> 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void solveInPlace(MatrixBase<Derived> &bAndX) const; 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LLT& compute(const MatrixType& matrix); 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the LLT decomposition matrix 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * TODO: document the storage layout 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const MatrixType& matrixLLT() const 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LLT is not initialized."); 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_matrix; 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType reconstructedMatrix() const; 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Reports whether previous computation was successful. 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns \c Success if computation was succesful, 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \c NumericalIssue if the matrix.appears to be negative. 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo info() const 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LLT is not initialized."); 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_info; 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index rows() const { return m_matrix.rows(); } 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index cols() const { return m_matrix.cols(); } 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename VectorType> 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LLT rankUpdate(const VectorType& vec, const RealScalar& sigma = 1); 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Used to compute and store L 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The strict upper part is not used and even not initialized. 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType m_matrix; 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo m_info; 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar, int UpLo> struct llt_inplace; 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, typename VectorType> 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstatic typename MatrixType::Index llt_rank_update_lower(MatrixType& mat, const VectorType& vec, const typename MatrixType::RealScalar& sigma) 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::ColXpr ColXpr; 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::remove_all<ColXpr>::type ColXprCleaned; 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename ColXprCleaned::SegmentReturnType ColXprSegment; 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,1> TempVectorType; 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename TempVectorType::SegmentReturnType TempVecSegment; 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int n = mat.cols(); 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(mat.rows()==n && vec.size()==n); 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath TempVectorType temp; 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(sigma>0) 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // This version is based on Givens rotations. 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // It is faster than the other one below, but only works for updates, 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // i.e., for sigma > 0 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp = sqrt(sigma) * vec; 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=0; i<n; ++i) 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiRotation<Scalar> g; 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath g.makeGivens(mat(i,i), -temp(i), &mat(i,i)); 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath int rs = n-i-1; 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(rs>0) 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColXprSegment x(mat.col(i).tail(rs)); 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath TempVecSegment y(temp.tail(rs)); 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath apply_rotation_in_the_plane(x, y, g); 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp = vec; 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar beta = 1; 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int j=0; j<n; ++j) 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar Ljj = real(mat.coeff(j,j)); 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar dj = abs2(Ljj); 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar wj = temp.coeff(j); 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar swj2 = sigma*abs2(wj); 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar gamma = dj*beta + swj2; 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar x = dj + swj2/beta; 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (x<=RealScalar(0)) 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return j; 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar nLjj = sqrt(x); 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat.coeffRef(j,j) = nLjj; 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath beta += swj2/dj; 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Update the terms of L 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rs = n-j-1; 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(rs) 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs); 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(gamma != 0) 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*conj(wj)/gamma)*temp.tail(rs); 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return -1; 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> struct llt_inplace<Scalar, Lower> 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static typename MatrixType::Index unblocked(MatrixType& mat) 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(mat.rows()==mat.cols()); 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = mat.rows(); 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index k = 0; k < size; ++k) 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rs = size-k-1; // remaining size 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1); 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k); 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k); 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar x = real(mat.coeff(k,k)); 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k>0) x -= A10.squaredNorm(); 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (x<=RealScalar(0)) 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return k; 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat.coeffRef(k,k) = x = sqrt(x); 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint(); 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (rs>0) A21 *= RealScalar(1)/x; 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return -1; 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static typename MatrixType::Index blocked(MatrixType& m) 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m.rows()==m.cols()); 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index size = m.rows(); 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(size<32) 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return unblocked(m); 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index blockSize = size/8; 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath blockSize = (blockSize/16)*16; 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128)); 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index k=0; k<size; k+=blockSize) 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // partition the matrix: 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // A00 | - | - 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // lu = A10 | A11 | - 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // A20 | A21 | A22 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index bs = (std::min)(blockSize, size-k); 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rs = size - k - bs; 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<MatrixType,Dynamic,Dynamic> A11(m,k, k, bs,bs); 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k, rs,bs); 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs); 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index ret; 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if((ret=unblocked(A11))>=0) return k+ret; 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21); 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return -1; 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType, typename VectorType> 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma) 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Eigen::internal::llt_rank_update_lower(mat, vec, sigma); 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> struct llt_inplace<Scalar, Upper> 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat) 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Transpose<MatrixType> matt(mat); 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return llt_inplace<Scalar, Lower>::unblocked(matt); 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat) 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Transpose<MatrixType> matt(mat); 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return llt_inplace<Scalar, Lower>::blocked(matt); 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType, typename VectorType> 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma) 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Transpose<MatrixType> matt(mat); 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma); 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> struct LLT_Traits<MatrixType,Lower> 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef const TriangularView<const MatrixType, Lower> MatrixL; 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU; 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline MatrixL getL(const MatrixType& m) { return m; } 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static bool inplace_decomposition(MatrixType& m) 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; } 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> struct LLT_Traits<MatrixType,Upper> 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL; 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef const TriangularView<const MatrixType, Upper> MatrixU; 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); } 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline MatrixU getU(const MatrixType& m) { return m; } 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static bool inplace_decomposition(MatrixType& m) 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; } 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns a reference to *this 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include TutorialLinAlgComputeTwice.cpp 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude TutorialLinAlgComputeTwice.out 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int _UpLo> 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathLLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a) 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(a.rows()==a.cols()); 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = a.rows(); 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix.resize(size, size); 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix = a; 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool ok = Traits::inplace_decomposition(m_matrix); 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = ok ? Success : NumericalIssue; 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Performs a rank one update (or dowdate) of the current decomposition. 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * If A = LL^* before the rank one update, 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * then after it we have LL^* = A + sigma * v v^* where \a v must be a vector 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * of same dimension. 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int _UpLo> 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename VectorType> 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathLLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma) 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType); 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(v.size()==m_matrix.cols()); 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized); 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0) 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = NumericalIssue; 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = Success; 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int UpLo, typename Rhs> 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct solve_retval<LLT<_MatrixType, UpLo>, Rhs> 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : solve_retval_base<LLT<_MatrixType, UpLo>, Rhs> 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef LLT<_MatrixType,UpLo> LLTType; 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs) 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Dest> void evalTo(Dest& dst) const 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dst = rhs(); 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dec().solveInPlace(dst); 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \internal use x = llt_object.solve(x); 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This is the \em in-place version of solve(). 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param bAndX represents both the right-hand side matrix b and result x. 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD. 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This version avoids a copy when the right hand side matrix b is not 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * needed anymore. 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa LLT::solve(), MatrixBase::llt() 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int _UpLo> 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LLT is not initialized."); 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_matrix.rows()==bAndX.rows()); 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matrixL().solveInPlace(bAndX); 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matrixU().solveInPlace(bAndX); 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \returns the matrix represented by the decomposition, 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * i.e., it returns the product: L L^*. 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function is provided for debug purpose. */ 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int _UpLo> 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathMatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LLT is not initialized."); 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return matrixL() * matrixL().adjoint().toDenseMatrix(); 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \cholesky_module 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns the LLT decomposition of \c *this 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathinline const LLT<typename MatrixBase<Derived>::PlainObject> 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathMatrixBase<Derived>::llt() const 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return LLT<PlainObject>(derived()); 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \cholesky_module 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns the LLT decomposition of \c *this 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, unsigned int UpLo> 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathinline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo> 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathSelfAdjointView<MatrixType, UpLo>::llt() const 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return LLT<PlainObject,UpLo>(m_matrix); 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_LLT_H 489