* add a HouseholderSequence class (not good enough yet for Triadiagonalization and HessenbergDecomposition)

* rework a bit AnyMatrixBase, and mobe it to a separate file

9 files changed, 328 insertions, 95 deletions

diff --git a/Eigen/Core b/Eigen/Core
index 854f737d6..42eb363a9 100644
--- a/Eigen/Core
+++ b/Eigen/Core
@@ -143,6 +143,7 @@ namespace Eigen {
 
 #include "src/Core/Functors.h"
 #include "src/Core/MatrixBase.h"
+#include "src/Core/AnyMatrixBase.h"
 #include "src/Core/Coeffs.h"
 
 #ifndef EIGEN_PARSED_BY_DOXYGEN // work around Doxygen bug triggered by Assign.h r814874
diff --git a/Eigen/Householder b/Eigen/Householder
index ba06bd8fb..ef3e61373 100644
--- a/Eigen/Householder
+++ b/Eigen/Householder
@@ -16,6 +16,7 @@ namespace Eigen {
   */
 
 #include "src/Householder/Householder.h"
+#include "src/Householder/HouseholderSequence.h"
 
 } // namespace Eigen
 
diff --git a/Eigen/src/Core/AnyMatrixBase.h b/Eigen/src/Core/AnyMatrixBase.h
new file mode 100644
index 000000000..cd354d7b1
--- /dev/null
+++ b/Eigen/src/Core/AnyMatrixBase.h
@@ -0,0 +1,153 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_ANYMATRIXBASE_H
+#define EIGEN_ANYMATRIXBASE_H
+
+
+/** Common base class for all classes T such that MatrixBase has an operator=(T) and a constructor MatrixBase(T).
+  *
+  * In other words, an AnyMatrixBase object is an object that can be copied into a MatrixBase.
+  *
+  * Besides MatrixBase-derived classes, this also includes special matrix classes such as diagonal matrices, etc.
+  *
+  * Notice that this class is trivial, it is only used to disambiguate overloaded functions.
+  */
+template<typename Derived> struct AnyMatrixBase
+{
+  typedef typename ei_plain_matrix_type<Derived>::type PlainMatrixType;
+
+  Derived& derived() { return *static_cast<Derived*>(this); }
+  const Derived& derived() const { return *static_cast<const Derived*>(this); }
+
+  /** \returns the number of rows. \sa cols(), RowsAtCompileTime */
+  inline int rows() const { return derived().rows(); }
+  /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/
+  inline int cols() const { return derived().cols(); }
+
+  /** \internal Don't use it, but do the equivalent: \code dst = *this; \endcode */
+  template<typename Dest> inline void evalTo(Dest& dst) const
+  { derived().evalTo(dst); }
+
+  /** \internal Don't use it, but do the equivalent: \code dst += *this; \endcode */
+  template<typename Dest> inline void addToDense(Dest& dst) const
+  {
+    // This is the default implementation,
+    // derived class can reimplement it in a more optimized way.
+    typename Dest::PlainMatrixType res(rows(),cols());
+    evalTo(res);
+    dst += res;
+  }
+
+  /** \internal Don't use it, but do the equivalent: \code dst -= *this; \endcode */
+  template<typename Dest> inline void subToDense(Dest& dst) const
+  {
+    // This is the default implementation,
+    // derived class can reimplement it in a more optimized way.
+    typename Dest::PlainMatrixType res(rows(),cols());
+    evalTo(res);
+    dst -= res;
+  }
+
+  /** \internal Don't use it, but do the equivalent: \code dst.applyOnTheRight(*this); \endcode */
+  template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
+  {
+    // This is the default implementation,
+    // derived class can reimplement it in a more optimized way.
+    dst = dst * this->derived();
+  }
+
+  /** \internal Don't use it, but do the equivalent: \code dst.applyOnTheLeft(*this); \endcode */
+  template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
+  {
+    // This is the default implementation,
+    // derived class can reimplement it in a more optimized way.
+    dst = this->derived() * dst;
+  }
+
+};
+
+/***************************************************************************
+* Implementation of matrix base methods
+***************************************************************************/
+
+/** Copies the generic expression \a other into *this. \returns a reference to *this.
+  * The expression must provide a (templated) evalToDense(Derived& dst) const function
+  * which does the actual job. In practice, this allows any user to write its own
+  * special matrix without having to modify MatrixBase */
+template<typename Derived>
+template<typename OtherDerived>
+Derived& MatrixBase<Derived>::operator=(const AnyMatrixBase<OtherDerived> &other)
+{
+  other.derived().evalTo(derived());
+  return derived();
+}
+
+template<typename Derived>
+template<typename OtherDerived>
+Derived& MatrixBase<Derived>::operator+=(const AnyMatrixBase<OtherDerived> &other)
+{
+  other.derived().addToDense(derived());
+  return derived();
+}
+
+template<typename Derived>
+template<typename OtherDerived>
+Derived& MatrixBase<Derived>::operator-=(const AnyMatrixBase<OtherDerived> &other)
+{
+  other.derived().subToDense(derived());
+  return derived();
+}
+
+/** replaces \c *this by \c *this * \a other.
+  *
+  * \returns a reference to \c *this
+  */
+template<typename Derived>
+template<typename OtherDerived>
+inline Derived&
+MatrixBase<Derived>::operator*=(const AnyMatrixBase<OtherDerived> &other)
+{
+  other.derived().applyThisOnTheRight(derived());
+  return derived();
+}
+
+/** replaces \c *this by \c *this * \a other. It is equivalent to MatrixBase::operator*=() */
+template<typename Derived>
+template<typename OtherDerived>
+inline void MatrixBase<Derived>::applyOnTheRight(const AnyMatrixBase<OtherDerived> &other)
+{
+  other.derived().applyThisOnTheRight(derived());
+}
+
+/** replaces \c *this by \c *this * \a other. */
+template<typename Derived>
+template<typename OtherDerived>
+inline void MatrixBase<Derived>::applyOnTheLeft(const AnyMatrixBase<OtherDerived> &other)
+{
+  other.derived().applyThisOnTheLeft(derived());
+}
+
+#endif // EIGEN_ANYMATRIXBASE_H
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index f267aa34d..4835f167c 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -26,44 +26,6 @@
 #ifndef EIGEN_MATRIXBASE_H
 #define EIGEN_MATRIXBASE_H
 
-
-/** Common base class for all classes T such that MatrixBase has an operator=(T) and a constructor MatrixBase(T).
-  *
-  * In other words, an AnyMatrixBase object is an object that can be copied into a MatrixBase.
-  *
-  * Besides MatrixBase-derived classes, this also includes special matrix classes such as diagonal matrices, etc.
-  *
-  * Notice that this class is trivial, it is only used to disambiguate overloaded functions.
-  */
-template<typename Derived> struct AnyMatrixBase
-{
-  typedef typename ei_plain_matrix_type<Derived>::type PlainMatrixType;
-
-  Derived& derived() { return *static_cast<Derived*>(this); }
-  const Derived& derived() const { return *static_cast<const Derived*>(this); }
-  /** \returns the number of rows. \sa cols(), RowsAtCompileTime */
-  inline int rows() const { return derived().rows(); }
-  /** \returns the number of columns. \sa rows(), ColsAtCompileTime*/
-  inline int cols() const { return derived().cols(); }
-
-  template<typename Dest> inline void evalTo(Dest& dst) const
-  { derived().evalTo(dst); }
-
-  template<typename Dest> inline void addToDense(Dest& dst) const
-  {
-    typename Dest::PlainMatrixType res(rows(),cols());
-    evalToDense(res);
-    dst += res;
-  }
-
-  template<typename Dest> inline void subToDense(Dest& dst) const
-  {
-    typename Dest::PlainMatrixType res(rows(),cols());
-    evalToDense(res);
-    dst -= res;
-  }
-};
-
 /** \class MatrixBase
   *
   * \brief Base class for all matrices, vectors, and expressions
@@ -96,7 +58,6 @@ template<typename Derived> class MatrixBase
 #endif // not EIGEN_PARSED_BY_DOXYGEN
 {
   public:
-
 #ifndef EIGEN_PARSED_BY_DOXYGEN
     using ei_special_scalar_op_base<Derived,typename ei_traits<Derived>::Scalar,
                 typename NumTraits<typename ei_traits<Derived>::Scalar>::Real>::operator*;
@@ -301,21 +262,14 @@ template<typename Derived> class MatrixBase
       */
     Derived& operator=(const MatrixBase& other);
 
-    /** Copies the generic expression \a other into *this. \returns a reference to *this.
-      * The expression must provide a (templated) evalToDense(Derived& dst) const function
-      * which does the actual job. In practice, this allows any user to write its own
-      * special matrix without having to modify MatrixBase */
     template<typename OtherDerived>
-    Derived& operator=(const AnyMatrixBase<OtherDerived> &other)
-    { other.derived().evalToDense(derived()); return derived(); }
+    Derived& operator=(const AnyMatrixBase<OtherDerived> &other);
 
     template<typename OtherDerived>
-    Derived& operator+=(const AnyMatrixBase<OtherDerived> &other)
-    { other.derived().addToDense(derived()); return derived(); }
+    Derived& operator+=(const AnyMatrixBase<OtherDerived> &other);
 
     template<typename OtherDerived>
-    Derived& operator-=(const AnyMatrixBase<OtherDerived> &other)
-    { other.derived().subToDense(derived()); return derived(); }
+    Derived& operator-=(const AnyMatrixBase<OtherDerived> &other);
 
     template<typename OtherDerived,typename OtherEvalType>
     Derived& operator=(const ReturnByValue<OtherDerived,OtherEvalType>& func);
@@ -436,6 +390,12 @@ template<typename Derived> class MatrixBase
     template<typename OtherDerived>
     Derived& operator*=(const AnyMatrixBase<OtherDerived>& other);
 
+    template<typename OtherDerived>
+    void applyOnTheLeft(const AnyMatrixBase<OtherDerived>& other);
+
+    template<typename OtherDerived>
+    void applyOnTheRight(const AnyMatrixBase<OtherDerived>& other);
+
     template<typename DiagonalDerived>
     const DiagonalProduct<Derived, DiagonalDerived, DiagonalOnTheRight>
     operator*(const DiagonalBase<DiagonalDerived> &diagonal) const;
diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h
index e7227d4f6..7f0c2df6e 100644
--- a/Eigen/src/Core/Product.h
+++ b/Eigen/src/Core/Product.h
@@ -434,18 +434,4 @@ MatrixBase<Derived>::operator*(const MatrixBase<OtherDerived> &other) const
   return typename ProductReturnType<Derived,OtherDerived>::Type(derived(), other.derived());
 }
 
