* HouseholderSequence:

  * be aware of number of actual householder vectors
    (optimization in non-full-rank case, no behavior change)
  * fix applyThisOnTheRight, it was using k instead of actual_k
* QR: rename matrixQ() to householderQ() where applicable

7 files changed, 43 insertions, 39 deletions

diff --git a/Eigen/src/Householder/HouseholderSequence.h b/Eigen/src/Householder/HouseholderSequence.h
index afeb758b3..1159f03ed 100644
--- a/Eigen/src/Householder/HouseholderSequence.h
+++ b/Eigen/src/Householder/HouseholderSequence.h
@@ -69,31 +69,35 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
 
     typedef HouseholderSequence<VectorsType,
       typename ei_meta_if<NumTraits<Scalar>::IsComplex,
-        typename ei_cleantype<typename CoeffsType::ConjugateReturnType>::type,
+        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)
+      : m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(v.diagonalSize())
+    {}
+
+    HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors)
+      : m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(actualVectors)
     {}
 
     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); }
+    { return HouseholderSequence(m_vectors, m_coeffs, !m_trans, m_actualVectors); }
 
     ConjugateReturnType conjugate() const
-    { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans); }
+    { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), m_trans, m_actualVectors); }
 
     ConjugateReturnType adjoint() const
-    { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans); }
+    { return ConjugateReturnType(m_vectors, m_coeffs.conjugate(), !m_trans, m_actualVectors); }
 
     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 vecs = m_actualVectors;
       int length = m_vectors.rows();
       dst.setIdentity();
       Matrix<Scalar,1,DestType::RowsAtCompileTime> temp(dst.rows());
@@ -111,22 +115,22 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
     /** \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());
+      int vecs = m_actualVectors; // number of householder vectors
+      int length = m_vectors.rows(); // size of the largest householder vector
+      Matrix<Scalar,1,Dest::RowsAtCompileTime> 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));
+        dst.corner(BottomRight, dst.rows(), length-actual_k)
+           .applyHouseholderOnTheRight(m_vectors.col(actual_k).end(length-actual_k-1), m_coeffs.coeff(actual_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
+      int vecs = m_actualVectors; // 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)
       {
@@ -156,6 +160,7 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
     typename VectorsType::Nested m_vectors;
     typename CoeffsType::Nested m_coeffs;
     bool m_trans;
+    int m_actualVectors;
 
 private:
   HouseholderSequence& operator=(const HouseholderSequence&);
@@ -167,4 +172,10 @@ HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsTyp
   return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans);
 }
 
+template<typename VectorsType, typename CoeffsType>
+HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors)
+{
+  return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans, actualVectors);
+}
+
 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H
diff --git a/Eigen/src/QR/ColPivHouseholderQR.h b/Eigen/src/QR/ColPivHouseholderQR.h
index 9c1e56dc6..1b2472a99 100644
--- a/Eigen/src/QR/ColPivHouseholderQR.h
+++ b/Eigen/src/QR/ColPivHouseholderQR.h
@@ -105,7 +105,7 @@ template<typename _MatrixType> class ColPivHouseholderQR
       return ei_solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
     }
 
-    HouseholderSequenceType matrixQ(void) const;
+    HouseholderSequenceType householderQ(void) const;
 
     /** \returns a reference to the matrix where the Householder QR decomposition is stored
       */
@@ -447,7 +447,8 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
     c.applyOnTheLeft(householderSequence(
       dec().matrixQR(),
       dec().hCoeffs(),
-      true
+      true,
+      dec().nonzeroPivots()
     ));
 
     dec().matrixQR()
@@ -470,10 +471,11 @@ struct ei_solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
 
 /** \returns the matrix Q as a sequence of householder transformations */
 template<typename MatrixType>
-typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>::matrixQ() const
+typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
+  ::householderQ() const
 {
   ei_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
-  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate(), false);
+  return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate(), false, m_nonzero_pivots);
 }
 
 #endif // EIGEN_HIDE_HEAVY_CODE
