source: pacpussensors/trunk/Vislab/lib3dv/eigen/Eigen/src/LU/FullPivLU.h@ 136

Last change on this file since 136 was 136, checked in by ldecherf, 7 years ago

Doc
File size: 27.6 KB
Line 
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
13namespace 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 */
46template<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
393template<typename MatrixType>
394FullPivLU<MatrixType>::FullPivLU()
395 : m_isInitialized(false), m_usePrescribedThreshold(false)
396{
397}
398
399template<typename MatrixType>
400FullPivLU<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
411template<typename MatrixType>
412FullPivLU<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
424template<typename MatrixType>
425FullPivLU<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
514template<typename MatrixType>
515typename 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. */
525template<typename MatrixType>
526MatrixType 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
549namespace internal {
550template<typename _MatrixType>
551struct 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
636template<typename _MatrixType>
637struct 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
674template<typename _MatrixType, typename Rhs>
675struct 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 */
742template<typename Derived>
743inline const FullPivLU<typename MatrixBase<Derived>::PlainObject>
744MatrixBase<Derived>::fullPivLu() const
745{
746 return FullPivLU<PlainObject>(eval());
747}
748
749} // end namespace Eigen
750
751#endif // EIGEN_LU_H
Note: See TracBrowser for help on using the repository browser.