-
-
-/** replaces \c *this by \c *this * \a other.
-  *
-  * \returns a reference to \c *this
-  */
-template<typename Derived>
-template<typename OtherDerived>
-inline Derived &
-MatrixBase<Derived>::operator*=(const AnyMatrixBase<OtherDerived> &other)
-{
-  return derived() = derived() * other.derived();
-}
-
 #endif // EIGEN_PRODUCT_H
diff --git a/Eigen/src/Core/TriangularMatrix.h b/Eigen/src/Core/TriangularMatrix.h
index b0362f20c..17726bca3 100644
--- a/Eigen/src/Core/TriangularMatrix.h
+++ b/Eigen/src/Core/TriangularMatrix.h
@@ -91,9 +91,9 @@ template<typename Derived> class TriangularBase : public AnyMatrixBase<Derived>
     #endif // not EIGEN_PARSED_BY_DOXYGEN
 
     template<typename DenseDerived>
-    void evalToDense(MatrixBase<DenseDerived> &other) const;
+    void evalTo(MatrixBase<DenseDerived> &other) const;
     template<typename DenseDerived>
-    void evalToDenseLazy(MatrixBase<DenseDerived> &other) const;
+    void evalToLazy(MatrixBase<DenseDerived> &other) const;
 
   protected:
 
