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{ 1937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::ColXpr ColXpr; 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename internal::remove_all<ColXpr>::type ColXprCleaned; 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename ColXprCleaned::SegmentReturnType ColXprSegment; 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar,Dynamic,1> TempVectorType; 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename TempVectorType::SegmentReturnType TempVecSegment; 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index n = mat.cols(); 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(mat.rows()==n && vec.size()==n); 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath TempVectorType temp; 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(sigma>0) 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // This version is based on Givens rotations. 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // It is faster than the other one below, but only works for updates, 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // i.e., for sigma > 0 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp = sqrt(sigma) * vec; 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(Index i=0; i<n; ++i) 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath JacobiRotation<Scalar> g; 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath g.makeGivens(mat(i,i), -temp(i), &mat(i,i)); 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index rs = n-i-1; 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(rs>0) 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColXprSegment x(mat.col(i).tail(rs)); 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath TempVecSegment y(temp.tail(rs)); 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath apply_rotation_in_the_plane(x, y, g); 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp = vec; 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar beta = 1; 2337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for(Index j=0; j<n; ++j) 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 2357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar Ljj = numext::real(mat.coeff(j,j)); 2367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar dj = numext::abs2(Ljj); 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar wj = temp.coeff(j); 2387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar swj2 = sigma*numext::abs2(wj); 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar gamma = dj*beta + swj2; 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar x = dj + swj2/beta; 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (x<=RealScalar(0)) 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return j; 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar nLjj = sqrt(x); 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat.coeffRef(j,j) = nLjj; 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath beta += swj2/dj; 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Update the terms of L 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rs = n-j-1; 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(rs) 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp.tail(rs) -= (wj/Ljj) * mat.col(j).tail(rs); 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(gamma != 0) 2547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mat.col(j).tail(rs) = (nLjj/Ljj) * mat.col(j).tail(rs) + (nLjj * sigma*numext::conj(wj)/gamma)*temp.tail(rs); 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return -1; 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> struct llt_inplace<Scalar, Lower> 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static typename MatrixType::Index unblocked(MatrixType& mat) 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 2677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(mat.rows()==mat.cols()); 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = mat.rows(); 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index k = 0; k < size; ++k) 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rs = size-k-1; // remaining size 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1); 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k); 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k); 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar x = numext::real(mat.coeff(k,k)); 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k>0) x -= A10.squaredNorm(); 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (x<=RealScalar(0)) 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return k; 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat.coeffRef(k,k) = x = sqrt(x); 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k>0 && rs>0) A21.noalias() -= A20 * A10.adjoint(); 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (rs>0) A21 *= RealScalar(1)/x; 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return -1; 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static typename MatrixType::Index blocked(MatrixType& m) 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m.rows()==m.cols()); 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index size = m.rows(); 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(size<32) 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return unblocked(m); 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index blockSize = size/8; 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath blockSize = (blockSize/16)*16; 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath blockSize = (std::min)((std::max)(blockSize,Index(8)), Index(128)); 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index k=0; k<size; k+=blockSize) 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // partition the matrix: 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // A00 | - | - 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // lu = A10 | A11 | - 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // A20 | A21 | A22 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index bs = (std::min)(blockSize, size-k); 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rs = size - k - bs; 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<MatrixType,Dynamic,Dynamic> A11(m,k, k, bs,bs); 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<MatrixType,Dynamic,Dynamic> A21(m,k+bs,k, rs,bs); 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<MatrixType,Dynamic,Dynamic> A22(m,k+bs,k+bs,rs,rs); 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index ret; 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if((ret=unblocked(A11))>=0) return k+ret; 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(rs>0) A11.adjoint().template triangularView<Upper>().template solveInPlace<OnTheRight>(A21); 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(rs>0) A22.template selfadjointView<Lower>().rankUpdate(A21,-1); // bottleneck 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return -1; 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType, typename VectorType> 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma) 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Eigen::internal::llt_rank_update_lower(mat, vec, sigma); 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> struct llt_inplace<Scalar, Upper> 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static EIGEN_STRONG_INLINE typename MatrixType::Index unblocked(MatrixType& mat) 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Transpose<MatrixType> matt(mat); 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return llt_inplace<Scalar, Lower>::unblocked(matt); 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType> 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static EIGEN_STRONG_INLINE typename MatrixType::Index blocked(MatrixType& mat) 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Transpose<MatrixType> matt(mat); 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return llt_inplace<Scalar, Lower>::blocked(matt); 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType, typename VectorType> 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static typename MatrixType::Index rankUpdate(MatrixType& mat, const VectorType& vec, const RealScalar& sigma) 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Transpose<MatrixType> matt(mat); 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return llt_inplace<Scalar, Lower>::rankUpdate(matt, vec.conjugate(), sigma); 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> struct LLT_Traits<MatrixType,Lower> 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef const TriangularView<const MatrixType, Lower> MatrixL; 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef const TriangularView<const typename MatrixType::AdjointReturnType, Upper> MatrixU; 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline MatrixL getL(const MatrixType& m) { return m; } 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static bool inplace_decomposition(MatrixType& m) 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { return llt_inplace<typename MatrixType::Scalar, Lower>::blocked(m)==-1; } 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> struct LLT_Traits<MatrixType,Upper> 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef const TriangularView<const typename MatrixType::AdjointReturnType, Lower> MatrixL; 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef const TriangularView<const MatrixType, Upper> MatrixU; 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); } 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline MatrixU getU(const MatrixType& m) { return m; } 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static bool inplace_decomposition(MatrixType& m) 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { return llt_inplace<typename MatrixType::Scalar, Upper>::blocked(m)==-1; } 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Computes / recomputes the Cholesky decomposition A = LL^* = U^*U of \a matrix 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns a reference to *this 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include TutorialLinAlgComputeTwice.cpp 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude TutorialLinAlgComputeTwice.out 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int _UpLo> 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathLLT<MatrixType,_UpLo>& LLT<MatrixType,_UpLo>::compute(const MatrixType& a) 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(a.rows()==a.cols()); 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = a.rows(); 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix.resize(size, size); 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix = a; 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool ok = Traits::inplace_decomposition(m_matrix); 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = ok ? Success : NumericalIssue; 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Performs a rank one update (or dowdate) of the current decomposition. 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * If A = LL^* before the rank one update, 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * then after it we have LL^* = A + sigma * v v^* where \a v must be a vector 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * of same dimension. 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int _UpLo> 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename VectorType> 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathLLT<_MatrixType,_UpLo> LLT<_MatrixType,_UpLo>::rankUpdate(const VectorType& v, const RealScalar& sigma) 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_STATIC_ASSERT_VECTOR_ONLY(VectorType); 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(v.size()==m_matrix.cols()); 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized); 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(internal::llt_inplace<typename MatrixType::Scalar, UpLo>::rankUpdate(m_matrix,v,sigma)>=0) 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = NumericalIssue; 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_info = Success; 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int UpLo, typename Rhs> 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct solve_retval<LLT<_MatrixType, UpLo>, Rhs> 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : solve_retval_base<LLT<_MatrixType, UpLo>, Rhs> 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef LLT<_MatrixType,UpLo> LLTType; 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_MAKE_SOLVE_HELPERS(LLTType,Rhs) 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Dest> void evalTo(Dest& dst) const 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dst = rhs(); 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dec().solveInPlace(dst); 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \internal use x = llt_object.solve(x); 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This is the \em in-place version of solve(). 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param bAndX represents both the right-hand side matrix b and result x. 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD. 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This version avoids a copy when the right hand side matrix b is not 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * needed anymore. 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa LLT::solve(), MatrixBase::llt() 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int _UpLo> 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid LLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LLT is not initialized."); 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_matrix.rows()==bAndX.rows()); 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matrixL().solveInPlace(bAndX); 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matrixU().solveInPlace(bAndX); 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \returns the matrix represented by the decomposition, 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * i.e., it returns the product: L L^*. 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function is provided for debug purpose. */ 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int _UpLo> 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathMatrixType LLT<MatrixType,_UpLo>::reconstructedMatrix() const 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LLT is not initialized."); 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return matrixL() * matrixL().adjoint().toDenseMatrix(); 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \cholesky_module 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns the LLT decomposition of \c *this 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathinline const LLT<typename MatrixBase<Derived>::PlainObject> 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathMatrixBase<Derived>::llt() const 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return LLT<PlainObject>(derived()); 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \cholesky_module 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns the LLT decomposition of \c *this 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, unsigned int UpLo> 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathinline const LLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo> 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathSelfAdjointView<MatrixType, UpLo>::llt() const 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return LLT<PlainObject,UpLo>(m_matrix); 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 489c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_LLT_H 491