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-2011 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Keir Mierle <mierle@gmail.com> 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2011 Timothy E. Holy <tim.holy@gmail.com > 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_LDLT_H 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_LDLT_H 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename MatrixType, int UpLo> struct LDLT_Traits; 207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // PositiveSemiDef means positive semi-definite and non-zero; same for NegativeSemiDef 227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez enum SignMatrix { PositiveSemiDef, NegativeSemiDef, ZeroSign, Indefinite }; 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \ingroup Cholesky_Module 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class LDLT 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Robust Cholesky decomposition of a matrix with pivoting 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param MatrixType the type of the matrix of which to compute the LDL^T Cholesky decomposition 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper. 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The other triangular part won't be read. 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Perform a robust Cholesky decomposition of a positive semidefinite or negative semidefinite 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix \f$ A \f$ such that \f$ A = P^TLDL^*P \f$, where P is a permutation matrix, L 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * is lower triangular with a unit diagonal and D is a diagonal matrix. 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The decomposition uses pivoting to ensure stability, so that L will have 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * zeros in the bottom right rank(A) - n submatrix. Avoiding the square root 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * on D also stabilizes the computation. 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * decomposition to determine whether a system of equations has a solution. 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::ldlt(), class LLT 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int _UpLo> class LDLT 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowsAtCompileTime = MatrixType::RowsAtCompileTime, 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColsAtCompileTime = MatrixType::ColsAtCompileTime, 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Options = MatrixType::Options & ~RowMajorBit, // these are the options for the TmpMatrixType, we need a ColMajor matrix here! 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath UpLo = _UpLo 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 typedef Matrix<Scalar, RowsAtCompileTime, 1, Options, MaxRowsAtCompileTime, 1> TmpMatrixType; 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType; 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType; 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef internal::LDLT_Traits<MatrixType,UpLo> Traits; 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Default Constructor. 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The default constructor is useful in cases in which the user intends to 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * perform decompositions via LDLT::compute(const MatrixType&). 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LDLT() 767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : m_matrix(), 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_transpositions(), 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_sign(internal::ZeroSign), 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_isInitialized(false) 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez {} 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Default Constructor with memory preallocation 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Like the default constructor but with preallocation of the internal data 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * according to the specified problem \a size. 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa LDLT() 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LDLT(Index size) 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_matrix(size, size), 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_transpositions(size), 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_temporary(size), 927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_sign(internal::ZeroSign), 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false) 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor with decomposition 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This calculates the decomposition for the input \a matrix. 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa LDLT(Index size) 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LDLT(const MatrixType& matrix) 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_matrix(matrix.rows(), matrix.cols()), 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_transpositions(matrix.rows()), 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_temporary(matrix.rows()), 1057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_sign(internal::ZeroSign), 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false) 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix); 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** Clear any existing decomposition 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa rankUpdate(w,sigma) 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void setZero() 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = false; 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns a view of the upper triangular matrix U */ 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline typename Traits::MatrixU matrixU() const 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LDLT is not initialized."); 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Traits::getU(m_matrix); 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns a view of the lower triangular matrix L */ 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline typename Traits::MatrixL matrixL() const 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LDLT is not initialized."); 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Traits::getL(m_matrix); 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the permutation matrix P as a transposition sequence. 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const TranspositionType& transpositionsP() const 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LDLT is not initialized."); 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_transpositions; 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the coefficients of the diagonal matrix D */ 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Diagonal<const MatrixType> vectorD() const 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LDLT is not initialized."); 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_matrix.diagonal(); 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns true if the matrix is positive (semidefinite) */ 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline bool isPositive() const 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LDLT is not initialized."); 1527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_sign == internal::PositiveSemiDef || m_sign == internal::ZeroSign; 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath #ifdef EIGEN2_SUPPORT 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline bool isPositiveDefinite() const 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return isPositive(); 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath #endif 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns true if the matrix is negative (semidefinite) */ 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline bool isNegative(void) const 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LDLT is not initialized."); 1667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_sign == internal::NegativeSemiDef || m_sign == internal::ZeroSign; 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns a solution x of \f$ A x = b \f$ using the current decomposition of A. 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function also supports in-place solves using the syntax <tt>x = decompositionObject.solve(x)</tt> . 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \note_about_checking_solutions 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * More precisely, this method solves \f$ A x = b \f$ using the decomposition \f$ A = P^T L D L^* P \f$ 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * by solving the systems \f$ P^T y_1 = b \f$, \f$ L y_2 = y_1 \f$, \f$ D y_3 = y_2 \f$, 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ L^* y_4 = y_3 \f$ and \f$ P x = y_4 \f$ in succession. If the matrix \f$ A \f$ is singular, then 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ D \f$ will also be singular (all the other matrices are invertible). In that case, the 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * least-square solution of \f$ D y_3 = y_2 \f$ is computed. This does not mean that this function 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * computes the least-square solution of \f$ A x = b \f$ is \f$ A \f$ is singular. 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::ldlt() 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Rhs> 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const internal::solve_retval<LDLT, Rhs> 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath solve(const MatrixBase<Rhs>& b) const 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LDLT is not initialized."); 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_matrix.rows()==b.rows() 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath && "LDLT::solve(): invalid number of rows of the right hand side matrix b"); 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return internal::solve_retval<LDLT, Rhs>(*this, b.derived()); 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath #ifdef EIGEN2_SUPPORT 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename OtherDerived, typename ResultType> 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *result = this->solve(b); 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath #endif 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Derived> 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool solveInPlace(MatrixBase<Derived> &bAndX) const; 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath LDLT& compute(const MatrixType& matrix); 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename Derived> 2097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez LDLT& rankUpdate(const MatrixBase<Derived>& w, const RealScalar& alpha=1); 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \returns the internal LDLT decomposition matrix 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * TODO: document the storage layout 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline const MatrixType& matrixLDLT() const 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LDLT is not initialized."); 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_matrix; 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType reconstructedMatrix() const; 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index rows() const { return m_matrix.rows(); } 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline Index cols() const { return m_matrix.cols(); } 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Reports whether previous computation was successful. 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns \c Success if computation was succesful, 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \c NumericalIssue if the matrix.appears to be negative. 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo info() const 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LDLT is not initialized."); 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return Success; 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \internal 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Used to compute and store the Cholesky decomposition A = L D L^* = U^* D U. 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The strict upper part is used during the decomposition, the strict lower 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * part correspond to the coefficients of L (its diagonal is equal to 1 and 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * is not stored), and the diagonal entries correspond to D. 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType m_matrix; 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath TranspositionType m_transpositions; 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath TmpMatrixType m_temporary; 2487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::SignMatrix m_sign; 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<int UpLo> struct ldlt_inplace; 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<> struct ldlt_inplace<Lower> 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType, typename TranspositionType, typename Workspace> 2597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign) 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 2617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(mat.rows()==mat.cols()); 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = mat.rows(); 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (size <= 1) 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath transpositions.setIdentity(); 2717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (numext::real(mat.coeff(0,0)) > 0) sign = PositiveSemiDef; 2727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (numext::real(mat.coeff(0,0)) < 0) sign = NegativeSemiDef; 2737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else sign = ZeroSign; 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index k = 0; k < size; ++k) 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Find largest diagonal element 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index index_of_biggest_in_corner; 2817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner); 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath index_of_biggest_in_corner += k; 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath transpositions.coeffRef(k) = index_of_biggest_in_corner; 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(k != index_of_biggest_in_corner) 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // apply the transposition while taking care to consider only 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // the lower triangular part 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k)); 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s)); 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner)); 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(int i=k+1;i<index_of_biggest_in_corner;++i) 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar tmp = mat.coeffRef(i,k); 2967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mat.coeffRef(i,k) = numext::conj(mat.coeffRef(index_of_biggest_in_corner,i)); 2977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mat.coeffRef(index_of_biggest_in_corner,i) = numext::conj(tmp); 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(NumTraits<Scalar>::IsComplex) 3007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mat.coeffRef(index_of_biggest_in_corner,k) = numext::conj(mat.coeff(index_of_biggest_in_corner,k)); 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // partition the matrix: 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // A00 | - | - 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // lu = A10 | A11 | - 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // A20 | A21 | A22 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rs = size - k - 1; 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1); 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k); 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k); 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(k>0) 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez temp.head(k) = mat.diagonal().real().head(k).asDiagonal() * A10.adjoint(); 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat.coeffRef(k,k) -= (A10 * temp.head(k)).value(); 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(rs>0) 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath A21.noalias() -= A20 * temp.head(k); 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 3197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // In some previous versions of Eigen (e.g., 3.2.1), the scaling was omitted if the pivot 3217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // was smaller than the cutoff value. However, soince LDLT is not rank-revealing 3227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // we should only make sure we do not introduce INF or NaN values. 3237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // LAPACK also uses 0 as the cutoff value. 3247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar realAkk = numext::real(mat.coeffRef(k,k)); 3257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if((rs>0) && (abs(realAkk) > RealScalar(0))) 3267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez A21 /= realAkk; 3277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (sign == PositiveSemiDef) { 3297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (realAkk < 0) sign = Indefinite; 3307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } else if (sign == NegativeSemiDef) { 3317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (realAkk > 0) sign = Indefinite; 3327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } else if (sign == ZeroSign) { 3337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (realAkk > 0) sign = PositiveSemiDef; 3347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (realAkk < 0) sign = NegativeSemiDef; 3357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Reference for the algorithm: Davis and Hager, "Multiple Rank 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Modifications of a Sparse Cholesky Factorization" (Algorithm 1) 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Trivial rearrangements of their computations (Timothy E. Holy) 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // allow their algorithm to work for rank-1 updates even if the 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // original matrix is not of full rank. 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Here only rank-1 updates are implemented, to reduce the 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // requirement for intermediate storage and improve accuracy 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType, typename WDerived> 3497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static bool updateInPlace(MatrixType& mat, MatrixBase<WDerived>& w, const typename MatrixType::RealScalar& sigma=1) 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using numext::isfinite; 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::RealScalar RealScalar; 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = mat.rows(); 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(mat.cols() == size && w.size()==size); 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar alpha = 1; 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Apply the update 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j = 0; j < size; j++) 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Check for termination due to an original decomposition of low-rank 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (!(isfinite)(alpha)) 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath break; 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Update the diagonal terms 3697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar dj = numext::real(mat.coeff(j,j)); 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar wj = w.coeff(j); 3717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar swj2 = sigma*numext::abs2(wj); 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar gamma = dj*alpha + swj2; 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath mat.coeffRef(j,j) += swj2/alpha; 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath alpha += swj2/dj; 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Update the terms of L 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rs = size-j-1; 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath w.tail(rs) -= wj * mat.col(j).tail(rs); 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(gamma != 0) 3827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez mat.col(j).tail(rs) += (sigma*numext::conj(wj)/gamma)*w.tail(rs); 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType> 3887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, const typename MatrixType::RealScalar& sigma=1) 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Apply the permutation to the input w 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp = transpositions * w; 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return ldlt_inplace<Lower>::updateInPlace(mat,tmp,sigma); 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<> struct ldlt_inplace<Upper> 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType, typename TranspositionType, typename Workspace> 4007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign) 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Transpose<MatrixType> matt(mat); 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign); 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType> 4077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w, const typename MatrixType::RealScalar& sigma=1) 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Transpose<MatrixType> matt(mat); 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return ldlt_inplace<Lower>::update(matt, transpositions, tmp, w.conjugate(), sigma); 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> struct LDLT_Traits<MatrixType,Lower> 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef const TriangularView<const MatrixType, UnitLower> MatrixL; 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitUpper> MatrixU; 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline MatrixL getL(const MatrixType& m) { return m; } 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> struct LDLT_Traits<MatrixType,Upper> 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitLower> MatrixL; 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef const TriangularView<const MatrixType, UnitUpper> MatrixU; 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline MatrixL getL(const MatrixType& m) { return m.adjoint(); } 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static inline MatrixU getU(const MatrixType& m) { return m; } 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int _UpLo> 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathLDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const MatrixType& a) 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(a.rows()==a.cols()); 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = a.rows(); 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix = a; 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_transpositions.resize(size); 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = false; 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_temporary.resize(size); 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign); 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Update the LDLT decomposition: given A = L D L^T, efficiently compute the decomposition of A + sigma w w^T. 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param w a vector to be incorporated into the decomposition. 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param sigma a scalar, +1 for updates and -1 for "downdates," which correspond to removing previously-added column vectors. Optional; default value is +1. 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa setZero() 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int _UpLo> 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 4597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezLDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::rankUpdate(const MatrixBase<Derived>& w, const typename NumTraits<typename MatrixType::Scalar>::Real& sigma) 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = w.rows(); 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_isInitialized) 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_matrix.rows()==size); 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix.resize(size,size); 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix.setZero(); 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_transpositions.resize(size); 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = 0; i < size; i++) 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_transpositions.coeffRef(i) = i; 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_temporary.resize(size); 4747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_sign = sigma>=0 ? internal::PositiveSemiDef : internal::NegativeSemiDef; 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::ldlt_inplace<UpLo>::update(m_matrix, m_transpositions, m_temporary, w, sigma); 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType, int _UpLo, typename Rhs> 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct solve_retval<LDLT<_MatrixType,_UpLo>, Rhs> 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : solve_retval_base<LDLT<_MatrixType,_UpLo>, Rhs> 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef LDLT<_MatrixType,_UpLo> LDLTType; 489c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_MAKE_SOLVE_HELPERS(LDLTType,Rhs) 490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 491c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template<typename Dest> void evalTo(Dest& dst) const 492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(rhs().rows() == dec().matrixLDLT().rows()); 494c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // dst = P b 495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dst = dec().transpositionsP() * rhs(); 496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 497c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // dst = L^-1 (P b) 498c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dec().matrixL().solveInPlace(dst); 499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // dst = D^-1 (L^-1 P b) 501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // more precisely, use pseudo-inverse of D (see bug 241) 502c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using std::abs; 503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using std::max; 504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename LDLTType::MatrixType MatrixType; 505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename LDLTType::Scalar Scalar; 506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename LDLTType::RealScalar RealScalar; 5077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const typename Diagonal<const MatrixType>::RealReturnType vectorD(dec().vectorD()); 5087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // In some previous versions, tolerance was set to the max of 1/highest and the maximal diagonal entry * epsilon 5097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // as motivated by LAPACK's xGELSS: 5107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // RealScalar tolerance = (max)(vectorD.array().abs().maxCoeff() *NumTraits<RealScalar>::epsilon(),RealScalar(1) / NumTraits<RealScalar>::highest()); 5117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // However, LDLT is not rank revealing, and so adjusting the tolerance wrt to the highest 5127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // diagonal element is not well justified and to numerical issues in some cases. 5137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Moreover, Lapack's xSYTRS routines use 0 for the tolerance. 5147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealScalar tolerance = RealScalar(1) / NumTraits<RealScalar>::highest(); 5157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 516c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = 0; i < vectorD.size(); ++i) { 517c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(abs(vectorD(i)) > tolerance) 5187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dst.row(i) /= vectorD(i); 519c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 5207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez dst.row(i).setZero(); 521c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 522c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 523c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // dst = L^-T (D^-1 L^-1 P b) 524c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dec().matrixU().solveInPlace(dst); 525c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // dst = P^-1 (L^-T D^-1 L^-1 P b) = A^-1 b 527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath dst = dec().transpositionsP().transpose() * dst; 528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 529c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 530c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 532c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \internal use x = ldlt_object.solve(x); 533c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 534c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This is the \em in-place version of solve(). 535c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param bAndX represents both the right-hand side matrix b and result x. 537c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD. 539c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This version avoids a copy when the right hand side matrix b is not 541c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * needed anymore. 542c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 543c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa LDLT::solve(), MatrixBase::ldlt() 544c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 545c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType,int _UpLo> 546c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 547c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathbool LDLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const 548c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 549c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LDLT is not initialized."); 5507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_matrix.rows() == bAndX.rows()); 551c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 552c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bAndX = this->solve(bAndX); 553c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 554c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return true; 555c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 556c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 557c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \returns the matrix represented by the decomposition, 558c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * i.e., it returns the product: P^T L D L^* P. 559c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function is provided for debug purpose. */ 560c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, int _UpLo> 561c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathMatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const 562c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 563c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "LDLT is not initialized."); 564c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = m_matrix.rows(); 565c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType res(size,size); 566c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 567c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // P 568c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath res.setIdentity(); 569c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath res = transpositionsP() * res; 570c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // L^* P 571c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath res = matrixU() * res; 572c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // D(L^*P) 5737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez res = vectorD().real().asDiagonal() * res; 574c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // L(DL^*P) 575c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath res = matrixL() * res; 576c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // P^T (LDL^*P) 577c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath res = transpositionsP().transpose() * res; 578c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 579c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return res; 580c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 581c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 582c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \cholesky_module 583c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns the Cholesky decomposition with full pivoting without square root of \c *this 584c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 585c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType, unsigned int UpLo> 586c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathinline const LDLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo> 587c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathSelfAdjointView<MatrixType, UpLo>::ldlt() const 588c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 589c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return LDLT<PlainObject,UpLo>(m_matrix); 590c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 591c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 592c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \cholesky_module 593c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns the Cholesky decomposition with full pivoting without square root of \c *this 594c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 595c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 596c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathinline const LDLT<typename MatrixBase<Derived>::PlainObject> 597c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathMatrixBase<Derived>::ldlt() const 598c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 599c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return LDLT<PlainObject>(derived()); 600c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 601c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 602c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 603c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 604c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_LDLT_H 605