1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud (at) 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 /* 11 12 NOTE: the _symbolic, and _numeric functions has been adapted from 13 the LDL library: 14 15 LDL Copyright (c) 2005 by Timothy A. Davis. All Rights Reserved. 16 17 LDL License: 18 19 Your use or distribution of LDL or any modified version of 20 LDL implies that you agree to this License. 21 22 This library is free software; you can redistribute it and/or 23 modify it under the terms of the GNU Lesser General Public 24 License as published by the Free Software Foundation; either 25 version 2.1 of the License, or (at your option) any later version. 26 27 This library is distributed in the hope that it will be useful, 28 but WITHOUT ANY WARRANTY; without even the implied warranty of 29 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 30 Lesser General Public License for more details. 31 32 You should have received a copy of the GNU Lesser General Public 33 License along with this library; if not, write to the Free Software 34 Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 35 USA 36 37 Permission is hereby granted to use or copy this program under the 38 terms of the GNU LGPL, provided that the Copyright, this License, 39 and the Availability of the original version is retained on all copies. 40 User documentation of any code that uses this code or any modified 41 version of this code must cite the Copyright, this License, the 42 Availability note, and "Used by permission." Permission to modify 43 the code and to distribute modified code is granted, provided the 44 Copyright, this License, and the Availability note are retained, 45 and a notice that the code was modified is included. 46 */ 47 48 #include "../Core/util/NonMPL2.h" 49 50 #ifndef EIGEN_SIMPLICIAL_CHOLESKY_H 51 #define EIGEN_SIMPLICIAL_CHOLESKY_H 52 53 namespace Eigen { 54 55 enum SimplicialCholeskyMode { 56 SimplicialCholeskyLLT, 57 SimplicialCholeskyLDLT 58 }; 59 60 /** \ingroup SparseCholesky_Module 61 * \brief A direct sparse Cholesky factorizations 62 * 63 * These classes provide LL^T and LDL^T Cholesky factorizations of sparse matrices that are 64 * selfadjoint and positive definite. The factorization allows for solving A.X = B where 65 * X and B can be either dense or sparse. 66 * 67 * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization 68 * such that the factorized matrix is P A P^-1. 69 * 70 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 71 * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower 72 * or Upper. Default is Lower. 73 * 74 */ 75 template<typename Derived> 76 class SimplicialCholeskyBase : internal::noncopyable 77 { 78 public: 79 typedef typename internal::traits<Derived>::MatrixType MatrixType; 80 enum { UpLo = internal::traits<Derived>::UpLo }; 81 typedef typename MatrixType::Scalar Scalar; 82 typedef typename MatrixType::RealScalar RealScalar; 83 typedef typename MatrixType::Index Index; 84 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; 85 typedef Matrix<Scalar,Dynamic,1> VectorType; 86 87 public: 88 89 /** Default constructor */ 90 SimplicialCholeskyBase() 91 : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1) 92 {} 93 94 SimplicialCholeskyBase(const MatrixType& matrix) 95 : m_info(Success), m_isInitialized(false), m_shiftOffset(0), 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 succesful, 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 solution x of \f$ A x = b \f$ using the current decomposition of A. 122 * 123 * \sa compute() 124 */ 125 template<typename Rhs> 126 inline const internal::solve_retval<SimplicialCholeskyBase, Rhs> 127 solve(const MatrixBase<Rhs>& b) const 128 { 129 eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized."); 130 eigen_assert(rows()==b.rows() 131 && "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b"); 132 return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived()); 133 } 134 135 /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A. 136 * 137 * \sa compute() 138 */ 139 template<typename Rhs> 140 inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs> 141 solve(const SparseMatrixBase<Rhs>& b) const 142 { 143 eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized."); 144 eigen_assert(rows()==b.rows() 145 && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b"); 146 return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived()); 147 } 148 149 /** \returns the permutation P 150 * \sa permutationPinv() */ 151 const PermutationMatrix<Dynamic,Dynamic,Index>& permutationP() const 152 { return m_P; } 153 154 /** \returns the inverse P^-1 of the permutation P 155 * \sa permutationP() */ 156 const PermutationMatrix<Dynamic,Dynamic,Index>& permutationPinv() const 157 { return m_Pinv; } 158 159 /** Sets the shift parameters that will be used to adjust the diagonal coefficients during the numerical factorization. 160 * 161 * During the numerical factorization, the diagonal coefficients are transformed by the following linear model:\n 162 * \c d_ii = \a offset + \a scale * \c d_ii 163 * 164 * The default is the identity transformation with \a offset=0, and \a scale=1. 165 * 166 * \returns a reference to \c *this. 167 */ 168 Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1) 169 { 170 m_shiftOffset = offset; 171 m_shiftScale = scale; 172 return derived(); 173 } 174 175 #ifndef EIGEN_PARSED_BY_DOXYGEN 176 /** \internal */ 177 template<typename Stream> 178 void dumpMemory(Stream& s) 179 { 180 int total = 0; 181 s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n"; 182 s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n"; 183 s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 184 s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 185 s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 186 s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 187 s << " TOTAL: " << (total>> 20) << "Mb" << "\n"; 188 } 189 190 /** \internal */ 191 template<typename Rhs,typename Dest> 192 void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const 193 { 194 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 195 eigen_assert(m_matrix.rows()==b.rows()); 196 197 if(m_info!=Success) 198 return; 199 200 if(m_P.size()>0) 201 dest = m_P * b; 202 else 203 dest = b; 204 205 if(m_matrix.nonZeros()>0) // otherwise L==I 206 derived().matrixL().solveInPlace(dest); 207 208 if(m_diag.size()>0) 209 dest = m_diag.asDiagonal().inverse() * dest; 210 211 if (m_matrix.nonZeros()>0) // otherwise U==I 212 derived().matrixU().solveInPlace(dest); 213 214 if(m_P.size()>0) 215 dest = m_Pinv * dest; 216 } 217 218 /** \internal */ 219 template<typename Rhs, typename DestScalar, int DestOptions, typename DestIndex> 220 void _solve_sparse(const Rhs& b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const 221 { 222 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 223 eigen_assert(m_matrix.rows()==b.rows()); 224 225 // we process the sparse rhs per block of NbColsAtOnce columns temporarily stored into a dense matrix. 226 static const int NbColsAtOnce = 4; 227 int rhsCols = b.cols(); 228 int size = b.rows(); 229 Eigen::Matrix<DestScalar,Dynamic,Dynamic> tmp(size,rhsCols); 230 for(int k=0; k<rhsCols; k+=NbColsAtOnce) 231 { 232 int actualCols = std::min<int>(rhsCols-k, NbColsAtOnce); 233 tmp.leftCols(actualCols) = b.middleCols(k,actualCols); 234 tmp.leftCols(actualCols) = derived().solve(tmp.leftCols(actualCols)); 235 dest.middleCols(k,actualCols) = tmp.leftCols(actualCols).sparseView(); 236 } 237 } 238 239 #endif // EIGEN_PARSED_BY_DOXYGEN 240 241 protected: 242 243 /** Computes the sparse Cholesky decomposition of \a matrix */ 244 template<bool DoLDLT> 245 void compute(const MatrixType& matrix) 246 { 247 eigen_assert(matrix.rows()==matrix.cols()); 248 Index size = matrix.cols(); 249 CholMatrixType ap(size,size); 250 ordering(matrix, ap); 251 analyzePattern_preordered(ap, DoLDLT); 252 factorize_preordered<DoLDLT>(ap); 253 } 254 255 template<bool DoLDLT> 256 void factorize(const MatrixType& a) 257 { 258 eigen_assert(a.rows()==a.cols()); 259 int size = a.cols(); 260 CholMatrixType ap(size,size); 261 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); 262 factorize_preordered<DoLDLT>(ap); 263 } 264 265 template<bool DoLDLT> 266 void factorize_preordered(const CholMatrixType& a); 267 268 void analyzePattern(const MatrixType& a, bool doLDLT) 269 { 270 eigen_assert(a.rows()==a.cols()); 271 int size = a.cols(); 272 CholMatrixType ap(size,size); 273 ordering(a, ap); 274 analyzePattern_preordered(ap,doLDLT); 275 } 276 void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT); 277 278 void ordering(const MatrixType& a, CholMatrixType& ap); 279 280 /** keeps off-diagonal entries; drops diagonal entries */ 281 struct keep_diag { 282 inline bool operator() (const Index& row, const Index& col, const Scalar&) const 283 { 284 return row!=col; 285 } 286 }; 287 288 mutable ComputationInfo m_info; 289 bool m_isInitialized; 290 bool m_factorizationIsOk; 291 bool m_analysisIsOk; 292 293 CholMatrixType m_matrix; 294 VectorType m_diag; // the diagonal coefficients (LDLT mode) 295 VectorXi m_parent; // elimination tree 296 VectorXi m_nonZerosPerCol; 297 PermutationMatrix<Dynamic,Dynamic,Index> m_P; // the permutation 298 PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // the inverse permutation 299 300 RealScalar m_shiftOffset; 301 RealScalar m_shiftScale; 302 }; 303 304 template<typename _MatrixType, int _UpLo = Lower> class SimplicialLLT; 305 template<typename _MatrixType, int _UpLo = Lower> class SimplicialLDLT; 306 template<typename _MatrixType, int _UpLo = Lower> class SimplicialCholesky; 307 308 namespace internal { 309 310 template<typename _MatrixType, int _UpLo> struct traits<SimplicialLLT<_MatrixType,_UpLo> > 311 { 312 typedef _MatrixType MatrixType; 313 enum { UpLo = _UpLo }; 314 typedef typename MatrixType::Scalar Scalar; 315 typedef typename MatrixType::Index Index; 316 typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType; 317 typedef SparseTriangularView<CholMatrixType, Eigen::Lower> MatrixL; 318 typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU; 319 static inline MatrixL getL(const MatrixType& m) { return m; } 320 static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } 321 }; 322 323 template<typename _MatrixType,int _UpLo> struct traits<SimplicialLDLT<_MatrixType,_UpLo> > 324 { 325 typedef _MatrixType MatrixType; 326 enum { UpLo = _UpLo }; 327 typedef typename MatrixType::Scalar Scalar; 328 typedef typename MatrixType::Index Index; 329 typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType; 330 typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower> MatrixL; 331 typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU; 332 static inline MatrixL getL(const MatrixType& m) { return m; } 333 static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } 334 }; 335 336 template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_MatrixType,_UpLo> > 337 { 338 typedef _MatrixType MatrixType; 339 enum { UpLo = _UpLo }; 340 }; 341 342 } 343 344 /** \ingroup SparseCholesky_Module 345 * \class SimplicialLLT 346 * \brief A direct sparse LLT Cholesky factorizations 347 * 348 * This class provides a LL^T Cholesky factorizations of sparse matrices that are 349 * selfadjoint and positive definite. The factorization allows for solving A.X = B where 350 * X and B can be either dense or sparse. 351 * 352 * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization 353 * such that the factorized matrix is P A P^-1. 354 * 355 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 356 * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower 357 * or Upper. Default is Lower. 358 * 359 * \sa class SimplicialLDLT 360 */ 361 template<typename _MatrixType, int _UpLo> 362 class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo> > 363 { 364 public: 365 typedef _MatrixType MatrixType; 366 enum { UpLo = _UpLo }; 367 typedef SimplicialCholeskyBase<SimplicialLLT> Base; 368 typedef typename MatrixType::Scalar Scalar; 369 typedef typename MatrixType::RealScalar RealScalar; 370 typedef typename MatrixType::Index Index; 371 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; 372 typedef Matrix<Scalar,Dynamic,1> VectorType; 373 typedef internal::traits<SimplicialLLT> Traits; 374 typedef typename Traits::MatrixL MatrixL; 375 typedef typename Traits::MatrixU MatrixU; 376 public: 377 /** Default constructor */ 378 SimplicialLLT() : Base() {} 379 /** Constructs and performs the LLT factorization of \a matrix */ 380 SimplicialLLT(const MatrixType& matrix) 381 : Base(matrix) {} 382 383 /** \returns an expression of the factor L */ 384 inline const MatrixL matrixL() const { 385 eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized"); 386 return Traits::getL(Base::m_matrix); 387 } 388 389 /** \returns an expression of the factor U (= L^*) */ 390 inline const MatrixU matrixU() const { 391 eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized"); 392 return Traits::getU(Base::m_matrix); 393 } 394 395 /** Computes the sparse Cholesky decomposition of \a matrix */ 396 SimplicialLLT& compute(const MatrixType& matrix) 397 { 398 Base::template compute<false>(matrix); 399 return *this; 400 } 401 402 /** Performs a symbolic decomposition on the sparcity of \a matrix. 403 * 404 * This function is particularly useful when solving for several problems having the same structure. 405 * 406 * \sa factorize() 407 */ 408 void analyzePattern(const MatrixType& a) 409 { 410 Base::analyzePattern(a, false); 411 } 412 413 /** Performs a numeric decomposition of \a matrix 414 * 415 * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. 416 * 417 * \sa analyzePattern() 418 */ 419 void factorize(const MatrixType& a) 420 { 421 Base::template factorize<false>(a); 422 } 423 424 /** \returns the determinant of the underlying matrix from the current factorization */ 425 Scalar determinant() const 426 { 427 Scalar detL = Base::m_matrix.diagonal().prod(); 428 return internal::abs2(detL); 429 } 430 }; 431 432 /** \ingroup SparseCholesky_Module 433 * \class SimplicialLDLT 434 * \brief A direct sparse LDLT Cholesky factorizations without square root. 435 * 436 * This class provides a LDL^T Cholesky factorizations without square root of sparse matrices that are 437 * selfadjoint and positive definite. The factorization allows for solving A.X = B where 438 * X and B can be either dense or sparse. 439 * 440 * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization 441 * such that the factorized matrix is P A P^-1. 442 * 443 * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<> 444 * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower 445 * or Upper. Default is Lower. 446 * 447 * \sa class SimplicialLLT 448 */ 449 template<typename _MatrixType, int _UpLo> 450 class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo> > 451 { 452 public: 453 typedef _MatrixType MatrixType; 454 enum { UpLo = _UpLo }; 455 typedef SimplicialCholeskyBase<SimplicialLDLT> Base; 456 typedef typename MatrixType::Scalar Scalar; 457 typedef typename MatrixType::RealScalar RealScalar; 458 typedef typename MatrixType::Index Index; 459 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; 460 typedef Matrix<Scalar,Dynamic,1> VectorType; 461 typedef internal::traits<SimplicialLDLT> Traits; 462 typedef typename Traits::MatrixL MatrixL; 463 typedef typename Traits::MatrixU MatrixU; 464 public: 465 /** Default constructor */ 466 SimplicialLDLT() : Base() {} 467 468 /** Constructs and performs the LLT factorization of \a matrix */ 469 SimplicialLDLT(const MatrixType& matrix) 470 : Base(matrix) {} 471 472 /** \returns a vector expression of the diagonal D */ 473 inline const VectorType vectorD() const { 474 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); 475 return Base::m_diag; 476 } 477 /** \returns an expression of the factor L */ 478 inline const MatrixL matrixL() const { 479 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); 480 return Traits::getL(Base::m_matrix); 481 } 482 483 /** \returns an expression of the factor U (= L^*) */ 484 inline const MatrixU matrixU() const { 485 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); 486 return Traits::getU(Base::m_matrix); 487 } 488 489 /** Computes the sparse Cholesky decomposition of \a matrix */ 490 SimplicialLDLT& compute(const MatrixType& matrix) 491 { 492 Base::template compute<true>(matrix); 493 return *this; 494 } 495 496 /** Performs a symbolic decomposition on the sparcity of \a matrix. 497 * 498 * This function is particularly useful when solving for several problems having the same structure. 499 * 500 * \sa factorize() 501 */ 502 void analyzePattern(const MatrixType& a) 503 { 504 Base::analyzePattern(a, true); 505 } 506 507 /** Performs a numeric decomposition of \a matrix 508 * 509 * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. 510 * 511 * \sa analyzePattern() 512 */ 513 void factorize(const MatrixType& a) 514 { 515 Base::template factorize<true>(a); 516 } 517 518 /** \returns the determinant of the underlying matrix from the current factorization */ 519 Scalar determinant() const 520 { 521 return Base::m_diag.prod(); 522 } 523 }; 524 525 /** \deprecated use SimplicialLDLT or class SimplicialLLT 526 * \ingroup SparseCholesky_Module 527 * \class SimplicialCholesky 528 * 529 * \sa class SimplicialLDLT, class SimplicialLLT 530 */ 531 template<typename _MatrixType, int _UpLo> 532 class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo> > 533 { 534 public: 535 typedef _MatrixType MatrixType; 536 enum { UpLo = _UpLo }; 537 typedef SimplicialCholeskyBase<SimplicialCholesky> Base; 538 typedef typename MatrixType::Scalar Scalar; 539 typedef typename MatrixType::RealScalar RealScalar; 540 typedef typename MatrixType::Index Index; 541 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; 542 typedef Matrix<Scalar,Dynamic,1> VectorType; 543 typedef internal::traits<SimplicialCholesky> Traits; 544 typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits; 545 typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits; 546 public: 547 SimplicialCholesky() : Base(), m_LDLT(true) {} 548 549 SimplicialCholesky(const MatrixType& matrix) 550 : Base(), m_LDLT(true) 551 { 552 compute(matrix); 553 } 554 555 SimplicialCholesky& setMode(SimplicialCholeskyMode mode) 556 { 557 switch(mode) 558 { 559 case SimplicialCholeskyLLT: 560 m_LDLT = false; 561 break; 562 case SimplicialCholeskyLDLT: 563 m_LDLT = true; 564 break; 565 default: 566 break; 567 } 568 569 return *this; 570 } 571 572 inline const VectorType vectorD() const { 573 eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized"); 574 return Base::m_diag; 575 } 576 inline const CholMatrixType rawMatrix() const { 577 eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized"); 578 return Base::m_matrix; 579 } 580 581 /** Computes the sparse Cholesky decomposition of \a matrix */ 582 SimplicialCholesky& compute(const MatrixType& matrix) 583 { 584 if(m_LDLT) 585 Base::template compute<true>(matrix); 586 else 587 Base::template compute<false>(matrix); 588 return *this; 589 } 590 591 /** Performs a symbolic decomposition on the sparcity of \a matrix. 592 * 593 * This function is particularly useful when solving for several problems having the same structure. 594 * 595 * \sa factorize() 596 */ 597 void analyzePattern(const MatrixType& a) 598 { 599 Base::analyzePattern(a, m_LDLT); 600 } 601 602 /** Performs a numeric decomposition of \a matrix 603 * 604 * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed. 605 * 606 * \sa analyzePattern() 607 */ 608 void factorize(const MatrixType& a) 609 { 610 if(m_LDLT) 611 Base::template factorize<true>(a); 612 else 613 Base::template factorize<false>(a); 614 } 615 616 /** \internal */ 617 template<typename Rhs,typename Dest> 618 void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const 619 { 620 eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 621 eigen_assert(Base::m_matrix.rows()==b.rows()); 622 623 if(Base::m_info!=Success) 624 return; 625 626 if(Base::m_P.size()>0) 627 dest = Base::m_P * b; 628 else 629 dest = b; 630 631 if(Base::m_matrix.nonZeros()>0) // otherwise L==I 632 { 633 if(m_LDLT) 634 LDLTTraits::getL(Base::m_matrix).solveInPlace(dest); 635 else 636 LLTTraits::getL(Base::m_matrix).solveInPlace(dest); 637 } 638 639 if(Base::m_diag.size()>0) 640 dest = Base::m_diag.asDiagonal().inverse() * dest; 641 642 if (Base::m_matrix.nonZeros()>0) // otherwise I==I 643 { 644 if(m_LDLT) 645 LDLTTraits::getU(Base::m_matrix).solveInPlace(dest); 646 else 647 LLTTraits::getU(Base::m_matrix).solveInPlace(dest); 648 } 649 650 if(Base::m_P.size()>0) 651 dest = Base::m_Pinv * dest; 652 } 653 654 Scalar determinant() const 655 { 656 if(m_LDLT) 657 { 658 return Base::m_diag.prod(); 659 } 660 else 661 { 662 Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod(); 663 return internal::abs2(detL); 664 } 665 } 666 667 protected: 668 bool m_LDLT; 669 }; 670 671 template<typename Derived> 672 void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixType& ap) 673 { 674 eigen_assert(a.rows()==a.cols()); 675 const Index size = a.rows(); 676 // TODO allows to configure the permutation 677 // Note that amd compute the inverse permutation 678 { 679 CholMatrixType C; 680 C = a.template selfadjointView<UpLo>(); 681 // remove diagonal entries: 682 // seems not to be needed 683 // C.prune(keep_diag()); 684 internal::minimum_degree_ordering(C, m_Pinv); 685 } 686 687 if(m_Pinv.size()>0) 688 m_P = m_Pinv.inverse(); 689 else 690 m_P.resize(0); 691 692 ap.resize(size,size); 693 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); 694 } 695 696 template<typename Derived> 697 void SimplicialCholeskyBase<Derived>::analyzePattern_preordered(const CholMatrixType& ap, bool doLDLT) 698 { 699 const Index size = ap.rows(); 700 m_matrix.resize(size, size); 701 m_parent.resize(size); 702 m_nonZerosPerCol.resize(size); 703 704 ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0); 705 706 for(Index k = 0; k < size; ++k) 707 { 708 /* L(k,:) pattern: all nodes reachable in etree from nz in A(0:k-1,k) */ 709 m_parent[k] = -1; /* parent of k is not yet known */ 710 tags[k] = k; /* mark node k as visited */ 711 m_nonZerosPerCol[k] = 0; /* count of nonzeros in column k of L */ 712 for(typename CholMatrixType::InnerIterator it(ap,k); it; ++it) 713 { 714 Index i = it.index(); 715 if(i < k) 716 { 717 /* follow path from i to root of etree, stop at flagged node */ 718 for(; tags[i] != k; i = m_parent[i]) 719 { 720 /* find parent of i if not yet determined */ 721 if (m_parent[i] == -1) 722 m_parent[i] = k; 723 m_nonZerosPerCol[i]++; /* L (k,i) is nonzero */ 724 tags[i] = k; /* mark i as visited */ 725 } 726 } 727 } 728 } 729 730 /* construct Lp index array from m_nonZerosPerCol column counts */ 731 Index* Lp = m_matrix.outerIndexPtr(); 732 Lp[0] = 0; 733 for(Index k = 0; k < size; ++k) 734 Lp[k+1] = Lp[k] + m_nonZerosPerCol[k] + (doLDLT ? 0 : 1); 735 736 m_matrix.resizeNonZeros(Lp[size]); 737 738 m_isInitialized = true; 739 m_info = Success; 740 m_analysisIsOk = true; 741 m_factorizationIsOk = false; 742 } 743 744 745 template<typename Derived> 746 template<bool DoLDLT> 747 void SimplicialCholeskyBase<Derived>::factorize_preordered(const CholMatrixType& ap) 748 { 749 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); 750 eigen_assert(ap.rows()==ap.cols()); 751 const Index size = ap.rows(); 752 eigen_assert(m_parent.size()==size); 753 eigen_assert(m_nonZerosPerCol.size()==size); 754 755 const Index* Lp = m_matrix.outerIndexPtr(); 756 Index* Li = m_matrix.innerIndexPtr(); 757 Scalar* Lx = m_matrix.valuePtr(); 758 759 ei_declare_aligned_stack_constructed_variable(Scalar, y, size, 0); 760 ei_declare_aligned_stack_constructed_variable(Index, pattern, size, 0); 761 ei_declare_aligned_stack_constructed_variable(Index, tags, size, 0); 762 763 bool ok = true; 764 m_diag.resize(DoLDLT ? size : 0); 765 766 for(Index k = 0; k < size; ++k) 767 { 768 // compute nonzero pattern of kth row of L, in topological order 769 y[k] = 0.0; // Y(0:k) is now all zero 770 Index top = size; // stack for pattern is empty 771 tags[k] = k; // mark node k as visited 772 m_nonZerosPerCol[k] = 0; // count of nonzeros in column k of L 773 for(typename MatrixType::InnerIterator it(ap,k); it; ++it) 774 { 775 Index i = it.index(); 776 if(i <= k) 777 { 778 y[i] += internal::conj(it.value()); /* scatter A(i,k) into Y (sum duplicates) */ 779 Index len; 780 for(len = 0; tags[i] != k; i = m_parent[i]) 781 { 782 pattern[len++] = i; /* L(k,i) is nonzero */ 783 tags[i] = k; /* mark i as visited */ 784 } 785 while(len > 0) 786 pattern[--top] = pattern[--len]; 787 } 788 } 789 790 /* compute numerical values kth row of L (a sparse triangular solve) */ 791 792 RealScalar d = internal::real(y[k]) * m_shiftScale + m_shiftOffset; // get D(k,k), apply the shift function, and clear Y(k) 793 y[k] = 0.0; 794 for(; top < size; ++top) 795 { 796 Index i = pattern[top]; /* pattern[top:n-1] is pattern of L(:,k) */ 797 Scalar yi = y[i]; /* get and clear Y(i) */ 798 y[i] = 0.0; 799 800 /* the nonzero entry L(k,i) */ 801 Scalar l_ki; 802 if(DoLDLT) 803 l_ki = yi / m_diag[i]; 804 else 805 yi = l_ki = yi / Lx[Lp[i]]; 806 807 Index p2 = Lp[i] + m_nonZerosPerCol[i]; 808 Index p; 809 for(p = Lp[i] + (DoLDLT ? 0 : 1); p < p2; ++p) 810 y[Li[p]] -= internal::conj(Lx[p]) * yi; 811 d -= internal::real(l_ki * internal::conj(yi)); 812 Li[p] = k; /* store L(k,i) in column form of L */ 813 Lx[p] = l_ki; 814 ++m_nonZerosPerCol[i]; /* increment count of nonzeros in col i */ 815 } 816 if(DoLDLT) 817 { 818 m_diag[k] = d; 819 if(d == RealScalar(0)) 820 { 821 ok = false; /* failure, D(k,k) is zero */ 822 break; 823 } 824 } 825 else 826 { 827 Index p = Lp[k] + m_nonZerosPerCol[k]++; 828 Li[p] = k ; /* store L(k,k) = sqrt (d) in column k */ 829 if(d <= RealScalar(0)) { 830 ok = false; /* failure, matrix is not positive definite */ 831 break; 832 } 833 Lx[p] = internal::sqrt(d) ; 834 } 835 } 836 837 m_info = ok ? Success : NumericalIssue; 838 m_factorizationIsOk = true; 839 } 840 841 namespace internal { 842 843 template<typename Derived, typename Rhs> 844 struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs> 845 : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs> 846 { 847 typedef SimplicialCholeskyBase<Derived> Dec; 848 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 849 850 template<typename Dest> void evalTo(Dest& dst) const 851 { 852 dec().derived()._solve(rhs(),dst); 853 } 854 }; 855 856 template<typename Derived, typename Rhs> 857 struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs> 858 : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs> 859 { 860 typedef SimplicialCholeskyBase<Derived> Dec; 861 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) 862 863 template<typename Dest> void evalTo(Dest& dst) const 864 { 865 dec().derived()._solve_sparse(rhs(),dst); 866 } 867 }; 868 869 } // end namespace internal 870 871 } // end namespace Eigen 872 873 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H 874
