add new Householder module

7 files changed, 223 insertions, 3 deletions

diff --git a/CMakeLists.txt b/CMakeLists.txt
index 15c21f591..85e8490f1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -8,8 +8,7 @@ string(REGEX MATCH "define *EIGEN_MAJOR_VERSION ([0-9]*)" _eigen2_major_version_
 set(EIGEN2_MAJOR_VERSION "${CMAKE_MATCH_1}")
 string(REGEX MATCH "define *EIGEN_MINOR_VERSION ([0-9]*)" _eigen2_minor_version_match "${_eigen2_version_header}")
 set(EIGEN2_MINOR_VERSION "${CMAKE_MATCH_1}")
-# TODO find a way to automatically remove the unstable suffix
-set(EIGEN_VERSION_NUMBER ${EIGEN2_WORLD_VERSION}.${EIGEN2_MAJOR_VERSION}.${EIGEN2_MINOR_VERSION}-unstable)
+set(EIGEN_VERSION_NUMBER ${EIGEN2_WORLD_VERSION}.${EIGEN2_MAJOR_VERSION}.${EIGEN2_MINOR_VERSION})
 
 # if the mercurial program is absent, this will leave the EIGEN_HG_REVISION string empty,
 # but won't stop CMake.
diff --git a/Eigen/Householder b/Eigen/Householder
new file mode 100644
index 000000000..ba06bd8fb
--- /dev/null
+++ b/Eigen/Householder
@@ -0,0 +1,24 @@
+#ifndef EIGEN_HOUSEHOLDER_MODULE_H
+#define EIGEN_HOUSEHOLDER_MODULE_H
+
+#include "Core"
+
+#include "src/Core/util/DisableMSVCWarnings.h"
+
+namespace Eigen {
+
+/** \defgroup Householder_Module Householder module
+  * This module provides householder transformations.
+  *
+  * \code
+  * #include <Eigen/Householder>
+  * \endcode
+  */
+
+#include "src/Householder/Householder.h"
+
+} // namespace Eigen
+
+#include "src/Core/util/EnableMSVCWarnings.h"
+
+#endif // EIGEN_HOUSEHOLDER_MODULE_H
diff --git a/Eigen/src/Core/MatrixBase.h b/Eigen/src/Core/MatrixBase.h
index b881c09c3..eb29c0ebe 100644
--- a/Eigen/src/Core/MatrixBase.h
+++ b/Eigen/src/Core/MatrixBase.h
@@ -40,6 +40,7 @@ template<typename Derived> struct AnyMatrixBase
   Derived& derived() { return *static_cast<Derived*>(this); }
   const Derived& derived() const { return *static_cast<const Derived*>(this); }
 };
