diff --git a/Eigen/LU b/Eigen/LU index 64dcdee60..d46bacf4b 100644 --- a/Eigen/LU +++ b/Eigen/LU @@ -25,6 +25,7 @@ #include "src/misc/Kernel.h" #include "src/misc/Image.h" +#include "src/misc/RankRevealingBase.h" // IWYU pragma: begin_exports #include "src/LU/FullPivLU.h" diff --git a/Eigen/QR b/Eigen/QR index c38b453b0..78c05632e 100644 --- a/Eigen/QR +++ b/Eigen/QR @@ -31,6 +31,8 @@ * \endcode */ +#include "src/misc/RankRevealingBase.h" + // IWYU pragma: begin_exports #include "src/QR/HouseholderQR.h" #include "src/QR/FullPivHouseholderQR.h" diff --git a/Eigen/src/LU/FullPivLU.h b/Eigen/src/LU/FullPivLU.h index 74bdcf34c..098b16fe2 100644 --- a/Eigen/src/LU/FullPivLU.h +++ b/Eigen/src/LU/FullPivLU.h @@ -60,11 +60,23 @@ struct traits > : traits * \sa MatrixBase::fullPivLu(), MatrixBase::determinant(), MatrixBase::inverse() */ template -class FullPivLU : public SolverBase > { +class FullPivLU : public SolverBase >, + public RankRevealingBase > { public: typedef MatrixType_ MatrixType; typedef SolverBase Base; + typedef RankRevealingBase RankRevealingBase_; friend class SolverBase; + friend class RankRevealingBase; + using RankRevealingBase_::dimensionOfKernel; + using RankRevealingBase_::isInjective; + using RankRevealingBase_::isInvertible; + using RankRevealingBase_::isSurjective; + using RankRevealingBase_::maxPivot; + using RankRevealingBase_::nonzeroPivots; + using RankRevealingBase_::rank; + using RankRevealingBase_::setThreshold; + using RankRevealingBase_::threshold; EIGEN_GENERIC_PUBLIC_INTERFACE(FullPivLU) enum { @@ -148,23 +160,6 @@ class FullPivLU : public SolverBase > return m_lu; } - /** \returns the number of nonzero pivots in the LU decomposition. - * Here nonzero is meant in the exact sense, not in a fuzzy sense. - * So that notion isn't really intrinsically interesting, but it is - * still useful when implementing algorithms. - * - * \sa rank() - */ - inline Index nonzeroPivots() const { - eigen_assert(m_isInitialized && "LU is not initialized."); - return m_nonzero_pivots; - } - - /** \returns the absolute value of the biggest pivot, i.e. the biggest - * diagonal coefficient of U. - */ - RealScalar maxPivot() const { return m_maxpivot; } - /** \returns the permutation matrix P * * \sa permutationQ() @@ -278,114 +273,10 @@ class FullPivLU : public SolverBase > */ typename internal::traits::Scalar determinant() const; - /** Allows to prescribe a threshold to be used by certain methods, such as rank(), - * who need to determine when pivots are to be considered nonzero. This is not used for the - * LU decomposition itself. - * - * When it needs to get the threshold value, Eigen calls threshold(). By default, this - * uses a formula to automatically determine a reasonable threshold. - * Once you have called the present method setThreshold(const RealScalar&), - * your value is used instead. - * - * \param threshold The new value to use as the threshold. - * - * A pivot will be considered nonzero if its absolute value is strictly greater than - * \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$ - * where maxpivot is the biggest pivot. - * - * If you want to come back to the default behavior, call setThreshold(Default_t) - */ - FullPivLU& setThreshold(const RealScalar& threshold) { - m_usePrescribedThreshold = true; - m_prescribedThreshold = threshold; - return *this; - } - - /** Allows to come back to the default behavior, letting Eigen use its default formula for - * determining the threshold. - * - * You should pass the special object Eigen::Default as parameter here. - * \code lu.setThreshold(Eigen::Default); \endcode - * - * See the documentation of setThreshold(const RealScalar&). - */ - FullPivLU& setThreshold(Default_t) { - m_usePrescribedThreshold = false; - return *this; - } - - /** Returns the threshold that will be used by certain methods such as rank(). - * - * See the documentation of setThreshold(const RealScalar&). - */ - RealScalar threshold() const { - eigen_assert(m_isInitialized || m_usePrescribedThreshold); - return m_usePrescribedThreshold ? m_prescribedThreshold - // Higham's backward error bound for Gaussian elimination with - // complete pivoting (Theorem 9.4) is ||ΔA||₂ ≤ c·min(m,n)·u·||A||₂. - // The factor of 4 covers the constant c. - : NumTraits::epsilon() * RealScalar(4 * m_lu.diagonalSize()); - } - - /** \returns the rank of the matrix of which *this is the LU decomposition. - * - * \note This method has to determine which pivots should be considered nonzero. - * For that, it uses the threshold value that you can control by calling - * setThreshold(const RealScalar&). - */ - inline Index rank() const { + /** \returns the absolute value of the i-th pivot coefficient (for RankRevealingBase). */ + RealScalar pivotCoeff(Index i) const { using std::abs; - eigen_assert(m_isInitialized && "LU is not initialized."); - RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); - Index result = 0; - for (Index i = 0; i < m_nonzero_pivots; ++i) result += (abs(m_lu.coeff(i, i)) > premultiplied_threshold); - return result; - } - - /** \returns the dimension of the kernel of the matrix of which *this is the LU decomposition. - * - * \note This method has to determine which pivots should be considered nonzero. - * For that, it uses the threshold value that you can control by calling - * setThreshold(const RealScalar&). - */ - inline Index dimensionOfKernel() const { - eigen_assert(m_isInitialized && "LU is not initialized."); - return cols() - rank(); - } - - /** \returns true if the matrix of which *this is the LU decomposition represents an injective - * linear map, i.e. has trivial kernel; false otherwise. - * - * \note This method has to determine which pivots should be considered nonzero. - * For that, it uses the threshold value that you can control by calling - * setThreshold(const RealScalar&). - */ - inline bool isInjective() const { - eigen_assert(m_isInitialized && "LU is not initialized."); - return rank() == cols(); - } - - /** \returns true if the matrix of which *this is the LU decomposition represents a surjective - * linear map; false otherwise. - * - * \note This method has to determine which pivots should be considered nonzero. - * For that, it uses the threshold value that you can control by calling - * setThreshold(const RealScalar&). - */ - inline bool isSurjective() const { - eigen_assert(m_isInitialized && "LU is not initialized."); - return rank() == rows(); - } - - /** \returns true if the matrix of which *this is the LU decomposition is invertible. - * - * \note This method has to determine which pivots should be considered nonzero. - * For that, it uses the threshold value that you can control by calling - * setThreshold(const RealScalar&). - */ - inline bool isInvertible() const { - eigen_assert(m_isInitialized && "LU is not initialized."); - return isInjective() && (m_lu.rows() == m_lu.cols()); + return abs(m_lu.coeff(i, i)); } /** \returns the inverse of the matrix of which *this is the LU decomposition. @@ -424,15 +315,13 @@ class FullPivLU : public SolverBase > PermutationQType m_q; IntColVectorType m_rowsTranspositions; IntRowVectorType m_colsTranspositions; - Index m_nonzero_pivots; RealScalar m_l1_norm; - RealScalar m_maxpivot, m_prescribedThreshold; signed char m_det_pq; - bool m_isInitialized, m_usePrescribedThreshold; + bool m_isInitialized; }; template -FullPivLU::FullPivLU() : m_isInitialized(false), m_usePrescribedThreshold(false) {} +FullPivLU::FullPivLU() : m_isInitialized(false) {} template FullPivLU::FullPivLU(Index rows, Index cols) @@ -441,8 +330,7 @@ FullPivLU::FullPivLU(Index rows, Index cols) m_q(cols), m_rowsTranspositions(rows), m_colsTranspositions(cols), - m_isInitialized(false), - m_usePrescribedThreshold(false) {} + m_isInitialized(false) {} template template @@ -452,8 +340,7 @@ FullPivLU::FullPivLU(const EigenBase& m m_q(matrix.cols()), m_rowsTranspositions(matrix.rows()), m_colsTranspositions(matrix.cols()), - m_isInitialized(false), - m_usePrescribedThreshold(false) { + m_isInitialized(false) { compute(matrix.derived()); } @@ -465,8 +352,7 @@ FullPivLU::FullPivLU(EigenBase& matrix) m_q(matrix.cols()), m_rowsTranspositions(matrix.rows()), m_colsTranspositions(matrix.cols()), - m_isInitialized(false), - m_usePrescribedThreshold(false) { + m_isInitialized(false) { computeInPlace(); } @@ -487,8 +373,8 @@ void FullPivLU::computeInPlace() { m_colsTranspositions.resize(m_lu.cols()); Index number_of_transpositions = 0; // number of NONTRIVIAL transpositions, i.e. m_rowsTranspositions[i]!=i - m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) - m_maxpivot = RealScalar(0); + this->m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) + this->m_maxpivot = RealScalar(0); for (Index k = 0; k < size; ++k) { // First, we need to find the pivot. @@ -507,7 +393,7 @@ void FullPivLU::computeInPlace() { if (numext::is_exactly_zero(biggest_in_corner)) { // before exiting, make sure to initialize the still uninitialized transpositions // in a sane state without destroying what we already have. - m_nonzero_pivots = k; + this->m_nonzero_pivots = k; for (Index i = k; i < size; ++i) { m_rowsTranspositions.coeffRef(i) = internal::convert_index(i); m_colsTranspositions.coeffRef(i) = internal::convert_index(i); @@ -517,7 +403,7 @@ void FullPivLU::computeInPlace() { RealScalar abs_pivot = internal::abs_knowing_score()( m_lu(row_of_biggest_in_corner, col_of_biggest_in_corner), biggest_in_corner); - if (abs_pivot > m_maxpivot) m_maxpivot = abs_pivot; + if (abs_pivot > this->m_maxpivot) this->m_maxpivot = abs_pivot; // Now that we've found the pivot, we need to apply the row/col swaps to // bring it to the location (k,k). diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h index 369507851..738656473 100644 --- a/Eigen/src/QR/ColPivHouseholderQR.h +++ b/Eigen/src/QR/ColPivHouseholderQR.h @@ -51,11 +51,23 @@ struct traits> : traits -class ColPivHouseholderQR : public SolverBase> { +class ColPivHouseholderQR : public SolverBase>, + public RankRevealingBase> { public: typedef MatrixType_ MatrixType; typedef SolverBase Base; + typedef RankRevealingBase RankRevealingBase_; friend class SolverBase; + friend class RankRevealingBase; + using RankRevealingBase_::dimensionOfKernel; + using RankRevealingBase_::isInjective; + using RankRevealingBase_::isInvertible; + using RankRevealingBase_::isSurjective; + using RankRevealingBase_::maxPivot; + using RankRevealingBase_::nonzeroPivots; + using RankRevealingBase_::rank; + using RankRevealingBase_::setThreshold; + using RankRevealingBase_::threshold; typedef PermutationIndex_ PermutationIndex; EIGEN_GENERIC_PUBLIC_INTERFACE(ColPivHouseholderQR) @@ -82,7 +94,6 @@ class ColPivHouseholderQR : public SolverBase premultiplied_threshold); - return result; - } - - /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition. - * - * \note This method has to determine which pivots should be considered nonzero. - * For that, it uses the threshold value that you can control by calling - * setThreshold(const RealScalar&). - */ - inline Index dimensionOfKernel() const { - eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - return cols() - rank(); - } - - /** \returns true if the matrix of which *this is the QR decomposition represents an injective - * linear map, i.e. has trivial kernel; false otherwise. - * - * \note This method has to determine which pivots should be considered nonzero. - * For that, it uses the threshold value that you can control by calling - * setThreshold(const RealScalar&). - */ - inline bool isInjective() const { - eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - return rank() == cols(); - } - - /** \returns true if the matrix of which *this is the QR decomposition represents a surjective - * linear map; false otherwise. - * - * \note This method has to determine which pivots should be considered nonzero. - * For that, it uses the threshold value that you can control by calling - * setThreshold(const RealScalar&). - */ - inline bool isSurjective() const { - eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - return rank() == rows(); - } - - /** \returns true if the matrix of which *this is the QR decomposition is invertible. - * - * \note This method has to determine which pivots should be considered nonzero. - * For that, it uses the threshold value that you can control by calling - * setThreshold(const RealScalar&). - */ - inline bool isInvertible() const { - eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - return isInjective() && isSurjective(); + return abs(m_qr.coeff(i, i)); } /** \returns the inverse of the matrix of which *this is the QR decomposition. @@ -332,72 +287,6 @@ class ColPivHouseholderQR : public SolverBase::epsilon() * RealScalar(4 * m_qr.diagonalSize()); - } - - /** \returns the number of nonzero pivots in the QR decomposition. - * Here nonzero is meant in the exact sense, not in a fuzzy sense. - * So that notion isn't really intrinsically interesting, but it is - * still useful when implementing algorithms. - * - * \sa rank() - */ - inline Index nonzeroPivots() const { - eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); - return m_nonzero_pivots; - } - - /** \returns the absolute value of the biggest pivot, i.e. the biggest - * diagonal coefficient of R. - */ - RealScalar maxPivot() const { return m_maxpivot; } - /** \brief Reports whether the QR factorization was successful. * * \note This function always returns \c Success. It is provided for compatibility @@ -431,9 +320,7 @@ class ColPivHouseholderQR : public SolverBase::computeInPlace() { numext::abs2(m_colNormsUpdated.maxCoeff() * NumTraits::epsilon()) / RealScalar(rows); RealScalar norm_downdate_threshold = numext::sqrt(NumTraits::epsilon()); - m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) - m_maxpivot = RealScalar(0); + this->m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) + this->m_maxpivot = RealScalar(0); for (Index k = 0; k < size; ++k) { // first, we look up in our table m_colNormsUpdated which column has the biggest norm @@ -526,7 +413,8 @@ void ColPivHouseholderQR::computeInPlace() { // Track the number of meaningful pivots but do not stop the decomposition to make // sure that the initial matrix is properly reproduced. See bug 941. - if (m_nonzero_pivots == size && biggest_col_sq_norm < threshold_helper * RealScalar(rows - k)) m_nonzero_pivots = k; + if (this->m_nonzero_pivots == size && biggest_col_sq_norm < threshold_helper * RealScalar(rows - k)) + this->m_nonzero_pivots = k; // apply the transposition to the columns m_colsTranspositions.coeffRef(k) = static_cast(biggest_col_index); @@ -545,7 +433,7 @@ void ColPivHouseholderQR::computeInPlace() { m_qr.coeffRef(k, k) = beta; // remember the maximum absolute value of diagonal coefficients - if (abs(beta) > m_maxpivot) m_maxpivot = abs(beta); + if (abs(beta) > this->m_maxpivot) this->m_maxpivot = abs(beta); // apply the householder transformation m_qr.bottomRightCorner(rows - k, cols - k - 1) diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h index 40f9047f5..547bad087 100644 --- a/Eigen/src/QR/FullPivHouseholderQR.h +++ b/Eigen/src/QR/FullPivHouseholderQR.h @@ -60,11 +60,23 @@ struct traits -class FullPivHouseholderQR : public SolverBase > { +class FullPivHouseholderQR : public SolverBase >, + public RankRevealingBase > { public: typedef MatrixType_ MatrixType; typedef SolverBase Base; + typedef RankRevealingBase RankRevealingBase_; friend class SolverBase; + friend class RankRevealingBase; + using RankRevealingBase_::dimensionOfKernel; + using RankRevealingBase_::isInjective; + using RankRevealingBase_::isInvertible; + using RankRevealingBase_::isSurjective; + using RankRevealingBase_::maxPivot; + using RankRevealingBase_::nonzeroPivots; + using RankRevealingBase_::rank; + using RankRevealingBase_::setThreshold; + using RankRevealingBase_::threshold; typedef PermutationIndex_ PermutationIndex; EIGEN_GENERIC_PUBLIC_INTERFACE(FullPivHouseholderQR) @@ -105,8 +117,7 @@ class FullPivHouseholderQR : public SolverBase premultiplied_threshold); - return result; - } - - /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition. - * - * \note This method has to determine which pivots should be considered nonzero. - * For that, it uses the threshold value that you can control by calling - * setThreshold(const RealScalar&). - */ - inline Index dimensionOfKernel() const { - eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); - return cols() - rank(); - } - - /** \returns true if the matrix of which *this is the QR decomposition represents an injective - * linear map, i.e. has trivial kernel; false otherwise. - * - * \note This method has to determine which pivots should be considered nonzero. - * For that, it uses the threshold value that you can control by calling - * setThreshold(const RealScalar&). - */ - inline bool isInjective() const { - eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); - return rank() == cols(); - } - - /** \returns true if the matrix of which *this is the QR decomposition represents a surjective - * linear map; false otherwise. - * - * \note This method has to determine which pivots should be considered nonzero. - * For that, it uses the threshold value that you can control by calling - * setThreshold(const RealScalar&). - */ - inline bool isSurjective() const { - eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); - return rank() == rows(); - } - - /** \returns true if the matrix of which *this is the QR decomposition is invertible. - * - * \note This method has to determine which pivots should be considered nonzero. - * For that, it uses the threshold value that you can control by calling - * setThreshold(const RealScalar&). - */ - inline bool isInvertible() const { - eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized."); - return isInjective() && isSurjective(); + return abs(m_qr.coeff(i, i)); } /** \returns the inverse of the matrix of which *this is the QR decomposition. @@ -353,72 +306,6 @@ class FullPivHouseholderQR : public SolverBase::epsilon() * RealScalar(4 * m_qr.diagonalSize()); - } - - /** \returns the number of nonzero pivots in the QR decomposition. - * Here nonzero is meant in the exact sense, not in a fuzzy sense. - * So that notion isn't really intrinsically interesting, but it is - * still useful when implementing algorithms. - * - * \sa rank() - */ - inline Index nonzeroPivots() const { - eigen_assert(m_isInitialized && "LU is not initialized."); - return m_nonzero_pivots; - } - - /** \returns the absolute value of the biggest pivot, i.e. the biggest - * diagonal coefficient of U. - */ - RealScalar maxPivot() const { return m_maxpivot; } - #ifndef EIGEN_PARSED_BY_DOXYGEN template void _solve_impl(const RhsType& rhs, DstType& dst) const; @@ -438,9 +325,7 @@ class FullPivHouseholderQR : public SolverBase::computeInPlace() { RealScalar biggest(0); - m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) - m_maxpivot = RealScalar(0); + this->m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) + this->m_maxpivot = RealScalar(0); for (Index k = 0; k < size; ++k) { Index row_of_biggest_in_corner, col_of_biggest_in_corner; @@ -532,7 +417,7 @@ void FullPivHouseholderQR::computeInPlace() { // if the corner is negligible, then we have less than full rank, and we can finish early if (internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision)) { - m_nonzero_pivots = k; + this->m_nonzero_pivots = k; for (Index i = k; i < size; i++) { m_rows_transpositions.coeffRef(i) = internal::convert_index(i); m_cols_transpositions.coeffRef(i) = internal::convert_index(i); @@ -557,7 +442,7 @@ void FullPivHouseholderQR::computeInPlace() { m_qr.coeffRef(k, k) = beta; // remember the maximum absolute value of diagonal coefficients - if (abs(beta) > m_maxpivot) m_maxpivot = abs(beta); + if (abs(beta) > this->m_maxpivot) this->m_maxpivot = abs(beta); m_qr.bottomRightCorner(rows - k, cols - k - 1) .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows - k - 1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k + 1)); diff --git a/Eigen/src/misc/RankRevealingBase.h b/Eigen/src/misc/RankRevealingBase.h new file mode 100644 index 000000000..1729c912f --- /dev/null +++ b/Eigen/src/misc/RankRevealingBase.h @@ -0,0 +1,178 @@ +// This file is part of Eigen, a lightweight C++ template library +// for linear algebra. +// +// This Source Code Form is subject to the terms of the Mozilla +// Public License v. 2.0. If a copy of the MPL was not distributed +// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#ifndef EIGEN_RANK_REVEALING_BASE_H +#define EIGEN_RANK_REVEALING_BASE_H + +// IWYU pragma: private +#include "./InternalHeaderCheck.h" + +namespace Eigen { + +/** \brief CRTP mixin providing threshold management, rank computation, and rank-derived queries + * for rank-revealing decompositions (FullPivLU, ColPivHouseholderQR, FullPivHouseholderQR). + * + * \tparam Derived the concrete decomposition class (CRTP parameter) + * + * The derived class must provide: + * - rows(), cols() (inherited from SolverBase) + * - m_isInitialized (bool member, also used by SolverBase) + * - pivotCoeff(Index i) returning the absolute value of the i-th pivot + */ +template +class RankRevealingBase { + public: + typedef typename internal::traits::Scalar Scalar; + typedef typename NumTraits::Real RealScalar; + + RankRevealingBase() + : m_usePrescribedThreshold(false), + m_prescribedThreshold(RealScalar(0)), + m_maxpivot(RealScalar(0)), + m_nonzero_pivots(0) {} + + /** Allows to prescribe a threshold to be used by certain methods, such as rank(), + * who need to determine when pivots are to be considered nonzero. This is not used for the + * decomposition itself. + * + * When it needs to get the threshold value, Eigen calls threshold(). By default, this + * uses a formula to automatically determine a reasonable threshold. + * Once you have called the present method setThreshold(const RealScalar&), + * your value is used instead. + * + * \param threshold The new value to use as the threshold. + * + * A pivot will be considered nonzero if its absolute value is strictly greater than + * \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$ + * where maxpivot is the biggest pivot. + * + * If you want to come back to the default behavior, call setThreshold(Default_t) + */ + Derived& setThreshold(const RealScalar& threshold) { + m_usePrescribedThreshold = true; + m_prescribedThreshold = threshold; + return self(); + } + + /** Allows to come back to the default behavior, letting Eigen use its default formula for + * determining the threshold. + * + * You should pass the special object Eigen::Default as parameter here. + * \code dec.setThreshold(Eigen::Default); \endcode + * + * See the documentation of setThreshold(const RealScalar&). + */ + Derived& setThreshold(Default_t) { + m_usePrescribedThreshold = false; + return self(); + } + + /** Returns the threshold that will be used by certain methods such as rank(). + * + * See the documentation of setThreshold(const RealScalar&). + */ + RealScalar threshold() const { + eigen_assert(self().m_isInitialized || m_usePrescribedThreshold); + // Higham's backward error bound: ||ΔA||₂ ≤ c·min(m,n)·u·||A||₂. + // The factor of 4 covers the constant c. + return m_usePrescribedThreshold + ? m_prescribedThreshold + : NumTraits::epsilon() * RealScalar(4 * (std::min)(self().rows(), self().cols())); + } + + /** \returns the rank of the matrix of which *this is the decomposition. + * + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). + */ + inline Index rank() const { + using std::abs; + eigen_assert(self().m_isInitialized && "Decomposition is not initialized."); + RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); + Index result = 0; + for (Index i = 0; i < m_nonzero_pivots; ++i) result += (self().pivotCoeff(i) > premultiplied_threshold); + return result; + } + + /** \returns the dimension of the kernel of the matrix of which *this is the decomposition. + * + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). + */ + inline Index dimensionOfKernel() const { + eigen_assert(self().m_isInitialized && "Decomposition is not initialized."); + return self().cols() - rank(); + } + + /** \returns true if the matrix of which *this is the decomposition represents an injective + * linear map, i.e. has trivial kernel; false otherwise. + * + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). + */ + inline bool isInjective() const { + eigen_assert(self().m_isInitialized && "Decomposition is not initialized."); + return rank() == self().cols(); + } + + /** \returns true if the matrix of which *this is the decomposition represents a surjective + * linear map; false otherwise. + * + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). + */ + inline bool isSurjective() const { + eigen_assert(self().m_isInitialized && "Decomposition is not initialized."); + return rank() == self().rows(); + } + + /** \returns true if the matrix of which *this is the decomposition is invertible. + * + * \note This method has to determine which pivots should be considered nonzero. + * For that, it uses the threshold value that you can control by calling + * setThreshold(const RealScalar&). + */ + inline bool isInvertible() const { + eigen_assert(self().m_isInitialized && "Decomposition is not initialized."); + return isInjective() && isSurjective(); + } + + /** \returns the number of nonzero pivots in the decomposition. + * Here nonzero is meant in the exact sense, not in a fuzzy sense. + * So that notion isn't really intrinsically interesting, but it is + * still useful when implementing algorithms. + * + * \sa rank() + */ + inline Index nonzeroPivots() const { + eigen_assert(self().m_isInitialized && "Decomposition is not initialized."); + return m_nonzero_pivots; + } + + /** \returns the absolute value of the biggest pivot, i.e. the biggest + * diagonal coefficient of U (or R). + */ + RealScalar maxPivot() const { return m_maxpivot; } + + protected: + bool m_usePrescribedThreshold; + RealScalar m_prescribedThreshold; + RealScalar m_maxpivot; + Index m_nonzero_pivots; + + private: + Derived& self() { return static_cast(*this); } + const Derived& self() const { return static_cast(*this); } +}; + +} // end namespace Eigen + +#endif // EIGEN_RANK_REVEALING_BASE_H