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
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058 template<typename _MatrixType> class FullPivLU
00059 {
00060 public:
00061 typedef _MatrixType MatrixType;
00062 enum {
00063 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00064 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00065 Options = MatrixType::Options,
00066 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00067 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00068 };
00069 typedef typename MatrixType::Scalar Scalar;
00070 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
00071 typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
00072 typedef typename MatrixType::Index Index;
00073 typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
00074 typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
00075 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
00076 typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType;
00077
00078
00079
00080
00081
00082
00083
00084 FullPivLU();
00085
00086
00087
00088
00089
00090
00091
00092 FullPivLU(Index rows, Index cols);
00093
00094
00095
00096
00097
00098
00099 FullPivLU(const MatrixType& matrix);
00100
00101
00102
00103
00104
00105
00106
00107
00108 FullPivLU& compute(const MatrixType& matrix);
00109
00110
00111
00112
00113
00114
00115
00116 inline const MatrixType& matrixLU() const
00117 {
00118 eigen_assert(m_isInitialized && "LU is not initialized.");
00119 return m_lu;
00120 }
00121
00122
00123
00124
00125
00126
00127
00128
00129 inline Index nonzeroPivots() const
00130 {
00131 eigen_assert(m_isInitialized && "LU is not initialized.");
00132 return m_nonzero_pivots;
00133 }
00134
00135
00136
00137
00138 RealScalar maxPivot() const { return m_maxpivot; }
00139
00140
00141
00142
00143
00144 inline const PermutationPType& permutationP() const
00145 {
00146 eigen_assert(m_isInitialized && "LU is not initialized.");
00147 return m_p;
00148 }
00149
00150
00151
00152
00153
00154 inline const PermutationQType& permutationQ() const
00155 {
00156 eigen_assert(m_isInitialized && "LU is not initialized.");
00157 return m_q;
00158 }
00159
00160
00161
00162
00163
00164
00165
00166
00167
00168
00169
00170
00171
00172
00173
00174 inline const internal::kernel_retval<FullPivLU> kernel() const
00175 {
00176 eigen_assert(m_isInitialized && "LU is not initialized.");
00177 return internal::kernel_retval<FullPivLU>(*this);
00178 }
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199 inline const internal::image_retval<FullPivLU>
00200 image(const MatrixType& originalMatrix) const
00201 {
00202 eigen_assert(m_isInitialized && "LU is not initialized.");
00203 return internal::image_retval<FullPivLU>(*this, originalMatrix);
00204 }
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225 template<typename Rhs>
00226 inline const internal::solve_retval<FullPivLU, Rhs>
00227 solve(const MatrixBase<Rhs>& b) const
00228 {
00229 eigen_assert(m_isInitialized && "LU is not initialized.");
00230 return internal::solve_retval<FullPivLU, Rhs>(*this, b.derived());
00231 }
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248 typename internal::traits<MatrixType>::Scalar determinant() const;
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267 FullPivLU& setThreshold(const RealScalar& threshold)
00268 {
00269 m_usePrescribedThreshold = true;
00270 m_prescribedThreshold = threshold;
00271 return *this;
00272 }
00273
00274
00275
00276
00277
00278
00279
00280
00281
00282 FullPivLU& setThreshold(Default_t)
00283 {
00284 m_usePrescribedThreshold = false;
00285 }
00286
00287
00288
00289
00290
00291 RealScalar threshold() const
00292 {
00293 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00294 return m_usePrescribedThreshold ? m_prescribedThreshold
00295
00296
00297 : NumTraits<Scalar>::epsilon() * m_lu.diagonalSize();
00298 }
00299
00300
00301
00302
00303
00304
00305
00306 inline Index rank() const
00307 {
00308 eigen_assert(m_isInitialized && "LU is not initialized.");
00309 RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
00310 Index result = 0;
00311 for(Index i = 0; i < m_nonzero_pivots; ++i)
00312 result += (internal::abs(m_lu.coeff(i,i)) > premultiplied_threshold);
00313 return result;
00314 }
00315
00316
00317
00318
00319
00320
00321
00322 inline Index dimensionOfKernel() const
00323 {
00324 eigen_assert(m_isInitialized && "LU is not initialized.");
00325 return cols() - rank();
00326 }
00327
00328
00329
00330
00331
00332
00333
00334
00335 inline bool isInjective() const
00336 {
00337 eigen_assert(m_isInitialized && "LU is not initialized.");
00338 return rank() == cols();
00339 }
00340
00341
00342
00343
00344
00345
00346
00347
00348 inline bool isSurjective() const
00349 {
00350 eigen_assert(m_isInitialized && "LU is not initialized.");
00351 return rank() == rows();
00352 }
00353
00354
00355
00356
00357
00358
00359
00360 inline bool isInvertible() const
00361 {
00362 eigen_assert(m_isInitialized && "LU is not initialized.");
00363 return isInjective() && (m_lu.rows() == m_lu.cols());
00364 }
00365
00366
00367
00368
00369
00370
00371
00372
00373 inline const internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType> inverse() const
00374 {
00375 eigen_assert(m_isInitialized && "LU is not initialized.");
00376 eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the inverse of a non-square matrix!");
00377 return internal::solve_retval<FullPivLU,typename MatrixType::IdentityReturnType>
00378 (*this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
00379 }
00380
00381 MatrixType reconstructedMatrix() const;
00382
00383 inline Index rows() const { return m_lu.rows(); }
00384 inline Index cols() const { return m_lu.cols(); }
00385
00386 protected:
00387 MatrixType m_lu;
00388 PermutationPType m_p;
00389 PermutationQType m_q;
00390 IntColVectorType m_rowsTranspositions;
00391 IntRowVectorType m_colsTranspositions;
00392 Index m_det_pq, m_nonzero_pivots;
00393 RealScalar m_maxpivot, m_prescribedThreshold;
00394 bool m_isInitialized, m_usePrescribedThreshold;
00395 };
00396
00397 template<typename MatrixType>
00398 FullPivLU<MatrixType>::FullPivLU()
00399 : m_isInitialized(false), m_usePrescribedThreshold(false)
00400 {
00401 }
00402
00403 template<typename MatrixType>
00404 FullPivLU<MatrixType>::FullPivLU(Index rows, Index cols)
00405 : m_lu(rows, cols),
00406 m_p(rows),
00407 m_q(cols),
00408 m_rowsTranspositions(rows),
00409 m_colsTranspositions(cols),
00410 m_isInitialized(false),
00411 m_usePrescribedThreshold(false)
00412 {
00413 }
00414
00415 template<typename MatrixType>
00416 FullPivLU<MatrixType>::FullPivLU(const MatrixType& matrix)
00417 : m_lu(matrix.rows(), matrix.cols()),
00418 m_p(matrix.rows()),
00419 m_q(matrix.cols()),
00420 m_rowsTranspositions(matrix.rows()),
00421 m_colsTranspositions(matrix.cols()),
00422 m_isInitialized(false),
00423 m_usePrescribedThreshold(false)
00424 {
00425 compute(matrix);
00426 }
00427
00428 template<typename MatrixType>
00429 FullPivLU<MatrixType>& FullPivLU<MatrixType>::compute(const MatrixType& matrix)
00430 {
00431 m_isInitialized = true;
00432 m_lu = matrix;
00433
00434 const Index size = matrix.diagonalSize();
00435 const Index rows = matrix.rows();
00436 const Index cols = matrix.cols();
00437
00438
00439
00440 m_rowsTranspositions.resize(matrix.rows());
00441 m_colsTranspositions.resize(matrix.cols());
00442 Index number_of_transpositions = 0;
00443
00444 m_nonzero_pivots = size;
00445 m_maxpivot = RealScalar(0);
00446 RealScalar cutoff(0);
00447
00448 for(Index k = 0; k < size; ++k)
00449 {
00450
00451
00452
00453 Index row_of_biggest_in_corner, col_of_biggest_in_corner;
00454 RealScalar biggest_in_corner;
00455 biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
00456 .cwiseAbs()
00457 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
00458 row_of_biggest_in_corner += k;
00459 col_of_biggest_in_corner += k;
00460
00461
00462 if(k == 0) cutoff = biggest_in_corner * NumTraits<Scalar>::epsilon();
00463
00464
00465
00466
00467
00468 if(biggest_in_corner < cutoff)
00469 {
00470
00471
00472 m_nonzero_pivots = k;
00473 for(Index i = k; i < size; ++i)
00474 {
00475 m_rowsTranspositions.coeffRef(i) = i;
00476 m_colsTranspositions.coeffRef(i) = i;
00477 }
00478 break;
00479 }
00480
00481 if(biggest_in_corner > m_maxpivot) m_maxpivot = biggest_in_corner;
00482
00483
00484
00485
00486 m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
00487 m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
00488 if(k != row_of_biggest_in_corner) {
00489 m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
00490 ++number_of_transpositions;
00491 }
00492 if(k != col_of_biggest_in_corner) {
00493 m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
00494 ++number_of_transpositions;
00495 }
00496
00497
00498
00499
00500 if(k<rows-1)
00501 m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
00502 if(k<size-1)
00503 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);
00504 }
00505
00506
00507
00508
00509 m_p.setIdentity(rows);
00510 for(Index k = size-1; k >= 0; --k)
00511 m_p.applyTranspositionOnTheRight(k, m_rowsTranspositions.coeff(k));
00512
00513 m_q.setIdentity(cols);
00514 for(Index k = 0; k < size; ++k)
00515 m_q.applyTranspositionOnTheRight(k, m_colsTranspositions.coeff(k));
00516
00517 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00518 return *this;
00519 }
00520
00521 template<typename MatrixType>
00522 typename internal::traits<MatrixType>::Scalar FullPivLU<MatrixType>::determinant() const
00523 {
00524 eigen_assert(m_isInitialized && "LU is not initialized.");
00525 eigen_assert(m_lu.rows() == m_lu.cols() && "You can't take the determinant of a non-square matrix!");
00526 return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
00527 }
00528
00529
00530
00531
00532 template<typename MatrixType>
00533 MatrixType FullPivLU<MatrixType>::reconstructedMatrix() const
00534 {
00535 eigen_assert(m_isInitialized && "LU is not initialized.");
00536 const Index smalldim = std::min(m_lu.rows(), m_lu.cols());
00537
00538 MatrixType res(m_lu.rows(),m_lu.cols());
00539
00540 res = m_lu.leftCols(smalldim)
00541 .template triangularView<UnitLower>().toDenseMatrix()
00542 * m_lu.topRows(smalldim)
00543 .template triangularView<Upper>().toDenseMatrix();
00544
00545
00546 res = m_p.inverse() * res;
00547
00548
00549 res = res * m_q.inverse();
00550
00551 return res;
00552 }
00553
00554
00555
00556 namespace internal {
00557 template<typename _MatrixType>
00558 struct kernel_retval<FullPivLU<_MatrixType> >
00559 : kernel_retval_base<FullPivLU<_MatrixType> >
00560 {
00561 EIGEN_MAKE_KERNEL_HELPERS(FullPivLU<_MatrixType>)
00562
00563 enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
00564 MatrixType::MaxColsAtCompileTime,
00565 MatrixType::MaxRowsAtCompileTime)
00566 };
00567
00568 template<typename Dest> void evalTo(Dest& dst) const
00569 {
00570 const Index cols = dec().matrixLU().cols(), dimker = cols - rank();
00571 if(dimker == 0)
00572 {
00573
00574
00575
00576 dst.setZero();
00577 return;
00578 }
00579
00580
00581
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00592
00593
00594
00595
00596 Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
00597 RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
00598 Index p = 0;
00599 for(Index i = 0; i < dec().nonzeroPivots(); ++i)
00600 if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
00601 pivots.coeffRef(p++) = i;
00602 eigen_internal_assert(p == rank());
00603
00604
00605
00606
00607
00608 Matrix<typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
00609 MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
00610 m(dec().matrixLU().block(0, 0, rank(), cols));
00611 for(Index i = 0; i < rank(); ++i)
00612 {
00613 if(i) m.row(i).head(i).setZero();
00614 m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i);
00615 }
00616 m.block(0, 0, rank(), rank());
00617 m.block(0, 0, rank(), rank()).template triangularView<StrictlyLower>().setZero();
00618 for(Index i = 0; i < rank(); ++i)
00619 m.col(i).swap(m.col(pivots.coeff(i)));
00620
00621
00622
00623
00624 m.topLeftCorner(rank(), rank())
00625 .template triangularView<Upper>().solveInPlace(
00626 m.topRightCorner(rank(), dimker)
00627 );
00628
00629
00630 for(Index i = rank()-1; i >= 0; --i)
00631 m.col(i).swap(m.col(pivots.coeff(i)));
00632
00633
00634 for(Index i = 0; i < rank(); ++i) dst.row(dec().permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker);
00635 for(Index i = rank(); i < cols; ++i) dst.row(dec().permutationQ().indices().coeff(i)).setZero();
00636 for(Index k = 0; k < dimker; ++k) dst.coeffRef(dec().permutationQ().indices().coeff(rank()+k), k) = Scalar(1);
00637 }
00638 };
00639
00640
00641
00642 template<typename _MatrixType>
00643 struct image_retval<FullPivLU<_MatrixType> >
00644 : image_retval_base<FullPivLU<_MatrixType> >
00645 {
00646 EIGEN_MAKE_IMAGE_HELPERS(FullPivLU<_MatrixType>)
00647
00648 enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
00649 MatrixType::MaxColsAtCompileTime,
00650 MatrixType::MaxRowsAtCompileTime)
00651 };
00652
00653 template<typename Dest> void evalTo(Dest& dst) const
00654 {
00655 if(rank() == 0)
00656 {
00657
00658
00659
00660 dst.setZero();
00661 return;
00662 }
00663
00664 Matrix<Index, Dynamic, 1, 0, MaxSmallDimAtCompileTime, 1> pivots(rank());
00665 RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
00666 Index p = 0;
00667 for(Index i = 0; i < dec().nonzeroPivots(); ++i)
00668 if(abs(dec().matrixLU().coeff(i,i)) > premultiplied_threshold)
00669 pivots.coeffRef(p++) = i;
00670 eigen_internal_assert(p == rank());
00671
00672 for(Index i = 0; i < rank(); ++i)
00673 dst.col(i) = originalMatrix().col(dec().permutationQ().indices().coeff(pivots.coeff(i)));
00674 }
00675 };
00676
00677
00678
00679 template<typename _MatrixType, typename Rhs>
00680 struct solve_retval<FullPivLU<_MatrixType>, Rhs>
00681 : solve_retval_base<FullPivLU<_MatrixType>, Rhs>
00682 {
00683 EIGEN_MAKE_SOLVE_HELPERS(FullPivLU<_MatrixType>,Rhs)
00684
00685 template<typename Dest> void evalTo(Dest& dst) const
00686 {
00687
00688
00689
00690
00691
00692
00693
00694
00695 const Index rows = dec().rows(), cols = dec().cols(),
00696 nonzero_pivots = dec().nonzeroPivots();
00697 eigen_assert(rhs().rows() == rows);
00698 const Index smalldim = std::min(rows, cols);
00699
00700 if(nonzero_pivots == 0)
00701 {
00702 dst.setZero();
00703 return;
00704 }
00705
00706 typename Rhs::PlainObject c(rhs().rows(), rhs().cols());
00707
00708
00709 c = dec().permutationP() * rhs();
00710
00711
00712 dec().matrixLU()
00713 .topLeftCorner(smalldim,smalldim)
00714 .template triangularView<UnitLower>()
00715 .solveInPlace(c.topRows(smalldim));
00716 if(rows>cols)
00717 {
00718 c.bottomRows(rows-cols)
00719 -= dec().matrixLU().bottomRows(rows-cols)
00720 * c.topRows(cols);
00721 }
00722
00723
00724 dec().matrixLU()
00725 .topLeftCorner(nonzero_pivots, nonzero_pivots)
00726 .template triangularView<Upper>()
00727 .solveInPlace(c.topRows(nonzero_pivots));
00728
00729
00730 for(Index i = 0; i < nonzero_pivots; ++i)
00731 dst.row(dec().permutationQ().indices().coeff(i)) = c.row(i);
00732 for(Index i = nonzero_pivots; i < dec().matrixLU().cols(); ++i)
00733 dst.row(dec().permutationQ().indices().coeff(i)).setZero();
00734 }
00735 };
00736
00737 }
00738
00739
00740
00741
00742
00743
00744
00745
00746
00747 template<typename Derived>
00748 inline const FullPivLU<typename MatrixBase<Derived>::PlainObject>
00749 MatrixBase<Derived>::fullPivLu() const
00750 {
00751 return FullPivLU<PlainObject>(eval());
00752 }
00753
00754 #endif // EIGEN_LU_H