+
 /** Common base class for all classes T such that there are overloaded operator* allowing to
   * multiply a MatrixBase by a T on both sides.
   *
@@ -749,6 +750,19 @@ template<typename Derived> class MatrixBase
     template<typename Derived1, typename Derived2>
     Derived& lazyAssign(const SparseProduct<Derived1,Derived2,DenseTimeSparseProduct>& product);
 
+////////// Householder module ///////////
+
+    template<typename EssentialPart>
+    void makeHouseholder(EssentialPart *essential,
+                         RealScalar *beta) const;
+    template<typename EssentialPart>
+    void applyHouseholderOnTheLeft(const EssentialPart& essential,
+                                   const RealScalar& beta);
+    template<typename EssentialPart>
+    void applyHouseholderOnTheRight(const EssentialPart& essential,
+                                    const RealScalar& beta);
+    
+
     #ifdef EIGEN_MATRIXBASE_PLUGIN
     #include EIGEN_MATRIXBASE_PLUGIN
     #endif
diff --git a/Eigen/src/Householder/CMakeLists.txt b/Eigen/src/Householder/CMakeLists.txt
new file mode 100644
index 000000000..ce4937db0
--- /dev/null
+++ b/Eigen/src/Householder/CMakeLists.txt
@@ -0,0 +1,6 @@
+FILE(GLOB Eigen_Householder_SRCS "*.h")
+
+INSTALL(FILES
+  ${Eigen_Householder_SRCS}
+  DESTINATION ${INCLUDE_INSTALL_DIR}/Eigen/src/Householder COMPONENT Devel
+  )
diff --git a/Eigen/src/Householder/Householder.h b/Eigen/src/Householder/Householder.h
new file mode 100644
index 000000000..b5571f40d
--- /dev/null
+++ b/Eigen/src/Householder/Householder.h
@@ -0,0 +1,88 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 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
+// 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_H
+#define EIGEN_HOUSEHOLDER_H
+
+template<int n> struct ei_decrement_size
+{
+  enum {
+    ret = (n==1 || n==Dynamic) ? n : n-1
+  };
+};
+
+template<typename Derived>
+template<typename EssentialPart>
+void MatrixBase<Derived>::makeHouseholder(
+  EssentialPart *essential,
+  RealScalar *beta) const
+{
+  EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart)
+  RealScalar _squaredNorm = squaredNorm();
+  Scalar c0;
+  if(ei_abs2(coeff(0)) <= ei_abs2(precision<Scalar>()) * _squaredNorm)
+  {
+    c0 = ei_sqrt(_squaredNorm);
+  }
+  else
+  {
+    Scalar sign = coeff(0) / ei_abs(coeff(0));
+    c0 = coeff(0) + sign * ei_sqrt(_squaredNorm);
+  }
+  *essential = end(size()-1) / c0; // FIXME take advantage of fixed size
+  const RealScalar c0abs2 = ei_abs2(c0);
+  *beta = RealScalar(2) * c0abs2 / (c0abs2 + _squaredNorm - ei_abs2(coeff(0)));
+}
+
+template<typename Derived>
+template<typename EssentialPart>
+void MatrixBase<Derived>::applyHouseholderOnTheLeft(
+  const EssentialPart& essential,
+  const RealScalar& beta)
+{
+  Matrix<Scalar, 1, ColsAtCompileTime, PlainMatrixType::Options, 1, MaxColsAtCompileTime> tmp(cols());
+  tmp = row(0) + essential.adjoint() * block(1,0,rows()-1,cols());
+  // FIXME take advantage of fixed size
+  // FIXME play with lazy()
+  // FIXME maybe not a good idea to use matrix product
+  row(0) -= beta * tmp;
+  block(1,0,rows()-1,cols()) -= beta * essential * tmp;
+}
+
+template<typename Derived>
+template<typename EssentialPart>
+void MatrixBase<Derived>::applyHouseholderOnTheRight(
+  const EssentialPart& essential,
+  const RealScalar& beta)
+{
+  Matrix<Scalar, RowsAtCompileTime, 1, PlainMatrixType::Options, MaxRowsAtCompileTime, 1> tmp(rows());
+  tmp = col(0) + block(0,1,rows(),cols()-1) * essential.conjugate();
+  // FIXME take advantage of fixed size
+  // FIXME play with lazy()
+  // FIXME maybe not a good idea to use matrix product
+  col(0) -= beta * tmp;
+  block(0,1,rows(),cols()-1) -= beta * tmp * essential.transpose();
+}
+
+#endif // EIGEN_HOUSEHOLDER_H
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index 2f1aa21a3..9c65b9c51 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -149,7 +149,7 @@ ei_add_test(sparse_basic)
 ei_add_test(sparse_product)
 ei_add_test(sparse_solvers " " "${SPARSE_LIBS}")
 ei_add_test(umeyama)
-
+ei_add_test(householder)
 
 ei_add_property(EIGEN_TESTING_SUMMARY "CXX:               ${CMAKE_CXX_COMPILER}\n")
 if(CMAKE_COMPILER_IS_GNUCXX)
diff --git a/test/householder.cpp b/test/householder.cpp
new file mode 100644
index 000000000..d8b11c9f1
--- /dev/null
+++ b/test/householder.cpp
@@ -0,0 +1,89 @@
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2009 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
+// 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/>.
+
+#include "main.h"
+#include <Eigen/Householder>
+
+template<typename MatrixType> void householder(const MatrixType& m)
+{
+  /* this test covers the following files:
+     Householder.h
+  */
+  int rows = m.rows();
+  int cols = m.cols();
+
+  typedef typename MatrixType::Scalar Scalar;
+  typedef typename NumTraits<Scalar>::Real RealScalar;
+  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, 1> VectorType;
+  typedef Matrix<Scalar, ei_decrement_size<MatrixType::RowsAtCompileTime>::ret, 1> EssentialVectorType;
+  typedef Matrix<Scalar, MatrixType::RowsAtCompileTime, MatrixType::RowsAtCompileTime> SquareMatrixType;
+  
+  RealScalar beta;
+  EssentialVectorType essential;
+
+  VectorType v1 = VectorType::Random(rows), v2;
+  v2 = v1;
+  v1.makeHouseholder(&essential, &beta);
+  v1.applyHouseholderOnTheLeft(essential,beta);
+  
+  VERIFY_IS_APPROX(v1.norm(), v2.norm());
+  VERIFY_IS_MUCH_SMALLER_THAN(v1.end(rows-1).norm(), v1.norm());
+  v1 = VectorType::Random(rows);
+  v2 = v1;
+  v1.applyHouseholderOnTheLeft(essential,beta);
+  VERIFY_IS_APPROX(v1.norm(), v2.norm());
+  
+  MatrixType m1(rows, cols),
+             m2(rows, cols);
+
+  v1 = VectorType::Random(rows);
+  m1.colwise() = v1;
+  m2 = m1;
+  m1.col(0).makeHouseholder(&essential, &beta);
+  m1.applyHouseholderOnTheLeft(essential,beta);
+  VERIFY_IS_APPROX(m1.norm(), m2.norm());
+  VERIFY_IS_MUCH_SMALLER_THAN(m1.block(1,0,rows-1,cols).norm(), m1.norm());
+  
+  v1 = VectorType::Random(rows);
+  SquareMatrixType m3(rows,rows), m4(rows,rows);
+  m3.rowwise() = v1.transpose();
+  m4 = m3;
+  m3.row(0).makeHouseholder(&essential, &beta);
+  m3.applyHouseholderOnTheRight(essential,beta);
+  VERIFY_IS_APPROX(m3.norm(), m4.norm());
+  VERIFY_IS_MUCH_SMALLER_THAN(m3.block(0,1,rows,rows-1).norm(), m3.norm());
+}
+
+void test_householder()
+{
+  for(int i = 0; i < g_repeat; i++) {
+    CALL_SUBTEST( householder(Matrix<double,2,2>()) );
+    CALL_SUBTEST( householder(Matrix<float,2,3>()) );
+    CALL_SUBTEST( householder(Matrix<double,3,5>()) );
+    CALL_SUBTEST( householder(Matrix<float,4,4>()) );
+    CALL_SUBTEST( householder(MatrixXd(10,12)) );
+    CALL_SUBTEST( householder(MatrixXcf(16,17)) );
+  }
+
+}