SimplicialCholesky.h (24216B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2012 Gael Guennebaud <gael.guennebaud@inria.fr> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_SIMPLICIAL_CHOLESKY_H 11 #define EIGEN_SIMPLICIAL_CHOLESKY_H 12 13 namespace Eigen { 14 15 enum SimplicialCholeskyMode { 16 SimplicialCholeskyLLT, 17 SimplicialCholeskyLDLT 18 }; 19 20 namespace internal { 21 template<typename CholMatrixType, typename InputMatrixType> 22 struct simplicial_cholesky_grab_input { 23 typedef CholMatrixType const * ConstCholMatrixPtr; 24 static void run(const InputMatrixType& input, ConstCholMatrixPtr &pmat, CholMatrixType &tmp) 25 { 26 tmp = input; 27 pmat = &tmp; 28 } 29 }; 30 31 template<typename MatrixType> 32 struct simplicial_cholesky_grab_input<MatrixType,MatrixType> { 33 typedef MatrixType const * ConstMatrixPtr; 34 static void run(const MatrixType& input, ConstMatrixPtr &pmat, MatrixType &/*tmp*/) 35 { 36 pmat = &input; 37 } 38 }; 39 } // end namespace internal 40 41 /** \ingroup SparseCholesky_Module 42 * \brief A base class for direct sparse Cholesky factorizations 43 * 44 * This is a base class for LL^T and LDL^T Cholesky factorizations of sparse matrices that are 45 * selfadjoint and positive definite. These factorizations allow for solving A.X = B where 46 * X and B can be either dense or sparse. 47 * 48 * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization 49 * such that the factorized matrix is P A P^-1. 50 * 51 * \tparam Derived the type of the derived class, that is the actual factorization type. 52 * 53 */ 54 template<typename Derived> 55 class SimplicialCholeskyBase : public SparseSolverBase<Derived> 56 { 57 typedef SparseSolverBase<Derived> Base; 58 using Base::m_isInitialized; 59 60 public: 61 typedef typename internal::traits<Derived>::MatrixType MatrixType; 62 typedef typename internal::traits<Derived>::OrderingType OrderingType; 63 enum { UpLo = internal::traits<Derived>::UpLo }; 64 typedef typename MatrixType::Scalar Scalar; 65 typedef typename MatrixType::RealScalar RealScalar; 66 typedef typename MatrixType::StorageIndex StorageIndex; 67 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType; 68 typedef CholMatrixType const * ConstCholMatrixPtr; 69 typedef Matrix<Scalar,Dynamic,1> VectorType; 70 typedef Matrix<StorageIndex,Dynamic,1> VectorI; 71 72 enum { 73 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 74 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 75 }; 76 77 public: 78 79 using Base::derived; 80 81 /** Default constructor */ 82 SimplicialCholeskyBase() 83 : m_info(Success), 84 m_factorizationIsOk(false), 85 m_analysisIsOk(false), 86 m_shiftOffset(0), 87 m_shiftScale(1) 88 {} 89 90 explicit SimplicialCholeskyBase(const MatrixType& matrix) 91 : m_info(Success), 92 m_factorizationIsOk(false), 93 m_analysisIsOk(false), 94 m_shiftOffset(0), 95 m_shiftScale(1) 96 { 97 derived().compute(matrix); 98 } 99 100 ~SimplicialCholeskyBase() 101 { 102 } 103 104 Derived& derived() { return *static_cast<Derived*>(this); } 105 const Derived& derived() const { return *static_cast<const Derived*>(this); } 106 107 inline Index cols() const { return m_matrix.cols(); } 108 inline Index rows() const { return m_matrix.rows(); } 109 110 /** \brief Reports whether previous computation was successful. 111 * 112 * \returns \c Success if computation was successful, 113 * \c NumericalIssue if the matrix.appears to be negative. 114 */ 115 ComputationInfo info() const 116 { 117 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 118 return m_info; 119 } 120 121 /** \returns the permutation P 122 * \sa permutationPinv() */ 123 const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& permutationP() const 124 { return m_P; } 125 126 /** \returns the inverse P^-1 of the permutation P 127 * \sa permutationP() */ 128 const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& permutationPinv() const 129 { return m_Pinv; } 130 131 /** Sets the shift parameters that will be used to adjust the diagonal coefficients during the numerical factorization. 132 * 133 * During the numerical factorization, the diagonal coefficients are transformed by the following linear model:\n 134 * \c d_ii = \a offset + \a scale * \c d_ii 135 * 136 * The default is the identity transformation with \a offset=0, and \a scale=1. 137 * 138 * \returns a reference to \c *this. 139 */ 140 Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1) 141 { 142 m_shiftOffset = offset; 143 m_shiftScale = scale; 144 return derived(); 145 } 146 147 #ifndef EIGEN_PARSED_BY_DOXYGEN 148 /** \internal */ 149 template<typename Stream> 150 void dumpMemory(Stream& s) 151 { 152 int total = 0; 153 s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n"; 154 s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n"; 155 s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 156 s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 157 s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 158 s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 159 s << " TOTAL: " << (total>> 20) << "Mb" << "\n"; 160 } 161 162 /** \internal */ 163 template<typename Rhs,typename Dest> 164 void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const 165 { 166 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 167 eigen_assert(m_matrix.rows()==b.rows()); 168 169 if(m_info!=Success) 170 return; 171 172 if(m_P.size()>0) 173 dest = m_P * b; 174 else 175 dest = b; 176 177 if(m_matrix.nonZeros()>0) // otherwise L==I 178 derived().matrixL().solveInPlace(dest); 179 180 if(m_diag.size()>0) 181 dest = m_diag.asDiagonal().inverse() * dest; 182 183 if (m_matrix.nonZeros()>0) // otherwise U==I 184 derived().matrixU().solveInPlace(dest); 185 186 if(m_P.size()>0) 187 dest = m_Pinv * dest; 188 } 189 190 template<typename Rhs,typename Dest> 191 void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const 192 { 193 internal::solve_sparse_through_dense_panels(derived(), b, dest); 194 } 195 196 #endif // EIGEN_PARSED_BY_DOXYGEN 197 198 protected: 199 200 /** Computes the sparse Cholesky decomposition of \a matrix */ 201 template<bool DoLDLT> 202 void compute(const MatrixType& matrix) 203 { 204 eigen_assert(matrix.rows()==matrix.cols()); 205 Index size = matrix.cols(); 206 CholMatrixType tmp(size,size); 207 ConstCholMatrixPtr pmat; 208 ordering(matrix, pmat, tmp); 209 analyzePattern_preordered(*pmat, DoLDLT); 210 factorize_preordered<DoLDLT>(*pmat); 211 } 212 213 template<bool DoLDLT> 214 void factorize(const MatrixType& a) 215 { 216 eigen_assert(a.rows()==a.cols()); 217 Index size = a.cols(); 218 CholMatrixType tmp(size,size); 219 ConstCholMatrixPtr pmat; 220 221 if(m_P.size() == 0 && (int(UpLo) & int(Upper)) == Upper) 222 { 223 // If there is no ordering, try to directly use the input matrix without any copy 224 internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp); 225 } 226 else 227 { 228 tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); 229 pmat = &tmp; 230 } 231 232 factorize_preordered<DoLDLT>(*pmat); 233 } 234 235 template<bool DoLDLT> 236 void factorize_preordered(const CholMatrixType& a); 237 238 void analyzePattern(const MatrixType& a, bool doLDLT) 239 { 240 eigen_assert(a.rows()==a.cols()); 241 Index size = a.cols(); 242 CholMatrixType tmp(size,size); 243 ConstCholMatrixPtr pmat; 244 ordering(a, pmat, tmp); 245 analyzePattern_preordered(*pmat,doLDLT); 246 } 247 void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT); 248 249 void ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap); 250 251 /** keeps off-diagonal entries; drops diagonal entries */ 252 struct keep_diag { 253 inline bool operator() (const Index& row, const Index& col, const Scalar&) const 254 { 255 return row!=col; 256 } 257 }; 258 259 mutable ComputationInfo m_info; 260 bool m_factorizationIsOk; 261 bool m_analysisIsOk; 262 263 CholMatrixType m_matrix; 264 VectorType m_diag; // the diagonal coefficients (LDLT mode) 265 VectorI m_parent; // elimination tree 266 VectorI m_nonZerosPerCol; 267 PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_P; // the permutation 268 PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_Pinv; // the inverse permutation 269 270 RealScalar m_shiftOffset; 271 RealScalar m_shiftScale; 272 }; 273 274 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLLT; 275 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLDLT; 276 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialCholesky; 277 278 namespace internal { 279 280 template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialLLT<_MatrixType,_UpLo,_Ordering> > 281 { 282 typedef _MatrixType MatrixType; 283 typedef _Ordering OrderingType; 284 enum { UpLo = _UpLo }; 285 typedef typename MatrixType::Scalar Scalar; 286 typedef typename MatrixType::StorageIndex StorageIndex; 287 typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType; 288 typedef TriangularView<const CholMatrixType, Eigen::Lower> MatrixL; 289 typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU; 290 static inline MatrixL getL(const CholMatrixType& m) { return MatrixL(m); } 291 static inline MatrixU getU(const CholMatrixType& m) { return MatrixU(m.adjoint()); } 292 }; 293 294 template<typename _MatrixType,int _UpLo, typename _Ordering> struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> > 295 { 296 typedef _MatrixType MatrixType; 297 typedef _Ordering OrderingType; 298 enum { UpLo = _UpLo }; 299 typedef typename MatrixType::Scalar Scalar; 300 typedef typename MatrixType::StorageIndex StorageIndex; 301 typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType; 302 typedef TriangularView<const CholMatrixType, Eigen::UnitLower> MatrixL; 303 typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU; 304 static inline MatrixL getL(const CholMatrixType& m) { return MatrixL(m); } 305 static inline MatrixU getU(const CholMatrixType& m) { return MatrixU(m.adjoint()); } 306 }; 307 308 template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> > 309 { 310 typedef _MatrixType MatrixType; 311 typedef _Ordering OrderingType; 312 enum { UpLo = _UpLo }; 313 }; 314 315 } 316 317 /** \ingroup SparseCholesky_Module 318 * \class SimplicialLLT 319 * \brief A direct sparse LLT Cholesky factorizations 320 * 321 * This class provides a LL^T Cholesky factorizations of sparse matrices that are 322 * selfadjoint and positive definite. The factorization allows for solving A.X = B where 323 * X and B can be either dense or sparse. 324 * 325 * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization 326 * such that the factorized matrix is P A P^-1. 327 * 328 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 329 * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower 330 * or Upper. Default is Lower. 331 * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<> 332 * 333 * \implsparsesolverconcept 334 * 335 * \sa class SimplicialLDLT, class AMDOrdering, class NaturalOrdering 336 */ 337 template<typename _MatrixType, int _UpLo, typename _Ordering> 338 class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> > 339 { 340 public: 341 typedef _MatrixType MatrixType; 342 enum { UpLo = _UpLo }; 343 typedef SimplicialCholeskyBase<SimplicialLLT> Base; 344 typedef typename MatrixType::Scalar Scalar; 345 typedef typename MatrixType::RealScalar RealScalar; 346 typedef typename MatrixType::StorageIndex StorageIndex; 347 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; 348 typedef Matrix<Scalar,Dynamic,1> VectorType; 349 typedef internal::traits<SimplicialLLT> Traits; 350 typedef typename Traits::MatrixL MatrixL; 351 typedef typename Traits::MatrixU MatrixU; 352 public: 353 /** Default constructor */ 354 SimplicialLLT() : Base() {} 355 /** Constructs and performs the LLT factorization of \a matrix */ 356 explicit SimplicialLLT(const MatrixType& matrix) 357 : Base(matrix) {} 358 359 /** \returns an expression of the factor L */ 360 inline const MatrixL matrixL() const { 361 eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized"); 362 return Traits::getL(Base::m_matrix); 363 } 364 365 /** \returns an expression of the factor U (= L^*) */ 366 inline const MatrixU matrixU() const { 367 eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized"); 368 return Traits::getU(Base::m_matrix); 369 } 370 371 /** Computes the sparse Cholesky decomposition of \a matrix */ 372 SimplicialLLT& compute(const MatrixType& matrix) 373 { 374 Base::template compute<false>(matrix); 375 return *this; 376 } 377 378 /** Performs a symbolic decomposition on the sparcity of \a matrix. 379 * 380 * This function is particularly useful when solving for several problems having the same structure. 381 * 382 * \sa factorize() 383 */ 384 void analyzePattern(const MatrixType& a) 385 { 386 Base::analyzePattern(a, false); 387 } 388 389 /** Performs a numeric decomposition of \a matrix 390 * 391 * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. 392 * 393 * \sa analyzePattern() 394 */ 395 void factorize(const MatrixType& a) 396 { 397 Base::template factorize<false>(a); 398 } 399 400 /** \returns the determinant of the underlying matrix from the current factorization */ 401 Scalar determinant() const 402 { 403 Scalar detL = Base::m_matrix.diagonal().prod(); 404 return numext::abs2(detL); 405 } 406 }; 407 408 /** \ingroup SparseCholesky_Module 409 * \class SimplicialLDLT 410 * \brief A direct sparse LDLT Cholesky factorizations without square root. 411 * 412 * This class provides a LDL^T Cholesky factorizations without square root of sparse matrices that are 413 * selfadjoint and positive definite. The factorization allows for solving A.X = B where 414 * X and B can be either dense or sparse. 415 * 416 * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization 417 * such that the factorized matrix is P A P^-1. 418 * 419 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 420 * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower 421 * or Upper. Default is Lower. 422 * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<> 423 * 424 * \implsparsesolverconcept 425 * 426 * \sa class SimplicialLLT, class AMDOrdering, class NaturalOrdering 427 */ 428 template<typename _MatrixType, int _UpLo, typename _Ordering> 429 class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> > 430 { 431 public: 432 typedef _MatrixType MatrixType; 433 enum { UpLo = _UpLo }; 434 typedef SimplicialCholeskyBase<SimplicialLDLT> Base; 435 typedef typename MatrixType::Scalar Scalar; 436 typedef typename MatrixType::RealScalar RealScalar; 437 typedef typename MatrixType::StorageIndex StorageIndex; 438 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType; 439 typedef Matrix<Scalar,Dynamic,1> VectorType; 440 typedef internal::traits<SimplicialLDLT> Traits; 441 typedef typename Traits::MatrixL MatrixL; 442 typedef typename Traits::MatrixU MatrixU; 443 public: 444 /** Default constructor */ 445 SimplicialLDLT() : Base() {} 446 447 /** Constructs and performs the LLT factorization of \a matrix */ 448 explicit SimplicialLDLT(const MatrixType& matrix) 449 : Base(matrix) {} 450 451 /** \returns a vector expression of the diagonal D */ 452 inline const VectorType vectorD() const { 453 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); 454 return Base::m_diag; 455 } 456 /** \returns an expression of the factor L */ 457 inline const MatrixL matrixL() const { 458 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); 459 return Traits::getL(Base::m_matrix); 460 } 461 462 /** \returns an expression of the factor U (= L^*) */ 463 inline const MatrixU matrixU() const { 464 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); 465 return Traits::getU(Base::m_matrix); 466 } 467 468 /** Computes the sparse Cholesky decomposition of \a matrix */ 469 SimplicialLDLT& compute(const MatrixType& matrix) 470 { 471 Base::template compute<true>(matrix); 472 return *this; 473 } 474 475 /** Performs a symbolic decomposition on the sparcity of \a matrix. 476 * 477 * This function is particularly useful when solving for several problems having the same structure. 478 * 479 * \sa factorize() 480 */ 481 void analyzePattern(const MatrixType& a) 482 { 483 Base::analyzePattern(a, true); 484 } 485 486 /** Performs a numeric decomposition of \a matrix 487 * 488 * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. 489 * 490 * \sa analyzePattern() 491 */ 492 void factorize(const MatrixType& a) 493 { 494 Base::template factorize<true>(a); 495 } 496 497 /** \returns the determinant of the underlying matrix from the current factorization */ 498 Scalar determinant() const 499 { 500 return Base::m_diag.prod(); 501 } 502 }; 503 504 /** \deprecated use SimplicialLDLT or class SimplicialLLT 505 * \ingroup SparseCholesky_Module 506 * \class SimplicialCholesky 507 * 508 * \sa class SimplicialLDLT, class SimplicialLLT 509 */ 510 template<typename _MatrixType, int _UpLo, typename _Ordering> 511 class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> > 512 { 513 public: 514 typedef _MatrixType MatrixType; 515 enum { UpLo = _UpLo }; 516 typedef SimplicialCholeskyBase<SimplicialCholesky> Base; 517 typedef typename MatrixType::Scalar Scalar; 518 typedef typename MatrixType::RealScalar RealScalar; 519 typedef typename MatrixType::StorageIndex StorageIndex; 520 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType; 521 typedef Matrix<Scalar,Dynamic,1> VectorType; 522 typedef internal::traits<SimplicialCholesky> Traits; 523 typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits; 524 typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits; 525 public: 526 SimplicialCholesky() : Base(), m_LDLT(true) {} 527 528 explicit SimplicialCholesky(const MatrixType& matrix) 529 : Base(), m_LDLT(true) 530 { 531 compute(matrix); 532 } 533 534 SimplicialCholesky& setMode(SimplicialCholeskyMode mode) 535 { 536 switch(mode) 537 { 538 case SimplicialCholeskyLLT: 539 m_LDLT = false; 540 break; 541 case SimplicialCholeskyLDLT: 542 m_LDLT = true; 543 break; 544 default: 545 break; 546 } 547 548 return *this; 549 } 550 551 inline const VectorType vectorD() const { 552 eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized"); 553 return Base::m_diag; 554 } 555 inline const CholMatrixType rawMatrix() const { 556 eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized"); 557 return Base::m_matrix; 558 } 559 560 /** Computes the sparse Cholesky decomposition of \a matrix */ 561 SimplicialCholesky& compute(const MatrixType& matrix) 562 { 563 if(m_LDLT) 564 Base::template compute<true>(matrix); 565 else 566 Base::template compute<false>(matrix); 567 return *this; 568 } 569 570 /** Performs a symbolic decomposition on the sparcity of \a matrix. 571 * 572 * This function is particularly useful when solving for several problems having the same structure. 573 * 574 * \sa factorize() 575 */ 576 void analyzePattern(const MatrixType& a) 577 { 578 Base::analyzePattern(a, m_LDLT); 579 } 580 581 /** Performs a numeric decomposition of \a matrix 582 * 583 * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. 584 * 585 * \sa analyzePattern() 586 */ 587 void factorize(const MatrixType& a) 588 { 589 if(m_LDLT) 590 Base::template factorize<true>(a); 591 else 592 Base::template factorize<false>(a); 593 } 594 595 /** \internal */ 596 template<typename Rhs,typename Dest> 597 void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const 598 { 599 eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 600 eigen_assert(Base::m_matrix.rows()==b.rows()); 601 602 if(Base::m_info!=Success) 603 return; 604 605 if(Base::m_P.size()>0) 606 dest = Base::m_P * b; 607 else 608 dest = b; 609 610 if(Base::m_matrix.nonZeros()>0) // otherwise L==I 611 { 612 if(m_LDLT) 613 LDLTTraits::getL(Base::m_matrix).solveInPlace(dest); 614 else 615 LLTTraits::getL(Base::m_matrix).solveInPlace(dest); 616 } 617 618 if(Base::m_diag.size()>0) 619 dest = Base::m_diag.real().asDiagonal().inverse() * dest; 620 621 if (Base::m_matrix.nonZeros()>0) // otherwise I==I 622 { 623 if(m_LDLT) 624 LDLTTraits::getU(Base::m_matrix).solveInPlace(dest); 625 else 626 LLTTraits::getU(Base::m_matrix).solveInPlace(dest); 627 } 628 629 if(Base::m_P.size()>0) 630 dest = Base::m_Pinv * dest; 631 } 632 633 /** \internal */ 634 template<typename Rhs,typename Dest> 635 void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const 636 { 637 internal::solve_sparse_through_dense_panels(*this, b, dest); 638 } 639 640 Scalar determinant() const 641 { 642 if(m_LDLT) 643 { 644 return Base::m_diag.prod(); 645 } 646 else 647 { 648 Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod(); 649 return numext::abs2(detL); 650 } 651 } 652 653 protected: 654 bool m_LDLT; 655 }; 656 657 template<typename Derived> 658 void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap) 659 { 660 eigen_assert(a.rows()==a.cols()); 661 const Index size = a.rows(); 662 pmat = ≈ 663 // Note that ordering methods compute the inverse permutation 664 if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value) 665 { 666 { 667 CholMatrixType C; 668 C = a.template selfadjointView<UpLo>(); 669 670 OrderingType ordering; 671 ordering(C,m_Pinv); 672 } 673 674 if(m_Pinv.size()>0) m_P = m_Pinv.inverse(); 675 else m_P.resize(0); 676 677 ap.resize(size,size); 678 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); 679 } 680 else 681 { 682 m_Pinv.resize(0); 683 m_P.resize(0); 684 if(int(UpLo)==int(Lower) || MatrixType::IsRowMajor) 685 { 686 // we have to transpose the lower part to to the upper one 687 ap.resize(size,size); 688 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>(); 689 } 690 else 691 internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap); 692 } 693 } 694 695 } // end namespace Eigen 696 697 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H