CholmodSupport.h
1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
5//
6// This Source Code Form is subject to the terms of the Mozilla
7// Public License v. 2.0. If a copy of the MPL was not distributed
8// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9
10#ifndef EIGEN_CHOLMODSUPPORT_H
11#define EIGEN_CHOLMODSUPPORT_H
12
13namespace Eigen {
14
15namespace internal {
16
17template<typename Scalar, typename CholmodType>
18void cholmod_configure_matrix(CholmodType& mat)
19{
20 if (internal::is_same<Scalar,float>::value)
21 {
22 mat.xtype = CHOLMOD_REAL;
23 mat.dtype = CHOLMOD_SINGLE;
24 }
25 else if (internal::is_same<Scalar,double>::value)
26 {
27 mat.xtype = CHOLMOD_REAL;
28 mat.dtype = CHOLMOD_DOUBLE;
29 }
30 else if (internal::is_same<Scalar,std::complex<float> >::value)
31 {
32 mat.xtype = CHOLMOD_COMPLEX;
33 mat.dtype = CHOLMOD_SINGLE;
34 }
35 else if (internal::is_same<Scalar,std::complex<double> >::value)
36 {
37 mat.xtype = CHOLMOD_COMPLEX;
38 mat.dtype = CHOLMOD_DOUBLE;
39 }
40 else
41 {
42 eigen_assert(false && "Scalar type not supported by CHOLMOD");
43 }
44}
45
46} // namespace internal
47
51template<typename _Scalar, int _Options, typename _Index>
52cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_Index>& mat)
53{
54 typedef SparseMatrix<_Scalar,_Options,_Index> MatrixType;
55 cholmod_sparse res;
56 res.nzmax = mat.nonZeros();
57 res.nrow = mat.rows();;
58 res.ncol = mat.cols();
59 res.p = mat.outerIndexPtr();
60 res.i = mat.innerIndexPtr();
61 res.x = mat.valuePtr();
62 res.sorted = 1;
63 if(mat.isCompressed())
64 {
65 res.packed = 1;
66 }
67 else
68 {
69 res.packed = 0;
70 res.nz = mat.innerNonZeroPtr();
71 }
72
73 res.dtype = 0;
74 res.stype = -1;
75
76 if (internal::is_same<_Index,int>::value)
77 {
78 res.itype = CHOLMOD_INT;
79 }
80 else
81 {
82 eigen_assert(false && "Index type different than int is not supported yet");
83 }
84
85 // setup res.xtype
86 internal::cholmod_configure_matrix<_Scalar>(res);
87
88 res.stype = 0;
89
90 return res;
91}
92
93template<typename _Scalar, int _Options, typename _Index>
94const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat)
95{
96 cholmod_sparse res = viewAsCholmod(mat.const_cast_derived());
97 return res;
98}
99
102template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo>
103cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat)
104{
105 cholmod_sparse res = viewAsCholmod(mat.matrix().const_cast_derived());
106
107 if(UpLo==Upper) res.stype = 1;
108 if(UpLo==Lower) res.stype = -1;
109
110 return res;
111}
112
115template<typename Derived>
116cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat)
117{
118 EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
119 typedef typename Derived::Scalar Scalar;
120
121 cholmod_dense res;
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();
127 res.z = 0;
128
129 internal::cholmod_configure_matrix<Scalar>(res);
130
131 return res;
132}
133
136template<typename Scalar, int Flags, typename Index>
137MappedSparseMatrix<Scalar,Flags,Index> viewAsEigen(cholmod_sparse& cm)
138{
139 return MappedSparseMatrix<Scalar,Flags,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) );
142}
143
144enum CholmodMode {
145 CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
146};
147
148
154template<typename _MatrixType, int _UpLo, typename Derived>
155class CholmodBase : internal::noncopyable
156{
157 public:
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;
164
165 public:
166
167 CholmodBase()
168 : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
169 {
170 cholmod_start(&m_cholmod);
171 }
172
173 CholmodBase(const MatrixType& matrix)
174 : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
175 {
176 cholmod_start(&m_cholmod);
177 compute(matrix);
178 }
179
180 ~CholmodBase()
181 {
182 if(m_cholmodFactor)
183 cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
184 cholmod_finish(&m_cholmod);
185 }
186
187 inline Index cols() const { return m_cholmodFactor->n; }
188 inline Index rows() const { return m_cholmodFactor->n; }
189
190 Derived& derived() { return *static_cast<Derived*>(this); }
191 const Derived& derived() const { return *static_cast<const Derived*>(this); }
192
198 ComputationInfo info() const
199 {
200 eigen_assert(m_isInitialized && "Decomposition is not initialized.");
201 return m_info;
202 }
203
205 Derived& compute(const MatrixType& matrix)
206 {
207 analyzePattern(matrix);
208 factorize(matrix);
209 return derived();
210 }
211
216 template<typename Rhs>
217 inline const internal::solve_retval<CholmodBase, Rhs>
218 solve(const MatrixBase<Rhs>& b) const
219 {
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());
224 }
225
230 template<typename Rhs>
231 inline const internal::sparse_solve_retval<CholmodBase, Rhs>
232 solve(const SparseMatrixBase<Rhs>& b) const
233 {
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());
238 }
239
246 void analyzePattern(const MatrixType& matrix)
247 {
248 if(m_cholmodFactor)
249 {
250 cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
251 m_cholmodFactor = 0;
252 }
253 cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
254 m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
255
256 this->m_isInitialized = true;
257 this->m_info = Success;
258 m_analysisIsOk = true;
259 m_factorizationIsOk = false;
260 }
261
268 void factorize(const MatrixType& matrix)
269 {
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);
273
274 this->m_info = Success;
275 m_factorizationIsOk = true;
276 }
277
280 cholmod_common& cholmod() { return m_cholmod; }
281
282 #ifndef EIGEN_PARSED_BY_DOXYGEN
284 template<typename Rhs,typename Dest>
285 void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
286 {
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());
290
291 // note: cd stands for Cholmod Dense
292 cholmod_dense b_cd = viewAsCholmod(b.const_cast_derived());
293 cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
294 if(!x_cd)
295 {
296 this->m_info = NumericalIssue;
297 }
298 // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.)
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);
301 }
302
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
306 {
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());
310
311 // note: cs stands for Cholmod Sparse
312 cholmod_sparse b_cs = viewAsCholmod(b);
313 cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
314 if(!x_cs)
315 {
316 this->m_info = NumericalIssue;
317 }
318 // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.)
319 dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
320 cholmod_free_sparse(&x_cs, &m_cholmod);
321 }
322 #endif // EIGEN_PARSED_BY_DOXYGEN
323
324 template<typename Stream>
325 void dumpMemory(Stream& s)
326 {}
327
328 protected:
329 mutable cholmod_common m_cholmod;
330 cholmod_factor* m_cholmodFactor;
331 mutable ComputationInfo m_info;
332 bool m_isInitialized;
333 int m_factorizationIsOk;
334 int m_analysisIsOk;
335};
336
355template<typename _MatrixType, int _UpLo = Lower>
356class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> >
357{
358 typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT> Base;
359 using Base::m_cholmod;
360
361 public:
362
363 typedef _MatrixType MatrixType;
364
365 CholmodSimplicialLLT() : Base() { init(); }
366
367 CholmodSimplicialLLT(const MatrixType& matrix) : Base()
368 {
369 init();
370 compute(matrix);
371 }
372
373 ~CholmodSimplicialLLT() {}
374 protected:
375 void init()
376 {
377 m_cholmod.final_asis = 0;
378 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
379 m_cholmod.final_ll = 1;
380 }
381};
382
383
402template<typename _MatrixType, int _UpLo = Lower>
403class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> >
404{
405 typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT> Base;
406 using Base::m_cholmod;
407
408 public:
409
410 typedef _MatrixType MatrixType;
411
412 CholmodSimplicialLDLT() : Base() { init(); }
413
414 CholmodSimplicialLDLT(const MatrixType& matrix) : Base()
415 {
416 init();
417 compute(matrix);
418 }
419
420 ~CholmodSimplicialLDLT() {}
421 protected:
422 void init()
423 {
424 m_cholmod.final_asis = 1;
425 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
426 }
427};
428
447template<typename _MatrixType, int _UpLo = Lower>
448class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> >
449{
450 typedef CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT> Base;
451 using Base::m_cholmod;
452
453 public:
454
455 typedef _MatrixType MatrixType;
456
457 CholmodSupernodalLLT() : Base() { init(); }
458
459 CholmodSupernodalLLT(const MatrixType& matrix) : Base()
460 {
461 init();
462 compute(matrix);
463 }
464
465 ~CholmodSupernodalLLT() {}
466 protected:
467 void init()
468 {
469 m_cholmod.final_asis = 1;
470 m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
471 }
472};
473
494template<typename _MatrixType, int _UpLo = Lower>
495class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> >
496{
497 typedef CholmodBase<_MatrixType, _UpLo, CholmodDecomposition> Base;
498 using Base::m_cholmod;
499
500 public:
501
502 typedef _MatrixType MatrixType;
503
504 CholmodDecomposition() : Base() { init(); }
505
506 CholmodDecomposition(const MatrixType& matrix) : Base()
507 {
508 init();
509 compute(matrix);
510 }
511
512 ~CholmodDecomposition() {}
513
514 void setMode(CholmodMode mode)
515 {
516 switch(mode)
517 {
518 case CholmodAuto:
519 m_cholmod.final_asis = 1;
520 m_cholmod.supernodal = CHOLMOD_AUTO;
521 break;
522 case CholmodSimplicialLLt:
523 m_cholmod.final_asis = 0;
524 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
525 m_cholmod.final_ll = 1;
526 break;
527 case CholmodSupernodalLLt:
528 m_cholmod.final_asis = 1;
529 m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
530 break;
531 case CholmodLDLt:
532 m_cholmod.final_asis = 1;
533 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
534 break;
535 default:
536 break;
537 }
538 }
539 protected:
540 void init()
541 {
542 m_cholmod.final_asis = 1;
543 m_cholmod.supernodal = CHOLMOD_AUTO;
544 }
545};
546
547namespace internal {
548
549template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
550struct solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
551 : solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
552{
553 typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
554 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
555
556 template<typename Dest> void evalTo(Dest& dst) const
557 {
558 dec()._solve(rhs(),dst);
559 }
560};
561
562template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
563struct sparse_solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
564 : sparse_solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
565{
566 typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
567 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
568
569 template<typename Dest> void evalTo(Dest& dst) const
570 {
571 dec()._solve(rhs(),dst);
572 }
573};
574
575} // end namespace internal
576
577} // end namespace Eigen
578
579#endif // EIGEN_CHOLMODSUPPORT_H
Eigen::CholmodBase
The base class for the direct Cholesky factorization of Cholmod.
Definition CholmodSupport.h:156
Eigen::CholmodBase::compute
Derived & compute(const MatrixType &matrix)
Definition CholmodSupport.h:205
Eigen::CholmodBase::factorize
void factorize(const MatrixType &matrix)
Definition CholmodSupport.h:268
Eigen::CholmodBase::cholmod
cholmod_common & cholmod()
Definition CholmodSupport.h:280
Eigen::CholmodBase::solve
const internal::solve_retval< CholmodBase, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition CholmodSupport.h:218
Eigen::CholmodBase::analyzePattern
void analyzePattern(const MatrixType &matrix)
Definition CholmodSupport.h:246
Eigen::CholmodBase::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition CholmodSupport.h:198
Eigen::CholmodBase::solve
const internal::sparse_solve_retval< CholmodBase, Rhs > solve(const SparseMatrixBase< Rhs > &b) const
Definition CholmodSupport.h:232
Eigen::MappedSparseMatrix
Sparse matrix.
Definition MappedSparseMatrix.h:33
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition MatrixBase.h:50
Eigen::SparseMatrixBase
Base class of any sparse matrices or sparse expressions.
Definition SparseMatrixBase.h:27
Eigen::SparseMatrix
A versatible sparse matrix representation.
Definition SparseMatrix.h:87
Eigen::SparseSelfAdjointView
Pseudo expression to manipulate a triangular sparse matrix as a selfadjoint matrix.
Definition SparseSelfAdjointView.h:51
Eigen::ComputationInfo
ComputationInfo
Definition Constants.h:367
Eigen::NumericalIssue
@ NumericalIssue
Definition Constants.h:371
Eigen::Success
@ Success
Definition Constants.h:369
Eigen::Upper
@ Upper
Definition Constants.h:164
Eigen::Lower
@ Lower
Definition Constants.h:162
Eigen::RowMajorBit
const unsigned int RowMajorBit
Definition Constants.h:48
Eigen::internal
Definition LDLT.h:18
Eigen::EigenBase::rows
Index rows() const
Definition EigenBase.h:44
Eigen::EigenBase::derived
Derived & derived()
Definition EigenBase.h:34