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>
57 res.nrow = mat.
rows();;
58 res.ncol = mat.
cols();
76 if (internal::is_same<_Index,int>::value)
78 res.itype = CHOLMOD_INT;
82 eigen_assert(
false &&
"Index type different than int is not supported yet");
86 internal::cholmod_configure_matrix<_Scalar>(res);
93 template<
typename _Scalar,
int _Options,
typename _Index>
94 const cholmod_sparse
viewAsCholmod(
const SparseMatrix<_Scalar,_Options,_Index>& mat)
102 template<
typename _Scalar,
int _Options,
typename _Index,
unsigned int UpLo>
105 cholmod_sparse res =
viewAsCholmod(mat.matrix().const_cast_derived());
107 if(UpLo==
Upper) res.stype = 1;
108 if(UpLo==
Lower) res.stype = -1;
115 template<
typename Derived>
118 EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&
RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
119 typedef typename Derived::Scalar Scalar;
122 res.nrow = mat.rows();
123 res.ncol = mat.cols();
124 res.nzmax = res.nrow * res.ncol;
125 res.d = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
126 res.x = mat.derived().data();
129 internal::cholmod_configure_matrix<Scalar>(res);
136 template<
typename Scalar,
int Flags,
typename Index>
140 (cm.nrow, cm.ncol,
reinterpret_cast<Index*
>(cm.p)[cm.ncol],
141 reinterpret_cast<Index*>(cm.p),
reinterpret_cast<Index*
>(cm.i),reinterpret_cast<Scalar*>(cm.x) );
145 CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
154 template<
typename _MatrixType,
int _UpLo,
typename Derived>
158 typedef _MatrixType MatrixType;
159 enum { UpLo = _UpLo };
160 typedef typename MatrixType::Scalar Scalar;
161 typedef typename MatrixType::RealScalar RealScalar;
162 typedef MatrixType CholMatrixType;
163 typedef typename MatrixType::Index Index;
168 : m_cholmodFactor(0), m_info(
Success), m_isInitialized(
false)
170 cholmod_start(&m_cholmod);
174 : m_cholmodFactor(0), m_info(
Success), m_isInitialized(
false)
176 cholmod_start(&m_cholmod);
183 cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
184 cholmod_finish(&m_cholmod);
187 inline Index cols()
const {
return m_cholmodFactor->n; }
188 inline Index rows()
const {
return m_cholmodFactor->n; }
190 Derived& derived() {
return *
static_cast<Derived*
>(
this); }
191 const Derived& derived()
const {
return *
static_cast<const Derived*
>(
this); }
200 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
216 template<
typename Rhs>
217 inline const internal::solve_retval<CholmodBase, Rhs>
220 eigen_assert(m_isInitialized &&
"LLT is not initialized.");
221 eigen_assert(rows()==b.rows()
222 &&
"CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
223 return internal::solve_retval<CholmodBase, Rhs>(*
this, b.derived());
230 template<
typename Rhs>
231 inline const internal::sparse_solve_retval<CholmodBase, Rhs>
234 eigen_assert(m_isInitialized &&
"LLT is not initialized.");
235 eigen_assert(rows()==b.
rows()
236 &&
"CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
237 return internal::sparse_solve_retval<CholmodBase, Rhs>(*
this, b.
derived());
250 cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
253 cholmod_sparse A =
viewAsCholmod(matrix.template selfadjointView<UpLo>());
254 m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
256 this->m_isInitialized =
true;
258 m_analysisIsOk =
true;
259 m_factorizationIsOk =
false;
270 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
271 cholmod_sparse A =
viewAsCholmod(matrix.template selfadjointView<UpLo>());
272 cholmod_factorize(&A, 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_assert(size==b.rows());
293 cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
299 dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
300 cholmod_free_dense(&x_cd, &m_cholmod);
304 template<
typename RhsScalar,
int RhsOptions,
typename RhsIndex,
typename DestScalar,
int DestOptions,
typename DestIndex>
305 void _solve(
const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest)
const
307 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
308 const Index size = m_cholmodFactor->n;
309 eigen_assert(size==b.rows());
313 cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
319 dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
320 cholmod_free_sparse(&x_cs, &m_cholmod);
322 #endif // EIGEN_PARSED_BY_DOXYGEN
324 template<
typename Stream>
325 void dumpMemory(Stream& s)
329 mutable cholmod_common m_cholmod;
330 cholmod_factor* m_cholmodFactor;
332 bool m_isInitialized;
333 int m_factorizationIsOk;
355 template<
typename _MatrixType,
int _UpLo = Lower>
359 using Base::m_cholmod;
363 typedef _MatrixType MatrixType;
377 m_cholmod.final_asis = 0;
378 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
379 m_cholmod.final_ll = 1;
402 template<
typename _MatrixType,
int _UpLo = Lower>
406 using Base::m_cholmod;
410 typedef _MatrixType MatrixType;
424 m_cholmod.final_asis = 1;
425 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
447 template<
typename _MatrixType,
int _UpLo = Lower>
451 using Base::m_cholmod;
455 typedef _MatrixType MatrixType;
469 m_cholmod.final_asis = 1;
470 m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
494 template<
typename _MatrixType,
int _UpLo = Lower>
498 using Base::m_cholmod;
502 typedef _MatrixType MatrixType;
514 void setMode(CholmodMode mode)
519 m_cholmod.final_asis = 1;
520 m_cholmod.supernodal = CHOLMOD_AUTO;
522 case CholmodSimplicialLLt:
523 m_cholmod.final_asis = 0;
524 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
525 m_cholmod.final_ll = 1;
527 case CholmodSupernodalLLt:
528 m_cholmod.final_asis = 1;
529 m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
532 m_cholmod.final_asis = 1;
533 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
542 m_cholmod.final_asis = 1;
543 m_cholmod.supernodal = CHOLMOD_AUTO;
549 template<
typename _MatrixType,
int _UpLo,
typename Derived,
typename Rhs>
550 struct solve_retval<
CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
551 : solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
554 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
556 template<typename Dest>
void evalTo(Dest& dst)
const
558 dec()._solve(rhs(),dst);
562 template<
typename _MatrixType,
int _UpLo,
typename Derived,
typename Rhs>
563 struct sparse_solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
564 : sparse_solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
566 typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
567 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
569 template<typename Dest>
void evalTo(Dest& dst)
const
571 dec()._solve(rhs(),dst);
579 #endif // EIGEN_CHOLMODSUPPORT_H