![]() |
Eigen
3.3.3
|
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