QR: add isInjective(), isSurjective(),

    mark isFullRank() deprecated,
    add solve() (mix of Keir's patch and LU::solve())
=> there is big problem with complex which are not working

1 files changed, 119 insertions, 11 deletions

diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h
index 13d49a4a3..0ea839d25 100644
--- a/Eigen/src/QR/QR.h
+++ b/Eigen/src/QR/QR.h
@@ -55,12 +55,64 @@ template<typename MatrixType> class QR
     {
       _compute(matrix);
     }
-
-    /** \returns whether or not the matrix is of full rank */
-    bool isFullRank() const { return rank() == std::min(m_qr.rows(),m_qr.cols()); }
     
+    /** \deprecated use isInjective()
+      * \returns whether or not the matrix is of full rank
+      *
+      * \note Since the rank is computed only once, i.e. the first time it is needed, this
+      *       method almost does not perform any further computation.
+      */
+    bool isFullRank() const EIGEN_DEPRECATED { return rank() == m_qr.cols(); }
+    
+    /** \returns the rank of the matrix of which *this is the QR decomposition.
+      *
+      * \note Since the rank is computed only once, i.e. the first time it is needed, this
+      *       method almost does not perform any further computation.
+      */
     int rank() const;
+    
+    /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition.
+      *
+      * \note Since the rank is computed only once, i.e. the first time it is needed, this
+      *       method almost does not perform any further computation.
+      */
+    inline int dimensionOfKernel() const
+    {
+      return m_qr.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 Since the rank is computed only once, i.e. the first time it is needed, this
+      *       method almost does not perform any further computation.
+      */
+    inline bool isInjective() const
+    {
+      return rank() == m_qr.cols();
+    }
+    
+    /** \returns true if the matrix of which *this is the QR decomposition represents a surjective
+      *          linear map; false otherwise.
+      *
+      * \note Since the rank is computed only once, i.e. the first time it is needed, this
+      *       method almost does not perform any further computation.
+      */
+    inline bool isSurjective() const
+    {
+      return rank() == m_qr.rows();
+    }
 
+    /** \returns true if the matrix of which *this is the QR decomposition is invertible.
+      *
+      * \note Since the rank is computed only once, i.e. the first time it is needed, this
+      *       method almost does not perform any further computation.
+      */
+    inline bool isInvertible() const
+    {
+      return isInjective() && isSurjective();
+    }
+    
     /** \returns a read-only expression of the matrix R of the actual the QR decomposition */
     const Part<NestByValue<MatrixRBlockType>, UpperTriangular>
     matrixR(void) const
@@ -69,6 +121,32 @@ template<typename MatrixType> class QR
       return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template part<UpperTriangular>();
     }
 
+    /** This method finds a solution x to the equation Ax=b, where A is the matrix of which
+      * *this is the QR decomposition, if any exists.
+      *
+      * \param b the right-hand-side of the equation to solve.
+      *
+      * \param result a pointer to the vector/matrix in which to store the solution, if any exists.
+      *          Resized if necessary, so that result->rows()==A.cols() and result->cols()==b.cols().
+      *          If no solution exists, *result is left with undefined coefficients.
+      *
+      * \returns true if any solution exists, false if no solution exists.
+      *
+      * \note If there exist more than one solution, this method will arbitrarily choose one.
+      *       If you need a complete analysis of the space of solutions, take the one solution obtained
+      *       by this method and add to it elements of the kernel, as determined by kernel().
+      *
+      * \note The case where b is a matrix is not yet implemented. Also, this
+      *       code is space inefficient.
+      *
+      * Example: \include QR_solve.cpp
+      * Output: \verbinclude QR_solve.out
+      *
+      * \sa MatrixBase::solveTriangular(), kernel(), computeKernel(), inverse(), computeInverse()
+      */
+    template<typename OtherDerived, typename ResultType>
+    bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
+
     MatrixType matrixQ(void) const;
 
   private:
@@ -88,12 +166,11 @@ int QR<MatrixType>::rank() const
 {
   if (!m_rankIsUptodate)
   {
-    RealScalar maxCoeff = m_qr.diagonal().maxCoeff();
-    int n = std::min(m_qr.rows(),m_qr.cols());
-    m_rank = n;
-    for (int i=0; i<n; ++i)
-      if (ei_isMuchSmallerThan(m_qr.diagonal().coeff(i), maxCoeff))
-        --m_rank;
+    RealScalar maxCoeff = m_qr.diagonal().cwise().abs().maxCoeff();
+    int n = m_qr.cols();
+    m_rank = 0;
+    while(m_rank<n && !ei_isMuchSmallerThan(m_qr.diagonal().coeff(m_rank), maxCoeff))
+      ++m_rank;
     m_rankIsUptodate = true;
   }
   return m_rank;
@@ -132,7 +209,7 @@ void QR<MatrixType>::_compute(const MatrixType& matrix)
         m_hCoeffs.coeffRef(k) = 0;
       }
     }
-    else if ( (!ei_isMuchSmallerThan(beta=m_qr.col(k).end(remainingSize-1).squaredNorm(),static_cast<Scalar>(1))) || ei_imag(v0)==0 )
+    else if ( (!ei_isMuchSmallerThan(beta=m_qr.col(k).end(remainingSize-1).squaredNorm(),static_cast<Scalar>(1))) ) // FIXME what about ei_imag(v0) ??
     {
       // form k-th Householder vector
       beta = ei_sqrt(ei_abs2(v0)+beta);
@@ -160,9 +237,40 @@ void QR<MatrixType>::_compute(const MatrixType& matrix)
   }
 }
 
+template<typename MatrixType>
+template<typename OtherDerived, typename ResultType>
+bool QR<MatrixType>::solve(
+  const MatrixBase<OtherDerived>& b,
+  ResultType *result
+) const
+{
+  const int rows = m_qr.rows();
+  ei_assert(b.rows() == rows);
+  result->resize(rows, b.cols());
+  
+  // TODO(keir): There is almost certainly a faster way to multiply by
+  // Q^T without explicitly forming matrixQ(). Investigate.
+  *result = matrixQ().transpose()*b;
+  
+  if(!isSurjective())
+  {
+    // is result is in the image of R ?
+    RealScalar biggest_in_res = result->corner(TopLeft, m_rank, result->cols()).cwise().abs().maxCoeff();
+    for(int col = 0; col < result->cols(); ++col)
+      for(int row = m_rank; row < result->rows(); ++row)
+        if(!ei_isMuchSmallerThan(result->coeff(row,col), biggest_in_res))
+          return false;
+  }
+  m_qr.corner(TopLeft, m_rank, m_rank)
+      .template marked<UpperTriangular>()
+      .solveTriangularInPlace(result->corner(TopLeft, m_rank, result->cols()));
+
+  return true;
+}
+
 /** \returns the matrix Q */
 template<typename MatrixType>
-MatrixType QR<MatrixType>::matrixQ(void) const
+MatrixType QR<MatrixType>::matrixQ() const
 {
   // compute the product Q_0 Q_1 ... Q_n-1,
   // where Q_k is the k-th Householder transformation I - h_k v_k v_k'