$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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 cholmod_sparse res; 00055 res.nzmax = mat.nonZeros(); 00056 res.nrow = mat.rows();; 00057 res.ncol = mat.cols(); 00058 res.p = mat.outerIndexPtr(); 00059 res.i = mat.innerIndexPtr(); 00060 res.x = mat.valuePtr(); 00061 res.z = 0; 00062 res.sorted = 1; 00063 if(mat.isCompressed()) 00064 { 00065 res.packed = 1; 00066 res.nz = 0; 00067 } 00068 else 00069 { 00070 res.packed = 0; 00071 res.nz = mat.innerNonZeroPtr(); 00072 } 00073 00074 res.dtype = 0; 00075 res.stype = -1; 00076 00077 if (internal::is_same<_Index,int>::value) 00078 { 00079 res.itype = CHOLMOD_INT; 00080 } 00081 else if (internal::is_same<_Index,UF_long>::value) 00082 { 00083 res.itype = CHOLMOD_LONG; 00084 } 00085 else 00086 { 00087 eigen_assert(false && "Index type not supported yet"); 00088 } 00089 00090 // setup res.xtype 00091 internal::cholmod_configure_matrix<_Scalar>(res); 00092 00093 res.stype = 0; 00094 00095 return res; 00096 } 00097 00098 template<typename _Scalar, int _Options, typename _Index> 00099 const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat) 00100 { 00101 cholmod_sparse res = viewAsCholmod(mat.const_cast_derived()); 00102 return res; 00103 } 00104 00107 template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo> 00108 cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat) 00109 { 00110 cholmod_sparse res = viewAsCholmod(mat.matrix().const_cast_derived()); 00111 00112 if(UpLo==Upper) res.stype = 1; 00113 if(UpLo==Lower) res.stype = -1; 00114 00115 return res; 00116 } 00117 00120 template<typename Derived> 00121 cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat) 00122 { 00123 EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES); 00124 typedef typename Derived::Scalar Scalar; 00125 00126 cholmod_dense res; 00127 res.nrow = mat.rows(); 00128 res.ncol = mat.cols(); 00129 res.nzmax = res.nrow * res.ncol; 00130 res.d = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride(); 00131 res.x = (void*)(mat.derived().data()); 00132 res.z = 0; 00133 00134 internal::cholmod_configure_matrix<Scalar>(res); 00135 00136 return res; 00137 } 00138 00141 template<typename Scalar, int Flags, typename Index> 00142 MappedSparseMatrix<Scalar,Flags,Index> viewAsEigen(cholmod_sparse& cm) 00143 { 00144 return MappedSparseMatrix<Scalar,Flags,Index> 00145 (cm.nrow, cm.ncol, static_cast<Index*>(cm.p)[cm.ncol], 00146 static_cast<Index*>(cm.p), static_cast<Index*>(cm.i),static_cast<Scalar*>(cm.x) ); 00147 } 00148 00149 enum CholmodMode { 00150 CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt 00151 }; 00152 00153 00159 template<typename _MatrixType, int _UpLo, typename Derived> 00160 class CholmodBase : internal::noncopyable 00161 { 00162 public: 00163 typedef _MatrixType MatrixType; 00164 enum { UpLo = _UpLo }; 00165 typedef typename MatrixType::Scalar Scalar; 00166 typedef typename MatrixType::RealScalar RealScalar; 00167 typedef MatrixType CholMatrixType; 00168 typedef typename MatrixType::Index Index; 00169 00170 public: 00171 00172 CholmodBase() 00173 : m_cholmodFactor(0), m_info(Success), m_isInitialized(false) 00174 { 00175 m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0); 00176 cholmod_start(&m_cholmod); 00177 } 00178 00179 CholmodBase(const MatrixType& matrix) 00180 : m_cholmodFactor(0), m_info(Success), m_isInitialized(false) 00181 { 00182 m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0); 00183 cholmod_start(&m_cholmod); 00184 compute(matrix); 00185 } 00186 00187 ~CholmodBase() 00188 { 00189 if(m_cholmodFactor) 00190 cholmod_free_factor(&m_cholmodFactor, &m_cholmod); 00191 cholmod_finish(&m_cholmod); 00192 } 00193 00194 inline Index cols() const { return m_cholmodFactor->n; } 00195 inline Index rows() const { return m_cholmodFactor->n; } 00196 00197 Derived& derived() { return *static_cast<Derived*>(this); } 00198 const Derived& derived() const { return *static_cast<const Derived*>(this); } 00199 00205 ComputationInfo info() const 00206 { 00207 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 00208 return m_info; 00209 } 00210 00212 Derived& compute(const MatrixType& matrix) 00213 { 00214 analyzePattern(matrix); 00215 factorize(matrix); 00216 return derived(); 00217 } 00218 00223 template<typename Rhs> 00224 inline const internal::solve_retval<CholmodBase, Rhs> 00225 solve(const MatrixBase<Rhs>& b) const 00226 { 00227 eigen_assert(m_isInitialized && "LLT is not initialized."); 00228 eigen_assert(rows()==b.rows() 00229 && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b"); 00230 return internal::solve_retval<CholmodBase, Rhs>(*this, b.derived()); 00231 } 00232 00237 template<typename Rhs> 00238 inline const internal::sparse_solve_retval<CholmodBase, Rhs> 00239 solve(const SparseMatrixBase<Rhs>& b) const 00240 { 00241 eigen_assert(m_isInitialized && "LLT is not initialized."); 00242 eigen_assert(rows()==b.rows() 00243 && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b"); 00244 return internal::sparse_solve_retval<CholmodBase, Rhs>(*this, b.derived()); 00245 } 00246 00253 void analyzePattern(const MatrixType& matrix) 00254 { 00255 if(m_cholmodFactor) 00256 { 00257 cholmod_free_factor(&m_cholmodFactor, &m_cholmod); 00258 m_cholmodFactor = 0; 00259 } 00260 cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>()); 00261 m_cholmodFactor = cholmod_analyze(&A, &m_cholmod); 00262 00263 this->m_isInitialized = true; 00264 this->m_info = Success; 00265 m_analysisIsOk = true; 00266 m_factorizationIsOk = false; 00267 } 00268 00275 void factorize(const MatrixType& matrix) 00276 { 00277 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()"); 00278 cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>()); 00279 cholmod_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod); 00280 00281 // If the factorization failed, minor is the column at which it did. On success minor == n. 00282 this->m_info = (m_cholmodFactor->minor == m_cholmodFactor->n ? Success : NumericalIssue); 00283 m_factorizationIsOk = true; 00284 } 00285 00288 cholmod_common& cholmod() { return m_cholmod; } 00289 00290 #ifndef EIGEN_PARSED_BY_DOXYGEN 00291 00292 template<typename Rhs,typename Dest> 00293 void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const 00294 { 00295 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 00296 const Index size = m_cholmodFactor->n; 00297 EIGEN_UNUSED_VARIABLE(size); 00298 eigen_assert(size==b.rows()); 00299 00300 // note: cd stands for Cholmod Dense 00301 Rhs& b_ref(b.const_cast_derived()); 00302 cholmod_dense b_cd = viewAsCholmod(b_ref); 00303 cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod); 00304 if(!x_cd) 00305 { 00306 this->m_info = NumericalIssue; 00307 } 00308 // TODO optimize this copy by swapping when possible (be careful with alignment, etc.) 00309 dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols()); 00310 cholmod_free_dense(&x_cd, &m_cholmod); 00311 } 00312 00314 template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex> 00315 void _solve(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const 00316 { 00317 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()"); 00318 const Index size = m_cholmodFactor->n; 00319 EIGEN_UNUSED_VARIABLE(size); 00320 eigen_assert(size==b.rows()); 00321 00322 // note: cs stands for Cholmod Sparse 00323 cholmod_sparse b_cs = viewAsCholmod(b); 00324 cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod); 00325 if(!x_cs) 00326 { 00327 this->m_info = NumericalIssue; 00328 } 00329 // TODO optimize this copy by swapping when possible (be careful with alignment, etc.) 00330 dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs); 00331 cholmod_free_sparse(&x_cs, &m_cholmod); 00332 } 00333 #endif // EIGEN_PARSED_BY_DOXYGEN 00334 00335 00345 Derived& setShift(const RealScalar& offset) 00346 { 00347 m_shiftOffset[0] = offset; 00348 return derived(); 00349 } 00350 00351 template<typename Stream> 00352 void dumpMemory(Stream& /*s*/) 00353 {} 00354 00355 protected: 00356 mutable cholmod_common m_cholmod; 00357 cholmod_factor* m_cholmodFactor; 00358 RealScalar m_shiftOffset[2]; 00359 mutable ComputationInfo m_info; 00360 bool m_isInitialized; 00361 int m_factorizationIsOk; 00362 int m_analysisIsOk; 00363 }; 00364 00383 template<typename _MatrixType, int _UpLo = Lower> 00384 class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> > 00385 { 00386 typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT> Base; 00387 using Base::m_cholmod; 00388 00389 public: 00390 00391 typedef _MatrixType MatrixType; 00392 00393 CholmodSimplicialLLT() : Base() { init(); } 00394 00395 CholmodSimplicialLLT(const MatrixType& matrix) : Base() 00396 { 00397 init(); 00398 compute(matrix); 00399 } 00400 00401 ~CholmodSimplicialLLT() {} 00402 protected: 00403 void init() 00404 { 00405 m_cholmod.final_asis = 0; 00406 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL; 00407 m_cholmod.final_ll = 1; 00408 } 00409 }; 00410 00411 00430 template<typename _MatrixType, int _UpLo = Lower> 00431 class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> > 00432 { 00433 typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT> Base; 00434 using Base::m_cholmod; 00435 00436 public: 00437 00438 typedef _MatrixType MatrixType; 00439 00440 CholmodSimplicialLDLT() : Base() { init(); } 00441 00442 CholmodSimplicialLDLT(const MatrixType& matrix) : Base() 00443 { 00444 init(); 00445 compute(matrix); 00446 } 00447 00448 ~CholmodSimplicialLDLT() {} 00449 protected: 00450 void init() 00451 { 00452 m_cholmod.final_asis = 1; 00453 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL; 00454 } 00455 }; 00456 00475 template<typename _MatrixType, int _UpLo = Lower> 00476 class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> > 00477 { 00478 typedef CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT> Base; 00479 using Base::m_cholmod; 00480 00481 public: 00482 00483 typedef _MatrixType MatrixType; 00484 00485 CholmodSupernodalLLT() : Base() { init(); } 00486 00487 CholmodSupernodalLLT(const MatrixType& matrix) : Base() 00488 { 00489 init(); 00490 compute(matrix); 00491 } 00492 00493 ~CholmodSupernodalLLT() {} 00494 protected: 00495 void init() 00496 { 00497 m_cholmod.final_asis = 1; 00498 m_cholmod.supernodal = CHOLMOD_SUPERNODAL; 00499 } 00500 }; 00501 00522 template<typename _MatrixType, int _UpLo = Lower> 00523 class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> > 00524 { 00525 typedef CholmodBase<_MatrixType, _UpLo, CholmodDecomposition> Base; 00526 using Base::m_cholmod; 00527 00528 public: 00529 00530 typedef _MatrixType MatrixType; 00531 00532 CholmodDecomposition() : Base() { init(); } 00533 00534 CholmodDecomposition(const MatrixType& matrix) : Base() 00535 { 00536 init(); 00537 compute(matrix); 00538 } 00539 00540 ~CholmodDecomposition() {} 00541 00542 void setMode(CholmodMode mode) 00543 { 00544 switch(mode) 00545 { 00546 case CholmodAuto: 00547 m_cholmod.final_asis = 1; 00548 m_cholmod.supernodal = CHOLMOD_AUTO; 00549 break; 00550 case CholmodSimplicialLLt: 00551 m_cholmod.final_asis = 0; 00552 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL; 00553 m_cholmod.final_ll = 1; 00554 break; 00555 case CholmodSupernodalLLt: 00556 m_cholmod.final_asis = 1; 00557 m_cholmod.supernodal = CHOLMOD_SUPERNODAL; 00558 break; 00559 case CholmodLDLt: 00560 m_cholmod.final_asis = 1; 00561 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL; 00562 break; 00563 default: 00564 break; 00565 } 00566 } 00567 protected: 00568 void init() 00569 { 00570 m_cholmod.final_asis = 1; 00571 m_cholmod.supernodal = CHOLMOD_AUTO; 00572 } 00573 }; 00574 00575 namespace internal { 00576 00577 template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs> 00578 struct solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs> 00579 : solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs> 00580 { 00581 typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec; 00582 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 00583 00584 template<typename Dest> void evalTo(Dest& dst) const 00585 { 00586 dec()._solve(rhs(),dst); 00587 } 00588 }; 00589 00590 template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs> 00591 struct sparse_solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs> 00592 : sparse_solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs> 00593 { 00594 typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec; 00595 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) 00596 00597 template<typename Dest> void evalTo(Dest& dst) const 00598 { 00599 dec()._solve(rhs(),dst); 00600 } 00601 }; 00602 00603 } // end namespace internal 00604 00605 } // end namespace Eigen 00606 00607 #endif // EIGEN_CHOLMODSUPPORT_H