Add blocking inside HouseholderQR based on code from Vincent Lejeune.

This is all internal stuff for now.

Eigen/Householder

@@ -17,6 +17,7 @@ namespace Eigen {
 
 #include "src/Householder/Householder.h"
 #include "src/Householder/HouseholderSequence.h"
+#include "src/Householder/BlockHouseholder.h"
 
 } // namespace Eigen
 

Eigen/src/Householder/BlockHouseholder.h

@@ -0,0 +1,76 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2010 Vincent Lejeune
+// Copyright (C) 2010 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_BLOCK_HOUSEHOLDER_H
+#define EIGEN_BLOCK_HOUSEHOLDER_H
+
+// This file contains some helper function to deal with block householder reflectors
+
+/** \internal */
+template<typename TriangularFactorType,typename VectorsType,typename CoeffsType>
+void ei_make_block_householder_triangular_factor(TriangularFactorType& triFactor, const VectorsType& vectors, const CoeffsType& hCoeffs)
+{
+  typedef typename TriangularFactorType::Index Index;
+  typedef typename VectorsType::Scalar Scalar;
+  const Index nbVecs = vectors.cols();
+  ei_assert(triFactor.rows() == nbVecs && triFactor.cols() == nbVecs && vectors.rows()>=nbVecs);
+
+  for(Index i = 0; i < nbVecs; i++)
+  {
+    Index rs = vectors.rows() - i;
+    Scalar Vii = vectors(i,i);
+    vectors.const_cast_derived().coeffRef(i,i) = Scalar(1);
+    triFactor.col(i).head(i).noalias() = vectors.block(i, 0, rs, i).adjoint()
+                                       * vectors.col(i).tail(rs);
+    triFactor.col(i).head(i) *= -hCoeffs(i);
+    vectors.const_cast_derived().coeffRef(i, i) = Vii;
+    // FIXME add .noalias() once the triangular product can work inplace
+    triFactor.col(i).head(i) = triFactor.block(0,0,i,i).template triangularView<Upper>()
+                             * triFactor.col(i).head(i);
+    triFactor(i,i) = hCoeffs(i);
+  }
+}
+
+/** \internal */
+template<typename MatrixType,typename VectorsType,typename CoeffsType>
+void ei_apply_block_householder_on_the_left(MatrixType& mat, const VectorsType& vectors, const CoeffsType& hCoeffs)
+{
+  typedef typename MatrixType::Index Index;
+  enum { TFactorSize = MatrixType::ColsAtCompileTime };
+  Index nbVecs = vectors.cols();
+  Matrix<typename MatrixType::Scalar, TFactorSize, TFactorSize> T(nbVecs,nbVecs);
+  ei_make_block_householder_triangular_factor(T, vectors, hCoeffs);
+
+  const TriangularView<VectorsType, UnitLower>& V(vectors);
+
+  // A -= V T V^* A
+  Matrix<typename MatrixType::Scalar,Dynamic,Dynamic> tmp = V.adjoint() * mat;
+  // FIXME add .noalias() once the triangular product can work inplace
+  tmp = T.template triangularView<Upper>().adjoint() * tmp;
+  mat.noalias() -= V * tmp;
+}
+
+
+#endif // EIGEN_BLOCK_HOUSEHOLDER_H

Eigen/src/QR/HouseholderQR.h

@@ -1,8 +1,9 @@
 // This file is part of Eigen, a lightweight C++ template library
 // for linear algebra.
 //
-// Copyright (C) 2008-2009 Gael Guennebaud <g.gael@free.fr>
+// Copyright (C) 2008-2010 Gael Guennebaud <g.gael@free.fr>
 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
+// Copyright (C) 2010 Vincent Lejeune
 //
 // Eigen is free software; you can redistribute it and/or
 // modify it under the terms of the GNU Lesser General Public
@@ -229,6 +230,59 @@ void ei_householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typena
   }
 }
 
+/** \internal */
+template<typename MatrixQR, typename HCoeffs>
+void ei_householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs,
+                                       typename MatrixQR::Index maxBlockSize=32,
+                                       typename MatrixQR::Scalar* tempData = 0)
+{
+  typedef typename MatrixQR::Index Index;
+  typedef typename MatrixQR::Scalar Scalar;
+  typedef typename MatrixQR::RealScalar RealScalar;
+  typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
+
+  Index rows = mat.rows();
+  Index cols = mat.cols();
+  Index size = std::min(rows, cols);
+
+  typedef Matrix<Scalar,Dynamic,1,0,MatrixQR::MaxColsAtCompileTime,1> TempType;
+  TempType tempVector;
+  if(tempData==0)
+  {
+    tempVector.resize(cols);
+    tempData = tempVector.data();
+  }
+
+  Index blockSize = std::min(maxBlockSize,size);
+
+  int k = 0;
+  for (k = 0; k < size; k += blockSize)
+  {
+    Index bs = std::min(size-k,blockSize);  // actual size of the block
+    Index tcols = cols - k - bs;            // trailing columns
+    Index brows = rows-k;                   // rows of the block
+
+    // partition the matrix:
+    //        A00 | A01 | A02
+    // mat  = A10 | A11 | A12
+    //        A20 | A21 | A22
+    // and performs the qr dec of [A11^T A12^T]^T
+    // and update [A21^T A22^T]^T using level 3 operations.
+    // Finally, the algorithm continue on A22
+
+    BlockType A11_21 = mat.block(k,k,brows,bs);
+    Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
+
+    ei_householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
+
+    if(tcols)
+    {
+      BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
+      ei_apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint());
+    }
+  }
+}
+
 template<typename MatrixType>
 HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
 {
@@ -241,7 +295,7 @@ HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType&
 
   m_temp.resize(cols);
 
-  ei_householder_qr_inplace_unblocked(m_qr, m_hCoeffs, m_temp.data());
+  ei_householder_qr_inplace_blocked(m_qr, m_hCoeffs, 48, m_temp.data());
 
   m_isInitialized = true;
   return *this;

test/qr.cpp

@@ -36,10 +36,10 @@ template<typename MatrixType> void qr(const MatrixType& m)
 
   MatrixType a = MatrixType::Random(rows,cols);
   HouseholderQR<MatrixType> qrOfA(a);
-  
+
   MatrixQType q = qrOfA.householderQ();
   VERIFY_IS_UNITARY(q);
-  
+
   MatrixType r = qrOfA.matrixQR().template triangularView<Upper>();
   VERIFY_IS_APPROX(a, qrOfA.householderQ() * r);
 }
@@ -112,8 +112,8 @@ template<typename MatrixType> void qr_verify_assert()
 void test_qr()
 {
   for(int i = 0; i < g_repeat; i++) {
-   CALL_SUBTEST_1( qr(MatrixXf(47,40)) );
+   CALL_SUBTEST_1( qr(MatrixXf(ei_random<int>(1,200),ei_random<int>(1,200))) );
-   CALL_SUBTEST_2( qr(MatrixXcd(17,7)) );
+   CALL_SUBTEST_2( qr(MatrixXcd(ei_random<int>(1,200),ei_random<int>(1,200))) );
    CALL_SUBTEST_3(( qr_fixedsize<Matrix<float,3,4>, 2 >() ));
    CALL_SUBTEST_4(( qr_fixedsize<Matrix<double,6,2>, 4 >() ));
    CALL_SUBTEST_5(( qr_fixedsize<Matrix<double,2,5>, 7 >() ));