![]() |
Eigen
3.3.3
|
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2008-2012 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_SIMPLICIAL_CHOLESKY_H 00011 #define EIGEN_SIMPLICIAL_CHOLESKY_H 00012 00013 namespace Eigen { 00014 00015 enum SimplicialCholeskyMode { 00016 SimplicialCholeskyLLT, 00017 SimplicialCholeskyLDLT 00018 }; 00019 00020 namespace internal { 00021 template<typename CholMatrixType, typename InputMatrixType> 00022 struct simplicial_cholesky_grab_input { 00023 typedef CholMatrixType const * ConstCholMatrixPtr; 00024 static void run(const InputMatrixType& input, ConstCholMatrixPtr &pmat, CholMatrixType &tmp) 00025 { 00026 tmp = input; 00027 pmat = &tmp; 00028 } 00029 }; 00030 00031 template<typename MatrixType> 00032 struct simplicial_cholesky_grab_input<MatrixType,MatrixType> { 00033 typedef MatrixType const * ConstMatrixPtr; 00034 static void run(const MatrixType& input, ConstMatrixPtr &pmat, MatrixType &/*tmp*/) 00035 { 00036 pmat = &input; 00037 } 00038 }; 00039 } // end namespace internal 00040 00054 template<typename Derived> 00055 class SimplicialCholeskyBase : public SparseSolverBase<Derived> 00056 { 00057 typedef SparseSolverBase<Derived> Base; 00058 using Base::m_isInitialized; 00059 00060 public: 00061 typedef typename internal::traits<Derived>::MatrixType MatrixType; 00062 typedef typename internal::traits<Derived>::OrderingType OrderingType; 00063 enum { UpLo = internal::traits<Derived>::UpLo }; 00064 typedef typename MatrixType::Scalar Scalar; 00065 typedef typename MatrixType::RealScalar RealScalar; 00066 typedef typename MatrixType::StorageIndex StorageIndex; 00067 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType; 00068 typedef CholMatrixType const * ConstCholMatrixPtr; 00069 typedef Matrix<Scalar,Dynamic,1> VectorType; 00070 typedef Matrix<StorageIndex,Dynamic,1> VectorI; 00071 00072 enum { 00073 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00074 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00075 }; 00076 00077 public: 00078 00079 using Base::derived; 00080 00082 SimplicialCholeskyBase() 00083 : m_info(Success), m_shiftOffset(0), m_shiftScale(1) 00084 {} 00085 00086 explicit SimplicialCholeskyBase(const MatrixType& matrix) 00087 : m_info(Success), m_shiftOffset(0), m_shiftScale(1) 00088 { 00089 derived().compute(matrix); 00090 } 00091 00092 ~SimplicialCholeskyBase() 00093 { 00094 } 00095 00096 Derived& derived() { return *static_cast<Derived*>(this); } 00097 const Derived& derived() const { return *static_cast<const Derived*>(this); } 00098 00099 inline Index cols() const { return m_matrix.cols(); } 00100 inline Index rows() const { return m_matrix.rows(); } 00101 00107 ComputationInfo info() const 00108 { 00109 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 00110 return m_info; 00111 } 00112 00115 const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& permutationP() const 00116 { return m_P; } 00117 00120 const PermutationMatrix<Dynamic,Dynamic,StorageIndex>& permutationPinv() const 00121 { return m_Pinv; } 00122 00132 Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1) 00133 { 00134 m_shiftOffset = offset; 00135 m_shiftScale = scale; 00136 return derived(); 00137 } 00138 00139 #ifndef EIGEN_PARSED_BY_DOXYGEN 00140 00141 template<typename Stream> 00142 void dumpMemory(Stream& s) 00143 { 00144 int total = 0; 00145 s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n"; 00146 s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n"; 00147 s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 00148 s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 00149 s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 00150 s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 00151 s << " TOTAL: " << (total>> 20) << "Mb" << "\n"; 00152 } 00153 00155 template<typename Rhs,typename Dest> 00156 void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const 00157 { 00158 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 00159 eigen_assert(m_matrix.rows()==b.rows()); 00160 00161 if(m_info!=Success) 00162 return; 00163 00164 if(m_P.size()>0) 00165 dest = m_P * b; 00166 else 00167 dest = b; 00168 00169 if(m_matrix.nonZeros()>0) // otherwise L==I 00170 derived().matrixL().solveInPlace(dest); 00171 00172 if(m_diag.size()>0) 00173 dest = m_diag.asDiagonal().inverse() * dest; 00174 00175 if (m_matrix.nonZeros()>0) // otherwise U==I 00176 derived().matrixU().solveInPlace(dest); 00177 00178 if(m_P.size()>0) 00179 dest = m_Pinv * dest; 00180 } 00181 00182 template<typename Rhs,typename Dest> 00183 void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const 00184 { 00185 internal::solve_sparse_through_dense_panels(derived(), b, dest); 00186 } 00187 00188 #endif // EIGEN_PARSED_BY_DOXYGEN 00189 00190 protected: 00191 00193 template<bool DoLDLT> 00194 void compute(const MatrixType& matrix) 00195 { 00196 eigen_assert(matrix.rows()==matrix.cols()); 00197 Index size = matrix.cols(); 00198 CholMatrixType tmp(size,size); 00199 ConstCholMatrixPtr pmat; 00200 ordering(matrix, pmat, tmp); 00201 analyzePattern_preordered(*pmat, DoLDLT); 00202 factorize_preordered<DoLDLT>(*pmat); 00203 } 00204 00205 template<bool DoLDLT> 00206 void factorize(const MatrixType& a) 00207 { 00208 eigen_assert(a.rows()==a.cols()); 00209 Index size = a.cols(); 00210 CholMatrixType tmp(size,size); 00211 ConstCholMatrixPtr pmat; 00212 00213 if(m_P.size()==0 && (UpLo&Upper)==Upper) 00214 { 00215 // If there is no ordering, try to directly use the input matrix without any copy 00216 internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp); 00217 } 00218 else 00219 { 00220 tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); 00221 pmat = &tmp; 00222 } 00223 00224 factorize_preordered<DoLDLT>(*pmat); 00225 } 00226 00227 template<bool DoLDLT> 00228 void factorize_preordered(const CholMatrixType& a); 00229 00230 void analyzePattern(const MatrixType& a, bool doLDLT) 00231 { 00232 eigen_assert(a.rows()==a.cols()); 00233 Index size = a.cols(); 00234 CholMatrixType tmp(size,size); 00235 ConstCholMatrixPtr pmat; 00236 ordering(a, pmat, tmp); 00237 analyzePattern_preordered(*pmat,doLDLT); 00238 } 00239 void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT); 00240 00241 void ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap); 00242 00244 struct keep_diag { 00245 inline bool operator() (const Index& row, const Index& col, const Scalar&) const 00246 { 00247 return row!=col; 00248 } 00249 }; 00250 00251 mutable ComputationInfo m_info; 00252 bool m_factorizationIsOk; 00253 bool m_analysisIsOk; 00254 00255 CholMatrixType m_matrix; 00256 VectorType m_diag; // the diagonal coefficients (LDLT mode) 00257 VectorI m_parent; // elimination tree 00258 VectorI m_nonZerosPerCol; 00259 PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_P; // the permutation 00260 PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_Pinv; // the inverse permutation 00261 00262 RealScalar m_shiftOffset; 00263 RealScalar m_shiftScale; 00264 }; 00265 00266 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLLT; 00267 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLDLT; 00268 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialCholesky; 00269 00270 namespace internal { 00271 00272 template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialLLT<_MatrixType,_UpLo,_Ordering> > 00273 { 00274 typedef _MatrixType MatrixType; 00275 typedef _Ordering OrderingType; 00276 enum { UpLo = _UpLo }; 00277 typedef typename MatrixType::Scalar Scalar; 00278 typedef typename MatrixType::StorageIndex StorageIndex; 00279 typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType; 00280 typedef TriangularView<const CholMatrixType, Eigen::Lower> MatrixL; 00281 typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU; 00282 static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); } 00283 static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); } 00284 }; 00285 00286 template<typename _MatrixType,int _UpLo, typename _Ordering> struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> > 00287 { 00288 typedef _MatrixType MatrixType; 00289 typedef _Ordering OrderingType; 00290 enum { UpLo = _UpLo }; 00291 typedef typename MatrixType::Scalar Scalar; 00292 typedef typename MatrixType::StorageIndex StorageIndex; 00293 typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType; 00294 typedef TriangularView<const CholMatrixType, Eigen::UnitLower> MatrixL; 00295 typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU; 00296 static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); } 00297 static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); } 00298 }; 00299 00300 template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> > 00301 { 00302 typedef _MatrixType MatrixType; 00303 typedef _Ordering OrderingType; 00304 enum { UpLo = _UpLo }; 00305 }; 00306 00307 } 00308 00329 template<typename _MatrixType, int _UpLo, typename _Ordering> 00330 class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> > 00331 { 00332 public: 00333 typedef _MatrixType MatrixType; 00334 enum { UpLo = _UpLo }; 00335 typedef SimplicialCholeskyBase<SimplicialLLT> Base; 00336 typedef typename MatrixType::Scalar Scalar; 00337 typedef typename MatrixType::RealScalar RealScalar; 00338 typedef typename MatrixType::StorageIndex StorageIndex; 00339 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; 00340 typedef Matrix<Scalar,Dynamic,1> VectorType; 00341 typedef internal::traits<SimplicialLLT> Traits; 00342 typedef typename Traits::MatrixL MatrixL; 00343 typedef typename Traits::MatrixU MatrixU; 00344 public: 00346 SimplicialLLT() : Base() {} 00348 explicit SimplicialLLT(const MatrixType& matrix) 00349 : Base(matrix) {} 00350 00352 inline const MatrixL matrixL() const { 00353 eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized"); 00354 return Traits::getL(Base::m_matrix); 00355 } 00356 00358 inline const MatrixU matrixU() const { 00359 eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized"); 00360 return Traits::getU(Base::m_matrix); 00361 } 00362 00364 SimplicialLLT& compute(const MatrixType& matrix) 00365 { 00366 Base::template compute<false>(matrix); 00367 return *this; 00368 } 00369 00376 void analyzePattern(const MatrixType& a) 00377 { 00378 Base::analyzePattern(a, false); 00379 } 00380 00387 void factorize(const MatrixType& a) 00388 { 00389 Base::template factorize<false>(a); 00390 } 00391 00393 Scalar determinant() const 00394 { 00395 Scalar detL = Base::m_matrix.diagonal().prod(); 00396 return numext::abs2(detL); 00397 } 00398 }; 00399 00420 template<typename _MatrixType, int _UpLo, typename _Ordering> 00421 class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> > 00422 { 00423 public: 00424 typedef _MatrixType MatrixType; 00425 enum { UpLo = _UpLo }; 00426 typedef SimplicialCholeskyBase<SimplicialLDLT> Base; 00427 typedef typename MatrixType::Scalar Scalar; 00428 typedef typename MatrixType::RealScalar RealScalar; 00429 typedef typename MatrixType::StorageIndex StorageIndex; 00430 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType; 00431 typedef Matrix<Scalar,Dynamic,1> VectorType; 00432 typedef internal::traits<SimplicialLDLT> Traits; 00433 typedef typename Traits::MatrixL MatrixL; 00434 typedef typename Traits::MatrixU MatrixU; 00435 public: 00437 SimplicialLDLT() : Base() {} 00438 00440 explicit SimplicialLDLT(const MatrixType& matrix) 00441 : Base(matrix) {} 00442 00444 inline const VectorType vectorD() const { 00445 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); 00446 return Base::m_diag; 00447 } 00449 inline const MatrixL matrixL() const { 00450 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); 00451 return Traits::getL(Base::m_matrix); 00452 } 00453 00455 inline const MatrixU matrixU() const { 00456 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); 00457 return Traits::getU(Base::m_matrix); 00458 } 00459 00461 SimplicialLDLT& compute(const MatrixType& matrix) 00462 { 00463 Base::template compute<true>(matrix); 00464 return *this; 00465 } 00466 00473 void analyzePattern(const MatrixType& a) 00474 { 00475 Base::analyzePattern(a, true); 00476 } 00477 00484 void factorize(const MatrixType& a) 00485 { 00486 Base::template factorize<true>(a); 00487 } 00488 00490 Scalar determinant() const 00491 { 00492 return Base::m_diag.prod(); 00493 } 00494 }; 00495 00502 template<typename _MatrixType, int _UpLo, typename _Ordering> 00503 class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> > 00504 { 00505 public: 00506 typedef _MatrixType MatrixType; 00507 enum { UpLo = _UpLo }; 00508 typedef SimplicialCholeskyBase<SimplicialCholesky> Base; 00509 typedef typename MatrixType::Scalar Scalar; 00510 typedef typename MatrixType::RealScalar RealScalar; 00511 typedef typename MatrixType::StorageIndex StorageIndex; 00512 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType; 00513 typedef Matrix<Scalar,Dynamic,1> VectorType; 00514 typedef internal::traits<SimplicialCholesky> Traits; 00515 typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits; 00516 typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits; 00517 public: 00518 SimplicialCholesky() : Base(), m_LDLT(true) {} 00519 00520 explicit SimplicialCholesky(const MatrixType& matrix) 00521 : Base(), m_LDLT(true) 00522 { 00523 compute(matrix); 00524 } 00525 00526 SimplicialCholesky& setMode(SimplicialCholeskyMode mode) 00527 { 00528 switch(mode) 00529 { 00530 case SimplicialCholeskyLLT: 00531 m_LDLT = false; 00532 break; 00533 case SimplicialCholeskyLDLT: 00534 m_LDLT = true; 00535 break; 00536 default: 00537 break; 00538 } 00539 00540 return *this; 00541 } 00542 00543 inline const VectorType vectorD() const { 00544 eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized"); 00545 return Base::m_diag; 00546 } 00547 inline const CholMatrixType rawMatrix() const { 00548 eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized"); 00549 return Base::m_matrix; 00550 } 00551 00553 SimplicialCholesky& compute(const MatrixType& matrix) 00554 { 00555 if(m_LDLT) 00556 Base::template compute<true>(matrix); 00557 else 00558 Base::template compute<false>(matrix); 00559 return *this; 00560 } 00561 00568 void analyzePattern(const MatrixType& a) 00569 { 00570 Base::analyzePattern(a, m_LDLT); 00571 } 00572 00579 void factorize(const MatrixType& a) 00580 { 00581 if(m_LDLT) 00582 Base::template factorize<true>(a); 00583 else 00584 Base::template factorize<false>(a); 00585 } 00586 00588 template<typename Rhs,typename Dest> 00589 void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const 00590 { 00591 eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 00592 eigen_assert(Base::m_matrix.rows()==b.rows()); 00593 00594 if(Base::m_info!=Success) 00595 return; 00596 00597 if(Base::m_P.size()>0) 00598 dest = Base::m_P * b; 00599 else 00600 dest = b; 00601 00602 if(Base::m_matrix.nonZeros()>0) // otherwise L==I 00603 { 00604 if(m_LDLT) 00605 LDLTTraits::getL(Base::m_matrix).solveInPlace(dest); 00606 else 00607 LLTTraits::getL(Base::m_matrix).solveInPlace(dest); 00608 } 00609 00610 if(Base::m_diag.size()>0) 00611 dest = Base::m_diag.asDiagonal().inverse() * dest; 00612 00613 if (Base::m_matrix.nonZeros()>0) // otherwise I==I 00614 { 00615 if(m_LDLT) 00616 LDLTTraits::getU(Base::m_matrix).solveInPlace(dest); 00617 else 00618 LLTTraits::getU(Base::m_matrix).solveInPlace(dest); 00619 } 00620 00621 if(Base::m_P.size()>0) 00622 dest = Base::m_Pinv * dest; 00623 } 00624 00626 template<typename Rhs,typename Dest> 00627 void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const 00628 { 00629 internal::solve_sparse_through_dense_panels(*this, b, dest); 00630 } 00631 00632 Scalar determinant() const 00633 { 00634 if(m_LDLT) 00635 { 00636 return Base::m_diag.prod(); 00637 } 00638 else 00639 { 00640 Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod(); 00641 return numext::abs2(detL); 00642 } 00643 } 00644 00645 protected: 00646 bool m_LDLT; 00647 }; 00648 00649 template<typename Derived> 00650 void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap) 00651 { 00652 eigen_assert(a.rows()==a.cols()); 00653 const Index size = a.rows(); 00654 pmat = ≈ 00655 // Note that ordering methods compute the inverse permutation 00656 if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value) 00657 { 00658 { 00659 CholMatrixType C; 00660 C = a.template selfadjointView<UpLo>(); 00661 00662 OrderingType ordering; 00663 ordering(C,m_Pinv); 00664 } 00665 00666 if(m_Pinv.size()>0) m_P = m_Pinv.inverse(); 00667 else m_P.resize(0); 00668 00669 ap.resize(size,size); 00670 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); 00671 } 00672 else 00673 { 00674 m_Pinv.resize(0); 00675 m_P.resize(0); 00676 if(int(UpLo)==int(Lower) || MatrixType::IsRowMajor) 00677 { 00678 // we have to transpose the lower part to to the upper one 00679 ap.resize(size,size); 00680 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>(); 00681 } 00682 else 00683 internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap); 00684 } 00685 } 00686 00687 } // end namespace Eigen 00688 00689 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H