port the qr module to ei_solve_xxx.

5 files changed, 181 insertions, 174 deletions

diff --git a/Eigen/src/LU/PartialPivLU.h b/Eigen/src/LU/PartialPivLU.h
index 8f3b7dfc1..eeec3533f 100644
--- a/Eigen/src/LU/PartialPivLU.h
+++ b/Eigen/src/LU/PartialPivLU.h
@@ -198,8 +198,6 @@ PartialPivLU<MatrixType>::PartialPivLU(const MatrixType& matrix)
   compute(matrix);
 }
 
-
-
 /** This is the blocked version of ei_fullpivlu_unblocked() */
 template<typename Scalar, int StorageOrder>
 struct ei_partial_lu_impl
diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h
index 05287ff3c..a774fdd73 100644
--- a/Eigen/src/QR/ColPivHouseholderQR.h
+++ b/Eigen/src/QR/ColPivHouseholderQR.h
@@ -42,17 +42,17 @@
   *
   * \sa MatrixBase::colPivHouseholderQr()
   */
-template<typename MatrixType> class ColPivHouseholderQR
+template<typename _MatrixType> class ColPivHouseholderQR
 {
   public:
 
+    typedef _MatrixType MatrixType;
     enum {
       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
       Options = MatrixType::Options,
       DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime)
     };
-
     typedef typename MatrixType::Scalar Scalar;
     typedef typename MatrixType::RealScalar RealScalar;
     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
@@ -83,22 +83,27 @@ template<typename MatrixType> class ColPivHouseholderQR
     /** 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.
       *
-      * \returns \c true if a solution exists, \c false if no solution 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 a solution.
       *
       * \note The case where b is a matrix is not yet implemented. Also, this
       *       code is space inefficient.
       *
+      * \note_about_checking_solutions
+      *
+      * \note_about_arbitrary_choice_of_solution
+      *
       * Example: \include ColPivHouseholderQR_solve.cpp
       * Output: \verbinclude ColPivHouseholderQR_solve.out
       */
-    template<typename OtherDerived, typename ResultType>
-    bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
+    template<typename Rhs>
+    inline const ei_solve_return_value<ColPivHouseholderQR, Rhs>
+    solve(const MatrixBase<Rhs>& b) const
+    {
+      ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
+      return ei_solve_return_value<ColPivHouseholderQR, Rhs>(*this, b.derived());
+    }
 
     HouseholderSequenceType matrixQ(void) const;
 
@@ -204,36 +209,24 @@ template<typename MatrixType> class ColPivHouseholderQR
       return isInjective() && isSurjective();
     }
 
-    /** Computes the inverse of the matrix of which *this is the QR decomposition.
-      *
-      * \param result a pointer to the matrix into which to store the inverse. Resized if needed.
-      *
-      * \note If this matrix is not invertible, *result is left with undefined coefficients.
-      *       Use isInvertible() to first determine whether this matrix is invertible.
-      *
-      * \sa inverse()
-      */
-    inline void computeInverse(MatrixType *result) const
-    {
-      ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
-      ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the inverse of a non-square matrix!");
-      solve(MatrixType::Identity(m_qr.rows(), m_qr.cols()), result);
-    }
-
     /** \returns the inverse of the matrix of which *this is the QR decomposition.
       *
       * \note If this matrix is not invertible, the returned matrix has undefined coefficients.
       *       Use isInvertible() to first determine whether this matrix is invertible.
-      *
-      * \sa computeInverse()
       */
-    inline MatrixType inverse() const
+    inline const
+    ei_solve_return_value<ColPivHouseholderQR, NestByValue<typename MatrixType::IdentityReturnType> >
+    inverse() const
     {
-      MatrixType result;
-      computeInverse(&result);
-      return result;
+      ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
+      return ei_solve_return_value<ColPivHouseholderQR,NestByValue<typename MatrixType::IdentityReturnType> >
+               (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()).nestByValue());
     }
 
+    inline int rows() const { return m_qr.rows(); }
+    inline int cols() const { return m_qr.cols(); }
+    const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
+
   protected:
     MatrixType m_qr;
     HCoeffsType m_hCoeffs;
@@ -331,50 +324,56 @@ ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const
   return *this;
 }
 