@@ -546,23 +546,23 @@ void TriangularView<MatrixType, Mode>::lazyAssign(const TriangularBase<OtherDeri
   * If the matrix is triangular, the opposite part is set to zero. */
 template<typename Derived>
 template<typename DenseDerived>
-void TriangularBase<Derived>::evalToDense(MatrixBase<DenseDerived> &other) const
+void TriangularBase<Derived>::evalTo(MatrixBase<DenseDerived> &other) const
 {
   if(ei_traits<Derived>::Flags & EvalBeforeAssigningBit)
   {
     typename Derived::PlainMatrixType other_evaluated(rows(), cols());
-    evalToDenseLazy(other_evaluated);
+    evalToLazy(other_evaluated);
     other.derived().swap(other_evaluated);
   }
   else
-    evalToDenseLazy(other.derived());
+    evalToLazy(other.derived());
 }
 
 /** Assigns a triangular or selfadjoint matrix to a dense matrix.
   * If the matrix is triangular, the opposite part is set to zero. */
 template<typename Derived>
 template<typename DenseDerived>
-void TriangularBase<Derived>::evalToDenseLazy(MatrixBase<DenseDerived> &other) const
+void TriangularBase<Derived>::evalToLazy(MatrixBase<DenseDerived> &other) const
 {
   const bool unroll =   DenseDerived::SizeAtCompileTime * Derived::CoeffReadCost / 2
                      <= EIGEN_UNROLLING_LIMIT;
diff --git a/Eigen/src/Core/util/ForwardDeclarations.h b/Eigen/src/Core/util/ForwardDeclarations.h
index c5f27d80b..3f66738f0 100644
--- a/Eigen/src/Core/util/ForwardDeclarations.h
+++ b/Eigen/src/Core/util/ForwardDeclarations.h
@@ -123,6 +123,7 @@ template<typename MatrixType> class SVD;
 template<typename MatrixType, unsigned int Options = 0> class JacobiSVD;
 template<typename MatrixType, int UpLo = LowerTriangular> class LLT;
 template<typename MatrixType> class LDLT;
+template<typename VectorsType, typename CoeffsType> class HouseholderSequence;
 template<typename Scalar>     class PlanarRotation;
 
 // Geometry module:
diff --git a/Eigen/src/Householder/HouseholderSequence.h b/Eigen/src/Householder/HouseholderSequence.h
new file mode 100644
index 000000000..86395213b
--- /dev/null
+++ b/Eigen/src/Householder/HouseholderSequence.h
@@ -0,0 +1,146 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H
+#define EIGEN_HOUSEHOLDER_SEQUENCE_H
+
+/** \ingroup Householder_Module
+  * \householder_module
+  * \class HouseholderSequence
+  * \brief Represents a sequence of householder reflections with decreasing size
+  *
+  * This class represents a product sequence of householder reflections \f$ H = \Pi_0^{n-1} H_i \f$
+  * where \f$ H_i \f$ is the i-th householder transformation \f$ I - h_i v_i v_i^* \f$,
+  * \f$ v_i \f$ is the i-th householder vector \f$ [ 1, m_vectors(i+1,i), m_vectors(i+2,i), ...] \f$
+  * and \f$ h_i \f$ is the i-th householder coefficient \c m_coeffs[i].
+  *
+  * Typical usages are listed below, where H is a HouseholderSequence:
+  * \code
+  * A.applyOnTheRight(H);             // A = A * H
+  * A.applyOnTheLeft(H);              // A = H * A
+  * A.applyOnTheRight(H.adjoint());   // A = A * H^*
+  * A.applyOnTheLeft(H.adjoint());    // A = H^* * A
+  * MatrixXd Q = H;                   // conversion to a dense matrix
+  * \endcode
+  * In addition to the adjoint, you can also apply the inverse (=adjoint), the transpose, and the conjugate.
+  *
+  * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
+  */
+
+template<typename VectorsType, typename CoeffsType>
+struct ei_traits<HouseholderSequence<VectorsType,CoeffsType> >
+{
+  typedef typename VectorsType::Scalar Scalar;
+  enum {
+    RowsAtCompileTime = ei_traits<VectorsType>::RowsAtCompileTime,
+    ColsAtCompileTime = ei_traits<VectorsType>::RowsAtCompileTime,
+    MaxRowsAtCompileTime = ei_traits<VectorsType>::MaxRowsAtCompileTime,
+    MaxColsAtCompileTime = ei_traits<VectorsType>::MaxRowsAtCompileTime,
+    Flags = 0
+  };
+};
+
+template<typename VectorsType, typename CoeffsType> class HouseholderSequence
+  : public AnyMatrixBase<HouseholderSequence<VectorsType,CoeffsType> >
+{
+    typedef typename VectorsType::Scalar Scalar;
+  public:
+
+    typedef HouseholderSequence<VectorsType,
+      typename ei_meta_if<NumTraits<Scalar>::IsComplex,
+        NestByValue<typename ei_cleantype<typename CoeffsType::ConjugateReturnType>::type >,
+        CoeffsType>::ret> ConjugateReturnType;
+
+    HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans = false)
+      : m_vectors(v), m_coeffs(h), m_trans(trans)
+    {}
+
+    int rows() const { return m_vectors.rows(); }
+    int cols() const { return m_vectors.rows(); }
+
+    HouseholderSequence transpose() const
+    { return  HouseholderSequence(m_vectors, m_coeffs, !m_trans); }
+
+    ConjugateReturnType conjugate() const
+    { return  ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans); }
+
+    ConjugateReturnType adjoint() const
+    { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans); }
+
+    ConjugateReturnType inverse() const { return adjoint(); }
+
+    /** \internal */
+    template<typename DestType> void evalTo(DestType& dst) const
+    {
+      int vecs = std::min(m_vectors.cols(),m_vectors.rows());
+      int length = m_vectors.rows();
+      dst.setIdentity();
+      Matrix<Scalar,1,DestType::RowsAtCompileTime> temp(dst.rows());
+      for(int k = vecs-1; k >= 0; --k)
+      {
+        if(m_trans)
+          dst.corner(BottomRight, length-k, length-k)
+          .applyHouseholderOnTheRight(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(0));
+        else
+          dst.corner(BottomRight, length-k, length-k)
+            .applyHouseholderOnTheLeft(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(k));
+      }
+    }
+
+    /** \internal */
+    template<typename Dest> inline void applyThisOnTheRight(Dest& dst) const
+    {
+      int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors
+      int length = m_vectors.rows();                          // size of the largest householder vector
+      Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.rows());
+      for(int k = 0; k < vecs; ++k)
+      {
+        int actual_k = m_trans ? vecs-k-1 : k;
+        dst.corner(BottomRight, dst.rows(), length-k)
+           .applyHouseholderOnTheRight(m_vectors.col(k).end(length-k-1), m_coeffs.coeff(k), &temp.coeffRef(0));
+      }
+    }
+
+    /** \internal */
+    template<typename Dest> inline void applyThisOnTheLeft(Dest& dst) const
+    {
+      int vecs = std::min(m_vectors.cols(),m_vectors.rows()); // number of householder vectors
+      int length = m_vectors.rows();                          // size of the largest householder vector
+      Matrix<Scalar,1,Dest::ColsAtCompileTime> temp(dst.cols());
+      for(int k = 0; k < vecs; ++k)
+      {
+        int actual_k = m_trans ? k : vecs-k-1;
+        dst.corner(BottomRight, length-actual_k, dst.cols())
+           .applyHouseholderOnTheLeft(m_vectors.col(actual_k).end(length-actual_k-1), m_coeffs.coeff(actual_k), &temp.coeffRef(0));
+      }
+    }
+
+  protected:
+
+    typename VectorsType::Nested m_vectors;
+    typename CoeffsType::Nested m_coeffs;
+    bool m_trans;
+};
+
+#endif // EIGEN_HOUSEHOLDER_SEQUENCE_H
diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h
index a89305869..39edda80c 100644
--- a/Eigen/src/QR/HouseholderQR.h
+++ b/Eigen/src/QR/HouseholderQR.h
@@ -56,12 +56,13 @@ template<typename MatrixType> class HouseholderQR
       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;
     typedef Matrix<Scalar, DiagSizeAtCompileTime, 1> HCoeffsType;
     typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType;
