Eigen  3.2.10
 
Loading...
Searching...
No Matches
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> struct cholmod_configure_matrix;
18
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;
24 }
25};
26
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;
32 }
33};
34
35// Other scalar types are not yet suppotred by Cholmod
36// template<> struct cholmod_configure_matrix<float> {
37// template<typename CholmodType>
38// static void run(CholmodType& mat) {
39// mat.xtype = CHOLMOD_REAL;
40// mat.dtype = CHOLMOD_SINGLE;
41// }
42// };
43//
44// template<> struct cholmod_configure_matrix<std::complex<float> > {
45// template<typename CholmodType>
46// static void run(CholmodType& mat) {
47// mat.xtype = CHOLMOD_COMPLEX;
48// mat.dtype = CHOLMOD_SINGLE;
49// }
50// };
51
52} // namespace internal
53
57template<typename _Scalar, int _Options, typename _Index>
58cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_Index>& mat)
59{
60 cholmod_sparse res;
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();
67 res.z = 0;
68 res.sorted = 1;
69 if(mat.isCompressed())
70 {
71 res.packed = 1;
72 res.nz = 0;
73 }
74 else
75 {
76 res.packed = 0;
77 res.nz = mat.innerNonZeroPtr();
78 }
79
80 res.dtype = 0;
81 res.stype = -1;
82
83 if (internal::is_same<_Index,int>::value)
84 {
85 res.itype = CHOLMOD_INT;
86 }
87 else if (internal::is_same<_Index,SuiteSparse_long>::value)
88 {
89 res.itype = CHOLMOD_LONG;
90 }
91 else
92 {
93 eigen_assert(false && "Index type not supported yet");
94 }
95
96 // setup res.xtype
97 internal::cholmod_configure_matrix<_Scalar>::run(res);
98
99 res.stype = 0;
100
101 return res;
102}
103
104template<typename _Scalar, int _Options, typename _Index>
105const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat)
106{
107 cholmod_sparse res = viewAsCholmod(mat.const_cast_derived());
108 return res;
109}
110
113template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo>
114cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat)
115{
116 cholmod_sparse res = viewAsCholmod(mat.matrix().const_cast_derived());
117
118 if(UpLo==Upper) res.stype = 1;
119 if(UpLo==Lower) res.stype = -1;
120
121 return res;
122}
123
126template<typename Derived>
127cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat)
128{
129 EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
130 typedef typename Derived::Scalar Scalar;
131
132 cholmod_dense res;
133 res.nrow = mat.rows();
134 res.ncol = mat.cols();
135 res.nzmax = res.nrow * res.ncol;
136 res.d = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
137 res.x = (void*)(mat.derived().data());
138 res.z = 0;
139
140 internal::cholmod_configure_matrix<Scalar>::run(res);
141
142 return res;
143}
144
147template<typename Scalar, int Flags, typename Index>
148MappedSparseMatrix<Scalar,Flags,Index> viewAsEigen(cholmod_sparse& cm)
149{
150 return MappedSparseMatrix<Scalar,Flags,Index>
151 (cm.nrow, cm.ncol, static_cast<Index*>(cm.p)[cm.ncol],
152 static_cast<Index*>(cm.p), static_cast<Index*>(cm.i),static_cast<Scalar*>(cm.x) );
153}
154
155enum CholmodMode {
156 CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
157};
158
159
165template<typename _MatrixType, int _UpLo, typename Derived>
166class CholmodBase : internal::noncopyable
167{
168 public:
169 typedef _MatrixType MatrixType;
170 enum { UpLo = _UpLo };
171 typedef typename MatrixType::Scalar Scalar;
172 typedef typename MatrixType::RealScalar RealScalar;
173 typedef MatrixType CholMatrixType;
174 typedef typename MatrixType::Index Index;
175
176 public:
177
178 CholmodBase()
179 : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
180 {
181 EIGEN_STATIC_ASSERT((internal::is_same<double,RealScalar>::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY);
182 m_shiftOffset[0] = m_shiftOffset[1] = 0.0;
183 cholmod_start(&m_cholmod);
184 }
185
186 CholmodBase(const MatrixType& matrix)
187 : m_cholmodFactor(0), m_info(Success), m_isInitialized(false)
188 {
189 EIGEN_STATIC_ASSERT((internal::is_same<double,RealScalar>::value), CHOLMOD_SUPPORTS_DOUBLE_PRECISION_ONLY);
190 m_shiftOffset[0] = m_shiftOffset[1] = 0.0;
191 cholmod_start(&m_cholmod);
192 compute(matrix);
193 }
194
195 ~CholmodBase()
196 {
197 if(m_cholmodFactor)
198 cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
199 cholmod_finish(&m_cholmod);
200 }
201
202 inline Index cols() const { return m_cholmodFactor->n; }
203 inline Index rows() const { return m_cholmodFactor->n; }
204
205 Derived& derived() { return *static_cast<Derived*>(this); }
206 const Derived& derived() const { return *static_cast<const Derived*>(this); }
207
213 ComputationInfo info() const
214 {
215 eigen_assert(m_isInitialized && "Decomposition is not initialized.");
216 return m_info;
217 }
218
220 Derived& compute(const MatrixType& matrix)
221 {
222 analyzePattern(matrix);
223 factorize(matrix);
224 return derived();
225 }
226
231 template<typename Rhs>
232 inline const internal::solve_retval<CholmodBase, Rhs>
233 solve(const MatrixBase<Rhs>& b) const
234 {
235 eigen_assert(m_isInitialized && "LLT is not initialized.");
236 eigen_assert(rows()==b.rows()
237 && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
238 return internal::solve_retval<CholmodBase, Rhs>(*this, b.derived());
239 }
240
245 template<typename Rhs>
246 inline const internal::sparse_solve_retval<CholmodBase, Rhs>
247 solve(const SparseMatrixBase<Rhs>& b) const
248 {
249 eigen_assert(m_isInitialized && "LLT is not initialized.");
250 eigen_assert(rows()==b.rows()
251 && "CholmodDecomposition::solve(): invalid number of rows of the right hand side matrix b");
252 return internal::sparse_solve_retval<CholmodBase, Rhs>(*this, b.derived());
253 }
254
261 void analyzePattern(const MatrixType& matrix)
262 {
263 if(m_cholmodFactor)
264 {
265 cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
266 m_cholmodFactor = 0;
267 }
268 cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
269 m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
270
271 this->m_isInitialized = true;
272 this->m_info = Success;
273 m_analysisIsOk = true;
274 m_factorizationIsOk = false;
275 }
276
283 void factorize(const MatrixType& matrix)
284 {
285 eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
286 cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
287 cholmod_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod);
288
289 // If the factorization failed, minor is the column at which it did. On success minor == n.
290 this->m_info = (m_cholmodFactor->minor == m_cholmodFactor->n ? Success : NumericalIssue);
291 m_factorizationIsOk = true;
292 }
293
296 cholmod_common& cholmod() { return m_cholmod; }
297
298 #ifndef EIGEN_PARSED_BY_DOXYGEN
300 template<typename Rhs,typename Dest>
301 void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
302 {
303 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
304 const Index size = m_cholmodFactor->n;
305 EIGEN_UNUSED_VARIABLE(size);
306 eigen_assert(size==b.rows());
307
308 // note: cd stands for Cholmod Dense
309 Rhs& b_ref(b.const_cast_derived());
310 cholmod_dense b_cd = viewAsCholmod(b_ref);
311 cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
312 if(!x_cd)
313 {
314 this->m_info = NumericalIssue;
315 }
316 // TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
317 dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
318 cholmod_free_dense(&x_cd, &m_cholmod);
319 }
320
322 template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex>
323 void _solve(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
324 {
325 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
326 const Index size = m_cholmodFactor->n;
327 EIGEN_UNUSED_VARIABLE(size);
328 eigen_assert(size==b.rows());
329
330 // note: cs stands for Cholmod Sparse
331 cholmod_sparse b_cs = viewAsCholmod(b);
332 cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
333 if(!x_cs)
334 {
335 this->m_info = NumericalIssue;
336 }
337 // TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
338 dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
339 cholmod_free_sparse(&x_cs, &m_cholmod);
340 }
341 #endif // EIGEN_PARSED_BY_DOXYGEN
342
343
353 Derived& setShift(const RealScalar& offset)
354 {
355 m_shiftOffset[0] = double(offset);
356 return derived();
357 }
358
359 template<typename Stream>
360 void dumpMemory(Stream& /*s*/)
361 {}
362
363 protected:
364 mutable cholmod_common m_cholmod;
365 cholmod_factor* m_cholmodFactor;
366 double m_shiftOffset[2];
367 mutable ComputationInfo m_info;
368 bool m_isInitialized;
369 int m_factorizationIsOk;
370 int m_analysisIsOk;
371};
372
393template<typename _MatrixType, int _UpLo = Lower>
394class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> >
395{
396 typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT> Base;
397 using Base::m_cholmod;
398
399 public:
400
401 typedef _MatrixType MatrixType;
402
403 CholmodSimplicialLLT() : Base() { init(); }
404
405 CholmodSimplicialLLT(const MatrixType& matrix) : Base()
406 {
407 init();
408 Base::compute(matrix);
409 }
410
411 ~CholmodSimplicialLLT() {}
412 protected:
413 void init()
414 {
415 m_cholmod.final_asis = 0;
416 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
417 m_cholmod.final_ll = 1;
418 }
419};
420
421
442template<typename _MatrixType, int _UpLo = Lower>
443class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> >
444{
445 typedef CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT> Base;
446 using Base::m_cholmod;
447
448 public:
449
450 typedef _MatrixType MatrixType;
451
452 CholmodSimplicialLDLT() : Base() { init(); }
453
454 CholmodSimplicialLDLT(const MatrixType& matrix) : Base()
455 {
456 init();
457 Base::compute(matrix);
458 }
459
460 ~CholmodSimplicialLDLT() {}
461 protected:
462 void init()
463 {
464 m_cholmod.final_asis = 1;
465 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
466 }
467};
468
489template<typename _MatrixType, int _UpLo = Lower>
490class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> >
491{
492 typedef CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT> Base;
493 using Base::m_cholmod;
494
495 public:
496
497 typedef _MatrixType MatrixType;
498
499 CholmodSupernodalLLT() : Base() { init(); }
500
501 CholmodSupernodalLLT(const MatrixType& matrix) : Base()
502 {
503 init();
504 Base::compute(matrix);
505 }
506
507 ~CholmodSupernodalLLT() {}
508 protected:
509 void init()
510 {
511 m_cholmod.final_asis = 1;
512 m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
513 }
514};
515
538template<typename _MatrixType, int _UpLo = Lower>
539class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> >
540{
541 typedef CholmodBase<_MatrixType, _UpLo, CholmodDecomposition> Base;
542 using Base::m_cholmod;
543
544 public:
545
546 typedef _MatrixType MatrixType;
547
548 CholmodDecomposition() : Base() { init(); }
549
550 CholmodDecomposition(const MatrixType& matrix) : Base()
551 {
552 init();
553 Base::compute(matrix);
554 }
555
556 ~CholmodDecomposition() {}
557
558 void setMode(CholmodMode mode)
559 {
560 switch(mode)
561 {
562 case CholmodAuto:
563 m_cholmod.final_asis = 1;
564 m_cholmod.supernodal = CHOLMOD_AUTO;
565 break;
566 case CholmodSimplicialLLt:
567 m_cholmod.final_asis = 0;
568 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
569 m_cholmod.final_ll = 1;
570 break;
571 case CholmodSupernodalLLt:
572 m_cholmod.final_asis = 1;
573 m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
574 break;
575 case CholmodLDLt:
576 m_cholmod.final_asis = 1;
577 m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
578 break;
579 default:
580 break;
581 }
582 }
583 protected:
584 void init()
585 {
586 m_cholmod.final_asis = 1;
587 m_cholmod.supernodal = CHOLMOD_AUTO;
588 }
589};
590
591namespace internal {
592
593template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
594struct solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
595 : solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
596{
597 typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
598 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
599
600 template<typename Dest> void evalTo(Dest& dst) const
601 {
602 dec()._solve(rhs(),dst);
603 }
604};
605
606template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
607struct sparse_solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
608 : sparse_solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
609{
610 typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
611 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
612
613 template<typename Dest> void evalTo(Dest& dst) const
614 {
615 dec()._solve(rhs(),dst);
616 }
617};
618
619} // end namespace internal
620
621} // end namespace Eigen
622
623#endif // EIGEN_CHOLMODSUPPORT_H
Eigen::CholmodBase
The base class for the direct Cholesky factorization of Cholmod.
Definition CholmodSupport.h:167
Eigen::CholmodBase::compute
Derived & compute(const MatrixType &matrix)
Definition CholmodSupport.h:220
Eigen::CholmodBase::factorize
void factorize(const MatrixType &matrix)
Definition CholmodSupport.h:283
Eigen::CholmodBase::cholmod
cholmod_common & cholmod()
Definition CholmodSupport.h:296
Eigen::CholmodBase::setShift
Derived & setShift(const RealScalar &offset)
Definition CholmodSupport.h:353
Eigen::CholmodBase::solve
const internal::solve_retval< CholmodBase, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition CholmodSupport.h:233
Eigen::CholmodBase::analyzePattern
void analyzePattern(const MatrixType &matrix)
Definition CholmodSupport.h:261
Eigen::CholmodBase::info
ComputationInfo info() const
Reports whether previous computation was successful.
Definition CholmodSupport.h:213
Eigen::CholmodBase::solve
const internal::sparse_solve_retval< CholmodBase, Rhs > solve(const SparseMatrixBase< Rhs > &b) const
Definition CholmodSupport.h:247
Eigen::MappedSparseMatrix
Sparse matrix.
Definition MappedSparseMatrix.h:33
Eigen::MatrixBase
Base class for all dense matrices, vectors, and expressions.
Definition MatrixBase.h:50
Eigen::PlainObjectBase< Matrix< _Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols > >::Map
static ConstMapType Map(const Scalar *data)
Definition PlainObjectBase.h:505
Eigen::SparseMatrixBase
Base class of any sparse matrices or sparse expressions.
Definition SparseMatrixBase.h:34
Eigen::SparseMatrixBase::rows
Index rows() const
Definition SparseMatrixBase.h:160
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:374
Eigen::NumericalIssue
@ NumericalIssue
Definition Constants.h:378
Eigen::Success
@ Success
Definition Constants.h:376
Eigen::Upper
@ Upper
Definition Constants.h:169
Eigen::Lower
@ Lower
Definition Constants.h:167
Eigen::RowMajorBit
const unsigned int RowMajorBit
Definition Constants.h:53
Eigen::internal
Definition LDLT.h:18
Eigen::EigenBase::derived
Derived & derived()
Definition EigenBase.h:34