Introduce third template parameter to HouseholderSequence: int Side.

When it's OnTheRight, we read householder vectors as rows above the diagonal.
With unit test. The use case will be bidiagonalization.

Eigen/src/Core/util/ForwardDeclarations.h

@@ -134,7 +134,7 @@ template<typename MatrixType> class SVD;
 template<typename MatrixType, unsigned int Options = 0> class JacobiSVD;
 template<typename MatrixType, int UpLo = Lower> class LLT;
 template<typename MatrixType> class LDLT;
-template<typename VectorsType, typename CoeffsType> class HouseholderSequence;
+template<typename VectorsType, typename CoeffsType, int Side=OnTheLeft> class HouseholderSequence;
 template<typename Scalar>     class PlanarRotation;
 
 // Geometry module:

Eigen/src/Householder/Householder.h

@@ -1,7 +1,7 @@
 // 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) 2010 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

Eigen/src/Householder/HouseholderSequence.h

@@ -2,6 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2009 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -48,31 +49,61 @@
   * \sa MatrixBase::applyOnTheLeft(), MatrixBase::applyOnTheRight()
   */
 
-template<typename VectorsType, typename CoeffsType>
+template<typename VectorsType, typename CoeffsType, int Side>
-struct ei_traits<HouseholderSequence<VectorsType,CoeffsType> >
+struct ei_traits<HouseholderSequence<VectorsType,CoeffsType,Side> >
 {
   typedef typename VectorsType::Scalar Scalar;
   enum {
-    RowsAtCompileTime = ei_traits<VectorsType>::RowsAtCompileTime,
+    RowsAtCompileTime = Side==OnTheLeft ? ei_traits<VectorsType>::RowsAtCompileTime
-    ColsAtCompileTime = ei_traits<VectorsType>::RowsAtCompileTime,
+                                        : ei_traits<VectorsType>::ColsAtCompileTime,
-    MaxRowsAtCompileTime = ei_traits<VectorsType>::MaxRowsAtCompileTime,
+    ColsAtCompileTime = RowsAtCompileTime,
-    MaxColsAtCompileTime = ei_traits<VectorsType>::MaxRowsAtCompileTime,
+    MaxRowsAtCompileTime = Side==OnTheLeft ? ei_traits<VectorsType>::MaxRowsAtCompileTime
+                                           : ei_traits<VectorsType>::MaxColsAtCompileTime,
+    MaxColsAtCompileTime = MaxRowsAtCompileTime,
     Flags = 0
   };
 };
 
-template<typename VectorsType, typename CoeffsType> class HouseholderSequence
+template<typename VectorsType, typename CoeffsType, int Side>
-  : public AnyMatrixBase<HouseholderSequence<VectorsType,CoeffsType> >
+struct ei_hseq_side_dependent_impl
+{
+  typedef Block<VectorsType, Dynamic, 1> EssentialVectorType;
+  typedef HouseholderSequence<VectorsType, CoeffsType, OnTheLeft> HouseholderSequenceType;
+  static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, int k)
+  {
+    const int start = k+1+h.m_shift;
+    return Block<VectorsType,Dynamic,1>(h.m_vectors, start, k, h.rows()-start, 1);
+  }
+};
+
+template<typename VectorsType, typename CoeffsType>
+struct ei_hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight>
+{
+  typedef Transpose<Block<VectorsType, 1, Dynamic> > EssentialVectorType;
+  typedef HouseholderSequence<VectorsType, CoeffsType, OnTheRight> HouseholderSequenceType;
+  static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, int k)
+  {
+    const int start = k+1+h.m_shift;
+    return Block<VectorsType,1,Dynamic>(h.m_vectors, k, start, 1, h.rows()-start).transpose();
+  }
+};
+
+template<typename VectorsType, typename CoeffsType, int Side> class HouseholderSequence
+  : public AnyMatrixBase<HouseholderSequence<VectorsType,CoeffsType,Side> >
 {
     typedef typename VectorsType::Scalar Scalar;
-    typedef Block<VectorsType, Dynamic, 1> EssentialVectorType;
+    typedef typename ei_hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::EssentialVectorType
+            EssentialVectorType;
     
   public:
 
-    typedef HouseholderSequence<VectorsType,
+    typedef HouseholderSequence<
+      VectorsType,
       typename ei_meta_if<NumTraits<Scalar>::IsComplex,
         typename ei_cleantype<typename CoeffsType::ConjugateReturnType>::type,
-        CoeffsType>::ret> ConjugateReturnType;
+        CoeffsType>::ret,
+      Side
+    > ConjugateReturnType;
 
     HouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans = false)
       : m_vectors(v), m_coeffs(h), m_trans(trans), m_actualVectors(v.diagonalSize()),
