$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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 00035 template<typename Derived> 00036 class SimplicialCholeskyBase : internal::noncopyable 00037 { 00038 public: 00039 typedef typename internal::traits<Derived>::MatrixType MatrixType; 00040 typedef typename internal::traits<Derived>::OrderingType OrderingType; 00041 enum { UpLo = internal::traits<Derived>::UpLo }; 00042 typedef typename MatrixType::Scalar Scalar; 00043 typedef typename MatrixType::RealScalar RealScalar; 00044 typedef typename MatrixType::Index Index; 00045 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; 00046 typedef Matrix<Scalar,Dynamic,1> VectorType; 00047 00048 public: 00049 00051 SimplicialCholeskyBase() 00052 : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1) 00053 {} 00054 00055 SimplicialCholeskyBase(const MatrixType& matrix) 00056 : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1) 00057 { 00058 derived().compute(matrix); 00059 } 00060 00061 ~SimplicialCholeskyBase() 00062 { 00063 } 00064 00065 Derived& derived() { return *static_cast<Derived*>(this); } 00066 const Derived& derived() const { return *static_cast<const Derived*>(this); } 00067 00068 inline Index cols() const { return m_matrix.cols(); } 00069 inline Index rows() const { return m_matrix.rows(); } 00070 00076 ComputationInfo info() const 00077 { 00078 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 00079 return m_info; 00080 } 00081 00086 template<typename Rhs> 00087 inline const internal::solve_retval<SimplicialCholeskyBase, Rhs> 00088 solve(const MatrixBase<Rhs>& b) const 00089 { 00090 eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized."); 00091 eigen_assert(rows()==b.rows() 00092 && "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b"); 00093 return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived()); 00094 } 00095 00100 template<typename Rhs> 00101 inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs> 00102 solve(const SparseMatrixBase<Rhs>& b) const 00103 { 00104 eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized."); 00105 eigen_assert(rows()==b.rows() 00106 && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b"); 00107 return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived()); 00108 } 00109 00112 const PermutationMatrix<Dynamic,Dynamic,Index>& permutationP() const 00113 { return m_P; } 00114 00117 const PermutationMatrix<Dynamic,Dynamic,Index>& permutationPinv() const 00118 { return m_Pinv; } 00119 00129 Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1) 00130 { 00131 m_shiftOffset = offset; 00132 m_shiftScale = scale; 00133 return derived(); 00134 } 00135 00136 #ifndef EIGEN_PARSED_BY_DOXYGEN 00137 00138 template<typename Stream> 00139 void dumpMemory(Stream& s) 00140 { 00141 int total = 0; 00142 s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n"; 00143 s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n"; 00144 s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 00145 s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 00146 s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 00147 s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n"; 00148 s << " TOTAL: " << (total>> 20) << "Mb" << "\n"; 00149 } 00150 00152 template<typename Rhs,typename Dest> 00153 void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const 00154 { 00155 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 00156 eigen_assert(m_matrix.rows()==b.rows()); 00157 00158 if(m_info!=Success) 00159 return; 00160 00161 if(m_P.size()>0) 00162 dest = m_P * b; 00163 else 00164 dest = b; 00165 00166 if(m_matrix.nonZeros()>0) // otherwise L==I 00167 derived().matrixL().solveInPlace(dest); 00168 00169 if(m_diag.size()>0) 00170 dest = m_diag.asDiagonal().inverse() * dest; 00171 00172 if (m_matrix.nonZeros()>0) // otherwise U==I 00173 derived().matrixU().solveInPlace(dest); 00174 00175 if(m_P.size()>0) 00176 dest = m_Pinv * dest; 00177 } 00178 00179 #endif // EIGEN_PARSED_BY_DOXYGEN 00180 00181 protected: 00182 00184 template<bool DoLDLT> 00185 void compute(const MatrixType& matrix) 00186 { 00187 eigen_assert(matrix.rows()==matrix.cols()); 00188 Index size = matrix.cols(); 00189 CholMatrixType ap(size,size); 00190 ordering(matrix, ap); 00191 analyzePattern_preordered(ap, DoLDLT); 00192 factorize_preordered<DoLDLT>(ap); 00193 } 00194 00195 template<bool DoLDLT> 00196 void factorize(const MatrixType& a) 00197 { 00198 eigen_assert(a.rows()==a.cols()); 00199 int size = a.cols(); 00200 CholMatrixType ap(size,size); 00201 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); 00202 factorize_preordered<DoLDLT>(ap); 00203 } 00204 00205 template<bool DoLDLT> 00206 void factorize_preordered(const CholMatrixType& a); 00207 00208 void analyzePattern(const MatrixType& a, bool doLDLT) 00209 { 00210 eigen_assert(a.rows()==a.cols()); 00211 int size = a.cols(); 00212 CholMatrixType ap(size,size); 00213 ordering(a, ap); 00214 analyzePattern_preordered(ap,doLDLT); 00215 } 00216 void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT); 00217 00218 void ordering(const MatrixType& a, CholMatrixType& ap); 00219 00221 struct keep_diag { 00222 inline bool operator() (const Index& row, const Index& col, const Scalar&) const 00223 { 00224 return row!=col; 00225 } 00226 }; 00227 00228 mutable ComputationInfo m_info; 00229 bool m_isInitialized; 00230 bool m_factorizationIsOk; 00231 bool m_analysisIsOk; 00232 00233 CholMatrixType m_matrix; 00234 VectorType m_diag; // the diagonal coefficients (LDLT mode) 00235 VectorXi m_parent; // elimination tree 00236 VectorXi m_nonZerosPerCol; 00237 PermutationMatrix<Dynamic,Dynamic,Index> m_P; // the permutation 00238 PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // the inverse permutation 00239 00240 RealScalar m_shiftOffset; 00241 RealScalar m_shiftScale; 00242 }; 00243 00244 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLLT; 00245 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLDLT; 00246 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialCholesky; 00247 00248 namespace internal { 00249 00250 template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialLLT<_MatrixType,_UpLo,_Ordering> > 00251 { 00252 typedef _MatrixType MatrixType; 00253 typedef _Ordering OrderingType; 00254 enum { UpLo = _UpLo }; 00255 typedef typename MatrixType::Scalar Scalar; 00256 typedef typename MatrixType::Index Index; 00257 typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType; 00258 typedef SparseTriangularView<CholMatrixType, Eigen::Lower> MatrixL; 00259 typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU; 00260 static inline MatrixL getL(const MatrixType& m) { return m; } 00261 static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } 00262 }; 00263 00264 template<typename _MatrixType,int _UpLo, typename _Ordering> struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> > 00265 { 00266 typedef _MatrixType MatrixType; 00267 typedef _Ordering OrderingType; 00268 enum { UpLo = _UpLo }; 00269 typedef typename MatrixType::Scalar Scalar; 00270 typedef typename MatrixType::Index Index; 00271 typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType; 00272 typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower> MatrixL; 00273 typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU; 00274 static inline MatrixL getL(const MatrixType& m) { return m; } 00275 static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); } 00276 }; 00277 00278 template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> > 00279 { 00280 typedef _MatrixType MatrixType; 00281 typedef _Ordering OrderingType; 00282 enum { UpLo = _UpLo }; 00283 }; 00284 00285 } 00286 00305 template<typename _MatrixType, int _UpLo, typename _Ordering> 00306 class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> > 00307 { 00308 public: 00309 typedef _MatrixType MatrixType; 00310 enum { UpLo = _UpLo }; 00311 typedef SimplicialCholeskyBase<SimplicialLLT> Base; 00312 typedef typename MatrixType::Scalar Scalar; 00313 typedef typename MatrixType::RealScalar RealScalar; 00314 typedef typename MatrixType::Index Index; 00315 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; 00316 typedef Matrix<Scalar,Dynamic,1> VectorType; 00317 typedef internal::traits<SimplicialLLT> Traits; 00318 typedef typename Traits::MatrixL MatrixL; 00319 typedef typename Traits::MatrixU MatrixU; 00320 public: 00322 SimplicialLLT() : Base() {} 00324 SimplicialLLT(const MatrixType& matrix) 00325 : Base(matrix) {} 00326 00328 inline const MatrixL matrixL() const { 00329 eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized"); 00330 return Traits::getL(Base::m_matrix); 00331 } 00332 00334 inline const MatrixU matrixU() const { 00335 eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized"); 00336 return Traits::getU(Base::m_matrix); 00337 } 00338 00340 SimplicialLLT& compute(const MatrixType& matrix) 00341 { 00342 Base::template compute<false>(matrix); 00343 return *this; 00344 } 00345 00352 void analyzePattern(const MatrixType& a) 00353 { 00354 Base::analyzePattern(a, false); 00355 } 00356 00363 void factorize(const MatrixType& a) 00364 { 00365 Base::template factorize<false>(a); 00366 } 00367 00369 Scalar determinant() const 00370 { 00371 Scalar detL = Base::m_matrix.diagonal().prod(); 00372 return numext::abs2(detL); 00373 } 00374 }; 00375 00394 template<typename _MatrixType, int _UpLo, typename _Ordering> 00395 class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> > 00396 { 00397 public: 00398 typedef _MatrixType MatrixType; 00399 enum { UpLo = _UpLo }; 00400 typedef SimplicialCholeskyBase<SimplicialLDLT> Base; 00401 typedef typename MatrixType::Scalar Scalar; 00402 typedef typename MatrixType::RealScalar RealScalar; 00403 typedef typename MatrixType::Index Index; 00404 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; 00405 typedef Matrix<Scalar,Dynamic,1> VectorType; 00406 typedef internal::traits<SimplicialLDLT> Traits; 00407 typedef typename Traits::MatrixL MatrixL; 00408 typedef typename Traits::MatrixU MatrixU; 00409 public: 00411 SimplicialLDLT() : Base() {} 00412 00414 SimplicialLDLT(const MatrixType& matrix) 00415 : Base(matrix) {} 00416 00418 inline const VectorType vectorD() const { 00419 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); 00420 return Base::m_diag; 00421 } 00423 inline const MatrixL matrixL() const { 00424 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); 00425 return Traits::getL(Base::m_matrix); 00426 } 00427 00429 inline const MatrixU matrixU() const { 00430 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized"); 00431 return Traits::getU(Base::m_matrix); 00432 } 00433 00435 SimplicialLDLT& compute(const MatrixType& matrix) 00436 { 00437 Base::template compute<true>(matrix); 00438 return *this; 00439 } 00440 00447 void analyzePattern(const MatrixType& a) 00448 { 00449 Base::analyzePattern(a, true); 00450 } 00451 00458 void factorize(const MatrixType& a) 00459 { 00460 Base::template factorize<true>(a); 00461 } 00462 00464 Scalar determinant() const 00465 { 00466 return Base::m_diag.prod(); 00467 } 00468 }; 00469 00476 template<typename _MatrixType, int _UpLo, typename _Ordering> 00477 class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> > 00478 { 00479 public: 00480 typedef _MatrixType MatrixType; 00481 enum { UpLo = _UpLo }; 00482 typedef SimplicialCholeskyBase<SimplicialCholesky> Base; 00483 typedef typename MatrixType::Scalar Scalar; 00484 typedef typename MatrixType::RealScalar RealScalar; 00485 typedef typename MatrixType::Index Index; 00486 typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType; 00487 typedef Matrix<Scalar,Dynamic,1> VectorType; 00488 typedef internal::traits<SimplicialCholesky> Traits; 00489 typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits; 00490 typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits; 00491 public: 00492 SimplicialCholesky() : Base(), m_LDLT(true) {} 00493 00494 SimplicialCholesky(const MatrixType& matrix) 00495 : Base(), m_LDLT(true) 00496 { 00497 compute(matrix); 00498 } 00499 00500 SimplicialCholesky& setMode(SimplicialCholeskyMode mode) 00501 { 00502 switch(mode) 00503 { 00504 case SimplicialCholeskyLLT: 00505 m_LDLT = false; 00506 break; 00507 case SimplicialCholeskyLDLT: 00508 m_LDLT = true; 00509 break; 00510 default: 00511 break; 00512 } 00513 00514 return *this; 00515 } 00516 00517 inline const VectorType vectorD() const { 00518 eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized"); 00519 return Base::m_diag; 00520 } 00521 inline const CholMatrixType rawMatrix() const { 00522 eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized"); 00523 return Base::m_matrix; 00524 } 00525 00527 SimplicialCholesky& compute(const MatrixType& matrix) 00528 { 00529 if(m_LDLT) 00530 Base::template compute<true>(matrix); 00531 else 00532 Base::template compute<false>(matrix); 00533 return *this; 00534 } 00535 00542 void analyzePattern(const MatrixType& a) 00543 { 00544 Base::analyzePattern(a, m_LDLT); 00545 } 00546 00553 void factorize(const MatrixType& a) 00554 { 00555 if(m_LDLT) 00556 Base::template factorize<true>(a); 00557 else 00558 Base::template factorize<false>(a); 00559 } 00560 00562 template<typename Rhs,typename Dest> 00563 void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const 00564 { 00565 eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 00566 eigen_assert(Base::m_matrix.rows()==b.rows()); 00567 00568 if(Base::m_info!=Success) 00569 return; 00570 00571 if(Base::m_P.size()>0) 00572 dest = Base::m_P * b; 00573 else 00574 dest = b; 00575 00576 if(Base::m_matrix.nonZeros()>0) // otherwise L==I 00577 { 00578 if(m_LDLT) 00579 LDLTTraits::getL(Base::m_matrix).solveInPlace(dest); 00580 else 00581 LLTTraits::getL(Base::m_matrix).solveInPlace(dest); 00582 } 00583 00584 if(Base::m_diag.size()>0) 00585 dest = Base::m_diag.asDiagonal().inverse() * dest; 00586 00587 if (Base::m_matrix.nonZeros()>0) // otherwise I==I 00588 { 00589 if(m_LDLT) 00590 LDLTTraits::getU(Base::m_matrix).solveInPlace(dest); 00591 else 00592 LLTTraits::getU(Base::m_matrix).solveInPlace(dest); 00593 } 00594 00595 if(Base::m_P.size()>0) 00596 dest = Base::m_Pinv * dest; 00597 } 00598 00599 Scalar determinant() const 00600 { 00601 if(m_LDLT) 00602 { 00603 return Base::m_diag.prod(); 00604 } 00605 else 00606 { 00607 Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod(); 00608 return numext::abs2(detL); 00609 } 00610 } 00611 00612 protected: 00613 bool m_LDLT; 00614 }; 00615 00616 template<typename Derived> 00617 void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixType& ap) 00618 { 00619 eigen_assert(a.rows()==a.cols()); 00620 const Index size = a.rows(); 00621 // Note that amd compute the inverse permutation 00622 { 00623 CholMatrixType C; 00624 C = a.template selfadjointView<UpLo>(); 00625 00626 OrderingType ordering; 00627 ordering(C,m_Pinv); 00628 } 00629 00630 if(m_Pinv.size()>0) 00631 m_P = m_Pinv.inverse(); 00632 else 00633 m_P.resize(0); 00634 00635 ap.resize(size,size); 00636 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P); 00637 } 00638 00639 namespace internal { 00640 00641 template<typename Derived, typename Rhs> 00642 struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs> 00643 : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs> 00644 { 00645 typedef SimplicialCholeskyBase<Derived> Dec; 00646 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 00647 00648 template<typename Dest> void evalTo(Dest& dst) const 00649 { 00650 dec().derived()._solve(rhs(),dst); 00651 } 00652 }; 00653 00654 template<typename Derived, typename Rhs> 00655 struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs> 00656 : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs> 00657 { 00658 typedef SimplicialCholeskyBase<Derived> Dec; 00659 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) 00660 00661 template<typename Dest> void evalTo(Dest& dst) const 00662 { 00663 this->defaultEvalTo(dst); 00664 } 00665 }; 00666 00667 } // end namespace internal 00668 00669 } // end namespace Eigen 00670 00671 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H