Eigen  3.3.3
SimplicialCholesky.h
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 = &ap;
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
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends