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
00026 #ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
00027 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
00028
00029 namespace Eigen {
00030
00031 namespace internal {
00032
00033 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType;
00034
00035 template<typename MatrixType>
00036 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
00037 {
00038 typedef typename MatrixType::PlainObject ReturnType;
00039 };
00040
00041 }
00042
00064 template<typename _MatrixType> class FullPivHouseholderQR
00065 {
00066 public:
00067
00068 typedef _MatrixType MatrixType;
00069 enum {
00070 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00071 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00072 Options = MatrixType::Options,
00073 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00074 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00075 };
00076 typedef typename MatrixType::Scalar Scalar;
00077 typedef typename MatrixType::RealScalar RealScalar;
00078 typedef typename MatrixType::Index Index;
00079 typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
00080 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00081 typedef Matrix<Index, 1, ColsAtCompileTime, RowMajor, 1, MaxColsAtCompileTime> IntRowVectorType;
00082 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
00083 typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
00084 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00085 typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
00086
00092 FullPivHouseholderQR()
00093 : m_qr(),
00094 m_hCoeffs(),
00095 m_rows_transpositions(),
00096 m_cols_transpositions(),
00097 m_cols_permutation(),
00098 m_temp(),
00099 m_isInitialized(false),
00100 m_usePrescribedThreshold(false) {}
00101
00108 FullPivHouseholderQR(Index rows, Index cols)
00109 : m_qr(rows, cols),
00110 m_hCoeffs((std::min)(rows,cols)),
00111 m_rows_transpositions(rows),
00112 m_cols_transpositions(cols),
00113 m_cols_permutation(cols),
00114 m_temp((std::min)(rows,cols)),
00115 m_isInitialized(false),
00116 m_usePrescribedThreshold(false) {}
00117
00118 FullPivHouseholderQR(const MatrixType& matrix)
00119 : m_qr(matrix.rows(), matrix.cols()),
00120 m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
00121 m_rows_transpositions(matrix.rows()),
00122 m_cols_transpositions(matrix.cols()),
00123 m_cols_permutation(matrix.cols()),
00124 m_temp((std::min)(matrix.rows(), matrix.cols())),
00125 m_isInitialized(false),
00126 m_usePrescribedThreshold(false)
00127 {
00128 compute(matrix);
00129 }
00130
00148 template<typename Rhs>
00149 inline const internal::solve_retval<FullPivHouseholderQR, Rhs>
00150 solve(const MatrixBase<Rhs>& b) const
00151 {
00152 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00153 return internal::solve_retval<FullPivHouseholderQR, Rhs>(*this, b.derived());
00154 }
00155
00158 MatrixQReturnType matrixQ(void) const;
00159
00162 const MatrixType& matrixQR() const
00163 {
00164 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00165 return m_qr;
00166 }
00167
00168 FullPivHouseholderQR& compute(const MatrixType& matrix);
00169
00170 const PermutationType& colsPermutation() const
00171 {
00172 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00173 return m_cols_permutation;
00174 }
00175
00176 const IntColVectorType& rowsTranspositions() const
00177 {
00178 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00179 return m_rows_transpositions;
00180 }
00181
00195 typename MatrixType::RealScalar absDeterminant() const;
00196
00209 typename MatrixType::RealScalar logAbsDeterminant() const;
00210
00217 inline Index rank() const
00218 {
00219 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00220 RealScalar premultiplied_threshold = internal::abs(m_maxpivot) * threshold();
00221 Index result = 0;
00222 for(Index i = 0; i < m_nonzero_pivots; ++i)
00223 result += (internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold);
00224 return result;
00225 }
00226
00233 inline Index dimensionOfKernel() const
00234 {
00235 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00236 return cols() - rank();
00237 }
00238
00246 inline bool isInjective() const
00247 {
00248 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00249 return rank() == cols();
00250 }
00251
00259 inline bool isSurjective() const
00260 {
00261 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00262 return rank() == rows();
00263 }
00264
00271 inline bool isInvertible() const
00272 {
00273 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00274 return isInjective() && isSurjective();
00275 }
00276 inline const
00282 internal::solve_retval<FullPivHouseholderQR, typename MatrixType::IdentityReturnType>
00283 inverse() const
00284 {
00285 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00286 return internal::solve_retval<FullPivHouseholderQR,typename MatrixType::IdentityReturnType>
00287 (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
00288 }
00289
00290 inline Index rows() const { return m_qr.rows(); }
00291 inline Index cols() const { return m_qr.cols(); }
00292 const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00293
00311 FullPivHouseholderQR& setThreshold(const RealScalar& threshold)
00312 {
00313 m_usePrescribedThreshold = true;
00314 m_prescribedThreshold = threshold;
00315 return *this;
00316 }
00317
00326 FullPivHouseholderQR& setThreshold(Default_t)
00327 {
00328 m_usePrescribedThreshold = false;
00329 return *this;
00330 }
00331
00336 RealScalar threshold() const
00337 {
00338 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00339 return m_usePrescribedThreshold ? m_prescribedThreshold
00340
00341
00342 : NumTraits<Scalar>::epsilon() * m_qr.diagonalSize();
00343 }
00344
00352 inline Index nonzeroPivots() const
00353 {
00354 eigen_assert(m_isInitialized && "LU is not initialized.");
00355 return m_nonzero_pivots;
00356 }
00357
00361 RealScalar maxPivot() const { return m_maxpivot; }
00362
00363 protected:
00364 MatrixType m_qr;
00365 HCoeffsType m_hCoeffs;
00366 IntColVectorType m_rows_transpositions;
00367 IntRowVectorType m_cols_transpositions;
00368 PermutationType m_cols_permutation;
00369 RowVectorType m_temp;
00370 bool m_isInitialized, m_usePrescribedThreshold;
00371 RealScalar m_prescribedThreshold, m_maxpivot;
00372 Index m_nonzero_pivots;
00373 RealScalar m_precision;
00374 Index m_det_pq;
00375 };
00376
00377 template<typename MatrixType>
00378 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
00379 {
00380 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00381 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00382 return internal::abs(m_qr.diagonal().prod());
00383 }
00384
00385 template<typename MatrixType>
00386 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
00387 {
00388 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00389 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00390 return m_qr.diagonal().cwiseAbs().array().log().sum();
00391 }
00392
00393 template<typename MatrixType>
00394 FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00395 {
00396 Index rows = matrix.rows();
00397 Index cols = matrix.cols();
00398 Index size = (std::min)(rows,cols);
00399
00400 m_qr = matrix;
00401 m_hCoeffs.resize(size);
00402
00403 m_temp.resize(cols);
00404
00405 m_precision = NumTraits<Scalar>::epsilon() * size;
00406
00407 m_rows_transpositions.resize(matrix.rows());
00408 m_cols_transpositions.resize(matrix.cols());
00409 Index number_of_transpositions = 0;
00410
00411 RealScalar biggest(0);
00412
00413 m_nonzero_pivots = size;
00414 m_maxpivot = RealScalar(0);
00415
00416 for (Index k = 0; k < size; ++k)
00417 {
00418 Index row_of_biggest_in_corner, col_of_biggest_in_corner;
00419 RealScalar biggest_in_corner;
00420
00421 biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k)
00422 .cwiseAbs()
00423 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
00424 row_of_biggest_in_corner += k;
00425 col_of_biggest_in_corner += k;
00426 if(k==0) biggest = biggest_in_corner;
00427
00428
00429 if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
00430 {
00431 m_nonzero_pivots = k;
00432 for(Index i = k; i < size; i++)
00433 {
00434 m_rows_transpositions.coeffRef(i) = i;
00435 m_cols_transpositions.coeffRef(i) = i;
00436 m_hCoeffs.coeffRef(i) = Scalar(0);
00437 }
00438 break;
00439 }
00440
00441 m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
00442 m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
00443 if(k != row_of_biggest_in_corner) {
00444 m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
00445 ++number_of_transpositions;
00446 }
00447 if(k != col_of_biggest_in_corner) {
00448 m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
00449 ++number_of_transpositions;
00450 }
00451
00452 RealScalar beta;
00453 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
00454 m_qr.coeffRef(k,k) = beta;
00455
00456
00457 if(internal::abs(beta) > m_maxpivot) m_maxpivot = internal::abs(beta);
00458
00459 m_qr.bottomRightCorner(rows-k, cols-k-1)
00460 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
00461 }
00462
00463 m_cols_permutation.setIdentity(cols);
00464 for(Index k = 0; k < size; ++k)
00465 m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
00466
00467 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00468 m_isInitialized = true;
00469
00470 return *this;
00471 }
00472
00473 namespace internal {
00474
00475 template<typename _MatrixType, typename Rhs>
00476 struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
00477 : solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs>
00478 {
00479 EIGEN_MAKE_SOLVE_HELPERS(FullPivHouseholderQR<_MatrixType>,Rhs)
00480
00481 template<typename Dest> void evalTo(Dest& dst) const
00482 {
00483 const Index rows = dec().rows(), cols = dec().cols();
00484 eigen_assert(rhs().rows() == rows);
00485
00486
00487
00488 if(dec().rank()==0)
00489 {
00490 dst.setZero();
00491 return;
00492 }
00493
00494 typename Rhs::PlainObject c(rhs());
00495
00496 Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs().cols());
00497 for (Index k = 0; k < dec().rank(); ++k)
00498 {
00499 Index remainingSize = rows-k;
00500 c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k)));
00501 c.bottomRightCorner(remainingSize, rhs().cols())
00502 .applyHouseholderOnTheLeft(dec().matrixQR().col(k).tail(remainingSize-1),
00503 dec().hCoeffs().coeff(k), &temp.coeffRef(0));
00504 }
00505
00506 if(!dec().isSurjective())
00507 {
00508
00509 RealScalar biggest_in_upper_part_of_c = c.topRows( dec().rank() ).cwiseAbs().maxCoeff();
00510 RealScalar biggest_in_lower_part_of_c = c.bottomRows(rows-dec().rank()).cwiseAbs().maxCoeff();
00511
00512 const RealScalar m_precision = NumTraits<Scalar>::epsilon() * (std::min)(rows,cols);
00513
00514 if(!internal::isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
00515 return;
00516 }
00517 dec().matrixQR()
00518 .topLeftCorner(dec().rank(), dec().rank())
00519 .template triangularView<Upper>()
00520 .solveInPlace(c.topRows(dec().rank()));
00521
00522 for(Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
00523 for(Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
00524 }
00525 };
00526
00533 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
00534 : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
00535 {
00536 public:
00537 typedef typename MatrixType::Index Index;
00538 typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
00539 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00540 typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
00541 MatrixType::MaxRowsAtCompileTime> WorkVectorType;
00542
00543 FullPivHouseholderQRMatrixQReturnType(const MatrixType& qr,
00544 const HCoeffsType& hCoeffs,
00545 const IntColVectorType& rowsTranspositions)
00546 : m_qr(qr),
00547 m_hCoeffs(hCoeffs),
00548 m_rowsTranspositions(rowsTranspositions)
00549 {}
00550
00551 template <typename ResultType>
00552 void evalTo(ResultType& result) const
00553 {
00554 const Index rows = m_qr.rows();
00555 WorkVectorType workspace(rows);
00556 evalTo(result, workspace);
00557 }
00558
00559 template <typename ResultType>
00560 void evalTo(ResultType& result, WorkVectorType& workspace) const
00561 {
00562
00563
00564
00565 const Index rows = m_qr.rows();
00566 const Index cols = m_qr.cols();
00567 const Index size = (std::min)(rows, cols);
00568 workspace.resize(rows);
00569 result.setIdentity(rows, rows);
00570 for (Index k = size-1; k >= 0; k--)
00571 {
00572 result.block(k, k, rows-k, rows-k)
00573 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), internal::conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
00574 result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
00575 }
00576 }
00577
00578 Index rows() const { return m_qr.rows(); }
00579 Index cols() const { return m_qr.rows(); }
00580
00581 protected:
00582 typename MatrixType::Nested m_qr;
00583 typename HCoeffsType::Nested m_hCoeffs;
00584 typename IntColVectorType::Nested m_rowsTranspositions;
00585 };
00586
00587 }
00588
00589 template<typename MatrixType>
00590 inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const
00591 {
00592 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
00593 return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
00594 }
00595
00600 template<typename Derived>
00601 const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
00602 MatrixBase<Derived>::fullPivHouseholderQr() const
00603 {
00604 return FullPivHouseholderQR<PlainObject>(eval());
00605 }
00606
00607 }
00608
00609 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H