Main MRPT website > C++ reference
MRPT logo

PartialPivLU.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 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00006 //
00007 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
00025 
00026 #ifndef EIGEN_PARTIALLU_H
00027 #define EIGEN_PARTIALLU_H
00028 
00029 /** \ingroup LU_Module
00030   *
00031   * \class PartialPivLU
00032   *
00033   * \brief LU decomposition of a matrix with partial pivoting, and related features
00034   *
00035   * \param MatrixType the type of the matrix of which we are computing the LU decomposition
00036   *
00037   * This class represents a LU decomposition of a \b square \b invertible matrix, with partial pivoting: the matrix A
00038   * is decomposed as A = PLU where L is unit-lower-triangular, U is upper-triangular, and P
00039   * is a permutation matrix.
00040   *
00041   * Typically, partial pivoting LU decomposition is only considered numerically stable for square invertible
00042   * matrices. Thus LAPACK's dgesv and dgesvx require the matrix to be square and invertible. The present class
00043   * does the same. It will assert that the matrix is square, but it won't (actually it can't) check that the
00044   * matrix is invertible: it is your task to check that you only use this decomposition on invertible matrices.
00045   *
00046   * The guaranteed safe alternative, working for all matrices, is the full pivoting LU decomposition, provided
00047   * by class FullPivLU.
00048   *
00049   * This is \b not a rank-revealing LU decomposition. Many features are intentionally absent from this class,
00050   * such as rank computation. If you need these features, use class FullPivLU.
00051   *
00052   * This LU decomposition is suitable to invert invertible matrices. It is what MatrixBase::inverse() uses
00053   * in the general case.
00054   * On the other hand, it is \b not suitable to determine whether a given matrix is invertible.
00055   *
00056   * The data of the LU decomposition can be directly accessed through the methods matrixLU(), permutationP().
00057   *
00058   * \sa MatrixBase::partialPivLu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse(), class FullPivLU
00059   */
00060 template<typename _MatrixType> class PartialPivLU
00061 {
00062   public:
00063 
00064     typedef _MatrixType MatrixType;
00065     enum {
00066       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00067       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00068       Options = MatrixType::Options,
00069       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00070       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00071     };
00072     typedef typename MatrixType::Scalar Scalar;
00073     typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00074     typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
00075     typedef typename MatrixType::Index Index;
00076     typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
00077     typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
00078 
00079 
00080     /**
00081     * \brief Default Constructor.
00082     *
00083     * The default constructor is useful in cases in which the user intends to
00084     * perform decompositions via PartialPivLU::compute(const MatrixType&).
00085     */
00086     PartialPivLU();
00087 
00088     /** \brief Default Constructor with memory preallocation
00089       *
00090       * Like the default constructor but with preallocation of the internal data
00091       * according to the specified problem \a size.
00092       * \sa PartialPivLU()
00093       */
00094     PartialPivLU(Index size);
00095 
00096     /** Constructor.
00097       *
00098       * \param matrix the matrix of which to compute the LU decomposition.
00099       *
00100       * \warning The matrix should have full rank (e.g. if it's square, it should be invertible).
00101       * If you need to deal with non-full rank, use class FullPivLU instead.
00102       */
00103     PartialPivLU(const MatrixType& matrix);
00104 
00105     PartialPivLU& compute(const MatrixType& matrix);
00106 
00107     /** \returns the LU decomposition matrix: the upper-triangular part is U, the
00108       * unit-lower-triangular part is L (at least for square matrices; in the non-square
00109       * case, special care is needed, see the documentation of class FullPivLU).
00110       *
00111       * \sa matrixL(), matrixU()
00112       */
00113     inline const MatrixType& matrixLU() const
00114     {
00115       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00116       return m_lu;
00117     }
00118 
00119     /** \returns the permutation matrix P.
00120       */
00121     inline const PermutationType& permutationP() const
00122     {
00123       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00124       return m_p;
00125     }
00126 
00127     /** This method returns the solution x to the equation Ax=b, where A is the matrix of which
00128       * *this is the LU decomposition.
00129       *
00130       * \param b the right-hand-side of the equation to solve. Can be a vector or a matrix,
00131       *          the only requirement in order for the equation to make sense is that
00132       *          b.rows()==A.rows(), where A is the matrix of which *this is the LU decomposition.
00133       *
00134       * \returns the solution.
00135       *
00136       * Example: \include PartialPivLU_solve.cpp
00137       * Output: \verbinclude PartialPivLU_solve.out
00138       *
00139       * Since this PartialPivLU class assumes anyway that the matrix A is invertible, the solution
00140       * theoretically exists and is unique regardless of b.
00141       *
00142       * \sa TriangularView::solve(), inverse(), computeInverse()
00143       */
00144     template<typename Rhs>
00145     inline const internal::solve_retval<PartialPivLU, Rhs>
00146     solve(const MatrixBase<Rhs>& b) const
00147     {
00148       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00149       return internal::solve_retval<PartialPivLU, Rhs>(*this, b.derived());
00150     }
00151 
00152     /** \returns the inverse of the matrix of which *this is the LU decomposition.
00153       *
00154       * \warning The matrix being decomposed here is assumed to be invertible. If you need to check for
00155       *          invertibility, use class FullPivLU instead.
00156       *
00157       * \sa MatrixBase::inverse(), LU::inverse()
00158       */
00159     inline const internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType> inverse() const
00160     {
00161       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00162       return internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType>
00163                (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
00164     }
00165 
00166     /** \returns the determinant of the matrix of which
00167       * *this is the LU decomposition. It has only linear complexity
00168       * (that is, O(n) where n is the dimension of the square matrix)
00169       * as the LU decomposition has already been computed.
00170       *
00171       * \note For fixed-size matrices of size up to 4, MatrixBase::determinant() offers
00172       *       optimized paths.
00173       *
00174       * \warning a determinant can be very big or small, so for matrices
00175       * of large enough dimension, there is a risk of overflow/underflow.
00176       *
00177       * \sa MatrixBase::determinant()
00178       */
00179     typename internal::traits<MatrixType>::Scalar determinant() const;
00180 
00181     MatrixType reconstructedMatrix() const;
00182 
00183     inline Index rows() const { return m_lu.rows(); }
00184     inline Index cols() const { return m_lu.cols(); }
00185 
00186   protected:
00187     MatrixType m_lu;
00188     PermutationType m_p;
00189     TranspositionType m_rowsTranspositions;
00190     Index m_det_p;
00191     bool m_isInitialized;
00192 };
00193 
00194 template<typename MatrixType>
00195 PartialPivLU<MatrixType>::PartialPivLU()
00196   : m_lu(),
00197     m_p(),
00198     m_rowsTranspositions(),
00199     m_det_p(0),
00200     m_isInitialized(false)
00201 {
00202 }
00203 
00204 template<typename MatrixType>
00205 PartialPivLU<MatrixType>::PartialPivLU(Index size)
00206   : m_lu(size, size),
00207     m_p(size),
00208     m_rowsTranspositions(size),
00209     m_det_p(0),
00210     m_isInitialized(false)
00211 {
00212 }
00213 
00214 template<typename MatrixType>
00215 PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix)
00216   : m_lu(matrix.rows(), matrix.rows()),
00217     m_p(matrix.rows()),
00218     m_rowsTranspositions(matrix.rows()),
00219     m_det_p(0),
00220     m_isInitialized(false)
00221 {
00222   compute(matrix);
00223 }
00224 
00225 namespace internal {
00226 
00227 /** \internal This is the blocked version of fullpivlu_unblocked() */
00228 template<typename Scalar, int StorageOrder>
00229 struct partial_lu_impl
00230 {
00231   // FIXME add a stride to Map, so that the following mapping becomes easier,
00232   // another option would be to create an expression being able to automatically
00233   // warp any Map, Matrix, and Block expressions as a unique type, but since that's exactly
00234   // a Map + stride, why not adding a stride to Map, and convenient ctors from a Matrix,
00235   // and Block.
00236   typedef Map<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > MapLU;
00237   typedef Block<MapLU, Dynamic, Dynamic> MatrixType;
00238   typedef Block<MatrixType,Dynamic,Dynamic> BlockType;
00239   typedef typename MatrixType::RealScalar RealScalar;
00240   typedef typename MatrixType::Index Index;
00241 
00242   /** \internal performs the LU decomposition in-place of the matrix \a lu
00243     * using an unblocked algorithm.
00244     *
00245     * In addition, this function returns the row transpositions in the
00246     * vector \a row_transpositions which must have a size equal to the number
00247     * of columns of the matrix \a lu, and an integer \a nb_transpositions
00248     * which returns the actual number of transpositions.
00249     *
00250     * \returns false if some pivot is exactly zero, in which case the matrix is left with
00251     *          undefined coefficients (to avoid generating inf/nan values). Returns true
00252     *          otherwise.
00253     */
00254   static bool unblocked_lu(MatrixType& lu, Index* row_transpositions, Index& nb_transpositions)
00255   {
00256     const Index rows = lu.rows();
00257     const Index size = std::min(lu.rows(),lu.cols());
00258     nb_transpositions = 0;
00259     for(Index k = 0; k < size; ++k)
00260     {
00261       Index row_of_biggest_in_col;
00262       RealScalar biggest_in_corner
00263         = lu.col(k).tail(rows-k).cwiseAbs().maxCoeff(&row_of_biggest_in_col);
00264       row_of_biggest_in_col += k;
00265 
00266       if(biggest_in_corner == 0) // the pivot is exactly zero: the matrix is singular
00267       {
00268         // end quickly, avoid generating inf/nan values. Although in this unblocked_lu case
00269         // the result is still valid, there's no need to boast about it because
00270         // the blocked_lu code can't guarantee the same.
00271         // before exiting, make sure to initialize the still uninitialized row_transpositions
00272         // in a sane state without destroying what we already have.
00273         for(Index i = k; i < size; i++)
00274           row_transpositions[i] = i;
00275         return false;
00276       }
00277 
00278       row_transpositions[k] = row_of_biggest_in_col;
00279 
00280       if(k != row_of_biggest_in_col)
00281       {
00282         lu.row(k).swap(lu.row(row_of_biggest_in_col));
00283         ++nb_transpositions;
00284       }
00285 
00286       if(k<rows-1)
00287       {
00288         Index rrows = rows-k-1;
00289         Index rsize = size-k-1;
00290         lu.col(k).tail(rrows) /= lu.coeff(k,k);
00291         lu.bottomRightCorner(rrows,rsize).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rsize);
00292       }
00293     }
00294     return true;
00295   }
00296 
00297   /** \internal performs the LU decomposition in-place of the matrix represented
00298     * by the variables \a rows, \a cols, \a lu_data, and \a lu_stride using a
00299     * recursive, blocked algorithm.
00300     *
00301     * In addition, this function returns the row transpositions in the
00302     * vector \a row_transpositions which must have a size equal to the number
00303     * of columns of the matrix \a lu, and an integer \a nb_transpositions
00304     * which returns the actual number of transpositions.
00305     *
00306     * \returns false if some pivot is exactly zero, in which case the matrix is left with
00307     *          undefined coefficients (to avoid generating inf/nan values). Returns true
00308     *          otherwise.
00309     *
00310     * \note This very low level interface using pointers, etc. is to:
00311     *   1 - reduce the number of instanciations to the strict minimum
00312     *   2 - avoid infinite recursion of the instanciations with Block<Block<Block<...> > >
00313     */
00314   static bool blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, Index* row_transpositions, Index& nb_transpositions, Index maxBlockSize=256)
00315   {
00316     MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols);
00317     MatrixType lu(lu1,0,0,rows,cols);
00318 
00319     const Index size = std::min(rows,cols);
00320 
00321     // if the matrix is too small, no blocking:
00322     if(size<=16)
00323     {
00324       return unblocked_lu(lu, row_transpositions, nb_transpositions);
00325     }
00326 
00327     // automatically adjust the number of subdivisions to the size
00328     // of the matrix so that there is enough sub blocks:
00329     Index blockSize;
00330     {
00331       blockSize = size/8;
00332       blockSize = (blockSize/16)*16;
00333       blockSize = std::min(std::max(blockSize,Index(8)), maxBlockSize);
00334     }
00335 
00336     nb_transpositions = 0;
00337     for(Index k = 0; k < size; k+=blockSize)
00338     {
00339       Index bs = std::min(size-k,blockSize); // actual size of the block
00340       Index trows = rows - k - bs; // trailing rows
00341       Index tsize = size - k - bs; // trailing size
00342 
00343       // partition the matrix:
00344       //                          A00 | A01 | A02
00345       // lu  = A_0 | A_1 | A_2 =  A10 | A11 | A12
00346       //                          A20 | A21 | A22
00347       BlockType A_0(lu,0,0,rows,k);
00348       BlockType A_2(lu,0,k+bs,rows,tsize);
00349       BlockType A11(lu,k,k,bs,bs);
00350       BlockType A12(lu,k,k+bs,bs,tsize);
00351       BlockType A21(lu,k+bs,k,trows,bs);
00352       BlockType A22(lu,k+bs,k+bs,trows,tsize);
00353 
00354       Index nb_transpositions_in_panel;
00355       // recursively call the blocked LU algorithm on [A11^T A21^T]^T
00356       // with a very small blocking size:
00357       if(!blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
00358                    row_transpositions+k, nb_transpositions_in_panel, 16))
00359       {
00360         // end quickly with undefined coefficients, just avoid generating inf/nan values.
00361         // before exiting, make sure to initialize the still uninitialized row_transpositions
00362         // in a sane state without destroying what we already have.
00363         for(Index i=k; i<size; ++i)
00364           row_transpositions[i] = i;
00365         return false;
00366       }
00367       nb_transpositions += nb_transpositions_in_panel;
00368 
00369       // update permutations and apply them to A_0
00370       for(Index i=k; i<k+bs; ++i)
00371       {
00372         Index piv = (row_transpositions[i] += k);
00373         A_0.row(i).swap(A_0.row(piv));
00374       }
00375 
00376       if(trows)
00377       {
00378         // apply permutations to A_2
00379         for(Index i=k;i<k+bs; ++i)
00380           A_2.row(i).swap(A_2.row(row_transpositions[i]));
00381 
00382         // A12 = A11^-1 A12
00383         A11.template triangularView<UnitLower>().solveInPlace(A12);
00384 
00385         A22.noalias() -= A21 * A12;
00386       }
00387     }
00388     return true;
00389   }
00390 };
00391 
00392 /** \internal performs the LU decomposition with partial pivoting in-place.
00393   */
00394 template<typename MatrixType, typename TranspositionType>
00395 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename MatrixType::Index& nb_transpositions)
00396 {
00397   eigen_assert(lu.cols() == row_transpositions.size());
00398   eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
00399 
00400   partial_lu_impl
00401     <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor>
00402     ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
00403 }
00404 
00405 } // end namespace internal
00406 
00407 template<typename MatrixType>
00408 PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& matrix)
00409 {
00410   m_lu = matrix;
00411 
00412   eigen_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
00413   const Index size = matrix.rows();
00414 
00415   m_rowsTranspositions.resize(size);
00416 
00417   Index nb_transpositions;
00418   internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
00419   m_det_p = (nb_transpositions%2) ? -1 : 1;
00420 
00421   m_p = m_rowsTranspositions;
00422 
00423   m_isInitialized = true;
00424   return *this;
00425 }
00426 
00427 template<typename MatrixType>
00428 typename internal::traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const
00429 {
00430   eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
00431   return Scalar(m_det_p) * m_lu.diagonal().prod();
00432 }
00433 
00434 /** \returns the matrix represented by the decomposition,
00435  * i.e., it returns the product: P^{-1} L U.
00436  * This function is provided for debug purpose. */
00437 template<typename MatrixType>
00438 MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
00439 {
00440   eigen_assert(m_isInitialized && "LU is not initialized.");
00441   // LU
00442   MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
00443                  * m_lu.template triangularView<Upper>();
00444 
00445   // P^{-1}(LU)
00446   res = m_p.inverse() * res;
00447 
00448   return res;
00449 }
00450 
00451 /***** Implementation of solve() *****************************************************/
00452 
00453 namespace internal {
00454 
00455 template<typename _MatrixType, typename Rhs>
00456 struct solve_retval<PartialPivLU<_MatrixType>, Rhs>
00457   : solve_retval_base<PartialPivLU<_MatrixType>, Rhs>
00458 {
00459   EIGEN_MAKE_SOLVE_HELPERS(PartialPivLU<_MatrixType>,Rhs)
00460 
00461   template<typename Dest> void evalTo(Dest& dst) const
00462   {
00463     /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
00464     * So we proceed as follows:
00465     * Step 1: compute c = Pb.
00466     * Step 2: replace c by the solution x to Lx = c.
00467     * Step 3: replace c by the solution x to Ux = c.
00468     */
00469 
00470     eigen_assert(rhs().rows() == dec().matrixLU().rows());
00471 
00472     // Step 1
00473     dst = dec().permutationP() * rhs();
00474 
00475     // Step 2
00476     dec().matrixLU().template triangularView<UnitLower>().solveInPlace(dst);
00477 
00478     // Step 3
00479     dec().matrixLU().template triangularView<Upper>().solveInPlace(dst);
00480   }
00481 };
00482 
00483 } // end namespace internal
00484 
00485 /******** MatrixBase methods *******/
00486 
00487 /** \lu_module
00488   *
00489   * \return the partial-pivoting LU decomposition of \c *this.
00490   *
00491   * \sa class PartialPivLU
00492   */
00493 template<typename Derived>
00494 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
00495 MatrixBase<Derived>::partialPivLu() const
00496 {
00497   return PartialPivLU<PlainObject>(eval());
00498 }
00499 
00500 /** \lu_module
00501   *
00502   * Synonym of partialPivLu().
00503   *
00504   * \return the partial-pivoting LU decomposition of \c *this.
00505   *
00506   * \sa class PartialPivLU
00507   */
00508 template<typename Derived>
00509 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
00510 MatrixBase<Derived>::lu() const
00511 {
00512   return PartialPivLU<PlainObject>(eval());
00513 }
00514 
00515 #endif // EIGEN_PARTIALLU_H



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN: at Sat Mar 26 06:40:17 UTC 2011