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 MatrixType m_lu; 175 PermutationType m_p; 176 TranspositionType m_rowsTranspositions; 177 Index m_det_p; 178 bool m_isInitialized; 179 }; 180 181 template<typename MatrixType> 182 PartialPivLU<MatrixType>::PartialPivLU() 183 : m_lu(), 184 m_p(), 185 m_rowsTranspositions(), 186 m_det_p(0), 187 m_isInitialized(false) 188 { 189 } 190 191 template<typename MatrixType> 192 PartialPivLU<MatrixType>::PartialPivLU(Index size) 193 : m_lu(size, size), 194 m_p(size), 195 m_rowsTranspositions(size), 196 m_det_p(0), 197 m_isInitialized(false) 198 { 199 } 200 201 template<typename MatrixType> 202 PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix) 203 : m_lu(matrix.rows(), matrix.rows()), 204 m_p(matrix.rows()), 205 m_rowsTranspositions(matrix.rows()), 206 m_det_p(0), 207 m_isInitialized(false) 208 { 209 compute(matrix); 210 } 211 212 namespace internal { 213 214 /** \internal This is the blocked version of fullpivlu_unblocked() */ 215 template<typename Scalar, int StorageOrder, typename PivIndex> 216 struct partial_lu_impl 217 { 218 // FIXME add a stride to Map, so that the following mapping becomes easier, 219 // another option would be to create an expression being able to automatically 220 // warp any Map, Matrix, and Block expressions as a unique type, but since that's exactly 221 // a Map + stride, why not adding a stride to Map, and convenient ctors from a Matrix, 222 // and Block. 223 typedef Map<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > MapLU; 224 typedef Block<MapLU, Dynamic, Dynamic> MatrixType; 225 typedef Block<MatrixType,Dynamic,Dynamic> BlockType; 226 typedef typename MatrixType::RealScalar RealScalar; 227 typedef typename MatrixType::Index Index; 228 229 /** \internal performs the LU decomposition in-place of the matrix \a lu 230 * using an unblocked algorithm. 231 * 232 * In addition, this function returns the row transpositions in the 233 * vector \a row_transpositions which must have a size equal to the number 234 * of columns of the matrix \a lu, and an integer \a nb_transpositions 235 * which returns the actual number of transpositions. 236 * 237 * \returns The index of the first pivot which is exactly zero if any, or a negative number otherwise. 238 */ 239 static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions) 240 { 241 const Index rows = lu.rows(); 242 const Index cols = lu.cols(); 243 const Index size = (std::min)(rows,cols); 244 nb_transpositions = 0; 245 int first_zero_pivot = -1; 246 for(Index k = 0; k < size; ++k) 247 { 248 Index rrows = rows-k-1; 249 Index rcols = cols-k-1; 250 251 Index row_of_biggest_in_col; 252 RealScalar biggest_in_corner 253 = lu.col(k).tail(rows-k).cwiseAbs().maxCoeff(&row_of_biggest_in_col); 254 row_of_biggest_in_col += k; 255 256 row_transpositions[k] = row_of_biggest_in_col; 257 258 if(biggest_in_corner != RealScalar(0)) 259 { 260 if(k != row_of_biggest_in_col) 261 { 262 lu.row(k).swap(lu.row(row_of_biggest_in_col)); 263 ++nb_transpositions; 264 } 265 266 // FIXME shall we introduce a safe quotient expression in cas 1/lu.coeff(k,k) 267 // overflow but not the actual quotient? 268 lu.col(k).tail(rrows) /= lu.coeff(k,k); 269 } 270 else if(first_zero_pivot==-1) 271 { 272 // the pivot is exactly zero, we record the index of the first pivot which is exactly 0, 273 // and continue the factorization such we still have A = PLU 274 first_zero_pivot = k; 275 } 276 277 if(k<rows-1) 278 lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols); 279 } 280 return first_zero_pivot; 281 } 282 283 /** \internal performs the LU decomposition in-place of the matrix represented 284 * by the variables \a rows, \a cols, \a lu_data, and \a lu_stride using a 285 * recursive, blocked algorithm. 286 * 287 * In addition, this function returns the row transpositions in the 288 * vector \a row_transpositions which must have a size equal to the number 289 * of columns of the matrix \a lu, and an integer \a nb_transpositions 290 * which returns the actual number of transpositions. 291 * 292 * \returns The index of the first pivot which is exactly zero if any, or a negative number otherwise. 293 * 294 * \note This very low level interface using pointers, etc. is to: 295 * 1 - reduce the number of instanciations to the strict minimum 296 * 2 - avoid infinite recursion of the instanciations with Block<Block<Block<...> > > 297 */ 298 static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256) 299 { 300 MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols); 301 MatrixType lu(lu1,0,0,rows,cols); 302 303 const Index size = (std::min)(rows,cols); 304 305 // if the matrix is too small, no blocking: 306 if(size<=16) 307 { 308 return unblocked_lu(lu, row_transpositions, nb_transpositions); 309 } 310 311 // automatically adjust the number of subdivisions to the size 312 // of the matrix so that there is enough sub blocks: 313 Index blockSize; 314 { 315 blockSize = size/8; 316 blockSize = (blockSize/16)*16; 317 blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize); 318 } 319 320 nb_transpositions = 0; 321 int first_zero_pivot = -1; 322 for(Index k = 0; k < size; k+=blockSize) 323 { 324 Index bs = (std::min)(size-k,blockSize); // actual size of the block 325 Index trows = rows - k - bs; // trailing rows 326 Index tsize = size - k - bs; // trailing size 327 328 // partition the matrix: 329 // A00 | A01 | A02 330 // lu = A_0 | A_1 | A_2 = A10 | A11 | A12 331 // A20 | A21 | A22 332 BlockType A_0(lu,0,0,rows,k); 333 BlockType A_2(lu,0,k+bs,rows,tsize); 334 BlockType A11(lu,k,k,bs,bs); 335 BlockType A12(lu,k,k+bs,bs,tsize); 336 BlockType A21(lu,k+bs,k,trows,bs); 337 BlockType A22(lu,k+bs,k+bs,trows,tsize); 338 339 PivIndex nb_transpositions_in_panel; 340 // recursively call the blocked LU algorithm on [A11^T A21^T]^T 341 // with a very small blocking size: 342 Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride, 343 row_transpositions+k, nb_transpositions_in_panel, 16); 344 if(ret>=0 && first_zero_pivot==-1) 345 first_zero_pivot = k+ret; 346 347 nb_transpositions += nb_transpositions_in_panel; 348 // update permutations and apply them to A_0 349 for(Index i=k; i<k+bs; ++i) 350 { 351 Index piv = (row_transpositions[i] += k); 352 A_0.row(i).swap(A_0.row(piv)); 353 } 354 355 if(trows) 356 { 357 // apply permutations to A_2 358 for(Index i=k;i<k+bs; ++i) 359 A_2.row(i).swap(A_2.row(row_transpositions[i])); 360 361 // A12 = A11^-1 A12 362 A11.template triangularView<UnitLower>().solveInPlace(A12); 363 364 A22.noalias() -= A21 * A12; 365 } 366 } 367 return first_zero_pivot; 368 } 369 }; 370 371 /** \internal performs the LU decomposition with partial pivoting in-place. 372 */ 373 template<typename MatrixType, typename TranspositionType> 374 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::Index& nb_transpositions) 375 { 376 eigen_assert(lu.cols() == row_transpositions.size()); 377 eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1); 378 379 partial_lu_impl 380 <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor, typename TranspositionType::Index> 381 ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions); 382 } 383 384 } // end namespace internal 385 386 template<typename MatrixType> 387 PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(const MatrixType& matrix) 388 { 389 m_lu = matrix; 390 391 eigen_assert(matrix.rows() == matrix.cols() && "PartialPivLU is only for square (and moreover invertible) matrices"); 392 const Index size = matrix.rows(); 393 394 m_rowsTranspositions.resize(size); 395 396 typename TranspositionType::Index nb_transpositions; 397 internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions); 398 m_det_p = (nb_transpositions%2) ? -1 : 1; 399 400 m_p = m_rowsTranspositions; 401 402 m_isInitialized = true; 403 return *this; 404 } 405 406 template<typename MatrixType> 407 typename internal::traits<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const 408 { 409 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 410 return Scalar(m_det_p) * m_lu.diagonal().prod(); 411 } 412 413 /** \returns the matrix represented by the decomposition, 414 * i.e., it returns the product: P^{-1} L U. 415 * This function is provided for debug purpose. */ 416 template<typename MatrixType> 417 MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const 418 { 419 eigen_assert(m_isInitialized && "LU is not initialized."); 420 // LU 421 MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix() 422 * m_lu.template triangularView<Upper>(); 423 424 // P^{-1}(LU) 425 res = m_p.inverse() * res; 426 427 return res; 428 } 429 430 /***** Implementation of solve() *****************************************************/ 431 432 namespace internal { 433 434 template<typename _MatrixType, typename Rhs> 435 struct solve_retval<PartialPivLU<_MatrixType>, Rhs> 436 : solve_retval_base<PartialPivLU<_MatrixType>, Rhs> 437 { 438 EIGEN_MAKE_SOLVE_HELPERS(PartialPivLU<_MatrixType>,Rhs) 439 440 template<typename Dest> void evalTo(Dest& dst) const 441 { 442 /* The decomposition PA = LU can be rewritten as A = P^{-1} L U. 443 * So we proceed as follows: 444 * Step 1: compute c = Pb. 445 * Step 2: replace c by the solution x to Lx = c. 446 * Step 3: replace c by the solution x to Ux = c. 447 */ 448 449 eigen_assert(rhs().rows() == dec().matrixLU().rows()); 450 451 // Step 1 452 dst = dec().permutationP() * rhs(); 453 454 // Step 2 455 dec().matrixLU().template triangularView<UnitLower>().solveInPlace(dst); 456 457 // Step 3 458 dec().matrixLU().template triangularView<Upper>().solveInPlace(dst); 459 } 460 }; 461 462 } // end namespace internal 463 464 /******** MatrixBase methods *******/ 465 466 /** \lu_module 467 * 468 * \return the partial-pivoting LU decomposition of \c *this. 469 * 470 * \sa class PartialPivLU 471 */ 472 template<typename Derived> 473 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject> 474 MatrixBase<Derived>::partialPivLu() const 475 { 476 return PartialPivLU<PlainObject>(eval()); 477 } 478 479 #if EIGEN2_SUPPORT_STAGE > STAGE20_RESOLVE_API_CONFLICTS 480 /** \lu_module 481 * 482 * Synonym of partialPivLu(). 483 * 484 * \return the partial-pivoting LU decomposition of \c *this. 485 * 486 * \sa class PartialPivLU 487 */ 488 template<typename Derived> 489 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject> 490 MatrixBase<Derived>::lu() const 491 { 492 return PartialPivLU<PlainObject>(eval()); 493 } 494 #endif 495 496 } // end namespace Eigen 497 498 #endif // EIGEN_PARTIALLU_H 499
