CholmodSupport.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-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_CHOLMODSUPPORT_H
00026 #define EIGEN_CHOLMODSUPPORT_H
00027 
00028 namespace Eigen { 
00029 
00030 namespace internal {
00031 
00032 template<typename Scalar, typename CholmodType>
00033 void cholmod_configure_matrix(CholmodType& mat)
00034 {
00035   if (internal::is_same<Scalar,float>::value)
00036   {
00037     mat.xtype = CHOLMOD_REAL;
00038     mat.dtype = CHOLMOD_SINGLE;
00039   }
00040   else if (internal::is_same<Scalar,double>::value)
00041   {
00042     mat.xtype = CHOLMOD_REAL;
00043     mat.dtype = CHOLMOD_DOUBLE;
00044   }
00045   else if (internal::is_same<Scalar,std::complex<float> >::value)
00046   {
00047     mat.xtype = CHOLMOD_COMPLEX;
00048     mat.dtype = CHOLMOD_SINGLE;
00049   }
00050   else if (internal::is_same<Scalar,std::complex<double> >::value)
00051   {
00052     mat.xtype = CHOLMOD_COMPLEX;
00053     mat.dtype = CHOLMOD_DOUBLE;
00054   }
00055   else
00056   {
00057     eigen_assert(false && "Scalar type not supported by CHOLMOD");
00058   }
00059 }
00060 
00061 } // namespace internal
00062 
00066 template<typename _Scalar, int _Options, typename _Index>
00067 cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_Index>& mat)
00068 {
00069   typedef SparseMatrix<_Scalar,_Options,_Index> MatrixType;
00070   cholmod_sparse res;
00071   res.nzmax   = mat.nonZeros();
00072   res.nrow    = mat.rows();;
00073   res.ncol    = mat.cols();
00074   res.p       = mat.outerIndexPtr();
00075   res.i       = mat.innerIndexPtr();
00076   res.x       = mat.valuePtr();
00077   res.sorted  = 1;
00078   if(mat.isCompressed())
00079   {
00080     res.packed  = 1;
00081   }
00082   else
00083   {
00084     res.packed  = 0;
00085     res.nz = mat.innerNonZeroPtr();
00086   }
00087 
00088   res.dtype   = 0;
00089   res.stype   = -1;
00090   
00091   if (internal::is_same<_Index,int>::value)
00092   {
00093     res.itype = CHOLMOD_INT;
00094   }
00095   else
00096   {
00097     eigen_assert(false && "Index type different than int is not supported yet");
00098   }
00099 
00100   // setup res.xtype
00101   internal::cholmod_configure_matrix<_Scalar>(res);
00102   
00103   res.stype = 0;
00104   
00105   return res;
00106 }
00107 
00108 template<typename _Scalar, int _Options, typename _Index>
00109 const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat)
00110 {
00111   cholmod_sparse res = viewAsCholmod(mat.const_cast_derived());
00112   return res;
00113 }
00114 
00117 template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo>
00118 cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat)
00119 {
00120   cholmod_sparse res = viewAsCholmod(mat.matrix().const_cast_derived());
00121   
00122   if(UpLo==Upper) res.stype =  1;
00123   if(UpLo==Lower) res.stype = -1;
00124 
00125   return res;
00126 }
00127 
00130 template<typename Derived>
00131 cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat)
00132 {
00133   EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00134   typedef typename Derived::Scalar Scalar;
00135 
00136   cholmod_dense res;
00137   res.nrow   = mat.rows();
00138   res.ncol   = mat.cols();
00139   res.nzmax  = res.nrow * res.ncol;
00140   res.d      = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
00141   res.x      = mat.derived().data();
00142   res.z      = 0;
00143 
00144   internal::cholmod_configure_matrix<Scalar>(res);
00145 
00146   return res;
00147 }
00148 
00151 template<typename Scalar, int Flags, typename Index>
00152 MappedSparseMatrix<Scalar,Flags,Index> viewAsEigen(cholmod_sparse& cm)
00153 {
00154   return MappedSparseMatrix<Scalar,Flags,Index>
00155          (cm.nrow, cm.ncol, reinterpret_cast<Index*>(cm.p)[cm.ncol],
00156           reinterpret_cast<Index*>(cm.p), reinterpret_cast<Index*>(cm.i),reinterpret_cast<Scalar*>(cm.x) );
00157 }
00158 
00159 enum CholmodMode {
00160   CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
00161 };
00162 
00163 
00169 template<typename _MatrixType, int _UpLo, typename Derived>
00170 class CholmodBase : internal::noncopyable
00171 {
00172   public:
00173     typedef _MatrixType MatrixType;
00174     enum { UpLo = _UpLo };
00175     typedef typename MatrixType::Scalar Scalar;
00176     typedef typename MatrixType::RealScalar RealScalar;
00177     typedef MatrixType CholMatrixType;
00178     typedef typename MatrixType::Index Index;
00179 
00180   public:
00181 
00182     CholmodBase()
00183       : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
00184     {
00185       cholmod_start(&m_cholmod);
00186     }
00187 
00188     CholmodBase(const MatrixType& matrix)
00189       : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
00190     {
00191       cholmod_start(&m_cholmod);
00192       compute(matrix);
00193     }
00194 
00195     ~CholmodBase()
00196     {
00197       if(m_cholmodFactor)
00198         cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
00199       cholmod_finish(&m_cholmod);
00200     }
00201     
00202     inline Index cols() const { return m_cholmodFactor->n; }
00203     inline Index rows() const { return m_cholmodFactor->n; }
00204     
00205     Derived& derived() { return *static_cast<Derived*>(this); }
00206     const Derived& derived() const { return *static_cast<const Derived*>(this); }
00207     
00213     ComputationInfo info() const
00214     {
00215       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00216       return m_info;
00217     }
00218 
00220     Derived& compute(const MatrixType& matrix)
00221     {
00222       analyzePattern(matrix);
00223       factorize(matrix);
00224       return derived();
00225     }
00226     
00231     template<typename Rhs>
00232     inline const internal::solve_retval<CholmodBase, Rhs>
00233     solve(const MatrixBase<Rhs>& b) const
00234     {
00235       eigen_assert(m_isInitialized && "LLT is not initialized.");
00236       eigen_assert(rows()==b.rows()
00237                 && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
00238       return internal::solve_retval<CholmodBase, Rhs>(*this, b.derived());
00239     }
00240     
00245     template<typename Rhs>
00246     inline const internal::sparse_solve_retval<CholmodBase, Rhs>
00247     solve(const SparseMatrixBase<Rhs>& b) const
00248     {
00249       eigen_assert(m_isInitialized && "LLT is not initialized.");
00250       eigen_assert(rows()==b.rows()
00251                 && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
00252       return internal::sparse_solve_retval<CholmodBase, Rhs>(*this, b.derived());
00253     }
00254     
00261     void analyzePattern(const MatrixType& matrix)
00262     {
00263       if(m_cholmodFactor)
00264       {
00265         cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
00266         m_cholmodFactor = 0;
00267       }
00268       cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
00269       m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
00270       
00271       this->m_isInitialized = true;
00272       this->m_info = Success;
00273       m_analysisIsOk = true;
00274       m_factorizationIsOk = false;
00275     }
00276     
00283     void factorize(const MatrixType& matrix)
00284     {
00285       eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00286       cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
00287       cholmod_factorize(&A, m_cholmodFactor, &m_cholmod);
00288       
00289       this->m_info = Success;
00290       m_factorizationIsOk = true;
00291     }
00292     
00295     cholmod_common& cholmod() { return m_cholmod; }
00296     
00297     #ifndef EIGEN_PARSED_BY_DOXYGEN
00298 
00299     template<typename Rhs,typename Dest>
00300     void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
00301     {
00302       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00303       const Index size = m_cholmodFactor->n;
00304       eigen_assert(size==b.rows());
00305 
00306       // note: cd stands for Cholmod Dense
00307       cholmod_dense b_cd = viewAsCholmod(b.const_cast_derived());
00308       cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
00309       if(!x_cd)
00310       {
00311         this->m_info = NumericalIssue;
00312       }
00313       // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.)
00314       dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
00315       cholmod_free_dense(&x_cd, &m_cholmod);
00316     }
00317     
00319     template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex>
00320     void _solve(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
00321     {
00322       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00323       const Index size = m_cholmodFactor->n;
00324       eigen_assert(size==b.rows());
00325 
00326       // note: cs stands for Cholmod Sparse
00327       cholmod_sparse b_cs = viewAsCholmod(b);
00328       cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
00329       if(!x_cs)
00330       {
00331         this->m_info = NumericalIssue;
00332       }
00333       // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.)
00334       dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
00335       cholmod_free_sparse(&x_cs, &m_cholmod);
00336     }
00337     #endif // EIGEN_PARSED_BY_DOXYGEN
00338     
00339     template<typename Stream>
00340     void dumpMemory(Stream& s)
00341     {}
00342     
00343   protected:
00344     mutable cholmod_common m_cholmod;
00345     cholmod_factor* m_cholmodFactor;
00346     mutable ComputationInfo m_info;
00347     bool m_isInitialized;
00348     int m_factorizationIsOk;
00349     int m_analysisIsOk;
00350 };
00351 
00370 template<typename _MatrixType, int _UpLo = Lower>
00371 class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> >
00372 {
00373     typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT> Base;
00374     using Base::m_cholmod;
00375     
00376   public:
00377     
00378     typedef _MatrixType MatrixType;
00379     
00380     CholmodSimplicialLLT() : Base() { init(); }
00381 
00382     CholmodSimplicialLLT(const MatrixType& matrix) : Base()
00383     {
00384       init();
00385       compute(matrix);
00386     }
00387 
00388     ~CholmodSimplicialLLT() {}
00389   protected:
00390     void init()
00391     {
00392       m_cholmod.final_asis = 0;
00393       m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00394       m_cholmod.final_ll = 1;
00395     }
00396 };
00397 
00398 
00417 template<typename _MatrixType, int _UpLo = Lower>
00418 class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> >
00419 {
00420     typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT> Base;
00421     using Base::m_cholmod;
00422     
00423   public:
00424     
00425     typedef _MatrixType MatrixType;
00426     
00427     CholmodSimplicialLDLT() : Base() { init(); }
00428 
00429     CholmodSimplicialLDLT(const MatrixType& matrix) : Base()
00430     {
00431       init();
00432       compute(matrix);
00433     }
00434 
00435     ~CholmodSimplicialLDLT() {}
00436   protected:
00437     void init()
00438     {
00439       m_cholmod.final_asis = 1;
00440       m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00441     }
00442 };
00443 
00462 template<typename _MatrixType, int _UpLo = Lower>
00463 class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> >
00464 {
00465     typedef CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT> Base;
00466     using Base::m_cholmod;
00467     
00468   public:
00469     
00470     typedef _MatrixType MatrixType;
00471     
00472     CholmodSupernodalLLT() : Base() { init(); }
00473 
00474     CholmodSupernodalLLT(const MatrixType& matrix) : Base()
00475     {
00476       init();
00477       compute(matrix);
00478     }
00479 
00480     ~CholmodSupernodalLLT() {}
00481   protected:
00482     void init()
00483     {
00484       m_cholmod.final_asis = 1;
00485       m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
00486     }
00487 };
00488 
00509 template<typename _MatrixType, int _UpLo = Lower>
00510 class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> >
00511 {
00512     typedef CholmodBase<_MatrixType, _UpLo, CholmodDecomposition> Base;
00513     using Base::m_cholmod;
00514     
00515   public:
00516     
00517     typedef _MatrixType MatrixType;
00518     
00519     CholmodDecomposition() : Base() { init(); }
00520 
00521     CholmodDecomposition(const MatrixType& matrix) : Base()
00522     {
00523       init();
00524       compute(matrix);
00525     }
00526 
00527     ~CholmodDecomposition() {}
00528     
00529     void setMode(CholmodMode mode)
00530     {
00531       switch(mode)
00532       {
00533         case CholmodAuto:
00534           m_cholmod.final_asis = 1;
00535           m_cholmod.supernodal = CHOLMOD_AUTO;
00536           break;
00537         case CholmodSimplicialLLt:
00538           m_cholmod.final_asis = 0;
00539           m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00540           m_cholmod.final_ll = 1;
00541           break;
00542         case CholmodSupernodalLLt:
00543           m_cholmod.final_asis = 1;
00544           m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
00545           break;
00546         case CholmodLDLt:
00547           m_cholmod.final_asis = 1;
00548           m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00549           break;
00550         default:
00551           break;
00552       }
00553     }
00554   protected:
00555     void init()
00556     {
00557       m_cholmod.final_asis = 1;
00558       m_cholmod.supernodal = CHOLMOD_AUTO;
00559     }
00560 };
00561 
00562 namespace internal {
00563   
00564 template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
00565 struct solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
00566   : solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
00567 {
00568   typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
00569   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00570 
00571   template<typename Dest> void evalTo(Dest& dst) const
00572   {
00573     dec()._solve(rhs(),dst);
00574   }
00575 };
00576 
00577 template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
00578 struct sparse_solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
00579   : sparse_solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
00580 {
00581   typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
00582   EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00583 
00584   template<typename Dest> void evalTo(Dest& dst) const
00585   {
00586     dec()._solve(rhs(),dst);
00587   }
00588 };
00589 
00590 } // end namespace internal
00591 
00592 } // end namespace Eigen
00593 
00594 #endif // EIGEN_CHOLMODSUPPORT_H