Please, help us to better know about our user community by answering the following short survey: https://forms.gle/wpyrxWi18ox9Z5ae9
Eigen  3.4.0
 
Loading...
Searching...
No Matches
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
13namespace Eigen {
14
15enum SimplicialCholeskyMode {
16 SimplicialCholeskyLLT,
17 SimplicialCholeskyLDLT
18};
19
20namespace 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
54template<typename Derived>
56{
58 using Base::m_isInitialized;
59
60 public:
61 typedef typename internal::traits<Derived>::MatrixType MatrixType;
62 typedef typename internal::traits<Derived>::OrderingType OrderingType;
63 enum { UpLo = internal::traits<Derived>::UpLo };
64 typedef typename MatrixType::Scalar Scalar;
65 typedef typename MatrixType::RealScalar RealScalar;
66 typedef typename MatrixType::StorageIndex StorageIndex;
68 typedef CholMatrixType const * ConstCholMatrixPtr;
71
72 enum {
73 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
74 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
75 };
76
77 public:
78
79 using Base::derived;
80
83 : m_info(Success),
84 m_factorizationIsOk(false),
85 m_analysisIsOk(false),
86 m_shiftOffset(0),
87 m_shiftScale(1)
88 {}
89
90 explicit SimplicialCholeskyBase(const MatrixType& matrix)
91 : m_info(Success),
92 m_factorizationIsOk(false),
93 m_analysisIsOk(false),
94 m_shiftOffset(0),
95 m_shiftScale(1)
96 {
97 derived().compute(matrix);
98 }
99
100 ~SimplicialCholeskyBase()
101 {
102 }
103
104 Derived& derived() { return *static_cast<Derived*>(this); }
105 const Derived& derived() const { return *static_cast<const Derived*>(this); }
106
107 inline Index cols() const { return m_matrix.cols(); }
108 inline Index rows() const { return m_matrix.rows(); }
109
116 {
117 eigen_assert(m_isInitialized && "Decomposition is not initialized.");
118 return m_info;
119 }
120
124 { return m_P; }
125
129 { return m_Pinv; }
130
140 Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1)
141 {
142 m_shiftOffset = offset;
143 m_shiftScale = scale;
144 return derived();
145 }
146
147#ifndef EIGEN_PARSED_BY_DOXYGEN
149 template<typename Stream>
150 void dumpMemory(Stream& s)
151 {
152 int total = 0;
153 s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
154 s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
155 s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
156 s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
157 s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
158 s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
159 s << " TOTAL: " << (total>> 20) << "Mb" << "\n";
160 }
161
163 template<typename Rhs,typename Dest>
164 void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
165 {
166 eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
167 eigen_assert(m_matrix.rows()==b.rows());
168
169 if(m_info!=Success)
170 return;
171
172 if(m_P.size()>0)
173 dest = m_P * b;
174 else
175 dest = b;
176
177 if(m_matrix.nonZeros()>0) // otherwise L==I
178 derived().matrixL().solveInPlace(dest);
179
180 if(m_diag.size()>0)
181 dest = m_diag.asDiagonal().inverse() * dest;
182
183 if (m_matrix.nonZeros()>0) // otherwise U==I
184 derived().matrixU().solveInPlace(dest);
185
186 if(m_P.size()>0)
187 dest = m_Pinv * dest;
188 }
189
190 template<typename Rhs,typename Dest>
191 void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const
192 {
193 internal::solve_sparse_through_dense_panels(derived(), b, dest);
194 }
195
196#endif // EIGEN_PARSED_BY_DOXYGEN
197
198 protected:
199
201 template<bool DoLDLT>
202 void compute(const MatrixType& matrix)
203 {
204 eigen_assert(matrix.rows()==matrix.cols());
205 Index size = matrix.cols();
206 CholMatrixType tmp(size,size);
207 ConstCholMatrixPtr pmat;
208 ordering(matrix, pmat, tmp);
209 analyzePattern_preordered(*pmat, DoLDLT);
210 factorize_preordered<DoLDLT>(*pmat);
211 }
212
213 template<bool DoLDLT>
214 void factorize(const MatrixType& a)
215 {
216 eigen_assert(a.rows()==a.cols());
217 Index size = a.cols();
218 CholMatrixType tmp(size,size);
219 ConstCholMatrixPtr pmat;
220
221 if(m_P.size() == 0 && (int(UpLo) & int(Upper)) == Upper)
222 {
223 // If there is no ordering, try to directly use the input matrix without any copy
224 internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, tmp);
225 }
226 else
227 {
228 tmp.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
229 pmat = &tmp;
230 }
231
232 factorize_preordered<DoLDLT>(*pmat);
233 }
234
235 template<bool DoLDLT>
236 void factorize_preordered(const CholMatrixType& a);
237
238 void analyzePattern(const MatrixType& a, bool doLDLT)
239 {
240 eigen_assert(a.rows()==a.cols());
241 Index size = a.cols();
242 CholMatrixType tmp(size,size);
243 ConstCholMatrixPtr pmat;
244 ordering(a, pmat, tmp);
245 analyzePattern_preordered(*pmat,doLDLT);
246 }
247 void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
248
249 void ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap);
250
252 struct keep_diag {
253 inline bool operator() (const Index& row, const Index& col, const Scalar&) const
254 {
255 return row!=col;
256 }
257 };
258
259 mutable ComputationInfo m_info;
260 bool m_factorizationIsOk;
261 bool m_analysisIsOk;
262
263 CholMatrixType m_matrix;
264 VectorType m_diag; // the diagonal coefficients (LDLT mode)
265 VectorI m_parent; // elimination tree
266 VectorI m_nonZerosPerCol;
268 PermutationMatrix<Dynamic,Dynamic,StorageIndex> m_Pinv; // the inverse permutation
269
270 RealScalar m_shiftOffset;
271 RealScalar m_shiftScale;
272};
273
274template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLLT;
275template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialLDLT;
276template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::StorageIndex> > class SimplicialCholesky;
277
278namespace internal {
279
280template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
281{
282 typedef _MatrixType MatrixType;
283 typedef _Ordering OrderingType;
284 enum { UpLo = _UpLo };
285 typedef typename MatrixType::Scalar Scalar;
286 typedef typename MatrixType::StorageIndex StorageIndex;
287 typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
288 typedef TriangularView<const CholMatrixType, Eigen::Lower> MatrixL;
289 typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
290 static inline MatrixL getL(const CholMatrixType& m) { return MatrixL(m); }
291 static inline MatrixU getU(const CholMatrixType& m) { return MatrixU(m.adjoint()); }
292};
293
294template<typename _MatrixType,int _UpLo, typename _Ordering> struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
295{
296 typedef _MatrixType MatrixType;
297 typedef _Ordering OrderingType;
298 enum { UpLo = _UpLo };
299 typedef typename MatrixType::Scalar Scalar;
300 typedef typename MatrixType::StorageIndex StorageIndex;
301 typedef SparseMatrix<Scalar, ColMajor, StorageIndex> CholMatrixType;
302 typedef TriangularView<const CholMatrixType, Eigen::UnitLower> MatrixL;
303 typedef TriangularView<const typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
304 static inline MatrixL getL(const CholMatrixType& m) { return MatrixL(m); }
305 static inline MatrixU getU(const CholMatrixType& m) { return MatrixU(m.adjoint()); }
306};
307
308template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
309{
310 typedef _MatrixType MatrixType;
311 typedef _Ordering OrderingType;
312 enum { UpLo = _UpLo };
313};
314
315}
316
337template<typename _MatrixType, int _UpLo, typename _Ordering>
338 class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
339{
340public:
341 typedef _MatrixType MatrixType;
342 enum { UpLo = _UpLo };
344 typedef typename MatrixType::Scalar Scalar;
345 typedef typename MatrixType::RealScalar RealScalar;
346 typedef typename MatrixType::StorageIndex StorageIndex;
349 typedef internal::traits<SimplicialLLT> Traits;
350 typedef typename Traits::MatrixL MatrixL;
351 typedef typename Traits::MatrixU MatrixU;
352public:
356 explicit SimplicialLLT(const MatrixType& matrix)
357 : Base(matrix) {}
358
360 inline const MatrixL matrixL() const {
361 eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
362 return Traits::getL(Base::m_matrix);
363 }
364
366 inline const MatrixU matrixU() const {
367 eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
368 return Traits::getU(Base::m_matrix);
369 }
370
372 SimplicialLLT& compute(const MatrixType& matrix)
373 {
374 Base::template compute<false>(matrix);
375 return *this;
376 }
377
384 void analyzePattern(const MatrixType& a)
385 {
386 Base::analyzePattern(a, false);
387 }
388
395 void factorize(const MatrixType& a)
396 {
397 Base::template factorize<false>(a);
398 }
399
401 Scalar determinant() const
402 {
403 Scalar detL = Base::m_matrix.diagonal().prod();
404 return numext::abs2(detL);
405 }
406};
407
428template<typename _MatrixType, int _UpLo, typename _Ordering>
429 class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
430{
431public:
432 typedef _MatrixType MatrixType;
433 enum { UpLo = _UpLo };
435 typedef typename MatrixType::Scalar Scalar;
436 typedef typename MatrixType::RealScalar RealScalar;
437 typedef typename MatrixType::StorageIndex StorageIndex;
440 typedef internal::traits<SimplicialLDLT> Traits;
441 typedef typename Traits::MatrixL MatrixL;
442 typedef typename Traits::MatrixU MatrixU;
443public:
446
448 explicit SimplicialLDLT(const MatrixType& matrix)
449 : Base(matrix) {}
450
452 inline const VectorType vectorD() const {
453 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
454 return Base::m_diag;
455 }
457 inline const MatrixL matrixL() const {
458 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
459 return Traits::getL(Base::m_matrix);
460 }
461
463 inline const MatrixU matrixU() const {
464 eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
465 return Traits::getU(Base::m_matrix);
466 }
467
469 SimplicialLDLT& compute(const MatrixType& matrix)
470 {
471 Base::template compute<true>(matrix);
472 return *this;
473 }
474
481 void analyzePattern(const MatrixType& a)
482 {
483 Base::analyzePattern(a, true);
484 }
485
492 void factorize(const MatrixType& a)
493 {
494 Base::template factorize<true>(a);
495 }
496
498 Scalar determinant() const
499 {
500 return Base::m_diag.prod();
501 }
502};
503
510template<typename _MatrixType, int _UpLo, typename _Ordering>
511 class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
512{
513public:
514 typedef _MatrixType MatrixType;
515 enum { UpLo = _UpLo };
517 typedef typename MatrixType::Scalar Scalar;
518 typedef typename MatrixType::RealScalar RealScalar;
519 typedef typename MatrixType::StorageIndex StorageIndex;
522 typedef internal::traits<SimplicialCholesky> Traits;
523 typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
524 typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
525 public:
526 SimplicialCholesky() : Base(), m_LDLT(true) {}
527
528 explicit SimplicialCholesky(const MatrixType& matrix)
529 : Base(), m_LDLT(true)
530 {
531 compute(matrix);
532 }
533
534 SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
535 {
536 switch(mode)
537 {
538 case SimplicialCholeskyLLT:
539 m_LDLT = false;
540 break;
541 case SimplicialCholeskyLDLT:
542 m_LDLT = true;
543 break;
544 default:
545 break;
546 }
547
548 return *this;
549 }
550
551 inline const VectorType vectorD() const {
552 eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
553 return Base::m_diag;
554 }
555 inline const CholMatrixType rawMatrix() const {
556 eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
557 return Base::m_matrix;
558 }
559
561 SimplicialCholesky& compute(const MatrixType& matrix)
562 {
563 if(m_LDLT)
564 Base::template compute<true>(matrix);
565 else
566 Base::template compute<false>(matrix);
567 return *this;
568 }
569
576 void analyzePattern(const MatrixType& a)
577 {
578 Base::analyzePattern(a, m_LDLT);
579 }
580
587 void factorize(const MatrixType& a)
588 {
589 if(m_LDLT)
590 Base::template factorize<true>(a);
591 else
592 Base::template factorize<false>(a);
593 }
594
596 template<typename Rhs,typename Dest>
597 void _solve_impl(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
598 {
599 eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
600 eigen_assert(Base::m_matrix.rows()==b.rows());
601
602 if(Base::m_info!=Success)
603 return;
604
605 if(Base::m_P.size()>0)
606 dest = Base::m_P * b;
607 else
608 dest = b;
609
610 if(Base::m_matrix.nonZeros()>0) // otherwise L==I
611 {
612 if(m_LDLT)
613 LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
614 else
615 LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
616 }
617
618 if(Base::m_diag.size()>0)
619 dest = Base::m_diag.real().asDiagonal().inverse() * dest;
620
621 if (Base::m_matrix.nonZeros()>0) // otherwise I==I
622 {
623 if(m_LDLT)
624 LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
625 else
626 LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
627 }
628
629 if(Base::m_P.size()>0)
630 dest = Base::m_Pinv * dest;
631 }
632
634 template<typename Rhs,typename Dest>
635 void _solve_impl(const SparseMatrixBase<Rhs> &b, SparseMatrixBase<Dest> &dest) const
636 {
637 internal::solve_sparse_through_dense_panels(*this, b, dest);
638 }
639
640 Scalar determinant() const
641 {
642 if(m_LDLT)
643 {
644 return Base::m_diag.prod();
645 }
646 else
647 {
648 Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
649 return numext::abs2(detL);
650 }
651 }
652
653 protected:
654 bool m_LDLT;
655};
656
657template<typename Derived>
658void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, ConstCholMatrixPtr &pmat, CholMatrixType& ap)
659{
660 eigen_assert(a.rows()==a.cols());
661 const Index size = a.rows();
662 pmat = &ap;
663 // Note that ordering methods compute the inverse permutation
664 if(!internal::is_same<OrderingType,NaturalOrdering<Index> >::value)
665 {
666 {
667 CholMatrixType C;
668 C = a.template selfadjointView<UpLo>();
669
670 OrderingType ordering;
671 ordering(C,m_Pinv);
672 }
673
674 if(m_Pinv.size()>0) m_P = m_Pinv.inverse();
675 else m_P.resize(0);
676
677 ap.resize(size,size);
678 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
679 }
680 else
681 {
682 m_Pinv.resize(0);
683 m_P.resize(0);
684 if(int(UpLo)==int(Lower) || MatrixType::IsRowMajor)
685 {
686 // we have to transpose the lower part to to the upper one
687 ap.resize(size,size);
688 ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>();
689 }
690 else
691 internal::simplicial_cholesky_grab_input<CholMatrixType,MatrixType>::run(a, pmat, ap);
692 }
693}
694
695} // end namespace Eigen
696
697#endif // EIGEN_SIMPLICIAL_CHOLESKY_H
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
const DiagonalWrapper< const Derived > asDiagonal() const
Definition: DiagonalMatrix.h:325
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Index size() const
Definition: PermutationMatrix.h:97
Permutation matrix.
Definition: PermutationMatrix.h:298
A base class for direct sparse Cholesky factorizations.
Definition: SimplicialCholesky.h:56
SimplicialCholeskyBase()
Definition: SimplicialCholesky.h:82
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SimplicialCholesky.h:115
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationP() const
Definition: SimplicialCholesky.h:123
Derived & setShift(const RealScalar &offset, const RealScalar &scale=1)
Definition: SimplicialCholesky.h:140
void compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:202
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationPinv() const
Definition: SimplicialCholesky.h:128
Definition: SimplicialCholesky.h:512
SimplicialCholesky & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:561
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:576
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:587
A direct sparse LDLT Cholesky factorizations without square root.
Definition: SimplicialCholesky.h:430
SimplicialLDLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:448
SimplicialLDLT()
Definition: SimplicialCholesky.h:445
SimplicialLDLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:469
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:492
Scalar determinant() const
Definition: SimplicialCholesky.h:498
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:481
const VectorType vectorD() const
Definition: SimplicialCholesky.h:452
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:457
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:463
A direct sparse LLT Cholesky factorizations.
Definition: SimplicialCholesky.h:339
const MatrixU matrixU() const
Definition: SimplicialCholesky.h:366
SimplicialLLT(const MatrixType &matrix)
Definition: SimplicialCholesky.h:356
SimplicialLLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:372
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:395
Scalar determinant() const
Definition: SimplicialCholesky.h:401
SimplicialLLT()
Definition: SimplicialCholesky.h:354
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:384
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:360
A versatible sparse matrix representation.
Definition: SparseMatrix.h:98
Index nonZeros() const
Definition: SparseCompressedBase.h:56
Index rows() const
Definition: SparseMatrix.h:138
Index cols() const
Definition: SparseMatrix.h:140
A base class for sparse solvers.
Definition: SparseSolverBase.h:68
ComputationInfo
Definition: Constants.h:440
@ Upper
Definition: Constants.h:211
@ Success
Definition: Constants.h:442
Namespace containing all symbols from the Eigen library.
Definition: Core:141
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
Definition: SimplicialCholesky.h:252