Eigen  3.3.3
CholmodSupport.h
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2010 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_CHOLMODSUPPORT_H
00011 #define EIGEN_CHOLMODSUPPORT_H
00012 
00013 namespace Eigen { 
00014 
00015 namespace internal {
00016 
00017 template<typename Scalar> struct cholmod_configure_matrix;
00018 
00019 template<> struct cholmod_configure_matrix<double> {
00020   template<typename CholmodType>
00021   static void run(CholmodType& mat) {
00022     mat.xtype = CHOLMOD_REAL;
00023     mat.dtype = CHOLMOD_DOUBLE;
00024   }
00025 };
00026 
00027 template<> struct cholmod_configure_matrix<std::complex<double> > {
00028   template<typename CholmodType>
00029   static void run(CholmodType& mat) {
00030     mat.xtype = CHOLMOD_COMPLEX;
00031     mat.dtype = CHOLMOD_DOUBLE;
00032   }
00033 };
00034 
00035 // Other scalar types are not yet suppotred by Cholmod
00036 // template<> struct cholmod_configure_matrix<float> {
00037 //   template<typename CholmodType>
00038 //   static void run(CholmodType& mat) {
00039 //     mat.xtype = CHOLMOD_REAL;
00040 //     mat.dtype = CHOLMOD_SINGLE;
00041 //   }
00042 // };
00043 //
00044 // template<> struct cholmod_configure_matrix<std::complex<float> > {
00045 //   template<typename CholmodType>
00046 //   static void run(CholmodType& mat) {
00047 //     mat.xtype = CHOLMOD_COMPLEX;
00048 //     mat.dtype = CHOLMOD_SINGLE;
00049 //   }
00050 // };
00051 
00052 } // namespace internal
00053 
00057 template<typename _Scalar, int _Options, typename _StorageIndex>
00058 cholmod_sparse viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_StorageIndex> > mat)
00059 {
00060   cholmod_sparse res;
00061   res.nzmax   = mat.nonZeros();
00062   res.nrow    = mat.rows();
00063   res.ncol    = mat.cols();
00064   res.p       = mat.outerIndexPtr();
00065   res.i       = mat.innerIndexPtr();
00066   res.x       = mat.valuePtr();
00067   res.z       = 0;
00068   res.sorted  = 1;
00069   if(mat.isCompressed())
00070   {
00071     res.packed  = 1;
00072     res.nz = 0;
00073   }
00074   else
00075   {
00076     res.packed  = 0;
00077     res.nz = mat.innerNonZeroPtr();
00078   }
00079 
00080   res.dtype   = 0;
00081   res.stype   = -1;
00082   
00083   if (internal::is_same<_StorageIndex,int>::value)
00084   {
00085     res.itype = CHOLMOD_INT;
00086   }
00087   else if (internal::is_same<_StorageIndex,long>::value)
00088   {
00089     res.itype = CHOLMOD_LONG;
00090   }
00091   else
00092   {
00093     eigen_assert(false && "Index type not supported yet");
00094   }
00095 
00096   // setup res.xtype
00097   internal::cholmod_configure_matrix<_Scalar>::run(res);
00098   
00099   res.stype = 0;
00100   
00101   return res;
00102 }
00103 
00104 template<typename _Scalar, int _Options, typename _Index>
00105 const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat)
00106 {
00107   cholmod_sparse res = viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_Index> >(mat.const_cast_derived()));
00108   return res;
00109 }
00110 
00111 template<typename _Scalar, int _Options, typename _Index>
00112 const cholmod_sparse viewAsCholmod(const SparseVector<_Scalar,_Options,_Index>& mat)
00113 {
00114   cholmod_sparse res = viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_Index> >(mat.const_cast_derived()));
00115   return res;
00116 }
00117 
00120 template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo>
00121 cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<const SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat)
00122 {
00123   cholmod_sparse res = viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_Index> >(mat.matrix().const_cast_derived()));
00124   
00125   if(UpLo==Upper) res.stype =  1;
00126   if(UpLo==Lower) res.stype = -1;
00127 
00128   return res;
00129 }
00130 
00133 template<typename Derived>
00134 cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat)
00135 {
00136   EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00137   typedef typename Derived::Scalar Scalar;
00138 
00139   cholmod_dense res;
00140   res.nrow   = mat.rows();
00141   res.ncol   = mat.cols();
00142   res.nzmax  = res.nrow * res.ncol;
00143   res.d      = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
00144   res.x      = (void*)(mat.derived().data());
00145   res.z      = 0;
00146 
00147   internal::cholmod_configure_matrix<Scalar>::run(res);
00148 
00149   return res;
00150 }
00151 
00154 template<typename Scalar, int Flags, typename StorageIndex>
00155 MappedSparseMatrix<Scalar,Flags,StorageIndex> viewAsEigen(cholmod_sparse& cm)
00156 {
00157   return MappedSparseMatrix<Scalar,Flags,StorageIndex>
00158          (cm.nrow, cm.ncol, static_cast<StorageIndex*>(cm.p)[cm.ncol],
00159           static_cast<StorageIndex*>(cm.p), static_cast<StorageIndex*>(cm.i),static_cast<Scalar*>(cm.x) );
00160 }
00161 
00162 enum CholmodMode {
00163   CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
00164 };
00165 
00166 
00172 template<typename _MatrixType, int _UpLo, typename Derived>
00173 class CholmodBase : public SparseSolverBase<Derived>
00174 {
00175   protected:
00176     typedef SparseSolverBase<Derived> Base;
00177     using Base::derived;
00178     using Base::m_isInitialized;
00179   public:
00180     typedef _MatrixType MatrixType;
00181     enum { UpLo = _UpLo };
00182     typedef typename MatrixType::Scalar Scalar;
00183     typedef typename MatrixType::RealScalar RealScalar;
00184     typedef MatrixType CholMatrixType;
00185     typedef typename MatrixType::StorageIndex StorageIndex;
00186     enum {
00187       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00188       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00189     };
00190 
00191   public:
00192 
00193     CholmodBase()
00194       : m_cholmodFactor(0), m_info(Success), m_factorizationIsOk(false), m_analysisIsOk(false)
00195     {
00196       EIGEN_STATIC_ASSERT((internal::is_same<double,RealScalar>::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY);
00197       m_shiftOffset[0] = m_shiftOffset[1] = 0.0;
00198       cholmod_start(&m_cholmod);
00199     }
00200 
00201     explicit CholmodBase(const MatrixType& matrix)
00202       : m_cholmodFactor(0), m_info(Success), m_factorizationIsOk(false), m_analysisIsOk(false)
00203     {
00204       EIGEN_STATIC_ASSERT((internal::is_same<double,RealScalar>::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY);
00205       m_shiftOffset[0] = m_shiftOffset[1] = 0.0;
00206       cholmod_start(&m_cholmod);
00207       compute(matrix);
00208     }
00209 
00210     ~CholmodBase()
00211     {
00212       if(m_cholmodFactor)
00213         cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
00214       cholmod_finish(&m_cholmod);
00215     }
00216     
00217     inline StorageIndex cols() const { return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
00218     inline StorageIndex rows() const { return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
00219     
00225     ComputationInfo info() const
00226     {
00227       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00228       return m_info;
00229     }
00230 
00232     Derived& compute(const MatrixType& matrix)
00233     {
00234       analyzePattern(matrix);
00235       factorize(matrix);
00236       return derived();
00237     }
00238     
00245     void analyzePattern(const MatrixType& matrix)
00246     {
00247       if(m_cholmodFactor)
00248       {
00249         cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
00250         m_cholmodFactor = 0;
00251       }
00252       cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
00253       m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
00254       
00255       this->m_isInitialized = true;
00256       this->m_info = Success;
00257       m_analysisIsOk = true;
00258       m_factorizationIsOk = false;
00259     }
00260     
00267     void factorize(const MatrixType& matrix)
00268     {
00269       eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00270       cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
00271       cholmod_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod);
00272 
00273       // If the factorization failed, minor is the column at which it did. On success minor == n.
00274       this->m_info = (m_cholmodFactor->minor == m_cholmodFactor->n ? Success : NumericalIssue);
00275       m_factorizationIsOk = true;
00276     }
00277     
00280     cholmod_common& cholmod() { return m_cholmod; }
00281     
00282     #ifndef EIGEN_PARSED_BY_DOXYGEN
00283 
00284     template<typename Rhs,typename Dest>
00285     void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
00286     {
00287       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00288       const Index size = m_cholmodFactor->n;
00289       EIGEN_UNUSED_VARIABLE(size);
00290       eigen_assert(size==b.rows());
00291       
00292       // Cholmod needs column-major stoarge without inner-stride, which corresponds to the default behavior of Ref.
00293       Ref<const Matrix<typename Rhs::Scalar,Dynamic,Dynamic,ColMajor> > b_ref(b.derived());
00294 
00295       cholmod_dense b_cd = viewAsCholmod(b_ref);
00296       cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
00297       if(!x_cd)
00298       {
00299         this->m_info = NumericalIssue;
00300         return;
00301       }
00302       // TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
00303       dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
00304       cholmod_free_dense(&x_cd, &m_cholmod);
00305     }
00306     
00308     template<typename RhsDerived, typename DestDerived>
00309     void _solve_impl(const SparseMatrixBase<RhsDerived> &b, SparseMatrixBase<DestDerived> &dest) const
00310     {
00311       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00312       const Index size = m_cholmodFactor->n;
00313       EIGEN_UNUSED_VARIABLE(size);
00314       eigen_assert(size==b.rows());
00315 
00316       // note: cs stands for Cholmod Sparse
00317       Ref<SparseMatrix<typename RhsDerived::Scalar,ColMajor,typename RhsDerived::StorageIndex> > b_ref(b.const_cast_derived());
00318       cholmod_sparse b_cs = viewAsCholmod(b_ref);
00319       cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
00320       if(!x_cs)
00321       {
00322         this->m_info = NumericalIssue;
00323         return;
00324       }
00325       // TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
00326       dest.derived() = viewAsEigen<typename DestDerived::Scalar,ColMajor,typename DestDerived::StorageIndex>(*x_cs);
00327       cholmod_free_sparse(&x_cs, &m_cholmod);
00328     }
00329     #endif // EIGEN_PARSED_BY_DOXYGEN
00330     
00331     
00341     Derived& setShift(const RealScalar& offset)
00342     {
00343       m_shiftOffset[0] = double(offset);
00344       return derived();
00345     }
00346     
00348     Scalar determinant() const
00349     {
00350       using std::exp;
00351       return exp(logDeterminant());
00352     }
00353 
00355     Scalar logDeterminant() const
00356     {
00357       using std::log;
00358       using numext::real;
00359       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00360 
00361       RealScalar logDet = 0;
00362       Scalar *x = static_cast<Scalar*>(m_cholmodFactor->x);
00363       if (m_cholmodFactor->is_super)
00364       {
00365         // Supernodal factorization stored as a packed list of dense column-major blocs,
00366         // as described by the following structure:
00367 
00368         // super[k] == index of the first column of the j-th super node
00369         StorageIndex *super = static_cast<StorageIndex*>(m_cholmodFactor->super);
00370         // pi[k] == offset to the description of row indices
00371         StorageIndex *pi = static_cast<StorageIndex*>(m_cholmodFactor->pi);
00372         // px[k] == offset to the respective dense block
00373         StorageIndex *px = static_cast<StorageIndex*>(m_cholmodFactor->px);
00374 
00375         Index nb_super_nodes = m_cholmodFactor->nsuper;
00376         for (Index k=0; k < nb_super_nodes; ++k)
00377         {
00378           StorageIndex ncols = super[k + 1] - super[k];
00379           StorageIndex nrows = pi[k + 1] - pi[k];
00380 
00381           Map<const Array<Scalar,1,Dynamic>, 0, InnerStride<> > sk(x + px[k], ncols, InnerStride<>(nrows+1));
00382           logDet += sk.real().log().sum();
00383         }
00384       }
00385       else
00386       {
00387         // Simplicial factorization stored as standard CSC matrix.
00388         StorageIndex *p = static_cast<StorageIndex*>(m_cholmodFactor->p);
00389         Index size = m_cholmodFactor->n;
00390         for (Index k=0; k<size; ++k)
00391           logDet += log(real( x[p[k]] ));
00392       }
00393       if (m_cholmodFactor->is_ll)
00394         logDet *= 2.0;
00395       return logDet;
00396     };
00397 
00398     template<typename Stream>
00399     void dumpMemory(Stream& /*s*/)
00400     {}
00401     
00402   protected:
00403     mutable cholmod_common m_cholmod;
00404     cholmod_factor* m_cholmodFactor;
00405     double m_shiftOffset[2];
00406     mutable ComputationInfo m_info;
00407     int m_factorizationIsOk;
00408     int m_analysisIsOk;
00409 };
00410 
00433 template<typename _MatrixType, int _UpLo = Lower>
00434 class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> >
00435 {
00436     typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT> Base;
00437     using Base::m_cholmod;
00438     
00439   public:
00440     
00441     typedef _MatrixType MatrixType;
00442     
00443     CholmodSimplicialLLT() : Base() { init(); }
00444 
00445     CholmodSimplicialLLT(const MatrixType& matrix) : Base()
00446     {
00447       init();
00448       this->compute(matrix);
00449     }
00450 
00451     ~CholmodSimplicialLLT() {}
00452   protected:
00453     void init()
00454     {
00455       m_cholmod.final_asis = 0;
00456       m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00457       m_cholmod.final_ll = 1;
00458     }
00459 };
00460 
00461 
00484 template<typename _MatrixType, int _UpLo = Lower>
00485 class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> >
00486 {
00487     typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT> Base;
00488     using Base::m_cholmod;
00489     
00490   public:
00491     
00492     typedef _MatrixType MatrixType;
00493     
00494     CholmodSimplicialLDLT() : Base() { init(); }
00495 
00496     CholmodSimplicialLDLT(const MatrixType& matrix) : Base()
00497     {
00498       init();
00499       this->compute(matrix);
00500     }
00501 
00502     ~CholmodSimplicialLDLT() {}
00503   protected:
00504     void init()
00505     {
00506       m_cholmod.final_asis = 1;
00507       m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00508     }
00509 };
00510 
00533 template<typename _MatrixType, int _UpLo = Lower>
00534 class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> >
00535 {
00536     typedef CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT> Base;
00537     using Base::m_cholmod;
00538     
00539   public:
00540     
00541     typedef _MatrixType MatrixType;
00542     
00543     CholmodSupernodalLLT() : Base() { init(); }
00544 
00545     CholmodSupernodalLLT(const MatrixType& matrix) : Base()
00546     {
00547       init();
00548       this->compute(matrix);
00549     }
00550 
00551     ~CholmodSupernodalLLT() {}
00552   protected:
00553     void init()
00554     {
00555       m_cholmod.final_asis = 1;
00556       m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
00557     }
00558 };
00559 
00584 template<typename _MatrixType, int _UpLo = Lower>
00585 class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> >
00586 {
00587     typedef CholmodBase<_MatrixType, _UpLo, CholmodDecomposition> Base;
00588     using Base::m_cholmod;
00589     
00590   public:
00591     
00592     typedef _MatrixType MatrixType;
00593     
00594     CholmodDecomposition() : Base() { init(); }
00595 
00596     CholmodDecomposition(const MatrixType& matrix) : Base()
00597     {
00598       init();
00599       this->compute(matrix);
00600     }
00601 
00602     ~CholmodDecomposition() {}
00603     
00604     void setMode(CholmodMode mode)
00605     {
00606       switch(mode)
00607       {
00608         case CholmodAuto:
00609           m_cholmod.final_asis = 1;
00610           m_cholmod.supernodal = CHOLMOD_AUTO;
00611           break;
00612         case CholmodSimplicialLLt:
00613           m_cholmod.final_asis = 0;
00614           m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00615           m_cholmod.final_ll = 1;
00616           break;
00617         case CholmodSupernodalLLt:
00618           m_cholmod.final_asis = 1;
00619           m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
00620           break;
00621         case CholmodLDLt:
00622           m_cholmod.final_asis = 1;
00623           m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00624           break;
00625         default:
00626           break;
00627       }
00628     }
00629   protected:
00630     void init()
00631     {
00632       m_cholmod.final_asis = 1;
00633       m_cholmod.supernodal = CHOLMOD_AUTO;
00634     }
00635 };
00636 
00637 } // end namespace Eigen
00638 
00639 #endif // EIGEN_CHOLMODSUPPORT_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends