Eigen  3.2.93
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>
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,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>
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>
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>
143 {
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  enum {
174  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
175  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
176  };
177 
178  public:
179 
180  CholmodBase()
181  : m_cholmodFactor(0), m_info(Success), m_factorizationIsOk(false), m_analysisIsOk(false)
182  {
183  m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(0.0);
184  cholmod_start(&m_cholmod);
185  }
186 
187  explicit CholmodBase(const MatrixType& matrix)
188  : m_cholmodFactor(0), m_info(Success), m_factorizationIsOk(false), m_analysisIsOk(false)
189  {
190  m_shiftOffset[0] = m_shiftOffset[1] = RealScalar(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 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); }
204 
211  {
212  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
213  return m_info;
214  }
215 
217  Derived& compute(const MatrixType& matrix)
218  {
219  analyzePattern(matrix);
220  factorize(matrix);
221  return derived();
222  }
223 
230  void analyzePattern(const MatrixType& matrix)
231  {
232  if(m_cholmodFactor)
233  {
234  cholmod_free_factor(&m_cholmodFactor, &m_cholmod);
235  m_cholmodFactor = 0;
236  }
237  cholmod_sparse A = viewAsCholmod(matrix.template selfadjointView<UpLo>());
238  m_cholmodFactor = cholmod_analyze(&A, &m_cholmod);
239 
240  this->m_isInitialized = true;
241  this->m_info = Success;
242  m_analysisIsOk = true;
243  m_factorizationIsOk = false;
244  }
245 
252  void factorize(const MatrixType& matrix)
253  {
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);
257 
258  // If the factorization failed, minor is the column at which it did. On success minor == n.
259  this->m_info = (m_cholmodFactor->minor == m_cholmodFactor->n ? Success : NumericalIssue);
260  m_factorizationIsOk = true;
261  }
262 
265  cholmod_common& cholmod() { return m_cholmod; }
266 
267  #ifndef EIGEN_PARSED_BY_DOXYGEN
268 
269  template<typename Rhs,typename Dest>
270  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
271  {
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());
276 
277  // Cholmod needs column-major stoarge without inner-stride, which corresponds to the default behavior of Ref.
279 
280  cholmod_dense b_cd = viewAsCholmod(b_ref);
281  cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
282  if(!x_cd)
283  {
284  this->m_info = NumericalIssue;
285  return;
286  }
287  // TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
288  dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
289  cholmod_free_dense(&x_cd, &m_cholmod);
290  }
291 
293  template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex>
295  {
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());
300 
301  // note: cs stands for Cholmod Sparse
302  cholmod_sparse b_cs = viewAsCholmod(b);
303  cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
304  if(!x_cs)
305  {
306  this->m_info = NumericalIssue;
307  return;
308  }
309  // TODO optimize this copy by swapping when possible (be careful with alignment, etc.)
310  dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
311  cholmod_free_sparse(&x_cs, &m_cholmod);
312  }
313  #endif // EIGEN_PARSED_BY_DOXYGEN
314 
315 
325  Derived& setShift(const RealScalar& offset)
326  {
327  m_shiftOffset[0] = offset;
328  return derived();
329  }
330 
332  Scalar determinant() const
333  {
334  using std::exp;
335  return exp(logDeterminant());
336  }
337 
339  Scalar logDeterminant() const
340  {
341  using std::log;
342  using numext::real;
343  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
344 
345  RealScalar logDet = 0;
346  Scalar *x = static_cast<Scalar*>(m_cholmodFactor->x);
347  if (m_cholmodFactor->is_super)
348  {
349  // Supernodal factorization stored as a packed list of dense column-major blocs,
350  // as described by the following structure:
351 
352  // super[k] == index of the first column of the j-th super node
353  StorageIndex *super = static_cast<StorageIndex*>(m_cholmodFactor->super);
354  // pi[k] == offset to the description of row indices
355  StorageIndex *pi = static_cast<StorageIndex*>(m_cholmodFactor->pi);
356  // px[k] == offset to the respective dense block
357  StorageIndex *px = static_cast<StorageIndex*>(m_cholmodFactor->px);
358 
359  Index nb_super_nodes = m_cholmodFactor->nsuper;
360  for (Index k=0; k < nb_super_nodes; ++k)
361  {
362  StorageIndex ncols = super[k + 1] - super[k];
363  StorageIndex nrows = pi[k + 1] - pi[k];
364 
365  Map<const Array<Scalar,1,Dynamic>, 0, InnerStride<> > sk(x + px[k], ncols, InnerStride<>(nrows+1));
366  logDet += sk.real().log().sum();
367  }
368  }
369  else
370  {
371  // Simplicial factorization stored as standard CSC matrix.
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]] ));
376  }
377  if (m_cholmodFactor->is_ll)
378  logDet *= 2.0;
379  return logDet;
380  };
381 
382  template<typename Stream>
383  void dumpMemory(Stream& /*s*/)
384  {}
385 
386  protected:
387  mutable cholmod_common m_cholmod;
388  cholmod_factor* m_cholmodFactor;
389  RealScalar m_shiftOffset[2];
390  mutable ComputationInfo m_info;
391  int m_factorizationIsOk;
392  int m_analysisIsOk;
393 };
394 
415 template<typename _MatrixType, int _UpLo = Lower>
416 class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> >
417 {
419  using Base::m_cholmod;
420 
421  public:
422 
423  typedef _MatrixType MatrixType;
424 
425  CholmodSimplicialLLT() : Base() { init(); }
426 
427  CholmodSimplicialLLT(const MatrixType& matrix) : Base()
428  {
429  init();
430  this->compute(matrix);
431  }
432 
434  protected:
435  void init()
436  {
437  m_cholmod.final_asis = 0;
438  m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
439  m_cholmod.final_ll = 1;
440  }
441 };
442 
443 
464 template<typename _MatrixType, int _UpLo = Lower>
465 class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> >
466 {
468  using Base::m_cholmod;
469 
470  public:
471 
472  typedef _MatrixType MatrixType;
473 
474  CholmodSimplicialLDLT() : Base() { init(); }
475 
476  CholmodSimplicialLDLT(const MatrixType& matrix) : Base()
477  {
478  init();
479  this->compute(matrix);
480  }
481 
483  protected:
484  void init()
485  {
486  m_cholmod.final_asis = 1;
487  m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
488  }
489 };
490 
511 template<typename _MatrixType, int _UpLo = Lower>
512 class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> >
513 {
515  using Base::m_cholmod;
516 
517  public:
518 
519  typedef _MatrixType MatrixType;
520 
521  CholmodSupernodalLLT() : Base() { init(); }
522 
523  CholmodSupernodalLLT(const MatrixType& matrix) : Base()
524  {
525  init();
526  this->compute(matrix);
527  }
528 
530  protected:
531  void init()
532  {
533  m_cholmod.final_asis = 1;
534  m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
535  }
536 };
537 
560 template<typename _MatrixType, int _UpLo = Lower>
561 class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> >
562 {
564  using Base::m_cholmod;
565 
566  public:
567 
568  typedef _MatrixType MatrixType;
569 
570  CholmodDecomposition() : Base() { init(); }
571 
572  CholmodDecomposition(const MatrixType& matrix) : Base()
573  {
574  init();
575  this->compute(matrix);
576  }
577 
579 
580  void setMode(CholmodMode mode)
581  {
582  switch(mode)
583  {
584  case CholmodAuto:
585  m_cholmod.final_asis = 1;
586  m_cholmod.supernodal = CHOLMOD_AUTO;
587  break;
588  case CholmodSimplicialLLt:
589  m_cholmod.final_asis = 0;
590  m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
591  m_cholmod.final_ll = 1;
592  break;
593  case CholmodSupernodalLLt:
594  m_cholmod.final_asis = 1;
595  m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
596  break;
597  case CholmodLDLt:
598  m_cholmod.final_asis = 1;
599  m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
600  break;
601  default:
602  break;
603  }
604  }
605  protected:
606  void init()
607  {
608  m_cholmod.final_asis = 1;
609  m_cholmod.supernodal = CHOLMOD_AUTO;
610  }
611 };
612 
613 } // end namespace Eigen
614 
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)