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 internal {
00029 
00030 template<typename Scalar, typename CholmodType>
00031 void cholmod_configure_matrix(CholmodType& mat)
00032 {
00033   if (internal::is_same<Scalar,float>::value)
00034   {
00035     mat.xtype = CHOLMOD_REAL;
00036     mat.dtype = CHOLMOD_SINGLE;
00037   }
00038   else if (internal::is_same<Scalar,double>::value)
00039   {
00040     mat.xtype = CHOLMOD_REAL;
00041     mat.dtype = CHOLMOD_DOUBLE;
00042   }
00043   else if (internal::is_same<Scalar,std::complex<float> >::value)
00044   {
00045     mat.xtype = CHOLMOD_COMPLEX;
00046     mat.dtype = CHOLMOD_SINGLE;
00047   }
00048   else if (internal::is_same<Scalar,std::complex<double> >::value)
00049   {
00050     mat.xtype = CHOLMOD_COMPLEX;
00051     mat.dtype = CHOLMOD_DOUBLE;
00052   }
00053   else
00054   {
00055     eigen_assert(false && "Scalar type not supported by CHOLMOD");
00056   }
00057 }
00058 
00059 } // namespace internal
00060 
00064 template<typename _Scalar, int _Options, typename _Index>
00065 cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_Index>& mat)
00066 {
00067   typedef SparseMatrix<_Scalar,_Options,_Index> MatrixType;
00068   cholmod_sparse res;
00069   res.nzmax   = mat.nonZeros();
00070   res.nrow    = mat.rows();;
00071   res.ncol    = mat.cols();
00072   res.p       = mat._outerIndexPtr();
00073   res.i       = mat._innerIndexPtr();
00074   res.x       = mat._valuePtr();
00075   res.sorted  = 1;
00076   res.packed  = 1;
00077   res.dtype   = 0;
00078   res.stype   = -1;
00079   
00080   if (internal::is_same<_Index,int>::value)
00081   {
00082     res.itype = CHOLMOD_INT;
00083   }
00084   else
00085   {
00086     eigen_assert(false && "Index type different than int is not supported yet");
00087   }
00088 
00089   // setup res.xtype
00090   internal::cholmod_configure_matrix<_Scalar>(res);
00091   
00092   res.stype = 0;
00093   
00094   return res;
00095 }
00096 
00097 template<typename _Scalar, int _Options, typename _Index>
00098 const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat)
00099 {
00100   cholmod_sparse res = viewAsCholmod(mat.const_cast_derived());
00101   return res;
00102 }
00103 
00106 template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo>
00107 cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat)
00108 {
00109   cholmod_sparse res = viewAsCholmod(mat.matrix().const_cast_derived());
00110   
00111   if(UpLo==Upper) res.stype =  1;
00112   if(UpLo==Lower) res.stype = -1;
00113 
00114   return res;
00115 }
00116 
00119 template<typename Derived>
00120 cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat)
00121 {
00122   EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00123   typedef typename Derived::Scalar Scalar;
00124 
00125   cholmod_dense res;
00126   res.nrow   = mat.rows();
00127   res.ncol   = mat.cols();
00128   res.nzmax  = res.nrow * res.ncol;
00129   res.d      = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
00130   res.x      = mat.derived().data();
00131   res.z      = 0;
00132 
00133   internal::cholmod_configure_matrix<Scalar>(res);
00134 
00135   return res;
00136 }
00137 
00140 template<typename Scalar, int Flags, typename Index>
00141 MappedSparseMatrix<Scalar,Flags,Index> viewAsEigen(cholmod_sparse& cm)
00142 {
00143   return MappedSparseMatrix<Scalar,Flags,Index>
00144          (cm.nrow, cm.ncol, reinterpret_cast<Index*>(cm.p)[cm.ncol],
00145           reinterpret_cast<Index*>(cm.p), reinterpret_cast<Index*>(cm.i),reinterpret_cast<Scalar*>(cm.x) );
00146 }
00147 
00148 enum CholmodMode {
00149   CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
00150 };
00151 
00163 template<typename _MatrixType, int _UpLo = Lower>
00164 class CholmodDecomposition
00165 {
00166   public:
00167     typedef _MatrixType MatrixType;
00168     enum { UpLo = _UpLo };
00169     typedef typename MatrixType::Scalar Scalar;
00170     typedef typename MatrixType::RealScalar RealScalar;
00171     typedef MatrixType CholMatrixType;
00172     typedef typename MatrixType::Index Index;
00173 
00174   public:
00175 
00176     CholmodDecomposition()
00177       : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
00178     {
00179       cholmod_start(&m_cholmod);
00180       setMode(CholmodLDLt);
00181     }
00182 
00183     CholmodDecomposition(const MatrixType& matrix)
00184       : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
00185     {
00186       cholmod_start(&m_cholmod);
00187       compute(matrix);
00188     }
00189 
00190     ~CholmodDecomposition()
00191     {
00192       if(m_cholmodFactor)
00193         cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
00194       cholmod_finish(&m_cholmod);
00195     }
00196     
00197     inline Index cols() const { return m_cholmodFactor->n; }
00198     inline Index rows() const { return m_cholmodFactor->n; }
00199     
00200     void setMode(CholmodMode mode)
00201     {
00202       switch(mode)
00203       {
00204         case CholmodAuto:
00205           m_cholmod.final_asis = 1;
00206           m_cholmod.supernodal = CHOLMOD_AUTO;
00207           break;
00208         case CholmodSimplicialLLt:
00209           m_cholmod.final_asis = 0;
00210           m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00211           m_cholmod.final_ll = 1;
00212           break;
00213         case CholmodSupernodalLLt:
00214           m_cholmod.final_asis = 1;
00215           m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
00216           break;
00217         case CholmodLDLt:
00218           m_cholmod.final_asis = 1;
00219           m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
00220           break;
00221         default:
00222           break;
00223       }
00224     }
00225     
00231     ComputationInfo info() const
00232     {
00233       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00234       return m_info;
00235     }
00236 
00238     void compute(const MatrixType& matrix)
00239     {
00240       analyzePattern(matrix);
00241       factorize(matrix);
00242     }
00243     
00248     template<typename Rhs>
00249     inline const internal::solve_retval<CholmodDecomposition, Rhs>
00250     solve(const MatrixBase<Rhs>& b) const
00251     {
00252       eigen_assert(m_isInitialized && "LLT is not initialized.");
00253       eigen_assert(rows()==b.rows()
00254                 && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
00255       return internal::solve_retval<CholmodDecomposition, Rhs>(*this, b.derived());
00256     }
00257     
00262     template<typename Rhs>
00263     inline const internal::sparse_solve_retval<CholmodDecomposition, Rhs>
00264     solve(const SparseMatrixBase<Rhs>& b) const
00265     {
00266       eigen_assert(m_isInitialized && "LLT is not initialized.");
00267       eigen_assert(rows()==b.rows()
00268                 && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
00269       return internal::sparse_solve_retval<CholmodDecomposition, Rhs>(*this, b.derived());
00270     }
00271     
00278     void analyzePattern(const MatrixType& matrix)
00279     {
00280       if(m_cholmodFactor)
00281       {
00282         cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
00283         m_cholmodFactor = 0;
00284       }
00285       cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
00286       m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
00287       
00288       this->m_isInitialized = true;
00289       this->m_info = Success;
00290       m_analysisIsOk = true;
00291       m_factorizationIsOk = false;
00292     }
00293     
00300     void factorize(const MatrixType& matrix)
00301     {
00302       eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
00303       cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
00304       cholmod_factorize(&A, m_cholmodFactor, &m_cholmod);
00305       
00306       this->m_info = Success;
00307       m_factorizationIsOk = true;
00308     }
00309     
00312     cholmod_common& cholmod() { return m_cholmod; }
00313     
00314     #ifndef EIGEN_PARSED_BY_DOXYGEN
00315 
00316     template<typename Rhs,typename Dest>
00317     void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
00318     {
00319       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00320       const Index size = m_cholmodFactor->n;
00321       eigen_assert(size==b.rows());
00322 
00323       // note: cd stands for Cholmod Dense
00324       cholmod_dense b_cd = viewAsCholmod(b.const_cast_derived());
00325       cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
00326       if(!x_cd)
00327       {
00328         this->m_info = NumericalIssue;
00329       }
00330       // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.)
00331       dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
00332       cholmod_free_dense(&x_cd, &m_cholmod);
00333     }
00334     
00336     template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex>
00337     void _solve(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
00338     {
00339       eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
00340       const Index size = m_cholmodFactor->n;
00341       eigen_assert(size==b.rows());
00342 
00343       // note: cs stands for Cholmod Sparse
00344       cholmod_sparse b_cs = viewAsCholmod(b);
00345       cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
00346       if(!x_cs)
00347       {
00348         this->m_info = NumericalIssue;
00349       }
00350       // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.)
00351       dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
00352       cholmod_free_sparse(&x_cs, &m_cholmod);
00353     }
00354     #endif // EIGEN_PARSED_BY_DOXYGEN
00355     
00356     template<typename Stream>
00357     void dumpMemory(Stream& s)
00358     {}
00359 
00360   protected:
00361     mutable cholmod_common m_cholmod;
00362     cholmod_factor* m_cholmodFactor;
00363     mutable ComputationInfo m_info;
00364     bool m_isInitialized;
00365     int m_factorizationIsOk;
00366     int m_analysisIsOk;
00367 };
00368 
00369 namespace internal {
00370   
00371 template<typename _MatrixType, int _UpLo, typename Rhs>
00372 struct solve_retval<CholmodDecomposition<_MatrixType,_UpLo>, Rhs>
00373   : solve_retval_base<CholmodDecomposition<_MatrixType,_UpLo>, Rhs>
00374 {
00375   typedef CholmodDecomposition<_MatrixType,_UpLo> Dec;
00376   EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00377 
00378   template<typename Dest> void evalTo(Dest& dst) const
00379   {
00380     dec()._solve(rhs(),dst);
00381   }
00382 };
00383 
00384 template<typename _MatrixType, int _UpLo, typename Rhs>
00385 struct sparse_solve_retval<CholmodDecomposition<_MatrixType,_UpLo>, Rhs>
00386   : sparse_solve_retval_base<CholmodDecomposition<_MatrixType,_UpLo>, Rhs>
00387 {
00388   typedef CholmodDecomposition<_MatrixType,_UpLo> Dec;
00389   EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00390 
00391   template<typename Dest> void evalTo(Dest& dst) const
00392   {
00393     dec()._solve(rhs(),dst);
00394   }
00395 };
00396 
00397 }
00398 
00399 #endif // EIGEN_CHOLMODSUPPORT_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:31