Eigen  3.2.91
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 
20 namespace internal {
21  template<typename CholMatrixType, typename InputMatrixType>
22  struct simplicial_cholesky_grab_input {
23  typedef CholMatrixType const * ConstCholMatrixPtr;
24  static void run(const InputMatrixType& input, ConstCholMatrixPtr &pmat, CholMatrixType &tmp)
25  {
26  tmp = input;
27  pmat = &tmp;
28  }
29  };
30 
31  template<typename MatrixType>
32  struct simplicial_cholesky_grab_input<MatrixType,MatrixType> {
33  typedef MatrixType const * ConstMatrixPtr;
34  static void run(const MatrixType& input, ConstMatrixPtr &pmat, MatrixType &/*tmp*/)
35  {
36  pmat = &input;
37  }
38  };
39 } // end namespace internal
40 
56 template<typename Derived>
58 {
60  using Base::m_isInitialized;
61 
62  public:
63  typedef typename internal::traits<Derived>::MatrixType MatrixType;
64  typedef typename internal::traits<Derived>::OrderingType OrderingType;
65  enum { UpLo = internal::traits<Derived>::UpLo };
66  typedef typename MatrixType::Scalar Scalar;
67  typedef typename MatrixType::RealScalar RealScalar;
68  typedef typename MatrixType::StorageIndex StorageIndex;
70  typedef CholMatrixType const * ConstCholMatrixPtr;
73 
74  public:
75 
76  using Base::derived;
77 
80  : m_info(Success), m_shiftOffset(0), m_shiftScale(1)
81  {}
82 
83  explicit SimplicialCholeskyBase(const MatrixType& matrix)
84  : m_info(Success), m_shiftOffset(0), m_shiftScale(1)
85  {
86  derived().compute(matrix);
87  }
88 
90  {
91  }
92 
93  Derived& derived() { return *static_cast<Derived*>(this); }
94  const Derived& derived() const { return *static_cast<const Derived*>(this); }
95 
96  inline Index cols() const { return m_matrix.cols(); }
97  inline Index rows() const { return m_matrix.rows(); }
98 
105  {
106  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
107  return m_info;
108  }
109 
113  { return m_P; }
114 
118  { return m_Pinv; }
119 
129  Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1)
130  {
131  m_shiftOffset = offset;
132  m_shiftScale = scale;
133  return derived();
134  }
135 
136 #ifndef EIGEN_PARSED_BY_DOXYGEN
137 
138  template<typename Stream>
139  void dumpMemory(Stream& s)
140  {
141  int total = 0;
142  s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
143  s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
144  s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
145  s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
146  s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
147  s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
148  s << " TOTAL: " << (total>> 20) << "Mb" << "\n";
149  }
150 
152  template<typename Rhs,typename Dest>
153  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
154  {
155  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
156  eigen_assert(m_matrix.rows()==b.rows());
157 
158  if(m_info!=Success)
159  return;
160 
161  if(m_P.size()>0)
162  dest = m_P * b;
163  else
164  dest = b;
165 
166  if(m_matrix.nonZeros()>0) // otherwise L==I
167  derived().matrixL().solveInPlace(dest);
168 
169  if(m_diag.size()>0)
170  dest = m_diag.asDiagonal().inverse() * dest;
171 
172  if (m_matrix.nonZeros()>0) // otherwise U==I
173  derived().matrixU().solveInPlace(dest);
174 
175  if(m_P.size()>0)
176  dest = m_Pinv * dest;
177  }
178 
179  template<typename Rhs,typename Dest>
180  void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const
181  {
182  internal::solve_sparse_through_dense_panels(derived(), b, dest);
183  }
184 
185 #endif // EIGEN_PARSED_BY_DOXYGEN
186 
187  protected:
188 
190  template<bool DoLDLT>
191  void compute(const MatrixType& matrix)
192  {
193  eigen_assert(matrix.rows()==matrix.cols());
194  Index size = matrix.cols();
195  CholMatrixType tmp(size,size);
196  ConstCholMatrixPtr pmat;
197  ordering(matrix, pmat, tmp);
198  analyzePattern_preordered(*pmat, DoLDLT);
199  factorize_preordered<DoLDLT>(*pmat);
200  }
201 
202  template<bool DoLDLT>
203  void factorize(const MatrixType& a)
204  {
205  eigen_assert(a.rows()==a.cols());
206  Index size = a.cols();
207  CholMatrixType tmp(size,size);
208  ConstCholMatrixPtr pmat;
209 
210  if(m_P.size()==0 && (UpLo&Upper)==Upper)
211  {
212  // If there is no ordering, try to directly use the input matrix without any copy
213  internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp);
214  }
215  else
216  {
217  tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
218  pmat = &tmp;
219  }
220 
221  factorize_preordered<DoLDLT>(*pmat);
222  }
223 
224  template<bool DoLDLT>
225  void factorize_preordered(const CholMatrixType& a);
226 
227  void analyzePattern(const MatrixType& a, bool doLDLT)
228  {
229  eigen_assert(a.rows()==a.cols());
230  Index size = a.cols();
231  CholMatrixType tmp(size,size);
232  ConstCholMatrixPtr pmat;
233  ordering(a, pmat, tmp);
234  analyzePattern_preordered(*pmat,doLDLT);
235  }
236  void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
237 
238  void ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap);
239 
241  struct keep_diag {
242  inline bool operator() (const Index& row, const Index& col, const Scalar&) const
243  {
244  return row!=col;
245  }
246  };
247 
248  mutable ComputationInfo m_info;
249  bool m_factorizationIsOk;
250  bool m_analysisIsOk;
251 
252  CholMatrixType m_matrix;
253  VectorType m_diag; // the diagonal coefficients (LDLT mode)
254  VectorI m_parent; // elimination tree
255  VectorI m_nonZerosPerCol;
257  PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_Pinv; // the inverse permutation
258 
259  RealScalar m_shiftOffset;
260  RealScalar m_shiftScale;
261 };
262 
263 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLLT;
264 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLDLT;
265 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialCholesky;
266 
267 namespace internal {
268 
269 template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
270 {
271  typedef _MatrixType MatrixType;
272  typedef _Ordering OrderingType;
273  enum { UpLo = _UpLo };
274  typedef typename MatrixType::Scalar Scalar;
275  typedef typename MatrixType::StorageIndex StorageIndex;
276  typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
279  static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); }
280  static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
281 };
282 
283 template<typename _MatrixType,int _UpLo, typename _Ordering> struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
284 {
285  typedef _MatrixType MatrixType;
286  typedef _Ordering OrderingType;
287  enum { UpLo = _UpLo };
288  typedef typename MatrixType::Scalar Scalar;
289  typedef typename MatrixType::StorageIndex StorageIndex;
290  typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
291  typedef TriangularView<const CholMatrixType, Eigen::UnitLower> MatrixL;
292  typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
293  static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); }
294  static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
295 };
296 
297 template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
298 {
299  typedef _MatrixType MatrixType;
300  typedef _Ordering OrderingType;
301  enum { UpLo = _UpLo };
302 };
303 
304 }
305 
324 template<typename _MatrixType, int _UpLo, typename _Ordering>
325  class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
326 {
327 public:
328  typedef _MatrixType MatrixType;
329  enum { UpLo = _UpLo };
330  typedef SimplicialCholeskyBase<SimplicialLLT> Base;
331  typedef typename MatrixType::Scalar Scalar;
332  typedef typename MatrixType::RealScalar RealScalar;
333  typedef typename MatrixType::StorageIndex StorageIndex;
334  typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
335  typedef Matrix<Scalar,Dynamic,1> VectorType;
336  typedef internal::traits<SimplicialLLT> Traits;
337  typedef typename Traits::MatrixL MatrixL;
338  typedef typename Traits::MatrixU MatrixU;
339 public:
341  SimplicialLLT() : Base() {}
343  explicit SimplicialLLT(const MatrixType& matrix)
344  : Base(matrix) {}
345 
347  inline const MatrixL matrixL() const {
348  eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
349  return Traits::getL(Base::m_matrix);
350  }
351 
353  inline const MatrixU matrixU() const {
354  eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
355  return Traits::getU(Base::m_matrix);
356  }
357 
359  SimplicialLLT& compute(const MatrixType& matrix)
360  {
361  Base::template compute<false>(matrix);
362  return *this;
363  }
364 
371  void analyzePattern(const MatrixType& a)
372  {
373  Base::analyzePattern(a, false);
374  }
375 
382  void factorize(const MatrixType& a)
383  {
384  Base::template factorize<false>(a);
385  }
386 
388  Scalar determinant() const
389  {
390  Scalar detL = Base::m_matrix.diagonal().prod();
391  return numext::abs2(detL);
392  }
393 };
394 
413 template<typename _MatrixType, int _UpLo, typename _Ordering>
414  class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
415 {
416 public:
417  typedef _MatrixType MatrixType;
418  enum { UpLo = _UpLo };
419  typedef SimplicialCholeskyBase<SimplicialLDLT> Base;
420  typedef typename MatrixType::Scalar Scalar;
421  typedef typename MatrixType::RealScalar RealScalar;
422  typedef typename MatrixType::StorageIndex StorageIndex;
423  typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
424  typedef Matrix<Scalar,Dynamic,1> VectorType;
425  typedef internal::traits<SimplicialLDLT> Traits;
426  typedef typename Traits::MatrixL MatrixL;
427  typedef typename Traits::MatrixU MatrixU;
428 public:
430  SimplicialLDLT() : Base() {}
431 
433  explicit SimplicialLDLT(const MatrixType& matrix)
434  : Base(matrix) {}
435 
437  inline const VectorType vectorD() const {
438  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
439  return Base::m_diag;
440  }
442  inline const MatrixL matrixL() const {
443  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
444  return Traits::getL(Base::m_matrix);
445  }
446 
448  inline const MatrixU matrixU() const {
449  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
450  return Traits::getU(Base::m_matrix);
451  }
452 
454  SimplicialLDLT& compute(const MatrixType& matrix)
455  {
456  Base::template compute<true>(matrix);
457  return *this;
458  }
459 
466  void analyzePattern(const MatrixType& a)
467  {
468  Base::analyzePattern(a, true);
469  }
470 
477  void factorize(const MatrixType& a)
478  {
479  Base::template factorize<true>(a);
480  }
481 
483  Scalar determinant() const
484  {
485  return Base::m_diag.prod();
486  }
487 };
488 
495 template<typename _MatrixType, int _UpLo, typename _Ordering>
496  class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
497 {
498 public:
499  typedef _MatrixType MatrixType;
500  enum { UpLo = _UpLo };
501  typedef SimplicialCholeskyBase<SimplicialCholesky> Base;
502  typedef typename MatrixType::Scalar Scalar;
503  typedef typename MatrixType::RealScalar RealScalar;
504  typedef typename MatrixType::StorageIndex StorageIndex;
505  typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
506  typedef Matrix<Scalar,Dynamic,1> VectorType;
507  typedef internal::traits<SimplicialCholesky> Traits;
508  typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
509  typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
510  public:
511  SimplicialCholesky() : Base(), m_LDLT(true) {}
512 
513  explicit SimplicialCholesky(const MatrixType& matrix)
514  : Base(), m_LDLT(true)
515  {
516  compute(matrix);
517  }
518 
519  SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
520  {
521  switch(mode)
522  {
523  case SimplicialCholeskyLLT:
524  m_LDLT = false;
525  break;
526  case SimplicialCholeskyLDLT:
527  m_LDLT = true;
528  break;
529  default:
530  break;
531  }
532 
533  return *this;
534  }
535 
536  inline const VectorType vectorD() const {
537  eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
538  return Base::m_diag;
539  }
540  inline const CholMatrixType rawMatrix() const {
541  eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
542  return Base::m_matrix;
543  }
544 
546  SimplicialCholesky& compute(const MatrixType& matrix)
547  {
548  if(m_LDLT)
549  Base::template compute<true>(matrix);
550  else
551  Base::template compute<false>(matrix);
552  return *this;
553  }
554 
561  void analyzePattern(const MatrixType& a)
562  {
563  Base::analyzePattern(a, m_LDLT);
564  }
565 
572  void factorize(const MatrixType& a)
573  {
574  if(m_LDLT)
575  Base::template factorize<true>(a);
576  else
577  Base::template factorize<false>(a);
578  }
579 
581  template<typename Rhs,typename Dest>
582  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
583  {
584  eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
585  eigen_assert(Base::m_matrix.rows()==b.rows());
586 
587  if(Base::m_info!=Success)
588  return;
589 
590  if(Base::m_P.size()>0)
591  dest = Base::m_P * b;
592  else
593  dest = b;
594 
595  if(Base::m_matrix.nonZeros()>0) // otherwise L==I
596  {
597  if(m_LDLT)
598  LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
599  else
600  LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
601  }
602 
603  if(Base::m_diag.size()>0)
604  dest = Base::m_diag.asDiagonal().inverse() * dest;
605 
606  if (Base::m_matrix.nonZeros()>0) // otherwise I==I
607  {
608  if(m_LDLT)
609  LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
610  else
611  LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
612  }
613 
614  if(Base::m_P.size()>0)
615  dest = Base::m_Pinv * dest;
616  }
617 
619  template<typename Rhs,typename Dest>
620  void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const
621  {
622  internal::solve_sparse_through_dense_panels(*this, b, dest);
623  }
624 
625  Scalar determinant() const
626  {
627  if(m_LDLT)
628  {
629  return Base::m_diag.prod();
630  }
631  else
632  {
633  Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
634  return numext::abs2(detL);
635  }
636  }
637 
638  protected:
639  bool m_LDLT;
640 };
641 
642 template<typename Derived>
643 void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap)
644 {
645  eigen_assert(a.rows()==a.cols());
646  const Index size = a.rows();
647  pmat = &ap;
648  // Note that ordering methods compute the inverse permutation
649  if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value)
650  {
651  {
652  CholMatrixType C;
653  C = a.template selfadjointView<UpLo>();
654 
655  OrderingType ordering;
656  ordering(C,m_Pinv);
657  }
658 
659  if(m_Pinv.size()>0) m_P = m_Pinv.inverse();
660  else m_P.resize(0);
661 
662  ap.resize(size,size);
663  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
664  }
665  else
666  {
667  m_Pinv.resize(0);
668  m_P.resize(0);
669  if(int(UpLo)==int(Lower) || MatrixType::IsRowMajor)
670  {
671  // we have to transpose the lower part to to the upper one
672  ap.resize(size,size);
673  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>();
674  }
675  else
676  internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap);
677  }
678 }
679 
680 } // end namespace Eigen
681 
682 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H
SimplicialLDLT()
Definition: SimplicialCholesky.h:430
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:347
Index cols() const
Definition: SparseMatrix.h:132
SimplicialCholeskyBase()
Definition: SimplicialCholesky.h:79
Derived & setShift(const RealScalar &offset, const RealScalar &scale=1)
Definition: SimplicialCholesky.h:129
SimplicialLLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:359
Definition: Constants.h:196
A base class for sparse solvers.
Definition: SparseSolverBase.h:53
Index size() const
Definition: PermutationMatrix.h:110
Definition: LDLT.h:16
A direct sparse LDLT Cholesky factorizations without square root.
Definition: SimplicialCholesky.h:264
SimplicialLLT()
Definition: SimplicialCholesky.h:341
SimplicialLDLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:454
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:382
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:561
Definition: Constants.h:198
void compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:191
Scalar determinant() const
Definition: SimplicialCholesky.h:483
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationP() const
Definition: SimplicialCholesky.h:112
Definition: SimplicialCholesky.h:265
A direct sparse Cholesky factorizations.
Definition: SimplicialCholesky.h:57
SimplicialLLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:343
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:442
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:353
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:572
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SimplicialCholesky.h:104
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:448
Definition: Constants.h:424
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationPinv() const
Definition: SimplicialCholesky.h:117
Definition: Eigen_Colamd.h:54
SimplicialLDLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:433
Scalar determinant() const
Definition: SimplicialCholesky.h:388
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:477
Definition: SimplicialCholesky.h:241
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:186
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:371
A direct sparse LLT Cholesky factorizations.
Definition: SimplicialCholesky.h:263
const DiagonalWrapper< const Derived > asDiagonal() const
Definition: DiagonalMatrix.h:278
ComputationInfo
Definition: Constants.h:422
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const VectorType vectorD() const
Definition: SimplicialCholesky.h:437
Index rows() const
Definition: SparseMatrix.h:130
SimplicialCholesky & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:546
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:466