SimplicialCholesky.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-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 
18 };
19 
35 template<typename Derived>
37 {
38  public:
41  typedef typename MatrixType::Scalar Scalar;
42  typedef typename MatrixType::RealScalar RealScalar;
43  typedef typename MatrixType::Index Index;
46 
47  public:
48 
52  {}
53 
54  SimplicialCholeskyBase(const MatrixType& matrix)
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>
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>
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");
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 
231 
232  CholMatrixType m_matrix;
233  VectorType m_diag; // the diagonal coefficients (LDLT mode)
234  VectorXi m_parent; // elimination tree
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;
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;
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 };
307  typedef typename MatrixType::Scalar Scalar;
308  typedef typename MatrixType::RealScalar RealScalar;
309  typedef typename MatrixType::Index Index;
313  typedef typename Traits::MatrixL MatrixL;
314  typedef typename Traits::MatrixU MatrixU;
315 public:
317  SimplicialLLT() : Base() {}
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 };
395  typedef typename MatrixType::Scalar Scalar;
396  typedef typename MatrixType::RealScalar RealScalar;
397  typedef typename MatrixType::Index Index;
401  typedef typename Traits::MatrixL MatrixL;
402  typedef typename Traits::MatrixU MatrixU;
403 public:
405  SimplicialLDLT() : Base() {}
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 };
477  typedef typename MatrixType::Scalar Scalar;
478  typedef typename MatrixType::RealScalar RealScalar;
479  typedef typename MatrixType::Index Index;
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 
495  {
496  switch(mode)
497  {
499  m_LDLT = false;
500  break;
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>
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());
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>
639  : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
640 {
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>
652  : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
653 {
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
PermutationMatrix< Dynamic, Dynamic, Index > m_P
SimplicialCholesky(const MatrixType &matrix)
internal::traits< SimplicialLDLT< MatrixType, UpLo > > LDLTTraits
bool operator()(const Index &row, const Index &col, const Scalar &) const
void factorize(const MatrixType &a)
Index cols() const
Definition: SparseMatrix.h:121
const AdjointReturnType adjoint() const
SimplicialCholeskyBase< SimplicialCholesky > Base
SimplicialCholeskyBase< SimplicialLLT > Base
Derived & setShift(const RealScalar &offset, const RealScalar &scale=1)
void analyzePattern(const MatrixType &a)
const internal::solve_retval< SimplicialCholeskyBase, Rhs > solve(const MatrixBase< Rhs > &b) const
MatrixType::RealScalar RealScalar
internal::traits< SimplicialLLT > Traits
iterative scaling algorithm to equilibrate rows and column norms in matrices
Definition: matrix.hpp:471
SimplicialCholeskyBase(const MatrixType &matrix)
A direct sparse LDLT Cholesky factorizations without square root.
void analyzePattern(const MatrixType &a)
Index nonZeros() const
Definition: SparseMatrix.h:246
SimplicialCholesky & setMode(SimplicialCholeskyMode mode)
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs2_op< Scalar >, const Derived > abs2() const
void analyzePattern_preordered(const CholMatrixType &a, bool doLDLT)
SparseMatrix< Scalar, ColMajor, Index > CholMatrixType
SparseMatrix< Scalar, ColMajor, Index > CholMatrixType
Matrix< Scalar, Dynamic, 1 > VectorType
const VectorType vectorD() const
internal::traits< SimplicialLLT< MatrixType, UpLo > > LLTTraits
SimplicialLLT & compute(const MatrixType &matrix)
const internal::sparse_solve_retval< SimplicialCholeskyBase, Rhs > solve(const SparseMatrixBase< Rhs > &b) const
SparseMatrix< Scalar, ColMajor, Index > CholMatrixType
const MatrixU matrixU() const
Scalar determinant() const
SparseSymmetricPermutationProduct< Derived, Upper|Lower > twistedBy(const PermutationMatrix< Dynamic, Dynamic, Index > &perm) const
Matrix< Scalar, Dynamic, 1 > VectorType
const Derived & derived() const
void factorize(const MatrixType &a)
void compute(const MatrixType &matrix)
Base class of any sparse matrices or sparse expressions.
void ordering(const MatrixType &a, CholMatrixType &ap)
Matrix< Scalar, Dynamic, 1 > VectorType
void factorize(const MatrixType &a)
void factorize(const MatrixType &a)
A direct sparse Cholesky factorizations.
const MatrixL matrixL() const
const PermutationMatrix< Dynamic, Dynamic, Index > & permutationPinv() const
MatrixType::Scalar Scalar
internal::traits< SimplicialLDLT > Traits
void resize(Index rows, Index cols)
Definition: SparseMatrix.h:596
MatrixType::Index Index
void minimum_degree_ordering(SparseMatrix< Scalar, ColMajor, Index > &C, PermutationMatrix< Dynamic, Dynamic, Index > &perm)
Definition: Amd.h:91
void analyzePattern(const MatrixType &a, bool doLDLT)
ComputationInfo info() const
Reports whether previous computation was successful.
MatrixType::RealScalar RealScalar
#define EIGEN_MAKE_SPARSE_SOLVE_HELPERS(DecompositionType, Rhs)
Definition: SparseSolve.h:71
Transpose< PermutationBase > inverse() const
Scalar determinant() const
void rhs(const real_t *x, real_t *f)
void factorize_preordered(const CholMatrixType &a)
MatrixType::RealScalar RealScalar
const VectorType vectorD() const
MatrixType::Scalar Scalar
SimplicialCholesky & compute(const MatrixType &matrix)
SparseTriangularView< CholMatrixType, Eigen::Lower > MatrixL
MatrixType::Index Index
RowXpr row(Index i)
Definition: BlockMethods.h:725
internal::traits< Derived >::MatrixType MatrixType
const Derived & derived() const
SimplicialLDLT & compute(const MatrixType &matrix)
SparseTriangularView< typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper > MatrixU
void analyzePattern(const MatrixType &a)
internal::traits< SimplicialCholesky > Traits
SparseMatrix< Scalar, ColMajor, Index > CholMatrixType
void _solve(const MatrixBase< Rhs > &b, MatrixBase< Dest > &dest) const
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType, Rhs)
Definition: Solve.h:61
Expression of a diagonal/subdiagonal/superdiagonal in a matrix.
Definition: Diagonal.h:64
Matrix< Scalar, Dynamic, 1 > VectorType
MatrixType::RealScalar RealScalar
A direct sparse LLT Cholesky factorizations.
SparseTriangularView< typename CholMatrixType::AdjointReturnType, Eigen::Upper > MatrixU
SimplicialCholeskyBase< SimplicialLDLT > Base
ColXpr col(Index i)
Definition: BlockMethods.h:708
const PermutationMatrix< Dynamic, Dynamic, Index > & permutationP() const
const DiagonalWrapper< const Derived > asDiagonal() const
const MatrixL matrixL() const
#define eigen_assert(x)
void resize(Index newSize)
SimplicialLDLT(const MatrixType &matrix)
const CholMatrixType rawMatrix() const
const MatrixU matrixU() const
void _solve(const MatrixBase< Rhs > &b, MatrixBase< Dest > &dest) const
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
SparseTriangularView< CholMatrixType, Eigen::UnitLower > MatrixL
Index rows() const
Definition: SparseMatrix.h:119
PermutationMatrix< Dynamic, Dynamic, Index > m_Pinv
SimplicialLLT(const MatrixType &matrix)


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Mon Jun 10 2019 12:35:06