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 _StorageIndex>
56 res.nrow = mat.
rows();;
57 res.ncol = mat.
cols();
77 if (internal::is_same<_StorageIndex,int>::value)
79 res.itype = CHOLMOD_INT;
81 else if (internal::is_same<_StorageIndex,long>::value)
83 res.itype = CHOLMOD_LONG;
87 eigen_assert(
false &&
"Index type not supported yet");
91 internal::cholmod_configure_matrix<_Scalar>(res);
98 template<
typename _Scalar,
int _Options,
typename _Index>
101 cholmod_sparse res =
viewAsCholmod(mat.const_cast_derived());
107 template<
typename _Scalar,
int _Options,
typename _Index,
unsigned int UpLo>
110 cholmod_sparse res =
viewAsCholmod(mat.matrix().const_cast_derived());
112 if(UpLo==
Upper) res.stype = 1;
113 if(UpLo==
Lower) res.stype = -1;
120 template<
typename Derived>
123 EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&
RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
124 typedef typename Derived::Scalar Scalar;
127 res.nrow = mat.
rows();
128 res.ncol = mat.
cols();
129 res.nzmax = res.nrow * res.ncol;
130 res.d = Derived::IsVectorAtCompileTime ? mat.
derived().size() : mat.
derived().outerStride();
131 res.x = (
void*)(mat.
derived().data());
134 internal::cholmod_configure_matrix<Scalar>(res);
141 template<
typename Scalar,
int Flags,
typename StorageIndex>
145 (cm.nrow, cm.ncol,
static_cast<StorageIndex*
>(cm.p)[cm.ncol],
146 static_cast<StorageIndex*>(cm.p),
static_cast<StorageIndex*
>(cm.i),static_cast<Scalar*>(cm.x) );
150 CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
159 template<
typename _MatrixType,
int _UpLo,
typename Derived>
165 using Base::m_isInitialized;
167 typedef _MatrixType MatrixType;
168 enum { UpLo = _UpLo };
169 typedef typename MatrixType::Scalar Scalar;
170 typedef typename MatrixType::RealScalar RealScalar;
171 typedef MatrixType CholMatrixType;
172 typedef typename MatrixType::StorageIndex StorageIndex;
174 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
175 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
181 : m_cholmodFactor(0), m_info(
Success), m_factorizationIsOk(
false), m_analysisIsOk(
false)
183 m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
184 cholmod_start(&m_cholmod);
188 : m_cholmodFactor(0), m_info(
Success), m_factorizationIsOk(
false), m_analysisIsOk(
false)
190 m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
191 cholmod_start(&m_cholmod);
198 cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
199 cholmod_finish(&m_cholmod);
202 inline StorageIndex cols()
const {
return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
203 inline StorageIndex rows()
const {
return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
212 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
219 analyzePattern(matrix);
234 cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
237 cholmod_sparse A =
viewAsCholmod(matrix.template selfadjointView<UpLo>());
238 m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
240 this->m_isInitialized =
true;
242 m_analysisIsOk =
true;
243 m_factorizationIsOk =
false;
254 eigen_assert(m_analysisIsOk &&
"You must first call analyzePattern()");
255 cholmod_sparse A =
viewAsCholmod(matrix.template selfadjointView<UpLo>());
256 cholmod_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod);
260 m_factorizationIsOk =
true;
265 cholmod_common&
cholmod() {
return m_cholmod; }
267 #ifndef EIGEN_PARSED_BY_DOXYGEN 269 template<
typename Rhs,
typename Dest>
272 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
273 const Index size = m_cholmodFactor->n;
274 EIGEN_UNUSED_VARIABLE(size);
275 eigen_assert(size==b.
rows());
281 cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
289 cholmod_free_dense(&x_cd, &m_cholmod);
293 template<
typename RhsScalar,
int RhsOptions,
typename RhsIndex,
typename DestScalar,
int DestOptions,
typename DestIndex>
296 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
297 const Index size = m_cholmodFactor->n;
298 EIGEN_UNUSED_VARIABLE(size);
299 eigen_assert(size==b.
rows());
303 cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
310 dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
311 cholmod_free_sparse(&x_cs, &m_cholmod);
313 #endif // EIGEN_PARSED_BY_DOXYGEN 327 m_shiftOffset[0] = offset;
335 return exp(logDeterminant());
343 eigen_assert(m_factorizationIsOk &&
"The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
345 RealScalar logDet = 0;
346 Scalar *x =
static_cast<Scalar*
>(m_cholmodFactor->x);
347 if (m_cholmodFactor->is_super)
353 StorageIndex *super =
static_cast<StorageIndex*
>(m_cholmodFactor->super);
355 StorageIndex *pi =
static_cast<StorageIndex*
>(m_cholmodFactor->pi);
357 StorageIndex *px =
static_cast<StorageIndex*
>(m_cholmodFactor->px);
359 Index nb_super_nodes = m_cholmodFactor->nsuper;
360 for (
Index k=0; k < nb_super_nodes; ++k)
362 StorageIndex ncols = super[k + 1] - super[k];
363 StorageIndex nrows = pi[k + 1] - pi[k];
366 logDet += sk.real().log().sum();
372 StorageIndex *p =
static_cast<StorageIndex*
>(m_cholmodFactor->p);
373 Index size = m_cholmodFactor->n;
374 for (
Index k=0; k<size; ++k)
375 logDet +=
log(
real( x[p[k]] ));
377 if (m_cholmodFactor->is_ll)
382 template<
typename Stream>
383 void dumpMemory(Stream& )
387 mutable cholmod_common m_cholmod;
388 cholmod_factor* m_cholmodFactor;
389 RealScalar m_shiftOffset[2];
391 int m_factorizationIsOk;
415 template<
typename _MatrixType,
int _UpLo = Lower>
419 using Base::m_cholmod;
423 typedef _MatrixType MatrixType;
430 this->compute(matrix);
437 m_cholmod.final_asis = 0;
438 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
439 m_cholmod.final_ll = 1;
464 template<
typename _MatrixType,
int _UpLo = Lower>
468 using Base::m_cholmod;
472 typedef _MatrixType MatrixType;
479 this->compute(matrix);
486 m_cholmod.final_asis = 1;
487 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
511 template<
typename _MatrixType,
int _UpLo = Lower>
515 using Base::m_cholmod;
519 typedef _MatrixType MatrixType;
526 this->compute(matrix);
533 m_cholmod.final_asis = 1;
534 m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
560 template<
typename _MatrixType,
int _UpLo = Lower>
564 using Base::m_cholmod;
568 typedef _MatrixType MatrixType;
575 this->compute(matrix);
580 void setMode(CholmodMode mode)
585 m_cholmod.final_asis = 1;
586 m_cholmod.supernodal = CHOLMOD_AUTO;
588 case CholmodSimplicialLLt:
589 m_cholmod.final_asis = 0;
590 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
591 m_cholmod.final_ll = 1;
593 case CholmodSupernodalLLt:
594 m_cholmod.final_asis = 1;
595 m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
598 m_cholmod.final_asis = 1;
599 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
608 m_cholmod.final_asis = 1;
609 m_cholmod.supernodal = CHOLMOD_AUTO;
615 #endif // EIGEN_CHOLMODSUPPORT_H bool isCompressed() const
Definition: SparseCompressedBase.h:107
Scalar logDeterminant() const
Definition: CholmodSupport.h:339
const Scalar * valuePtr() const
Definition: SparseMatrix.h:143
A versatible sparse matrix representation.
Definition: SparseMatrix.h:92
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:88
void factorize(const MatrixType &matrix)
Definition: CholmodSupport.h:252
A base class for sparse solvers.
Definition: SparseSolverBase.h:53
MappedSparseMatrix< Scalar, Flags, StorageIndex > viewAsEigen(cholmod_sparse &cm)
Definition: CholmodSupport.h:142
Namespace containing all symbols from the Eigen library.
Definition: Core:271
Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix.
Definition: SparseSelfAdjointView.h:43
Index cols() const
Definition: SparseMatrix.h:133
A supernodal Cholesky (LLT) factorization and solver based on Cholmod.
Definition: CholmodSupport.h:512
Derived & derived()
Definition: EigenBase.h:44
Index nonZeros() const
Definition: SparseCompressedBase.h:56
const unsigned int RowMajorBit
Definition: Constants.h:61
cholmod_common & cholmod()
Definition: CholmodSupport.h:265
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: CholmodSupport.h:210
Scalar determinant() const
Definition: CholmodSupport.h:332
Definition: Constants.h:204
Index rows() const
Definition: EigenBase.h:58
const StorageIndex * innerNonZeroPtr() const
Definition: SparseMatrix.h:170
Definition: Constants.h:434
Convenience specialization of Stride to specify only an inner stride See class Map for some examples...
Definition: Stride.h:90
void analyzePattern(const MatrixType &matrix)
Definition: CholmodSupport.h:230
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: XprHelper.h:35
Index rows() const
Definition: SparseMatrix.h:131
A general Cholesky factorization and solver based on Cholmod.
Definition: CholmodSupport.h:561
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_real_op< typename Derived::Scalar >, const Derived > real(const Eigen::ArrayBase< Derived > &x)
Definition: Constants.h:206
The base class for the direct Cholesky factorization of Cholmod.
Definition: CholmodSupport.h:160
Definition: Constants.h:432
A matrix or vector expression mapping an existing expression.
Definition: Ref.h:184
A simplicial direct Cholesky (LLT) factorization and solver based on Cholmod.
Definition: CholmodSupport.h:416
const StorageIndex * innerIndexPtr() const
Definition: SparseMatrix.h:152
Definition: Eigen_Colamd.h:50
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_log_op< typename Derived::Scalar >, const Derived > log(const Eigen::ArrayBase< Derived > &x)
A simplicial direct Cholesky (LDLT) factorization and solver based on Cholmod.
Definition: CholmodSupport.h:465
const StorageIndex * outerIndexPtr() const
Definition: SparseMatrix.h:161
Derived & compute(const MatrixType &matrix)
Definition: CholmodSupport.h:217
Derived & setShift(const RealScalar &offset)
Definition: CholmodSupport.h:325
cholmod_sparse viewAsCholmod(SparseMatrix< _Scalar, _Options, _StorageIndex > &mat)
Definition: CholmodSupport.h:52
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:178
Index cols() const
Definition: EigenBase.h:61
ComputationInfo
Definition: Constants.h:430
Sparse matrix.
Definition: MappedSparseMatrix.h:32
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_exp_op< typename Derived::Scalar >, const Derived > exp(const Eigen::ArrayBase< Derived > &x)