CholmodSupport.h
Go to the documentation of this file.
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 // Eigen is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 3 of the License, or (at your option) any later version.
10 //
11 // Alternatively, you can redistribute it and/or
12 // modify it under the terms of the GNU General Public License as
13 // published by the Free Software Foundation; either version 2 of
14 // the License, or (at your option) any later version.
15 //
16 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
17 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
18 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
19 // GNU General Public License for more details.
20 //
21 // You should have received a copy of the GNU Lesser General Public
22 // License and a copy of the GNU General Public License along with
23 // Eigen. If not, see <http://www.gnu.org/licenses/>.
24 
25 #ifndef EIGEN_CHOLMODSUPPORT_H
26 #define EIGEN_CHOLMODSUPPORT_H
27 
28 namespace Eigen {
29 
30 namespace internal {
31 
32 template<typename Scalar, typename CholmodType>
33 void cholmod_configure_matrix(CholmodType& mat)
34 {
35  if (internal::is_same<Scalar,float>::value)
36  {
37  mat.xtype = CHOLMOD_REAL;
38  mat.dtype = CHOLMOD_SINGLE;
39  }
40  else if (internal::is_same<Scalar,double>::value)
41  {
42  mat.xtype = CHOLMOD_REAL;
43  mat.dtype = CHOLMOD_DOUBLE;
44  }
45  else if (internal::is_same<Scalar,std::complex<float> >::value)
46  {
47  mat.xtype = CHOLMOD_COMPLEX;
48  mat.dtype = CHOLMOD_SINGLE;
49  }
50  else if (internal::is_same<Scalar,std::complex<double> >::value)
51  {
52  mat.xtype = CHOLMOD_COMPLEX;
53  mat.dtype = CHOLMOD_DOUBLE;
54  }
55  else
56  {
57  eigen_assert(false && "Scalar type not supported by CHOLMOD");
58  }
59 }
60 
61 } // namespace internal
62 
66 template<typename _Scalar, int _Options, typename _Index>
68 {
69  typedef SparseMatrix<_Scalar,_Options,_Index> MatrixType;
70  cholmod_sparse res;
71  res.nzmax = mat.nonZeros();
72  res.nrow = mat.rows();;
73  res.ncol = mat.cols();
74  res.p = mat.outerIndexPtr();
75  res.i = mat.innerIndexPtr();
76  res.x = mat.valuePtr();
77  res.sorted = 1;
78  if(mat.isCompressed())
79  {
80  res.packed = 1;
81  }
82  else
83  {
84  res.packed = 0;
85  res.nz = mat.innerNonZeroPtr();
86  }
87 
88  res.dtype = 0;
89  res.stype = -1;
90 
91  if (internal::is_same<_Index,int>::value)
92  {
93  res.itype = CHOLMOD_INT;
94  }
95  else
96  {
97  eigen_assert(false && "Index type different than int is not supported yet");
98  }
99 
100  // setup res.xtype
101  internal::cholmod_configure_matrix<_Scalar>(res);
102 
103  res.stype = 0;
104 
105  return res;
106 }
107 
108 template<typename _Scalar, int _Options, typename _Index>
109 const cholmod_sparse viewAsCholmod(const SparseMatrix<_Scalar,_Options,_Index>& mat)
110 {
111  cholmod_sparse res = viewAsCholmod(mat.const_cast_derived());
112  return res;
113 }
114 
117 template<typename _Scalar, int _Options, typename _Index, unsigned int UpLo>
119 {
120  cholmod_sparse res = viewAsCholmod(mat.matrix().const_cast_derived());
121 
122  if(UpLo==Upper) res.stype = 1;
123  if(UpLo==Lower) res.stype = -1;
124 
125  return res;
126 }
127 
130 template<typename Derived>
132 {
133  EIGEN_STATIC_ASSERT((internal::traits<Derived>::Flags&RowMajorBit)==0,THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
134  typedef typename Derived::Scalar Scalar;
135 
136  cholmod_dense res;
137  res.nrow = mat.rows();
138  res.ncol = mat.cols();
139  res.nzmax = res.nrow * res.ncol;
140  res.d = Derived::IsVectorAtCompileTime ? mat.derived().size() : mat.derived().outerStride();
141  res.x = mat.derived().data();
142  res.z = 0;
143 
144  internal::cholmod_configure_matrix<Scalar>(res);
145 
146  return res;
147 }
148 
151 template<typename Scalar, int Flags, typename Index>
153 {
155  (cm.nrow, cm.ncol, reinterpret_cast<Index*>(cm.p)[cm.ncol],
156  reinterpret_cast<Index*>(cm.p), reinterpret_cast<Index*>(cm.i),reinterpret_cast<Scalar*>(cm.x) );
157 }
158 
161 };
162 
163 
169 template<typename _MatrixType, int _UpLo, typename Derived>
170 class CholmodBase : internal::noncopyable
171 {
172  public:
173  typedef _MatrixType MatrixType;
174  enum { UpLo = _UpLo };
175  typedef typename MatrixType::Scalar Scalar;
176  typedef typename MatrixType::RealScalar RealScalar;
178  typedef typename MatrixType::Index Index;
179 
180  public:
181 
184  {
185  cholmod_start(&m_cholmod);
186  }
187 
188  CholmodBase(const MatrixType& matrix)
190  {
191  cholmod_start(&m_cholmod);
192  compute(matrix);
193  }
194 
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 
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(&A, m_cholmodFactor, &m_cholmod);
288 
289  this->m_info = Success;
290  m_factorizationIsOk = true;
291  }
292 
296 
297  #ifndef EIGEN_PARSED_BY_DOXYGEN
298 
299  template<typename Rhs,typename Dest>
300  void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
301  {
302  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
303  const Index size = m_cholmodFactor->n;
304  eigen_assert(size==b.rows());
305 
306  // note: cd stands for Cholmod Dense
307  cholmod_dense b_cd = viewAsCholmod(b.const_cast_derived());
308  cholmod_dense* x_cd = cholmod_solve(CHOLMOD_A, m_cholmodFactor, &b_cd, &m_cholmod);
309  if(!x_cd)
310  {
311  this->m_info = NumericalIssue;
312  }
313  // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.)
314  dest = Matrix<Scalar,Dest::RowsAtCompileTime,Dest::ColsAtCompileTime>::Map(reinterpret_cast<Scalar*>(x_cd->x),b.rows(),b.cols());
315  cholmod_free_dense(&x_cd, &m_cholmod);
316  }
317 
319  template<typename RhsScalar, int RhsOptions, typename RhsIndex, typename DestScalar, int DestOptions, typename DestIndex>
320  void _solve(const SparseMatrix<RhsScalar,RhsOptions,RhsIndex> &b, SparseMatrix<DestScalar,DestOptions,DestIndex> &dest) const
321  {
322  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
323  const Index size = m_cholmodFactor->n;
324  eigen_assert(size==b.rows());
325 
326  // note: cs stands for Cholmod Sparse
327  cholmod_sparse b_cs = viewAsCholmod(b);
328  cholmod_sparse* x_cs = cholmod_spsolve(CHOLMOD_A, m_cholmodFactor, &b_cs, &m_cholmod);
329  if(!x_cs)
330  {
331  this->m_info = NumericalIssue;
332  }
333  // TODO optimize this copy by swapping when possible (be carreful with alignment, etc.)
334  dest = viewAsEigen<DestScalar,DestOptions,DestIndex>(*x_cs);
335  cholmod_free_sparse(&x_cs, &m_cholmod);
336  }
337  #endif // EIGEN_PARSED_BY_DOXYGEN
338 
339  template<typename Stream>
340  void dumpMemory(Stream& s)
341  {}
342 
343  protected:
345  cholmod_factor* m_cholmodFactor;
350 };
351 
370 template<typename _MatrixType, int _UpLo = Lower>
371 class CholmodSimplicialLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLLT<_MatrixType, _UpLo> >
372 {
374  using Base::m_cholmod;
375 
376  public:
377 
378  typedef _MatrixType MatrixType;
379 
381 
383  {
384  init();
385  compute(matrix);
386  }
387 
389  protected:
390  void init()
391  {
392  m_cholmod.final_asis = 0;
393  m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
394  m_cholmod.final_ll = 1;
395  }
396 };
397 
398 
417 template<typename _MatrixType, int _UpLo = Lower>
418 class CholmodSimplicialLDLT : public CholmodBase<_MatrixType, _UpLo, CholmodSimplicialLDLT<_MatrixType, _UpLo> >
419 {
421  using Base::m_cholmod;
422 
423  public:
424 
425  typedef _MatrixType MatrixType;
426 
428 
430  {
431  init();
432  compute(matrix);
433  }
434 
436  protected:
437  void init()
438  {
439  m_cholmod.final_asis = 1;
440  m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
441  }
442 };
443 
462 template<typename _MatrixType, int _UpLo = Lower>
463 class CholmodSupernodalLLT : public CholmodBase<_MatrixType, _UpLo, CholmodSupernodalLLT<_MatrixType, _UpLo> >
464 {
466  using Base::m_cholmod;
467 
468  public:
469 
470  typedef _MatrixType MatrixType;
471 
473 
475  {
476  init();
477  compute(matrix);
478  }
479 
481  protected:
482  void init()
483  {
484  m_cholmod.final_asis = 1;
485  m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
486  }
487 };
488 
509 template<typename _MatrixType, int _UpLo = Lower>
510 class CholmodDecomposition : public CholmodBase<_MatrixType, _UpLo, CholmodDecomposition<_MatrixType, _UpLo> >
511 {
513  using Base::m_cholmod;
514 
515  public:
516 
517  typedef _MatrixType MatrixType;
518 
520 
522  {
523  init();
524  compute(matrix);
525  }
526 
528 
529  void setMode(CholmodMode mode)
530  {
531  switch(mode)
532  {
533  case CholmodAuto:
534  m_cholmod.final_asis = 1;
535  m_cholmod.supernodal = CHOLMOD_AUTO;
536  break;
538  m_cholmod.final_asis = 0;
539  m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
540  m_cholmod.final_ll = 1;
541  break;
543  m_cholmod.final_asis = 1;
544  m_cholmod.supernodal = CHOLMOD_SUPERNODAL;
545  break;
546  case CholmodLDLt:
547  m_cholmod.final_asis = 1;
548  m_cholmod.supernodal = CHOLMOD_SIMPLICIAL;
549  break;
550  default:
551  break;
552  }
553  }
554  protected:
555  void init()
556  {
557  m_cholmod.final_asis = 1;
558  m_cholmod.supernodal = CHOLMOD_AUTO;
559  }
560 };
561 
562 namespace internal {
563 
564 template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
565 struct solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
566  : solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
567 {
568  typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
569  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
570 
571  template<typename Dest> void evalTo(Dest& dst) const
572  {
573  dec()._solve(rhs(),dst);
574  }
575 };
576 
577 template<typename _MatrixType, int _UpLo, typename Derived, typename Rhs>
578 struct sparse_solve_retval<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
579  : sparse_solve_retval_base<CholmodBase<_MatrixType,_UpLo,Derived>, Rhs>
580 {
581  typedef CholmodBase<_MatrixType,_UpLo,Derived> Dec;
583 
584  template<typename Dest> void evalTo(Dest& dst) const
585  {
586  dec()._solve(rhs(),dst);
587  }
588 };
589 
590 } // end namespace internal
591 
592 } // end namespace Eigen
593 
594 #endif // EIGEN_CHOLMODSUPPORT_H