diff --git a/Eigen/src/QR/HouseholderQR.h b/Eigen/src/QR/HouseholderQR.h
index 8d842d129..95496b943 100644
--- a/Eigen/src/QR/HouseholderQR.h
+++ b/Eigen/src/QR/HouseholderQR.h
@@ -104,11 +104,10 @@ template<typename _MatrixType> class HouseholderQR
       ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
       return ei_solve_retval<HouseholderQR, Rhs>(*this, b.derived());
     }
-    
-    MatrixQType matrixQ() const;
 
-    HouseholderSequenceType matrixQAsHouseholderSequence() const
+    HouseholderSequenceType householderQ() const
     {
+      ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
       return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
     }
 
@@ -240,14 +239,6 @@ struct ei_solve_retval<HouseholderQR<_MatrixType>, Rhs>
   }
 };
 
-/** \returns the matrix Q */
-template<typename MatrixType>
-typename HouseholderQR<MatrixType>::MatrixQType HouseholderQR<MatrixType>::matrixQ() const
-{
-  ei_assert(m_isInitialized && "HouseholderQR is not initialized.");
-  return matrixQAsHouseholderSequence();
-}
-
 #endif // EIGEN_HIDE_HEAVY_CODE
 
 /** \return the Householder QR decomposition of \c *this.
diff --git a/test/geo_hyperplane.cpp b/test/geo_hyperplane.cpp
index 3cf5655c2..1fd1e281a 100644
--- a/test/geo_hyperplane.cpp
+++ b/test/geo_hyperplane.cpp
@@ -64,7 +64,7 @@ template<typename HyperplaneType> void hyperplane(const HyperplaneType& _plane)
   // transform
   if (!NumTraits<Scalar>::IsComplex)
   {
-    MatrixType rot = MatrixType::Random(dim,dim).householderQr().matrixQ();
+    MatrixType rot = MatrixType::Random(dim,dim).householderQr().householderQ();
     DiagonalMatrix<Scalar,HyperplaneType::AmbientDimAtCompileTime> scaling(VectorType::Random());
     Translation<Scalar,HyperplaneType::AmbientDimAtCompileTime> translation(VectorType::Random());
 
diff --git a/test/main.h b/test/main.h
index 3fded75c8..9624838ff 100644
--- a/test/main.h
+++ b/test/main.h
@@ -360,7 +360,7 @@ void createRandomMatrixOfRank(int desired_rank, int rows, int cols, MatrixType&
 
   HouseholderQR<MatrixAType> qra(a);
   HouseholderQR<MatrixBType> qrb(b);
-  m = qra.matrixQ() * d * qrb.matrixQ();
+  m = qra.householderQ() * d * qrb.householderQ();
 }
 
 } // end namespace Eigen
diff --git a/test/qr.cpp b/test/qr.cpp
index 90b5c4446..3848ce0a5 100644
--- a/test/qr.cpp
+++ b/test/qr.cpp
@@ -38,13 +38,13 @@ template<typename MatrixType> void qr(const MatrixType& m)
   HouseholderQR<MatrixType> qrOfA(a);
   MatrixType r = qrOfA.matrixQR();
   
-  MatrixQType q = qrOfA.matrixQ();
+  MatrixQType q = qrOfA.householderQ();
   VERIFY_IS_UNITARY(q);
   
   // FIXME need better way to construct trapezoid
   for(int i = 0; i < rows; i++) for(int j = 0; j < cols; j++) if(i>j) r(i,j) = Scalar(0);
 
-  VERIFY_IS_APPROX(a, qrOfA.matrixQ() * r);
+  VERIFY_IS_APPROX(a, qrOfA.householderQ() * r);
 }
 
 template<typename MatrixType, int Cols2> void qr_fixedsize()
@@ -58,7 +58,7 @@ template<typename MatrixType, int Cols2> void qr_fixedsize()
   // FIXME need better way to construct trapezoid
   for(int i = 0; i < Rows; i++) for(int j = 0; j < Cols; j++) if(i>j) r(i,j) = Scalar(0);
 
-  VERIFY_IS_APPROX(m1, qr.matrixQ() * r);
+  VERIFY_IS_APPROX(m1, qr.householderQ() * r);
 
   Matrix<Scalar,Cols,Cols2> m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
   Matrix<Scalar,Rows,Cols2> m3 = m1*m2;
@@ -93,7 +93,7 @@ template<typename MatrixType> void qr_invertible()
   m1.setZero();
   for(int i = 0; i < size; i++) m1(i,i) = ei_random<Scalar>();
   RealScalar absdet = ei_abs(m1.diagonal().prod());
-  m3 = qr.matrixQ(); // get a unitary
+  m3 = qr.householderQ(); // get a unitary
   m1 = m3 * m1 * m3;
   qr.compute(m1);
   VERIFY_IS_APPROX(absdet, qr.absDeterminant());
@@ -107,7 +107,7 @@ template<typename MatrixType> void qr_verify_assert()
   HouseholderQR<MatrixType> qr;
   VERIFY_RAISES_ASSERT(qr.matrixQR())
   VERIFY_RAISES_ASSERT(qr.solve(tmp))
-  VERIFY_RAISES_ASSERT(qr.matrixQ())
+  VERIFY_RAISES_ASSERT(qr.householderQ())
   VERIFY_RAISES_ASSERT(qr.absDeterminant())
   VERIFY_RAISES_ASSERT(qr.logAbsDeterminant())
 }
diff --git a/test/qr_colpivoting.cpp b/test/qr_colpivoting.cpp
index 48b6de3f5..8b56cd296 100644
--- a/test/qr_colpivoting.cpp
+++ b/test/qr_colpivoting.cpp
@@ -44,7 +44,7 @@ template<typename MatrixType> void qr()
   VERIFY(!qr.isInvertible());
   VERIFY(!qr.isSurjective());
 
-  MatrixQType q = qr.matrixQ();
+  MatrixQType q = qr.householderQ();
   VERIFY_IS_UNITARY(q);
 
   MatrixType r = qr.matrixQR().template triangularView<UpperTriangular>();
@@ -73,7 +73,7 @@ template<typename MatrixType, int Cols2> void qr_fixedsize()
   VERIFY(!qr.isSurjective());
 
   Matrix<Scalar,Rows,Cols> r = qr.matrixQR().template triangularView<UpperTriangular>();
-  Matrix<Scalar,Rows,Cols> c = qr.matrixQ() * r * qr.colsPermutation().inverse();
+  Matrix<Scalar,Rows,Cols> c = qr.householderQ() * r * qr.colsPermutation().inverse();
   VERIFY_IS_APPROX(m1, c);
 
   Matrix<Scalar,Cols,Cols2> m2 = Matrix<Scalar,Cols,Cols2>::Random(Cols,Cols2);
@@ -109,7 +109,7 @@ template<typename MatrixType> void qr_invertible()
   m1.setZero();
   for(int i = 0; i < size; i++) m1(i,i) = ei_random<Scalar>();
   RealScalar absdet = ei_abs(m1.diagonal().prod());
-  m3 = qr.matrixQ(); // get a unitary
+  m3 = qr.householderQ(); // get a unitary
   m1 = m3 * m1 * m3;
   qr.compute(m1);
   VERIFY_IS_APPROX(absdet, qr.absDeterminant());
@@ -123,7 +123,7 @@ template<typename MatrixType> void qr_verify_assert()
   ColPivHouseholderQR<MatrixType> qr;
   VERIFY_RAISES_ASSERT(qr.matrixQR())
   VERIFY_RAISES_ASSERT(qr.solve(tmp))
-  VERIFY_RAISES_ASSERT(qr.matrixQ())
+  VERIFY_RAISES_ASSERT(qr.householderQ())
   VERIFY_RAISES_ASSERT(qr.dimensionOfKernel())
   VERIFY_RAISES_ASSERT(qr.isInjective())
   VERIFY_RAISES_ASSERT(qr.isSurjective())