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 | /** \ingroup SparseCholesky_Module
|
---|
21 | * \brief A direct sparse Cholesky factorizations
|
---|
22 | *
|
---|
23 | * These classes provide LL^T and LDL^T Cholesky factorizations of sparse matrices that are
|
---|
24 | * selfadjoint and positive definite. The factorization allows for solving A.X = B where
|
---|
25 | * X and B can be either dense or sparse.
|
---|
26 | *
|
---|
27 | * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
|
---|
28 | * such that the factorized matrix is P A P^-1.
|
---|
29 | *
|
---|
30 | * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
|
---|
31 | * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
|
---|
32 | * or Upper. Default is Lower.
|
---|
33 | *
|
---|
34 | */
|
---|
35 | template<typename Derived>
|
---|
36 | class SimplicialCholeskyBase : internal::noncopyable
|
---|
37 | {
|
---|
38 | public:
|
---|
39 | typedef typename internal::traits<Derived>::MatrixType MatrixType;
|
---|
40 | typedef typename internal::traits<Derived>::OrderingType OrderingType;
|
---|
41 | enum { UpLo = internal::traits<Derived>::UpLo };
|
---|
42 | typedef typename MatrixType::Scalar Scalar;
|
---|
43 | typedef typename MatrixType::RealScalar RealScalar;
|
---|
44 | typedef typename MatrixType::Index Index;
|
---|
45 | typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
|
---|
46 | typedef Matrix<Scalar,Dynamic,1> VectorType;
|
---|
47 |
|
---|
48 | public:
|
---|
49 |
|
---|
50 | /** Default constructor */
|
---|
51 | SimplicialCholeskyBase()
|
---|
52 | : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
|
---|
53 | {}
|
---|
54 |
|
---|
55 | SimplicialCholeskyBase(const MatrixType& matrix)
|
---|
56 | : m_info(Success), m_isInitialized(false), m_shiftOffset(0), m_shiftScale(1)
|
---|
57 | {
|
---|
58 | derived().compute(matrix);
|
---|
59 | }
|
---|
60 |
|
---|
61 | ~SimplicialCholeskyBase()
|
---|
62 | {
|
---|
63 | }
|
---|
64 |
|
---|
65 | Derived& derived() { return *static_cast<Derived*>(this); }
|
---|
66 | const Derived& derived() const { return *static_cast<const Derived*>(this); }
|
---|
67 |
|
---|
68 | inline Index cols() const { return m_matrix.cols(); }
|
---|
69 | inline Index rows() const { return m_matrix.rows(); }
|
---|
70 |
|
---|
71 | /** \brief Reports whether previous computation was successful.
|
---|
72 | *
|
---|
73 | * \returns \c Success if computation was succesful,
|
---|
74 | * \c NumericalIssue if the matrix.appears to be negative.
|
---|
75 | */
|
---|
76 | ComputationInfo info() const
|
---|
77 | {
|
---|
78 | eigen_assert(m_isInitialized && "Decomposition is not initialized.");
|
---|
79 | return m_info;
|
---|
80 | }
|
---|
81 |
|
---|
82 | /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
|
---|
83 | *
|
---|
84 | * \sa compute()
|
---|
85 | */
|
---|
86 | template<typename Rhs>
|
---|
87 | inline const internal::solve_retval<SimplicialCholeskyBase, Rhs>
|
---|
88 | solve(const MatrixBase<Rhs>& b) const
|
---|
89 | {
|
---|
90 | eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
|
---|
91 | eigen_assert(rows()==b.rows()
|
---|
92 | && "SimplicialCholeskyBase::solve(): invalid number of rows of the right hand side matrix b");
|
---|
93 | return internal::solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
|
---|
94 | }
|
---|
95 |
|
---|
96 | /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
|
---|
97 | *
|
---|
98 | * \sa compute()
|
---|
99 | */
|
---|
100 | template<typename Rhs>
|
---|
101 | inline const internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>
|
---|
102 | solve(const SparseMatrixBase<Rhs>& b) const
|
---|
103 | {
|
---|
104 | eigen_assert(m_isInitialized && "Simplicial LLT or LDLT is not initialized.");
|
---|
105 | eigen_assert(rows()==b.rows()
|
---|
106 | && "SimplicialCholesky::solve(): invalid number of rows of the right hand side matrix b");
|
---|
107 | return internal::sparse_solve_retval<SimplicialCholeskyBase, Rhs>(*this, b.derived());
|
---|
108 | }
|
---|
109 |
|
---|
110 | /** \returns the permutation P
|
---|
111 | * \sa permutationPinv() */
|
---|
112 | const PermutationMatrix<Dynamic,Dynamic,Index>& permutationP() const
|
---|
113 | { return m_P; }
|
---|
114 |
|
---|
115 | /** \returns the inverse P^-1 of the permutation P
|
---|
116 | * \sa permutationP() */
|
---|
117 | const PermutationMatrix<Dynamic,Dynamic,Index>& permutationPinv() const
|
---|
118 | { return m_Pinv; }
|
---|
119 |
|
---|
120 | /** Sets the shift parameters that will be used to adjust the diagonal coefficients during the numerical factorization.
|
---|
121 | *
|
---|
122 | * During the numerical factorization, the diagonal coefficients are transformed by the following linear model:\n
|
---|
123 | * \c d_ii = \a offset + \a scale * \c d_ii
|
---|
124 | *
|
---|
125 | * The default is the identity transformation with \a offset=0, and \a scale=1.
|
---|
126 | *
|
---|
127 | * \returns a reference to \c *this.
|
---|
128 | */
|
---|
129 | Derived& setShift(const RealScalar& offset, const RealScalar& scale = 1)
|
---|
130 | {
|
---|
131 | m_shiftOffset = offset;
|
---|
132 | m_shiftScale = scale;
|
---|
133 | return derived();
|
---|
134 | }
|
---|
135 |
|
---|
136 | #ifndef EIGEN_PARSED_BY_DOXYGEN
|
---|
137 | /** \internal */
|
---|
138 | template<typename Stream>
|
---|
139 | void dumpMemory(Stream& s)
|
---|
140 | {
|
---|
141 | int total = 0;
|
---|
142 | s << " L: " << ((total+=(m_matrix.cols()+1) * sizeof(int) + m_matrix.nonZeros()*(sizeof(int)+sizeof(Scalar))) >> 20) << "Mb" << "\n";
|
---|
143 | s << " diag: " << ((total+=m_diag.size() * sizeof(Scalar)) >> 20) << "Mb" << "\n";
|
---|
144 | s << " tree: " << ((total+=m_parent.size() * sizeof(int)) >> 20) << "Mb" << "\n";
|
---|
145 | s << " nonzeros: " << ((total+=m_nonZerosPerCol.size() * sizeof(int)) >> 20) << "Mb" << "\n";
|
---|
146 | s << " perm: " << ((total+=m_P.size() * sizeof(int)) >> 20) << "Mb" << "\n";
|
---|
147 | s << " perm^-1: " << ((total+=m_Pinv.size() * sizeof(int)) >> 20) << "Mb" << "\n";
|
---|
148 | s << " TOTAL: " << (total>> 20) << "Mb" << "\n";
|
---|
149 | }
|
---|
150 |
|
---|
151 | /** \internal */
|
---|
152 | template<typename Rhs,typename Dest>
|
---|
153 | void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
|
---|
154 | {
|
---|
155 | eigen_assert(m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
|
---|
156 | eigen_assert(m_matrix.rows()==b.rows());
|
---|
157 |
|
---|
158 | if(m_info!=Success)
|
---|
159 | return;
|
---|
160 |
|
---|
161 | if(m_P.size()>0)
|
---|
162 | dest = m_P * b;
|
---|
163 | else
|
---|
164 | dest = b;
|
---|
165 |
|
---|
166 | if(m_matrix.nonZeros()>0) // otherwise L==I
|
---|
167 | derived().matrixL().solveInPlace(dest);
|
---|
168 |
|
---|
169 | if(m_diag.size()>0)
|
---|
170 | dest = m_diag.asDiagonal().inverse() * dest;
|
---|
171 |
|
---|
172 | if (m_matrix.nonZeros()>0) // otherwise U==I
|
---|
173 | derived().matrixU().solveInPlace(dest);
|
---|
174 |
|
---|
175 | if(m_P.size()>0)
|
---|
176 | dest = m_Pinv * dest;
|
---|
177 | }
|
---|
178 |
|
---|
179 | #endif // EIGEN_PARSED_BY_DOXYGEN
|
---|
180 |
|
---|
181 | protected:
|
---|
182 |
|
---|
183 | /** Computes the sparse Cholesky decomposition of \a matrix */
|
---|
184 | template<bool DoLDLT>
|
---|
185 | void compute(const MatrixType& matrix)
|
---|
186 | {
|
---|
187 | eigen_assert(matrix.rows()==matrix.cols());
|
---|
188 | Index size = matrix.cols();
|
---|
189 | CholMatrixType ap(size,size);
|
---|
190 | ordering(matrix, ap);
|
---|
191 | analyzePattern_preordered(ap, DoLDLT);
|
---|
192 | factorize_preordered<DoLDLT>(ap);
|
---|
193 | }
|
---|
194 |
|
---|
195 | template<bool DoLDLT>
|
---|
196 | void factorize(const MatrixType& a)
|
---|
197 | {
|
---|
198 | eigen_assert(a.rows()==a.cols());
|
---|
199 | int size = a.cols();
|
---|
200 | CholMatrixType ap(size,size);
|
---|
201 | ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
|
---|
202 | factorize_preordered<DoLDLT>(ap);
|
---|
203 | }
|
---|
204 |
|
---|
205 | template<bool DoLDLT>
|
---|
206 | void factorize_preordered(const CholMatrixType& a);
|
---|
207 |
|
---|
208 | void analyzePattern(const MatrixType& a, bool doLDLT)
|
---|
209 | {
|
---|
210 | eigen_assert(a.rows()==a.cols());
|
---|
211 | int size = a.cols();
|
---|
212 | CholMatrixType ap(size,size);
|
---|
213 | ordering(a, ap);
|
---|
214 | analyzePattern_preordered(ap,doLDLT);
|
---|
215 | }
|
---|
216 | void analyzePattern_preordered(const CholMatrixType& a, bool doLDLT);
|
---|
217 |
|
---|
218 | void ordering(const MatrixType& a, CholMatrixType& ap);
|
---|
219 |
|
---|
220 | /** keeps off-diagonal entries; drops diagonal entries */
|
---|
221 | struct keep_diag {
|
---|
222 | inline bool operator() (const Index& row, const Index& col, const Scalar&) const
|
---|
223 | {
|
---|
224 | return row!=col;
|
---|
225 | }
|
---|
226 | };
|
---|
227 |
|
---|
228 | mutable ComputationInfo m_info;
|
---|
229 | bool m_isInitialized;
|
---|
230 | bool m_factorizationIsOk;
|
---|
231 | bool m_analysisIsOk;
|
---|
232 |
|
---|
233 | CholMatrixType m_matrix;
|
---|
234 | VectorType m_diag; // the diagonal coefficients (LDLT mode)
|
---|
235 | VectorXi m_parent; // elimination tree
|
---|
236 | VectorXi m_nonZerosPerCol;
|
---|
237 | PermutationMatrix<Dynamic,Dynamic,Index> m_P; // the permutation
|
---|
238 | PermutationMatrix<Dynamic,Dynamic,Index> m_Pinv; // the inverse permutation
|
---|
239 |
|
---|
240 | RealScalar m_shiftOffset;
|
---|
241 | RealScalar m_shiftScale;
|
---|
242 | };
|
---|
243 |
|
---|
244 | template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLLT;
|
---|
245 | template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialLDLT;
|
---|
246 | template<typename _MatrixType, int _UpLo = Lower, typename _Ordering = AMDOrdering<typename _MatrixType::Index> > class SimplicialCholesky;
|
---|
247 |
|
---|
248 | namespace internal {
|
---|
249 |
|
---|
250 | template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
|
---|
251 | {
|
---|
252 | typedef _MatrixType MatrixType;
|
---|
253 | typedef _Ordering OrderingType;
|
---|
254 | enum { UpLo = _UpLo };
|
---|
255 | typedef typename MatrixType::Scalar Scalar;
|
---|
256 | typedef typename MatrixType::Index Index;
|
---|
257 | typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
|
---|
258 | typedef SparseTriangularView<CholMatrixType, Eigen::Lower> MatrixL;
|
---|
259 | typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::Upper> MatrixU;
|
---|
260 | static inline MatrixL getL(const MatrixType& m) { return m; }
|
---|
261 | static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
|
---|
262 | };
|
---|
263 |
|
---|
264 | template<typename _MatrixType,int _UpLo, typename _Ordering> struct traits<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
|
---|
265 | {
|
---|
266 | typedef _MatrixType MatrixType;
|
---|
267 | typedef _Ordering OrderingType;
|
---|
268 | enum { UpLo = _UpLo };
|
---|
269 | typedef typename MatrixType::Scalar Scalar;
|
---|
270 | typedef typename MatrixType::Index Index;
|
---|
271 | typedef SparseMatrix<Scalar, ColMajor, Index> CholMatrixType;
|
---|
272 | typedef SparseTriangularView<CholMatrixType, Eigen::UnitLower> MatrixL;
|
---|
273 | typedef SparseTriangularView<typename CholMatrixType::AdjointReturnType, Eigen::UnitUpper> MatrixU;
|
---|
274 | static inline MatrixL getL(const MatrixType& m) { return m; }
|
---|
275 | static inline MatrixU getU(const MatrixType& m) { return m.adjoint(); }
|
---|
276 | };
|
---|
277 |
|
---|
278 | template<typename _MatrixType, int _UpLo, typename _Ordering> struct traits<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
|
---|
279 | {
|
---|
280 | typedef _MatrixType MatrixType;
|
---|
281 | typedef _Ordering OrderingType;
|
---|
282 | enum { UpLo = _UpLo };
|
---|
283 | };
|
---|
284 |
|
---|
285 | }
|
---|
286 |
|
---|
287 | /** \ingroup SparseCholesky_Module
|
---|
288 | * \class SimplicialLLT
|
---|
289 | * \brief A direct sparse LLT Cholesky factorizations
|
---|
290 | *
|
---|
291 | * This class provides a LL^T Cholesky factorizations of sparse matrices that are
|
---|
292 | * selfadjoint and positive definite. The factorization allows for solving A.X = B where
|
---|
293 | * X and B can be either dense or sparse.
|
---|
294 | *
|
---|
295 | * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
|
---|
296 | * such that the factorized matrix is P A P^-1.
|
---|
297 | *
|
---|
298 | * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
|
---|
299 | * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
|
---|
300 | * or Upper. Default is Lower.
|
---|
301 | * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<>
|
---|
302 | *
|
---|
303 | * \sa class SimplicialLDLT, class AMDOrdering, class NaturalOrdering
|
---|
304 | */
|
---|
305 | template<typename _MatrixType, int _UpLo, typename _Ordering>
|
---|
306 | class SimplicialLLT : public SimplicialCholeskyBase<SimplicialLLT<_MatrixType,_UpLo,_Ordering> >
|
---|
307 | {
|
---|
308 | public:
|
---|
309 | typedef _MatrixType MatrixType;
|
---|
310 | enum { UpLo = _UpLo };
|
---|
311 | typedef SimplicialCholeskyBase<SimplicialLLT> Base;
|
---|
312 | typedef typename MatrixType::Scalar Scalar;
|
---|
313 | typedef typename MatrixType::RealScalar RealScalar;
|
---|
314 | typedef typename MatrixType::Index Index;
|
---|
315 | typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
|
---|
316 | typedef Matrix<Scalar,Dynamic,1> VectorType;
|
---|
317 | typedef internal::traits<SimplicialLLT> Traits;
|
---|
318 | typedef typename Traits::MatrixL MatrixL;
|
---|
319 | typedef typename Traits::MatrixU MatrixU;
|
---|
320 | public:
|
---|
321 | /** Default constructor */
|
---|
322 | SimplicialLLT() : Base() {}
|
---|
323 | /** Constructs and performs the LLT factorization of \a matrix */
|
---|
324 | SimplicialLLT(const MatrixType& matrix)
|
---|
325 | : Base(matrix) {}
|
---|
326 |
|
---|
327 | /** \returns an expression of the factor L */
|
---|
328 | inline const MatrixL matrixL() const {
|
---|
329 | eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
|
---|
330 | return Traits::getL(Base::m_matrix);
|
---|
331 | }
|
---|
332 |
|
---|
333 | /** \returns an expression of the factor U (= L^*) */
|
---|
334 | inline const MatrixU matrixU() const {
|
---|
335 | eigen_assert(Base::m_factorizationIsOk && "Simplicial LLT not factorized");
|
---|
336 | return Traits::getU(Base::m_matrix);
|
---|
337 | }
|
---|
338 |
|
---|
339 | /** Computes the sparse Cholesky decomposition of \a matrix */
|
---|
340 | SimplicialLLT& compute(const MatrixType& matrix)
|
---|
341 | {
|
---|
342 | Base::template compute<false>(matrix);
|
---|
343 | return *this;
|
---|
344 | }
|
---|
345 |
|
---|
346 | /** Performs a symbolic decomposition on the sparcity of \a matrix.
|
---|
347 | *
|
---|
348 | * This function is particularly useful when solving for several problems having the same structure.
|
---|
349 | *
|
---|
350 | * \sa factorize()
|
---|
351 | */
|
---|
352 | void analyzePattern(const MatrixType& a)
|
---|
353 | {
|
---|
354 | Base::analyzePattern(a, false);
|
---|
355 | }
|
---|
356 |
|
---|
357 | /** Performs a numeric decomposition of \a matrix
|
---|
358 | *
|
---|
359 | * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
|
---|
360 | *
|
---|
361 | * \sa analyzePattern()
|
---|
362 | */
|
---|
363 | void factorize(const MatrixType& a)
|
---|
364 | {
|
---|
365 | Base::template factorize<false>(a);
|
---|
366 | }
|
---|
367 |
|
---|
368 | /** \returns the determinant of the underlying matrix from the current factorization */
|
---|
369 | Scalar determinant() const
|
---|
370 | {
|
---|
371 | Scalar detL = Base::m_matrix.diagonal().prod();
|
---|
372 | return numext::abs2(detL);
|
---|
373 | }
|
---|
374 | };
|
---|
375 |
|
---|
376 | /** \ingroup SparseCholesky_Module
|
---|
377 | * \class SimplicialLDLT
|
---|
378 | * \brief A direct sparse LDLT Cholesky factorizations without square root.
|
---|
379 | *
|
---|
380 | * This class provides a LDL^T Cholesky factorizations without square root of sparse matrices that are
|
---|
381 | * selfadjoint and positive definite. The factorization allows for solving A.X = B where
|
---|
382 | * X and B can be either dense or sparse.
|
---|
383 | *
|
---|
384 | * In order to reduce the fill-in, a symmetric permutation P is applied prior to the factorization
|
---|
385 | * such that the factorized matrix is P A P^-1.
|
---|
386 | *
|
---|
387 | * \tparam _MatrixType the type of the sparse matrix A, it must be a SparseMatrix<>
|
---|
388 | * \tparam _UpLo the triangular part that will be used for the computations. It can be Lower
|
---|
389 | * or Upper. Default is Lower.
|
---|
390 | * \tparam _Ordering The ordering method to use, either AMDOrdering<> or NaturalOrdering<>. Default is AMDOrdering<>
|
---|
391 | *
|
---|
392 | * \sa class SimplicialLLT, class AMDOrdering, class NaturalOrdering
|
---|
393 | */
|
---|
394 | template<typename _MatrixType, int _UpLo, typename _Ordering>
|
---|
395 | class SimplicialLDLT : public SimplicialCholeskyBase<SimplicialLDLT<_MatrixType,_UpLo,_Ordering> >
|
---|
396 | {
|
---|
397 | public:
|
---|
398 | typedef _MatrixType MatrixType;
|
---|
399 | enum { UpLo = _UpLo };
|
---|
400 | typedef SimplicialCholeskyBase<SimplicialLDLT> Base;
|
---|
401 | typedef typename MatrixType::Scalar Scalar;
|
---|
402 | typedef typename MatrixType::RealScalar RealScalar;
|
---|
403 | typedef typename MatrixType::Index Index;
|
---|
404 | typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
|
---|
405 | typedef Matrix<Scalar,Dynamic,1> VectorType;
|
---|
406 | typedef internal::traits<SimplicialLDLT> Traits;
|
---|
407 | typedef typename Traits::MatrixL MatrixL;
|
---|
408 | typedef typename Traits::MatrixU MatrixU;
|
---|
409 | public:
|
---|
410 | /** Default constructor */
|
---|
411 | SimplicialLDLT() : Base() {}
|
---|
412 |
|
---|
413 | /** Constructs and performs the LLT factorization of \a matrix */
|
---|
414 | SimplicialLDLT(const MatrixType& matrix)
|
---|
415 | : Base(matrix) {}
|
---|
416 |
|
---|
417 | /** \returns a vector expression of the diagonal D */
|
---|
418 | inline const VectorType vectorD() const {
|
---|
419 | eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
|
---|
420 | return Base::m_diag;
|
---|
421 | }
|
---|
422 | /** \returns an expression of the factor L */
|
---|
423 | inline const MatrixL matrixL() const {
|
---|
424 | eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
|
---|
425 | return Traits::getL(Base::m_matrix);
|
---|
426 | }
|
---|
427 |
|
---|
428 | /** \returns an expression of the factor U (= L^*) */
|
---|
429 | inline const MatrixU matrixU() const {
|
---|
430 | eigen_assert(Base::m_factorizationIsOk && "Simplicial LDLT not factorized");
|
---|
431 | return Traits::getU(Base::m_matrix);
|
---|
432 | }
|
---|
433 |
|
---|
434 | /** Computes the sparse Cholesky decomposition of \a matrix */
|
---|
435 | SimplicialLDLT& compute(const MatrixType& matrix)
|
---|
436 | {
|
---|
437 | Base::template compute<true>(matrix);
|
---|
438 | return *this;
|
---|
439 | }
|
---|
440 |
|
---|
441 | /** Performs a symbolic decomposition on the sparcity of \a matrix.
|
---|
442 | *
|
---|
443 | * This function is particularly useful when solving for several problems having the same structure.
|
---|
444 | *
|
---|
445 | * \sa factorize()
|
---|
446 | */
|
---|
447 | void analyzePattern(const MatrixType& a)
|
---|
448 | {
|
---|
449 | Base::analyzePattern(a, true);
|
---|
450 | }
|
---|
451 |
|
---|
452 | /** Performs a numeric decomposition of \a matrix
|
---|
453 | *
|
---|
454 | * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
|
---|
455 | *
|
---|
456 | * \sa analyzePattern()
|
---|
457 | */
|
---|
458 | void factorize(const MatrixType& a)
|
---|
459 | {
|
---|
460 | Base::template factorize<true>(a);
|
---|
461 | }
|
---|
462 |
|
---|
463 | /** \returns the determinant of the underlying matrix from the current factorization */
|
---|
464 | Scalar determinant() const
|
---|
465 | {
|
---|
466 | return Base::m_diag.prod();
|
---|
467 | }
|
---|
468 | };
|
---|
469 |
|
---|
470 | /** \deprecated use SimplicialLDLT or class SimplicialLLT
|
---|
471 | * \ingroup SparseCholesky_Module
|
---|
472 | * \class SimplicialCholesky
|
---|
473 | *
|
---|
474 | * \sa class SimplicialLDLT, class SimplicialLLT
|
---|
475 | */
|
---|
476 | template<typename _MatrixType, int _UpLo, typename _Ordering>
|
---|
477 | class SimplicialCholesky : public SimplicialCholeskyBase<SimplicialCholesky<_MatrixType,_UpLo,_Ordering> >
|
---|
478 | {
|
---|
479 | public:
|
---|
480 | typedef _MatrixType MatrixType;
|
---|
481 | enum { UpLo = _UpLo };
|
---|
482 | typedef SimplicialCholeskyBase<SimplicialCholesky> Base;
|
---|
483 | typedef typename MatrixType::Scalar Scalar;
|
---|
484 | typedef typename MatrixType::RealScalar RealScalar;
|
---|
485 | typedef typename MatrixType::Index Index;
|
---|
486 | typedef SparseMatrix<Scalar,ColMajor,Index> CholMatrixType;
|
---|
487 | typedef Matrix<Scalar,Dynamic,1> VectorType;
|
---|
488 | typedef internal::traits<SimplicialCholesky> Traits;
|
---|
489 | typedef internal::traits<SimplicialLDLT<MatrixType,UpLo> > LDLTTraits;
|
---|
490 | typedef internal::traits<SimplicialLLT<MatrixType,UpLo> > LLTTraits;
|
---|
491 | public:
|
---|
492 | SimplicialCholesky() : Base(), m_LDLT(true) {}
|
---|
493 |
|
---|
494 | SimplicialCholesky(const MatrixType& matrix)
|
---|
495 | : Base(), m_LDLT(true)
|
---|
496 | {
|
---|
497 | compute(matrix);
|
---|
498 | }
|
---|
499 |
|
---|
500 | SimplicialCholesky& setMode(SimplicialCholeskyMode mode)
|
---|
501 | {
|
---|
502 | switch(mode)
|
---|
503 | {
|
---|
504 | case SimplicialCholeskyLLT:
|
---|
505 | m_LDLT = false;
|
---|
506 | break;
|
---|
507 | case SimplicialCholeskyLDLT:
|
---|
508 | m_LDLT = true;
|
---|
509 | break;
|
---|
510 | default:
|
---|
511 | break;
|
---|
512 | }
|
---|
513 |
|
---|
514 | return *this;
|
---|
515 | }
|
---|
516 |
|
---|
517 | inline const VectorType vectorD() const {
|
---|
518 | eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
|
---|
519 | return Base::m_diag;
|
---|
520 | }
|
---|
521 | inline const CholMatrixType rawMatrix() const {
|
---|
522 | eigen_assert(Base::m_factorizationIsOk && "Simplicial Cholesky not factorized");
|
---|
523 | return Base::m_matrix;
|
---|
524 | }
|
---|
525 |
|
---|
526 | /** Computes the sparse Cholesky decomposition of \a matrix */
|
---|
527 | SimplicialCholesky& compute(const MatrixType& matrix)
|
---|
528 | {
|
---|
529 | if(m_LDLT)
|
---|
530 | Base::template compute<true>(matrix);
|
---|
531 | else
|
---|
532 | Base::template compute<false>(matrix);
|
---|
533 | return *this;
|
---|
534 | }
|
---|
535 |
|
---|
536 | /** Performs a symbolic decomposition on the sparcity of \a matrix.
|
---|
537 | *
|
---|
538 | * This function is particularly useful when solving for several problems having the same structure.
|
---|
539 | *
|
---|
540 | * \sa factorize()
|
---|
541 | */
|
---|
542 | void analyzePattern(const MatrixType& a)
|
---|
543 | {
|
---|
544 | Base::analyzePattern(a, m_LDLT);
|
---|
545 | }
|
---|
546 |
|
---|
547 | /** Performs a numeric decomposition of \a matrix
|
---|
548 | *
|
---|
549 | * The given matrix must has the same sparcity than the matrix on which the symbolic decomposition has been performed.
|
---|
550 | *
|
---|
551 | * \sa analyzePattern()
|
---|
552 | */
|
---|
553 | void factorize(const MatrixType& a)
|
---|
554 | {
|
---|
555 | if(m_LDLT)
|
---|
556 | Base::template factorize<true>(a);
|
---|
557 | else
|
---|
558 | Base::template factorize<false>(a);
|
---|
559 | }
|
---|
560 |
|
---|
561 | /** \internal */
|
---|
562 | template<typename Rhs,typename Dest>
|
---|
563 | void _solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
|
---|
564 | {
|
---|
565 | eigen_assert(Base::m_factorizationIsOk && "The decomposition is not in a valid state for solving, you must first call either compute() or symbolic()/numeric()");
|
---|
566 | eigen_assert(Base::m_matrix.rows()==b.rows());
|
---|
567 |
|
---|
568 | if(Base::m_info!=Success)
|
---|
569 | return;
|
---|
570 |
|
---|
571 | if(Base::m_P.size()>0)
|
---|
572 | dest = Base::m_P * b;
|
---|
573 | else
|
---|
574 | dest = b;
|
---|
575 |
|
---|
576 | if(Base::m_matrix.nonZeros()>0) // otherwise L==I
|
---|
577 | {
|
---|
578 | if(m_LDLT)
|
---|
579 | LDLTTraits::getL(Base::m_matrix).solveInPlace(dest);
|
---|
580 | else
|
---|
581 | LLTTraits::getL(Base::m_matrix).solveInPlace(dest);
|
---|
582 | }
|
---|
583 |
|
---|
584 | if(Base::m_diag.size()>0)
|
---|
585 | dest = Base::m_diag.asDiagonal().inverse() * dest;
|
---|
586 |
|
---|
587 | if (Base::m_matrix.nonZeros()>0) // otherwise I==I
|
---|
588 | {
|
---|
589 | if(m_LDLT)
|
---|
590 | LDLTTraits::getU(Base::m_matrix).solveInPlace(dest);
|
---|
591 | else
|
---|
592 | LLTTraits::getU(Base::m_matrix).solveInPlace(dest);
|
---|
593 | }
|
---|
594 |
|
---|
595 | if(Base::m_P.size()>0)
|
---|
596 | dest = Base::m_Pinv * dest;
|
---|
597 | }
|
---|
598 |
|
---|
599 | Scalar determinant() const
|
---|
600 | {
|
---|
601 | if(m_LDLT)
|
---|
602 | {
|
---|
603 | return Base::m_diag.prod();
|
---|
604 | }
|
---|
605 | else
|
---|
606 | {
|
---|
607 | Scalar detL = Diagonal<const CholMatrixType>(Base::m_matrix).prod();
|
---|
608 | return numext::abs2(detL);
|
---|
609 | }
|
---|
610 | }
|
---|
611 |
|
---|
612 | protected:
|
---|
613 | bool m_LDLT;
|
---|
614 | };
|
---|
615 |
|
---|
616 | template<typename Derived>
|
---|
617 | void SimplicialCholeskyBase<Derived>::ordering(const MatrixType& a, CholMatrixType& ap)
|
---|
618 | {
|
---|
619 | eigen_assert(a.rows()==a.cols());
|
---|
620 | const Index size = a.rows();
|
---|
621 | // Note that amd compute the inverse permutation
|
---|
622 | {
|
---|
623 | CholMatrixType C;
|
---|
624 | C = a.template selfadjointView<UpLo>();
|
---|
625 |
|
---|
626 | OrderingType ordering;
|
---|
627 | ordering(C,m_Pinv);
|
---|
628 | }
|
---|
629 |
|
---|
630 | if(m_Pinv.size()>0)
|
---|
631 | m_P = m_Pinv.inverse();
|
---|
632 | else
|
---|
633 | m_P.resize(0);
|
---|
634 |
|
---|
635 | ap.resize(size,size);
|
---|
636 | ap.template selfadjointView<Upper>() = a.template selfadjointView<UpLo>().twistedBy(m_P);
|
---|
637 | }
|
---|
638 |
|
---|
639 | namespace internal {
|
---|
640 |
|
---|
641 | template<typename Derived, typename Rhs>
|
---|
642 | struct solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
|
---|
643 | : solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
|
---|
644 | {
|
---|
645 | typedef SimplicialCholeskyBase<Derived> Dec;
|
---|
646 | EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
|
---|
647 |
|
---|
648 | template<typename Dest> void evalTo(Dest& dst) const
|
---|
649 | {
|
---|
650 | dec().derived()._solve(rhs(),dst);
|
---|
651 | }
|
---|
652 | };
|
---|
653 |
|
---|
654 | template<typename Derived, typename Rhs>
|
---|
655 | struct sparse_solve_retval<SimplicialCholeskyBase<Derived>, Rhs>
|
---|
656 | : sparse_solve_retval_base<SimplicialCholeskyBase<Derived>, Rhs>
|
---|
657 | {
|
---|
658 | typedef SimplicialCholeskyBase<Derived> Dec;
|
---|
659 | EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
|
---|
660 |
|
---|
661 | template<typename Dest> void evalTo(Dest& dst) const
|
---|
662 | {
|
---|
663 | this->defaultEvalTo(dst);
|
---|
664 | }
|
---|
665 | };
|
---|
666 |
|
---|
667 | } // end namespace internal
|
---|
668 |
|
---|
669 | } // end namespace Eigen
|
---|
670 |
|
---|
671 | #endif // EIGEN_SIMPLICIAL_CHOLESKY_H
|
---|