00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025 #ifndef EIGEN_LU_H
00026 #define EIGEN_LU_H
00027
00028 namespace Eigen {
00029
00060 template<typename _MatrixType> class FullPivLU
00061 {
00062 public:
00063 typedef _MatrixType MatrixType;
00064 enum {
00065 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00066 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00067 Options = MatrixType::Options,
00068 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00069 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00070 };
00071 typedef typename MatrixType::Scalar Scalar;
00072 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00073 typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
00074 typedef typename MatrixType::Index Index;
00075 typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
00076 typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
00077 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
00078 typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType;
00079
00086 FullPivLU();
00087
00094 FullPivLU(Index rows, Index cols);
00095
00101 FullPivLU(const MatrixType& matrix);
00102
00110 FullPivLU& compute(const MatrixType& matrix);
00111
00118 inline const MatrixType& matrixLU() const
00119 {
00120 eigen_assert(m_isInitialized && "LU is not initialized.");
00121 return m_lu;
00122 }
00123
00131 inline Index nonzeroPivots() const
00132 {
00133 eigen_assert(m_isInitialized && "LU is not initialized.");
00134 return m_nonzero_pivots;
00135 }
00136
00140 RealScalar maxPivot() const { return m_maxpivot; }
00141
00146 inline const PermutationPType& permutationP() const
00147 {
00148 eigen_assert(m_isInitialized && "LU is not initialized.");
00149 return m_p;
00150 }
00151
00156 inline const PermutationQType& permutationQ() const
00157 {
00158 eigen_assert(m_isInitialized && "LU is not initialized.");
00159 return m_q;
00160 }
00161
00176 inline const internal::kernel_retval<FullPivLU> kernel() const
00177 {
00178 eigen_assert(m_isInitialized && "LU is not initialized.");
00179 return internal::kernel_retval<FullPivLU>(*this);
00180 }
00181
00201 inline const internal::image_retval<FullPivLU>
00202 image(const MatrixType& originalMatrix) const
00203 {
00204 eigen_assert(m_isInitialized && "LU is not initialized.");
00205 return internal::image_retval<FullPivLU>(*this, originalMatrix);
00206 }
00207
00227 template<typename Rhs>
00228 inline const internal::solve_retval<FullPivLU, Rhs>
00229 solve(const MatrixBase<Rhs>& b) const
00230 {
00231 eigen_assert(m_isInitialized && "LU is not initialized.");
00232 return internal::solve_retval<FullPivLU, Rhs>(*this, b.derived());
00233 }
00234
00250 typename internal::traits<MatrixType>::Scalar determinant() const;
00251
00269 FullPivLU& setThreshold(const RealScalar& threshold)
00270 {
00271 m_usePrescribedThreshold = true;
00272 m_prescribedThreshold = threshold;
00273 return *this;
00274 }
00275
00284 FullPivLU& setThreshold(Default_t)
00285 {
00286 m_usePrescribedThreshold = false;
00287 return *this;
00288 }
00289
00294 RealScalar threshold() const
00295 {
00296 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00297 return m_usePrescribedThreshold ? m_prescribedThreshold
00298
00299
00300 : NumTraits<Scalar>::epsilon() * m_lu.diagonalSize();
00301 }
00302
00309 inline Index rank() const
00310 {
00311 eigen_assert(m_isInitialized && "LU is not initialized.");
00312 RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
00313 Index result = 0;
00314 for(Index i = 0; i < m_nonzero_pivots; ++i)
00315 result += (internal::abs(m_lu.coeff(i,i)) > premultiplied_threshold);
00316 return result;
00317 }
00318
00325 inline Index dimensionOfKernel() const
00326 {
00327 eigen_assert(m_isInitialized && "LU is not initialized.");
00328 return cols() - rank();
00329 }
00330
00338 inline bool isInjective() const
00339 {
00340 eigen_assert(m_isInitialized && "LU is not initialized.");
00341 return rank() == cols();
00342 }
00343
00351 inline bool isSurjective() const
00352 {
00353 eigen_assert(m_isInitialized && "LU is not initialized.");
00354 return rank() == rows();
00355 }
00356
00363 inline bool isInvertible() const
00364 {
00365 eigen_assert(m_isInitialized && "LU is not initialized.");
00366 return isInjective() && (m_lu.rows() == m_lu.cols());
00367 }
00368
00376 inline const internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType> inverse() const
00377 {
00378 eigen_assert(m_isInitialized && "LU is not initialized.");
00379 eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
00380 return internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType>
00381 (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
00382 }
00383
00384 MatrixType reconstructedMatrix() const;
00385
00386 inline Index rows() const { return m_lu.rows(); }
00387 inline Index cols() const { return m_lu.cols(); }
00388
00389 protected:
00390 MatrixType m_lu;
00391 PermutationPType m_p;
00392 PermutationQType m_q;
00393 IntColVectorType m_rowsTranspositions;
00394 IntRowVectorType m_colsTranspositions;
00395 Index m_det_pq, m_nonzero_pivots;
00396 RealScalar m_maxpivot, m_prescribedThreshold;
00397 bool m_isInitialized, m_usePrescribedThreshold;
00398 };
00399
00400 template<typename MatrixType>
00401 FullPivLU<MatrixType>::FullPivLU()
00402 : m_isInitialized(false), m_usePrescribedThreshold(false)
00403 {
00404 }
00405
00406 template<typename MatrixType>
00407 FullPivLU<MatrixType>::FullPivLU(Index rows, Index cols)
00408 : m_lu(rows, cols),
00409 m_p(rows),
00410 m_q(cols),
00411 m_rowsTranspositions(rows),
00412 m_colsTranspositions(cols),
00413 m_isInitialized(false),
00414 m_usePrescribedThreshold(false)
00415 {
00416 }
00417
00418 template<typename MatrixType>
00419 FullPivLU<MatrixType>::FullPivLU(const MatrixType& matrix)
00420 : m_lu(matrix.rows(), matrix.cols()),
00421 m_p(matrix.rows()),
00422 m_q(matrix.cols()),
00423 m_rowsTranspositions(matrix.rows()),
00424 m_colsTranspositions(matrix.cols()),
00425 m_isInitialized(false),
00426 m_usePrescribedThreshold(false)
00427 {
00428 compute(matrix);
00429 }
00430
00431 template<typename MatrixType>
00432 FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
00433 {
00434 m_isInitialized = true;
00435 m_lu = matrix;
00436
00437 const Index size = matrix.diagonalSize();
00438 const Index rows = matrix.rows();
00439 const Index cols = matrix.cols();
00440
00441
00442
00443 m_rowsTranspositions.resize(matrix.rows());
00444 m_colsTranspositions.resize(matrix.cols());
00445 Index number_of_transpositions = 0;
00446
00447 m_nonzero_pivots = size;
00448 m_maxpivot = RealScalar(0);
00449
00450 for(Index k = 0; k < size; ++k)
00451 {
00452
00453
00454
00455 Index row_of_biggest_in_corner, col_of_biggest_in_corner;
00456 RealScalar biggest_in_corner;
00457 biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
00458 .cwiseAbs()
00459 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
00460 row_of_biggest_in_corner += k;
00461 col_of_biggest_in_corner += k;
00462
00463 if(biggest_in_corner==RealScalar(0))
00464 {
00465
00466
00467 m_nonzero_pivots = k;
00468 for(Index i = k; i < size; ++i)
00469 {
00470 m_rowsTranspositions.coeffRef(i) = i;
00471 m_colsTranspositions.coeffRef(i) = i;
00472 }
00473 break;
00474 }
00475
00476 if(biggest_in_corner > m_maxpivot) m_maxpivot = biggest_in_corner;
00477
00478
00479
00480
00481 m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
00482 m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
00483 if(k != row_of_biggest_in_corner) {
00484 m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
00485 ++number_of_transpositions;
00486 }
00487 if(k != col_of_biggest_in_corner) {
00488 m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
00489 ++number_of_transpositions;
00490 }
00491
00492
00493
00494
00495 if(k<rows-1)
00496 m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
00497 if(k<size-1)
00498 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);
00499 }
00500
00501
00502
00503
00504 m_p.setIdentity(rows);
00505 for(Index k = size-1; k >= 0; --k)
00506 m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
00507
00508 m_q.setIdentity(cols);
00509 for(Index k = 0; k < size; ++k)
00510 m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
00511
00512 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00513 return *this;
00514 }
00515
00516 template<typename MatrixType>
00517 typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
00518 {
00519 eigen_assert(m_isInitialized && "LU is not initialized.");
00520 eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
00521 return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
00522 }
00523
00527 template<typename MatrixType>
00528 MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
00529 {
00530 eigen_assert(m_isInitialized && "LU is not initialized.");
00531 const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols());
00532
00533 MatrixType res(m_lu.rows(),m_lu.cols());
00534
00535 res = m_lu.leftCols(smalldim)
00536 .template triangularView<UnitLower>().toDenseMatrix()
00537 * m_lu.topRows(smalldim)
00538 .template triangularView<Upper>().toDenseMatrix();
00539
00540
00541 res = m_p.inverse() * res;
00542
00543
00544 res = res * m_q.inverse();
00545
00546 return res;
00547 }
00548
00549
00550
00551 namespace internal {
00552 template<typename _MatrixType>
00553 struct kernel_retval<FullPivLU<_MatrixType> >
00554 : kernel_retval_base<FullPivLU<_MatrixType> >
00555 {
00556 EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<_MatrixType>)
00557
00558 enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
00559 MatrixType::MaxColsAtCompileTime,
00560 MatrixType::MaxRowsAtCompileTime)
00561 };
00562
00563 template<typename Dest> void evalTo(Dest& dst) const
00564 {
00565 const Index cols = dec().matrixLU().cols(), dimker = cols - rank();
00566 if(dimker == 0)
00567 {
00568
00569
00570
00571 dst.setZero();
00572 return;
00573 }
00574
00575
00576
00577
00578
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591 Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
00592 RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
00593 Index p = 0;
00594 for(Index i = 0; i < dec().nonzeroPivots(); ++i)
00595 if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
00596 pivots.coeffRef(p++) = i;
00597 eigen_internal_assert(p == rank());
00598
00599
00600
00601
00602
00603 Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
00604 MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
00605 m(dec().matrixLU().block(0, 0, rank(), cols));
00606 for(Index i = 0; i < rank(); ++i)
00607 {
00608 if(i) m.row(i).head(i).setZero();
00609 m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i);
00610 }
00611 m.block(0, 0, rank(), rank());
00612 m.block(0, 0, rank(), rank()).template triangularView<StrictlyLower>().setZero();
00613 for(Index i = 0; i < rank(); ++i)
00614 m.col(i).swap(m.col(pivots.coeff(i)));
00615
00616
00617
00618
00619 m.topLeftCorner(rank(), rank())
00620 .template triangularView<Upper>().solveInPlace(
00621 m.topRightCorner(rank(), dimker)
00622 );
00623
00624
00625 for(Index i = rank()-1; i >= 0; --i)
00626 m.col(i).swap(m.col(pivots.coeff(i)));
00627
00628
00629 for(Index i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker);
00630 for(Index i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero();
00631 for(Index k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1);
00632 }
00633 };
00634
00635
00636
00637 template<typename _MatrixType>
00638 struct image_retval<FullPivLU<_MatrixType> >
00639 : image_retval_base<FullPivLU<_MatrixType> >
00640 {
00641 EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<_MatrixType>)
00642
00643 enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
00644 MatrixType::MaxColsAtCompileTime,
00645 MatrixType::MaxRowsAtCompileTime)
00646 };
00647
00648 template<typename Dest> void evalTo(Dest& dst) const
00649 {
00650 if(rank() == 0)
00651 {
00652
00653
00654
00655 dst.setZero();
00656 return;
00657 }
00658
00659 Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
00660 RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
00661 Index p = 0;
00662 for(Index i = 0; i < dec().nonzeroPivots(); ++i)
00663 if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
00664 pivots.coeffRef(p++) = i;
00665 eigen_internal_assert(p == rank());
00666
00667 for(Index i = 0; i < rank(); ++i)
00668 dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i)));
00669 }
00670 };
00671
00672
00673
00674 template<typename _MatrixType, typename Rhs>
00675 struct solve_retval<FullPivLU<_MatrixType>, Rhs>
00676 : solve_retval_base<FullPivLU<_MatrixType>, Rhs>
00677 {
00678 EIGEN_MAKE_SOLVE_HELPERS(FullPivLU<_MatrixType>,Rhs)
00679
00680 template<typename Dest> void evalTo(Dest& dst) const
00681 {
00682
00683
00684
00685
00686
00687
00688
00689
00690 const Index rows = dec().rows(), cols = dec().cols(),
00691 nonzero_pivots = dec().nonzeroPivots();
00692 eigen_assert(rhs().rows() == rows);
00693 const Index smalldim = (std::min)(rows, cols);
00694
00695 if(nonzero_pivots == 0)
00696 {
00697 dst.setZero();
00698 return;
00699 }
00700
00701 typename Rhs::PlainObject c(rhs().rows(), rhs().cols());
00702
00703
00704 c = dec().permutationP() * rhs();
00705
00706
00707 dec().matrixLU()
00708 .topLeftCorner(smalldim,smalldim)
00709 .template triangularView<UnitLower>()
00710 .solveInPlace(c.topRows(smalldim));
00711 if(rows>cols)
00712 {
00713 c.bottomRows(rows-cols)
00714 -= dec().matrixLU().bottomRows(rows-cols)
00715 * c.topRows(cols);
00716 }
00717
00718
00719 dec().matrixLU()
00720 .topLeftCorner(nonzero_pivots, nonzero_pivots)
00721 .template triangularView<Upper>()
00722 .solveInPlace(c.topRows(nonzero_pivots));
00723
00724
00725 for(Index i = 0; i < nonzero_pivots; ++i)
00726 dst.row(dec().permutationQ().indices().coeff(i)) = c.row(i);
00727 for(Index i = nonzero_pivots; i < dec().matrixLU().cols(); ++i)
00728 dst.row(dec().permutationQ().indices().coeff(i)).setZero();
00729 }
00730 };
00731
00732 }
00733
00734
00735
00742 template<typename Derived>
00743 inline const FullPivLU<typename MatrixBase<Derived>::PlainObject>
00744 MatrixBase<Derived>::fullPivLu() const
00745 {
00746 return FullPivLU<PlainObject>(eval());
00747 }
00748
00749 }
00750
00751 #endif // EIGEN_LU_H