Eigen  3.2.92
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  enum {
75  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
76  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
77  };
78 
79  public:
80 
81  using Base::derived;
82 
85  : m_info(Success), m_shiftOffset(0), m_shiftScale(1)
86  {}
87 
88  explicit SimplicialCholeskyBase(const MatrixType& matrix)
89  : m_info(Success), m_shiftOffset(0), m_shiftScale(1)
90  {
91  derived().compute(matrix);
92  }
93 
95  {
96  }
97 
98  Derived& derived() { return *static_cast<Derived*>(this); }
99  const Derived& derived() const { return *static_cast<const Derived*>(this); }
100 
101  inline Index cols() const { return m_matrix.cols(); }
102  inline Index rows() const { return m_matrix.rows(); }
103 
110  {
111  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
112  return m_info;
113  }
114 
118  { return m_P; }
119 
123  { return m_Pinv; }
124 
134  Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1)
135  {
136  m_shiftOffset = offset;
137  m_shiftScale = scale;
138  return derived();
139  }
140 
141 #ifndef EIGEN_PARSED_BY_DOXYGEN
142 
143  template<typename Stream>
144  void dumpMemory(Stream& s)
145  {
146  int total = 0;
147  s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
148  s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
149  s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
150  s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
151  s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
152  s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
153  s << " TOTAL: " << (total>> 20) << "Mb" << "\n";
154  }
155 
157  template<typename Rhs,typename Dest>
158  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
159  {
160  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
161  eigen_assert(m_matrix.rows()==b.rows());
162 
163  if(m_info!=Success)
164  return;
165 
166  if(m_P.size()>0)
167  dest = m_P * b;
168  else
169  dest = b;
170 
171  if(m_matrix.nonZeros()>0) // otherwise L==I
172  derived().matrixL().solveInPlace(dest);
173 
174  if(m_diag.size()>0)
175  dest = m_diag.asDiagonal().inverse() * dest;
176 
177  if (m_matrix.nonZeros()>0) // otherwise U==I
178  derived().matrixU().solveInPlace(dest);
179 
180  if(m_P.size()>0)
181  dest = m_Pinv * dest;
182  }
183 
184  template<typename Rhs,typename Dest>
185  void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const
186  {
187  internal::solve_sparse_through_dense_panels(derived(), b, dest);
188  }
189 
190 #endif // EIGEN_PARSED_BY_DOXYGEN
191 
192  protected:
193 
195  template<bool DoLDLT>
196  void compute(const MatrixType& matrix)
197  {
198  eigen_assert(matrix.rows()==matrix.cols());
199  Index size = matrix.cols();
200  CholMatrixType tmp(size,size);
201  ConstCholMatrixPtr pmat;
202  ordering(matrix, pmat, tmp);
203  analyzePattern_preordered(*pmat, DoLDLT);
204  factorize_preordered<DoLDLT>(*pmat);
205  }
206 
207  template<bool DoLDLT>
208  void factorize(const MatrixType& a)
209  {
210  eigen_assert(a.rows()==a.cols());
211  Index size = a.cols();
212  CholMatrixType tmp(size,size);
213  ConstCholMatrixPtr pmat;
214 
215  if(m_P.size()==0 && (UpLo&Upper)==Upper)
216  {
217  // If there is no ordering, try to directly use the input matrix without any copy
218  internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp);
219  }
220  else
221  {
222  tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
223  pmat = &tmp;
224  }
225 
226  factorize_preordered<DoLDLT>(*pmat);
227  }
228 
229  template<bool DoLDLT>
230  void factorize_preordered(const CholMatrixType& a);
231 
232  void analyzePattern(const MatrixType& a, bool doLDLT)
233  {
234  eigen_assert(a.rows()==a.cols());
235  Index size = a.cols();
236  CholMatrixType tmp(size,size);
237  ConstCholMatrixPtr pmat;
238  ordering(a, pmat, tmp);
239  analyzePattern_preordered(*pmat,doLDLT);
240  }
241  void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
242 
243  void ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap);
244 
246  struct keep_diag {
247  inline bool operator() (const Index& row, const Index& col, const Scalar&) const
248  {
249  return row!=col;
250  }
251  };
252 
253  mutable ComputationInfo m_info;
254  bool m_factorizationIsOk;
255  bool m_analysisIsOk;
256 
257  CholMatrixType m_matrix;
258  VectorType m_diag; // the diagonal coefficients (LDLT mode)
259  VectorI m_parent; // elimination tree
260  VectorI m_nonZerosPerCol;
262  PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_Pinv; // the inverse permutation
263 
264  RealScalar m_shiftOffset;
265  RealScalar m_shiftScale;
266 };
267 
268 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLLT;
269 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLDLT;
270 template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialCholesky;
271 
272 namespace internal {
273 
274 template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
275 {
276  typedef _MatrixType MatrixType;
277  typedef _Ordering OrderingType;
278  enum { UpLo = _UpLo };
279  typedef typename MatrixType::Scalar Scalar;
280  typedef typename MatrixType::StorageIndex StorageIndex;
281  typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
284  static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); }
285  static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
286 };
287 
288 template<typename _MatrixType,int _UpLo, typename _Ordering> struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
289 {
290  typedef _MatrixType MatrixType;
291  typedef _Ordering OrderingType;
292  enum { UpLo = _UpLo };
293  typedef typename MatrixType::Scalar Scalar;
294  typedef typename MatrixType::StorageIndex StorageIndex;
295  typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
296  typedef TriangularView<const CholMatrixType, Eigen::UnitLower> MatrixL;
297  typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
298  static inline MatrixL getL(const MatrixType& m) { return MatrixL(m); }
299  static inline MatrixU getU(const MatrixType& m) { return MatrixU(m.adjoint()); }
300 };
301 
302 template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
303 {
304  typedef _MatrixType MatrixType;
305  typedef _Ordering OrderingType;
306  enum { UpLo = _UpLo };
307 };
308 
309 }
310 
331 template<typename _MatrixType, int _UpLo, typename _Ordering>
332  class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
333 {
334 public:
335  typedef _MatrixType MatrixType;
336  enum { UpLo = _UpLo };
337  typedef SimplicialCholeskyBase<SimplicialLLT> Base;
338  typedef typename MatrixType::Scalar Scalar;
339  typedef typename MatrixType::RealScalar RealScalar;
340  typedef typename MatrixType::StorageIndex StorageIndex;
341  typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
342  typedef Matrix<Scalar,Dynamic,1> VectorType;
343  typedef internal::traits<SimplicialLLT> Traits;
344  typedef typename Traits::MatrixL MatrixL;
345  typedef typename Traits::MatrixU MatrixU;
346 public:
348  SimplicialLLT() : Base() {}
350  explicit SimplicialLLT(const MatrixType& matrix)
351  : Base(matrix) {}
352 
354  inline const MatrixL matrixL() const {
355  eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
356  return Traits::getL(Base::m_matrix);
357  }
358 
360  inline const MatrixU matrixU() const {
361  eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
362  return Traits::getU(Base::m_matrix);
363  }
364 
366  SimplicialLLT& compute(const MatrixType& matrix)
367  {
368  Base::template compute<false>(matrix);
369  return *this;
370  }
371 
378  void analyzePattern(const MatrixType& a)
379  {
380  Base::analyzePattern(a, false);
381  }
382 
389  void factorize(const MatrixType& a)
390  {
391  Base::template factorize<false>(a);
392  }
393 
395  Scalar determinant() const
396  {
397  Scalar detL = Base::m_matrix.diagonal().prod();
398  return numext::abs2(detL);
399  }
400 };
401 
422 template<typename _MatrixType, int _UpLo, typename _Ordering>
423  class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
424 {
425 public:
426  typedef _MatrixType MatrixType;
427  enum { UpLo = _UpLo };
428  typedef SimplicialCholeskyBase<SimplicialLDLT> Base;
429  typedef typename MatrixType::Scalar Scalar;
430  typedef typename MatrixType::RealScalar RealScalar;
431  typedef typename MatrixType::StorageIndex StorageIndex;
432  typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
433  typedef Matrix<Scalar,Dynamic,1> VectorType;
434  typedef internal::traits<SimplicialLDLT> Traits;
435  typedef typename Traits::MatrixL MatrixL;
436  typedef typename Traits::MatrixU MatrixU;
437 public:
439  SimplicialLDLT() : Base() {}
440 
442  explicit SimplicialLDLT(const MatrixType& matrix)
443  : Base(matrix) {}
444 
446  inline const VectorType vectorD() const {
447  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
448  return Base::m_diag;
449  }
451  inline const MatrixL matrixL() const {
452  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
453  return Traits::getL(Base::m_matrix);
454  }
455 
457  inline const MatrixU matrixU() const {
458  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
459  return Traits::getU(Base::m_matrix);
460  }
461 
463  SimplicialLDLT& compute(const MatrixType& matrix)
464  {
465  Base::template compute<true>(matrix);
466  return *this;
467  }
468 
475  void analyzePattern(const MatrixType& a)
476  {
477  Base::analyzePattern(a, true);
478  }
479 
486  void factorize(const MatrixType& a)
487  {
488  Base::template factorize<true>(a);
489  }
490 
492  Scalar determinant() const
493  {
494  return Base::m_diag.prod();
495  }
496 };
497 
504 template<typename _MatrixType, int _UpLo, typename _Ordering>
505  class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
506 {
507 public:
508  typedef _MatrixType MatrixType;
509  enum { UpLo = _UpLo };
510  typedef SimplicialCholeskyBase<SimplicialCholesky> Base;
511  typedef typename MatrixType::Scalar Scalar;
512  typedef typename MatrixType::RealScalar RealScalar;
513  typedef typename MatrixType::StorageIndex StorageIndex;
514  typedef SparseMatrix<Scalar,ColMajor,StorageIndex> CholMatrixType;
515  typedef Matrix<Scalar,Dynamic,1> VectorType;
516  typedef internal::traits<SimplicialCholesky> Traits;
517  typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
518  typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
519  public:
520  SimplicialCholesky() : Base(), m_LDLT(true) {}
521 
522  explicit SimplicialCholesky(const MatrixType& matrix)
523  : Base(), m_LDLT(true)
524  {
525  compute(matrix);
526  }
527 
528  SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
529  {
530  switch(mode)
531  {
532  case SimplicialCholeskyLLT:
533  m_LDLT = false;
534  break;
535  case SimplicialCholeskyLDLT:
536  m_LDLT = true;
537  break;
538  default:
539  break;
540  }
541 
542  return *this;
543  }
544 
545  inline const VectorType vectorD() const {
546  eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
547  return Base::m_diag;
548  }
549  inline const CholMatrixType rawMatrix() const {
550  eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
551  return Base::m_matrix;
552  }
553 
555  SimplicialCholesky& compute(const MatrixType& matrix)
556  {
557  if(m_LDLT)
558  Base::template compute<true>(matrix);
559  else
560  Base::template compute<false>(matrix);
561  return *this;
562  }
563 
570  void analyzePattern(const MatrixType& a)
571  {
572  Base::analyzePattern(a, m_LDLT);
573  }
574 
581  void factorize(const MatrixType& a)
582  {
583  if(m_LDLT)
584  Base::template factorize<true>(a);
585  else
586  Base::template factorize<false>(a);
587  }
588 
590  template<typename Rhs,typename Dest>
591  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
592  {
593  eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
594  eigen_assert(Base::m_matrix.rows()==b.rows());
595 
596  if(Base::m_info!=Success)
597  return;
598 
599  if(Base::m_P.size()>0)
600  dest = Base::m_P * b;
601  else
602  dest = b;
603 
604  if(Base::m_matrix.nonZeros()>0) // otherwise L==I
605  {
606  if(m_LDLT)
607  LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
608  else
609  LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
610  }
611 
612  if(Base::m_diag.size()>0)
613  dest = Base::m_diag.asDiagonal().inverse() * dest;
614 
615  if (Base::m_matrix.nonZeros()>0) // otherwise I==I
616  {
617  if(m_LDLT)
618  LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
619  else
620  LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
621  }
622 
623  if(Base::m_P.size()>0)
624  dest = Base::m_Pinv * dest;
625  }
626 
628  template<typename Rhs,typename Dest>
629  void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const
630  {
631  internal::solve_sparse_through_dense_panels(*this, b, dest);
632  }
633 
634  Scalar determinant() const
635  {
636  if(m_LDLT)
637  {
638  return Base::m_diag.prod();
639  }
640  else
641  {
642  Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
643  return numext::abs2(detL);
644  }
645  }
646 
647  protected:
648  bool m_LDLT;
649 };
650 
651 template<typename Derived>
652 void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap)
653 {
654  eigen_assert(a.rows()==a.cols());
655  const Index size = a.rows();
656  pmat = &ap;
657  // Note that ordering methods compute the inverse permutation
658  if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value)
659  {
660  {
661  CholMatrixType C;
662  C = a.template selfadjointView<UpLo>();
663 
664  OrderingType ordering;
665  ordering(C,m_Pinv);
666  }
667 
668  if(m_Pinv.size()>0) m_P = m_Pinv.inverse();
669  else m_P.resize(0);
670 
671  ap.resize(size,size);
672  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
673  }
674  else
675  {
676  m_Pinv.resize(0);
677  m_P.resize(0);
678  if(int(UpLo)==int(Lower) || MatrixType::IsRowMajor)
679  {
680  // we have to transpose the lower part to to the upper one
681  ap.resize(size,size);
682  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>();
683  }
684  else
685  internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap);
686  }
687 }
688 
689 } // end namespace Eigen
690 
691 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H
SimplicialLDLT()
Definition: SimplicialCholesky.h:439
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:354
Index cols() const
Definition: SparseMatrix.h:133
SimplicialCholeskyBase()
Definition: SimplicialCholesky.h:84
Derived & setShift(const RealScalar &offset, const RealScalar &scale=1)
Definition: SimplicialCholesky.h:134
SimplicialLLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:366
A base class for sparse solvers.
Definition: SparseSolverBase.h:53
Index size() const
Definition: PermutationMatrix.h:109
Definition: LDLT.h:16
A direct sparse LDLT Cholesky factorizations without square root.
Definition: SimplicialCholesky.h:269
SimplicialLLT()
Definition: SimplicialCholesky.h:348
SimplicialLDLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:463
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:389
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:570
void compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:196
Definition: Constants.h:204
Scalar determinant() const
Definition: SimplicialCholesky.h:492
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationP() const
Definition: SimplicialCholesky.h:117
Definition: SimplicialCholesky.h:270
A direct sparse Cholesky factorizations.
Definition: SimplicialCholesky.h:57
SimplicialLLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:350
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:451
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:360
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:581
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SimplicialCholesky.h:109
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:457
Definition: Constants.h:432
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationPinv() const
Definition: SimplicialCholesky.h:122
Definition: Eigen_Colamd.h:54
SimplicialLDLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:442
Scalar determinant() const
Definition: SimplicialCholesky.h:395
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:486
Definition: SimplicialCholesky.h:246
Expression of a triangular part in a matrix.
Definition: TriangularMatrix.h:186
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:378
A direct sparse LLT Cholesky factorizations.
Definition: SimplicialCholesky.h:268
const DiagonalWrapper< const Derived > asDiagonal() const
Definition: DiagonalMatrix.h:278
Definition: Constants.h:206
ComputationInfo
Definition: Constants.h:430
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const VectorType vectorD() const
Definition: SimplicialCholesky.h:446
Index rows() const
Definition: SparseMatrix.h:131
SimplicialCholesky & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:555
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:475