-template<typename MatrixType>
-template<typename OtherDerived, typename ResultType>
-bool ColPivHouseholderQR<MatrixType>::solve(
-  const MatrixBase<OtherDerived>& b,
-  ResultType *result
-) const
+template<typename MatrixType, typename Rhs, typename Dest>
+struct ei_solve_impl<ColPivHouseholderQR<MatrixType>, Rhs, Dest>
+  : ei_solve_return_value<ColPivHouseholderQR<MatrixType>, Rhs>
 {
-  ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
-  result->resize(m_qr.cols(), b.cols());
-  if(m_rank==0)
+  void evalTo(Dest& dst) const
   {
-    if(b.squaredNorm() == RealScalar(0))
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::RealScalar RealScalar;
+    const ColPivHouseholderQR<MatrixType>& dec = this->m_dec;
+    const Rhs& rhs = this->m_rhs;
+    const int rows = dec.rows(), cols = dec.cols();
+    dst.resize(cols, rhs.cols());
+    ei_assert(rhs.rows() == rows);
+
+    // FIXME introduce nonzeroPivots() and use it here. and more generally,
+    // make the same improvements in this dec as in FullPivLU.
+    if(dec.rank()==0)
     {
-      result->setZero();
-      return true;
+      dst.setZero();
+      return;
     }
-    else return false;
-  }
 
-  const int rows = m_qr.rows();
-  ei_assert(b.rows() == rows);
+    typename Rhs::PlainMatrixType c(rhs);
 
-  typename OtherDerived::PlainMatrixType c(b);
+    // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
+    c.applyOnTheLeft(makeHouseholderSequence(
+      dec.matrixQR().corner(TopLeft,rows,dec.rank()),
+      dec.hCoeffs().start(dec.rank())).transpose()
+    );
 
-  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
-  c.applyOnTheLeft(makeHouseholderSequence(m_qr.corner(TopLeft,rows,m_rank), m_hCoeffs.start(m_rank)).transpose());
-
-  if(!isSurjective())
-  {
-    // is c is in the image of R ?
-    RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff();
-    RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-m_rank, c.cols()).cwise().abs().maxCoeff();
-    if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision*4))
-      return false;
-  }
+    if(!dec.isSurjective())
+    {
+      // is c is in the image of R ?
+      RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, dec.rank(), c.cols()).cwise().abs().maxCoeff();
+      RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-dec.rank(), c.cols()).cwise().abs().maxCoeff();
+      // FIXME brain dead
+      const RealScalar m_precision = epsilon<Scalar>() * std::min(rows,cols);
+      if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision*4))
+        return;
+    }
 
-  m_qr.corner(TopLeft, m_rank, m_rank)
-      .template triangularView<UpperTriangular>()
-      .solveInPlace(c.corner(TopLeft, m_rank, c.cols()));
+    dec.matrixQR()
+       .corner(TopLeft, dec.rank(), dec.rank())
+       .template triangularView<UpperTriangular>()
+       .solveInPlace(c.corner(TopLeft, dec.rank(), c.cols()));
 
-  for(int i = 0; i < m_rank; ++i) result->row(m_cols_permutation.coeff(i)) = c.row(i);
-  for(int i = m_rank; i < m_qr.cols(); ++i) result->row(m_cols_permutation.coeff(i)).setZero();
-  return true;
-}
+    for(int i = 0; i < dec.rank(); ++i) dst.row(dec.colsPermutation().coeff(i)) = c.row(i);
+    for(int i = dec.rank(); i < cols; ++i) dst.row(dec.colsPermutation().coeff(i)).setZero();
+  }
+};
 
 /** \returns the matrix Q as a sequence of householder transformations */
 template<typename MatrixType>
diff --git a/Eigen/src/QR/FullPivHouseholderQR.h b/Eigen/src/QR/FullPivHouseholderQR.h
index 07ec343a5..36ec71b95 100644
--- a/Eigen/src/QR/FullPivHouseholderQR.h
+++ b/Eigen/src/QR/FullPivHouseholderQR.h
@@ -42,17 +42,17 @@
   *
   * \sa MatrixBase::fullPivHouseholderQr()
   */
-template<typename MatrixType> class FullPivHouseholderQR
+template<typename _MatrixType> class FullPivHouseholderQR
 {
   public:
 
+    typedef _MatrixType MatrixType;
     enum {
       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
       Options = MatrixType::Options,
       DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime)
     };
