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
void factorize(const MatrixType &matrix)
Definition: CholmodSupport.h:272
Index rows() const
Definition: SparseMatrix.h:119
void analyzePattern(const MatrixType &matrix)
Definition: CholmodSupport.h:250
Index cols() const
Definition: SparseMatrix.h:121
const Scalar * valuePtr() const
Definition: SparseMatrix.h:131
A versatible sparse matrix representation.
Definition: SparseMatrix.h:85
bool isCompressed() const
Definition: SparseMatrix.h:116
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: CholmodSupport.h:202
Definition: Constants.h:378
const internal::sparse_solve_retval< CholmodBase, Rhs > solve(const SparseMatrixBase< Rhs > &b) const
Definition: CholmodSupport.h:236
Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix.
Definition: SparseSelfAdjointView.h:49
const Index * innerNonZeroPtr() const
Definition: SparseMatrix.h:158
A supernodal Cholesky (LLT) factorization and solver based on Cholmod.
Definition: CholmodSupport.h:473
const Index * outerIndexPtr() const
Definition: SparseMatrix.h:149
MappedSparseMatrix< Scalar, Flags, Index > viewAsEigen(cholmod_sparse &cm)
Definition: CholmodSupport.h:140
Index nonZeros() const
Definition: SparseMatrix.h:246
Base class of any sparse matrices or sparse expressions.
Definition: SparseMatrixBase.h:26
Derived & setShift(const RealScalar &offset)
Definition: CholmodSupport.h:342
Definition: Constants.h:169
cholmod_sparse viewAsCholmod(SparseMatrix< _Scalar, _Options, _Index > &mat)
Definition: CholmodSupport.h:52
Derived & derived()
Definition: EigenBase.h:34
Definition: Constants.h:167
A general Cholesky factorization and solver based on Cholmod.
Definition: CholmodSupport.h:520
const internal::solve_retval< CholmodBase, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: CholmodSupport.h:222
The base class for the direct Cholesky factorization of Cholmod.
Definition: CholmodSupport.h:158
A simplicial direct Cholesky (LLT) factorization and solver based on Cholmod.
Definition: CholmodSupport.h:381
A simplicial direct Cholesky (LDLT) factorization and solver based on Cholmod.
Definition: CholmodSupport.h:428
Definition: Constants.h:376
const unsigned int RowMajorBit
Definition: Constants.h:53
Index rows() const
Definition: SparseMatrixBase.h:150
ComputationInfo
Definition: Constants.h:374
Sparse matrix.
Definition: MappedSparseMatrix.h:31
Derived & compute(const MatrixType &matrix)
Definition: CholmodSupport.h:209
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
cholmod_common & cholmod()
Definition: CholmodSupport.h:285
const Index * innerIndexPtr() const
Definition: SparseMatrix.h:140