10 #ifndef EIGEN_CHOLMODSUPPORT_H
11 #define EIGEN_CHOLMODSUPPORT_H
17 template<
typename Scalar>
struct cholmod_configure_matrix;
19 template<>
struct cholmod_configure_matrix<double> {
20 template<
typename CholmodType>
21 static void run(CholmodType& mat) {
22 mat.xtype = CHOLMOD_REAL;
23 mat.dtype = CHOLMOD_DOUBLE;
27 template<>
struct cholmod_configure_matrix<std::complex<double> > {
28 template<
typename CholmodType>
29 static void run(CholmodType& mat) {
30 mat.xtype = CHOLMOD_COMPLEX;
31 mat.dtype = CHOLMOD_DOUBLE;
57 template<
typename _Scalar,
int _Options,
typename _StorageIndex>
61 res.nzmax = mat.nonZeros();
62 res.nrow = mat.rows();
63 res.ncol = mat.cols();
64 res.p = mat.outerIndexPtr();
65 res.i = mat.innerIndexPtr();
66 res.x = mat.valuePtr();
69 if(mat.isCompressed())
77 res.nz = mat.innerNonZeroPtr();
83 if (internal::is_same<_StorageIndex,int>::value)
85 res.itype = CHOLMOD_INT;
87 else if (internal::is_same<_StorageIndex,long>::value)
89 res.itype = CHOLMOD_LONG;
93 eigen_assert(
false &&
"Index type not supported yet");
97 internal::cholmod_configure_matrix<_Scalar>::run(res);
104 template<
typename _Scalar,
int _Options,
typename _Index>
105 const cholmod_sparse
viewAsCholmod(
const SparseMatrix<_Scalar,_Options,_Index>& mat)
107 cholmod_sparse res =
viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_Index> >(mat.const_cast_derived()));
111 template<
typename _Scalar,
int _Options,
typename _Index>
112 const cholmod_sparse
viewAsCholmod(
const SparseVector<_Scalar,_Options,_Index>& mat)
114 cholmod_sparse res =
viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_Index> >(mat.const_cast_derived()));
120 template<
typename _Scalar,
int _Options,
typename _Index,
unsigned int UpLo>
125 if(UpLo==
Upper) res.stype = 1;
126 if(UpLo==
Lower) res.stype = -1;
133 template<
typename Derived>
136 EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&
RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
137 typedef typename Derived::Scalar Scalar;
140 res.nrow = mat.
rows();
141 res.ncol = mat.
cols();
142 res.nzmax = res.nrow * res.ncol;
143 res.d = Derived::IsVectorAtCompileTime ? mat.
derived().size() : mat.
derived().outerStride();
144 res.x = (
void*)(mat.
derived().data());
147 internal::cholmod_configure_matrix<Scalar>::run(res);
154 template<
typename Scalar,
int Flags,
typename StorageIndex>
158 (cm.nrow, cm.ncol, static_cast<StorageIndex*>(cm.p)[cm.ncol],
159 static_cast<StorageIndex*>(cm.p), static_cast<StorageIndex*>(cm.i),static_cast<Scalar*>(cm.x) );
163 CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
172 template<
typename _MatrixType,
int _UpLo,
typename Derived>
178 using Base::m_isInitialized;
180 typedef _MatrixType MatrixType;
181 enum { UpLo = _UpLo };
182 typedef typename MatrixType::Scalar Scalar;
183 typedef typename MatrixType::RealScalar RealScalar;
184 typedef MatrixType CholMatrixType;
185 typedef typename MatrixType::StorageIndex StorageIndex;
187 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
188 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
194 : m_cholmodFactor(0), m_info(
Success), m_factorizationIsOk(
false), m_analysisIsOk(
false)
196 EIGEN_STATIC_ASSERT((internal::is_same<double,RealScalar>::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY);
197 m_shiftOffset[0] = m_shiftOffset[1] = 0.0;
198 cholmod_start(&m_cholmod);
202 : m_cholmodFactor(0), m_info(
Success), m_factorizationIsOk(
false), m_analysisIsOk(
false)
204 EIGEN_STATIC_ASSERT((internal::is_same<double,RealScalar>::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY);
205 m_shiftOffset[0] = m_shiftOffset[1] = 0.0;
206 cholmod_start(&m_cholmod);
213 cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
214 cholmod_finish(&m_cholmod);
217 inline StorageIndex cols()
const {
return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
218 inline StorageIndex rows()
const {
return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
227 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
249 cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
252 cholmod_sparse A =
viewAsCholmod(matrix.template selfadjointView<UpLo>());
253 m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
255 this->m_isInitialized =
true;
257 m_analysisIsOk =
true;
258 m_factorizationIsOk =
false;
269 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
270 cholmod_sparse A =
viewAsCholmod(matrix.template selfadjointView<UpLo>());
271 cholmod_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod);
275 m_factorizationIsOk =
true;
280 cholmod_common&
cholmod() {
return m_cholmod; }
282 #ifndef EIGEN_PARSED_BY_DOXYGEN
284 template<
typename Rhs,
typename Dest>
287 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
288 const Index size = m_cholmodFactor->n;
289 EIGEN_UNUSED_VARIABLE(size);
290 eigen_assert(size==b.
rows());
296 cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
304 cholmod_free_dense(&x_cd, &m_cholmod);
308 template<
typename RhsDerived,
typename DestDerived>
309 void _solve_impl(
const SparseMatrixBase<RhsDerived> &b, SparseMatrixBase<DestDerived> &dest)
const
311 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
312 const Index size = m_cholmodFactor->n;
313 EIGEN_UNUSED_VARIABLE(size);
314 eigen_assert(size==b.rows());
317 Ref<SparseMatrix<typename RhsDerived::Scalar,ColMajor,typename RhsDerived::StorageIndex> > b_ref(b.const_cast_derived());
319 cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
326 dest.derived() = viewAsEigen<typename DestDerived::Scalar,ColMajor,typename DestDerived::StorageIndex>(*x_cs);
327 cholmod_free_sparse(&x_cs, &m_cholmod);
329 #endif // EIGEN_PARSED_BY_DOXYGEN
343 m_shiftOffset[0] = double(offset);
359 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
361 RealScalar logDet = 0;
362 Scalar *x = static_cast<Scalar*>(m_cholmodFactor->x);
363 if (m_cholmodFactor->is_super)
369 StorageIndex *super = static_cast<StorageIndex*>(m_cholmodFactor->super);
371 StorageIndex *pi = static_cast<StorageIndex*>(m_cholmodFactor->pi);
373 StorageIndex *px = static_cast<StorageIndex*>(m_cholmodFactor->px);
375 Index nb_super_nodes = m_cholmodFactor->nsuper;
376 for (
Index k=0; k < nb_super_nodes; ++k)
378 StorageIndex ncols = super[k + 1] - super[k];
379 StorageIndex nrows = pi[k + 1] - pi[k];
382 logDet += sk.real().log().sum();
388 StorageIndex *p = static_cast<StorageIndex*>(m_cholmodFactor->p);
389 Index size = m_cholmodFactor->n;
390 for (
Index k=0; k<size; ++k)
391 logDet +=
log(
real( x[p[k]] ));
393 if (m_cholmodFactor->is_ll)
398 template<
typename Stream>
399 void dumpMemory(Stream& )
403 mutable cholmod_common m_cholmod;
404 cholmod_factor* m_cholmodFactor;
405 double m_shiftOffset[2];
407 int m_factorizationIsOk;
433 template<
typename _MatrixType,
int _UpLo = Lower>
437 using Base::m_cholmod;
441 typedef _MatrixType MatrixType;
455 m_cholmod.final_asis = 0;
456 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
457 m_cholmod.final_ll = 1;
484 template<
typename _MatrixType,
int _UpLo = Lower>
488 using Base::m_cholmod;
492 typedef _MatrixType MatrixType;
506 m_cholmod.final_asis = 1;
507 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
533 template<
typename _MatrixType,
int _UpLo = Lower>
537 using Base::m_cholmod;
541 typedef _MatrixType MatrixType;
555 m_cholmod.final_asis = 1;
556 m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
584 template<
typename _MatrixType,
int _UpLo = Lower>
588 using Base::m_cholmod;
592 typedef _MatrixType MatrixType;
604 void setMode(CholmodMode mode)
609 m_cholmod.final_asis = 1;
610 m_cholmod.supernodal = CHOLMOD_AUTO;
612 case CholmodSimplicialLLt:
613 m_cholmod.final_asis = 0;
614 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
615 m_cholmod.final_ll = 1;
617 case CholmodSupernodalLLt:
618 m_cholmod.final_asis = 1;
619 m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
622 m_cholmod.final_asis = 1;
623 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
632 m_cholmod.final_asis = 1;
633 m_cholmod.supernodal = CHOLMOD_AUTO;
639 #endif // EIGEN_CHOLMODSUPPORT_H