10 #ifndef EIGEN_CHOLMODSUPPORT_H
11 #define EIGEN_CHOLMODSUPPORT_H
17 template<
typename Scalar,
typename CholmodType>
18 void cholmod_configure_matrix(CholmodType& mat)
20 if (internal::is_same<Scalar,float>::value)
22 mat.xtype = CHOLMOD_REAL;
23 mat.dtype = CHOLMOD_SINGLE;
25 else if (internal::is_same<Scalar,double>::value)
27 mat.xtype = CHOLMOD_REAL;
28 mat.dtype = CHOLMOD_DOUBLE;
30 else if (internal::is_same<Scalar,std::complex<float> >::value)
32 mat.xtype = CHOLMOD_COMPLEX;
33 mat.dtype = CHOLMOD_SINGLE;
35 else if (internal::is_same<Scalar,std::complex<double> >::value)
37 mat.xtype = CHOLMOD_COMPLEX;
38 mat.dtype = CHOLMOD_DOUBLE;
42 eigen_assert(
false &&
"Scalar type not supported by CHOLMOD");
51 template<
typename _Scalar,
int _Options,
typename _Index>
56 res.nrow = mat.
rows();;
57 res.ncol = mat.
cols();
75 if (internal::is_same<_Index,int>::value)
77 res.itype = CHOLMOD_INT;
79 else if (internal::is_same<_Index,UF_long>::value)
81 res.itype = CHOLMOD_LONG;
85 eigen_assert(
false &&
"Index type not supported yet");
89 internal::cholmod_configure_matrix<_Scalar>(res);
96 template<
typename _Scalar,
int _Options,
typename _Index>
97 const cholmod_sparse
viewAsCholmod(
const SparseMatrix<_Scalar,_Options,_Index>& mat)
105 template<
typename _Scalar,
int _Options,
typename _Index,
unsigned int UpLo>
108 cholmod_sparse res =
viewAsCholmod(mat.matrix().const_cast_derived());
110 if(UpLo==
Upper) res.stype = 1;
111 if(UpLo==
Lower) res.stype = -1;
118 template<
typename Derived>
121 EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&
RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
122 typedef typename Derived::Scalar Scalar;
125 res.nrow = mat.rows();
126 res.ncol = mat.cols();
127 res.nzmax = res.nrow * res.ncol;
128 res.d = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
129 res.x = (
void*)(mat.derived().data());
132 internal::cholmod_configure_matrix<Scalar>(res);
139 template<
typename Scalar,
int Flags,
typename Index>
143 (cm.nrow, cm.ncol,
static_cast<Index*
>(cm.p)[cm.ncol],
144 static_cast<Index*>(cm.p),
static_cast<Index*
>(cm.i),static_cast<Scalar*>(cm.x) );
148 CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
157 template<
typename _MatrixType,
int _UpLo,
typename Derived>
161 typedef _MatrixType MatrixType;
162 enum { UpLo = _UpLo };
163 typedef typename MatrixType::Scalar Scalar;
164 typedef typename MatrixType::RealScalar RealScalar;
165 typedef MatrixType CholMatrixType;
166 typedef typename MatrixType::Index Index;
171 : m_cholmodFactor(0), m_info(
Success), m_isInitialized(
false)
173 cholmod_start(&m_cholmod);
177 : m_cholmodFactor(0), m_info(
Success), m_isInitialized(
false)
179 m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
180 cholmod_start(&m_cholmod);
187 cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
188 cholmod_finish(&m_cholmod);
191 inline Index cols()
const {
return m_cholmodFactor->n; }
192 inline Index rows()
const {
return m_cholmodFactor->n; }
194 Derived& derived() {
return *
static_cast<Derived*
>(
this); }
195 const Derived& derived()
const {
return *
static_cast<const Derived*
>(
this); }
204 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
220 template<
typename Rhs>
221 inline const internal::solve_retval<CholmodBase, Rhs>
224 eigen_assert(m_isInitialized &&
"LLT is not initialized.");
225 eigen_assert(rows()==b.rows()
226 &&
"CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
227 return internal::solve_retval<CholmodBase, Rhs>(*
this, b.derived());
234 template<
typename Rhs>
235 inline const internal::sparse_solve_retval<CholmodBase, Rhs>
238 eigen_assert(m_isInitialized &&
"LLT is not initialized.");
239 eigen_assert(rows()==b.
rows()
240 &&
"CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
241 return internal::sparse_solve_retval<CholmodBase, Rhs>(*
this, b.
derived());
254 cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
257 cholmod_sparse A =
viewAsCholmod(matrix.template selfadjointView<UpLo>());
258 m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
260 this->m_isInitialized =
true;
262 m_analysisIsOk =
true;
263 m_factorizationIsOk =
false;
274 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
275 cholmod_sparse A =
viewAsCholmod(matrix.template selfadjointView<UpLo>());
276 cholmod_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod);
280 m_factorizationIsOk =
true;
285 cholmod_common&
cholmod() {
return m_cholmod; }
287 #ifndef EIGEN_PARSED_BY_DOXYGEN
289 template<
typename Rhs,
typename Dest>
292 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
293 const Index size = m_cholmodFactor->n;
294 EIGEN_UNUSED_VARIABLE(size);
295 eigen_assert(size==b.rows());
298 Rhs& b_ref(b.const_cast_derived());
300 cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
306 dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
307 cholmod_free_dense(&x_cd, &m_cholmod);
311 template<
typename RhsScalar,
int RhsOptions,
typename RhsIndex,
typename DestScalar,
int DestOptions,
typename DestIndex>
312 void _solve(
const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest)
const
314 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
315 const Index size = m_cholmodFactor->n;
316 EIGEN_UNUSED_VARIABLE(size);
317 eigen_assert(size==b.rows());
321 cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
327 dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
328 cholmod_free_sparse(&x_cs, &m_cholmod);
330 #endif // EIGEN_PARSED_BY_DOXYGEN
344 m_shiftOffset[0] = offset;
348 template<
typename Stream>
349 void dumpMemory(Stream& )
353 mutable cholmod_common m_cholmod;
354 cholmod_factor* m_cholmodFactor;
355 RealScalar m_shiftOffset[2];
357 bool m_isInitialized;
358 int m_factorizationIsOk;
380 template<
typename _MatrixType,
int _UpLo = Lower>
384 using Base::m_cholmod;
388 typedef _MatrixType MatrixType;
402 m_cholmod.final_asis = 0;
403 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
404 m_cholmod.final_ll = 1;
427 template<
typename _MatrixType,
int _UpLo = Lower>
431 using Base::m_cholmod;
435 typedef _MatrixType MatrixType;
449 m_cholmod.final_asis = 1;
450 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
472 template<
typename _MatrixType,
int _UpLo = Lower>
476 using Base::m_cholmod;
480 typedef _MatrixType MatrixType;
494 m_cholmod.final_asis = 1;
495 m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
519 template<
typename _MatrixType,
int _UpLo = Lower>
523 using Base::m_cholmod;
527 typedef _MatrixType MatrixType;
539 void setMode(CholmodMode mode)
544 m_cholmod.final_asis = 1;
545 m_cholmod.supernodal = CHOLMOD_AUTO;
547 case CholmodSimplicialLLt:
548 m_cholmod.final_asis = 0;
549 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
550 m_cholmod.final_ll = 1;
552 case CholmodSupernodalLLt:
553 m_cholmod.final_asis = 1;
554 m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
557 m_cholmod.final_asis = 1;
558 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
567 m_cholmod.final_asis = 1;
568 m_cholmod.supernodal = CHOLMOD_AUTO;
574 template<
typename _MatrixType,
int _UpLo,
typename Derived,
typename Rhs>
575 struct solve_retval<
CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
576 : solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
579 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
581 template<typename Dest>
void evalTo(Dest& dst)
const
583 dec()._solve(rhs(),dst);
587 template<
typename _MatrixType,
int _UpLo,
typename Derived,
typename Rhs>
588 struct sparse_solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
589 : sparse_solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
591 typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
592 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
594 template<typename Dest>
void evalTo(Dest& dst)
const
596 dec()._solve(rhs(),dst);
604 #endif // EIGEN_CHOLMODSUPPORT_H