-
     typedef typename MatrixType::Scalar Scalar;
     typedef typename MatrixType::RealScalar RealScalar;
     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime> MatrixQType;
@@ -78,22 +78,27 @@ template<typename MatrixType> class FullPivHouseholderQR
     /** 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.
       *
-      * \returns \c true if a solution exists, \c false if no solution 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 a solution.
       *
       * \note The case where b is a matrix is not yet implemented. Also, this
       *       code is space inefficient.
       *
+      * \note_about_checking_solutions
+      *
+      * \note_about_arbitrary_choice_of_solution
+      *
       * Example: \include FullPivHouseholderQR_solve.cpp
       * Output: \verbinclude FullPivHouseholderQR_solve.out
       */
-    template<typename OtherDerived, typename ResultType>
-    bool solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
+    template<typename Rhs>
+    inline const ei_solve_return_value<FullPivHouseholderQR, Rhs>
+    solve(const MatrixBase<Rhs>& b) const
+    {
+      ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
+      return ei_solve_return_value<FullPivHouseholderQR, Rhs>(*this, b.derived());
+    }
 
     MatrixQType matrixQ(void) const;
 
@@ -205,36 +210,23 @@ template<typename MatrixType> class FullPivHouseholderQR
       return isInjective() && isSurjective();
     }
 
-    /** Computes the inverse of the matrix of which *this is the QR decomposition.
-      *
-      * \param result a pointer to the matrix into which to store the inverse. Resized if needed.
-      *
-      * \note If this matrix is not invertible, *result is left with undefined coefficients.
-      *       Use isInvertible() to first determine whether this matrix is invertible.
-      *
-      * \sa inverse()
-      */
-    inline void computeInverse(MatrixType *result) const
-    {
-      ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
-      ei_assert(m_qr.rows() == m_qr.cols() && "You can't take the inverse of a non-square matrix!");
-      solve(MatrixType::Identity(m_qr.rows(), m_qr.cols()), result);
-    }
-
     /** \returns the inverse of the matrix of which *this is the QR decomposition.
       *
       * \note If this matrix is not invertible, the returned matrix has undefined coefficients.
       *       Use isInvertible() to first determine whether this matrix is invertible.
-      *
-      * \sa computeInverse()
-      */
-    inline MatrixType inverse() const
+      */    inline const
+    ei_solve_return_value<FullPivHouseholderQR, NestByValue<typename MatrixType::IdentityReturnType> >
+    inverse() const
     {
-      MatrixType result;
-      computeInverse(&result);
-      return result;
+      ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
+      return ei_solve_return_value<FullPivHouseholderQR,NestByValue<typename MatrixType::IdentityReturnType> >
+               (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()).nestByValue());
     }
 
+    inline int rows() const { return m_qr.rows(); }
+    inline int cols() const { return m_qr.cols(); }
+    const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
+
   protected:
     MatrixType m_qr;
     HCoeffsType m_hCoeffs;
@@ -340,56 +332,59 @@ FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(cons
   return *this;
 }
 
-template<typename MatrixType>
-template<typename OtherDerived, typename ResultType>
-bool FullPivHouseholderQR<MatrixType>::solve(
-  const MatrixBase<OtherDerived>& b,
-  ResultType *result
-) const
+template<typename MatrixType, typename Rhs, typename Dest>
+struct ei_solve_impl<FullPivHouseholderQR<MatrixType>, Rhs, Dest>
+  : ei_solve_return_value<FullPivHouseholderQR<MatrixType>, Rhs>
 {
-  ei_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
-  result->resize(m_qr.cols(), b.cols());
-  if(m_rank==0)
+  void evalTo(Dest& dst) const
   {
-    if(b.squaredNorm() == RealScalar(0))
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::RealScalar RealScalar;
+    const FullPivHouseholderQR<MatrixType>& dec = this->m_dec;
+    const Rhs& rhs = this->m_rhs;
+    const int rows = dec.rows(), cols = dec.cols();
+    dst.resize(cols, rhs.cols());
+    ei_assert(rhs.rows() == rows);
+
+    // FIXME introduce nonzeroPivots() and use it here. and more generally,
+    // make the same improvements in this dec as in FullPivLU.
+    if(dec.rank()==0)
     {
-      result->setZero();
-      return true;
+      dst.setZero();
+      return;
     }
-    else return false;
-  }
 
-  const int rows = m_qr.rows();
-  const int cols = b.cols();
-  ei_assert(b.rows() == rows);
+    typename Rhs::PlainMatrixType c(rhs);
 
-  typename OtherDerived::PlainMatrixType c(b);
+    Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs.cols());
+    for (int k = 0; k < dec.rank(); ++k)
+    {
+      int remainingSize = rows-k;
+      c.row(k).swap(c.row(dec.rowsTranspositions().coeff(k)));
+      c.corner(BottomRight, remainingSize, rhs.cols())
+       .applyHouseholderOnTheLeft(dec.matrixQR().col(k).end(remainingSize-1),
+                                  dec.hCoeffs().coeff(k), &temp.coeffRef(0));
+    }
 
