LDLT.h (24934B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr> 5 // Copyright (C) 2009 Keir Mierle <mierle@gmail.com> 6 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> 7 // Copyright (C) 2011 Timothy E. Holy <tim.holy@gmail.com > 8 // 9 // This Source Code Form is subject to the terms of the Mozilla 10 // Public License v. 2.0. If a copy of the MPL was not distributed 11 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 12 13 #ifndef EIGEN_LDLT_H 14 #define EIGEN_LDLT_H 15 16 namespace Eigen { 17 18 namespace internal { 19 template<typename _MatrixType, int _UpLo> struct traits<LDLT<_MatrixType, _UpLo> > 20 : traits<_MatrixType> 21 { 22 typedef MatrixXpr XprKind; 23 typedef SolverStorage StorageKind; 24 typedef int StorageIndex; 25 enum { Flags = 0 }; 26 }; 27 28 template<typename MatrixType, int UpLo> struct LDLT_Traits; 29 30 // PositiveSemiDef means positive semi-definite and non-zero; same for NegativeSemiDef 31 enum SignMatrix { PositiveSemiDef, NegativeSemiDef, ZeroSign, Indefinite }; 32 } 33 34 /** \ingroup Cholesky_Module 35 * 36 * \class LDLT 37 * 38 * \brief Robust Cholesky decomposition of a matrix with pivoting 39 * 40 * \tparam _MatrixType the type of the matrix of which to compute the LDL^T Cholesky decomposition 41 * \tparam _UpLo the triangular part that will be used for the decompositon: Lower (default) or Upper. 42 * The other triangular part won't be read. 43 * 44 * Perform a robust Cholesky decomposition of a positive semidefinite or negative semidefinite 45 * matrix \f$ A \f$ such that \f$ A = P^TLDL^*P \f$, where P is a permutation matrix, L 46 * is lower triangular with a unit diagonal and D is a diagonal matrix. 47 * 48 * The decomposition uses pivoting to ensure stability, so that D will have 49 * zeros in the bottom right rank(A) - n submatrix. Avoiding the square root 50 * on D also stabilizes the computation. 51 * 52 * Remember that Cholesky decompositions are not rank-revealing. Also, do not use a Cholesky 53 * decomposition to determine whether a system of equations has a solution. 54 * 55 * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. 56 * 57 * \sa MatrixBase::ldlt(), SelfAdjointView::ldlt(), class LLT 58 */ 59 template<typename _MatrixType, int _UpLo> class LDLT 60 : public SolverBase<LDLT<_MatrixType, _UpLo> > 61 { 62 public: 63 typedef _MatrixType MatrixType; 64 typedef SolverBase<LDLT> Base; 65 friend class SolverBase<LDLT>; 66 67 EIGEN_GENERIC_PUBLIC_INTERFACE(LDLT) 68 enum { 69 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 70 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime, 71 UpLo = _UpLo 72 }; 73 typedef Matrix<Scalar, RowsAtCompileTime, 1, 0, MaxRowsAtCompileTime, 1> TmpMatrixType; 74 75 typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType; 76 typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType; 77 78 typedef internal::LDLT_Traits<MatrixType,UpLo> Traits; 79 80 /** \brief Default Constructor. 81 * 82 * The default constructor is useful in cases in which the user intends to 83 * perform decompositions via LDLT::compute(const MatrixType&). 84 */ 85 LDLT() 86 : m_matrix(), 87 m_transpositions(), 88 m_sign(internal::ZeroSign), 89 m_isInitialized(false) 90 {} 91 92 /** \brief Default Constructor with memory preallocation 93 * 94 * Like the default constructor but with preallocation of the internal data 95 * according to the specified problem \a size. 96 * \sa LDLT() 97 */ 98 explicit LDLT(Index size) 99 : m_matrix(size, size), 100 m_transpositions(size), 101 m_temporary(size), 102 m_sign(internal::ZeroSign), 103 m_isInitialized(false) 104 {} 105 106 /** \brief Constructor with decomposition 107 * 108 * This calculates the decomposition for the input \a matrix. 109 * 110 * \sa LDLT(Index size) 111 */ 112 template<typename InputType> 113 explicit LDLT(const EigenBase<InputType>& matrix) 114 : m_matrix(matrix.rows(), matrix.cols()), 115 m_transpositions(matrix.rows()), 116 m_temporary(matrix.rows()), 117 m_sign(internal::ZeroSign), 118 m_isInitialized(false) 119 { 120 compute(matrix.derived()); 121 } 122 123 /** \brief Constructs a LDLT factorization from a given matrix 124 * 125 * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. 126 * 127 * \sa LDLT(const EigenBase&) 128 */ 129 template<typename InputType> 130 explicit LDLT(EigenBase<InputType>& matrix) 131 : m_matrix(matrix.derived()), 132 m_transpositions(matrix.rows()), 133 m_temporary(matrix.rows()), 134 m_sign(internal::ZeroSign), 135 m_isInitialized(false) 136 { 137 compute(matrix.derived()); 138 } 139 140 /** Clear any existing decomposition 141 * \sa rankUpdate(w,sigma) 142 */ 143 void setZero() 144 { 145 m_isInitialized = false; 146 } 147 148 /** \returns a view of the upper triangular matrix U */ 149 inline typename Traits::MatrixU matrixU() const 150 { 151 eigen_assert(m_isInitialized && "LDLT is not initialized."); 152 return Traits::getU(m_matrix); 153 } 154 155 /** \returns a view of the lower triangular matrix L */ 156 inline typename Traits::MatrixL matrixL() const 157 { 158 eigen_assert(m_isInitialized && "LDLT is not initialized."); 159 return Traits::getL(m_matrix); 160 } 161 162 /** \returns the permutation matrix P as a transposition sequence. 163 */ 164 inline const TranspositionType& transpositionsP() const 165 { 166 eigen_assert(m_isInitialized && "LDLT is not initialized."); 167 return m_transpositions; 168 } 169 170 /** \returns the coefficients of the diagonal matrix D */ 171 inline Diagonal<const MatrixType> vectorD() const 172 { 173 eigen_assert(m_isInitialized && "LDLT is not initialized."); 174 return m_matrix.diagonal(); 175 } 176 177 /** \returns true if the matrix is positive (semidefinite) */ 178 inline bool isPositive() const 179 { 180 eigen_assert(m_isInitialized && "LDLT is not initialized."); 181 return m_sign == internal::PositiveSemiDef || m_sign == internal::ZeroSign; 182 } 183 184 /** \returns true if the matrix is negative (semidefinite) */ 185 inline bool isNegative(void) const 186 { 187 eigen_assert(m_isInitialized && "LDLT is not initialized."); 188 return m_sign == internal::NegativeSemiDef || m_sign == internal::ZeroSign; 189 } 190 191 #ifdef EIGEN_PARSED_BY_DOXYGEN 192 /** \returns a solution x of \f$ A x = b \f$ using the current decomposition of A. 193 * 194 * This function also supports in-place solves using the syntax <tt>x = decompositionObject.solve(x)</tt> . 195 * 196 * \note_about_checking_solutions 197 * 198 * More precisely, this method solves \f$ A x = b \f$ using the decomposition \f$ A = P^T L D L^* P \f$ 199 * 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$, 200 * \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 201 * \f$ D \f$ will also be singular (all the other matrices are invertible). In that case, the 202 * least-square solution of \f$ D y_3 = y_2 \f$ is computed. This does not mean that this function 203 * computes the least-square solution of \f$ A x = b \f$ if \f$ A \f$ is singular. 204 * 205 * \sa MatrixBase::ldlt(), SelfAdjointView::ldlt() 206 */ 207 template<typename Rhs> 208 inline const Solve<LDLT, Rhs> 209 solve(const MatrixBase<Rhs>& b) const; 210 #endif 211 212 template<typename Derived> 213 bool solveInPlace(MatrixBase<Derived> &bAndX) const; 214 215 template<typename InputType> 216 LDLT& compute(const EigenBase<InputType>& matrix); 217 218 /** \returns an estimate of the reciprocal condition number of the matrix of 219 * which \c *this is the LDLT decomposition. 220 */ 221 RealScalar rcond() const 222 { 223 eigen_assert(m_isInitialized && "LDLT is not initialized."); 224 return internal::rcond_estimate_helper(m_l1_norm, *this); 225 } 226 227 template <typename Derived> 228 LDLT& rankUpdate(const MatrixBase<Derived>& w, const RealScalar& alpha=1); 229 230 /** \returns the internal LDLT decomposition matrix 231 * 232 * TODO: document the storage layout 233 */ 234 inline const MatrixType& matrixLDLT() const 235 { 236 eigen_assert(m_isInitialized && "LDLT is not initialized."); 237 return m_matrix; 238 } 239 240 MatrixType reconstructedMatrix() const; 241 242 /** \returns the adjoint of \c *this, that is, a const reference to the decomposition itself as the underlying matrix is self-adjoint. 243 * 244 * This method is provided for compatibility with other matrix decompositions, thus enabling generic code such as: 245 * \code x = decomposition.adjoint().solve(b) \endcode 246 */ 247 const LDLT& adjoint() const { return *this; }; 248 249 EIGEN_DEVICE_FUNC inline EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT { return m_matrix.rows(); } 250 EIGEN_DEVICE_FUNC inline EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT { return m_matrix.cols(); } 251 252 /** \brief Reports whether previous computation was successful. 253 * 254 * \returns \c Success if computation was successful, 255 * \c NumericalIssue if the factorization failed because of a zero pivot. 256 */ 257 ComputationInfo info() const 258 { 259 eigen_assert(m_isInitialized && "LDLT is not initialized."); 260 return m_info; 261 } 262 263 #ifndef EIGEN_PARSED_BY_DOXYGEN 264 template<typename RhsType, typename DstType> 265 void _solve_impl(const RhsType &rhs, DstType &dst) const; 266 267 template<bool Conjugate, typename RhsType, typename DstType> 268 void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const; 269 #endif 270 271 protected: 272 273 static void check_template_parameters() 274 { 275 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 276 } 277 278 /** \internal 279 * Used to compute and store the Cholesky decomposition A = L D L^* = U^* D U. 280 * The strict upper part is used during the decomposition, the strict lower 281 * part correspond to the coefficients of L (its diagonal is equal to 1 and 282 * is not stored), and the diagonal entries correspond to D. 283 */ 284 MatrixType m_matrix; 285 RealScalar m_l1_norm; 286 TranspositionType m_transpositions; 287 TmpMatrixType m_temporary; 288 internal::SignMatrix m_sign; 289 bool m_isInitialized; 290 ComputationInfo m_info; 291 }; 292 293 namespace internal { 294 295 template<int UpLo> struct ldlt_inplace; 296 297 template<> struct ldlt_inplace<Lower> 298 { 299 template<typename MatrixType, typename TranspositionType, typename Workspace> 300 static bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign) 301 { 302 using std::abs; 303 typedef typename MatrixType::Scalar Scalar; 304 typedef typename MatrixType::RealScalar RealScalar; 305 typedef typename TranspositionType::StorageIndex IndexType; 306 eigen_assert(mat.rows()==mat.cols()); 307 const Index size = mat.rows(); 308 bool found_zero_pivot = false; 309 bool ret = true; 310 311 if (size <= 1) 312 { 313 transpositions.setIdentity(); 314 if(size==0) sign = ZeroSign; 315 else if (numext::real(mat.coeff(0,0)) > static_cast<RealScalar>(0) ) sign = PositiveSemiDef; 316 else if (numext::real(mat.coeff(0,0)) < static_cast<RealScalar>(0)) sign = NegativeSemiDef; 317 else sign = ZeroSign; 318 return true; 319 } 320 321 for (Index k = 0; k < size; ++k) 322 { 323 // Find largest diagonal element 324 Index index_of_biggest_in_corner; 325 mat.diagonal().tail(size-k).cwiseAbs().maxCoeff(&index_of_biggest_in_corner); 326 index_of_biggest_in_corner += k; 327 328 transpositions.coeffRef(k) = IndexType(index_of_biggest_in_corner); 329 if(k != index_of_biggest_in_corner) 330 { 331 // apply the transposition while taking care to consider only 332 // the lower triangular part 333 Index s = size-index_of_biggest_in_corner-1; // trailing size after the biggest element 334 mat.row(k).head(k).swap(mat.row(index_of_biggest_in_corner).head(k)); 335 mat.col(k).tail(s).swap(mat.col(index_of_biggest_in_corner).tail(s)); 336 std::swap(mat.coeffRef(k,k),mat.coeffRef(index_of_biggest_in_corner,index_of_biggest_in_corner)); 337 for(Index i=k+1;i<index_of_biggest_in_corner;++i) 338 { 339 Scalar tmp = mat.coeffRef(i,k); 340 mat.coeffRef(i,k) = numext::conj(mat.coeffRef(index_of_biggest_in_corner,i)); 341 mat.coeffRef(index_of_biggest_in_corner,i) = numext::conj(tmp); 342 } 343 if(NumTraits<Scalar>::IsComplex) 344 mat.coeffRef(index_of_biggest_in_corner,k) = numext::conj(mat.coeff(index_of_biggest_in_corner,k)); 345 } 346 347 // partition the matrix: 348 // A00 | - | - 349 // lu = A10 | A11 | - 350 // A20 | A21 | A22 351 Index rs = size - k - 1; 352 Block<MatrixType,Dynamic,1> A21(mat,k+1,k,rs,1); 353 Block<MatrixType,1,Dynamic> A10(mat,k,0,1,k); 354 Block<MatrixType,Dynamic,Dynamic> A20(mat,k+1,0,rs,k); 355 356 if(k>0) 357 { 358 temp.head(k) = mat.diagonal().real().head(k).asDiagonal() * A10.adjoint(); 359 mat.coeffRef(k,k) -= (A10 * temp.head(k)).value(); 360 if(rs>0) 361 A21.noalias() -= A20 * temp.head(k); 362 } 363 364 // In some previous versions of Eigen (e.g., 3.2.1), the scaling was omitted if the pivot 365 // was smaller than the cutoff value. However, since LDLT is not rank-revealing 366 // we should only make sure that we do not introduce INF or NaN values. 367 // Remark that LAPACK also uses 0 as the cutoff value. 368 RealScalar realAkk = numext::real(mat.coeffRef(k,k)); 369 bool pivot_is_valid = (abs(realAkk) > RealScalar(0)); 370 371 if(k==0 && !pivot_is_valid) 372 { 373 // The entire diagonal is zero, there is nothing more to do 374 // except filling the transpositions, and checking whether the matrix is zero. 375 sign = ZeroSign; 376 for(Index j = 0; j<size; ++j) 377 { 378 transpositions.coeffRef(j) = IndexType(j); 379 ret = ret && (mat.col(j).tail(size-j-1).array()==Scalar(0)).all(); 380 } 381 return ret; 382 } 383 384 if((rs>0) && pivot_is_valid) 385 A21 /= realAkk; 386 else if(rs>0) 387 ret = ret && (A21.array()==Scalar(0)).all(); 388 389 if(found_zero_pivot && pivot_is_valid) ret = false; // factorization failed 390 else if(!pivot_is_valid) found_zero_pivot = true; 391 392 if (sign == PositiveSemiDef) { 393 if (realAkk < static_cast<RealScalar>(0)) sign = Indefinite; 394 } else if (sign == NegativeSemiDef) { 395 if (realAkk > static_cast<RealScalar>(0)) sign = Indefinite; 396 } else if (sign == ZeroSign) { 397 if (realAkk > static_cast<RealScalar>(0)) sign = PositiveSemiDef; 398 else if (realAkk < static_cast<RealScalar>(0)) sign = NegativeSemiDef; 399 } 400 } 401 402 return ret; 403 } 404 405 // Reference for the algorithm: Davis and Hager, "Multiple Rank 406 // Modifications of a Sparse Cholesky Factorization" (Algorithm 1) 407 // Trivial rearrangements of their computations (Timothy E. Holy) 408 // allow their algorithm to work for rank-1 updates even if the 409 // original matrix is not of full rank. 410 // Here only rank-1 updates are implemented, to reduce the 411 // requirement for intermediate storage and improve accuracy 412 template<typename MatrixType, typename WDerived> 413 static bool updateInPlace(MatrixType& mat, MatrixBase<WDerived>& w, const typename MatrixType::RealScalar& sigma=1) 414 { 415 using numext::isfinite; 416 typedef typename MatrixType::Scalar Scalar; 417 typedef typename MatrixType::RealScalar RealScalar; 418 419 const Index size = mat.rows(); 420 eigen_assert(mat.cols() == size && w.size()==size); 421 422 RealScalar alpha = 1; 423 424 // Apply the update 425 for (Index j = 0; j < size; j++) 426 { 427 // Check for termination due to an original decomposition of low-rank 428 if (!(isfinite)(alpha)) 429 break; 430 431 // Update the diagonal terms 432 RealScalar dj = numext::real(mat.coeff(j,j)); 433 Scalar wj = w.coeff(j); 434 RealScalar swj2 = sigma*numext::abs2(wj); 435 RealScalar gamma = dj*alpha + swj2; 436 437 mat.coeffRef(j,j) += swj2/alpha; 438 alpha += swj2/dj; 439 440 441 // Update the terms of L 442 Index rs = size-j-1; 443 w.tail(rs) -= wj * mat.col(j).tail(rs); 444 if(gamma != 0) 445 mat.col(j).tail(rs) += (sigma*numext::conj(wj)/gamma)*w.tail(rs); 446 } 447 return true; 448 } 449 450 template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType> 451 static bool update(MatrixType& mat, const TranspositionType& transpositions, Workspace& tmp, const WType& w, const typename MatrixType::RealScalar& sigma=1) 452 { 453 // Apply the permutation to the input w 454 tmp = transpositions * w; 455 456 return ldlt_inplace<Lower>::updateInPlace(mat,tmp,sigma); 457 } 458 }; 459 460 template<> struct ldlt_inplace<Upper> 461 { 462 template<typename MatrixType, typename TranspositionType, typename Workspace> 463 static EIGEN_STRONG_INLINE bool unblocked(MatrixType& mat, TranspositionType& transpositions, Workspace& temp, SignMatrix& sign) 464 { 465 Transpose<MatrixType> matt(mat); 466 return ldlt_inplace<Lower>::unblocked(matt, transpositions, temp, sign); 467 } 468 469 template<typename MatrixType, typename TranspositionType, typename Workspace, typename WType> 470 static EIGEN_STRONG_INLINE bool update(MatrixType& mat, TranspositionType& transpositions, Workspace& tmp, WType& w, const typename MatrixType::RealScalar& sigma=1) 471 { 472 Transpose<MatrixType> matt(mat); 473 return ldlt_inplace<Lower>::update(matt, transpositions, tmp, w.conjugate(), sigma); 474 } 475 }; 476 477 template<typename MatrixType> struct LDLT_Traits<MatrixType,Lower> 478 { 479 typedef const TriangularView<const MatrixType, UnitLower> MatrixL; 480 typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitUpper> MatrixU; 481 static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); } 482 static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); } 483 }; 484 485 template<typename MatrixType> struct LDLT_Traits<MatrixType,Upper> 486 { 487 typedef const TriangularView<const typename MatrixType::AdjointReturnType, UnitLower> MatrixL; 488 typedef const TriangularView<const MatrixType, UnitUpper> MatrixU; 489 static inline MatrixL getL(const MatrixType& m) { return MatrixL(m.adjoint()); } 490 static inline MatrixU getU(const MatrixType& m) { return MatrixU(m); } 491 }; 492 493 } // end namespace internal 494 495 /** Compute / recompute the LDLT decomposition A = L D L^* = U^* D U of \a matrix 496 */ 497 template<typename MatrixType, int _UpLo> 498 template<typename InputType> 499 LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::compute(const EigenBase<InputType>& a) 500 { 501 check_template_parameters(); 502 503 eigen_assert(a.rows()==a.cols()); 504 const Index size = a.rows(); 505 506 m_matrix = a.derived(); 507 508 // Compute matrix L1 norm = max abs column sum. 509 m_l1_norm = RealScalar(0); 510 // TODO move this code to SelfAdjointView 511 for (Index col = 0; col < size; ++col) { 512 RealScalar abs_col_sum; 513 if (_UpLo == Lower) 514 abs_col_sum = m_matrix.col(col).tail(size - col).template lpNorm<1>() + m_matrix.row(col).head(col).template lpNorm<1>(); 515 else 516 abs_col_sum = m_matrix.col(col).head(col).template lpNorm<1>() + m_matrix.row(col).tail(size - col).template lpNorm<1>(); 517 if (abs_col_sum > m_l1_norm) 518 m_l1_norm = abs_col_sum; 519 } 520 521 m_transpositions.resize(size); 522 m_isInitialized = false; 523 m_temporary.resize(size); 524 m_sign = internal::ZeroSign; 525 526 m_info = internal::ldlt_inplace<UpLo>::unblocked(m_matrix, m_transpositions, m_temporary, m_sign) ? Success : NumericalIssue; 527 528 m_isInitialized = true; 529 return *this; 530 } 531 532 /** Update the LDLT decomposition: given A = L D L^T, efficiently compute the decomposition of A + sigma w w^T. 533 * \param w a vector to be incorporated into the decomposition. 534 * \param sigma a scalar, +1 for updates and -1 for "downdates," which correspond to removing previously-added column vectors. Optional; default value is +1. 535 * \sa setZero() 536 */ 537 template<typename MatrixType, int _UpLo> 538 template<typename Derived> 539 LDLT<MatrixType,_UpLo>& LDLT<MatrixType,_UpLo>::rankUpdate(const MatrixBase<Derived>& w, const typename LDLT<MatrixType,_UpLo>::RealScalar& sigma) 540 { 541 typedef typename TranspositionType::StorageIndex IndexType; 542 const Index size = w.rows(); 543 if (m_isInitialized) 544 { 545 eigen_assert(m_matrix.rows()==size); 546 } 547 else 548 { 549 m_matrix.resize(size,size); 550 m_matrix.setZero(); 551 m_transpositions.resize(size); 552 for (Index i = 0; i < size; i++) 553 m_transpositions.coeffRef(i) = IndexType(i); 554 m_temporary.resize(size); 555 m_sign = sigma>=0 ? internal::PositiveSemiDef : internal::NegativeSemiDef; 556 m_isInitialized = true; 557 } 558 559 internal::ldlt_inplace<UpLo>::update(m_matrix, m_transpositions, m_temporary, w, sigma); 560 561 return *this; 562 } 563 564 #ifndef EIGEN_PARSED_BY_DOXYGEN 565 template<typename _MatrixType, int _UpLo> 566 template<typename RhsType, typename DstType> 567 void LDLT<_MatrixType,_UpLo>::_solve_impl(const RhsType &rhs, DstType &dst) const 568 { 569 _solve_impl_transposed<true>(rhs, dst); 570 } 571 572 template<typename _MatrixType,int _UpLo> 573 template<bool Conjugate, typename RhsType, typename DstType> 574 void LDLT<_MatrixType,_UpLo>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const 575 { 576 // dst = P b 577 dst = m_transpositions * rhs; 578 579 // dst = L^-1 (P b) 580 // dst = L^-*T (P b) 581 matrixL().template conjugateIf<!Conjugate>().solveInPlace(dst); 582 583 // dst = D^-* (L^-1 P b) 584 // dst = D^-1 (L^-*T P b) 585 // more precisely, use pseudo-inverse of D (see bug 241) 586 using std::abs; 587 const typename Diagonal<const MatrixType>::RealReturnType vecD(vectorD()); 588 // In some previous versions, tolerance was set to the max of 1/highest (or rather numeric_limits::min()) 589 // and the maximal diagonal entry * epsilon as motivated by LAPACK's xGELSS: 590 // RealScalar tolerance = numext::maxi(vecD.array().abs().maxCoeff() * NumTraits<RealScalar>::epsilon(),RealScalar(1) / NumTraits<RealScalar>::highest()); 591 // However, LDLT is not rank revealing, and so adjusting the tolerance wrt to the highest 592 // diagonal element is not well justified and leads to numerical issues in some cases. 593 // Moreover, Lapack's xSYTRS routines use 0 for the tolerance. 594 // Using numeric_limits::min() gives us more robustness to denormals. 595 RealScalar tolerance = (std::numeric_limits<RealScalar>::min)(); 596 for (Index i = 0; i < vecD.size(); ++i) 597 { 598 if(abs(vecD(i)) > tolerance) 599 dst.row(i) /= vecD(i); 600 else 601 dst.row(i).setZero(); 602 } 603 604 // dst = L^-* (D^-* L^-1 P b) 605 // dst = L^-T (D^-1 L^-*T P b) 606 matrixL().transpose().template conjugateIf<Conjugate>().solveInPlace(dst); 607 608 // dst = P^T (L^-* D^-* L^-1 P b) = A^-1 b 609 // dst = P^-T (L^-T D^-1 L^-*T P b) = A^-1 b 610 dst = m_transpositions.transpose() * dst; 611 } 612 #endif 613 614 /** \internal use x = ldlt_object.solve(x); 615 * 616 * This is the \em in-place version of solve(). 617 * 618 * \param bAndX represents both the right-hand side matrix b and result x. 619 * 620 * \returns true always! If you need to check for existence of solutions, use another decomposition like LU, QR, or SVD. 621 * 622 * This version avoids a copy when the right hand side matrix b is not 623 * needed anymore. 624 * 625 * \sa LDLT::solve(), MatrixBase::ldlt() 626 */ 627 template<typename MatrixType,int _UpLo> 628 template<typename Derived> 629 bool LDLT<MatrixType,_UpLo>::solveInPlace(MatrixBase<Derived> &bAndX) const 630 { 631 eigen_assert(m_isInitialized && "LDLT is not initialized."); 632 eigen_assert(m_matrix.rows() == bAndX.rows()); 633 634 bAndX = this->solve(bAndX); 635 636 return true; 637 } 638 639 /** \returns the matrix represented by the decomposition, 640 * i.e., it returns the product: P^T L D L^* P. 641 * This function is provided for debug purpose. */ 642 template<typename MatrixType, int _UpLo> 643 MatrixType LDLT<MatrixType,_UpLo>::reconstructedMatrix() const 644 { 645 eigen_assert(m_isInitialized && "LDLT is not initialized."); 646 const Index size = m_matrix.rows(); 647 MatrixType res(size,size); 648 649 // P 650 res.setIdentity(); 651 res = transpositionsP() * res; 652 // L^* P 653 res = matrixU() * res; 654 // D(L^*P) 655 res = vectorD().real().asDiagonal() * res; 656 // L(DL^*P) 657 res = matrixL() * res; 658 // P^T (LDL^*P) 659 res = transpositionsP().transpose() * res; 660 661 return res; 662 } 663 664 /** \cholesky_module 665 * \returns the Cholesky decomposition with full pivoting without square root of \c *this 666 * \sa MatrixBase::ldlt() 667 */ 668 template<typename MatrixType, unsigned int UpLo> 669 inline const LDLT<typename SelfAdjointView<MatrixType, UpLo>::PlainObject, UpLo> 670 SelfAdjointView<MatrixType, UpLo>::ldlt() const 671 { 672 return LDLT<PlainObject,UpLo>(m_matrix); 673 } 674 675 /** \cholesky_module 676 * \returns the Cholesky decomposition with full pivoting without square root of \c *this 677 * \sa SelfAdjointView::ldlt() 678 */ 679 template<typename Derived> 680 inline const LDLT<typename MatrixBase<Derived>::PlainObject> 681 MatrixBase<Derived>::ldlt() const 682 { 683 return LDLT<PlainObject>(derived()); 684 } 685 686 } // end namespace Eigen 687 688 #endif // EIGEN_LDLT_H