@@ -85,14 +116,13 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
     {
     }
 
-    int rows() const { return m_vectors.rows(); }
+    int rows() const { return Side==OnTheLeft ? m_vectors.rows() : m_vectors.cols(); }
-    int cols() const { return m_vectors.rows(); }
+    int cols() const { return rows(); }
 
     const EssentialVectorType essentialVector(int k) const
     {
-      ei_assert(k>= 0 && k < m_actualVectors);
+      ei_assert(k >= 0 && k < m_actualVectors);
-      const int start = k+1+m_shift;
+      return ei_hseq_side_dependent_impl<VectorsType,CoeffsType,Side>::essentialVector(*this, k);
-      return Block<VectorsType,Dynamic,1>(m_vectors, start, k, rows()-start, 1);
     }
 
     HouseholderSequence transpose() const
@@ -164,6 +194,8 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
       return res;
     }
 
+    template<typename _VectorsType, typename _CoeffsType, int _Side> friend struct ei_hseq_side_dependent_impl;
+
   protected:
     typename VectorsType::Nested m_vectors;
     typename CoeffsType::Nested m_coeffs;
@@ -175,13 +207,25 @@ template<typename VectorsType, typename CoeffsType> class HouseholderSequence
 template<typename VectorsType, typename CoeffsType>
 HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h, bool trans=false)
 {
-  return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans);
+  return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h, trans);
 }
 
 template<typename VectorsType, typename CoeffsType>
 HouseholderSequence<VectorsType,CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors, int shift)
 {
-  return HouseholderSequence<VectorsType,CoeffsType>(v, h, trans, actualVectors, shift);
+  return HouseholderSequence<VectorsType,CoeffsType,OnTheLeft>(v, h, trans, actualVectors, shift);
+}
+
+template<typename VectorsType, typename CoeffsType>
+HouseholderSequence<VectorsType,CoeffsType> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans=false)
+{
+  return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h, trans);
+}
+
+template<typename VectorsType, typename CoeffsType>
+HouseholderSequence<VectorsType,CoeffsType> rightHouseholderSequence(const VectorsType& v, const CoeffsType& h, bool trans, int actualVectors, int shift)
+{
+  return HouseholderSequence<VectorsType,CoeffsType,OnTheRight>(v, h, trans, actualVectors, shift);
 }
 
 #endif // EIGEN_HOUSEHOLDER_SEQUENCE_H

test/householder.cpp

@@ -106,27 +106,15 @@ template<typename MatrixType> void householder(const MatrixType& m)
   m5.block(shift,0,brows,cols).template triangularView<StrictlyLower>().setZero();
   VERIFY_IS_APPROX(hseq * m5, m1); // test applying hseq directly
   m3 = hseq;
-  VERIFY_IS_APPROX(m3*m5, m1); // test evaluating hseq to a dense matrix, then applying
+  VERIFY_IS_APPROX(m3 * m5, m1); // test evaluating hseq to a dense matrix, then applying
 
-#if 0
   // test householder sequence on the right with a shift
 
-  TMatrixType tm1 = m1.transpose();
   TMatrixType tm2 = m2.transpose();
-  
+  HouseholderSequence<TMatrixType, HCoeffsVectorType, OnTheRight> rhseq(tm2, hc, false, hc.size(), shift);
-  int bcols = cols - shift;
+  VERIFY_IS_APPROX(rhseq * m5, m1); // test applying rhseq directly
-  VBlockMatrixType vbm = 
+  m3 = rhseq;
-  HouseholderQR<HBlockMatrixType> qr(hbm);
+  VERIFY_IS_APPROX(m3 * m5, m1); // test evaluating rhseq to a dense matrix, then applying
-  m2 = m1;
-  m2.block(shift,0,brows,cols) = qr.matrixQR();
-  HCoeffsVectorType hc = qr.hCoeffs().conjugate();
-  HouseholderSequence<MatrixType, HCoeffsVectorType> hseq(m2, hc, false, hc.size(), shift);
-  MatrixType m5 = m2;
-  m5.block(shift,0,brows,cols).template triangularView<StrictlyLower>().setZero();
-  VERIFY_IS_APPROX(hseq * m5, m1); // test applying hseq directly
-  m3 = hseq;
-  VERIFY_IS_APPROX(m3*m5, m1); // test evaluating hseq to a dense matrix, then applying
-#endif
 }
 
 void test_householder()