Eigen  3.4.90 (git rev 67eeba6e720c5745abc77ae6c92ce0a44aa7b7ae)
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 #include "./InternalHeaderCheck.h"
14 
15 namespace Eigen {
16 
17 enum SimplicialCholeskyMode {
18  SimplicialCholeskyLLT,
19  SimplicialCholeskyLDLT
20 };
21 
22 namespace internal {
23  template<typename CholMatrixType, typename InputMatrixType>
24  struct simplicial_cholesky_grab_input {
25  typedef CholMatrixType const * ConstCholMatrixPtr;
26  static void run(const InputMatrixType& input, ConstCholMatrixPtr &pmat, CholMatrixType &tmp)
27  {
28  tmp = input;
29  pmat = &tmp;
30  }
31  };
32 
33  template<typename MatrixType>
34  struct simplicial_cholesky_grab_input<MatrixType,MatrixType> {
35  typedef MatrixType const * ConstMatrixPtr;
36  static void run(const MatrixType& input, ConstMatrixPtr &pmat, MatrixType &/*tmp*/)
37  {
38  pmat = &input;
39  }
40  };
41 } // end namespace internal
42 
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),
86  m_factorizationIsOk(false),
87  m_analysisIsOk(false),
88  m_shiftOffset(0),
89  m_shiftScale(1)
90  {}
91 
92  explicit SimplicialCholeskyBase(const MatrixType& matrix)
93  : m_info(Success),
94  m_factorizationIsOk(false),
95  m_analysisIsOk(false),
96  m_shiftOffset(0),
97  m_shiftScale(1)
98  {
99  derived().compute(matrix);
100  }
101 
103  {
104  }
105 
106  Derived& derived() { return *static_cast<Derived*>(this); }
107  const Derived& derived() const { return *static_cast<const Derived*>(this); }
108 
109  inline Index cols() const { return m_matrix.cols(); }
110  inline Index rows() const { return m_matrix.rows(); }
111 
118  {
119  eigen_assert(m_isInitialized && "Decomposition is not initialized.");
120  return m_info;
121  }
122 
126  { return m_P; }
127 
131  { return m_Pinv; }
132 
142  Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1)
143  {
144  m_shiftOffset = offset;
145  m_shiftScale = scale;
146  return derived();
147  }
148 
149 #ifndef EIGEN_PARSED_BY_DOXYGEN
151  template<typename Stream>
152  void dumpMemory(Stream& s)
153  {
154  int total = 0;
155  s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
156  s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
157  s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
158  s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
159  s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
160  s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
161  s << " TOTAL: " << (total>> 20) << "Mb" << "\n";
162  }
163 
165  template<typename Rhs,typename Dest>
166  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
167  {
168  eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
169  eigen_assert(m_matrix.rows()==b.rows());
170 
171  if(m_info!=Success)
172  return;
173 
174  if(m_P.size()>0)
175  dest = m_P * b;
176  else
177  dest = b;
178 
179  if(m_matrix.nonZeros()>0) // otherwise L==I
180  derived().matrixL().solveInPlace(dest);
181 
182  if(m_diag.size()>0)
183  dest = m_diag.asDiagonal().inverse() * dest;
184 
185  if (m_matrix.nonZeros()>0) // otherwise U==I
186  derived().matrixU().solveInPlace(dest);
187 
188  if(m_P.size()>0)
189  dest = m_Pinv * dest;
190  }
191 
192  template<typename Rhs,typename Dest>
193  void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const
194  {
195  internal::solve_sparse_through_dense_panels(derived(), b, dest);
196  }
197 
198 #endif // EIGEN_PARSED_BY_DOXYGEN
199 
200  protected:
201 
203  template<bool DoLDLT>
204  void compute(const MatrixType& matrix)
205  {
206  eigen_assert(matrix.rows()==matrix.cols());
207  Index size = matrix.cols();
208  CholMatrixType tmp(size,size);
209  ConstCholMatrixPtr pmat;
210  ordering(matrix, pmat, tmp);
211  analyzePattern_preordered(*pmat, DoLDLT);
212  factorize_preordered<DoLDLT>(*pmat);
213  }
214 
215  template<bool DoLDLT>
216  void factorize(const MatrixType& a)
217  {
218  eigen_assert(a.rows()==a.cols());
219  Index size = a.cols();
220  CholMatrixType tmp(size,size);
221  ConstCholMatrixPtr pmat;
222 
223  if(m_P.size() == 0 && (int(UpLo) & int(Upper)) == Upper)
224  {
225  // If there is no ordering, try to directly use the input matrix without any copy
226  internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp);
227  }
228  else
229  {
230  tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
231  pmat = &tmp;
232  }
233 
234  factorize_preordered<DoLDLT>(*pmat);
235  }
236 
237  template<bool DoLDLT>
238  void factorize_preordered(const CholMatrixType& a);
239 
240  void analyzePattern(const MatrixType& a, bool doLDLT)
241  {
242  eigen_assert(a.rows()==a.cols());
243  Index size = a.cols();
244  CholMatrixType tmp(size,size);
245  ConstCholMatrixPtr pmat;
246  ordering(a, pmat, tmp);
247  analyzePattern_preordered(*pmat,doLDLT);
248  }
249  void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
250 
251  void ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap);
252 
254  struct keep_diag {
255  inline bool operator() (const Index& row, const Index& col, const Scalar&) const
256  {
257  return row!=col;
258  }
259  };
260 
261  mutable ComputationInfo m_info;
262  bool m_factorizationIsOk;
263  bool m_analysisIsOk;
264 
265  CholMatrixType m_matrix;
266  VectorType m_diag; // the diagonal coefficients (LDLT mode)
267  VectorI m_parent; // elimination tree
268  VectorI m_nonZerosPerCol;
270  PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_Pinv; // the inverse permutation
271 
272  RealScalar m_shiftOffset;
273  RealScalar m_shiftScale;
274 };
275 
276 template<typename MatrixType_, int UpLo_ = Lower, typename Ordering_ = AMDOrdering<typename MatrixType_::StorageIndex> > class SimplicialLLT;
277 template<typename MatrixType_, int UpLo_ = Lower, typename Ordering_ = AMDOrdering<typename MatrixType_::StorageIndex> > class SimplicialLDLT;
278 template<typename MatrixType_, int UpLo_ = Lower, typename Ordering_ = AMDOrdering<typename MatrixType_::StorageIndex> > class SimplicialCholesky;
279 
280 namespace internal {
281 
282 template<typename MatrixType_, int UpLo_, typename Ordering_> struct traits<SimplicialLLT<MatrixType_,UpLo_,Ordering_> >
283 {
284  typedef MatrixType_ MatrixType;
285  typedef Ordering_ OrderingType;
286  enum { UpLo = UpLo_ };
287  typedef typename MatrixType::Scalar Scalar;
288  typedef typename MatrixType::StorageIndex StorageIndex;
289  typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
290  typedef TriangularView<const CholMatrixType, Eigen::Lower> MatrixL;
291  typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
292  static inline MatrixL getL(const CholMatrixType& m) { return MatrixL(m); }
293  static inline MatrixU getU(const CholMatrixType& m) { return MatrixU(m.adjoint()); }
294 };
295 
296 template<typename MatrixType_,int UpLo_, typename Ordering_> struct traits<SimplicialLDLT<MatrixType_,UpLo_,Ordering_> >
297 {
298  typedef MatrixType_ MatrixType;
299  typedef Ordering_ OrderingType;
300  enum { UpLo = UpLo_ };
301  typedef typename MatrixType::Scalar Scalar;
302  typedef typename MatrixType::StorageIndex StorageIndex;
303  typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
304  typedef TriangularView<const CholMatrixType, Eigen::UnitLower> MatrixL;
305  typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
306  static inline MatrixL getL(const CholMatrixType& m) { return MatrixL(m); }
307  static inline MatrixU getU(const CholMatrixType& m) { return MatrixU(m.adjoint()); }
308 };
309 
310 template<typename MatrixType_, int UpLo_, typename Ordering_> struct traits<SimplicialCholesky<MatrixType_,UpLo_,Ordering_> >
311 {
312  typedef MatrixType_ MatrixType;
313  typedef Ordering_ OrderingType;
314  enum { UpLo = UpLo_ };
315 };
316 
317 }
318 
339 template<typename MatrixType_, int UpLo_, typename Ordering_>
340  class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<MatrixType_,UpLo_,Ordering_> >
341 {
342 public:
343  typedef MatrixType_ MatrixType;
344  enum { UpLo = UpLo_ };
346  typedef typename MatrixType::Scalar Scalar;
347  typedef typename MatrixType::RealScalar RealScalar;
348  typedef typename MatrixType::StorageIndex StorageIndex;
351  typedef internal::traits<SimplicialLLT> Traits;
352  typedef typename Traits::MatrixL MatrixL;
353  typedef typename Traits::MatrixU MatrixU;
354 public:
358  explicit SimplicialLLT(const MatrixType& matrix)
359  : Base(matrix) {}
360 
362  inline const MatrixL matrixL() const {
363  eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
364  return Traits::getL(Base::m_matrix);
365  }
366 
368  inline const MatrixU matrixU() const {
369  eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
370  return Traits::getU(Base::m_matrix);
371  }
372 
374  SimplicialLLT& compute(const MatrixType& matrix)
375  {
376  Base::template compute<false>(matrix);
377  return *this;
378  }
379 
386  void analyzePattern(const MatrixType& a)
387  {
388  Base::analyzePattern(a, false);
389  }
390 
397  void factorize(const MatrixType& a)
398  {
399  Base::template factorize<false>(a);
400  }
401 
403  Scalar determinant() const
404  {
405  Scalar detL = Base::m_matrix.diagonal().prod();
406  return numext::abs2(detL);
407  }
408 };
409 
430 template<typename MatrixType_, int UpLo_, typename Ordering_>
431  class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<MatrixType_,UpLo_,Ordering_> >
432 {
433 public:
434  typedef MatrixType_ MatrixType;
435  enum { UpLo = UpLo_ };
437  typedef typename MatrixType::Scalar Scalar;
438  typedef typename MatrixType::RealScalar RealScalar;
439  typedef typename MatrixType::StorageIndex StorageIndex;
442  typedef internal::traits<SimplicialLDLT> Traits;
443  typedef typename Traits::MatrixL MatrixL;
444  typedef typename Traits::MatrixU MatrixU;
445 public:
448 
450  explicit SimplicialLDLT(const MatrixType& matrix)
451  : Base(matrix) {}
452 
454  inline const VectorType vectorD() const {
455  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
456  return Base::m_diag;
457  }
459  inline const MatrixL matrixL() const {
460  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
461  return Traits::getL(Base::m_matrix);
462  }
463 
465  inline const MatrixU matrixU() const {
466  eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
467  return Traits::getU(Base::m_matrix);
468  }
469 
471  SimplicialLDLT& compute(const MatrixType& matrix)
472  {
473  Base::template compute<true>(matrix);
474  return *this;
475  }
476 
483  void analyzePattern(const MatrixType& a)
484  {
485  Base::analyzePattern(a, true);
486  }
487 
494  void factorize(const MatrixType& a)
495  {
496  Base::template factorize<true>(a);
497  }
498 
500  Scalar determinant() const
501  {
502  return Base::m_diag.prod();
503  }
504 };
505 
512 template<typename MatrixType_, int UpLo_, typename Ordering_>
513  class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<MatrixType_,UpLo_,Ordering_> >
514 {
515 public:
516  typedef MatrixType_ MatrixType;
517  enum { UpLo = UpLo_ };
519  typedef typename MatrixType::Scalar Scalar;
520  typedef typename MatrixType::RealScalar RealScalar;
521  typedef typename MatrixType::StorageIndex StorageIndex;
524  typedef internal::traits<SimplicialCholesky> Traits;
525  typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
526  typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
527  public:
528  SimplicialCholesky() : Base(), m_LDLT(true) {}
529 
530  explicit SimplicialCholesky(const MatrixType& matrix)
531  : Base(), m_LDLT(true)
532  {
533  compute(matrix);
534  }
535 
536  SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
537  {
538  switch(mode)
539  {
540  case SimplicialCholeskyLLT:
541  m_LDLT = false;
542  break;
543  case SimplicialCholeskyLDLT:
544  m_LDLT = true;
545  break;
546  default:
547  break;
548  }
549 
550  return *this;
551  }
552 
553  inline const VectorType vectorD() const {
554  eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
555  return Base::m_diag;
556  }
557  inline const CholMatrixType rawMatrix() const {
558  eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
559  return Base::m_matrix;
560  }
561 
563  SimplicialCholesky& compute(const MatrixType& matrix)
564  {
565  if(m_LDLT)
566  Base::template compute<true>(matrix);
567  else
568  Base::template compute<false>(matrix);
569  return *this;
570  }
571 
578  void analyzePattern(const MatrixType& a)
579  {
580  Base::analyzePattern(a, m_LDLT);
581  }
582 
589  void factorize(const MatrixType& a)
590  {
591  if(m_LDLT)
592  Base::template factorize<true>(a);
593  else
594  Base::template factorize<false>(a);
595  }
596 
598  template<typename Rhs,typename Dest>
599  void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
600  {
601  eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
602  eigen_assert(Base::m_matrix.rows()==b.rows());
603 
604  if(Base::m_info!=Success)
605  return;
606 
607  if(Base::m_P.size()>0)
608  dest = Base::m_P * b;
609  else
610  dest = b;
611 
612  if(Base::m_matrix.nonZeros()>0) // otherwise L==I
613  {
614  if(m_LDLT)
615  LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
616  else
617  LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
618  }
619 
620  if(Base::m_diag.size()>0)
621  dest = Base::m_diag.real().asDiagonal().inverse() * dest;
622 
623  if (Base::m_matrix.nonZeros()>0) // otherwise I==I
624  {
625  if(m_LDLT)
626  LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
627  else
628  LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
629  }
630 
631  if(Base::m_P.size()>0)
632  dest = Base::m_Pinv * dest;
633  }
634 
636  template<typename Rhs,typename Dest>
637  void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const
638  {
639  internal::solve_sparse_through_dense_panels(*this, b, dest);
640  }
641 
642  Scalar determinant() const
643  {
644  if(m_LDLT)
645  {
646  return Base::m_diag.prod();
647  }
648  else
649  {
650  Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
651  return numext::abs2(detL);
652  }
653  }
654 
655  protected:
656  bool m_LDLT;
657 };
658 
659 template<typename Derived>
660 void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap)
661 {
662  eigen_assert(a.rows()==a.cols());
663  const Index size = a.rows();
664  pmat = &ap;
665  // Note that ordering methods compute the inverse permutation
666  if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value)
667  {
668  {
669  CholMatrixType C;
670  C = a.template selfadjointView<UpLo>();
671 
672  OrderingType ordering;
673  ordering(C,m_Pinv);
674  }
675 
676  if(m_Pinv.size()>0) m_P = m_Pinv.inverse();
677  else m_P.resize(0);
678 
679  ap.resize(size,size);
680  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
681  }
682  else
683  {
684  m_Pinv.resize(0);
685  m_P.resize(0);
686  if(int(UpLo)==int(Lower) || MatrixType::IsRowMajor)
687  {
688  // we have to transpose the lower part to to the upper one
689  ap.resize(size,size);
690  ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>();
691  }
692  else
693  internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap);
694  }
695 }
696 
697 } // end namespace Eigen
698 
699 #endif // EIGEN_SIMPLICIAL_CHOLESKY_H
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:62
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
const DiagonalWrapper< const Derived > asDiagonal() const
Definition: DiagonalMatrix.h:327
Index size() const
Definition: PermutationMatrix.h:99
A base class for direct sparse Cholesky factorizations.
Definition: SimplicialCholesky.h:58
SimplicialCholeskyBase()
Definition: SimplicialCholesky.h:84
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationPinv() const
Definition: SimplicialCholesky.h:130
Derived & setShift(const RealScalar &offset, const RealScalar &scale=1)
Definition: SimplicialCholesky.h:142
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SimplicialCholesky.h:117
void compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:204
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationP() const
Definition: SimplicialCholesky.h:125
Definition: SimplicialCholesky.h:514
SimplicialCholesky & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:563
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:589
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:578
A direct sparse LDLT Cholesky factorizations without square root.
Definition: SimplicialCholesky.h:432
Scalar determinant() const
Definition: SimplicialCholesky.h:500
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:459
SimplicialLDLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:450
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:465
SimplicialLDLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:471
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:483
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:494
const VectorType vectorD() const
Definition: SimplicialCholesky.h:454
SimplicialLDLT()
Definition: SimplicialCholesky.h:447
A direct sparse LLT Cholesky factorizations.
Definition: SimplicialCholesky.h:341
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:362
Scalar determinant() const
Definition: SimplicialCholesky.h:403
SimplicialLLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:374
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:368
SimplicialLLT()
Definition: SimplicialCholesky.h:356
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:386
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:397
SimplicialLLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:358
Index nonZeros() const
Definition: SparseCompressedBase.h:58
Index cols() const
Definition: SparseMatrix.h:142
Index rows() const
Definition: SparseMatrix.h:140
A base class for sparse solvers.
Definition: SparseSolverBase.h:70
ComputationInfo
Definition: Constants.h:442
@ Lower
Definition: Constants.h:211
@ Upper
Definition: Constants.h:213
@ Success
Definition: Constants.h:444
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
Definition: SimplicialCholesky.h:254