00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
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 }
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
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
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
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
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
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 }
00591
00592 }
00593
00594 #endif // EIGEN_CHOLMODSUPPORT_H