Eigen  3.2.91
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 
13 namespace Eigen {
14 
15 namespace internal {
16 
17 template<typename Scalar, typename CholmodType>
18 void 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 
51 template<typename _Scalar, int _Options, typename _StorageIndex>
52 cholmod_sparse viewAsCholmod(SparseMatrix<_Scalar,_Options,_StorageIndex>& mat)
53 {
54  cholmod_sparse res;
55  res.nzmax = mat.nonZeros();
56  res.nrow = mat.rows();;
57  res.ncol = mat.cols();
58  res.p = mat.outerIndexPtr();
59  res.i = mat.innerIndexPtr();
60  res.x = mat.valuePtr();
61  res.z = 0;
62  res.sorted = 1;
63  if(mat.isCompressed())
64  {
65  res.packed = 1;
66  res.nz = 0;
67  }
68  else
69  {
70  res.packed = 0;
71  res.nz = mat.innerNonZeroPtr();
72  }
73 
74  res.dtype = 0;
75  res.stype = -1;
76 
77  if (internal::is_same<_StorageIndex,int>::value)
78  {
79  res.itype = CHOLMOD_INT;
80  }
81  else if (internal::is_same<_StorageIndex,SuiteSparse_long>::value)
82  {
83  res.itype = CHOLMOD_LONG;
84  }
85  else
86  {
87  eigen_assert(false && "Index type not supported yet");
88  }
89 
90  // setup res.xtype
91  internal::cholmod_configure_matrix<_Scalar>(res);
92 
93  res.stype = 0;
94 
95  return res;
96 }
97 
98 template<typename _Scalar, int _Options, typename _Index>
99 const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat)
100 {
101  cholmod_sparse res = viewAsCholmod(mat.const_cast_derived());
102  return res;
103 }
104 
107 template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo>
108 cholmod_sparse viewAsCholmod(const SparseSelfAdjointView<const SparseMatrix<_Scalar,_Options,_Index>, UpLo>& mat)
109 {
110  cholmod_sparse res = viewAsCholmod(mat.matrix().const_cast_derived());
111 
112  if(UpLo==Upper) res.stype = 1;
113  if(UpLo==Lower) res.stype = -1;
114 
115  return res;
116 }
117 
120 template<typename Derived>
121 cholmod_dense viewAsCholmod(MatrixBase<Derived>& mat)
122 {
123  EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
124  typedef typename Derived::Scalar Scalar;
125 
126  cholmod_dense res;
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());
132  res.z = 0;
133 
134  internal::cholmod_configure_matrix<Scalar>(res);
135 
136  return res;
137 }
138 
141 template<typename Scalar, int Flags, typename StorageIndex>
142 MappedSparseMatrix<Scalar,Flags,StorageIndex> viewAsEigen(cholmod_sparse& cm)
143 {
144  return MappedSparseMatrix<Scalar,Flags,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) );
147 }
148 
149 enum CholmodMode {
150  CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
151 };
152 
153 
159 template<typename _MatrixType, int _UpLo, typename Derived>
160 class CholmodBase : public SparseSolverBase<Derived>
161 {
162  protected:
164  using Base::derived;
165  using Base::m_isInitialized;
166  public:
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;
173 
174  public:
175 
176  CholmodBase()
177  : m_cholmodFactor(0), m_info(Success)
178  {
179  m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
180  cholmod_start(&m_cholmod);
181  }
182 
183  explicit CholmodBase(const MatrixType& matrix)
184  : m_cholmodFactor(0), m_info(Success)
185  {
186  m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
187  cholmod_start(&m_cholmod);
188  compute(matrix);
189  }
190 
191  ~CholmodBase()
192  {
193  if(m_cholmodFactor)
194  cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
195  cholmod_finish(&m_cholmod);
196  }
197 
198  inline StorageIndex cols() const { return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
199  inline StorageIndex rows() const { return internal::convert_index<StorageIndex, Index>(m_cholmodFactor->n); }
200 
207  {
208  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
209  return m_info;
210  }
211 
213  Derived& compute(const MatrixType& matrix)
214  {
215  analyzePattern(matrix);
216  factorize(matrix);
217  return derived();
218  }
219 
226  void analyzePattern(const MatrixType& matrix)
227  {
228  if(m_cholmodFactor)
229  {
230  cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
231  m_cholmodFactor = 0;
232  }
233  cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
234  m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
235 
236  this->m_isInitialized = true;
237  this->m_info = Success;
238  m_analysisIsOk = true;
239  m_factorizationIsOk = false;
240  }
241 
248  void factorize(const MatrixType& matrix)
249  {
250  eigen_assert(m_analysisIsOk && "You must first call analyzePattern()");
251  cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
252  cholmod_factorize_p(&A, m_shiftOffset, 0, 0, m_cholmodFactor, &m_cholmod);
253 
254  // If the factorization failed, minor is the column at which it did. On success minor == n.
255  this->m_info = (m_cholmodFactor->minor == m_cholmodFactor->n ? Success : NumericalIssue);
256  m_factorizationIsOk = true;
257  }
258 
261  cholmod_common& cholmod() { return m_cholmod; }
262 
263  #ifndef EIGEN_PARSED_BY_DOXYGEN
264 
265  template<typename Rhs,typename Dest>
266  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
267  {
268  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
269  const Index size = m_cholmodFactor->n;
270  EIGEN_UNUSED_VARIABLE(size);
271  eigen_assert(size==b.rows());
272 
273  // note: cd stands for Cholmod Dense
274  Rhs& b_ref(b.const_cast_derived());
275  cholmod_dense b_cd = viewAsCholmod(b_ref);
276  cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
277  if(!x_cd)
278  {
279  this->m_info = NumericalIssue;
280  return;
281  }
282  // TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
283  dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
284  cholmod_free_dense(&x_cd, &m_cholmod);
285  }
286 
288  template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex>
289  void _solve_impl(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
290  {
291  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
292  const Index size = m_cholmodFactor->n;
293  EIGEN_UNUSED_VARIABLE(size);
294  eigen_assert(size==b.rows());
295 
296  // note: cs stands for Cholmod Sparse
297  cholmod_sparse b_cs = viewAsCholmod(b);
298  cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
299  if(!x_cs)
300  {
301  this->m_info = NumericalIssue;
302  return;
303  }
304  // TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
305  dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
306  cholmod_free_sparse(&x_cs, &m_cholmod);
307  }
308  #endif // EIGEN_PARSED_BY_DOXYGEN
309 
310 
320  Derived& setShift(const RealScalar& offset)
321  {
322  m_shiftOffset[0] = offset;
323  return derived();
324  }
325 
326  template<typename Stream>
327  void dumpMemory(Stream& /*s*/)
328  {}
329 
330  protected:
331  mutable cholmod_common m_cholmod;
332  cholmod_factor* m_cholmodFactor;
333  RealScalar m_shiftOffset[2];
334  mutable ComputationInfo m_info;
335  int m_factorizationIsOk;
336  int m_analysisIsOk;
337 };
338 
357 template<typename _MatrixType, int _UpLo = Lower>
358 class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> >
359 {
361  using Base::m_cholmod;
362 
363  public:
364 
365  typedef _MatrixType MatrixType;
366 
367  CholmodSimplicialLLT() : Base() { init(); }
368 
369  CholmodSimplicialLLT(const MatrixType& matrix) : Base()
370  {
371  init();
372  this->compute(matrix);
373  }
374 
376  protected:
377  void init()
378  {
379  m_cholmod.final_asis = 0;
380  m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
381  m_cholmod.final_ll = 1;
382  }
383 };
384 
385 
404 template<typename _MatrixType, int _UpLo = Lower>
405 class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> >
406 {
408  using Base::m_cholmod;
409 
410  public:
411 
412  typedef _MatrixType MatrixType;
413 
414  CholmodSimplicialLDLT() : Base() { init(); }
415 
416  CholmodSimplicialLDLT(const MatrixType& matrix) : Base()
417  {
418  init();
419  this->compute(matrix);
420  }
421 
423  protected:
424  void init()
425  {
426  m_cholmod.final_asis = 1;
427  m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
428  }
429 };
430 
449 template<typename _MatrixType, int _UpLo = Lower>
450 class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> >
451 {
453  using Base::m_cholmod;
454 
455  public:
456 
457  typedef _MatrixType MatrixType;
458 
459  CholmodSupernodalLLT() : Base() { init(); }
460 
461  CholmodSupernodalLLT(const MatrixType& matrix) : Base()
462  {
463  init();
464  this->compute(matrix);
465  }
466 
468  protected:
469  void init()
470  {
471  m_cholmod.final_asis = 1;
472  m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
473  }
474 };
475 
496 template<typename _MatrixType, int _UpLo = Lower>
497 class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> >
498 {
500  using Base::m_cholmod;
501 
502  public:
503 
504  typedef _MatrixType MatrixType;
505 
506  CholmodDecomposition() : Base() { init(); }
507 
508  CholmodDecomposition(const MatrixType& matrix) : Base()
509  {
510  init();
511  this->compute(matrix);
512  }
513 
515 
516  void setMode(CholmodMode mode)
517  {
518  switch(mode)
519  {
520  case CholmodAuto:
521  m_cholmod.final_asis = 1;
522  m_cholmod.supernodal = CHOLMOD_AUTO;
523  break;
524  case CholmodSimplicialLLt:
525  m_cholmod.final_asis = 0;
526  m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
527  m_cholmod.final_ll = 1;
528  break;
529  case CholmodSupernodalLLt:
530  m_cholmod.final_asis = 1;
531  m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
532  break;
533  case CholmodLDLt:
534  m_cholmod.final_asis = 1;
535  m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
536  break;
537  default:
538  break;
539  }
540  }
541  protected:
542  void init()
543  {
544  m_cholmod.final_asis = 1;
545  m_cholmod.supernodal = CHOLMOD_AUTO;
546  }
547 };
548 
549 } // end namespace Eigen
550 
551 #endif // EIGEN_CHOLMODSUPPORT_H
Definition: Constants.h:196
void factorize(const MatrixType &matrix)
Definition: CholmodSupport.h:248
A base class for sparse solvers.
Definition: SparseSolverBase.h:53
Definition: LDLT.h:16
A supernodal Cholesky (LLT) factorization and solver based on Cholmod.
Definition: CholmodSupport.h:450
const unsigned int RowMajorBit
Definition: Constants.h:53
cholmod_common & cholmod()
Definition: CholmodSupport.h:261
Definition: Constants.h:198
Definition: Constants.h:426
void analyzePattern(const MatrixType &matrix)
Definition: CholmodSupport.h:226
A general Cholesky factorization and solver based on Cholmod.
Definition: CholmodSupport.h:497
The base class for the direct Cholesky factorization of Cholmod.
Definition: CholmodSupport.h:160
Definition: Constants.h:424
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: CholmodSupport.h:206
A simplicial direct Cholesky (LLT) factorization and solver based on Cholmod.
Definition: CholmodSupport.h:358
Definition: Eigen_Colamd.h:54
A simplicial direct Cholesky (LDLT) factorization and solver based on Cholmod.
Definition: CholmodSupport.h:405
Derived & compute(const MatrixType &matrix)
Definition: CholmodSupport.h:213
Derived & setShift(const RealScalar &offset)
Definition: CholmodSupport.h:320
ComputationInfo
Definition: Constants.h:422
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48