SimplicialCholesky.h
Go to the documentation of this file.
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     enum { UpLo = internal::traits<Derived>::UpLo };
00041     typedef typename MatrixType::Scalar Scalar;
00042     typedef typename MatrixType::RealScalar RealScalar;
00043     typedef typename MatrixType::Index Index;
00044     typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
00045     typedef Matrix<Scalar,Dynamic,1> VectorType;
00046 
00047   public:
00048 
00050     SimplicialCholeskyBase()
00051       : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
00052     {}
00053 
00054     SimplicialCholeskyBase(const MatrixType& matrix)
00055       : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
00056     {
00057       derived().compute(matrix);
00058     }
00059 
00060     ~SimplicialCholeskyBase()
00061     {
00062     }
00063 
00064     Derived& derived() { return *static_cast<Derived*>(this); }
00065     const Derived& derived() const { return *static_cast<const Derived*>(this); }
00066     
00067     inline Index cols() const { return m_matrix.cols(); }
00068     inline Index rows() const { return m_matrix.rows(); }
00069     
00075     ComputationInfo info() const
00076     {
00077       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00078       return m_info;
00079     }
00080     
00085     template<typename Rhs>
00086     inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
00087     solve(const MatrixBase<Rhs>& b) const
00088     {
00089       eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
00090       eigen_assert(rows()==b.rows()
00091                 && "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
00092       return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
00093     }
00094     
00099     template<typename Rhs>
00100     inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
00101     solve(const SparseMatrixBase<Rhs>& b) const
00102     {
00103       eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
00104       eigen_assert(rows()==b.rows()
00105                 && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
00106       return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
00107     }
00108     
00111     const PermutationMatrix<Dynamic,Dynamic,Index>& permutationP() const
00112     { return m_P; }
00113     
00116     const PermutationMatrix<Dynamic,Dynamic,Index>& permutationPinv() const
00117     { return m_Pinv; }
00118 
00128     Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1)
00129     {
00130       m_shiftOffset = offset;
00131       m_shiftScale = scale;
00132       return derived();
00133     }
00134 
00135 #ifndef EIGEN_PARSED_BY_DOXYGEN
00136 
00137     template<typename Stream>
00138     void dumpMemory(Stream& s)
00139     {
00140       int total = 0;
00141       s << "  L:        " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
00142       s << "  diag:     " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
00143       s << "  tree:     " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
00144       s << "  nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
00145       s << "  perm:     " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
00146       s << "  perm^-1:  " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
00147       s << "  TOTAL:    " << (total>> 20) << "Mb" << "\n";
00148     }
00149 
00151     template<typename Rhs,typename Dest>
00152     void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
00153     {
00154       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00155       eigen_assert(m_matrix.rows()==b.rows());
00156 
00157       if(m_info!=Success)
00158         return;
00159 
00160       if(m_P.size()>0)
00161         dest = m_P * b;
00162       else
00163         dest = b;
00164 
00165       if(m_matrix.nonZeros()>0) // otherwise L==I
00166         derived().matrixL().solveInPlace(dest);
00167 
00168       if(m_diag.size()>0)
00169         dest = m_diag.asDiagonal().inverse() * dest;
00170 
00171       if (m_matrix.nonZeros()>0) // otherwise U==I
00172         derived().matrixU().solveInPlace(dest);
00173 
00174       if(m_P.size()>0)
00175         dest = m_Pinv * dest;
00176     }
00177 
00178 #endif // EIGEN_PARSED_BY_DOXYGEN
00179 
00180   protected:
00181     
00183     template<bool DoLDLT>
00184     void compute(const MatrixType& matrix)
00185     {
00186       eigen_assert(matrix.rows()==matrix.cols());
00187       Index size = matrix.cols();
00188       CholMatrixType ap(size,size);
00189       ordering(matrix, ap);
00190       analyzePattern_preordered(ap, DoLDLT);
00191       factorize_preordered<DoLDLT>(ap);
00192     }
00193     
00194     template<bool DoLDLT>
00195     void factorize(const MatrixType& a)
00196     {
00197       eigen_assert(a.rows()==a.cols());
00198       int size = a.cols();
00199       CholMatrixType ap(size,size);
00200       ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
00201       factorize_preordered<DoLDLT>(ap);
00202     }
00203 
00204     template<bool DoLDLT>
00205     void factorize_preordered(const CholMatrixType& a);
00206 
00207     void analyzePattern(const MatrixType& a, bool doLDLT)
00208     {
00209       eigen_assert(a.rows()==a.cols());
00210       int size = a.cols();
00211       CholMatrixType ap(size,size);
00212       ordering(a, ap);
00213       analyzePattern_preordered(ap,doLDLT);
00214     }
00215     void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
00216     
00217     void ordering(const MatrixType& a, CholMatrixType& ap);
00218 
00220     struct keep_diag {
00221       inline bool operator() (const Index& row, const Index& col, const Scalar&) const
00222       {
00223         return row!=col;
00224       }
00225     };
00226 
00227     mutable ComputationInfo m_info;
00228     bool m_isInitialized;
00229     bool m_factorizationIsOk;
00230     bool m_analysisIsOk;
00231     
00232     CholMatrixType m_matrix;
00233     VectorType m_diag;                                // the diagonal coefficients (LDLT mode)
00234     VectorXi m_parent;                                // elimination tree
00235     VectorXi m_nonZerosPerCol;
00236     PermutationMatrix<Dynamic,Dynamic,Index> m_P;     // the permutation
00237     PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv;  // the inverse permutation
00238 
00239     RealScalar m_shiftOffset;
00240     RealScalar m_shiftScale;
00241 };
00242 
00243 template<typename _MatrixType, int _UpLo = Lower> class SimplicialLLT;
00244 template<typename _MatrixType, int _UpLo = Lower> class SimplicialLDLT;
00245 template<typename _MatrixType, int _UpLo = Lower> class SimplicialCholesky;
00246 
00247 namespace internal {
00248 
00249 template<typename _MatrixType, int _UpLo> struct traits<SimplicialLLT<_MatrixType,_UpLo> >
00250 {
00251   typedef _MatrixType MatrixType;
00252   enum { UpLo = _UpLo };
00253   typedef typename MatrixType::Scalar                         Scalar;
00254   typedef typename MatrixType::Index                          Index;
00255   typedef SparseMatrix<Scalar, ColMajor, Index>               CholMatrixType;
00256   typedef SparseTriangularView<CholMatrixType, Eigen::Lower>  MatrixL;
00257   typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper>   MatrixU;
00258   static inline MatrixL getL(const MatrixType& m) { return m; }
00259   static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00260 };
00261 
00262 template<typename _MatrixType,int _UpLo> struct traits<SimplicialLDLT<_MatrixType,_UpLo> >
00263 {
00264   typedef _MatrixType MatrixType;
00265   enum { UpLo = _UpLo };
00266   typedef typename MatrixType::Scalar                             Scalar;
00267   typedef typename MatrixType::Index                              Index;
00268   typedef SparseMatrix<Scalar, ColMajor, Index>                   CholMatrixType;
00269   typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower>  MatrixL;
00270   typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
00271   static inline MatrixL getL(const MatrixType& m) { return m; }
00272   static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
00273 };
00274 
00275 template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_MatrixType,_UpLo> >
00276 {
00277   typedef _MatrixType MatrixType;
00278   enum { UpLo = _UpLo };
00279 };
00280 
00281 }
00282 
00300 template<typename _MatrixType, int _UpLo>
00301     class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo> >
00302 {
00303 public:
00304     typedef _MatrixType MatrixType;
00305     enum { UpLo = _UpLo };
00306     typedef SimplicialCholeskyBase<SimplicialLLT> Base;
00307     typedef typename MatrixType::Scalar Scalar;
00308     typedef typename MatrixType::RealScalar RealScalar;
00309     typedef typename MatrixType::Index Index;
00310     typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
00311     typedef Matrix<Scalar,Dynamic,1> VectorType;
00312     typedef internal::traits<SimplicialLLT> Traits;
00313     typedef typename Traits::MatrixL  MatrixL;
00314     typedef typename Traits::MatrixU  MatrixU;
00315 public:
00317     SimplicialLLT() : Base() {}
00319     SimplicialLLT(const MatrixType& matrix)
00320         : Base(matrix) {}
00321 
00323     inline const MatrixL matrixL() const {
00324         eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
00325         return Traits::getL(Base::m_matrix);
00326     }
00327 
00329     inline const MatrixU matrixU() const {
00330         eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
00331         return Traits::getU(Base::m_matrix);
00332     }
00333     
00335     SimplicialLLT& compute(const MatrixType& matrix)
00336     {
00337       Base::template compute<false>(matrix);
00338       return *this;
00339     }
00340 
00347     void analyzePattern(const MatrixType& a)
00348     {
00349       Base::analyzePattern(a, false);
00350     }
00351 
00358     void factorize(const MatrixType& a)
00359     {
00360       Base::template factorize<false>(a);
00361     }
00362 
00364     Scalar determinant() const
00365     {
00366       Scalar detL = Base::m_matrix.diagonal().prod();
00367       return numext::abs2(detL);
00368     }
00369 };
00370 
00388 template<typename _MatrixType, int _UpLo>
00389     class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo> >
00390 {
00391 public:
00392     typedef _MatrixType MatrixType;
00393     enum { UpLo = _UpLo };
00394     typedef SimplicialCholeskyBase<SimplicialLDLT> Base;
00395     typedef typename MatrixType::Scalar Scalar;
00396     typedef typename MatrixType::RealScalar RealScalar;
00397     typedef typename MatrixType::Index Index;
00398     typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
00399     typedef Matrix<Scalar,Dynamic,1> VectorType;
00400     typedef internal::traits<SimplicialLDLT> Traits;
00401     typedef typename Traits::MatrixL  MatrixL;
00402     typedef typename Traits::MatrixU  MatrixU;
00403 public:
00405     SimplicialLDLT() : Base() {}
00406 
00408     SimplicialLDLT(const MatrixType& matrix)
00409         : Base(matrix) {}
00410 
00412     inline const VectorType vectorD() const {
00413         eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
00414         return Base::m_diag;
00415     }
00417     inline const MatrixL matrixL() const {
00418         eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
00419         return Traits::getL(Base::m_matrix);
00420     }
00421 
00423     inline const MatrixU matrixU() const {
00424         eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
00425         return Traits::getU(Base::m_matrix);
00426     }
00427 
00429     SimplicialLDLT& compute(const MatrixType& matrix)
00430     {
00431       Base::template compute<true>(matrix);
00432       return *this;
00433     }
00434     
00441     void analyzePattern(const MatrixType& a)
00442     {
00443       Base::analyzePattern(a, true);
00444     }
00445 
00452     void factorize(const MatrixType& a)
00453     {
00454       Base::template factorize<true>(a);
00455     }
00456 
00458     Scalar determinant() const
00459     {
00460       return Base::m_diag.prod();
00461     }
00462 };
00463 
00470 template<typename _MatrixType, int _UpLo>
00471     class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo> >
00472 {
00473 public:
00474     typedef _MatrixType MatrixType;
00475     enum { UpLo = _UpLo };
00476     typedef SimplicialCholeskyBase<SimplicialCholesky> Base;
00477     typedef typename MatrixType::Scalar Scalar;
00478     typedef typename MatrixType::RealScalar RealScalar;
00479     typedef typename MatrixType::Index Index;
00480     typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
00481     typedef Matrix<Scalar,Dynamic,1> VectorType;
00482     typedef internal::traits<SimplicialCholesky> Traits;
00483     typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
00484     typedef internal::traits<SimplicialLLT<MatrixType,UpLo>  > LLTTraits;
00485   public:
00486     SimplicialCholesky() : Base(), m_LDLT(true) {}
00487 
00488     SimplicialCholesky(const MatrixType& matrix)
00489       : Base(), m_LDLT(true)
00490     {
00491       compute(matrix);
00492     }
00493 
00494     SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
00495     {
00496       switch(mode)
00497       {
00498       case SimplicialCholeskyLLT:
00499         m_LDLT = false;
00500         break;
00501       case SimplicialCholeskyLDLT:
00502         m_LDLT = true;
00503         break;
00504       default:
00505         break;
00506       }
00507 
00508       return *this;
00509     }
00510 
00511     inline const VectorType vectorD() const {
00512         eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
00513         return Base::m_diag;
00514     }
00515     inline const CholMatrixType rawMatrix() const {
00516         eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
00517         return Base::m_matrix;
00518     }
00519     
00521     SimplicialCholesky& compute(const MatrixType& matrix)
00522     {
00523       if(m_LDLT)
00524         Base::template compute<true>(matrix);
00525       else
00526         Base::template compute<false>(matrix);
00527       return *this;
00528     }
00529 
00536     void analyzePattern(const MatrixType& a)
00537     {
00538       Base::analyzePattern(a, m_LDLT);
00539     }
00540 
00547     void factorize(const MatrixType& a)
00548     {
00549       if(m_LDLT)
00550         Base::template factorize<true>(a);
00551       else
00552         Base::template factorize<false>(a);
00553     }
00554 
00556     template<typename Rhs,typename Dest>
00557     void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
00558     {
00559       eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00560       eigen_assert(Base::m_matrix.rows()==b.rows());
00561 
00562       if(Base::m_info!=Success)
00563         return;
00564 
00565       if(Base::m_P.size()>0)
00566         dest = Base::m_P * b;
00567       else
00568         dest = b;
00569 
00570       if(Base::m_matrix.nonZeros()>0) // otherwise L==I
00571       {
00572         if(m_LDLT)
00573           LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
00574         else
00575           LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
00576       }
00577 
00578       if(Base::m_diag.size()>0)
00579         dest = Base::m_diag.asDiagonal().inverse() * dest;
00580 
00581       if (Base::m_matrix.nonZeros()>0) // otherwise I==I
00582       {
00583         if(m_LDLT)
00584           LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
00585         else
00586           LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
00587       }
00588 
00589       if(Base::m_P.size()>0)
00590         dest = Base::m_Pinv * dest;
00591     }
00592     
00593     Scalar determinant() const
00594     {
00595       if(m_LDLT)
00596       {
00597         return Base::m_diag.prod();
00598       }
00599       else
00600       {
00601         Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
00602         return numext::abs2(detL);
00603       }
00604     }
00605     
00606   protected:
00607     bool m_LDLT;
00608 };
00609 
00610 template<typename Derived>
00611 void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixType& ap)
00612 {
00613   eigen_assert(a.rows()==a.cols());
00614   const Index size = a.rows();
00615   // TODO allows to configure the permutation
00616   // Note that amd compute the inverse permutation
00617   {
00618     CholMatrixType C;
00619     C = a.template selfadjointView<UpLo>();
00620     // remove diagonal entries:
00621     // seems not to be needed
00622     // C.prune(keep_diag());
00623     internal::minimum_degree_ordering(C, m_Pinv);
00624   }
00625 
00626   if(m_Pinv.size()>0)
00627     m_P = m_Pinv.inverse();
00628   else
00629     m_P.resize(0);
00630 
00631   ap.resize(size,size);
00632   ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
00633 }
00634 
00635 namespace internal {
00636   
00637 template<typename Derived, typename Rhs>
00638 struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
00639   : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
00640 {
00641   typedef SimplicialCholeskyBase<Derived> Dec;
00642   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00643 
00644   template<typename Dest> void evalTo(Dest& dst) const
00645   {
00646     dec().derived()._solve(rhs(),dst);
00647   }
00648 };
00649 
00650 template<typename Derived, typename Rhs>
00651 struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
00652   : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
00653 {
00654   typedef SimplicialCholeskyBase<Derived> Dec;
00655   EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00656 
00657   template<typename Dest> void evalTo(Dest& dst) const
00658   {
00659     this->defaultEvalTo(dst);
00660   }
00661 };
00662 
00663 } // end namespace internal
00664 
00665 } // end namespace Eigen
00666 
00667 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:39:06