Eigen  3.2.0
 All Classes Namespaces Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SimplicialCholesky.h
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2008-2012 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_SIMPLICIAL_CHOLESKY_H
11 #define EIGEN_SIMPLICIAL_CHOLESKY_H
12 
13 namespace Eigen {
14 
15 enum SimplicialCholeskyMode {
16  SimplicialCholeskyLLT,
17  SimplicialCholeskyLDLT
18 };
19 
35 template<typename Derived>
36 class SimplicialCholeskyBase : internal::noncopyable
37 {
38  public:
39  typedef typename internal::traits<Derived>::MatrixType MatrixType;
40  enum { UpLo = internal::traits<Derived>::UpLo };
41  typedef typename MatrixType::Scalar Scalar;
42  typedef typename MatrixType::RealScalar RealScalar;
43  typedef typename MatrixType::Index Index;
46 
47  public:
48 
51  : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
52  {}
53 
54  SimplicialCholeskyBase(const MatrixType& matrix)
55  : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
56  {
57  derived().compute(matrix);
58  }
59 
61  {
62  }
63 
64  Derived& derived() { return *static_cast<Derived*>(this); }
65  const Derived& derived() const { return *static_cast<const Derived*>(this); }
66 
67  inline Index cols() const { return m_matrix.cols(); }
68  inline Index rows() const { return m_matrix.rows(); }
69 
76  {
77  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
78  return m_info;
79  }
80 
85  template<typename Rhs>
86  inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
87  solve(const MatrixBase<Rhs>& b) const
88  {
89  eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
90  eigen_assert(rows()==b.rows()
91  && "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
92  return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
93  }
94 
99  template<typename Rhs>
100  inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
101  solve(const SparseMatrixBase<Rhs>& b) const
102  {
103  eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
104  eigen_assert(rows()==b.rows()
105  && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
106  return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
107  }
108 
112  { return m_P; }
113 
117  { return m_Pinv; }
118 
128  Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1)
129  {
130  m_shiftOffset = offset;
131  m_shiftScale = scale;
132  return derived();
133  }
134 
135 #ifndef EIGEN_PARSED_BY_DOXYGEN
136 
137  template<typename Stream>
138  void dumpMemory(Stream& s)
139  {
140  int total = 0;
141  s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
142  s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
143  s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
144  s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
145  s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
146  s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
147  s << " TOTAL: " << (total>> 20) << "Mb" << "\n";
148  }
149 
151  template<typename Rhs,typename Dest>
152  void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
153  {
154  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
155  eigen_assert(m_matrix.rows()==b.rows());
156 
157  if(m_info!=Success)
158  return;
159 
160  if(m_P.size()>0)
161  dest = m_P * b;
162  else
163  dest = b;
164 
165  if(m_matrix.nonZeros()>0) // otherwise L==I
166  derived().matrixL().solveInPlace(dest);
167 
168  if(m_diag.size()>0)
169  dest = m_diag.asDiagonal().inverse() * dest;
170 
171  if (m_matrix.nonZeros()>0) // otherwise U==I
172  derived().matrixU().solveInPlace(dest);
173 
174  if(m_P.size()>0)
175  dest = m_Pinv * dest;
176  }
177 
178 #endif // EIGEN_PARSED_BY_DOXYGEN
179 
180  protected:
181 
183  template<bool DoLDLT>
184  void compute(const MatrixType& matrix)
185  {
186  eigen_assert(matrix.rows()==matrix.cols());
187  Index size = matrix.cols();
188  CholMatrixType ap(size,size);
189  ordering(matrix, ap);
190  analyzePattern_preordered(ap, DoLDLT);
191  factorize_preordered<DoLDLT>(ap);
192  }
193 
194  template<bool DoLDLT>
195  void factorize(const MatrixType& a)
196  {
197  eigen_assert(a.rows()==a.cols());
198  int size = a.cols();
199  CholMatrixType ap(size,size);
200  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
201  factorize_preordered<DoLDLT>(ap);
202  }
203 
204  template<bool DoLDLT>
205  void factorize_preordered(const CholMatrixType& a);
206 
207  void analyzePattern(const MatrixType& a, bool doLDLT)
208  {
209  eigen_assert(a.rows()==a.cols());
210  int size = a.cols();
211  CholMatrixType ap(size,size);
212  ordering(a, ap);
213  analyzePattern_preordered(ap,doLDLT);
214  }
215  void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
216 
217  void ordering(const MatrixType& a, CholMatrixType& ap);
218 
220  struct keep_diag {
221  inline bool operator() (const Index& row, const Index& col, const Scalar&) const
222  {
223  return row!=col;
224  }
225  };
226 
227  mutable ComputationInfo m_info;
228  bool m_isInitialized;
229  bool m_factorizationIsOk;
230  bool m_analysisIsOk;
231 
232  CholMatrixType m_matrix;
233  VectorType m_diag; // the diagonal coefficients (LDLT mode)
234  VectorXi m_parent; // elimination tree
235  VectorXi m_nonZerosPerCol;
236  PermutationMatrix<Dynamic,Dynamic,Index> m_P; // the permutation
237  PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // the inverse permutation
238 
239  RealScalar m_shiftOffset;
240  RealScalar m_shiftScale;
241 };
242 
243 template<typename _MatrixType, int _UpLo = Lower> class SimplicialLLT;
244 template<typename _MatrixType, int _UpLo = Lower> class SimplicialLDLT;
245 template<typename _MatrixType, int _UpLo = Lower> class SimplicialCholesky;
246 
247 namespace internal {
248 
249 template<typename _MatrixType, int _UpLo> struct traits<SimplicialLLT<_MatrixType,_UpLo> >
250 {
251  typedef _MatrixType MatrixType;
252  enum { UpLo = _UpLo };
253  typedef typename MatrixType::Scalar Scalar;
254  typedef typename MatrixType::Index Index;
255  typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
256  typedef SparseTriangularView<CholMatrixType, Eigen::Lower> MatrixL;
257  typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
258  static inline MatrixL getL(const MatrixType& m) { return m; }
259  static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
260 };
261 
262 template<typename _MatrixType,int _UpLo> struct traits<SimplicialLDLT<_MatrixType,_UpLo> >
263 {
264  typedef _MatrixType MatrixType;
265  enum { UpLo = _UpLo };
266  typedef typename MatrixType::Scalar Scalar;
267  typedef typename MatrixType::Index Index;
268  typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
269  typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower> MatrixL;
270  typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
271  static inline MatrixL getL(const MatrixType& m) { return m; }
272  static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
273 };
274 
275 template<typename _MatrixType, int _UpLo> struct traits<SimplicialCholesky<_MatrixType,_UpLo> >
276 {
277  typedef _MatrixType MatrixType;
278  enum { UpLo = _UpLo };
279 };
280 
281 }
282 
300 template<typename _MatrixType, int _UpLo>
301  class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo> >
302 {
303 public:
304  typedef _MatrixType MatrixType;
305  enum { UpLo = _UpLo };
306  typedef SimplicialCholeskyBase<SimplicialLLT> Base;
307  typedef typename MatrixType::Scalar Scalar;
308  typedef typename MatrixType::RealScalar RealScalar;
309  typedef typename MatrixType::Index Index;
310  typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
311  typedef Matrix<Scalar,Dynamic,1> VectorType;
312  typedef internal::traits<SimplicialLLT> Traits;
313  typedef typename Traits::MatrixL MatrixL;
314  typedef typename Traits::MatrixU MatrixU;
315 public:
319  SimplicialLLT(const MatrixType& matrix)
320  : Base(matrix) {}
321 
323  inline const MatrixL matrixL() const {
324  eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
325  return Traits::getL(Base::m_matrix);
326  }
327 
329  inline const MatrixU matrixU() const {
330  eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
331  return Traits::getU(Base::m_matrix);
332  }
333 
335  SimplicialLLT& compute(const MatrixType& matrix)
336  {
337  Base::template compute<false>(matrix);
338  return *this;
339  }
340 
347  void analyzePattern(const MatrixType& a)
348  {
349  Base::analyzePattern(a, false);
350  }
351 
358  void factorize(const MatrixType& a)
359  {
360  Base::template factorize<false>(a);
361  }
362 
364  Scalar determinant() const
365  {
366  Scalar detL = Base::m_matrix.diagonal().prod();
367  return numext::abs2(detL);
368  }
369 };
370 
388 template<typename _MatrixType, int _UpLo>
389  class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo> >
390 {
391 public:
392  typedef _MatrixType MatrixType;
393  enum { UpLo = _UpLo };
394  typedef SimplicialCholeskyBase<SimplicialLDLT> Base;
395  typedef typename MatrixType::Scalar Scalar;
396  typedef typename MatrixType::RealScalar RealScalar;
397  typedef typename MatrixType::Index Index;
398  typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
399  typedef Matrix<Scalar,Dynamic,1> VectorType;
400  typedef internal::traits<SimplicialLDLT> Traits;
401  typedef typename Traits::MatrixL MatrixL;
402  typedef typename Traits::MatrixU MatrixU;
403 public:
406 
408  SimplicialLDLT(const MatrixType& matrix)
409  : Base(matrix) {}
410 
412  inline const VectorType vectorD() const {
413  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
414  return Base::m_diag;
415  }
417  inline const MatrixL matrixL() const {
418  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
419  return Traits::getL(Base::m_matrix);
420  }
421 
423  inline const MatrixU matrixU() const {
424  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
425  return Traits::getU(Base::m_matrix);
426  }
427 
429  SimplicialLDLT& compute(const MatrixType& matrix)
430  {
431  Base::template compute<true>(matrix);
432  return *this;
433  }
434 
441  void analyzePattern(const MatrixType& a)
442  {
443  Base::analyzePattern(a, true);
444  }
445 
452  void factorize(const MatrixType& a)
453  {
454  Base::template factorize<true>(a);
455  }
456 
458  Scalar determinant() const
459  {
460  return Base::m_diag.prod();
461  }
462 };
463 
470 template<typename _MatrixType, int _UpLo>
471  class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo> >
472 {
473 public:
474  typedef _MatrixType MatrixType;
475  enum { UpLo = _UpLo };
476  typedef SimplicialCholeskyBase<SimplicialCholesky> Base;
477  typedef typename MatrixType::Scalar Scalar;
478  typedef typename MatrixType::RealScalar RealScalar;
479  typedef typename MatrixType::Index Index;
480  typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
481  typedef Matrix<Scalar,Dynamic,1> VectorType;
482  typedef internal::traits<SimplicialCholesky> Traits;
483  typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
484  typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
485  public:
486  SimplicialCholesky() : Base(), m_LDLT(true) {}
487 
488  SimplicialCholesky(const MatrixType& matrix)
489  : Base(), m_LDLT(true)
490  {
491  compute(matrix);
492  }
493 
494  SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
495  {
496  switch(mode)
497  {
498  case SimplicialCholeskyLLT:
499  m_LDLT = false;
500  break;
501  case SimplicialCholeskyLDLT:
502  m_LDLT = true;
503  break;
504  default:
505  break;
506  }
507 
508  return *this;
509  }
510 
511  inline const VectorType vectorD() const {
512  eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
513  return Base::m_diag;
514  }
515  inline const CholMatrixType rawMatrix() const {
516  eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
517  return Base::m_matrix;
518  }
519 
521  SimplicialCholesky& compute(const MatrixType& matrix)
522  {
523  if(m_LDLT)
524  Base::template compute<true>(matrix);
525  else
526  Base::template compute<false>(matrix);
527  return *this;
528  }
529 
536  void analyzePattern(const MatrixType& a)
537  {
538  Base::analyzePattern(a, m_LDLT);
539  }
540 
547  void factorize(const MatrixType& a)
548  {
549  if(m_LDLT)
550  Base::template factorize<true>(a);
551  else
552  Base::template factorize<false>(a);
553  }
554 
556  template<typename Rhs,typename Dest>
557  void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
558  {
559  eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
560  eigen_assert(Base::m_matrix.rows()==b.rows());
561 
562  if(Base::m_info!=Success)
563  return;
564 
565  if(Base::m_P.size()>0)
566  dest = Base::m_P * b;
567  else
568  dest = b;
569 
570  if(Base::m_matrix.nonZeros()>0) // otherwise L==I
571  {
572  if(m_LDLT)
573  LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
574  else
575  LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
576  }
577 
578  if(Base::m_diag.size()>0)
579  dest = Base::m_diag.asDiagonal().inverse() * dest;
580 
581  if (Base::m_matrix.nonZeros()>0) // otherwise I==I
582  {
583  if(m_LDLT)
584  LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
585  else
586  LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
587  }
588 
589  if(Base::m_P.size()>0)
590  dest = Base::m_Pinv * dest;
591  }
592 
593  Scalar determinant() const
594  {
595  if(m_LDLT)
596  {
597  return Base::m_diag.prod();
598  }
599  else
600  {
601  Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
602  return numext::abs2(detL);
603  }
604  }
605 
606  protected:
607  bool m_LDLT;
608 };
609 
610 template<typename Derived>
611 void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixType& ap)
612 {
613  eigen_assert(a.rows()==a.cols());
614  const Index size = a.rows();
615  // TODO allows to configure the permutation
616  // Note that amd compute the inverse permutation
617  {
618  CholMatrixType C;
619  C = a.template selfadjointView<UpLo>();
620  // remove diagonal entries:
621  // seems not to be needed
622  // C.prune(keep_diag());
623  internal::minimum_degree_ordering(C, m_Pinv);
624  }
625 
626  if(m_Pinv.size()>0)
627  m_P = m_Pinv.inverse();
628  else
629  m_P.resize(0);
630 
631  ap.resize(size,size);
632  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
633 }
634 
635 namespace internal {
636 
637 template<typename Derived, typename Rhs>
638 struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
639  : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
640 {
641  typedef SimplicialCholeskyBase<Derived> Dec;
642  EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
643 
644  template<typename Dest> void evalTo(Dest& dst) const
645  {
646  dec().derived()._solve(rhs(),dst);
647  }
648 };
649 
650 template<typename Derived, typename Rhs>
651 struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
652  : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
653 {
654  typedef SimplicialCholeskyBase<Derived> Dec;
655  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
656 
657  template<typename Dest> void evalTo(Dest& dst) const
658  {
659  this->defaultEvalTo(dst);
660  }
661 };
662 
663 } // end namespace internal
664 
665 } // end namespace Eigen
666 
667 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H