-  Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols);
-  for (int k = 0; k < m_rank; ++k)
-  {
-    int remainingSize = rows-k;
-    c.row(k).swap(c.row(m_rows_transpositions.coeff(k)));
-    c.corner(BottomRight, remainingSize, cols)
-     .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0));
-  }
+    if(!dec.isSurjective())
+    {
+      // is c is in the image of R ?
+      RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, dec.rank(), c.cols()).cwise().abs().maxCoeff();
+      RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-dec.rank(), c.cols()).cwise().abs().maxCoeff();
+      // FIXME brain dead
+      const RealScalar m_precision = epsilon<Scalar>() * std::min(rows,cols);
+      if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
+        return;
+    }
+    dec.matrixQR()
+       .corner(TopLeft, dec.rank(), dec.rank())
+       .template triangularView<UpperTriangular>()
+       .solveInPlace(c.corner(TopLeft, dec.rank(), c.cols()));
 
-  if(!isSurjective())
-  {
-    // is c is in the image of R ?
-    RealScalar biggest_in_upper_part_of_c = c.corner(TopLeft, m_rank, c.cols()).cwise().abs().maxCoeff();
-    RealScalar biggest_in_lower_part_of_c = c.corner(BottomLeft, rows-m_rank, c.cols()).cwise().abs().maxCoeff();
-    if(!ei_isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
-      return false;
+    for(int i = 0; i < dec.rank(); ++i) dst.row(dec.colsPermutation().coeff(i)) = c.row(i);
+    for(int i = dec.rank(); i < cols; ++i) dst.row(dec.colsPermutation().coeff(i)).setZero();
   }
-  m_qr.corner(TopLeft, m_rank, m_rank)
-      .template triangularView<UpperTriangular>()
-      .solveInPlace(c.corner(TopLeft, m_rank, c.cols()));
-
-  for(int i = 0; i < m_rank; ++i) result->row(m_cols_permutation.coeff(i)) = c.row(i);
-  for(int i = m_rank; i < m_qr.cols(); ++i) result->row(m_cols_permutation.coeff(i)).setZero();
-  return true;
-}
+};
 
 /** \returns the matrix Q */
 template<typename MatrixType>
diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h
index a32aa4eaf..6db0411d9 100644
--- a/Eigen/src/QR/HouseholderQR.h
+++ b/Eigen/src/QR/HouseholderQR.h
@@ -46,17 +46,17 @@
   *
   * \sa MatrixBase::householderQr()
   */
