Main MRPT website > C++ reference
MRPT logo

FullPivLU.h

Go to the documentation of this file.
00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
00005 //
00006 // Eigen is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 3 of the License, or (at your option) any later version.
00010 //
00011 // Alternatively, you can redistribute it and/or
00012 // modify it under the terms of the GNU General Public License as
00013 // published by the Free Software Foundation; either version 2 of
00014 // the License, or (at your option) any later version.
00015 //
00016 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00017 // // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00018 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00019 // GNU General Public License for more details.
00020 //
00021 // You should have received a copy of the GNU Lesser General Public
00022 // License and a copy of the GNU General Public License along with
00023 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00024 
00025 #ifndef EIGEN_LU_H
00026 #define EIGEN_LU_H
00027 
00028 /** \ingroup LU_Module
00029   *
00030   * \class FullPivLU
00031   *
00032   * \brief LU decomposition of a matrix with complete pivoting, and related features
00033   *
00034   * \param MatrixType the type of the matrix of which we are computing the LU decomposition
00035   *
00036   * This class represents a LU decomposition of any matrix, with complete pivoting: the matrix A
00037   * is decomposed as A = PLUQ where L is unit-lower-triangular, U is upper-triangular, and P and Q
00038   * are permutation matrices. This is a rank-revealing LU decomposition. The eigenvalues (diagonal
00039   * coefficients) of U are sorted in such a way that any zeros are at the end.
00040   *
00041   * This decomposition provides the generic approach to solving systems of linear equations, computing
00042   * the rank, invertibility, inverse, kernel, and determinant.
00043   *
00044   * This LU decomposition is very stable and well tested with large matrices. However there are use cases where the SVD
00045   * decomposition is inherently more stable and/or flexible. For example, when computing the kernel of a matrix,
00046   * working with the SVD allows to select the smallest singular values of the matrix, something that
00047   * the LU decomposition doesn't see.
00048   *
00049   * The data of the LU decomposition can be directly accessed through the methods matrixLU(),
00050   * permutationP(), permutationQ().
00051   *
00052   * As an exemple, here is how the original matrix can be retrieved:
00053   * \include class_FullPivLU.cpp
00054   * Output: \verbinclude class_FullPivLU.out
00055   *
00056   * \sa MatrixBase::fullPivLu(), MatrixBase::determinant(), MatrixBase::inverse()
00057   */
00058 template<typename _MatrixType> class FullPivLU
00059 {
00060   public:
00061     typedef _MatrixType MatrixType;
00062     enum {
00063       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00064       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00065       Options = MatrixType::Options,
00066       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00067       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00068     };
00069     typedef typename MatrixType::Scalar Scalar;
00070     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00071     typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
00072     typedef typename MatrixType::Index Index;
00073     typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
00074     typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
00075     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
00076     typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType;
00077 
00078     /**
00079       * \brief Default Constructor.
00080       *
00081       * The default constructor is useful in cases in which the user intends to
00082       * perform decompositions via LU::compute(const MatrixType&).
00083       */
00084     FullPivLU();
00085 
00086     /** \brief Default Constructor with memory preallocation
00087       *
00088       * Like the default constructor but with preallocation of the internal data
00089       * according to the specified problem \a size.
00090       * \sa FullPivLU()
00091       */
00092     FullPivLU(Index rows, Index cols);
00093 
00094     /** Constructor.
00095       *
00096       * \param matrix the matrix of which to compute the LU decomposition.
00097       *               It is required to be nonzero.
00098       */
00099     FullPivLU(const MatrixType& matrix);
00100 
00101     /** Computes the LU decomposition of the given matrix.
00102       *
00103       * \param matrix the matrix of which to compute the LU decomposition.
00104       *               It is required to be nonzero.
00105       *
00106       * \returns a reference to *this
00107       */
00108     FullPivLU& compute(const MatrixType& matrix);
00109 
00110     /** \returns the LU decomposition matrix: the upper-triangular part is U, the
00111       * unit-lower-triangular part is L (at least for square matrices; in the non-square
00112       * case, special care is needed, see the documentation of class FullPivLU).
00113       *
00114       * \sa matrixL(), matrixU()
00115       */
00116     inline const MatrixType& matrixLU() const
00117     {
00118       eigen_assert(m_isInitialized && "LU is not initialized.");
00119       return m_lu;
00120     }
00121 
00122     /** \returns the number of nonzero pivots in the LU decomposition.
00123       * Here nonzero is meant in the exact sense, not in a fuzzy sense.
00124       * So that notion isn't really intrinsically interesting, but it is
00125       * still useful when implementing algorithms.
00126       *
00127       * \sa rank()
00128       */
00129     inline Index nonzeroPivots() const
00130     {
00131       eigen_assert(m_isInitialized && "LU is not initialized.");
00132       return m_nonzero_pivots;
00133     }
00134 
00135     /** \returns the absolute value of the biggest pivot, i.e. the biggest
00136       *          diagonal coefficient of U.
00137       */
00138     RealScalar maxPivot() const { return m_maxpivot; }
00139 
00140     /** \returns the permutation matrix P
00141       *
00142       * \sa permutationQ()
00143       */
00144     inline const PermutationPType& permutationP() const
00145     {
00146       eigen_assert(m_isInitialized && "LU is not initialized.");
00147       return m_p;
00148     }
00149 
00150     /** \returns the permutation matrix Q
00151       *
00152       * \sa permutationP()
00153       */
00154     inline const PermutationQType& permutationQ() const
00155     {
00156       eigen_assert(m_isInitialized && "LU is not initialized.");
00157       return m_q;
00158     }
00159 
00160     /** \returns the kernel of the matrix, also called its null-space. The columns of the returned matrix
00161       * will form a basis of the kernel.
00162       *
00163       * \note If the kernel has dimension zero, then the returned matrix is a column-vector filled with zeros.
00164       *
00165       * \note This method has to determine which pivots should be considered nonzero.
00166       *       For that, it uses the threshold value that you can control by calling
00167       *       setThreshold(const RealScalar&).
00168       *
00169       * Example: \include FullPivLU_kernel.cpp
00170       * Output: \verbinclude FullPivLU_kernel.out
00171       *
00172       * \sa image()
00173       */
00174     inline const internal::kernel_retval<FullPivLU> kernel() const
00175     {
00176       eigen_assert(m_isInitialized && "LU is not initialized.");
00177       return internal::kernel_retval<FullPivLU>(*this);
00178     }
00179 
00180     /** \returns the image of the matrix, also called its column-space. The columns of the returned matrix
00181       * will form a basis of the kernel.
00182       *
00183       * \param originalMatrix the original matrix, of which *this is the LU decomposition.
00184       *                       The reason why it is needed to pass it here, is that this allows
00185       *                       a large optimization, as otherwise this method would need to reconstruct it
00186       *                       from the LU decomposition.
00187       *
00188       * \note If the image has dimension zero, then the returned matrix is a column-vector filled with zeros.
00189       *
00190       * \note This method has to determine which pivots should be considered nonzero.
00191       *       For that, it uses the threshold value that you can control by calling
00192       *       setThreshold(const RealScalar&).
00193       *
00194       * Example: \include FullPivLU_image.cpp
00195       * Output: \verbinclude FullPivLU_image.out
00196       *
00197       * \sa kernel()
00198       */
00199     inline const internal::image_retval<FullPivLU>
00200       image(const MatrixType& originalMatrix) const
00201     {
00202       eigen_assert(m_isInitialized && "LU is not initialized.");
00203       return internal::image_retval<FullPivLU>(*this, originalMatrix);
00204     }
00205 
00206     /** \return a solution x to the equation Ax=b, where A is the matrix of which
00207       * *this is the LU decomposition.
00208       *
00209       * \param b the right-hand-side of the equation to solve. Can be a vector or a matrix,
00210       *          the only requirement in order for the equation to make sense is that
00211       *          b.rows()==A.rows(), where A is the matrix of which *this is the LU decomposition.
00212       *
00213       * \returns a solution.
00214       *
00215       * \note_about_checking_solutions
00216       *
00217       * \note_about_arbitrary_choice_of_solution
00218       * \note_about_using_kernel_to_study_multiple_solutions
00219       *
00220       * Example: \include FullPivLU_solve.cpp
00221       * Output: \verbinclude FullPivLU_solve.out
00222       *
00223       * \sa TriangularView::solve(), kernel(), inverse()
00224       */
00225     template<typename Rhs>
00226     inline const internal::solve_retval<FullPivLU, Rhs>
00227     solve(const MatrixBase<Rhs>& b) const
00228     {
00229       eigen_assert(m_isInitialized && "LU is not initialized.");
00230       return internal::solve_retval<FullPivLU, Rhs>(*this, b.derived());
00231     }
00232 
00233     /** \returns the determinant of the matrix of which
00234       * *this is the LU decomposition. It has only linear complexity
00235       * (that is, O(n) where n is the dimension of the square matrix)
00236       * as the LU decomposition has already been computed.
00237       *
00238       * \note This is only for square matrices.
00239       *
00240       * \note For fixed-size matrices of size up to 4, MatrixBase::determinant() offers
00241       *       optimized paths.
00242       *
00243       * \warning a determinant can be very big or small, so for matrices
00244       * of large enough dimension, there is a risk of overflow/underflow.
00245       *
00246       * \sa MatrixBase::determinant()
00247       */
00248     typename internal::traits<MatrixType>::Scalar determinant() const;
00249 
00250     /** Allows to prescribe a threshold to be used by certain methods, such as rank(),
00251       * who need to determine when pivots are to be considered nonzero. This is not used for the
00252       * LU decomposition itself.
00253       *
00254       * When it needs to get the threshold value, Eigen calls threshold(). By default, this
00255       * uses a formula to automatically determine a reasonable threshold.
00256       * Once you have called the present method setThreshold(const RealScalar&),
00257       * your value is used instead.
00258       *
00259       * \param threshold The new value to use as the threshold.
00260       *
00261       * A pivot will be considered nonzero if its absolute value is strictly greater than
00262       *  \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$
00263       * where maxpivot is the biggest pivot.
00264       *
00265       * If you want to come back to the default behavior, call setThreshold(Default_t)
00266       */
00267     FullPivLU& setThreshold(const RealScalar& threshold)
00268     {
00269       m_usePrescribedThreshold = true;
00270       m_prescribedThreshold = threshold;
00271       return *this;
00272     }
00273 
00274     /** Allows to come back to the default behavior, letting Eigen use its default formula for
00275       * determining the threshold.
00276       *
00277       * You should pass the special object Eigen::Default as parameter here.
00278       * \code lu.setThreshold(Eigen::Default); \endcode
00279       *
00280       * See the documentation of setThreshold(const RealScalar&).
00281       */
00282     FullPivLU& setThreshold(Default_t)
00283     {
00284       m_usePrescribedThreshold = false;
00285     }
00286 
00287     /** Returns the threshold that will be used by certain methods such as rank().
00288       *
00289       * See the documentation of setThreshold(const RealScalar&).
00290       */
00291     RealScalar threshold() const
00292     {
00293       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00294       return m_usePrescribedThreshold ? m_prescribedThreshold
00295       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
00296       // and turns out to be identical to Higham's formula used already in LDLt.
00297                                       : NumTraits<Scalar>::epsilon() * m_lu.diagonalSize();
00298     }
00299 
00300     /** \returns the rank of the matrix of which *this is the LU decomposition.
00301       *
00302       * \note This method has to determine which pivots should be considered nonzero.
00303       *       For that, it uses the threshold value that you can control by calling
00304       *       setThreshold(const RealScalar&).
00305       */
00306     inline Index rank() const
00307     {
00308       eigen_assert(m_isInitialized && "LU is not initialized.");
00309       RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
00310       Index result = 0;
00311       for(Index i = 0; i < m_nonzero_pivots; ++i)
00312         result += (internal::abs(m_lu.coeff(i,i)) > premultiplied_threshold);
00313       return result;
00314     }
00315 
00316     /** \returns the dimension of the kernel of the matrix of which *this is the LU decomposition.
00317       *
00318       * \note This method has to determine which pivots should be considered nonzero.
00319       *       For that, it uses the threshold value that you can control by calling
00320       *       setThreshold(const RealScalar&).
00321       */
00322     inline Index dimensionOfKernel() const
00323     {
00324       eigen_assert(m_isInitialized && "LU is not initialized.");
00325       return cols() - rank();
00326     }
00327 
00328     /** \returns true if the matrix of which *this is the LU decomposition represents an injective
00329       *          linear map, i.e. has trivial kernel; false otherwise.
00330       *
00331       * \note This method has to determine which pivots should be considered nonzero.
00332       *       For that, it uses the threshold value that you can control by calling
00333       *       setThreshold(const RealScalar&).
00334       */
00335     inline bool isInjective() const
00336     {
00337       eigen_assert(m_isInitialized && "LU is not initialized.");
00338       return rank() == cols();
00339     }
00340 
00341     /** \returns true if the matrix of which *this is the LU decomposition represents a surjective
00342       *          linear map; false otherwise.
00343       *
00344       * \note This method has to determine which pivots should be considered nonzero.
00345       *       For that, it uses the threshold value that you can control by calling
00346       *       setThreshold(const RealScalar&).
00347       */
00348     inline bool isSurjective() const
00349     {
00350       eigen_assert(m_isInitialized && "LU is not initialized.");
00351       return rank() == rows();
00352     }
00353 
00354     /** \returns true if the matrix of which *this is the LU decomposition is invertible.
00355       *
00356       * \note This method has to determine which pivots should be considered nonzero.
00357       *       For that, it uses the threshold value that you can control by calling
00358       *       setThreshold(const RealScalar&).
00359       */
00360     inline bool isInvertible() const
00361     {
00362       eigen_assert(m_isInitialized && "LU is not initialized.");
00363       return isInjective() && (m_lu.rows() == m_lu.cols());
00364     }
00365 
00366     /** \returns the inverse of the matrix of which *this is the LU decomposition.
00367       *
00368       * \note If this matrix is not invertible, the returned matrix has undefined coefficients.
00369       *       Use isInvertible() to first determine whether this matrix is invertible.
00370       *
00371       * \sa MatrixBase::inverse()
00372       */
00373     inline const internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType> inverse() const
00374     {
00375       eigen_assert(m_isInitialized && "LU is not initialized.");
00376       eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
00377       return internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType>
00378                (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
00379     }
00380 
00381     MatrixType reconstructedMatrix() const;
00382 
00383     inline Index rows() const { return m_lu.rows(); }
00384     inline Index cols() const { return m_lu.cols(); }
00385 
00386   protected:
00387     MatrixType m_lu;
00388     PermutationPType m_p;
00389     PermutationQType m_q;
00390     IntColVectorType m_rowsTranspositions;
00391     IntRowVectorType m_colsTranspositions;
00392     Index m_det_pq, m_nonzero_pivots;
00393     RealScalar m_maxpivot, m_prescribedThreshold;
00394     bool m_isInitialized, m_usePrescribedThreshold;
00395 };
00396 
00397 template<typename MatrixType>
00398 FullPivLU<MatrixType>::FullPivLU()
00399   : m_isInitialized(false), m_usePrescribedThreshold(false)
00400 {
00401 }
00402 
00403 template<typename MatrixType>
00404 FullPivLU<MatrixType>::FullPivLU(Index rows, Index cols)
00405   : m_lu(rows, cols),
00406     m_p(rows),
00407     m_q(cols),
00408     m_rowsTranspositions(rows),
00409     m_colsTranspositions(cols),
00410     m_isInitialized(false),
00411     m_usePrescribedThreshold(false)
00412 {
00413 }
00414 
00415 template<typename MatrixType>
00416 FullPivLU<MatrixType>::FullPivLU(const MatrixType& matrix)
00417   : m_lu(matrix.rows(), matrix.cols()),
00418     m_p(matrix.rows()),
00419     m_q(matrix.cols()),
00420     m_rowsTranspositions(matrix.rows()),
00421     m_colsTranspositions(matrix.cols()),
00422     m_isInitialized(false),
00423     m_usePrescribedThreshold(false)
00424 {
00425   compute(matrix);
00426 }
00427 
00428 template<typename MatrixType>
00429 FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
00430 {
00431   m_isInitialized = true;
00432   m_lu = matrix;
00433 
00434   const Index size = matrix.diagonalSize();
00435   const Index rows = matrix.rows();
00436   const Index cols = matrix.cols();
00437 
00438   // will store the transpositions, before we accumulate them at the end.
00439   // can't accumulate on-the-fly because that will be done in reverse order for the rows.
00440   m_rowsTranspositions.resize(matrix.rows());
00441   m_colsTranspositions.resize(matrix.cols());
00442   Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i
00443 
00444   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
00445   m_maxpivot = RealScalar(0);
00446   RealScalar cutoff(0);
00447 
00448   for(Index k = 0; k < size; ++k)
00449   {
00450     // First, we need to find the pivot.
00451 
00452     // biggest coefficient in the remaining bottom-right corner (starting at row k, col k)
00453     Index row_of_biggest_in_corner, col_of_biggest_in_corner;
00454     RealScalar biggest_in_corner;
00455     biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
00456                         .cwiseAbs()
00457                         .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
00458     row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner,
00459     col_of_biggest_in_corner += k; // need to add k to them.
00460 
00461     // when k==0, biggest_in_corner is the biggest coeff absolute value in the original matrix
00462     if(k == 0) cutoff = biggest_in_corner * NumTraits<Scalar>::epsilon();
00463 
00464     // if the pivot (hence the corner) is "zero", terminate to avoid generating nan/inf values.
00465     // Notice that using an exact comparison (biggest_in_corner==0) here, as Golub-van Loan do in
00466     // their pseudo-code, results in numerical instability! The cutoff here has been validated
00467     // by running the unit test 'lu' with many repetitions.
00468     if(biggest_in_corner < cutoff)
00469     {
00470       // before exiting, make sure to initialize the still uninitialized transpositions
00471       // in a sane state without destroying what we already have.
00472       m_nonzero_pivots = k;
00473       for(Index i = k; i < size; ++i)
00474       {
00475         m_rowsTranspositions.coeffRef(i) = i;
00476         m_colsTranspositions.coeffRef(i) = i;
00477       }
00478       break;
00479     }
00480 
00481     if(biggest_in_corner > m_maxpivot) m_maxpivot = biggest_in_corner;
00482 
00483     // Now that we've found the pivot, we need to apply the row/col swaps to
00484     // bring it to the location (k,k).
00485 
00486     m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
00487     m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
00488     if(k != row_of_biggest_in_corner) {
00489       m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
00490       ++number_of_transpositions;
00491     }
00492     if(k != col_of_biggest_in_corner) {
00493       m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
00494       ++number_of_transpositions;
00495     }
00496 
00497     // Now that the pivot is at the right location, we update the remaining
00498     // bottom-right corner by Gaussian elimination.
00499 
00500     if(k<rows-1)
00501       m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
00502     if(k<size-1)
00503       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);
00504   }
00505 
00506   // the main loop is over, we still have to accumulate the transpositions to find the
00507   // permutations P and Q
00508 
00509   m_p.setIdentity(rows);
00510   for(Index k = size-1; k >= 0; --k)
00511     m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
00512 
00513   m_q.setIdentity(cols);
00514   for(Index k = 0; k < size; ++k)
00515     m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
00516 
00517   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00518   return *this;
00519 }
00520 
00521 template<typename MatrixType>
00522 typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
00523 {
00524   eigen_assert(m_isInitialized && "LU is not initialized.");
00525   eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
00526   return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
00527 }
00528 
00529 /** \returns the matrix represented by the decomposition,
00530  * i.e., it returns the product: P^{-1} L U Q^{-1}.
00531  * This function is provided for debug purpose. */
00532 template<typename MatrixType>
00533 MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
00534 {
00535   eigen_assert(m_isInitialized && "LU is not initialized.");
00536   const Index smalldim = std::min(m_lu.rows(), m_lu.cols());
00537   // LU
00538   MatrixType res(m_lu.rows(),m_lu.cols());
00539   // FIXME the .toDenseMatrix() should not be needed...
00540   res = m_lu.leftCols(smalldim)
00541             .template triangularView<UnitLower>().toDenseMatrix()
00542       * m_lu.topRows(smalldim)
00543             .template triangularView<Upper>().toDenseMatrix();
00544 
00545   // P^{-1}(LU)
00546   res = m_p.inverse() * res;
00547 
00548   // (P^{-1}LU)Q^{-1}
00549   res = res * m_q.inverse();
00550 
00551   return res;
00552 }
00553 
00554 /********* Implementation of kernel() **************************************************/
00555 
00556 namespace internal {
00557 template<typename _MatrixType>
00558 struct kernel_retval<FullPivLU<_MatrixType> >
00559   : kernel_retval_base<FullPivLU<_MatrixType> >
00560 {
00561   EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<_MatrixType>)
00562 
00563   enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
00564             MatrixType::MaxColsAtCompileTime,
00565             MatrixType::MaxRowsAtCompileTime)
00566   };
00567 
00568   template<typename Dest> void evalTo(Dest& dst) const
00569   {
00570     const Index cols = dec().matrixLU().cols(), dimker = cols - rank();
00571     if(dimker == 0)
00572     {
00573       // The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's
00574       // avoid crashing/asserting as that depends on floating point calculations. Let's
00575       // just return a single column vector filled with zeros.
00576       dst.setZero();
00577       return;
00578     }
00579 
00580     /* Let us use the following lemma:
00581       *
00582       * Lemma: If the matrix A has the LU decomposition PAQ = LU,
00583       * then Ker A = Q(Ker U).
00584       *
00585       * Proof: trivial: just keep in mind that P, Q, L are invertible.
00586       */
00587 
00588     /* Thus, all we need to do is to compute Ker U, and then apply Q.
00589       *
00590       * U is upper triangular, with eigenvalues sorted so that any zeros appear at the end.
00591       * Thus, the diagonal of U ends with exactly
00592       * dimKer zero's. Let us use that to construct dimKer linearly
00593       * independent vectors in Ker U.
00594       */
00595 
00596     Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
00597     RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
00598     Index p = 0;
00599     for(Index i = 0; i < dec().nonzeroPivots(); ++i)
00600       if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
00601         pivots.coeffRef(p++) = i;
00602     eigen_internal_assert(p == rank());
00603 
00604     // we construct a temporaty trapezoid matrix m, by taking the U matrix and
00605     // permuting the rows and cols to bring the nonnegligible pivots to the top of
00606     // the main diagonal. We need that to be able to apply our triangular solvers.
00607     // FIXME when we get triangularView-for-rectangular-matrices, this can be simplified
00608     Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
00609            MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
00610       m(dec().matrixLU().block(0, 0, rank(), cols));
00611     for(Index i = 0; i < rank(); ++i)
00612     {
00613       if(i) m.row(i).head(i).setZero();
00614       m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i);
00615     }
00616     m.block(0, 0, rank(), rank());
00617     m.block(0, 0, rank(), rank()).template triangularView<StrictlyLower>().setZero();
00618     for(Index i = 0; i < rank(); ++i)
00619       m.col(i).swap(m.col(pivots.coeff(i)));
00620 
00621     // ok, we have our trapezoid matrix, we can apply the triangular solver.
00622     // notice that the math behind this suggests that we should apply this to the
00623     // negative of the RHS, but for performance we just put the negative sign elsewhere, see below.
00624     m.topLeftCorner(rank(), rank())
00625      .template triangularView<Upper>().solveInPlace(
00626         m.topRightCorner(rank(), dimker)
00627       );
00628 
00629     // now we must undo the column permutation that we had applied!
00630     for(Index i = rank()-1; i >= 0; --i)
00631       m.col(i).swap(m.col(pivots.coeff(i)));
00632 
00633     // see the negative sign in the next line, that's what we were talking about above.
00634     for(Index i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker);
00635     for(Index i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero();
00636     for(Index k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1);
00637   }
00638 };
00639 
00640 /***** Implementation of image() *****************************************************/
00641 
00642 template<typename _MatrixType>
00643 struct image_retval<FullPivLU<_MatrixType> >
00644   : image_retval_base<FullPivLU<_MatrixType> >
00645 {
00646   EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<_MatrixType>)
00647 
00648   enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
00649             MatrixType::MaxColsAtCompileTime,
00650             MatrixType::MaxRowsAtCompileTime)
00651   };
00652 
00653   template<typename Dest> void evalTo(Dest& dst) const
00654   {
00655     if(rank() == 0)
00656     {
00657       // The Image is just {0}, so it doesn't have a basis properly speaking, but let's
00658       // avoid crashing/asserting as that depends on floating point calculations. Let's
00659       // just return a single column vector filled with zeros.
00660       dst.setZero();
00661       return;
00662     }
00663 
00664     Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
00665     RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
00666     Index p = 0;
00667     for(Index i = 0; i < dec().nonzeroPivots(); ++i)
00668       if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
00669         pivots.coeffRef(p++) = i;
00670     eigen_internal_assert(p == rank());
00671 
00672     for(Index i = 0; i < rank(); ++i)
00673       dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i)));
00674   }
00675 };
00676 
00677 /***** Implementation of solve() *****************************************************/
00678 
00679 template<typename _MatrixType, typename Rhs>
00680 struct solve_retval<FullPivLU<_MatrixType>, Rhs>
00681   : solve_retval_base<FullPivLU<_MatrixType>, Rhs>
00682 {
00683   EIGEN_MAKE_SOLVE_HELPERS(FullPivLU<_MatrixType>,Rhs)
00684 
00685   template<typename Dest> void evalTo(Dest& dst) const
00686   {
00687     /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}.
00688      * So we proceed as follows:
00689      * Step 1: compute c = P * rhs.
00690      * Step 2: replace c by the solution x to Lx = c. Exists because L is invertible.
00691      * Step 3: replace c by the solution x to Ux = c. May or may not exist.
00692      * Step 4: result = Q * c;
00693      */
00694 
00695     const Index rows = dec().rows(), cols = dec().cols(),
00696               nonzero_pivots = dec().nonzeroPivots();
00697     eigen_assert(rhs().rows() == rows);
00698     const Index smalldim = std::min(rows, cols);
00699 
00700     if(nonzero_pivots == 0)
00701     {
00702       dst.setZero();
00703       return;
00704     }
00705 
00706     typename Rhs::PlainObject c(rhs().rows(), rhs().cols());
00707 
00708     // Step 1
00709     c = dec().permutationP() * rhs();
00710 
00711     // Step 2
00712     dec().matrixLU()
00713         .topLeftCorner(smalldim,smalldim)
00714         .template triangularView<UnitLower>()
00715         .solveInPlace(c.topRows(smalldim));
00716     if(rows>cols)
00717     {
00718       c.bottomRows(rows-cols)
00719         -= dec().matrixLU().bottomRows(rows-cols)
00720          * c.topRows(cols);
00721     }
00722 
00723     // Step 3
00724     dec().matrixLU()
00725         .topLeftCorner(nonzero_pivots, nonzero_pivots)
00726         .template triangularView<Upper>()
00727         .solveInPlace(c.topRows(nonzero_pivots));
00728 
00729     // Step 4
00730     for(Index i = 0; i < nonzero_pivots; ++i)
00731       dst.row(dec().permutationQ().indices().coeff(i)) = c.row(i);
00732     for(Index i = nonzero_pivots; i < dec().matrixLU().cols(); ++i)
00733       dst.row(dec().permutationQ().indices().coeff(i)).setZero();
00734   }
00735 };
00736 
00737 } // end namespace internal
00738 
00739 /******* MatrixBase methods *****************************************************************/
00740 
00741 /** \lu_module
00742   *
00743   * \return the full-pivoting LU decomposition of \c *this.
00744   *
00745   * \sa class FullPivLU
00746   */
00747 template<typename Derived>
00748 inline const FullPivLU<typename MatrixBase<Derived>::PlainObject>
00749 MatrixBase<Derived>::fullPivLu() const
00750 {
00751   return FullPivLU<PlainObject>(eval());
00752 }
00753 
00754 #endif // EIGEN_LU_H



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN:exported at Tue Jan 25 21:56:31 UTC 2011