10#ifndef EIGEN_CHOLMODSUPPORT_H
11#define EIGEN_CHOLMODSUPPORT_H
17template<
typename Scalar>
struct cholmod_configure_matrix;
19template<>
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;
27template<>
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;
57template<
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,SuiteSparse_long>::value)
89 res.itype = CHOLMOD_LONG;
93 eigen_assert(
false &&
"Index type not supported yet");
97 internal::cholmod_configure_matrix<_Scalar>::run(res);
104template<
typename _Scalar,
int _Options,
typename _Index>
105const cholmod_sparse viewAsCholmod(
const SparseMatrix<_Scalar,_Options,_Index>& mat)
107 cholmod_sparse res = viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_Index> >(mat.const_cast_derived()));
111template<
typename _Scalar,
int _Options,
typename _Index>
112const cholmod_sparse
viewAsCholmod(
const SparseVector<_Scalar,_Options,_Index>& mat)
114 cholmod_sparse res =
viewAsCholmod(Ref<SparseMatrix<_Scalar,_Options,_Index> >(mat.const_cast_derived()));
120template<
typename _Scalar,
int _Options,
typename _Index,
unsigned int UpLo>
125 if(UpLo==
Upper) res.stype = 1;
126 if(UpLo==
Lower) res.stype = -1;
136template<
typename Derived>
139 EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&
RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
140 typedef typename Derived::Scalar Scalar;
143 res.nrow = mat.
rows();
144 res.ncol = mat.
cols();
145 res.nzmax = res.nrow * res.ncol;
146 res.d = Derived::IsVectorAtCompileTime ? mat.
derived().size() : mat.
derived().outerStride();
147 res.x = (
void*)(mat.
derived().data());
150 internal::cholmod_configure_matrix<Scalar>::run(res);
157template<
typename Scalar,
int Flags,
typename StorageIndex>
161 (cm.nrow, cm.ncol,
static_cast<StorageIndex*
>(cm.p)[cm.ncol],
162 static_cast<StorageIndex*
>(cm.p),
static_cast<StorageIndex*
>(cm.i),
static_cast<Scalar*
>(cm.x) );
169#define EIGEN_CHOLMOD_SPECIALIZE0(ret, name) \
170 template<typename _StorageIndex> inline ret cm_ ## name (cholmod_common &Common) { return cholmod_ ## name (&Common); } \
171 template<> inline ret cm_ ## name<SuiteSparse_long> (cholmod_common &Common) { return cholmod_l_ ## name (&Common); }
173#define EIGEN_CHOLMOD_SPECIALIZE1(ret, name, t1, a1) \
174 template<typename _StorageIndex> inline ret cm_ ## name (t1& a1, cholmod_common &Common) { return cholmod_ ## name (&a1, &Common); } \
175 template<> inline ret cm_ ## name<SuiteSparse_long> (t1& a1, cholmod_common &Common) { return cholmod_l_ ## name (&a1, &Common); }
177EIGEN_CHOLMOD_SPECIALIZE0(
int, start)
178EIGEN_CHOLMOD_SPECIALIZE0(
int, finish)
180EIGEN_CHOLMOD_SPECIALIZE1(
int, free_factor, cholmod_factor*, L)
181EIGEN_CHOLMOD_SPECIALIZE1(
int, free_dense, cholmod_dense*, X)
182EIGEN_CHOLMOD_SPECIALIZE1(
int, free_sparse, cholmod_sparse*, A)
184EIGEN_CHOLMOD_SPECIALIZE1(cholmod_factor*, analyze, cholmod_sparse, A)
186template<
typename _StorageIndex>
inline cholmod_dense* cm_solve (
int sys, cholmod_factor& L, cholmod_dense& B, cholmod_common &Common) {
return cholmod_solve (sys, &L, &B, &Common); }
187template<>
inline cholmod_dense* cm_solve<SuiteSparse_long> (
int sys, cholmod_factor& L, cholmod_dense& B, cholmod_common &Common) {
return cholmod_l_solve (sys, &L, &B, &Common); }
189template<
typename _StorageIndex>
inline cholmod_sparse* cm_spsolve (
int sys, cholmod_factor& L, cholmod_sparse& B, cholmod_common &Common) {
return cholmod_spsolve (sys, &L, &B, &Common); }
190template<>
inline cholmod_sparse* cm_spsolve<SuiteSparse_long> (
int sys, cholmod_factor& L, cholmod_sparse& B, cholmod_common &Common) {
return cholmod_l_spsolve (sys, &L, &B, &Common); }
192template<
typename _StorageIndex>
193inline int cm_factorize_p (cholmod_sparse* A,
double beta[2], _StorageIndex* fset, std::size_t fsize, cholmod_factor* L, cholmod_common &Common) {
return cholmod_factorize_p (A, beta, fset, fsize, L, &Common); }
195inline int cm_factorize_p<SuiteSparse_long> (cholmod_sparse* A,
double beta[2], SuiteSparse_long* fset, std::size_t fsize, cholmod_factor* L, cholmod_common &Common) {
return cholmod_l_factorize_p (A, beta, fset, fsize, L, &Common); }
197#undef EIGEN_CHOLMOD_SPECIALIZE0
198#undef EIGEN_CHOLMOD_SPECIALIZE1
204 CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
213template<
typename _MatrixType,
int _UpLo,
typename Derived>
219 using Base::m_isInitialized;
221 typedef _MatrixType MatrixType;
222 enum { UpLo = _UpLo };
223 typedef typename MatrixType::Scalar Scalar;
224 typedef typename MatrixType::RealScalar RealScalar;
225 typedef MatrixType CholMatrixType;
226 typedef typename MatrixType::StorageIndex StorageIndex;
228 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
229 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
235 : m_cholmodFactor(0), m_info(
Success), m_factorizationIsOk(
false), m_analysisIsOk(
false)
237 EIGEN_STATIC_ASSERT((internal::is_same<double,RealScalar>::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY);
238 m_shiftOffset[0] = m_shiftOffset[1] = 0.0;
239 internal::cm_start<StorageIndex>(m_cholmod);
243 : m_cholmodFactor(0), m_info(
Success), m_factorizationIsOk(
false), m_analysisIsOk(
false)
245 EIGEN_STATIC_ASSERT((internal::is_same<double,RealScalar>::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY);
246 m_shiftOffset[0] = m_shiftOffset[1] = 0.0;
247 internal::cm_start<StorageIndex>(m_cholmod);
254 internal::cm_free_factor<StorageIndex>(m_cholmodFactor, m_cholmod);
255 internal::cm_finish<StorageIndex>(m_cholmod);
258 inline StorageIndex cols()
const {
return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
259 inline StorageIndex rows()
const {
return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
268 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
290 internal::cm_free_factor<StorageIndex>(m_cholmodFactor, m_cholmod);
293 cholmod_sparse A =
viewAsCholmod(matrix.template selfadjointView<UpLo>());
294 m_cholmodFactor = internal::cm_analyze<StorageIndex>(A, m_cholmod);
296 this->m_isInitialized =
true;
298 m_analysisIsOk =
true;
299 m_factorizationIsOk =
false;
310 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
311 cholmod_sparse A =
viewAsCholmod(matrix.template selfadjointView<UpLo>());
312 internal::cm_factorize_p<StorageIndex>(&A, m_shiftOffset, 0, 0, m_cholmodFactor, m_cholmod);
316 m_factorizationIsOk =
true;
321 cholmod_common&
cholmod() {
return m_cholmod; }
323 #ifndef EIGEN_PARSED_BY_DOXYGEN
325 template<
typename Rhs,
typename Dest>
328 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
329 const Index size = m_cholmodFactor->n;
330 EIGEN_UNUSED_VARIABLE(size);
331 eigen_assert(size==b.
rows());
337 cholmod_dense* x_cd = internal::cm_solve<StorageIndex>(CHOLMOD_A, *m_cholmodFactor, b_cd, m_cholmod);
346 internal::cm_free_dense<StorageIndex>(x_cd, m_cholmod);
350 template<
typename RhsDerived,
typename DestDerived>
351 void _solve_impl(
const SparseMatrixBase<RhsDerived> &b, SparseMatrixBase<DestDerived> &dest)
const
353 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
354 const Index size = m_cholmodFactor->n;
355 EIGEN_UNUSED_VARIABLE(size);
356 eigen_assert(size==b.rows());
359 Ref<SparseMatrix<typename RhsDerived::Scalar,ColMajor,typename RhsDerived::StorageIndex> > b_ref(b.const_cast_derived());
361 cholmod_sparse* x_cs = internal::cm_spsolve<StorageIndex>(CHOLMOD_A, *m_cholmodFactor, b_cs, m_cholmod);
369 dest.derived() = viewAsEigen<typename DestDerived::Scalar,ColMajor,typename DestDerived::StorageIndex>(*x_cs);
370 internal::cm_free_sparse<StorageIndex>(x_cs, m_cholmod);
386 m_shiftOffset[0] = double(offset);
402 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
404 RealScalar logDet = 0;
405 Scalar *x =
static_cast<Scalar*
>(m_cholmodFactor->x);
406 if (m_cholmodFactor->is_super)
412 StorageIndex *super =
static_cast<StorageIndex*
>(m_cholmodFactor->super);
414 StorageIndex *pi =
static_cast<StorageIndex*
>(m_cholmodFactor->pi);
416 StorageIndex *px =
static_cast<StorageIndex*
>(m_cholmodFactor->px);
418 Index nb_super_nodes = m_cholmodFactor->nsuper;
419 for (
Index k=0; k < nb_super_nodes; ++k)
421 StorageIndex ncols = super[k + 1] - super[k];
422 StorageIndex nrows = pi[k + 1] - pi[k];
425 logDet += sk.real().log().sum();
431 StorageIndex *p =
static_cast<StorageIndex*
>(m_cholmodFactor->p);
432 Index size = m_cholmodFactor->n;
433 for (
Index k=0; k<size; ++k)
434 logDet += log(
real( x[p[k]] ));
436 if (m_cholmodFactor->is_ll)
441 template<
typename Stream>
442 void dumpMemory(Stream& )
446 mutable cholmod_common m_cholmod;
447 cholmod_factor* m_cholmodFactor;
448 double m_shiftOffset[2];
450 int m_factorizationIsOk;
476template<
typename _MatrixType,
int _UpLo = Lower>
480 using Base::m_cholmod;
484 typedef _MatrixType MatrixType;
498 m_cholmod.final_asis = 0;
499 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
500 m_cholmod.final_ll = 1;
527template<
typename _MatrixType,
int _UpLo = Lower>
531 using Base::m_cholmod;
535 typedef _MatrixType MatrixType;
549 m_cholmod.final_asis = 1;
550 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
576template<
typename _MatrixType,
int _UpLo = Lower>
580 using Base::m_cholmod;
584 typedef _MatrixType MatrixType;
598 m_cholmod.final_asis = 1;
599 m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
627template<
typename _MatrixType,
int _UpLo = Lower>
631 using Base::m_cholmod;
635 typedef _MatrixType MatrixType;
647 void setMode(CholmodMode mode)
652 m_cholmod.final_asis = 1;
653 m_cholmod.supernodal = CHOLMOD_AUTO;
655 case CholmodSimplicialLLt:
656 m_cholmod.final_asis = 0;
657 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
658 m_cholmod.final_ll = 1;
660 case CholmodSupernodalLLt:
661 m_cholmod.final_asis = 1;
662 m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
665 m_cholmod.final_asis = 1;
666 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
675 m_cholmod.final_asis = 1;
676 m_cholmod.supernodal = CHOLMOD_AUTO;
The base class for the direct Cholesky factorization of Cholmod.
Definition: CholmodSupport.h:215
cholmod_common & cholmod()
Definition: CholmodSupport.h:321
Derived & compute(const MatrixType &matrix)
Definition: CholmodSupport.h:273
Scalar logDeterminant() const
Definition: CholmodSupport.h:398
void analyzePattern(const MatrixType &matrix)
Definition: CholmodSupport.h:286
void factorize(const MatrixType &matrix)
Definition: CholmodSupport.h:308
Scalar determinant() const
Definition: CholmodSupport.h:391
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: CholmodSupport.h:266
Derived & setShift(const RealScalar &offset)
Definition: CholmodSupport.h:384
A general Cholesky factorization and solver based on Cholmod.
Definition: CholmodSupport.h:629
A simplicial direct Cholesky (LDLT) factorization and solver based on Cholmod.
Definition: CholmodSupport.h:529
A simplicial direct Cholesky (LLT) factorization and solver based on Cholmod.
Definition: CholmodSupport.h:478
A supernodal Cholesky (LLT) factorization and solver based on Cholmod.
Definition: CholmodSupport.h:578
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:63
Derived & derived()
Definition: EigenBase.h:46
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
Convenience specialization of Stride to specify only an inner stride See class Map for some examples.
Definition: Stride.h:96
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:96
Sparse matrix.
Definition: MappedSparseMatrix.h:34
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
static ConstMapType Map(const Scalar *data)
Definition: PlainObjectBase.h:644
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:283
A versatible sparse matrix representation.
Definition: SparseMatrix.h:98
Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix.
Definition: SparseSelfAdjointView.h:45
A base class for sparse solvers.
Definition: SparseSolverBase.h:68
ComputationInfo
Definition: Constants.h:440
@ Lower
Definition: Constants.h:209
@ Upper
Definition: Constants.h:211
@ NumericalIssue
Definition: Constants.h:444
@ Success
Definition: Constants.h:442
const unsigned int RowMajorBit
Definition: Constants.h:66
Namespace containing all symbols from the Eigen library.
Definition: Core:141
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_real_op< typename Derived::Scalar >, const Derived > real(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
cholmod_sparse viewAsCholmod(Ref< SparseMatrix< _Scalar, _Options, _StorageIndex > > mat)
Definition: CholmodSupport.h:58
MappedSparseMatrix< Scalar, Flags, StorageIndex > viewAsEigen(cholmod_sparse &cm)
Definition: CholmodSupport.h:158
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:233