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 _Index>
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 
93 template<typename _Scalar, int _Options, typename _Index>
94 const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat)
95 {
96  cholmod_sparse res = viewAsCholmod(mat.const_cast_derived());
97  return res;
98 }
99 
102 template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo>
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 
115 template<typename Derived>
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 
136 template<typename Scalar, int Flags, typename Index>
138 {
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 
144 enum CholmodMode {
145  CholmodAuto, CholmodSimplicialLLt, CholmodSupernodalLLt, CholmodLDLt
146 };
147 
148 
154 template<typename _MatrixType, int _UpLo, typename Derived>
155 class 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 
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
283 
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 
355 template<typename _MatrixType, int _UpLo = Lower>
356 class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> >
357 {
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 
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 
402 template<typename _MatrixType, int _UpLo = Lower>
403 class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> >
404 {
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 
421  protected:
422  void init()
423  {
424  m_cholmod.final_asis = 1;
425  m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
426  }
427 };
428 
447 template<typename _MatrixType, int _UpLo = Lower>
448 class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> >
449 {
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 
466  protected:
467  void init()
468  {
469  m_cholmod.final_asis = 1;
470  m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
471  }
472 };
473 
494 template<typename _MatrixType, int _UpLo = Lower>
495 class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> >
496 {
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 
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 
547 namespace internal {
548 
549 template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
550 struct solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
551  : solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
552 {
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 
562 template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
563 struct 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