+    typedef typename HouseholderSequence<MatrixQType,HCoeffsType>::ConjugateReturnType HouseholderSequenceType;
 
     /**
     * \brief Default Constructor.
@@ -97,7 +98,12 @@ template<typename MatrixType> class HouseholderQR
     template<typename OtherDerived, typename ResultType>
     void solve(const MatrixBase<OtherDerived>& b, ResultType *result) const;
 
-    MatrixQType matrixQ(void) const;
+    MatrixQType matrixQ() const;
+
+    HouseholderSequenceType matrixQAsHouseholderSequence() const
+    {
+      return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
+    }
 
     /** \returns a reference to the matrix where the Householder QR decomposition is stored
       * in a LAPACK-compatible way.
@@ -169,7 +175,7 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
   int rows = matrix.rows();
   int cols = matrix.cols();
   int size = std::min(rows,cols);
-  
+
   m_qr = matrix;
   m_hCoeffs.resize(size);
 
@@ -206,15 +212,7 @@ void HouseholderQR<MatrixType>::solve(
   result->resize(rows, cols);
 
   *result = b;
-  
-  Matrix<Scalar,1,MatrixType::ColsAtCompileTime> temp(cols);
-  for (int k = 0; k < cols; ++k)
-  {
-    int remainingSize = rows-k;
-
-    result->corner(BottomRight, remainingSize, cols)
-           .applyHouseholderOnTheLeft(m_qr.col(k).end(remainingSize-1), m_hCoeffs.coeff(k), &temp.coeffRef(0));
-  }
+  result->applyOnTheLeft(matrixQAsHouseholderSequence().inverse());
 
   const int rank = std::min(result->rows(), result->cols());
   m_qr.corner(TopLeft, rank, rank)
@@ -227,20 +225,7 @@ template<typename MatrixType>
 typename HouseholderQR<MatrixType>::MatrixQType HouseholderQR<MatrixType>::matrixQ() const
 {
   ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
-  // compute the product H'_0 H'_1 ... H'_n-1,
-  // where H_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();
-  int size = std::min(rows,cols);
-  MatrixQType res = MatrixQType::Identity(rows, rows);
-  Matrix<Scalar,1,MatrixType::RowsAtCompileTime> temp(rows);
-  for (int k = size-1; k >= 0; k--)
-  {
-    res.block(k, k, rows-k, rows-k)
-       .applyHouseholderOnTheLeft(m_qr.col(k).end(rows-k-1), ei_conj(m_hCoeffs.coeff(k)), &temp.coeffRef(k));
-  }
-  return res;
+  return matrixQAsHouseholderSequence();
 }
 
 #endif // EIGEN_HIDE_HEAVY_CODE