![]() |
Eigen
3.3.3
|
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 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_PARTIALLU_H 00012 #define EIGEN_PARTIALLU_H 00013 00014 namespace Eigen { 00015 00016 namespace internal { 00017 template<typename _MatrixType> struct traits<PartialPivLU<_MatrixType> > 00018 : traits<_MatrixType> 00019 { 00020 typedef MatrixXpr XprKind; 00021 typedef SolverStorage StorageKind; 00022 typedef traits<_MatrixType> BaseTraits; 00023 enum { 00024 Flags = BaseTraits::Flags & RowMajorBit, 00025 CoeffReadCost = Dynamic 00026 }; 00027 }; 00028 00029 template<typename T,typename Derived> 00030 struct enable_if_ref; 00031 // { 00032 // typedef Derived type; 00033 // }; 00034 00035 template<typename T,typename Derived> 00036 struct enable_if_ref<Ref<T>,Derived> { 00037 typedef Derived type; 00038 }; 00039 00040 } // end namespace internal 00041 00075 template<typename _MatrixType> class PartialPivLU 00076 : public SolverBase<PartialPivLU<_MatrixType> > 00077 { 00078 public: 00079 00080 typedef _MatrixType MatrixType; 00081 typedef SolverBase<PartialPivLU> Base; 00082 EIGEN_GENERIC_PUBLIC_INTERFACE(PartialPivLU) 00083 // FIXME StorageIndex defined in EIGEN_GENERIC_PUBLIC_INTERFACE should be int 00084 enum { 00085 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00086 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00087 }; 00088 typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType; 00089 typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType; 00090 typedef typename MatrixType::PlainObject PlainObject; 00091 00098 PartialPivLU(); 00099 00106 explicit PartialPivLU(Index size); 00107 00115 template<typename InputType> 00116 explicit PartialPivLU(const EigenBase<InputType>& matrix); 00117 00125 template<typename InputType> 00126 explicit PartialPivLU(EigenBase<InputType>& matrix); 00127 00128 template<typename InputType> 00129 PartialPivLU& compute(const EigenBase<InputType>& matrix) { 00130 m_lu = matrix.derived(); 00131 compute(); 00132 return *this; 00133 } 00134 00141 inline const MatrixType& matrixLU() const 00142 { 00143 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 00144 return m_lu; 00145 } 00146 00149 inline const PermutationType& permutationP() const 00150 { 00151 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 00152 return m_p; 00153 } 00154 00172 // FIXME this is a copy-paste of the base-class member to add the isInitialized assertion. 00173 template<typename Rhs> 00174 inline const Solve<PartialPivLU, Rhs> 00175 solve(const MatrixBase<Rhs>& b) const 00176 { 00177 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 00178 return Solve<PartialPivLU, Rhs>(*this, b.derived()); 00179 } 00180 00184 inline RealScalar rcond() const 00185 { 00186 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 00187 return internal::rcond_estimate_helper(m_l1_norm, *this); 00188 } 00189 00197 inline const Inverse<PartialPivLU> inverse() const 00198 { 00199 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 00200 return Inverse<PartialPivLU>(*this); 00201 } 00202 00216 Scalar determinant() const; 00217 00218 MatrixType reconstructedMatrix() const; 00219 00220 inline Index rows() const { return m_lu.rows(); } 00221 inline Index cols() const { return m_lu.cols(); } 00222 00223 #ifndef EIGEN_PARSED_BY_DOXYGEN 00224 template<typename RhsType, typename DstType> 00225 EIGEN_DEVICE_FUNC 00226 void _solve_impl(const RhsType &rhs, DstType &dst) const { 00227 /* The decomposition PA = LU can be rewritten as A = P^{-1} L U. 00228 * So we proceed as follows: 00229 * Step 1: compute c = Pb. 00230 * Step 2: replace c by the solution x to Lx = c. 00231 * Step 3: replace c by the solution x to Ux = c. 00232 */ 00233 00234 eigen_assert(rhs.rows() == m_lu.rows()); 00235 00236 // Step 1 00237 dst = permutationP() * rhs; 00238 00239 // Step 2 00240 m_lu.template triangularView<UnitLower>().solveInPlace(dst); 00241 00242 // Step 3 00243 m_lu.template triangularView<Upper>().solveInPlace(dst); 00244 } 00245 00246 template<bool Conjugate, typename RhsType, typename DstType> 00247 EIGEN_DEVICE_FUNC 00248 void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const { 00249 /* The decomposition PA = LU can be rewritten as A = P^{-1} L U. 00250 * So we proceed as follows: 00251 * Step 1: compute c = Pb. 00252 * Step 2: replace c by the solution x to Lx = c. 00253 * Step 3: replace c by the solution x to Ux = c. 00254 */ 00255 00256 eigen_assert(rhs.rows() == m_lu.cols()); 00257 00258 if (Conjugate) { 00259 // Step 1 00260 dst = m_lu.template triangularView<Upper>().adjoint().solve(rhs); 00261 // Step 2 00262 m_lu.template triangularView<UnitLower>().adjoint().solveInPlace(dst); 00263 } else { 00264 // Step 1 00265 dst = m_lu.template triangularView<Upper>().transpose().solve(rhs); 00266 // Step 2 00267 m_lu.template triangularView<UnitLower>().transpose().solveInPlace(dst); 00268 } 00269 // Step 3 00270 dst = permutationP().transpose() * dst; 00271 } 00272 #endif 00273 00274 protected: 00275 00276 static void check_template_parameters() 00277 { 00278 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00279 } 00280 00281 void compute(); 00282 00283 MatrixType m_lu; 00284 PermutationType m_p; 00285 TranspositionType m_rowsTranspositions; 00286 RealScalar m_l1_norm; 00287 signed char m_det_p; 00288 bool m_isInitialized; 00289 }; 00290 00291 template<typename MatrixType> 00292 PartialPivLU<MatrixType>::PartialPivLU() 00293 : m_lu(), 00294 m_p(), 00295 m_rowsTranspositions(), 00296 m_l1_norm(0), 00297 m_det_p(0), 00298 m_isInitialized(false) 00299 { 00300 } 00301 00302 template<typename MatrixType> 00303 PartialPivLU<MatrixType>::PartialPivLU(Index size) 00304 : m_lu(size, size), 00305 m_p(size), 00306 m_rowsTranspositions(size), 00307 m_l1_norm(0), 00308 m_det_p(0), 00309 m_isInitialized(false) 00310 { 00311 } 00312 00313 template<typename MatrixType> 00314 template<typename InputType> 00315 PartialPivLU<MatrixType>::PartialPivLU(const EigenBase<InputType>& matrix) 00316 : m_lu(matrix.rows(),matrix.cols()), 00317 m_p(matrix.rows()), 00318 m_rowsTranspositions(matrix.rows()), 00319 m_l1_norm(0), 00320 m_det_p(0), 00321 m_isInitialized(false) 00322 { 00323 compute(matrix.derived()); 00324 } 00325 00326 template<typename MatrixType> 00327 template<typename InputType> 00328 PartialPivLU<MatrixType>::PartialPivLU(EigenBase<InputType>& matrix) 00329 : m_lu(matrix.derived()), 00330 m_p(matrix.rows()), 00331 m_rowsTranspositions(matrix.rows()), 00332 m_l1_norm(0), 00333 m_det_p(0), 00334 m_isInitialized(false) 00335 { 00336 compute(); 00337 } 00338 00339 namespace internal { 00340 00342 template<typename Scalar, int StorageOrder, typename PivIndex> 00343 struct partial_lu_impl 00344 { 00345 // FIXME add a stride to Map, so that the following mapping becomes easier, 00346 // another option would be to create an expression being able to automatically 00347 // warp any Map, Matrix, and Block expressions as a unique type, but since that's exactly 00348 // a Map + stride, why not adding a stride to Map, and convenient ctors from a Matrix, 00349 // and Block. 00350 typedef Map<Matrix<Scalar, Dynamic, Dynamic, StorageOrder> > MapLU; 00351 typedef Block<MapLU, Dynamic, Dynamic> MatrixType; 00352 typedef Block<MatrixType,Dynamic,Dynamic> BlockType; 00353 typedef typename MatrixType::RealScalar RealScalar; 00354 00365 static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions) 00366 { 00367 typedef scalar_score_coeff_op<Scalar> Scoring; 00368 typedef typename Scoring::result_type Score; 00369 const Index rows = lu.rows(); 00370 const Index cols = lu.cols(); 00371 const Index size = (std::min)(rows,cols); 00372 nb_transpositions = 0; 00373 Index first_zero_pivot = -1; 00374 for(Index k = 0; k < size; ++k) 00375 { 00376 Index rrows = rows-k-1; 00377 Index rcols = cols-k-1; 00378 00379 Index row_of_biggest_in_col; 00380 Score biggest_in_corner 00381 = lu.col(k).tail(rows-k).unaryExpr(Scoring()).maxCoeff(&row_of_biggest_in_col); 00382 row_of_biggest_in_col += k; 00383 00384 row_transpositions[k] = PivIndex(row_of_biggest_in_col); 00385 00386 if(biggest_in_corner != Score(0)) 00387 { 00388 if(k != row_of_biggest_in_col) 00389 { 00390 lu.row(k).swap(lu.row(row_of_biggest_in_col)); 00391 ++nb_transpositions; 00392 } 00393 00394 // FIXME shall we introduce a safe quotient expression in cas 1/lu.coeff(k,k) 00395 // overflow but not the actual quotient? 00396 lu.col(k).tail(rrows) /= lu.coeff(k,k); 00397 } 00398 else if(first_zero_pivot==-1) 00399 { 00400 // the pivot is exactly zero, we record the index of the first pivot which is exactly 0, 00401 // and continue the factorization such we still have A = PLU 00402 first_zero_pivot = k; 00403 } 00404 00405 if(k<rows-1) 00406 lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols); 00407 } 00408 return first_zero_pivot; 00409 } 00410 00426 static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256) 00427 { 00428 MapLU lu1(lu_data,StorageOrder==RowMajor?rows:luStride,StorageOrder==RowMajor?luStride:cols); 00429 MatrixType lu(lu1,0,0,rows,cols); 00430 00431 const Index size = (std::min)(rows,cols); 00432 00433 // if the matrix is too small, no blocking: 00434 if(size<=16) 00435 { 00436 return unblocked_lu(lu, row_transpositions, nb_transpositions); 00437 } 00438 00439 // automatically adjust the number of subdivisions to the size 00440 // of the matrix so that there is enough sub blocks: 00441 Index blockSize; 00442 { 00443 blockSize = size/8; 00444 blockSize = (blockSize/16)*16; 00445 blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize); 00446 } 00447 00448 nb_transpositions = 0; 00449 Index first_zero_pivot = -1; 00450 for(Index k = 0; k < size; k+=blockSize) 00451 { 00452 Index bs = (std::min)(size-k,blockSize); // actual size of the block 00453 Index trows = rows - k - bs; // trailing rows 00454 Index tsize = size - k - bs; // trailing size 00455 00456 // partition the matrix: 00457 // A00 | A01 | A02 00458 // lu = A_0 | A_1 | A_2 = A10 | A11 | A12 00459 // A20 | A21 | A22 00460 BlockType A_0(lu,0,0,rows,k); 00461 BlockType A_2(lu,0,k+bs,rows,tsize); 00462 BlockType A11(lu,k,k,bs,bs); 00463 BlockType A12(lu,k,k+bs,bs,tsize); 00464 BlockType A21(lu,k+bs,k,trows,bs); 00465 BlockType A22(lu,k+bs,k+bs,trows,tsize); 00466 00467 PivIndex nb_transpositions_in_panel; 00468 // recursively call the blocked LU algorithm on [A11^T A21^T]^T 00469 // with a very small blocking size: 00470 Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride, 00471 row_transpositions+k, nb_transpositions_in_panel, 16); 00472 if(ret>=0 && first_zero_pivot==-1) 00473 first_zero_pivot = k+ret; 00474 00475 nb_transpositions += nb_transpositions_in_panel; 00476 // update permutations and apply them to A_0 00477 for(Index i=k; i<k+bs; ++i) 00478 { 00479 Index piv = (row_transpositions[i] += internal::convert_index<PivIndex>(k)); 00480 A_0.row(i).swap(A_0.row(piv)); 00481 } 00482 00483 if(trows) 00484 { 00485 // apply permutations to A_2 00486 for(Index i=k;i<k+bs; ++i) 00487 A_2.row(i).swap(A_2.row(row_transpositions[i])); 00488 00489 // A12 = A11^-1 A12 00490 A11.template triangularView<UnitLower>().solveInPlace(A12); 00491 00492 A22.noalias() -= A21 * A12; 00493 } 00494 } 00495 return first_zero_pivot; 00496 } 00497 }; 00498 00501 template<typename MatrixType, typename TranspositionType> 00502 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions, typename TranspositionType::StorageIndex& nb_transpositions) 00503 { 00504 eigen_assert(lu.cols() == row_transpositions.size()); 00505 eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1); 00506 00507 partial_lu_impl 00508 <typename MatrixType::Scalar, MatrixType::Flags&RowMajorBit?RowMajor:ColMajor, typename TranspositionType::StorageIndex> 00509 ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions); 00510 } 00511 00512 } // end namespace internal 00513 00514 template<typename MatrixType> 00515 void PartialPivLU<MatrixType>::compute() 00516 { 00517 check_template_parameters(); 00518 00519 // the row permutation is stored as int indices, so just to be sure: 00520 eigen_assert(m_lu.rows()<NumTraits<int>::highest()); 00521 00522 m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff(); 00523 00524 eigen_assert(m_lu.rows() == m_lu.cols() && "PartialPivLU is only for square (and moreover invertible) matrices"); 00525 const Index size = m_lu.rows(); 00526 00527 m_rowsTranspositions.resize(size); 00528 00529 typename TranspositionType::StorageIndex nb_transpositions; 00530 internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions); 00531 m_det_p = (nb_transpositions%2) ? -1 : 1; 00532 00533 m_p = m_rowsTranspositions; 00534 00535 m_isInitialized = true; 00536 } 00537 00538 template<typename MatrixType> 00539 typename PartialPivLU<MatrixType>::Scalar PartialPivLU<MatrixType>::determinant() const 00540 { 00541 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 00542 return Scalar(m_det_p) * m_lu.diagonal().prod(); 00543 } 00544 00548 template<typename MatrixType> 00549 MatrixType PartialPivLU<MatrixType>::reconstructedMatrix() const 00550 { 00551 eigen_assert(m_isInitialized && "LU is not initialized."); 00552 // LU 00553 MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix() 00554 * m_lu.template triangularView<Upper>(); 00555 00556 // P^{-1}(LU) 00557 res = m_p.inverse() * res; 00558 00559 return res; 00560 } 00561 00562 /***** Implementation details *****************************************************/ 00563 00564 namespace internal { 00565 00566 /***** Implementation of inverse() *****************************************************/ 00567 template<typename DstXprType, typename MatrixType> 00568 struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename PartialPivLU<MatrixType>::Scalar>, Dense2Dense> 00569 { 00570 typedef PartialPivLU<MatrixType> LuType; 00571 typedef Inverse<LuType> SrcXprType; 00572 static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename LuType::Scalar> &) 00573 { 00574 dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols())); 00575 } 00576 }; 00577 } // end namespace internal 00578 00579 /******** MatrixBase methods *******/ 00580 00587 template<typename Derived> 00588 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject> 00589 MatrixBase<Derived>::partialPivLu() const 00590 { 00591 return PartialPivLU<PlainObject>(eval()); 00592 } 00593 00602 template<typename Derived> 00603 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject> 00604 MatrixBase<Derived>::lu() const 00605 { 00606 return PartialPivLU<PlainObject>(eval()); 00607 } 00608 00609 } // end namespace Eigen 00610 00611 #endif // EIGEN_PARTIALLU_H