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 // 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, typename CholmodType>
00018 void cholmod_configure_matrix(CholmodType& mat)
00019 {
00020   if (internal::is_same<Scalar,float>::value)
00021   {
00022     mat.xtype = CHOLMOD_REAL;
00023     mat.dtype = CHOLMOD_SINGLE;
00024   }
00025   else if (internal::is_same<Scalar,double>::value)
00026   {
00027     mat.xtype = CHOLMOD_REAL;
00028     mat.dtype = CHOLMOD_DOUBLE;
00029   }
00030   else if (internal::is_same<Scalar,std::complex<float> >::value)
00031   {
00032     mat.xtype = CHOLMOD_COMPLEX;
00033     mat.dtype = CHOLMOD_SINGLE;
00034   }
00035   else if (internal::is_same<Scalar,std::complex<double> >::value)
00036   {
00037     mat.xtype = CHOLMOD_COMPLEX;
00038     mat.dtype = CHOLMOD_DOUBLE;
00039   }
00040   else
00041   {
00042     eigen_assert(false && "Scalar type not supported by CHOLMOD");
00043   }
00044 }
00045 
00046 } // namespace internal
00047 
00051 template<typename _Scalar, int _Options, typename _Index>
00052 cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_Index>& mat)
00053 {
00054   typedef SparseMatrix<_Scalar,_Options,_Index> MatrixType;
00055   cholmod_sparse res;
00056   res.nzmax   = mat.nonZeros();
00057   res.nrow    = mat.rows();;
00058   res.ncol    = mat.cols();
00059   res.p       = mat.outerIndexPtr();
00060   res.i       = mat.innerIndexPtr();
00061   res.x       = mat.valuePtr();
00062   res.sorted  = 1;
00063   if(mat.isCompressed())
00064   {
00065     res.packed  = 1;
00066   }
00067   else
00068   {
00069     res.packed  = 0;
00070     res.nz = mat.innerNonZeroPtr();
00071   }
00072 
00073   res.dtype   = 0;
00074   res.stype   = -1;
00075   
00076   if (internal::is_same<_Index,int>::value)
00077   {
00078     res.itype = CHOLMOD_INT;
00079   }
00080   else
00081   {
00082     eigen_assert(false && "Index type different than int is not supported yet");
00083   }
00084 
00085   // setup res.xtype
00086   internal::cholmod_configure_matrix<_Scalar>(res);
00087   
00088   res.stype = 0;
00089   
00090   return res;
00091 }
00092 
00093 template<typename _Scalar, int _Options, typename _Index>
00094 const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat)
00095 {
00096   cholmod_sparse res = viewAsCholmod(mat.const_cast_derived());
00097   return res;
00098 }
00099 
00102 template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo>
00103 cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat)
00104 {
00105   cholmod_sparse res = viewAsCholmod(mat.matrix().const_cast_derived());
00106   
00107   if(UpLo==Upper) res.stype =  1;
00108   if(UpLo==Lower) res.stype = -1;
00109 
00110   return res;
00111 }
00112 
00115 template<typename Derived>
00116 cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat)
00117 {
00118   EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00119   typedef typename Derived::Scalar Scalar;
00120 
00121   cholmod_dense res;
00122   res.nrow   = mat.rows();
00123   res.ncol   = mat.cols();
00124   res.nzmax  = res.nrow * res.ncol;
00125   res.d      = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
00126   res.x      = mat.derived().data();
00127   res.z      = 0;
00128 
00129   internal::cholmod_configure_matrix<Scalar>(res);
00130 
00131   return res;
00132 }
00133 
00136 template<typename Scalar, int Flags, typename Index>
00137 MappedSparseMatrix<Scalar,Flags,Index> viewAsEigen(cholmod_sparse& cm)
00138 {
00139   return MappedSparseMatrix<Scalar,Flags,Index>
00140          (cm.nrow, cm.ncol, reinterpret_cast<Index*>(cm.p)[cm.ncol],
00141           reinterpret_cast<Index*>(cm.p), reinterpret_cast<Index*>(cm.i),reinterpret_cast<Scalar*>(cm.x) );
00142 }
00143 
00144 enum CholmodMode {
00145   CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
00146 };
00147 
00148 
00154 template<typename _MatrixType, int _UpLo, typename Derived>
00155 class CholmodBase : internal::noncopyable
00156 {
00157   public:
00158     typedef _MatrixType MatrixType;
00159     enum { UpLo = _UpLo };
00160     typedef typename MatrixType::Scalar Scalar;
00161     typedef typename MatrixType::RealScalar RealScalar;
00162     typedef MatrixType CholMatrixType;
00163     typedef typename MatrixType::Index Index;
00164 
00165   public:
00166 
00167     CholmodBase()
00168       : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
00169     {
00170       cholmod_start(&m_cholmod);
00171     }
00172 
00173     CholmodBase(const MatrixType& matrix)
00174       : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
00175     {
00176       cholmod_start(&m_cholmod);
00177       compute(matrix);
00178     }
00179 
00180     ~CholmodBase()
00181     {
00182       if(m_cholmodFactor)
00183         cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
00184       cholmod_finish(&m_cholmod);
00185     }
00186     
00187     inline Index cols() const { return m_cholmodFactor->n; }
00188     inline Index rows() const { return m_cholmodFactor->n; }
00189     
00190     Derived& derived() { return *static_cast<Derived*>(this); }
00191     const Derived& derived() const { return *static_cast<const Derived*>(this); }
00192     
00198     ComputationInfo info() const
00199     {
00200       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00201       return m_info;
00202     }
00203 
00205     Derived& compute(const MatrixType& matrix)
00206     {
00207       analyzePattern(matrix);
00208       factorize(matrix);
00209       return derived();
00210     }
00211     
00216     template<typename Rhs>
00217     inline const internal::solve_retval<CholmodBase, Rhs>
00218     solve(const MatrixBase<Rhs>& b) const
00219     {
00220       eigen_assert(m_isInitialized && "LLT is not initialized.");
00221       eigen_assert(rows()==b.rows()
00222                 && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
00223       return internal::solve_retval<CholmodBase, Rhs>(*this, b.derived());
00224     }
00225     
00230     template<typename Rhs>
00231     inline const internal::sparse_solve_retval<CholmodBase, Rhs>
00232     solve(const SparseMatrixBase<Rhs>& b) const
00233     {
00234       eigen_assert(m_isInitialized && "LLT is not initialized.");
00235       eigen_assert(rows()==b.rows()
00236                 && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
00237       return internal::sparse_solve_retval<CholmodBase, Rhs>(*this, b.derived());
00238     }
00239     
00246     void analyzePattern(const MatrixType& matrix)
00247     {
00248       if(m_cholmodFactor)
00249       {
00250         cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
00251         m_cholmodFactor = 0;
00252       }
00253       cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
00254       m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
00255       
00256       this->m_isInitialized = true;
00257       this->m_info = Success;
00258       m_analysisIsOk = true;
00259       m_factorizationIsOk = false;
00260     }
00261     
00268     void factorize(const MatrixType& matrix)
00269     {
00270       eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00271       cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
00272       cholmod_factorize(&A, m_cholmodFactor, &m_cholmod);
00273       
00274       this->m_info = Success;
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(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_assert(size==b.rows());
00290 
00291       // note: cd stands for Cholmod Dense
00292       cholmod_dense b_cd = viewAsCholmod(b.const_cast_derived());
00293       cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
00294       if(!x_cd)
00295       {
00296         this->m_info = NumericalIssue;
00297       }
00298       // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.)
00299       dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
00300       cholmod_free_dense(&x_cd, &m_cholmod);
00301     }
00302     
00304     template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex>
00305     void _solve(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
00306     {
00307       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00308       const Index size = m_cholmodFactor->n;
00309       eigen_assert(size==b.rows());
00310 
00311       // note: cs stands for Cholmod Sparse
00312       cholmod_sparse b_cs = viewAsCholmod(b);
00313       cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
00314       if(!x_cs)
00315       {
00316         this->m_info = NumericalIssue;
00317       }
00318       // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.)
00319       dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
00320       cholmod_free_sparse(&x_cs, &m_cholmod);
00321     }
00322     #endif // EIGEN_PARSED_BY_DOXYGEN
00323     
00324     template<typename Stream>
00325     void dumpMemory(Stream& s)
00326     {}
00327     
00328   protected:
00329     mutable cholmod_common m_cholmod;
00330     cholmod_factor* m_cholmodFactor;
00331     mutable ComputationInfo m_info;
00332     bool m_isInitialized;
00333     int m_factorizationIsOk;
00334     int m_analysisIsOk;
00335 };
00336 
00355 template<typename _MatrixType, int _UpLo = Lower>
00356 class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> >
00357 {
00358     typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT> Base;
00359     using Base::m_cholmod;
00360     
00361   public:
00362     
00363     typedef _MatrixType MatrixType;
00364     
00365     CholmodSimplicialLLT() : Base() { init(); }
00366 
00367     CholmodSimplicialLLT(const MatrixType& matrix) : Base()
00368     {
00369       init();
00370       compute(matrix);
00371     }
00372 
00373     ~CholmodSimplicialLLT() {}
00374   protected:
00375     void init()
00376     {
00377       m_cholmod.final_asis = 0;
00378       m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00379       m_cholmod.final_ll = 1;
00380     }
00381 };
00382 
00383 
00402 template<typename _MatrixType, int _UpLo = Lower>
00403 class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> >
00404 {
00405     typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT> Base;
00406     using Base::m_cholmod;
00407     
00408   public:
00409     
00410     typedef _MatrixType MatrixType;
00411     
00412     CholmodSimplicialLDLT() : Base() { init(); }
00413 
00414     CholmodSimplicialLDLT(const MatrixType& matrix) : Base()
00415     {
00416       init();
00417       compute(matrix);
00418     }
00419 
00420     ~CholmodSimplicialLDLT() {}
00421   protected:
00422     void init()
00423     {
00424       m_cholmod.final_asis = 1;
00425       m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00426     }
00427 };
00428 
00447 template<typename _MatrixType, int _UpLo = Lower>
00448 class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> >
00449 {
00450     typedef CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT> Base;
00451     using Base::m_cholmod;
00452     
00453   public:
00454     
00455     typedef _MatrixType MatrixType;
00456     
00457     CholmodSupernodalLLT() : Base() { init(); }
00458 
00459     CholmodSupernodalLLT(const MatrixType& matrix) : Base()
00460     {
00461       init();
00462       compute(matrix);
00463     }
00464 
00465     ~CholmodSupernodalLLT() {}
00466   protected:
00467     void init()
00468     {
00469       m_cholmod.final_asis = 1;
00470       m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
00471     }
00472 };
00473 
00494 template<typename _MatrixType, int _UpLo = Lower>
00495 class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> >
00496 {
00497     typedef CholmodBase<_MatrixType, _UpLo, CholmodDecomposition> Base;
00498     using Base::m_cholmod;
00499     
00500   public:
00501     
00502     typedef _MatrixType MatrixType;
00503     
00504     CholmodDecomposition() : Base() { init(); }
00505 
00506     CholmodDecomposition(const MatrixType& matrix) : Base()
00507     {
00508       init();
00509       compute(matrix);
00510     }
00511 
00512     ~CholmodDecomposition() {}
00513     
00514     void setMode(CholmodMode mode)
00515     {
00516       switch(mode)
00517       {
00518         case CholmodAuto:
00519           m_cholmod.final_asis = 1;
00520           m_cholmod.supernodal = CHOLMOD_AUTO;
00521           break;
00522         case CholmodSimplicialLLt:
00523           m_cholmod.final_asis = 0;
00524           m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00525           m_cholmod.final_ll = 1;
00526           break;
00527         case CholmodSupernodalLLt:
00528           m_cholmod.final_asis = 1;
00529           m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
00530           break;
00531         case CholmodLDLt:
00532           m_cholmod.final_asis = 1;
00533           m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00534           break;
00535         default:
00536           break;
00537       }
00538     }
00539   protected:
00540     void init()
00541     {
00542       m_cholmod.final_asis = 1;
00543       m_cholmod.supernodal = CHOLMOD_AUTO;
00544     }
00545 };
00546 
00547 namespace internal {
00548   
00549 template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
00550 struct solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
00551   : solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
00552 {
00553   typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
00554   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00555 
00556   template<typename Dest> void evalTo(Dest& dst) const
00557   {
00558     dec()._solve(rhs(),dst);
00559   }
00560 };
00561 
00562 template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
00563 struct sparse_solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
00564   : sparse_solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
00565 {
00566   typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
00567   EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00568 
00569   template<typename Dest> void evalTo(Dest& dst) const
00570   {
00571     dec()._solve(rhs(),dst);
00572   }
00573 };
00574 
00575 } // end namespace internal
00576 
00577 } // end namespace Eigen
00578 
00579 #endif // EIGEN_CHOLMODSUPPORT_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:24:13