* rewrite of the QR decomposition:

  - works for complex
  - allows direct access to the matrix R
* removed the scale by the matrix dimensions in MatrixBase::isMuchSmallerThan(scalar)

3 files changed, 66 insertions, 45 deletions

diff --git a/Eigen/src/Core/Fuzzy.h b/Eigen/src/Core/Fuzzy.h
index e689fc913..0b2305c57 100644
--- a/Eigen/src/Core/Fuzzy.h
+++ b/Eigen/src/Core/Fuzzy.h
@@ -63,7 +63,10 @@ bool MatrixBase<Derived>::isApprox(
   * \note The fuzzy compares are done multiplicatively. A vector \f$ v \f$ is
   * considered to be much smaller than \f$ x \f$ within precision \f$ p \f$ if
   * \f[ \Vert v \Vert \leqslant p\,\vert x\vert. \f]
-  * For matrices, the comparison is done using the Hilbert-Schmidt norm.
+  *
+  * For matrices, the comparison is done using the Hilbert-Schmidt norm. For this reason,
+  * the value of the reference scalar \a other should come from the Hilbert-Schmidt norm
+  * of a reference matrix of same dimensions.
   *
   * \sa isApprox(), isMuchSmallerThan(const MatrixBase<OtherDerived>&, RealScalar) const
   */
@@ -73,7 +76,7 @@ bool MatrixBase<Derived>::isMuchSmallerThan(
   typename NumTraits<Scalar>::Real prec
 ) const
 {
-  return cwiseAbs2().sum() <= prec * prec * other * other * cols() * rows();
+  return cwiseAbs2().sum() <= prec * prec * other * other;
 }
 
 /** \returns \c true if the norm of \c *this is much smaller than the norm of \a other,
diff --git a/Eigen/src/QR/QR.h b/Eigen/src/QR/QR.h
index bd01cf71e..2bf6c72b2 100644
--- a/Eigen/src/QR/QR.h
+++ b/Eigen/src/QR/QR.h
@@ -32,12 +32,7 @@
   * \param MatrixType the type of the matrix of which we are computing the QR decomposition
   *
   * This class performs a QR decomposition using Householder transformations. The result is
-  * stored in a compact way.
-  *
-  * \todo add convenient method to direclty use the result in a compact way. First need to determine
-  * typical use cases though.
-  *
-  * \todo what about complex matrices ?
+  * stored in a compact way compatible with LAPACK.
   *
   * \sa MatrixBase::qr()
   */
@@ -46,20 +41,28 @@ template<typename MatrixType> class QR
   public:
 
     typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::RealScalar RealScalar;
+    typedef Block<MatrixType, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixRBlockType;
     typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, MatrixType::ColsAtCompileTime> MatrixTypeR;
     typedef Matrix<Scalar, MatrixType::ColsAtCompileTime, 1> VectorType;
 
     QR(const MatrixType& matrix)
       : m_qr(matrix.rows(), matrix.cols()),
-        m_norms(matrix.cols())
+        m_hCoeffs(matrix.cols())
     {
       _compute(matrix);
     }
 
     /** \returns whether or not the matrix is of full rank */
-    bool isFullRank() const { return ei_isMuchSmallerThan(m_norms.cwiseAbs().minCoeff(), Scalar(1)); }
+    bool isFullRank() const { return ei_isMuchSmallerThan(m_hCoeffs.cwiseAbs().minCoeff(), Scalar(1)); }
 
-    MatrixTypeR matrixR(void) const;
+    /** \returns a read-only expression of the matrix R of the actual the QR decomposition */
+    const Extract<NestByValue<MatrixRBlockType>, Upper>
+    matrixR(void) const
+    {
+      int cols = m_qr.cols();
+      return MatrixRBlockType(m_qr, 0, 0, cols, cols).nestByValue().template extract<Upper>();
+    }
 
     MatrixType matrixQ(void) const;
 
@@ -69,7 +72,7 @@ template<typename MatrixType> class QR
 
   protected:
     MatrixType m_qr;
-    VectorType m_norms;
+    VectorType m_hCoeffs;
 };
 
 template<typename MatrixType>
