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