-template<typename MatrixType> class HouseholderQR
+template<typename _MatrixType> class HouseholderQR
 {
   public:
 
+    typedef _MatrixType MatrixType;
     enum {
       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
       Options = MatrixType::Options,
       DiagSizeAtCompileTime = EIGEN_ENUM_MIN(ColsAtCompileTime,RowsAtCompileTime)
     };
-
     typedef typename MatrixType::Scalar Scalar;
     typedef typename MatrixType::RealScalar RealScalar;
     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, AutoAlign | (ei_traits<MatrixType>::Flags&RowMajorBit ? RowMajor : ColMajor)> MatrixQType;
@@ -85,19 +85,26 @@ template<typename MatrixType> class HouseholderQR
       *
       * \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 a solution.
       *
       * \note The case where b is a matrix is not yet implemented. Also, this
       *       code is space inefficient.
       *
+      * \note_about_checking_solutions
+      *
+      * \note_about_arbitrary_choice_of_solution
+      *
       * Example: \include HouseholderQR_solve.cpp
       * Output: \verbinclude HouseholderQR_solve.out
       */
-    template<typename OtherDerived, typename ResultType>
-    void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
-
+    template<typename Rhs>
+    inline const ei_solve_return_value<HouseholderQR, Rhs>
+    solve(const MatrixBase<Rhs>& b) const
+    {
+      ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
+      return ei_solve_return_value<HouseholderQR, Rhs>(*this, b.derived());
+    }
+    
     MatrixQType matrixQ() const;
 
     HouseholderSequenceType matrixQAsHouseholderSequence() const
@@ -145,6 +152,10 @@ template<typename MatrixType> class HouseholderQR
       */
     typename MatrixType::RealScalar logAbsDeterminant() const;
 
+    inline int rows() const { return m_qr.rows(); }
+    inline int cols() const { return m_qr.cols(); }
+    const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
+
   protected:
     MatrixType m_qr;
     HCoeffsType m_hCoeffs;
@@ -198,31 +209,36 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
   return *this;
 }
 
-template<typename MatrixType>
-template<typename OtherDerived, typename ResultType>
-void HouseholderQR<MatrixType>::solve(
-  const MatrixBase<OtherDerived>& b,
-  ResultType *result
-) const
+template<typename MatrixType, typename Rhs, typename Dest>
+struct ei_solve_impl<HouseholderQR<MatrixType>, Rhs, Dest>
+  : ei_solve_return_value<HouseholderQR<MatrixType>, Rhs>
 {
-  ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
-  result->derived().resize(m_qr.cols(), b.cols());
-  const int rows = m_qr.rows();
-  const int rank = std::min(m_qr.rows(), m_qr.cols());
-  ei_assert(b.rows() == rows);
-
-  typename OtherDerived::PlainMatrixType c(b);
-
-  // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
-  c.applyOnTheLeft(makeHouseholderSequence(m_qr.corner(TopLeft,rows,rank), m_hCoeffs.start(rank)).transpose());
-
-  m_qr.corner(TopLeft, rank, rank)
-      .template triangularView<UpperTriangular>()
-      .solveInPlace(c.corner(TopLeft, rank, c.cols()));
-
-  result->corner(TopLeft, rank, c.cols()) = c.corner(TopLeft,rank, c.cols());
-  result->corner(BottomLeft, result->rows()-rank, c.cols()).setZero();
-}
+  void evalTo(Dest& dst) const
+  {
+    const HouseholderQR<MatrixType>& dec = this->m_dec;
+    const Rhs& rhs = this->m_rhs;
+    const int rows = dec.rows(), cols = dec.cols();
+    dst.resize(cols, rhs.cols());
+    const int rank = std::min(rows, cols);
+    ei_assert(rhs.rows() == rows);
+
+    typename Rhs::PlainMatrixType c(rhs);
+
+    // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
+    c.applyOnTheLeft(makeHouseholderSequence(
+      dec.matrixQR().corner(TopLeft,rows,rank),
+      dec.hCoeffs().start(rank)).transpose()
+    );
+
+    dec.matrixQR()
+       .corner(TopLeft, rank, rank)
+       .template triangularView<UpperTriangular>()
+       .solveInPlace(c.corner(TopLeft, rank, c.cols()));
+
+    dst.corner(TopLeft, rank, c.cols()) = c.corner(TopLeft, rank, c.cols());
+    dst.corner(BottomLeft, cols-rank, c.cols()).setZero();
+  }
+};
 
 /** \returns the matrix Q */
 template<typename MatrixType>
diff --git a/Eigen/src/SVD/SVD.h b/Eigen/src/SVD/SVD.h
index b43123384..8ca425525 100644
--- a/Eigen/src/SVD/SVD.h
+++ b/Eigen/src/SVD/SVD.h
@@ -86,7 +86,6 @@ template<typename _MatrixType> class SVD
       * \note_about_checking_solutions
       *
       * \note_about_arbitrary_choice_of_solution
-      * \note_about_using_kernel_to_study_multiple_solutions
       *
       * \sa MatrixBase::svd(),
       */