@@ -83,58 +86,73 @@ void QR<MatrixType>::_compute(const MatrixType& matrix)
   {
     int remainingSize = rows-k;
 
-    Scalar nrm = m_qr.col(k).end(remainingSize).norm();
+    RealScalar beta;
+    Scalar v0 = m_qr.col(k).coeff(k);
 
-    if (nrm != Scalar(0))
+    if (remainingSize==1)
+    {
+      if (NumTraits<Scalar>::IsComplex)
+      {
+        // Householder transformation on the remaining single scalar
+        beta = ei_abs(v0);
+        if (ei_real(v0)>0)
+          beta = -beta;
+        m_qr.coeffRef(k,k) = beta;
+        m_hCoeffs.coeffRef(k) = (beta - v0) / beta;
+      }
+      else
+      {
+        m_hCoeffs.coeffRef(k) = 0;
+      }
+    }
+    else if ( (!ei_isMuchSmallerThan(beta=m_qr.col(k).end(remainingSize-1).norm2(),static_cast<Scalar>(1))) || ei_imag(v0)==0 )
     {
       // form k-th Householder vector
-      if (m_qr(k,k) < 0)
-        nrm = -nrm;
-
-      m_qr.col(k).end(rows-k) /= nrm;
-      m_qr(k,k) += 1.0;
-
-      // apply transformation to remaining columns
+      beta = ei_sqrt(ei_abs2(v0)+beta);
+      if (ei_real(v0)>=0.)
+        beta = -beta;
+      m_qr.col(k).end(remainingSize-1) /= v0-beta;
+      m_qr.coeffRef(k,k) = beta;
+      Scalar h = m_hCoeffs.coeffRef(k) = (beta - v0) / beta;
+
+      // apply the Householder transformation (I - h v v') to remaining columns, i.e.,
+      // R <- (I - h v v') * R   where v = [1,m_qr(k+1,k), m_qr(k+2,k), ...]
       int remainingCols = cols - k -1;
       if (remainingCols>0)
       {
-        m_qr.corner(BottomRight, remainingSize, remainingCols) -= (1./m_qr(k,k)) * m_qr.col(k).end(remainingSize)
-          * (m_qr.col(k).end(remainingSize).transpose() * m_qr.corner(BottomRight, remainingSize, remainingCols));
+        m_qr.coeffRef(k,k) = Scalar(1);
+        m_qr.corner(BottomRight, remainingSize, remainingCols) -= ei_conj(h) * m_qr.col(k).end(remainingSize)
+            * (m_qr.col(k).end(remainingSize).adjoint() * m_qr.corner(BottomRight, remainingSize, remainingCols));
+        m_qr.coeffRef(k,k) = beta;
       }
     }
-    m_norms[k] = -nrm;
+    else
+    {
+      m_hCoeffs.coeffRef(k) = 0;
+    }
   }
 }
 
-/** \returns the matrix R */
-template<typename MatrixType>
-typename QR<MatrixType>::MatrixTypeR QR<MatrixType>::matrixR(void) const
-{
-  int cols = m_qr.cols();
-  MatrixTypeR res = m_qr.block(0,0,cols,cols).template extract<StrictlyUpper>();
-  res.diagonal() = m_norms;
-  return res;
-}
-
 /** \returns the matrix Q */
 template<typename MatrixType>
 MatrixType QR<MatrixType>::matrixQ(void) 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'
+  // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
   int rows = m_qr.rows();
   int cols = m_qr.cols();
   MatrixType res = MatrixType::identity(rows, cols);
   for (int k = cols-1; k >= 0; k--)
   {
-    for (int j = k; j < cols; j++)
-    {
-      if (res(k,k) != Scalar(0))
-      {
-        int endLength = rows-k;
-        Scalar s = -(m_qr.col(k).end(endLength).transpose() * res.col(j).end(endLength))(0,0) / m_qr(k,k);
-
-        res.col(j).end(endLength) += s * m_qr.col(k).end(endLength);
-      }
-    }
+    // to make easier the computation of the transformation, let's temporarily
+    // overwrite m_qr(k,k) such that the end of m_qr.col(k) is exactly our Householder vector.
+    Scalar beta = m_qr.coeff(k,k);
+    m_qr.const_cast_derived().coeffRef(k,k) = 1;
+    int endLength = rows-k;
+    res.corner(BottomRight,endLength, cols-k) -= ((m_hCoeffs.coeff(k) * m_qr.col(k).end(endLength))
+      * (m_qr.col(k).end(endLength).adjoint() * res.corner(BottomRight,endLength, cols-k)).lazy()).lazy();
+    m_qr.const_cast_derived().coeffRef(k,k) = beta;
   }
   return res;
 }
diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/QR/Tridiagonalization.h
index b347b19a6..e5b412c46 100755
--- a/Eigen/src/QR/Tridiagonalization.h
+++ b/Eigen/src/QR/Tridiagonalization.h
@@ -226,7 +226,7 @@ Tridiagonalization<MatrixType>::matrixQ(void) const
     m_matrix.const_cast_derived().coeffRef(i+1,i) = 1;
 
     matQ.corner(BottomRight,n-i-1,n-i-1) -=
-      ((m_hCoeffs[i] * m_matrix.col(i).end(n-i-1)) *
+      ((m_hCoeffs.coeff(i) * m_matrix.col(i).end(n-i-1)) *
       (m_matrix.col(i).end(n-i-1).adjoint() * matQ.corner(BottomRight,n-i-1,n-i-1)).lazy()).lazy();
 
     m_matrix.const_cast_derived().coeffRef(i+1,i) = tmp;