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 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_PARTIALLU_H 12 #define EIGEN_PARTIALLU_H 13 14 namespace Eigen { 15 16 namespace internal { 17 template<typename _MatrixType> struct traits<PartialPivLU<_MatrixType> > 18 : traits<_MatrixType> 19 { 20 typedef MatrixXpr XprKind; 21 typedef SolverStorage StorageKind; 22 typedef int StorageIndex; 23 typedef traits<_MatrixType> BaseTraits; 24 enum { 25 Flags = BaseTraits::Flags & RowMajorBit, 26 CoeffReadCost = Dynamic 27 }; 28 }; 29 30 template<typename T,typename Derived> 31 struct enable_if_ref; 32 // { 33 // typedef Derived type; 34 // }; 35 36 template<typename T,typename Derived> 37 struct enable_if_ref<Ref<T>,Derived> { 38 typedef Derived type; 39 }; 40 41 } // end namespace internal 42 43 /** \ingroup LU_Module 44 * 45 * \class PartialPivLU 46 * 47 * \brief LU decomposition of a matrix with partial pivoting, and related features 48 * 49 * \tparam _MatrixType the type of the matrix of which we are computing the LU decomposition 50 * 51 * This class represents a LU decomposition of a \b square \b invertible matrix, with partial pivoting: the matrix A 52 * is decomposed as A = PLU where L is unit-lower-triangular, U is upper-triangular, and P 53 * is a permutation matrix. 54 * 55 * Typically, partial pivoting LU decomposition is only considered numerically stable for square invertible 56 * matrices. Thus LAPACK's dgesv and dgesvx require the matrix to be square and invertible. The present class 57 * does the same. It will assert that the matrix is square, but it won't (actually it can't) check that the 58 * matrix is invertible: it is your task to check that you only use this decomposition on invertible matrices. 59 * 60 * The guaranteed safe alternative, working for all matrices, is the full pivoting LU decomposition, provided 61 * by class FullPivLU. 62 * 63 * This is \b not a rank-revealing LU decomposition. Many features are intentionally absent from this class, 64 * such as rank computation. If you need these features, use class FullPivLU. 65 * 66 * This LU decomposition is suitable to invert invertible matrices. It is what MatrixBase::inverse() uses 67 * in the general case. 68 * On the other hand, it is \b not suitable to determine whether a given matrix is invertible. 69 * 70 * The data of the LU decomposition can be directly accessed through the methods matrixLU(), permutationP(). 71 * 72 * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. 73 * 74 * \sa MatrixBase::partialPivLu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse(), class FullPivLU 75 */ 76 template<typename _MatrixType> class PartialPivLU 77 : public SolverBase<PartialPivLU<_MatrixType> > 78 { 79 public: 80 81 typedef _MatrixType MatrixType; 82 typedef SolverBase<PartialPivLU> Base; 83 friend class SolverBase<PartialPivLU>; 84 85 EIGEN_GENERIC_PUBLIC_INTERFACE(PartialPivLU) 86 enum { 87 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 88 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 89 }; 90 typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType; 91 typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType; 92 typedef typename MatrixType::PlainObject PlainObject; 93 94 /** 95 * \brief Default Constructor. 96 * 97 * The default constructor is useful in cases in which the user intends to 98 * perform decompositions via PartialPivLU::compute(const MatrixType&). 99 */ 100 PartialPivLU(); 101 102 /** \brief Default Constructor with memory preallocation 103 * 104 * Like the default constructor but with preallocation of the internal data 105 * according to the specified problem \a size. 106 * \sa PartialPivLU() 107 */ 108 explicit PartialPivLU(Index size); 109 110 /** Constructor. 111 * 112 * \param matrix the matrix of which to compute the LU decomposition. 113 * 114 * \warning The matrix should have full rank (e.g. if it's square, it should be invertible). 115 * If you need to deal with non-full rank, use class FullPivLU instead. 116 */ 117 template<typename InputType> 118 explicit PartialPivLU(const EigenBase<InputType>& matrix); 119 120 /** Constructor for \link InplaceDecomposition inplace decomposition \endlink 121 * 122 * \param matrix the matrix of which to compute the LU decomposition. 123 * 124 * \warning The matrix should have full rank (e.g. if it's square, it should be invertible). 125 * If you need to deal with non-full rank, use class FullPivLU instead. 126 */ 127 template<typename InputType> 128 explicit PartialPivLU(EigenBase<InputType>& matrix); 129 130 template<typename InputType> 131 PartialPivLU& compute(const EigenBase<InputType>& matrix) { 132 m_lu = matrix.derived(); 133 compute(); 134 return *this; 135 } 136 137 /** \returns the LU decomposition matrix: the upper-triangular part is U, the 138 * unit-lower-triangular part is L (at least for square matrices; in the non-square 139 * case, special care is needed, see the documentation of class FullPivLU). 140 * 141 * \sa matrixL(), matrixU() 142 */ 143 inline const MatrixType& matrixLU() const 144 { 145 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 146 return m_lu; 147 } 148 149 /** \returns the permutation matrix P. 150 */ 151 inline const PermutationType& permutationP() const 152 { 153 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 154 return m_p; 155 } 156 157 #ifdef EIGEN_PARSED_BY_DOXYGEN 158 /** This method returns the solution x to the equation Ax=b, where A is the matrix of which 159 * *this is the LU decomposition. 160 * 161 * \param b the right-hand-side of the equation to solve. Can be a vector or a matrix, 162 * the only requirement in order for the equation to make sense is that 163 * b.rows()==A.rows(), where A is the matrix of which *this is the LU decomposition. 164 * 165 * \returns the solution. 166 * 167 * Example: \include PartialPivLU_solve.cpp 168 * Output: \verbinclude PartialPivLU_solve.out 169 * 170 * Since this PartialPivLU class assumes anyway that the matrix A is invertible, the solution 171 * theoretically exists and is unique regardless of b. 172 * 173 * \sa TriangularView::solve(), inverse(), computeInverse() 174 */ 175 template<typename Rhs> 176 inline const Solve<PartialPivLU, Rhs> 177 solve(const MatrixBase<Rhs>& b) const; 178 #endif 179 180 /** \returns an estimate of the reciprocal condition number of the matrix of which \c *this is 181 the LU decomposition. 182 */ 183 inline RealScalar rcond() const 184 { 185 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 186 return internal::rcond_estimate_helper(m_l1_norm, *this); 187 } 188 189 /** \returns the inverse of the matrix of which *this is the LU decomposition. 190 * 191 * \warning The matrix being decomposed here is assumed to be invertible. If you need to check for 192 * invertibility, use class FullPivLU instead. 193 * 194 * \sa MatrixBase::inverse(), LU::inverse() 195 */ 196 inline const Inverse<PartialPivLU> inverse() const 197 { 198 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 199 return Inverse<PartialPivLU>(*this); 200 } 201 202 /** \returns the determinant of the matrix of which 203 * *this is the LU decomposition. It has only linear complexity 204 * (that is, O(n) where n is the dimension of the square matrix) 205 * as the LU decomposition has already been computed. 206 * 207 * \note For fixed-size matrices of size up to 4, MatrixBase::determinant() offers 208 * optimized paths. 209 * 210 * \warning a determinant can be very big or small, so for matrices 211 * of large enough dimension, there is a risk of overflow/underflow. 212 * 213 * \sa MatrixBase::determinant() 214 */ 215 Scalar determinant() const; 216 217 MatrixType reconstructedMatrix() const; 218 219 EIGEN_CONSTEXPR inline Index rows() const EIGEN_NOEXCEPT { return m_lu.rows(); } 220 EIGEN_CONSTEXPR inline Index cols() const EIGEN_NOEXCEPT { return m_lu.cols(); } 221 222 #ifndef EIGEN_PARSED_BY_DOXYGEN 223 template<typename RhsType, typename DstType> 224 EIGEN_DEVICE_FUNC 225 void _solve_impl(const RhsType &rhs, DstType &dst) const { 226 /* The decomposition PA = LU can be rewritten as A = P^{-1} L U. 227 * So we proceed as follows: 228 * Step 1: compute c = Pb. 229 * Step 2: replace c by the solution x to Lx = c. 230 * Step 3: replace c by the solution x to Ux = c. 231 */ 232 233 // Step 1 234 dst = permutationP() * rhs; 235 236 // Step 2 237 m_lu.template triangularView<UnitLower>().solveInPlace(dst); 238 239 // Step 3 240 m_lu.template triangularView<Upper>().solveInPlace(dst); 241 } 242 243 template<bool Conjugate, typename RhsType, typename DstType> 244 EIGEN_DEVICE_FUNC 245 void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const { 246 /* The decomposition PA = LU can be rewritten as A^T = U^T L^T P. 247 * So we proceed as follows: 248 * Step 1: compute c as the solution to L^T c = b 249 * Step 2: replace c by the solution x to U^T x = c. 250 * Step 3: update c = P^-1 c. 251 */ 252 253 eigen_assert(rhs.rows() == m_lu.cols()); 254 255 // Step 1 256 dst = m_lu.template triangularView<Upper>().transpose() 257 .template conjugateIf<Conjugate>().solve(rhs); 258 // Step 2 259 m_lu.template triangularView<UnitLower>().transpose() 260 .template conjugateIf<Conjugate>().solveInPlace(dst); 261 // Step 3 262 dst = permutationP().transpose() * dst; 263 } 264 #endif 265 266 protected: 267 268 static void check_template_parameters() 269 { 270 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 271 } 272 273 void compute(); 274 275 MatrixType m_lu; 276 PermutationType m_p; 277 TranspositionType m_rowsTranspositions; 278 RealScalar m_l1_norm; 279 signed char m_det_p; 280 bool m_isInitialized; 281 }; 282 283 template<typename MatrixType> 284 PartialPivLU<MatrixType>::PartialPivLU() 285 : m_lu(), 286 m_p(), 287 m_rowsTranspositions(), 288 m_l1_norm(0), 289 m_det_p(0), 290 m_isInitialized(false) 291 { 292 } 293 294 template<typename MatrixType> 295 PartialPivLU<MatrixType>::PartialPivLU(Index size) 296 : m_lu(size, size), 297 m_p(size), 298 m_rowsTranspositions(size), 299 m_l1_norm(0), 300 m_det_p(0), 301 m_isInitialized(false) 302 { 303 } 304 305 template<typename MatrixType> 306 template<typename InputType> 307 PartialPivLU<MatrixType>::PartialPivLU(const EigenBase<InputType>& matrix) 308 : m_lu(matrix.rows(),matrix.cols()), 309 m_p(matrix.rows()), 310 m_rowsTranspositions(matrix.rows()), 311 m_l1_norm(0), 312 m_det_p(0), 313 m_isInitialized(false) 314 { 315 compute(matrix.derived()); 316 } 317 318 template<typename MatrixType> 319 template<typename InputType> 320 PartialPivLU<MatrixType>::PartialPivLU(EigenBase<InputType>& matrix) 321 : m_lu(matrix.derived()), 322 m_p(matrix.rows()), 323 m_rowsTranspositions(matrix.rows()), 324 m_l1_norm(0), 325 m_det_p(0), 326 m_isInitialized(false) 327 { 328 compute(); 329 } 330 331 namespace internal { 332 333 /** \internal This is the blocked version of fullpivlu_unblocked() */ 334 template<typename Scalar, int StorageOrder, typename PivIndex, int SizeAtCompileTime=Dynamic> 335 struct partial_lu_impl 336 { 337 static const int UnBlockedBound = 16; 338 static const bool UnBlockedAtCompileTime = SizeAtCompileTime!=Dynamic && SizeAtCompileTime<=UnBlockedBound; 339 static const int ActualSizeAtCompileTime = UnBlockedAtCompileTime ? SizeAtCompileTime : Dynamic; 340 // Remaining rows and columns at compile-time: 341 static const int RRows = SizeAtCompileTime==2 ? 1 : Dynamic; 342 static const int RCols = SizeAtCompileTime==2 ? 1 : Dynamic; 343 typedef Matrix<Scalar, ActualSizeAtCompileTime, ActualSizeAtCompileTime, StorageOrder> MatrixType; 344 typedef Ref<MatrixType> MatrixTypeRef; 345 typedef Ref<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > BlockType; 346 typedef typename MatrixType::RealScalar RealScalar; 347 348 /** \internal performs the LU decomposition in-place of the matrix \a lu 349 * using an unblocked algorithm. 350 * 351 * In addition, this function returns the row transpositions in the 352 * vector \a row_transpositions which must have a size equal to the number 353 * of columns of the matrix \a lu, and an integer \a nb_transpositions 354 * which returns the actual number of transpositions. 355 * 356 * \returns The index of the first pivot which is exactly zero if any, or a negative number otherwise. 357 */ 358 static Index unblocked_lu(MatrixTypeRef& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions) 359 { 360 typedef scalar_score_coeff_op<Scalar> Scoring; 361 typedef typename Scoring::result_type Score; 362 const Index rows = lu.rows(); 363 const Index cols = lu.cols(); 364 const Index size = (std::min)(rows,cols); 365 // For small compile-time matrices it is worth processing the last row separately: 366 // speedup: +100% for 2x2, +10% for others. 367 const Index endk = UnBlockedAtCompileTime ? size-1 : size; 368 nb_transpositions = 0; 369 Index first_zero_pivot = -1; 370 for(Index k = 0; k < endk; ++k) 371 { 372 int rrows = internal::convert_index<int>(rows-k-1); 373 int rcols = internal::convert_index<int>(cols-k-1); 374 375 Index row_of_biggest_in_col; 376 Score biggest_in_corner 377 = lu.col(k).tail(rows-k).unaryExpr(Scoring()).maxCoeff(&row_of_biggest_in_col); 378 row_of_biggest_in_col += k; 379 380 row_transpositions[k] = PivIndex(row_of_biggest_in_col); 381 382 if(biggest_in_corner != Score(0)) 383 { 384 if(k != row_of_biggest_in_col) 385 { 386 lu.row(k).swap(lu.row(row_of_biggest_in_col)); 387 ++nb_transpositions; 388 } 389 390 lu.col(k).tail(fix<RRows>(rrows)) /= lu.coeff(k,k); 391 } 392 else if(first_zero_pivot==-1) 393 { 394 // the pivot is exactly zero, we record the index of the first pivot which is exactly 0, 395 // and continue the factorization such we still have A = PLU 396 first_zero_pivot = k; 397 } 398 399 if(k<rows-1) 400 lu.bottomRightCorner(fix<RRows>(rrows),fix<RCols>(rcols)).noalias() -= lu.col(k).tail(fix<RRows>(rrows)) * lu.row(k).tail(fix<RCols>(rcols)); 401 } 402 403 // special handling of the last entry 404 if(UnBlockedAtCompileTime) 405 { 406 Index k = endk; 407 row_transpositions[k] = PivIndex(k); 408 if (Scoring()(lu(k, k)) == Score(0) && first_zero_pivot == -1) 409 first_zero_pivot = k; 410 } 411 412 return first_zero_pivot; 413 } 414 415 /** \internal performs the LU decomposition in-place of the matrix represented 416 * by the variables \a rows, \a cols, \a lu_data, and \a lu_stride using a 417 * recursive, blocked algorithm. 418 * 419 * In addition, this function returns the row transpositions in the 420 * vector \a row_transpositions which must have a size equal to the number 421 * of columns of the matrix \a lu, and an integer \a nb_transpositions 422 * which returns the actual number of transpositions. 423 * 424 * \returns The index of the first pivot which is exactly zero if any, or a negative number otherwise. 425 * 426 * \note This very low level interface using pointers, etc. is to: 427 * 1 - reduce the number of instantiations to the strict minimum 428 * 2 - avoid infinite recursion of the instantiations with Block<Block<Block<...> > > 429 */ 430 static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256) 431 { 432 MatrixTypeRef lu = MatrixType::Map(lu_data,rows, cols, OuterStride<>(luStride)); 433 434 const Index size = (std::min)(rows,cols); 435 436 // if the matrix is too small, no blocking: 437 if(UnBlockedAtCompileTime || size<=UnBlockedBound) 438 { 439 return unblocked_lu(lu, row_transpositions, nb_transpositions); 440 } 441 442 // automatically adjust the number of subdivisions to the size 443 // of the matrix so that there is enough sub blocks: 444 Index blockSize; 445 { 446 blockSize = size/8; 447 blockSize = (blockSize/16)*16; 448 blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize); 449 } 450 451 nb_transpositions = 0; 452 Index first_zero_pivot = -1; 453 for(Index k = 0; k < size; k+=blockSize) 454 { 455 Index bs = (std::min)(size-k,blockSize); // actual size of the block 456 Index trows = rows - k - bs; // trailing rows 457 Index tsize = size - k - bs; // trailing size 458 459 // partition the matrix: 460 // A00 | A01 | A02 461 // lu = A_0 | A_1 | A_2 = A10 | A11 | A12 462 // A20 | A21 | A22 463 BlockType A_0 = lu.block(0,0,rows,k); 464 BlockType A_2 = lu.block(0,k+bs,rows,tsize); 465 BlockType A11 = lu.block(k,k,bs,bs); 466 BlockType A12 = lu.block(k,k+bs,bs,tsize); 467 BlockType A21 = lu.block(k+bs,k,trows,bs); 468 BlockType A22 = lu.block(k+bs,k+bs,trows,tsize); 469 470 PivIndex nb_transpositions_in_panel; 471 // recursively call the blocked LU algorithm on [A11^T A21^T]^T 472 // with a very small blocking size: 473 Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride, 474 row_transpositions+k, nb_transpositions_in_panel, 16); 475 if(ret>=0 && first_zero_pivot==-1) 476 first_zero_pivot = k+ret; 477 478 nb_transpositions += nb_transpositions_in_panel; 479 // update permutations and apply them to A_0 480 for(Index i=k; i<k+bs; ++i) 481 { 482 Index piv = (row_transpositions[i] += internal::convert_index<PivIndex>(k)); 483 A_0.row(i).swap(A_0.row(piv)); 484 } 485 486 if(trows) 487 { 488 // apply permutations to A_2 489 for(Index i=k;i<k+bs; ++i) 490 A_2.row(i).swap(A_2.row(row_transpositions[i])); 491 492 // A12 = A11^-1 A12 493 A11.template triangularView<UnitLower>().solveInPlace(A12); 494 495 A22.noalias() -= A21 * A12; 496 } 497 } 498 return first_zero_pivot; 499 } 500 }; 501 502 /** \internal performs the LU decomposition with partial pivoting in-place. 503 */ 504 template<typename MatrixType, typename TranspositionType> 505 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::StorageIndex& nb_transpositions) 506 { 507 // Special-case of zero matrix. 508 if (lu.rows() == 0 || lu.cols() == 0) { 509 nb_transpositions = 0; 510 return; 511 } 512 eigen_assert(lu.cols() == row_transpositions.size()); 513 eigen_assert(row_transpositions.size() < 2 || (&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1); 514 515 partial_lu_impl 516 < typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor, 517 typename TranspositionType::StorageIndex, 518 EIGEN_SIZE_MIN_PREFER_FIXED(MatrixType::RowsAtCompileTime,MatrixType::ColsAtCompileTime)> 519 ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions); 520 } 521 522 } // end namespace internal 523 524 template<typename MatrixType> 525 void PartialPivLU<MatrixType>::compute() 526 { 527 check_template_parameters(); 528 529 // the row permutation is stored as int indices, so just to be sure: 530 eigen_assert(m_lu.rows()<NumTraits<int>::highest()); 531 532 if(m_lu.cols()>0) 533 m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff(); 534 else 535 m_l1_norm = RealScalar(0); 536 537 eigen_assert(m_lu.rows() == m_lu.cols() && "PartialPivLU is only for square (and moreover invertible) matrices"); 538 const Index size = m_lu.rows(); 539 540 m_rowsTranspositions.resize(size); 541 542 typename TranspositionType::StorageIndex nb_transpositions; 543 internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions); 544 m_det_p = (nb_transpositions%2) ? -1 : 1; 545 546 m_p = m_rowsTranspositions; 547 548 m_isInitialized = true; 549 } 550 551 template<typename MatrixType> 552 typename PartialPivLU<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const 553 { 554 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 555 return Scalar(m_det_p) * m_lu.diagonal().prod(); 556 } 557 558 /** \returns the matrix represented by the decomposition, 559 * i.e., it returns the product: P^{-1} L U. 560 * This function is provided for debug purpose. */ 561 template<typename MatrixType> 562 MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const 563 { 564 eigen_assert(m_isInitialized && "LU is not initialized."); 565 // LU 566 MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix() 567 * m_lu.template triangularView<Upper>(); 568 569 // P^{-1}(LU) 570 res = m_p.inverse() * res; 571 572 return res; 573 } 574 575 /***** Implementation details *****************************************************/ 576 577 namespace internal { 578 579 /***** Implementation of inverse() *****************************************************/ 580 template<typename DstXprType, typename MatrixType> 581 struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename PartialPivLU<MatrixType>::Scalar>, Dense2Dense> 582 { 583 typedef PartialPivLU<MatrixType> LuType; 584 typedef Inverse<LuType> SrcXprType; 585 static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename LuType::Scalar> &) 586 { 587 dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols())); 588 } 589 }; 590 } // end namespace internal 591 592 /******** MatrixBase methods *******/ 593 594 /** \lu_module 595 * 596 * \return the partial-pivoting LU decomposition of \c *this. 597 * 598 * \sa class PartialPivLU 599 */ 600 template<typename Derived> 601 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject> 602 MatrixBase<Derived>::partialPivLu() const 603 { 604 return PartialPivLU<PlainObject>(eval()); 605 } 606 607 /** \lu_module 608 * 609 * Synonym of partialPivLu(). 610 * 611 * \return the partial-pivoting LU decomposition of \c *this. 612 * 613 * \sa class PartialPivLU 614 */ 615 template<typename Derived> 616 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject> 617 MatrixBase<Derived>::lu() const 618 { 619 return PartialPivLU<PlainObject>(eval()); 620 } 621 622 } // end namespace Eigen 623 624 #endif // EIGEN_PARTIALLU_H 625