10 #ifndef EIGEN_CHOLMODSUPPORT_H
11 #define EIGEN_CHOLMODSUPPORT_H
13 #include "./InternalHeaderCheck.h"
19 template<
typename Scalar>
struct cholmod_configure_matrix;
21 template<>
struct cholmod_configure_matrix<double> {
22 template<
typename CholmodType>
23 static void run(CholmodType& mat) {
24 mat.xtype = CHOLMOD_REAL;
25 mat.dtype = CHOLMOD_DOUBLE;
29 template<>
struct cholmod_configure_matrix<std::complex<double> > {
30 template<
typename CholmodType>
31 static void run(CholmodType& mat) {
32 mat.xtype = CHOLMOD_COMPLEX;
33 mat.dtype = CHOLMOD_DOUBLE;
59 template<
typename Scalar_,
int Options_,
typename StorageIndex_>
63 res.nzmax = mat.nonZeros();
64 res.nrow = mat.rows();
65 res.ncol = mat.cols();
66 res.p = mat.outerIndexPtr();
67 res.i = mat.innerIndexPtr();
68 res.x = mat.valuePtr();
71 if(mat.isCompressed())
79 res.nz = mat.innerNonZeroPtr();
85 if (internal::is_same<StorageIndex_,int>::value)
87 res.itype = CHOLMOD_INT;
89 else if (internal::is_same<StorageIndex_,SuiteSparse_long>::value)
91 res.itype = CHOLMOD_LONG;
95 eigen_assert(
false &&
"Index type not supported yet");
99 internal::cholmod_configure_matrix<Scalar_>::run(res);
106 template<
typename Scalar_,
int Options_,
typename Index_>
107 const cholmod_sparse
viewAsCholmod(
const SparseMatrix<Scalar_,Options_,Index_>& mat)
109 cholmod_sparse res =
viewAsCholmod(Ref<SparseMatrix<Scalar_,Options_,Index_> >(mat.const_cast_derived()));
113 template<
typename Scalar_,
int Options_,
typename Index_>
114 const cholmod_sparse
viewAsCholmod(
const SparseVector<Scalar_,Options_,Index_>& mat)
116 cholmod_sparse res =
viewAsCholmod(Ref<SparseMatrix<Scalar_,Options_,Index_> >(mat.const_cast_derived()));
122 template<
typename Scalar_,
int Options_,
typename Index_,
unsigned int UpLo>
127 if(UpLo==
Upper) res.stype = 1;
128 if(UpLo==
Lower) res.stype = -1;
138 template<
typename Derived>
141 EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&
RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
142 typedef typename Derived::Scalar Scalar;
145 res.nrow = mat.
rows();
146 res.ncol = mat.
cols();
147 res.nzmax = res.nrow * res.ncol;
148 res.d = Derived::IsVectorAtCompileTime ? mat.
derived().size() : mat.
derived().outerStride();
149 res.x = (
void*)(mat.
derived().data());
152 internal::cholmod_configure_matrix<Scalar>::run(res);
159 template<
typename Scalar,
int Flags,
typename StorageIndex>
163 (cm.nrow, cm.ncol,
static_cast<StorageIndex*
>(cm.p)[cm.ncol],
164 static_cast<StorageIndex*
>(cm.p),
static_cast<StorageIndex*
>(cm.i),
static_cast<Scalar*
>(cm.x) );
171 #define EIGEN_CHOLMOD_SPECIALIZE0(ret, name) \
172 template<typename StorageIndex_> inline ret cm_ ## name (cholmod_common &Common) { return cholmod_ ## name (&Common); } \
173 template<> inline ret cm_ ## name<SuiteSparse_long> (cholmod_common &Common) { return cholmod_l_ ## name (&Common); }
175 #define EIGEN_CHOLMOD_SPECIALIZE1(ret, name, t1, a1) \
176 template<typename StorageIndex_> inline ret cm_ ## name (t1& a1, cholmod_common &Common) { return cholmod_ ## name (&a1, &Common); } \
177 template<> inline ret cm_ ## name<SuiteSparse_long> (t1& a1, cholmod_common &Common) { return cholmod_l_ ## name (&a1, &Common); }
179 EIGEN_CHOLMOD_SPECIALIZE0(
int, start)
180 EIGEN_CHOLMOD_SPECIALIZE0(
int, finish)
182 EIGEN_CHOLMOD_SPECIALIZE1(
int, free_factor, cholmod_factor*, L)
183 EIGEN_CHOLMOD_SPECIALIZE1(
int, free_dense, cholmod_dense*, X)
184 EIGEN_CHOLMOD_SPECIALIZE1(
int, free_sparse, cholmod_sparse*, A)
186 EIGEN_CHOLMOD_SPECIALIZE1(cholmod_factor*, analyze, cholmod_sparse, A)
188 template<
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); }
189 template<>
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); }
191 template<
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); }
192 template<>
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); }
194 template<
typename StorageIndex_>
195 inline 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); }
197 inline 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); }
199 #undef EIGEN_CHOLMOD_SPECIALIZE0
200 #undef EIGEN_CHOLMOD_SPECIALIZE1
206 CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
215 template<
typename MatrixType_,
int UpLo_,
typename Derived>
221 using Base::m_isInitialized;
223 typedef MatrixType_ MatrixType;
224 enum { UpLo = UpLo_ };
225 typedef typename MatrixType::Scalar Scalar;
226 typedef typename MatrixType::RealScalar RealScalar;
227 typedef MatrixType CholMatrixType;
228 typedef typename MatrixType::StorageIndex StorageIndex;
230 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
231 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
237 : m_cholmodFactor(0), m_info(
Success), m_factorizationIsOk(
false), m_analysisIsOk(
false)
239 EIGEN_STATIC_ASSERT((internal::is_same<double,RealScalar>::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY);
240 m_shiftOffset[0] = m_shiftOffset[1] = 0.0;
241 internal::cm_start<StorageIndex>(m_cholmod);
245 : m_cholmodFactor(0), m_info(
Success), m_factorizationIsOk(
false), m_analysisIsOk(
false)
247 EIGEN_STATIC_ASSERT((internal::is_same<double,RealScalar>::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY);
248 m_shiftOffset[0] = m_shiftOffset[1] = 0.0;
249 internal::cm_start<StorageIndex>(m_cholmod);
256 internal::cm_free_factor<StorageIndex>(m_cholmodFactor, m_cholmod);
257 internal::cm_finish<StorageIndex>(m_cholmod);
260 inline StorageIndex cols()
const {
return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
261 inline StorageIndex rows()
const {
return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
270 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
292 internal::cm_free_factor<StorageIndex>(m_cholmodFactor, m_cholmod);
295 cholmod_sparse A =
viewAsCholmod(matrix.template selfadjointView<UpLo>());
296 m_cholmodFactor = internal::cm_analyze<StorageIndex>(A, m_cholmod);
298 this->m_isInitialized =
true;
300 m_analysisIsOk =
true;
301 m_factorizationIsOk =
false;
312 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
313 cholmod_sparse A =
viewAsCholmod(matrix.template selfadjointView<UpLo>());
314 internal::cm_factorize_p<StorageIndex>(&A, m_shiftOffset, 0, 0, m_cholmodFactor, m_cholmod);
318 m_factorizationIsOk =
true;
323 cholmod_common&
cholmod() {
return m_cholmod; }
325 #ifndef EIGEN_PARSED_BY_DOXYGEN
327 template<
typename Rhs,
typename Dest>
330 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
331 const Index size = m_cholmodFactor->n;
332 EIGEN_UNUSED_VARIABLE(size);
333 eigen_assert(size==b.
rows());
339 cholmod_dense* x_cd = internal::cm_solve<StorageIndex>(CHOLMOD_A, *m_cholmodFactor, b_cd, m_cholmod);
348 internal::cm_free_dense<StorageIndex>(x_cd, m_cholmod);
352 template<
typename RhsDerived,
typename DestDerived>
353 void _solve_impl(
const SparseMatrixBase<RhsDerived> &b, SparseMatrixBase<DestDerived> &dest)
const
355 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
356 const Index size = m_cholmodFactor->n;
357 EIGEN_UNUSED_VARIABLE(size);
358 eigen_assert(size==b.rows());
361 Ref<SparseMatrix<typename RhsDerived::Scalar,ColMajor,typename RhsDerived::StorageIndex> > b_ref(b.const_cast_derived());
363 cholmod_sparse* x_cs = internal::cm_spsolve<StorageIndex>(CHOLMOD_A, *m_cholmodFactor, b_cs, m_cholmod);
371 dest.derived() = viewAsEigen<typename DestDerived::Scalar,ColMajor,typename DestDerived::StorageIndex>(*x_cs);
372 internal::cm_free_sparse<StorageIndex>(x_cs, m_cholmod);
388 m_shiftOffset[0] = double(offset);
404 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
406 RealScalar logDet = 0;
407 Scalar *x =
static_cast<Scalar*
>(m_cholmodFactor->x);
408 if (m_cholmodFactor->is_super)
414 StorageIndex *super =
static_cast<StorageIndex*
>(m_cholmodFactor->super);
416 StorageIndex *pi =
static_cast<StorageIndex*
>(m_cholmodFactor->pi);
418 StorageIndex *px =
static_cast<StorageIndex*
>(m_cholmodFactor->px);
420 Index nb_super_nodes = m_cholmodFactor->nsuper;
421 for (
Index k=0; k < nb_super_nodes; ++k)
423 StorageIndex ncols = super[k + 1] - super[k];
424 StorageIndex nrows = pi[k + 1] - pi[k];
427 logDet += sk.real().log().sum();
433 StorageIndex *p =
static_cast<StorageIndex*
>(m_cholmodFactor->p);
434 Index size = m_cholmodFactor->n;
435 for (
Index k=0; k<size; ++k)
436 logDet +=
log(
real( x[p[k]] ));
438 if (m_cholmodFactor->is_ll)
443 template<
typename Stream>
444 void dumpMemory(Stream& )
448 mutable cholmod_common m_cholmod;
449 cholmod_factor* m_cholmodFactor;
450 double m_shiftOffset[2];
452 int m_factorizationIsOk;
478 template<
typename MatrixType_,
int UpLo_ = Lower>
482 using Base::m_cholmod;
486 typedef MatrixType_ MatrixType;
500 m_cholmod.final_asis = 0;
501 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
502 m_cholmod.final_ll = 1;
529 template<
typename MatrixType_,
int UpLo_ = Lower>
533 using Base::m_cholmod;
537 typedef MatrixType_ MatrixType;
551 m_cholmod.final_asis = 1;
552 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
578 template<
typename MatrixType_,
int UpLo_ = Lower>
582 using Base::m_cholmod;
586 typedef MatrixType_ MatrixType;
600 m_cholmod.final_asis = 1;
601 m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
629 template<
typename MatrixType_,
int UpLo_ = Lower>
633 using Base::m_cholmod;
637 typedef MatrixType_ MatrixType;
649 void setMode(CholmodMode mode)
654 m_cholmod.final_asis = 1;
655 m_cholmod.supernodal = CHOLMOD_AUTO;
657 case CholmodSimplicialLLt:
658 m_cholmod.final_asis = 0;
659 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
660 m_cholmod.final_ll = 1;
662 case CholmodSupernodalLLt:
663 m_cholmod.final_asis = 1;
664 m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
667 m_cholmod.final_asis = 1;
668 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
677 m_cholmod.final_asis = 1;
678 m_cholmod.supernodal = CHOLMOD_AUTO;
The base class for the direct Cholesky factorization of Cholmod.
Definition: CholmodSupport.h:217
void factorize(const MatrixType &matrix)
Definition: CholmodSupport.h:310
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: CholmodSupport.h:268
cholmod_common & cholmod()
Definition: CholmodSupport.h:323
Scalar determinant() const
Definition: CholmodSupport.h:393
Derived & compute(const MatrixType &matrix)
Definition: CholmodSupport.h:275
Scalar logDeterminant() const
Definition: CholmodSupport.h:400
Derived & setShift(const RealScalar &offset)
Definition: CholmodSupport.h:386
void analyzePattern(const MatrixType &matrix)
Definition: CholmodSupport.h:288
A general Cholesky factorization and solver based on Cholmod.
Definition: CholmodSupport.h:631
A simplicial direct Cholesky (LDLT) factorization and solver based on Cholmod.
Definition: CholmodSupport.h:531
A simplicial direct Cholesky (LLT) factorization and solver based on Cholmod.
Definition: CholmodSupport.h:480
A supernodal Cholesky (LLT) factorization and solver based on Cholmod.
Definition: CholmodSupport.h:580
Derived & derived()
Definition: EigenBase.h:48
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:65
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:62
Convenience specialization of Stride to specify only an inner stride See class Map for some examples.
Definition: Stride.h:102
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:98
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
static ConstMapType Map(const Scalar *data)
Definition: PlainObjectBase.h:642
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:285
A versatible sparse matrix representation.
Definition: SparseMatrix.h:100
Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix.
Definition: SparseSelfAdjointView.h:47
A base class for sparse solvers.
Definition: SparseSolverBase.h:70
ComputationInfo
Definition: Constants.h:442
@ Lower
Definition: Constants.h:211
@ Upper
Definition: Constants.h:213
@ NumericalIssue
Definition: Constants.h:446
@ Success
Definition: Constants.h:444
const unsigned int RowMajorBit
Definition: Constants.h:68
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
cholmod_sparse viewAsCholmod(Ref< SparseMatrix< Scalar_, Options_, StorageIndex_ > > mat)
Definition: CholmodSupport.h:60
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_real_op< typename Derived::Scalar >, const Derived > real(const Eigen::ArrayBase< Derived > &x)
Map< SparseMatrix< Scalar, Flags, StorageIndex > > viewAsEigen(cholmod_sparse &cm)
Definition: CholmodSupport.h:160
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_exp_op< typename Derived::Scalar >, const Derived > exp(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_log_op< typename Derived::Scalar >, const Derived > log(const Eigen::ArrayBase< Derived > &x)
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231