src/LU/PartialPivLU.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2006-2009 Benoit Jacob <jacob.benoit.1@gmail.com>
0005 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_PARTIALLU_H
0012 #define EIGEN_PARTIALLU_H
0013
0014 namespace Eigen {
0015
0016 namespace internal {
0017 template<typename _MatrixType> struct traits<PartialPivLU<_MatrixType> >
0018  : traits<_MatrixType>
0019 {
0020   typedef MatrixXpr XprKind;
0021   typedef SolverStorage StorageKind;
0022   typedef int StorageIndex;
0023   typedef traits<_MatrixType> BaseTraits;
0024   enum {
0025     Flags = BaseTraits::Flags & RowMajorBit,
0026     CoeffReadCost = Dynamic
0027   };
0028 };
0029
0030 template<typename T,typename Derived>
0031 struct enable_if_ref;
0032 // {
0033 //   typedef Derived type;
0034 // };
0035
0036 template<typename T,typename Derived>
0037 struct enable_if_ref<Ref<T>,Derived> {
0038   typedef Derived type;
0039 };
0040
0041 } // end namespace internal
0042
0043 /** \ingroup LU_Module
0044   *
0045   * \class PartialPivLU
0046   *
0047   * \brief LU decomposition of a matrix with partial pivoting, and related features
0048   *
0049   * \tparam _MatrixType the type of the matrix of which we are computing the LU decomposition
0050   *
0051   * This class represents a LU decomposition of a \b square \b invertible matrix, with partial pivoting: the matrix A
0052   * is decomposed as A = PLU where L is unit-lower-triangular, U is upper-triangular, and P
0053   * is a permutation matrix.
0054   *
0055   * Typically, partial pivoting LU decomposition is only considered numerically stable for square invertible
0056   * matrices. Thus LAPACK's dgesv and dgesvx require the matrix to be square and invertible. The present class
0057   * does the same. It will assert that the matrix is square, but it won't (actually it can't) check that the
0058   * matrix is invertible: it is your task to check that you only use this decomposition on invertible matrices.
0059   *
0060   * The guaranteed safe alternative, working for all matrices, is the full pivoting LU decomposition, provided
0061   * by class FullPivLU.
0062   *
0063   * This is \b not a rank-revealing LU decomposition. Many features are intentionally absent from this class,
0064   * such as rank computation. If you need these features, use class FullPivLU.
0065   *
0066   * This LU decomposition is suitable to invert invertible matrices. It is what MatrixBase::inverse() uses
0067   * in the general case.
0068   * On the other hand, it is \b not suitable to determine whether a given matrix is invertible.
0069   *
0070   * The data of the LU decomposition can be directly accessed through the methods matrixLU(), permutationP().
0071   *
0072   * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
0073   *
0074   * \sa MatrixBase::partialPivLu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse(), class FullPivLU
0075   */
0076 template<typename _MatrixType> class PartialPivLU
0077   : public SolverBase<PartialPivLU<_MatrixType> >
0078 {
0079   public:
0080
0081     typedef _MatrixType MatrixType;
0082     typedef SolverBase<PartialPivLU> Base;
0083     friend class SolverBase<PartialPivLU>;
0084
0085     EIGEN_GENERIC_PUBLIC_INTERFACE(PartialPivLU)
0086     enum {
0087       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
0088       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
0089     };
0090     typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
0091     typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
0092     typedef typename MatrixType::PlainObject PlainObject;
0093
0094     /**
0095       * \brief Default Constructor.
0096       *
0097       * The default constructor is useful in cases in which the user intends to
0098       * perform decompositions via PartialPivLU::compute(const MatrixType&).
0099       */
0100     PartialPivLU();
0101
0102     /** \brief Default Constructor with memory preallocation
0103       *
0104       * Like the default constructor but with preallocation of the internal data
0105       * according to the specified problem \a size.
0106       * \sa PartialPivLU()
0107       */
0108     explicit PartialPivLU(Index size);
0109
0110     /** Constructor.
0111       *
0112       * \param matrix the matrix of which to compute the LU decomposition.
0113       *
0114       * \warning The matrix should have full rank (e.g. if it's square, it should be invertible).
0115       * If you need to deal with non-full rank, use class FullPivLU instead.
0116       */
0117     template<typename InputType>
0118     explicit PartialPivLU(const EigenBase<InputType>& matrix);
0119
0120     /** Constructor for \link InplaceDecomposition inplace decomposition \endlink
0121       *
0122       * \param matrix the matrix of which to compute the LU decomposition.
0123       *
0124       * \warning The matrix should have full rank (e.g. if it's square, it should be invertible).
0125       * If you need to deal with non-full rank, use class FullPivLU instead.
0126       */
0127     template<typename InputType>
0128     explicit PartialPivLU(EigenBase<InputType>& matrix);
0129
0130     template<typename InputType>
0131     PartialPivLU& compute(const EigenBase<InputType>& matrix) {
0132       m_lu = matrix.derived();
0133       compute();
0134       return *this;
0135     }
0136
0137     /** \returns the LU decomposition matrix: the upper-triangular part is U, the
0138       * unit-lower-triangular part is L (at least for square matrices; in the non-square
0139       * case, special care is needed, see the documentation of class FullPivLU).
0140       *
0141       * \sa matrixL(), matrixU()
0142       */
0143     inline const MatrixType& matrixLU() const
0144     {
0145       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
0146       return m_lu;
0147     }
0148
0149     /** \returns the permutation matrix P.
0150       */
0151     inline const PermutationType& permutationP() const
0152     {
0153       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
0154       return m_p;
0155     }
0156
0157     #ifdef EIGEN_PARSED_BY_DOXYGEN
0158     /** This method returns the solution x to the equation Ax=b, where A is the matrix of which
0159       * *this is the LU decomposition.
0160       *
0161       * \param b the right-hand-side of the equation to solve. Can be a vector or a matrix,
0162       *          the only requirement in order for the equation to make sense is that
0163       *          b.rows()==A.rows(), where A is the matrix of which *this is the LU decomposition.
0164       *
0165       * \returns the solution.
0166       *
0167       * Example: \include PartialPivLU_solve.cpp
0168       * Output: \verbinclude PartialPivLU_solve.out
0169       *
0170       * Since this PartialPivLU class assumes anyway that the matrix A is invertible, the solution
0171       * theoretically exists and is unique regardless of b.
0172       *
0173       * \sa TriangularView::solve(), inverse(), computeInverse()
0174       */
0175     template<typename Rhs>
0176     inline const Solve<PartialPivLU, Rhs>
0177     solve(const MatrixBase<Rhs>& b) const;
0178     #endif
0179
0180     /** \returns an estimate of the reciprocal condition number of the matrix of which \c *this is
0181         the LU decomposition.
0182       */
0183     inline RealScalar rcond() const
0184     {
0185       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
0186       return internal::rcond_estimate_helper(m_l1_norm, *this);
0187     }
0188
0189     /** \returns the inverse of the matrix of which *this is the LU decomposition.
0190       *
0191       * \warning The matrix being decomposed here is assumed to be invertible. If you need to check for
0192       *          invertibility, use class FullPivLU instead.
0193       *
0194       * \sa MatrixBase::inverse(), LU::inverse()
0195       */
0196     inline const Inverse<PartialPivLU> inverse() const
0197     {
0198       eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
0199       return Inverse<PartialPivLU>(*this);
0200     }
0201
0202     /** \returns the determinant of the matrix of which
0203       * *this is the LU decomposition. It has only linear complexity
0204       * (that is, O(n) where n is the dimension of the square matrix)
0205       * as the LU decomposition has already been computed.
0206       *
0207       * \note For fixed-size matrices of size up to 4, MatrixBase::determinant() offers
0208       *       optimized paths.
0209       *
0210       * \warning a determinant can be very big or small, so for matrices
0211       * of large enough dimension, there is a risk of overflow/underflow.
0212       *
0213       * \sa MatrixBase::determinant()
0214       */
0215     Scalar determinant() const;
0216
0217     MatrixType reconstructedMatrix() const;
0218
0219     EIGEN_CONSTEXPR inline Index rows() const EIGEN_NOEXCEPT { return m_lu.rows(); }
0220     EIGEN_CONSTEXPR inline Index cols() const EIGEN_NOEXCEPT { return m_lu.cols(); }
0221
0222     #ifndef EIGEN_PARSED_BY_DOXYGEN
0223     template<typename RhsType, typename DstType>
0224     EIGEN_DEVICE_FUNC
0225     void _solve_impl(const RhsType &rhs, DstType &dst) const {
0226      /* The decomposition PA = LU can be rewritten as A = P^{-1} L U.
0227       * So we proceed as follows:
0228       * Step 1: compute c = Pb.
0229       * Step 2: replace c by the solution x to Lx = c.
0230       * Step 3: replace c by the solution x to Ux = c.
0231       */
0232
0233       // Step 1
0234       dst = permutationP() * rhs;
0235
0236       // Step 2
0237       m_lu.template triangularView<UnitLower>().solveInPlace(dst);
0238
0239       // Step 3
0240       m_lu.template triangularView<Upper>().solveInPlace(dst);
0241     }
0242
0243     template<bool Conjugate, typename RhsType, typename DstType>
0244     EIGEN_DEVICE_FUNC
0245     void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const {
0246      /* The decomposition PA = LU can be rewritten as A^T = U^T L^T P.
0247       * So we proceed as follows:
0248       * Step 1: compute c as the solution to L^T c = b
0249       * Step 2: replace c by the solution x to U^T x = c.
0250       * Step 3: update  c = P^-1 c.
0251       */
0252
0253       eigen_assert(rhs.rows() == m_lu.cols());
0254
0255       // Step 1
0256       dst = m_lu.template triangularView<Upper>().transpose()
0257                 .template conjugateIf<Conjugate>().solve(rhs);
0258       // Step 2
0259       m_lu.template triangularView<UnitLower>().transpose()
0260           .template conjugateIf<Conjugate>().solveInPlace(dst);
0261       // Step 3
0262       dst = permutationP().transpose() * dst;
0263     }
0264     #endif
0265
0266   protected:
0267
0268     static void check_template_parameters()
0269     {
0270       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
0271     }
0272
0273     void compute();
0274
0275     MatrixType m_lu;
0276     PermutationType m_p;
0277     TranspositionType m_rowsTranspositions;
0278     RealScalar m_l1_norm;
0279     signed char m_det_p;
0280     bool m_isInitialized;
0281 };
0282
0283 template<typename MatrixType>
0284 PartialPivLU<MatrixType>::PartialPivLU()
0285   : m_lu(),
0286     m_p(),
0287     m_rowsTranspositions(),
0288     m_l1_norm(0),
0289     m_det_p(0),
0290     m_isInitialized(false)
0291 {
0292 }
0293
0294 template<typename MatrixType>
0295 PartialPivLU<MatrixType>::PartialPivLU(Index size)
0296   : m_lu(size, size),
0297     m_p(size),
0298     m_rowsTranspositions(size),
0299     m_l1_norm(0),
0300     m_det_p(0),
0301     m_isInitialized(false)
0302 {
0303 }
0304
0305 template<typename MatrixType>
0306 template<typename InputType>
0307 PartialPivLU<MatrixType>::PartialPivLU(const EigenBase<InputType>& matrix)
0308   : m_lu(matrix.rows(),matrix.cols()),
0309     m_p(matrix.rows()),
0310     m_rowsTranspositions(matrix.rows()),
0311     m_l1_norm(0),
0312     m_det_p(0),
0313     m_isInitialized(false)
0314 {
0315   compute(matrix.derived());
0316 }
0317
0318 template<typename MatrixType>
0319 template<typename InputType>
0320 PartialPivLU<MatrixType>::PartialPivLU(EigenBase<InputType>& matrix)
0321   : m_lu(matrix.derived()),
0322     m_p(matrix.rows()),
0323     m_rowsTranspositions(matrix.rows()),
0324     m_l1_norm(0),
0325     m_det_p(0),
0326     m_isInitialized(false)
0327 {
0328   compute();
0329 }
0330
0331 namespace internal {
0332
0333 /** \internal This is the blocked version of fullpivlu_unblocked() */
0334 template<typename Scalar, int StorageOrder, typename PivIndex, int SizeAtCompileTime=Dynamic>
0335 struct partial_lu_impl
0336 {
0337   static const int UnBlockedBound = 16;
0338   static const bool UnBlockedAtCompileTime = SizeAtCompileTime!=Dynamic && SizeAtCompileTime<=UnBlockedBound;
0339   static const int ActualSizeAtCompileTime = UnBlockedAtCompileTime ? SizeAtCompileTime : Dynamic;
0340   // Remaining rows and columns at compile-time:
0341   static const int RRows = SizeAtCompileTime==2 ? 1 : Dynamic;
0342   static const int RCols = SizeAtCompileTime==2 ? 1 : Dynamic;
0343   typedef Matrix<Scalar, ActualSizeAtCompileTime, ActualSizeAtCompileTime, StorageOrder> MatrixType;
0344   typedef Ref<MatrixType> MatrixTypeRef;
0345   typedef Ref<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > BlockType;
0346   typedef typename MatrixType::RealScalar RealScalar;
0347
0348   /** \internal performs the LU decomposition in-place of the matrix \a lu
0349     * using an unblocked algorithm.
0350     *
0351     * In addition, this function returns the row transpositions in the
0352     * vector \a row_transpositions which must have a size equal to the number
0353     * of columns of the matrix \a lu, and an integer \a nb_transpositions
0354     * which returns the actual number of transpositions.
0355     *
0356     * \returns The index of the first pivot which is exactly zero if any, or a negative number otherwise.
0357     */
0358   static Index unblocked_lu(MatrixTypeRef& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
0359   {
0360     typedef scalar_score_coeff_op<Scalar> Scoring;
0361     typedef typename Scoring::result_type Score;
0362     const Index rows = lu.rows();
0363     const Index cols = lu.cols();
0364     const Index size = (std::min)(rows,cols);
0365     // For small compile-time matrices it is worth processing the last row separately:
0366     //  speedup: +100% for 2x2, +10% for others.
0367     const Index endk = UnBlockedAtCompileTime ? size-1 : size;
0368     nb_transpositions = 0;
0369     Index first_zero_pivot = -1;
0370     for(Index k = 0; k < endk; ++k)
0371     {
0372       int rrows = internal::convert_index<int>(rows-k-1);
0373       int rcols = internal::convert_index<int>(cols-k-1);
0374
0375       Index row_of_biggest_in_col;
0376       Score biggest_in_corner
0377         = lu.col(k).tail(rows-k).unaryExpr(Scoring()).maxCoeff(&row_of_biggest_in_col);
0378       row_of_biggest_in_col += k;
0379
0380       row_transpositions[k] = PivIndex(row_of_biggest_in_col);
0381
0382       if(biggest_in_corner != Score(0))
0383       {
0384         if(k != row_of_biggest_in_col)
0385         {
0386           lu.row(k).swap(lu.row(row_of_biggest_in_col));
0387           ++nb_transpositions;
0388         }
0389
0390         lu.col(k).tail(fix<RRows>(rrows)) /= lu.coeff(k,k);
0391       }
0392       else if(first_zero_pivot==-1)
0393       {
0394         // the pivot is exactly zero, we record the index of the first pivot which is exactly 0,
0395         // and continue the factorization such we still have A = PLU
0396         first_zero_pivot = k;
0397       }
0398
0399       if(k<rows-1)
0400         lu.bottomRightCorner(fix<RRows>(rrows),fix<RCols>(rcols)).noalias() -= lu.col(k).tail(fix<RRows>(rrows)) * lu.row(k).tail(fix<RCols>(rcols));
0401     }
0402
0403     // special handling of the last entry
0404     if(UnBlockedAtCompileTime)
0405     {
0406       Index k = endk;
0407       row_transpositions[k] = PivIndex(k);
0408       if (Scoring()(lu(k, k)) == Score(0) && first_zero_pivot == -1)
0409         first_zero_pivot = k;
0410     }
0411
0412     return first_zero_pivot;
0413   }
0414
0415   /** \internal performs the LU decomposition in-place of the matrix represented
0416     * by the variables \a rows, \a cols, \a lu_data, and \a lu_stride using a
0417     * recursive, blocked algorithm.
0418     *
0419     * In addition, this function returns the row transpositions in the
0420     * vector \a row_transpositions which must have a size equal to the number
0421     * of columns of the matrix \a lu, and an integer \a nb_transpositions
0422     * which returns the actual number of transpositions.
0423     *
0424     * \returns The index of the first pivot which is exactly zero if any, or a negative number otherwise.
0425     *
0426     * \note This very low level interface using pointers, etc. is to:
0427     *   1 - reduce the number of instantiations to the strict minimum
0428     *   2 - avoid infinite recursion of the instantiations with Block<Block<Block<...> > >
0429     */
0430   static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256)
0431   {
0432     MatrixTypeRef lu = MatrixType::Map(lu_data,rows, cols, OuterStride<>(luStride));
0433
0434     const Index size = (std::min)(rows,cols);
0435
0436     // if the matrix is too small, no blocking:
0437     if(UnBlockedAtCompileTime || size<=UnBlockedBound)
0438     {
0439       return unblocked_lu(lu, row_transpositions, nb_transpositions);
0440     }
0441
0442     // automatically adjust the number of subdivisions to the size
0443     // of the matrix so that there is enough sub blocks:
0444     Index blockSize;
0445     {
0446       blockSize = size/8;
0447       blockSize = (blockSize/16)*16;
0448       blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
0449     }
0450
0451     nb_transpositions = 0;
0452     Index first_zero_pivot = -1;
0453     for(Index k = 0; k < size; k+=blockSize)
0454     {
0455       Index bs = (std::min)(size-k,blockSize); // actual size of the block
0456       Index trows = rows - k - bs; // trailing rows
0457       Index tsize = size - k - bs; // trailing size
0458
0459       // partition the matrix:
0460       //                          A00 | A01 | A02
0461       // lu  = A_0 | A_1 | A_2 =  A10 | A11 | A12
0462       //                          A20 | A21 | A22
0463       BlockType A_0 = lu.block(0,0,rows,k);
0464       BlockType A_2 = lu.block(0,k+bs,rows,tsize);
0465       BlockType A11 = lu.block(k,k,bs,bs);
0466       BlockType A12 = lu.block(k,k+bs,bs,tsize);
0467       BlockType A21 = lu.block(k+bs,k,trows,bs);
0468       BlockType A22 = lu.block(k+bs,k+bs,trows,tsize);
0469
0470       PivIndex nb_transpositions_in_panel;
0471       // recursively call the blocked LU algorithm on [A11^T A21^T]^T
0472       // with a very small blocking size:
0473       Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
0474                    row_transpositions+k, nb_transpositions_in_panel, 16);
0475       if(ret>=0 && first_zero_pivot==-1)
0476         first_zero_pivot = k+ret;
0477
0478       nb_transpositions += nb_transpositions_in_panel;
0479       // update permutations and apply them to A_0
0480       for(Index i=k; i<k+bs; ++i)
0481       {
0482         Index piv = (row_transpositions[i] += internal::convert_index<PivIndex>(k));
0483         A_0.row(i).swap(A_0.row(piv));
0484       }
0485
0486       if(trows)
0487       {
0488         // apply permutations to A_2
0489         for(Index i=k;i<k+bs; ++i)
0490           A_2.row(i).swap(A_2.row(row_transpositions[i]));
0491
0492         // A12 = A11^-1 A12
0493         A11.template triangularView<UnitLower>().solveInPlace(A12);
0494
0495         A22.noalias() -= A21 * A12;
0496       }
0497     }
0498     return first_zero_pivot;
0499   }
0500 };
0501
0502 /** \internal performs the LU decomposition with partial pivoting in-place.
0503   */
0504 template<typename MatrixType, typename TranspositionType>
0505 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::StorageIndex& nb_transpositions)
0506 {
0507   // Special-case of zero matrix.
0508   if (lu.rows() == 0 || lu.cols() == 0) {
0509     nb_transpositions = 0;
0510     return;
0511   }
0512   eigen_assert(lu.cols() == row_transpositions.size());
0513   eigen_assert(row_transpositions.size() < 2 || (&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
0514
0515   partial_lu_impl
0516     < typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor,
0517       typename TranspositionType::StorageIndex,
0518       EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime)>
0519     ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
0520 }
0521
0522 } // end namespace internal
0523
0524 template<typename MatrixType>
0525 void PartialPivLU<MatrixType>::compute()
0526 {
0527   check_template_parameters();
0528
0529   // the row permutation is stored as int indices, so just to be sure:
0530   eigen_assert(m_lu.rows()<NumTraits<int>::highest());
0531
0532   if(m_lu.cols()>0)
0533     m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
0534   else
0535     m_l1_norm = RealScalar(0);
0536
0537   eigen_assert(m_lu.rows() == m_lu.cols() && "PartialPivLU is only for square (and moreover invertible) matrices");
0538   const Index size = m_lu.rows();
0539
0540   m_rowsTranspositions.resize(size);
0541
0542   typename TranspositionType::StorageIndex nb_transpositions;
0543   internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
0544   m_det_p = (nb_transpositions%2) ? -1 : 1;
0545
0546   m_p = m_rowsTranspositions;
0547
0548   m_isInitialized = true;
0549 }
0550
0551 template<typename MatrixType>
0552 typename PartialPivLU<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const
0553 {
0554   eigen_assert(m_isInitialized && "PartialPivLU is not initialized.");
0555   return Scalar(m_det_p) * m_lu.diagonal().prod();
0556 }
0557
0558 /** \returns the matrix represented by the decomposition,
0559  * i.e., it returns the product: P^{-1} L U.
0560  * This function is provided for debug purpose. */
0561 template<typename MatrixType>
0562 MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const
0563 {
0564   eigen_assert(m_isInitialized && "LU is not initialized.");
0565   // LU
0566   MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
0567                  * m_lu.template triangularView<Upper>();
0568
0569   // P^{-1}(LU)
0570   res = m_p.inverse() * res;
0571
0572   return res;
0573 }
0574
0575 /***** Implementation details *****************************************************/
0576
0577 namespace internal {
0578
0579 /***** Implementation of inverse() *****************************************************/
0580 template<typename DstXprType, typename MatrixType>
0581 struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename PartialPivLU<MatrixType>::Scalar>, Dense2Dense>
0582 {
0583   typedef PartialPivLU<MatrixType> LuType;
0584   typedef Inverse<LuType> SrcXprType;
0585   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename LuType::Scalar> &)
0586   {
0587     dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
0588   }
0589 };
0590 } // end namespace internal
0591
0592 /******** MatrixBase methods *******/
0593
0594 /** \lu_module
0595   *
0596   * \return the partial-pivoting LU decomposition of \c *this.
0597   *
0598   * \sa class PartialPivLU
0599   */
0600 template<typename Derived>
0601 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
0602 MatrixBase<Derived>::partialPivLu() const
0603 {
0604   return PartialPivLU<PlainObject>(eval());
0605 }
0606
0607 /** \lu_module
0608   *
0609   * Synonym of partialPivLu().
0610   *
0611   * \return the partial-pivoting LU decomposition of \c *this.
0612   *
0613   * \sa class PartialPivLU
0614   */
0615 template<typename Derived>
0616 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
0617 MatrixBase<Derived>::lu() const
0618 {
0619   return PartialPivLU<PlainObject>(eval());
0620 }
0621
0622 } // end namespace Eigen
0623
0624 #endif // EIGEN_PARTIALLU_H