Eigen  3.4.90 (git rev a4098ac676528a83cfb73d4d26ce1b42ec05f47c)
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
15namespace Eigen {
16
17enum SimplicialCholeskyMode {
18 SimplicialCholeskyLLT,
19 SimplicialCholeskyLDLT
20};
21
22namespace 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
56template<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
102 ~SimplicialCholeskyBase()
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
276template<typename MatrixType_, int UpLo_ = Lower, typename Ordering_ = AMDOrdering<typename MatrixType_::StorageIndex> > class SimplicialLLT;
277template<typename MatrixType_, int UpLo_ = Lower, typename Ordering_ = AMDOrdering<typename MatrixType_::StorageIndex> > class SimplicialLDLT;
278template<typename MatrixType_, int UpLo_ = Lower, typename Ordering_ = AMDOrdering<typename MatrixType_::StorageIndex> > class SimplicialCholesky;
279
280namespace internal {
281
282template<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
296template<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
310template<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
339template<typename MatrixType_, int UpLo_, typename Ordering_>
340 class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<MatrixType_,UpLo_,Ordering_> >
341{
342public:
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;
354public:
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
430template<typename MatrixType_, int UpLo_, typename Ordering_>
431 class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<MatrixType_,UpLo_,Ordering_> >
432{
433public:
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;
445public:
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
512template<typename MatrixType_, int UpLo_, typename Ordering_>
513 class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<MatrixType_,UpLo_,Ordering_> >
514{
515public:
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
659template<typename Derived>
660void 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:325
Index size() const
Definition: PermutationMatrix.h:99
A base class for direct sparse Cholesky factorizations.
Definition: SimplicialCholesky.h:58
SimplicialCholeskyBase()
Definition: SimplicialCholesky.h:84
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SimplicialCholesky.h:117
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationP() const
Definition: SimplicialCholesky.h:125
Derived & setShift(const RealScalar &offset, const RealScalar &scale=1)
Definition: SimplicialCholesky.h:142
void compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:204
const PermutationMatrix< Dynamic, Dynamic, StorageIndex > & permutationPinv() const
Definition: SimplicialCholesky.h:130
Definition: SimplicialCholesky.h:514
void factorize(const MatrixType &a)
Definition: SimplicialCholesky.h:589
void analyzePattern(const MatrixType &a)
Definition: SimplicialCholesky.h:578
SimplicialCholesky & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:563
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
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 & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:471
SimplicialLDLT()
Definition: SimplicialCholesky.h:447
A direct sparse LLT Cholesky factorizations.
Definition: SimplicialCholesky.h:341
const MatrixL matrixL() const
Definition: SimplicialCholesky.h:362
SimplicialLLT & compute(const MatrixType &matrix)
Definition: SimplicialCholesky.h:374
Scalar determinant() const
Definition: SimplicialCholesky.h:403
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: B01_Experimental.dox:1
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
Definition: SimplicialCholesky.h:254