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 (at) gmail.com> 5 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud (at) 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 /** \ingroup LU_Module 17 * 18 * \class PartialPivLU 19 * 20 * \brief LU decomposition of a matrix with partial pivoting, and related features 21 * 22 * \param MatrixType the type of the matrix of which we are computing the LU decomposition 23 * 24 * This class represents a LU decomposition of a \b square \b invertible matrix, with partial pivoting: the matrix A 25 * is decomposed as A = PLU where L is unit-lower-triangular, U is upper-triangular, and P 26 * is a permutation matrix. 27 * 28 * Typically, partial pivoting LU decomposition is only considered numerically stable for square invertible 29 * matrices. Thus LAPACK's dgesv and dgesvx require the matrix to be square and invertible. The present class 30 * does the same. It will assert that the matrix is square, but it won't (actually it can't) check that the 31 * matrix is invertible: it is your task to check that you only use this decomposition on invertible matrices. 32 * 33 * The guaranteed safe alternative, working for all matrices, is the full pivoting LU decomposition, provided 34 * by class FullPivLU. 35 * 36 * This is \b not a rank-revealing LU decomposition. Many features are intentionally absent from this class, 37 * such as rank computation. If you need these features, use class FullPivLU. 38 * 39 * This LU decomposition is suitable to invert invertible matrices. It is what MatrixBase::inverse() uses 40 * in the general case. 41 * On the other hand, it is \b not suitable to determine whether a given matrix is invertible. 42 * 43 * The data of the LU decomposition can be directly accessed through the methods matrixLU(), permutationP(). 44 * 45 * \sa MatrixBase::partialPivLu(), MatrixBase::determinant(), MatrixBase::inverse(), MatrixBase::computeInverse(), class FullPivLU 46 */ 47 template<typename _MatrixType> class PartialPivLU 48 { 49 public: 50 51 typedef _MatrixType MatrixType; 52 enum { 53 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 54 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 55 Options = MatrixType::Options, 56 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 57 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 58 }; 59 typedef typename MatrixType::Scalar Scalar; 60 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar; 61 typedef typename internal::traits<MatrixType>::StorageKind StorageKind; 62 typedef typename MatrixType::Index Index; 63 typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType; 64 typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType; 65 66 67 /** 68 * \brief Default Constructor. 69 * 70 * The default constructor is useful in cases in which the user intends to 71 * perform decompositions via PartialPivLU::compute(const MatrixType&). 72 */ 73 PartialPivLU(); 74 75 /** \brief Default Constructor with memory preallocation 76 * 77 * Like the default constructor but with preallocation of the internal data 78 * according to the specified problem \a size. 79 * \sa PartialPivLU() 80 */ 81 PartialPivLU(Index size); 82 83 /** Constructor. 84 * 85 * \param matrix the matrix of which to compute the LU decomposition. 86 * 87 * \warning The matrix should have full rank (e.g. if it's square, it should be invertible). 88 * If you need to deal with non-full rank, use class FullPivLU instead. 89 */ 90 PartialPivLU(const MatrixType& matrix); 91 92 PartialPivLU& compute(const MatrixType& matrix); 93 94 /** \returns the LU decomposition matrix: the upper-triangular part is U, the 95 * unit-lower-triangular part is L (at least for square matrices; in the non-square 96 * case, special care is needed, see the documentation of class FullPivLU). 97 * 98 * \sa matrixL(), matrixU() 99 */ 100 inline const MatrixType& matrixLU() const 101 { 102 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 103 return m_lu; 104 } 105 106 /** \returns the permutation matrix P. 107 */ 108 inline const PermutationType& permutationP() const 109 { 110 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 111 return m_p; 112 } 113 114 /** This method returns the solution x to the equation Ax=b, where A is the matrix of which 115 * *this is the LU decomposition. 116 * 117 * \param b the right-hand-side of the equation to solve. Can be a vector or a matrix, 118 * the only requirement in order for the equation to make sense is that 119 * b.rows()==A.rows(), where A is the matrix of which *this is the LU decomposition. 120 * 121 * \returns the solution. 122 * 123 * Example: \include PartialPivLU_solve.cpp 124 * Output: \verbinclude PartialPivLU_solve.out 125 * 126 * Since this PartialPivLU class assumes anyway that the matrix A is invertible, the solution 127 * theoretically exists and is unique regardless of b. 128 * 129 * \sa TriangularView::solve(), inverse(), computeInverse() 130 */ 131 template<typename Rhs> 132 inline const internal::solve_retval<PartialPivLU, Rhs> 133 solve(const MatrixBase<Rhs>& b) const 134 { 135 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 136 return internal::solve_retval<PartialPivLU, Rhs>(*this, b.derived()); 137 } 138 139 /** \returns the inverse of the matrix of which *this is the LU decomposition. 140 * 141 * \warning The matrix being decomposed here is assumed to be invertible. If you need to check for 142 * invertibility, use class FullPivLU instead. 143 * 144 * \sa MatrixBase::inverse(), LU::inverse() 145 */ 146 inline const internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType> inverse() const 147 { 148 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 149 return internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType> 150 (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols())); 151 } 152 153 /** \returns the determinant of the matrix of which 154 * *this is the LU decomposition. It has only linear complexity 155 * (that is, O(n) where n is the dimension of the square matrix) 156 * as the LU decomposition has already been computed. 157 * 158 * \note For fixed-size matrices of size up to 4, MatrixBase::determinant() offers 159 * optimized paths. 160 * 161 * \warning a determinant can be very big or small, so for matrices 162 * of large enough dimension, there is a risk of overflow/underflow. 163 * 164 * \sa MatrixBase::determinant() 165 */ 166 typename internal::traits<MatrixType>::Scalar determinant() const; 167 168 MatrixType reconstructedMatrix() const; 169 170 inline Index rows() const { return m_lu.rows(); } 171 inline Index cols() const { return m_lu.cols(); } 172 173 protected: 174 175 static void check_template_parameters() 176 { 177 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 178 } 179 180 MatrixType m_lu; 181 PermutationType m_p; 182 TranspositionType m_rowsTranspositions; 183 Index m_det_p; 184 bool m_isInitialized; 185 }; 186 187 template<typename MatrixType> 188 PartialPivLU<MatrixType>::PartialPivLU() 189 : m_lu(), 190 m_p(), 191 m_rowsTranspositions(), 192 m_det_p(0), 193 m_isInitialized(false) 194 { 195 } 196 197 template<typename MatrixType> 198 PartialPivLU<MatrixType>::PartialPivLU(Index size) 199 : m_lu(size, size), 200 m_p(size), 201 m_rowsTranspositions(size), 202 m_det_p(0), 203 m_isInitialized(false) 204 { 205 } 206 207 template<typename MatrixType> 208 PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix) 209 : m_lu(matrix.rows(), matrix.rows()), 210 m_p(matrix.rows()), 211 m_rowsTranspositions(matrix.rows()), 212 m_det_p(0), 213 m_isInitialized(false) 214 { 215 compute(matrix); 216 } 217 218 namespace internal { 219 220 /** \internal This is the blocked version of fullpivlu_unblocked() */ 221 template<typename Scalar, int StorageOrder, typename PivIndex> 222 struct partial_lu_impl 223 { 224 // FIXME add a stride to Map, so that the following mapping becomes easier, 225 // another option would be to create an expression being able to automatically 226 // warp any Map, Matrix, and Block expressions as a unique type, but since that's exactly 227 // a Map + stride, why not adding a stride to Map, and convenient ctors from a Matrix, 228 // and Block. 229 typedef Map<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > MapLU; 230 typedef Block<MapLU, Dynamic, Dynamic> MatrixType; 231 typedef Block<MatrixType,Dynamic,Dynamic> BlockType; 232 typedef typename MatrixType::RealScalar RealScalar; 233 typedef typename MatrixType::Index Index; 234 235 /** \internal performs the LU decomposition in-place of the matrix \a lu 236 * using an unblocked algorithm. 237 * 238 * In addition, this function returns the row transpositions in the 239 * vector \a row_transpositions which must have a size equal to the number 240 * of columns of the matrix \a lu, and an integer \a nb_transpositions 241 * which returns the actual number of transpositions. 242 * 243 * \returns The index of the first pivot which is exactly zero if any, or a negative number otherwise. 244 */ 245 static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions) 246 { 247 const Index rows = lu.rows(); 248 const Index cols = lu.cols(); 249 const Index size = (std::min)(rows,cols); 250 nb_transpositions = 0; 251 Index first_zero_pivot = -1; 252 for(Index k = 0; k < size; ++k) 253 { 254 Index rrows = rows-k-1; 255 Index rcols = cols-k-1; 256 257 Index row_of_biggest_in_col; 258 RealScalar biggest_in_corner 259 = lu.col(k).tail(rows-k).cwiseAbs().maxCoeff(&row_of_biggest_in_col); 260 row_of_biggest_in_col += k; 261 262 row_transpositions[k] = PivIndex(row_of_biggest_in_col); 263 264 if(biggest_in_corner != RealScalar(0)) 265 { 266 if(k != row_of_biggest_in_col) 267 { 268 lu.row(k).swap(lu.row(row_of_biggest_in_col)); 269 ++nb_transpositions; 270 } 271 272 // FIXME shall we introduce a safe quotient expression in cas 1/lu.coeff(k,k) 273 // overflow but not the actual quotient? 274 lu.col(k).tail(rrows) /= lu.coeff(k,k); 275 } 276 else if(first_zero_pivot==-1) 277 { 278 // the pivot is exactly zero, we record the index of the first pivot which is exactly 0, 279 // and continue the factorization such we still have A = PLU 280 first_zero_pivot = k; 281 } 282 283 if(k<rows-1) 284 lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols); 285 } 286 return first_zero_pivot; 287 } 288 289 /** \internal performs the LU decomposition in-place of the matrix represented 290 * by the variables \a rows, \a cols, \a lu_data, and \a lu_stride using a 291 * recursive, blocked algorithm. 292 * 293 * In addition, this function returns the row transpositions in the 294 * vector \a row_transpositions which must have a size equal to the number 295 * of columns of the matrix \a lu, and an integer \a nb_transpositions 296 * which returns the actual number of transpositions. 297 * 298 * \returns The index of the first pivot which is exactly zero if any, or a negative number otherwise. 299 * 300 * \note This very low level interface using pointers, etc. is to: 301 * 1 - reduce the number of instanciations to the strict minimum 302 * 2 - avoid infinite recursion of the instanciations with Block<Block<Block<...> > > 303 */ 304 static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256) 305 { 306 MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols); 307 MatrixType lu(lu1,0,0,rows,cols); 308 309 const Index size = (std::min)(rows,cols); 310 311 // if the matrix is too small, no blocking: 312 if(size<=16) 313 { 314 return unblocked_lu(lu, row_transpositions, nb_transpositions); 315 } 316 317 // automatically adjust the number of subdivisions to the size 318 // of the matrix so that there is enough sub blocks: 319 Index blockSize; 320 { 321 blockSize = size/8; 322 blockSize = (blockSize/16)*16; 323 blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize); 324 } 325 326 nb_transpositions = 0; 327 Index first_zero_pivot = -1; 328 for(Index k = 0; k < size; k+=blockSize) 329 { 330 Index bs = (std::min)(size-k,blockSize); // actual size of the block 331 Index trows = rows - k - bs; // trailing rows 332 Index tsize = size - k - bs; // trailing size 333 334 // partition the matrix: 335 // A00 | A01 | A02 336 // lu = A_0 | A_1 | A_2 = A10 | A11 | A12 337 // A20 | A21 | A22 338 BlockType A_0(lu,0,0,rows,k); 339 BlockType A_2(lu,0,k+bs,rows,tsize); 340 BlockType A11(lu,k,k,bs,bs); 341 BlockType A12(lu,k,k+bs,bs,tsize); 342 BlockType A21(lu,k+bs,k,trows,bs); 343 BlockType A22(lu,k+bs,k+bs,trows,tsize); 344 345 PivIndex nb_transpositions_in_panel; 346 // recursively call the blocked LU algorithm on [A11^T A21^T]^T 347 // with a very small blocking size: 348 Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride, 349 row_transpositions+k, nb_transpositions_in_panel, 16); 350 if(ret>=0 && first_zero_pivot==-1) 351 first_zero_pivot = k+ret; 352 353 nb_transpositions += nb_transpositions_in_panel; 354 // update permutations and apply them to A_0 355 for(Index i=k; i<k+bs; ++i) 356 { 357 Index piv = (row_transpositions[i] += k); 358 A_0.row(i).swap(A_0.row(piv)); 359 } 360 361 if(trows) 362 { 363 // apply permutations to A_2 364 for(Index i=k;i<k+bs; ++i) 365 A_2.row(i).swap(A_2.row(row_transpositions[i])); 366 367 // A12 = A11^-1 A12 368 A11.template triangularView<UnitLower>().solveInPlace(A12); 369 370 A22.noalias() -= A21 * A12; 371 } 372 } 373 return first_zero_pivot; 374 } 375 }; 376 377 /** \internal performs the LU decomposition with partial pivoting in-place. 378 */ 379 template<typename MatrixType, typename TranspositionType> 380 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::Index& nb_transpositions) 381 { 382 eigen_assert(lu.cols() == row_transpositions.size()); 383 eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1); 384 385 partial_lu_impl 386 <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor, typename TranspositionType::Index> 387 ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions); 388 } 389 390 } // end namespace internal 391 392 template<typename MatrixType> 393 PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& matrix) 394 { 395 check_template_parameters(); 396 397 // the row permutation is stored as int indices, so just to be sure: 398 eigen_assert(matrix.rows()<NumTraits<int>::highest()); 399 400 m_lu = matrix; 401 402 eigen_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices"); 403 const Index size = matrix.rows(); 404 405 m_rowsTranspositions.resize(size); 406 407 typename TranspositionType::Index nb_transpositions; 408 internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions); 409 m_det_p = (nb_transpositions%2) ? -1 : 1; 410 411 m_p = m_rowsTranspositions; 412 413 m_isInitialized = true; 414 return *this; 415 } 416 417 template<typename MatrixType> 418 typename internal::traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const 419 { 420 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 421 return Scalar(m_det_p) * m_lu.diagonal().prod(); 422 } 423 424 /** \returns the matrix represented by the decomposition, 425 * i.e., it returns the product: P^{-1} L U. 426 * This function is provided for debug purpose. */ 427 template<typename MatrixType> 428 MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const 429 { 430 eigen_assert(m_isInitialized && "LU is not initialized."); 431 // LU 432 MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix() 433 * m_lu.template triangularView<Upper>(); 434 435 // P^{-1}(LU) 436 res = m_p.inverse() * res; 437 438 return res; 439 } 440 441 /***** Implementation of solve() *****************************************************/ 442 443 namespace internal { 444 445 template<typename _MatrixType, typename Rhs> 446 struct solve_retval<PartialPivLU<_MatrixType>, Rhs> 447 : solve_retval_base<PartialPivLU<_MatrixType>, Rhs> 448 { 449 EIGEN_MAKE_SOLVE_HELPERS(PartialPivLU<_MatrixType>,Rhs) 450 451 template<typename Dest> void evalTo(Dest& dst) const 452 { 453 /* The decomposition PA = LU can be rewritten as A = P^{-1} L U. 454 * So we proceed as follows: 455 * Step 1: compute c = Pb. 456 * Step 2: replace c by the solution x to Lx = c. 457 * Step 3: replace c by the solution x to Ux = c. 458 */ 459 460 eigen_assert(rhs().rows() == dec().matrixLU().rows()); 461 462 // Step 1 463 dst = dec().permutationP() * rhs(); 464 465 // Step 2 466 dec().matrixLU().template triangularView<UnitLower>().solveInPlace(dst); 467 468 // Step 3 469 dec().matrixLU().template triangularView<Upper>().solveInPlace(dst); 470 } 471 }; 472 473 } // end namespace internal 474 475 /******** MatrixBase methods *******/ 476 477 /** \lu_module 478 * 479 * \return the partial-pivoting LU decomposition of \c *this. 480 * 481 * \sa class PartialPivLU 482 */ 483 template<typename Derived> 484 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject> 485 MatrixBase<Derived>::partialPivLu() const 486 { 487 return PartialPivLU<PlainObject>(eval()); 488 } 489 490 #if EIGEN2_SUPPORT_STAGE > STAGE20_RESOLVE_API_CONFLICTS 491 /** \lu_module 492 * 493 * Synonym of partialPivLu(). 494 * 495 * \return the partial-pivoting LU decomposition of \c *this. 496 * 497 * \sa class PartialPivLU 498 */ 499 template<typename Derived> 500 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject> 501 MatrixBase<Derived>::lu() const 502 { 503 return PartialPivLU<PlainObject>(eval()); 504 } 505 #endif 506 507 } // end namespace Eigen 508 509 #endif // EIGEN_PARTIALLU_H 510
