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