1 | // This file is part of Eigen, a lightweight C++ template library
|
---|
2 | // for linear algebra.
|
---|
3 | //
|
---|
4 | // Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
|
---|
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_LU_H
|
---|
11 | #define EIGEN_LU_H
|
---|
12 |
|
---|
13 | namespace Eigen {
|
---|
14 |
|
---|
15 | /** \ingroup LU_Module
|
---|
16 | *
|
---|
17 | * \class FullPivLU
|
---|
18 | *
|
---|
19 | * \brief LU decomposition of a matrix with complete pivoting, and related features
|
---|
20 | *
|
---|
21 | * \param MatrixType the type of the matrix of which we are computing the LU decomposition
|
---|
22 | *
|
---|
23 | * This class represents a LU decomposition of any matrix, with complete pivoting: the matrix A is
|
---|
24 | * decomposed as \f$ A = P^{-1} L U Q^{-1} \f$ where L is unit-lower-triangular, U is
|
---|
25 | * upper-triangular, and P and Q are permutation matrices. This is a rank-revealing LU
|
---|
26 | * decomposition. The eigenvalues (diagonal coefficients) of U are sorted in such a way that any
|
---|
27 | * zeros are at the end.
|
---|
28 | *
|
---|
29 | * This decomposition provides the generic approach to solving systems of linear equations, computing
|
---|
30 | * the rank, invertibility, inverse, kernel, and determinant.
|
---|
31 | *
|
---|
32 | * This LU decomposition is very stable and well tested with large matrices. However there are use cases where the SVD
|
---|
33 | * decomposition is inherently more stable and/or flexible. For example, when computing the kernel of a matrix,
|
---|
34 | * working with the SVD allows to select the smallest singular values of the matrix, something that
|
---|
35 | * the LU decomposition doesn't see.
|
---|
36 | *
|
---|
37 | * The data of the LU decomposition can be directly accessed through the methods matrixLU(),
|
---|
38 | * permutationP(), permutationQ().
|
---|
39 | *
|
---|
40 | * As an exemple, here is how the original matrix can be retrieved:
|
---|
41 | * \include class_FullPivLU.cpp
|
---|
42 | * Output: \verbinclude class_FullPivLU.out
|
---|
43 | *
|
---|
44 | * \sa MatrixBase::fullPivLu(), MatrixBase::determinant(), MatrixBase::inverse()
|
---|
45 | */
|
---|
46 | template<typename _MatrixType> class FullPivLU
|
---|
47 | {
|
---|
48 | public:
|
---|
49 | typedef _MatrixType MatrixType;
|
---|
50 | enum {
|
---|
51 | RowsAtCompileTime = MatrixType::RowsAtCompileTime,
|
---|
52 | ColsAtCompileTime = MatrixType::ColsAtCompileTime,
|
---|
53 | Options = MatrixType::Options,
|
---|
54 | MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
|
---|
55 | MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
|
---|
56 | };
|
---|
57 | typedef typename MatrixType::Scalar Scalar;
|
---|
58 | typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
|
---|
59 | typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
|
---|
60 | typedef typename MatrixType::Index Index;
|
---|
61 | typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
|
---|
62 | typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
|
---|
63 | typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
|
---|
64 | typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType;
|
---|
65 |
|
---|
66 | /**
|
---|
67 | * \brief Default Constructor.
|
---|
68 | *
|
---|
69 | * The default constructor is useful in cases in which the user intends to
|
---|
70 | * perform decompositions via LU::compute(const MatrixType&).
|
---|
71 | */
|
---|
72 | FullPivLU();
|
---|
73 |
|
---|
74 | /** \brief Default Constructor with memory preallocation
|
---|
75 | *
|
---|
76 | * Like the default constructor but with preallocation of the internal data
|
---|
77 | * according to the specified problem \a size.
|
---|
78 | * \sa FullPivLU()
|
---|
79 | */
|
---|
80 | FullPivLU(Index rows, Index cols);
|
---|
81 |
|
---|
82 | /** Constructor.
|
---|
83 | *
|
---|
84 | * \param matrix the matrix of which to compute the LU decomposition.
|
---|
85 | * It is required to be nonzero.
|
---|
86 | */
|
---|
87 | FullPivLU(const MatrixType& matrix);
|
---|
88 |
|
---|
89 | /** Computes the LU decomposition of the given matrix.
|
---|
90 | *
|
---|
91 | * \param matrix the matrix of which to compute the LU decomposition.
|
---|
92 | * It is required to be nonzero.
|
---|
93 | *
|
---|
94 | * \returns a reference to *this
|
---|
95 | */
|
---|
96 | FullPivLU& compute(const MatrixType& matrix);
|
---|
97 |
|
---|
98 | /** \returns the LU decomposition matrix: the upper-triangular part is U, the
|
---|
99 | * unit-lower-triangular part is L (at least for square matrices; in the non-square
|
---|
100 | * case, special care is needed, see the documentation of class FullPivLU).
|
---|
101 | *
|
---|
102 | * \sa matrixL(), matrixU()
|
---|
103 | */
|
---|
104 | inline const MatrixType& matrixLU() const
|
---|
105 | {
|
---|
106 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
107 | return m_lu;
|
---|
108 | }
|
---|
109 |
|
---|
110 | /** \returns the number of nonzero pivots in the LU decomposition.
|
---|
111 | * Here nonzero is meant in the exact sense, not in a fuzzy sense.
|
---|
112 | * So that notion isn't really intrinsically interesting, but it is
|
---|
113 | * still useful when implementing algorithms.
|
---|
114 | *
|
---|
115 | * \sa rank()
|
---|
116 | */
|
---|
117 | inline Index nonzeroPivots() const
|
---|
118 | {
|
---|
119 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
120 | return m_nonzero_pivots;
|
---|
121 | }
|
---|
122 |
|
---|
123 | /** \returns the absolute value of the biggest pivot, i.e. the biggest
|
---|
124 | * diagonal coefficient of U.
|
---|
125 | */
|
---|
126 | RealScalar maxPivot() const { return m_maxpivot; }
|
---|
127 |
|
---|
128 | /** \returns the permutation matrix P
|
---|
129 | *
|
---|
130 | * \sa permutationQ()
|
---|
131 | */
|
---|
132 | inline const PermutationPType& permutationP() const
|
---|
133 | {
|
---|
134 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
135 | return m_p;
|
---|
136 | }
|
---|
137 |
|
---|
138 | /** \returns the permutation matrix Q
|
---|
139 | *
|
---|
140 | * \sa permutationP()
|
---|
141 | */
|
---|
142 | inline const PermutationQType& permutationQ() const
|
---|
143 | {
|
---|
144 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
145 | return m_q;
|
---|
146 | }
|
---|
147 |
|
---|
148 | /** \returns the kernel of the matrix, also called its null-space. The columns of the returned matrix
|
---|
149 | * will form a basis of the kernel.
|
---|
150 | *
|
---|
151 | * \note If the kernel has dimension zero, then the returned matrix is a column-vector filled with zeros.
|
---|
152 | *
|
---|
153 | * \note This method has to determine which pivots should be considered nonzero.
|
---|
154 | * For that, it uses the threshold value that you can control by calling
|
---|
155 | * setThreshold(const RealScalar&).
|
---|
156 | *
|
---|
157 | * Example: \include FullPivLU_kernel.cpp
|
---|
158 | * Output: \verbinclude FullPivLU_kernel.out
|
---|
159 | *
|
---|
160 | * \sa image()
|
---|
161 | */
|
---|
162 | inline const internal::kernel_retval<FullPivLU> kernel() const
|
---|
163 | {
|
---|
164 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
165 | return internal::kernel_retval<FullPivLU>(*this);
|
---|
166 | }
|
---|
167 |
|
---|
168 | /** \returns the image of the matrix, also called its column-space. The columns of the returned matrix
|
---|
169 | * will form a basis of the kernel.
|
---|
170 | *
|
---|
171 | * \param originalMatrix the original matrix, of which *this is the LU decomposition.
|
---|
172 | * The reason why it is needed to pass it here, is that this allows
|
---|
173 | * a large optimization, as otherwise this method would need to reconstruct it
|
---|
174 | * from the LU decomposition.
|
---|
175 | *
|
---|
176 | * \note If the image has dimension zero, then the returned matrix is a column-vector filled with zeros.
|
---|
177 | *
|
---|
178 | * \note This method has to determine which pivots should be considered nonzero.
|
---|
179 | * For that, it uses the threshold value that you can control by calling
|
---|
180 | * setThreshold(const RealScalar&).
|
---|
181 | *
|
---|
182 | * Example: \include FullPivLU_image.cpp
|
---|
183 | * Output: \verbinclude FullPivLU_image.out
|
---|
184 | *
|
---|
185 | * \sa kernel()
|
---|
186 | */
|
---|
187 | inline const internal::image_retval<FullPivLU>
|
---|
188 | image(const MatrixType& originalMatrix) const
|
---|
189 | {
|
---|
190 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
191 | return internal::image_retval<FullPivLU>(*this, originalMatrix);
|
---|
192 | }
|
---|
193 |
|
---|
194 | /** \return a solution x to the equation Ax=b, where A is the matrix of which
|
---|
195 | * *this is the LU decomposition.
|
---|
196 | *
|
---|
197 | * \param b the right-hand-side of the equation to solve. Can be a vector or a matrix,
|
---|
198 | * the only requirement in order for the equation to make sense is that
|
---|
199 | * b.rows()==A.rows(), where A is the matrix of which *this is the LU decomposition.
|
---|
200 | *
|
---|
201 | * \returns a solution.
|
---|
202 | *
|
---|
203 | * \note_about_checking_solutions
|
---|
204 | *
|
---|
205 | * \note_about_arbitrary_choice_of_solution
|
---|
206 | * \note_about_using_kernel_to_study_multiple_solutions
|
---|
207 | *
|
---|
208 | * Example: \include FullPivLU_solve.cpp
|
---|
209 | * Output: \verbinclude FullPivLU_solve.out
|
---|
210 | *
|
---|
211 | * \sa TriangularView::solve(), kernel(), inverse()
|
---|
212 | */
|
---|
213 | template<typename Rhs>
|
---|
214 | inline const internal::solve_retval<FullPivLU, Rhs>
|
---|
215 | solve(const MatrixBase<Rhs>& b) const
|
---|
216 | {
|
---|
217 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
218 | return internal::solve_retval<FullPivLU, Rhs>(*this, b.derived());
|
---|
219 | }
|
---|
220 |
|
---|
221 | /** \returns the determinant of the matrix of which
|
---|
222 | * *this is the LU decomposition. It has only linear complexity
|
---|
223 | * (that is, O(n) where n is the dimension of the square matrix)
|
---|
224 | * as the LU decomposition has already been computed.
|
---|
225 | *
|
---|
226 | * \note This is only for square matrices.
|
---|
227 | *
|
---|
228 | * \note For fixed-size matrices of size up to 4, MatrixBase::determinant() offers
|
---|
229 | * optimized paths.
|
---|
230 | *
|
---|
231 | * \warning a determinant can be very big or small, so for matrices
|
---|
232 | * of large enough dimension, there is a risk of overflow/underflow.
|
---|
233 | *
|
---|
234 | * \sa MatrixBase::determinant()
|
---|
235 | */
|
---|
236 | typename internal::traits<MatrixType>::Scalar determinant() const;
|
---|
237 |
|
---|
238 | /** Allows to prescribe a threshold to be used by certain methods, such as rank(),
|
---|
239 | * who need to determine when pivots are to be considered nonzero. This is not used for the
|
---|
240 | * LU decomposition itself.
|
---|
241 | *
|
---|
242 | * When it needs to get the threshold value, Eigen calls threshold(). By default, this
|
---|
243 | * uses a formula to automatically determine a reasonable threshold.
|
---|
244 | * Once you have called the present method setThreshold(const RealScalar&),
|
---|
245 | * your value is used instead.
|
---|
246 | *
|
---|
247 | * \param threshold The new value to use as the threshold.
|
---|
248 | *
|
---|
249 | * A pivot will be considered nonzero if its absolute value is strictly greater than
|
---|
250 | * \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$
|
---|
251 | * where maxpivot is the biggest pivot.
|
---|
252 | *
|
---|
253 | * If you want to come back to the default behavior, call setThreshold(Default_t)
|
---|
254 | */
|
---|
255 | FullPivLU& setThreshold(const RealScalar& threshold)
|
---|
256 | {
|
---|
257 | m_usePrescribedThreshold = true;
|
---|
258 | m_prescribedThreshold = threshold;
|
---|
259 | return *this;
|
---|
260 | }
|
---|
261 |
|
---|
262 | /** Allows to come back to the default behavior, letting Eigen use its default formula for
|
---|
263 | * determining the threshold.
|
---|
264 | *
|
---|
265 | * You should pass the special object Eigen::Default as parameter here.
|
---|
266 | * \code lu.setThreshold(Eigen::Default); \endcode
|
---|
267 | *
|
---|
268 | * See the documentation of setThreshold(const RealScalar&).
|
---|
269 | */
|
---|
270 | FullPivLU& setThreshold(Default_t)
|
---|
271 | {
|
---|
272 | m_usePrescribedThreshold = false;
|
---|
273 | return *this;
|
---|
274 | }
|
---|
275 |
|
---|
276 | /** Returns the threshold that will be used by certain methods such as rank().
|
---|
277 | *
|
---|
278 | * See the documentation of setThreshold(const RealScalar&).
|
---|
279 | */
|
---|
280 | RealScalar threshold() const
|
---|
281 | {
|
---|
282 | eigen_assert(m_isInitialized || m_usePrescribedThreshold);
|
---|
283 | return m_usePrescribedThreshold ? m_prescribedThreshold
|
---|
284 | // this formula comes from experimenting (see "LU precision tuning" thread on the list)
|
---|
285 | // and turns out to be identical to Higham's formula used already in LDLt.
|
---|
286 | : NumTraits<Scalar>::epsilon() * m_lu.diagonalSize();
|
---|
287 | }
|
---|
288 |
|
---|
289 | /** \returns the rank of the matrix of which *this is the LU decomposition.
|
---|
290 | *
|
---|
291 | * \note This method has to determine which pivots should be considered nonzero.
|
---|
292 | * For that, it uses the threshold value that you can control by calling
|
---|
293 | * setThreshold(const RealScalar&).
|
---|
294 | */
|
---|
295 | inline Index rank() const
|
---|
296 | {
|
---|
297 | using std::abs;
|
---|
298 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
299 | RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
|
---|
300 | Index result = 0;
|
---|
301 | for(Index i = 0; i < m_nonzero_pivots; ++i)
|
---|
302 | result += (abs(m_lu.coeff(i,i)) > premultiplied_threshold);
|
---|
303 | return result;
|
---|
304 | }
|
---|
305 |
|
---|
306 | /** \returns the dimension of the kernel of the matrix of which *this is the LU decomposition.
|
---|
307 | *
|
---|
308 | * \note This method has to determine which pivots should be considered nonzero.
|
---|
309 | * For that, it uses the threshold value that you can control by calling
|
---|
310 | * setThreshold(const RealScalar&).
|
---|
311 | */
|
---|
312 | inline Index dimensionOfKernel() const
|
---|
313 | {
|
---|
314 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
315 | return cols() - rank();
|
---|
316 | }
|
---|
317 |
|
---|
318 | /** \returns true if the matrix of which *this is the LU decomposition represents an injective
|
---|
319 | * linear map, i.e. has trivial kernel; false otherwise.
|
---|
320 | *
|
---|
321 | * \note This method has to determine which pivots should be considered nonzero.
|
---|
322 | * For that, it uses the threshold value that you can control by calling
|
---|
323 | * setThreshold(const RealScalar&).
|
---|
324 | */
|
---|
325 | inline bool isInjective() const
|
---|
326 | {
|
---|
327 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
328 | return rank() == cols();
|
---|
329 | }
|
---|
330 |
|
---|
331 | /** \returns true if the matrix of which *this is the LU decomposition represents a surjective
|
---|
332 | * linear map; false otherwise.
|
---|
333 | *
|
---|
334 | * \note This method has to determine which pivots should be considered nonzero.
|
---|
335 | * For that, it uses the threshold value that you can control by calling
|
---|
336 | * setThreshold(const RealScalar&).
|
---|
337 | */
|
---|
338 | inline bool isSurjective() const
|
---|
339 | {
|
---|
340 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
341 | return rank() == rows();
|
---|
342 | }
|
---|
343 |
|
---|
344 | /** \returns true if the matrix of which *this is the LU decomposition is invertible.
|
---|
345 | *
|
---|
346 | * \note This method has to determine which pivots should be considered nonzero.
|
---|
347 | * For that, it uses the threshold value that you can control by calling
|
---|
348 | * setThreshold(const RealScalar&).
|
---|
349 | */
|
---|
350 | inline bool isInvertible() const
|
---|
351 | {
|
---|
352 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
353 | return isInjective() && (m_lu.rows() == m_lu.cols());
|
---|
354 | }
|
---|
355 |
|
---|
356 | /** \returns the inverse of the matrix of which *this is the LU decomposition.
|
---|
357 | *
|
---|
358 | * \note If this matrix is not invertible, the returned matrix has undefined coefficients.
|
---|
359 | * Use isInvertible() to first determine whether this matrix is invertible.
|
---|
360 | *
|
---|
361 | * \sa MatrixBase::inverse()
|
---|
362 | */
|
---|
363 | inline const internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType> inverse() const
|
---|
364 | {
|
---|
365 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
366 | eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
|
---|
367 | return internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType>
|
---|
368 | (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
|
---|
369 | }
|
---|
370 |
|
---|
371 | MatrixType reconstructedMatrix() const;
|
---|
372 |
|
---|
373 | inline Index rows() const { return m_lu.rows(); }
|
---|
374 | inline Index cols() const { return m_lu.cols(); }
|
---|
375 |
|
---|
376 | protected:
|
---|
377 |
|
---|
378 | static void check_template_parameters()
|
---|
379 | {
|
---|
380 | EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
|
---|
381 | }
|
---|
382 |
|
---|
383 | MatrixType m_lu;
|
---|
384 | PermutationPType m_p;
|
---|
385 | PermutationQType m_q;
|
---|
386 | IntColVectorType m_rowsTranspositions;
|
---|
387 | IntRowVectorType m_colsTranspositions;
|
---|
388 | Index m_det_pq, m_nonzero_pivots;
|
---|
389 | RealScalar m_maxpivot, m_prescribedThreshold;
|
---|
390 | bool m_isInitialized, m_usePrescribedThreshold;
|
---|
391 | };
|
---|
392 |
|
---|
393 | template<typename MatrixType>
|
---|
394 | FullPivLU<MatrixType>::FullPivLU()
|
---|
395 | : m_isInitialized(false), m_usePrescribedThreshold(false)
|
---|
396 | {
|
---|
397 | }
|
---|
398 |
|
---|
399 | template<typename MatrixType>
|
---|
400 | FullPivLU<MatrixType>::FullPivLU(Index rows, Index cols)
|
---|
401 | : m_lu(rows, cols),
|
---|
402 | m_p(rows),
|
---|
403 | m_q(cols),
|
---|
404 | m_rowsTranspositions(rows),
|
---|
405 | m_colsTranspositions(cols),
|
---|
406 | m_isInitialized(false),
|
---|
407 | m_usePrescribedThreshold(false)
|
---|
408 | {
|
---|
409 | }
|
---|
410 |
|
---|
411 | template<typename MatrixType>
|
---|
412 | FullPivLU<MatrixType>::FullPivLU(const MatrixType& matrix)
|
---|
413 | : m_lu(matrix.rows(), matrix.cols()),
|
---|
414 | m_p(matrix.rows()),
|
---|
415 | m_q(matrix.cols()),
|
---|
416 | m_rowsTranspositions(matrix.rows()),
|
---|
417 | m_colsTranspositions(matrix.cols()),
|
---|
418 | m_isInitialized(false),
|
---|
419 | m_usePrescribedThreshold(false)
|
---|
420 | {
|
---|
421 | compute(matrix);
|
---|
422 | }
|
---|
423 |
|
---|
424 | template<typename MatrixType>
|
---|
425 | FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
|
---|
426 | {
|
---|
427 | check_template_parameters();
|
---|
428 |
|
---|
429 | // the permutations are stored as int indices, so just to be sure:
|
---|
430 | eigen_assert(matrix.rows()<=NumTraits<int>::highest() && matrix.cols()<=NumTraits<int>::highest());
|
---|
431 |
|
---|
432 | m_isInitialized = true;
|
---|
433 | m_lu = matrix;
|
---|
434 |
|
---|
435 | const Index size = matrix.diagonalSize();
|
---|
436 | const Index rows = matrix.rows();
|
---|
437 | const Index cols = matrix.cols();
|
---|
438 |
|
---|
439 | // will store the transpositions, before we accumulate them at the end.
|
---|
440 | // can't accumulate on-the-fly because that will be done in reverse order for the rows.
|
---|
441 | m_rowsTranspositions.resize(matrix.rows());
|
---|
442 | m_colsTranspositions.resize(matrix.cols());
|
---|
443 | Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
|
---|
444 |
|
---|
445 | m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
|
---|
446 | m_maxpivot = RealScalar(0);
|
---|
447 |
|
---|
448 | for(Index k = 0; k < size; ++k)
|
---|
449 | {
|
---|
450 | // First, we need to find the pivot.
|
---|
451 |
|
---|
452 | // biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
|
---|
453 | Index row_of_biggest_in_corner, col_of_biggest_in_corner;
|
---|
454 | RealScalar biggest_in_corner;
|
---|
455 | biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
|
---|
456 | .cwiseAbs()
|
---|
457 | .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
|
---|
458 | row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
|
---|
459 | col_of_biggest_in_corner += k; // need to add k to them.
|
---|
460 |
|
---|
461 | if(biggest_in_corner==RealScalar(0))
|
---|
462 | {
|
---|
463 | // before exiting, make sure to initialize the still uninitialized transpositions
|
---|
464 | // in a sane state without destroying what we already have.
|
---|
465 | m_nonzero_pivots = k;
|
---|
466 | for(Index i = k; i < size; ++i)
|
---|
467 | {
|
---|
468 | m_rowsTranspositions.coeffRef(i) = i;
|
---|
469 | m_colsTranspositions.coeffRef(i) = i;
|
---|
470 | }
|
---|
471 | break;
|
---|
472 | }
|
---|
473 |
|
---|
474 | if(biggest_in_corner > m_maxpivot) m_maxpivot = biggest_in_corner;
|
---|
475 |
|
---|
476 | // Now that we've found the pivot, we need to apply the row/col swaps to
|
---|
477 | // bring it to the location (k,k).
|
---|
478 |
|
---|
479 | m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
|
---|
480 | m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
|
---|
481 | if(k != row_of_biggest_in_corner) {
|
---|
482 | m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
|
---|
483 | ++number_of_transpositions;
|
---|
484 | }
|
---|
485 | if(k != col_of_biggest_in_corner) {
|
---|
486 | m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
|
---|
487 | ++number_of_transpositions;
|
---|
488 | }
|
---|
489 |
|
---|
490 | // Now that the pivot is at the right location, we update the remaining
|
---|
491 | // bottom-right corner by Gaussian elimination.
|
---|
492 |
|
---|
493 | if(k<rows-1)
|
---|
494 | m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
|
---|
495 | if(k<size-1)
|
---|
496 | m_lu.block(k+1,k+1,rows-k-1,cols-k-1).noalias() -= m_lu.col(k).tail(rows-k-1) * m_lu.row(k).tail(cols-k-1);
|
---|
497 | }
|
---|
498 |
|
---|
499 | // the main loop is over, we still have to accumulate the transpositions to find the
|
---|
500 | // permutations P and Q
|
---|
501 |
|
---|
502 | m_p.setIdentity(rows);
|
---|
503 | for(Index k = size-1; k >= 0; --k)
|
---|
504 | m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
|
---|
505 |
|
---|
506 | m_q.setIdentity(cols);
|
---|
507 | for(Index k = 0; k < size; ++k)
|
---|
508 | m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
|
---|
509 |
|
---|
510 | m_det_pq = (number_of_transpositions%2) ? -1 : 1;
|
---|
511 | return *this;
|
---|
512 | }
|
---|
513 |
|
---|
514 | template<typename MatrixType>
|
---|
515 | typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
|
---|
516 | {
|
---|
517 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
518 | eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
|
---|
519 | return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
|
---|
520 | }
|
---|
521 |
|
---|
522 | /** \returns the matrix represented by the decomposition,
|
---|
523 | * i.e., it returns the product: \f$ P^{-1} L U Q^{-1} \f$.
|
---|
524 | * This function is provided for debug purposes. */
|
---|
525 | template<typename MatrixType>
|
---|
526 | MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
|
---|
527 | {
|
---|
528 | eigen_assert(m_isInitialized && "LU is not initialized.");
|
---|
529 | const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols());
|
---|
530 | // LU
|
---|
531 | MatrixType res(m_lu.rows(),m_lu.cols());
|
---|
532 | // FIXME the .toDenseMatrix() should not be needed...
|
---|
533 | res = m_lu.leftCols(smalldim)
|
---|
534 | .template triangularView<UnitLower>().toDenseMatrix()
|
---|
535 | * m_lu.topRows(smalldim)
|
---|
536 | .template triangularView<Upper>().toDenseMatrix();
|
---|
537 |
|
---|
538 | // P^{-1}(LU)
|
---|
539 | res = m_p.inverse() * res;
|
---|
540 |
|
---|
541 | // (P^{-1}LU)Q^{-1}
|
---|
542 | res = res * m_q.inverse();
|
---|
543 |
|
---|
544 | return res;
|
---|
545 | }
|
---|
546 |
|
---|
547 | /********* Implementation of kernel() **************************************************/
|
---|
548 |
|
---|
549 | namespace internal {
|
---|
550 | template<typename _MatrixType>
|
---|
551 | struct kernel_retval<FullPivLU<_MatrixType> >
|
---|
552 | : kernel_retval_base<FullPivLU<_MatrixType> >
|
---|
553 | {
|
---|
554 | EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<_MatrixType>)
|
---|
555 |
|
---|
556 | enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
|
---|
557 | MatrixType::MaxColsAtCompileTime,
|
---|
558 | MatrixType::MaxRowsAtCompileTime)
|
---|
559 | };
|
---|
560 |
|
---|
561 | template<typename Dest> void evalTo(Dest& dst) const
|
---|
562 | {
|
---|
563 | using std::abs;
|
---|
564 | const Index cols = dec().matrixLU().cols(), dimker = cols - rank();
|
---|
565 | if(dimker == 0)
|
---|
566 | {
|
---|
567 | // The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's
|
---|
568 | // avoid crashing/asserting as that depends on floating point calculations. Let's
|
---|
569 | // just return a single column vector filled with zeros.
|
---|
570 | dst.setZero();
|
---|
571 | return;
|
---|
572 | }
|
---|
573 |
|
---|
574 | /* Let us use the following lemma:
|
---|
575 | *
|
---|
576 | * Lemma: If the matrix A has the LU decomposition PAQ = LU,
|
---|
577 | * then Ker A = Q(Ker U).
|
---|
578 | *
|
---|
579 | * Proof: trivial: just keep in mind that P, Q, L are invertible.
|
---|
580 | */
|
---|
581 |
|
---|
582 | /* Thus, all we need to do is to compute Ker U, and then apply Q.
|
---|
583 | *
|
---|
584 | * U is upper triangular, with eigenvalues sorted so that any zeros appear at the end.
|
---|
585 | * Thus, the diagonal of U ends with exactly
|
---|
586 | * dimKer zero's. Let us use that to construct dimKer linearly
|
---|
587 | * independent vectors in Ker U.
|
---|
588 | */
|
---|
589 |
|
---|
590 | Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
|
---|
591 | RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
|
---|
592 | Index p = 0;
|
---|
593 | for(Index i = 0; i < dec().nonzeroPivots(); ++i)
|
---|
594 | if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
|
---|
595 | pivots.coeffRef(p++) = i;
|
---|
596 | eigen_internal_assert(p == rank());
|
---|
597 |
|
---|
598 | // we construct a temporaty trapezoid matrix m, by taking the U matrix and
|
---|
599 | // permuting the rows and cols to bring the nonnegligible pivots to the top of
|
---|
600 | // the main diagonal. We need that to be able to apply our triangular solvers.
|
---|
601 | // FIXME when we get triangularView-for-rectangular-matrices, this can be simplified
|
---|
602 | Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
|
---|
603 | MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
|
---|
604 | m(dec().matrixLU().block(0, 0, rank(), cols));
|
---|
605 | for(Index i = 0; i < rank(); ++i)
|
---|
606 | {
|
---|
607 | if(i) m.row(i).head(i).setZero();
|
---|
608 | m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i);
|
---|
609 | }
|
---|
610 | m.block(0, 0, rank(), rank());
|
---|
611 | m.block(0, 0, rank(), rank()).template triangularView<StrictlyLower>().setZero();
|
---|
612 | for(Index i = 0; i < rank(); ++i)
|
---|
613 | m.col(i).swap(m.col(pivots.coeff(i)));
|
---|
614 |
|
---|
615 | // ok, we have our trapezoid matrix, we can apply the triangular solver.
|
---|
616 | // notice that the math behind this suggests that we should apply this to the
|
---|
617 | // negative of the RHS, but for performance we just put the negative sign elsewhere, see below.
|
---|
618 | m.topLeftCorner(rank(), rank())
|
---|
619 | .template triangularView<Upper>().solveInPlace(
|
---|
620 | m.topRightCorner(rank(), dimker)
|
---|
621 | );
|
---|
622 |
|
---|
623 | // now we must undo the column permutation that we had applied!
|
---|
624 | for(Index i = rank()-1; i >= 0; --i)
|
---|
625 | m.col(i).swap(m.col(pivots.coeff(i)));
|
---|
626 |
|
---|
627 | // see the negative sign in the next line, that's what we were talking about above.
|
---|
628 | for(Index i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker);
|
---|
629 | for(Index i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero();
|
---|
630 | for(Index k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1);
|
---|
631 | }
|
---|
632 | };
|
---|
633 |
|
---|
634 | /***** Implementation of image() *****************************************************/
|
---|
635 |
|
---|
636 | template<typename _MatrixType>
|
---|
637 | struct image_retval<FullPivLU<_MatrixType> >
|
---|
638 | : image_retval_base<FullPivLU<_MatrixType> >
|
---|
639 | {
|
---|
640 | EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<_MatrixType>)
|
---|
641 |
|
---|
642 | enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
|
---|
643 | MatrixType::MaxColsAtCompileTime,
|
---|
644 | MatrixType::MaxRowsAtCompileTime)
|
---|
645 | };
|
---|
646 |
|
---|
647 | template<typename Dest> void evalTo(Dest& dst) const
|
---|
648 | {
|
---|
649 | using std::abs;
|
---|
650 | if(rank() == 0)
|
---|
651 | {
|
---|
652 | // The Image is just {0}, so it doesn't have a basis properly speaking, but let's
|
---|
653 | // avoid crashing/asserting as that depends on floating point calculations. Let's
|
---|
654 | // just return a single column vector filled with zeros.
|
---|
655 | dst.setZero();
|
---|
656 | return;
|
---|
657 | }
|
---|
658 |
|
---|
659 | Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
|
---|
660 | RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
|
---|
661 | Index p = 0;
|
---|
662 | for(Index i = 0; i < dec().nonzeroPivots(); ++i)
|
---|
663 | if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
|
---|
664 | pivots.coeffRef(p++) = i;
|
---|
665 | eigen_internal_assert(p == rank());
|
---|
666 |
|
---|
667 | for(Index i = 0; i < rank(); ++i)
|
---|
668 | dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i)));
|
---|
669 | }
|
---|
670 | };
|
---|
671 |
|
---|
672 | /***** Implementation of solve() *****************************************************/
|
---|
673 |
|
---|
674 | template<typename _MatrixType, typename Rhs>
|
---|
675 | struct solve_retval<FullPivLU<_MatrixType>, Rhs>
|
---|
676 | : solve_retval_base<FullPivLU<_MatrixType>, Rhs>
|
---|
677 | {
|
---|
678 | EIGEN_MAKE_SOLVE_HELPERS(FullPivLU<_MatrixType>,Rhs)
|
---|
679 |
|
---|
680 | template<typename Dest> void evalTo(Dest& dst) const
|
---|
681 | {
|
---|
682 | /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
|
---|
683 | * So we proceed as follows:
|
---|
684 | * Step 1: compute c = P * rhs.
|
---|
685 | * Step 2: replace c by the solution x to Lx = c. Exists because L is invertible.
|
---|
686 | * Step 3: replace c by the solution x to Ux = c. May or may not exist.
|
---|
687 | * Step 4: result = Q * c;
|
---|
688 | */
|
---|
689 |
|
---|
690 | const Index rows = dec().rows(), cols = dec().cols(),
|
---|
691 | nonzero_pivots = dec().rank();
|
---|
692 | eigen_assert(rhs().rows() == rows);
|
---|
693 | const Index smalldim = (std::min)(rows, cols);
|
---|
694 |
|
---|
695 | if(nonzero_pivots == 0)
|
---|
696 | {
|
---|
697 | dst.setZero();
|
---|
698 | return;
|
---|
699 | }
|
---|
700 |
|
---|
701 | typename Rhs::PlainObject c(rhs().rows(), rhs().cols());
|
---|
702 |
|
---|
703 | // Step 1
|
---|
704 | c = dec().permutationP() * rhs();
|
---|
705 |
|
---|
706 | // Step 2
|
---|
707 | dec().matrixLU()
|
---|
708 | .topLeftCorner(smalldim,smalldim)
|
---|
709 | .template triangularView<UnitLower>()
|
---|
710 | .solveInPlace(c.topRows(smalldim));
|
---|
711 | if(rows>cols)
|
---|
712 | {
|
---|
713 | c.bottomRows(rows-cols)
|
---|
714 | -= dec().matrixLU().bottomRows(rows-cols)
|
---|
715 | * c.topRows(cols);
|
---|
716 | }
|
---|
717 |
|
---|
718 | // Step 3
|
---|
719 | dec().matrixLU()
|
---|
720 | .topLeftCorner(nonzero_pivots, nonzero_pivots)
|
---|
721 | .template triangularView<Upper>()
|
---|
722 | .solveInPlace(c.topRows(nonzero_pivots));
|
---|
723 |
|
---|
724 | // Step 4
|
---|
725 | for(Index i = 0; i < nonzero_pivots; ++i)
|
---|
726 | dst.row(dec().permutationQ().indices().coeff(i)) = c.row(i);
|
---|
727 | for(Index i = nonzero_pivots; i < dec().matrixLU().cols(); ++i)
|
---|
728 | dst.row(dec().permutationQ().indices().coeff(i)).setZero();
|
---|
729 | }
|
---|
730 | };
|
---|
731 |
|
---|
732 | } // end namespace internal
|
---|
733 |
|
---|
734 | /******* MatrixBase methods *****************************************************************/
|
---|
735 |
|
---|
736 | /** \lu_module
|
---|
737 | *
|
---|
738 | * \return the full-pivoting LU decomposition of \c *this.
|
---|
739 | *
|
---|
740 | * \sa class FullPivLU
|
---|
741 | */
|
---|
742 | template<typename Derived>
|
---|
743 | inline const FullPivLU<typename MatrixBase<Derived>::PlainObject>
|
---|
744 | MatrixBase<Derived>::fullPivLu() const
|
---|
745 | {
|
---|
746 | return FullPivLU<PlainObject>(eval());
|
---|
747 | }
|
---|
748 |
|
---|
749 | } // end namespace Eigen
|
---|
750 |
|
---|
751 | #endif // EIGEN_LU_H
|
---|