![]() |
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 // 00006 // This Source Code Form is subject to the terms of the Mozilla 00007 // Public License v. 2.0. If a copy of the MPL was not distributed 00008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00009 00010 #ifndef EIGEN_LU_H 00011 #define EIGEN_LU_H 00012 00013 namespace Eigen { 00014 00015 namespace internal { 00016 template<typename _MatrixType> struct traits<FullPivLU<_MatrixType> > 00017 : traits<_MatrixType> 00018 { 00019 typedef MatrixXpr XprKind; 00020 typedef SolverStorage StorageKind; 00021 enum { Flags = 0 }; 00022 }; 00023 00024 } // end namespace internal 00025 00059 template<typename _MatrixType> class FullPivLU 00060 : public SolverBase<FullPivLU<_MatrixType> > 00061 { 00062 public: 00063 typedef _MatrixType MatrixType; 00064 typedef SolverBase<FullPivLU> Base; 00065 00066 EIGEN_GENERIC_PUBLIC_INTERFACE(FullPivLU) 00067 // FIXME StorageIndex defined in EIGEN_GENERIC_PUBLIC_INTERFACE should be int 00068 enum { 00069 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 00070 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00071 }; 00072 typedef typename internal::plain_row_type<MatrixType, StorageIndex>::type IntRowVectorType; 00073 typedef typename internal::plain_col_type<MatrixType, StorageIndex>::type IntColVectorType; 00074 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType; 00075 typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType; 00076 typedef typename MatrixType::PlainObject PlainObject; 00077 00084 FullPivLU(); 00085 00092 FullPivLU(Index rows, Index cols); 00093 00099 template<typename InputType> 00100 explicit FullPivLU(const EigenBase<InputType>& matrix); 00101 00108 template<typename InputType> 00109 explicit FullPivLU(EigenBase<InputType>& matrix); 00110 00118 template<typename InputType> 00119 FullPivLU& compute(const EigenBase<InputType>& matrix) { 00120 m_lu = matrix.derived(); 00121 computeInPlace(); 00122 return *this; 00123 } 00124 00131 inline const MatrixType& matrixLU() const 00132 { 00133 eigen_assert(m_isInitialized && "LU is not initialized."); 00134 return m_lu; 00135 } 00136 00144 inline Index nonzeroPivots() const 00145 { 00146 eigen_assert(m_isInitialized && "LU is not initialized."); 00147 return m_nonzero_pivots; 00148 } 00149 00153 RealScalar maxPivot() const { return m_maxpivot; } 00154 00159 EIGEN_DEVICE_FUNC inline const PermutationPType& permutationP() const 00160 { 00161 eigen_assert(m_isInitialized && "LU is not initialized."); 00162 return m_p; 00163 } 00164 00169 inline const PermutationQType& permutationQ() const 00170 { 00171 eigen_assert(m_isInitialized && "LU is not initialized."); 00172 return m_q; 00173 } 00174 00189 inline const internal::kernel_retval<FullPivLU> kernel() const 00190 { 00191 eigen_assert(m_isInitialized && "LU is not initialized."); 00192 return internal::kernel_retval<FullPivLU>(*this); 00193 } 00194 00214 inline const internal::image_retval<FullPivLU> 00215 image(const MatrixType& originalMatrix) const 00216 { 00217 eigen_assert(m_isInitialized && "LU is not initialized."); 00218 return internal::image_retval<FullPivLU>(*this, originalMatrix); 00219 } 00220 00240 // FIXME this is a copy-paste of the base-class member to add the isInitialized assertion. 00241 template<typename Rhs> 00242 inline const Solve<FullPivLU, Rhs> 00243 solve(const MatrixBase<Rhs>& b) const 00244 { 00245 eigen_assert(m_isInitialized && "LU is not initialized."); 00246 return Solve<FullPivLU, Rhs>(*this, b.derived()); 00247 } 00248 00252 inline RealScalar rcond() const 00253 { 00254 eigen_assert(m_isInitialized && "PartialPivLU is not initialized."); 00255 return internal::rcond_estimate_helper(m_l1_norm, *this); 00256 } 00257 00273 typename internal::traits<MatrixType>::Scalar determinant() const; 00274 00292 FullPivLU& setThreshold(const RealScalar& threshold) 00293 { 00294 m_usePrescribedThreshold = true; 00295 m_prescribedThreshold = threshold; 00296 return *this; 00297 } 00298 00307 FullPivLU& setThreshold(Default_t) 00308 { 00309 m_usePrescribedThreshold = false; 00310 return *this; 00311 } 00312 00317 RealScalar threshold() const 00318 { 00319 eigen_assert(m_isInitialized || m_usePrescribedThreshold); 00320 return m_usePrescribedThreshold ? m_prescribedThreshold 00321 // this formula comes from experimenting (see "LU precision tuning" thread on the list) 00322 // and turns out to be identical to Higham's formula used already in LDLt. 00323 : NumTraits<Scalar>::epsilon() * m_lu.diagonalSize(); 00324 } 00325 00332 inline Index rank() const 00333 { 00334 using std::abs; 00335 eigen_assert(m_isInitialized && "LU is not initialized."); 00336 RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); 00337 Index result = 0; 00338 for(Index i = 0; i < m_nonzero_pivots; ++i) 00339 result += (abs(m_lu.coeff(i,i)) > premultiplied_threshold); 00340 return result; 00341 } 00342 00349 inline Index dimensionOfKernel() const 00350 { 00351 eigen_assert(m_isInitialized && "LU is not initialized."); 00352 return cols() - rank(); 00353 } 00354 00362 inline bool isInjective() const 00363 { 00364 eigen_assert(m_isInitialized && "LU is not initialized."); 00365 return rank() == cols(); 00366 } 00367 00375 inline bool isSurjective() const 00376 { 00377 eigen_assert(m_isInitialized && "LU is not initialized."); 00378 return rank() == rows(); 00379 } 00380 00387 inline bool isInvertible() const 00388 { 00389 eigen_assert(m_isInitialized && "LU is not initialized."); 00390 return isInjective() && (m_lu.rows() == m_lu.cols()); 00391 } 00392 00400 inline const Inverse<FullPivLU> inverse() const 00401 { 00402 eigen_assert(m_isInitialized && "LU is not initialized."); 00403 eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!"); 00404 return Inverse<FullPivLU>(*this); 00405 } 00406 00407 MatrixType reconstructedMatrix() const; 00408 00409 EIGEN_DEVICE_FUNC inline Index rows() const { return m_lu.rows(); } 00410 EIGEN_DEVICE_FUNC inline Index cols() const { return m_lu.cols(); } 00411 00412 #ifndef EIGEN_PARSED_BY_DOXYGEN 00413 template<typename RhsType, typename DstType> 00414 EIGEN_DEVICE_FUNC 00415 void _solve_impl(const RhsType &rhs, DstType &dst) const; 00416 00417 template<bool Conjugate, typename RhsType, typename DstType> 00418 EIGEN_DEVICE_FUNC 00419 void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const; 00420 #endif 00421 00422 protected: 00423 00424 static void check_template_parameters() 00425 { 00426 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00427 } 00428 00429 void computeInPlace(); 00430 00431 MatrixType m_lu; 00432 PermutationPType m_p; 00433 PermutationQType m_q; 00434 IntColVectorType m_rowsTranspositions; 00435 IntRowVectorType m_colsTranspositions; 00436 Index m_nonzero_pivots; 00437 RealScalar m_l1_norm; 00438 RealScalar m_maxpivot, m_prescribedThreshold; 00439 signed char m_det_pq; 00440 bool m_isInitialized, m_usePrescribedThreshold; 00441 }; 00442 00443 template<typename MatrixType> 00444 FullPivLU<MatrixType>::FullPivLU() 00445 : m_isInitialized(false), m_usePrescribedThreshold(false) 00446 { 00447 } 00448 00449 template<typename MatrixType> 00450 FullPivLU<MatrixType>::FullPivLU(Index rows, Index cols) 00451 : m_lu(rows, cols), 00452 m_p(rows), 00453 m_q(cols), 00454 m_rowsTranspositions(rows), 00455 m_colsTranspositions(cols), 00456 m_isInitialized(false), 00457 m_usePrescribedThreshold(false) 00458 { 00459 } 00460 00461 template<typename MatrixType> 00462 template<typename InputType> 00463 FullPivLU<MatrixType>::FullPivLU(const EigenBase<InputType>& matrix) 00464 : m_lu(matrix.rows(), matrix.cols()), 00465 m_p(matrix.rows()), 00466 m_q(matrix.cols()), 00467 m_rowsTranspositions(matrix.rows()), 00468 m_colsTranspositions(matrix.cols()), 00469 m_isInitialized(false), 00470 m_usePrescribedThreshold(false) 00471 { 00472 compute(matrix.derived()); 00473 } 00474 00475 template<typename MatrixType> 00476 template<typename InputType> 00477 FullPivLU<MatrixType>::FullPivLU(EigenBase<InputType>& matrix) 00478 : m_lu(matrix.derived()), 00479 m_p(matrix.rows()), 00480 m_q(matrix.cols()), 00481 m_rowsTranspositions(matrix.rows()), 00482 m_colsTranspositions(matrix.cols()), 00483 m_isInitialized(false), 00484 m_usePrescribedThreshold(false) 00485 { 00486 computeInPlace(); 00487 } 00488 00489 template<typename MatrixType> 00490 void FullPivLU<MatrixType>::computeInPlace() 00491 { 00492 check_template_parameters(); 00493 00494 // the permutations are stored as int indices, so just to be sure: 00495 eigen_assert(m_lu.rows()<=NumTraits<int>::highest() && m_lu.cols()<=NumTraits<int>::highest()); 00496 00497 m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff(); 00498 00499 const Index size = m_lu.diagonalSize(); 00500 const Index rows = m_lu.rows(); 00501 const Index cols = m_lu.cols(); 00502 00503 // will store the transpositions, before we accumulate them at the end. 00504 // can't accumulate on-the-fly because that will be done in reverse order for the rows. 00505 m_rowsTranspositions.resize(m_lu.rows()); 00506 m_colsTranspositions.resize(m_lu.cols()); 00507 Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i 00508 00509 m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) 00510 m_maxpivot = RealScalar(0); 00511 00512 for(Index k = 0; k < size; ++k) 00513 { 00514 // First, we need to find the pivot. 00515 00516 // biggest coefficient in the remaining bottom-right corner (starting at row k, col k) 00517 Index row_of_biggest_in_corner, col_of_biggest_in_corner; 00518 typedef internal::scalar_score_coeff_op<Scalar> Scoring; 00519 typedef typename Scoring::result_type Score; 00520 Score biggest_in_corner; 00521 biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k) 00522 .unaryExpr(Scoring()) 00523 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner); 00524 row_of_biggest_in_corner += k; // correct the values! since they were computed in the corner, 00525 col_of_biggest_in_corner += k; // need to add k to them. 00526 00527 if(biggest_in_corner==Score(0)) 00528 { 00529 // before exiting, make sure to initialize the still uninitialized transpositions 00530 // in a sane state without destroying what we already have. 00531 m_nonzero_pivots = k; 00532 for(Index i = k; i < size; ++i) 00533 { 00534 m_rowsTranspositions.coeffRef(i) = i; 00535 m_colsTranspositions.coeffRef(i) = i; 00536 } 00537 break; 00538 } 00539 00540 RealScalar abs_pivot = internal::abs_knowing_score<Scalar>()(m_lu(row_of_biggest_in_corner, col_of_biggest_in_corner), biggest_in_corner); 00541 if(abs_pivot > m_maxpivot) m_maxpivot = abs_pivot; 00542 00543 // Now that we've found the pivot, we need to apply the row/col swaps to 00544 // bring it to the location (k,k). 00545 00546 m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner; 00547 m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner; 00548 if(k != row_of_biggest_in_corner) { 00549 m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner)); 00550 ++number_of_transpositions; 00551 } 00552 if(k != col_of_biggest_in_corner) { 00553 m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner)); 00554 ++number_of_transpositions; 00555 } 00556 00557 // Now that the pivot is at the right location, we update the remaining 00558 // bottom-right corner by Gaussian elimination. 00559 00560 if(k<rows-1) 00561 m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k); 00562 if(k<size-1) 00563 m_lu.block(k+1,k+1,rows-k-1,cols-k-1).noalias() -= m_lu.col(k).tail(rows-k-1) * m_lu.row(k).tail(cols-k-1); 00564 } 00565 00566 // the main loop is over, we still have to accumulate the transpositions to find the 00567 // permutations P and Q 00568 00569 m_p.setIdentity(rows); 00570 for(Index k = size-1; k >= 0; --k) 00571 m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k)); 00572 00573 m_q.setIdentity(cols); 00574 for(Index k = 0; k < size; ++k) 00575 m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k)); 00576 00577 m_det_pq = (number_of_transpositions%2) ? -1 : 1; 00578 00579 m_isInitialized = true; 00580 } 00581 00582 template<typename MatrixType> 00583 typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const 00584 { 00585 eigen_assert(m_isInitialized && "LU is not initialized."); 00586 eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!"); 00587 return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod()); 00588 } 00589 00593 template<typename MatrixType> 00594 MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const 00595 { 00596 eigen_assert(m_isInitialized && "LU is not initialized."); 00597 const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols()); 00598 // LU 00599 MatrixType res(m_lu.rows(),m_lu.cols()); 00600 // FIXME the .toDenseMatrix() should not be needed... 00601 res = m_lu.leftCols(smalldim) 00602 .template triangularView<UnitLower>().toDenseMatrix() 00603 * m_lu.topRows(smalldim) 00604 .template triangularView<Upper>().toDenseMatrix(); 00605 00606 // P^{-1}(LU) 00607 res = m_p.inverse() * res; 00608 00609 // (P^{-1}LU)Q^{-1} 00610 res = res * m_q.inverse(); 00611 00612 return res; 00613 } 00614 00615 /********* Implementation of kernel() **************************************************/ 00616 00617 namespace internal { 00618 template<typename _MatrixType> 00619 struct kernel_retval<FullPivLU<_MatrixType> > 00620 : kernel_retval_base<FullPivLU<_MatrixType> > 00621 { 00622 EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<_MatrixType>) 00623 00624 enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED( 00625 MatrixType::MaxColsAtCompileTime, 00626 MatrixType::MaxRowsAtCompileTime) 00627 }; 00628 00629 template<typename Dest> void evalTo(Dest& dst) const 00630 { 00631 using std::abs; 00632 const Index cols = dec().matrixLU().cols(), dimker = cols - rank(); 00633 if(dimker == 0) 00634 { 00635 // The Kernel is just {0}, so it doesn't have a basis properly speaking, but let's 00636 // avoid crashing/asserting as that depends on floating point calculations. Let's 00637 // just return a single column vector filled with zeros. 00638 dst.setZero(); 00639 return; 00640 } 00641 00642 /* Let us use the following lemma: 00643 * 00644 * Lemma: If the matrix A has the LU decomposition PAQ = LU, 00645 * then Ker A = Q(Ker U). 00646 * 00647 * Proof: trivial: just keep in mind that P, Q, L are invertible. 00648 */ 00649 00650 /* Thus, all we need to do is to compute Ker U, and then apply Q. 00651 * 00652 * U is upper triangular, with eigenvalues sorted so that any zeros appear at the end. 00653 * Thus, the diagonal of U ends with exactly 00654 * dimKer zero's. Let us use that to construct dimKer linearly 00655 * independent vectors in Ker U. 00656 */ 00657 00658 Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank()); 00659 RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold(); 00660 Index p = 0; 00661 for(Index i = 0; i < dec().nonzeroPivots(); ++i) 00662 if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold) 00663 pivots.coeffRef(p++) = i; 00664 eigen_internal_assert(p == rank()); 00665 00666 // we construct a temporaty trapezoid matrix m, by taking the U matrix and 00667 // permuting the rows and cols to bring the nonnegligible pivots to the top of 00668 // the main diagonal. We need that to be able to apply our triangular solvers. 00669 // FIXME when we get triangularView-for-rectangular-matrices, this can be simplified 00670 Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options, 00671 MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime> 00672 m(dec().matrixLU().block(0, 0, rank(), cols)); 00673 for(Index i = 0; i < rank(); ++i) 00674 { 00675 if(i) m.row(i).head(i).setZero(); 00676 m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i); 00677 } 00678 m.block(0, 0, rank(), rank()); 00679 m.block(0, 0, rank(), rank()).template triangularView<StrictlyLower>().setZero(); 00680 for(Index i = 0; i < rank(); ++i) 00681 m.col(i).swap(m.col(pivots.coeff(i))); 00682 00683 // ok, we have our trapezoid matrix, we can apply the triangular solver. 00684 // notice that the math behind this suggests that we should apply this to the 00685 // negative of the RHS, but for performance we just put the negative sign elsewhere, see below. 00686 m.topLeftCorner(rank(), rank()) 00687 .template triangularView<Upper>().solveInPlace( 00688 m.topRightCorner(rank(), dimker) 00689 ); 00690 00691 // now we must undo the column permutation that we had applied! 00692 for(Index i = rank()-1; i >= 0; --i) 00693 m.col(i).swap(m.col(pivots.coeff(i))); 00694 00695 // see the negative sign in the next line, that's what we were talking about above. 00696 for(Index i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker); 00697 for(Index i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero(); 00698 for(Index k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1); 00699 } 00700 }; 00701 00702 /***** Implementation of image() *****************************************************/ 00703 00704 template<typename _MatrixType> 00705 struct image_retval<FullPivLU<_MatrixType> > 00706 : image_retval_base<FullPivLU<_MatrixType> > 00707 { 00708 EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<_MatrixType>) 00709 00710 enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED( 00711 MatrixType::MaxColsAtCompileTime, 00712 MatrixType::MaxRowsAtCompileTime) 00713 }; 00714 00715 template<typename Dest> void evalTo(Dest& dst) const 00716 { 00717 using std::abs; 00718 if(rank() == 0) 00719 { 00720 // The Image is just {0}, so it doesn't have a basis properly speaking, but let's 00721 // avoid crashing/asserting as that depends on floating point calculations. Let's 00722 // just return a single column vector filled with zeros. 00723 dst.setZero(); 00724 return; 00725 } 00726 00727 Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank()); 00728 RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold(); 00729 Index p = 0; 00730 for(Index i = 0; i < dec().nonzeroPivots(); ++i) 00731 if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold) 00732 pivots.coeffRef(p++) = i; 00733 eigen_internal_assert(p == rank()); 00734 00735 for(Index i = 0; i < rank(); ++i) 00736 dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i))); 00737 } 00738 }; 00739 00740 /***** Implementation of solve() *****************************************************/ 00741 00742 } // end namespace internal 00743 00744 #ifndef EIGEN_PARSED_BY_DOXYGEN 00745 template<typename _MatrixType> 00746 template<typename RhsType, typename DstType> 00747 void FullPivLU<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const 00748 { 00749 /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}. 00750 * So we proceed as follows: 00751 * Step 1: compute c = P * rhs. 00752 * Step 2: replace c by the solution x to Lx = c. Exists because L is invertible. 00753 * Step 3: replace c by the solution x to Ux = c. May or may not exist. 00754 * Step 4: result = Q * c; 00755 */ 00756 00757 const Index rows = this->rows(), 00758 cols = this->cols(), 00759 nonzero_pivots = this->rank(); 00760 eigen_assert(rhs.rows() == rows); 00761 const Index smalldim = (std::min)(rows, cols); 00762 00763 if(nonzero_pivots == 0) 00764 { 00765 dst.setZero(); 00766 return; 00767 } 00768 00769 typename RhsType::PlainObject c(rhs.rows(), rhs.cols()); 00770 00771 // Step 1 00772 c = permutationP() * rhs; 00773 00774 // Step 2 00775 m_lu.topLeftCorner(smalldim,smalldim) 00776 .template triangularView<UnitLower>() 00777 .solveInPlace(c.topRows(smalldim)); 00778 if(rows>cols) 00779 c.bottomRows(rows-cols) -= m_lu.bottomRows(rows-cols) * c.topRows(cols); 00780 00781 // Step 3 00782 m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots) 00783 .template triangularView<Upper>() 00784 .solveInPlace(c.topRows(nonzero_pivots)); 00785 00786 // Step 4 00787 for(Index i = 0; i < nonzero_pivots; ++i) 00788 dst.row(permutationQ().indices().coeff(i)) = c.row(i); 00789 for(Index i = nonzero_pivots; i < m_lu.cols(); ++i) 00790 dst.row(permutationQ().indices().coeff(i)).setZero(); 00791 } 00792 00793 template<typename _MatrixType> 00794 template<bool Conjugate, typename RhsType, typename DstType> 00795 void FullPivLU<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const 00796 { 00797 /* The decomposition PAQ = LU can be rewritten as A = P^{-1} L U Q^{-1}, 00798 * and since permutations are real and unitary, we can write this 00799 * as A^T = Q U^T L^T P, 00800 * So we proceed as follows: 00801 * Step 1: compute c = Q^T rhs. 00802 * Step 2: replace c by the solution x to U^T x = c. May or may not exist. 00803 * Step 3: replace c by the solution x to L^T x = c. 00804 * Step 4: result = P^T c. 00805 * If Conjugate is true, replace "^T" by "^*" above. 00806 */ 00807 00808 const Index rows = this->rows(), cols = this->cols(), 00809 nonzero_pivots = this->rank(); 00810 eigen_assert(rhs.rows() == cols); 00811 const Index smalldim = (std::min)(rows, cols); 00812 00813 if(nonzero_pivots == 0) 00814 { 00815 dst.setZero(); 00816 return; 00817 } 00818 00819 typename RhsType::PlainObject c(rhs.rows(), rhs.cols()); 00820 00821 // Step 1 00822 c = permutationQ().inverse() * rhs; 00823 00824 if (Conjugate) { 00825 // Step 2 00826 m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots) 00827 .template triangularView<Upper>() 00828 .adjoint() 00829 .solveInPlace(c.topRows(nonzero_pivots)); 00830 // Step 3 00831 m_lu.topLeftCorner(smalldim, smalldim) 00832 .template triangularView<UnitLower>() 00833 .adjoint() 00834 .solveInPlace(c.topRows(smalldim)); 00835 } else { 00836 // Step 2 00837 m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots) 00838 .template triangularView<Upper>() 00839 .transpose() 00840 .solveInPlace(c.topRows(nonzero_pivots)); 00841 // Step 3 00842 m_lu.topLeftCorner(smalldim, smalldim) 00843 .template triangularView<UnitLower>() 00844 .transpose() 00845 .solveInPlace(c.topRows(smalldim)); 00846 } 00847 00848 // Step 4 00849 PermutationPType invp = permutationP().inverse().eval(); 00850 for(Index i = 0; i < smalldim; ++i) 00851 dst.row(invp.indices().coeff(i)) = c.row(i); 00852 for(Index i = smalldim; i < rows; ++i) 00853 dst.row(invp.indices().coeff(i)).setZero(); 00854 } 00855 00856 #endif 00857 00858 namespace internal { 00859 00860 00861 /***** Implementation of inverse() *****************************************************/ 00862 template<typename DstXprType, typename MatrixType> 00863 struct Assignment<DstXprType, Inverse<FullPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivLU<MatrixType>::Scalar>, Dense2Dense> 00864 { 00865 typedef FullPivLU<MatrixType> LuType; 00866 typedef Inverse<LuType> SrcXprType; 00867 static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename MatrixType::Scalar> &) 00868 { 00869 dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols())); 00870 } 00871 }; 00872 } // end namespace internal 00873 00874 /******* MatrixBase methods *****************************************************************/ 00875 00882 template<typename Derived> 00883 inline const FullPivLU<typename MatrixBase<Derived>::PlainObject> 00884 MatrixBase<Derived>::fullPivLu() const 00885 { 00886 return FullPivLU<PlainObject>(eval()); 00887 } 00888 00889 } // end namespace Eigen 00890 00891 #endif // EIGEN_LU_H