[PATCH] store rank of original QR decomposition in CompleteOrthogonalDecomposition, so that if the additional householder transformations change the rank (as measured by QR's rank()) method, we still stick with the initially computed rank in the CODs pseudoInverse() method when computing the pseudoinverse

[Simon Rutishauser <simon.rutishauser@xxxxxxxxx>]

---
 Eigen/src/QR/CompleteOrthogonalDecomposition.h | 23 +++++++++++++++--------
 1 file changed, 15 insertions(+), 8 deletions(-)

diff --git a/Eigen/src/QR/CompleteOrthogonalDecomposition.h b/Eigen/src/QR/CompleteOrthogonalDecomposition.h
index 230d0d2..75abd6f 100644
--- a/Eigen/src/QR/CompleteOrthogonalDecomposition.h
+++ b/Eigen/src/QR/CompleteOrthogonalDecomposition.h
@@ -83,7 +83,7 @@ class CompleteOrthogonalDecomposition {
    * perform decompositions via
    * \c CompleteOrthogonalDecomposition::compute(const* MatrixType&).
    */
-  CompleteOrthogonalDecomposition() : m_cpqr(), m_zCoeffs(), m_temp() {}
+  CompleteOrthogonalDecomposition() : m_cpqr(), m_zCoeffs(), m_temp(), m_rank(0) {}
 
   /** \brief Default Constructor with memory preallocation
    *
@@ -92,7 +92,7 @@ class CompleteOrthogonalDecomposition {
    * \sa CompleteOrthogonalDecomposition()
    */
   CompleteOrthogonalDecomposition(Index rows, Index cols)
-      : m_cpqr(rows, cols), m_zCoeffs((std::min)(rows, cols)), m_temp(cols) {}
+      : m_cpqr(rows, cols), m_zCoeffs((std::min)(rows, cols)), m_temp(cols), m_rank(0) {}
 
   /** \brief Constructs a complete orthogonal decomposition from a given
    * matrix.
@@ -114,7 +114,8 @@ class CompleteOrthogonalDecomposition {
   explicit CompleteOrthogonalDecomposition(const EigenBase<InputType>& matrix)
       : m_cpqr(matrix.rows(), matrix.cols()),
         m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
-        m_temp(matrix.cols()) {
+        m_temp(matrix.cols()),
+        m_rank(0) {
     compute(matrix.derived());
   }
 
@@ -209,7 +210,11 @@ class CompleteOrthogonalDecomposition {
    * nonzero. For that, it uses the threshold value that you can control by
    * calling setThreshold(const RealScalar&).
    */
-  inline Index rank() const { return m_cpqr.rank(); }
+  inline Index rank() const
+  {
+    eigen_assert(m_cpqr.m_isInitialized && "Decomposition is not initialized.");
+    return m_rank;
+  }
 
   /** \returns the dimension of the kernel of the matrix of which *this is the
    * complete orthogonal decomposition.
@@ -227,7 +232,7 @@ class CompleteOrthogonalDecomposition {
    * nonzero. For that, it uses the threshold value that you can control by
    * calling setThreshold(const RealScalar&).
    */
-  inline bool isInjective() const { return m_cpqr.isInjective(); }
+  inline bool isInjective() const { return rank() == m_cpqr.cols(); }
 
   /** \returns true if the matrix of which *this is the decomposition represents
    * a surjective linear map; false otherwise.
@@ -236,7 +241,7 @@ class CompleteOrthogonalDecomposition {
    * nonzero. For that, it uses the threshold value that you can control by
    * calling setThreshold(const RealScalar&).
    */
-  inline bool isSurjective() const { return m_cpqr.isSurjective(); }
+  inline bool isSurjective() const { return rank() == m_cpqr.rows(); }
 
   /** \returns true if the matrix of which *this is the complete orthogonal
    * decomposition is invertible.
@@ -245,7 +250,7 @@ class CompleteOrthogonalDecomposition {
    * nonzero. For that, it uses the threshold value that you can control by
    * calling setThreshold(const RealScalar&).
    */
-  inline bool isInvertible() const { return m_cpqr.isInvertible(); }
+  inline bool isInvertible() const { return isInjective() && isSurjective(); }
 
   /** \returns the pseudo-inverse of the matrix of which *this is the complete
    * orthogonal decomposition.
@@ -362,6 +367,7 @@ class CompleteOrthogonalDecomposition {
   ColPivHouseholderQR<MatrixType> m_cpqr;
   HCoeffsType m_zCoeffs;
   RowVectorType m_temp;
+  Index m_rank;
 };
 
 template <typename MatrixType>
@@ -395,7 +401,8 @@ CompleteOrthogonalDecomposition<MatrixType>& CompleteOrthogonalDecomposition<
   // Compute the column pivoted QR factorization A P = Q R.
   m_cpqr.compute(matrix);
 
-  const Index rank = m_cpqr.rank();
+  m_rank = m_cpqr.rank();
+  const Index rank = m_rank;
   const Index cols = matrix.cols();
   const Index rows = matrix.rows();
   m_zCoeffs.resize((std::min)(rows, cols));
-- 
2.8.1


--nextPart1997099.Ihaho3VsFr--