add new Householder module
diff --git a/CMakeLists.txt b/CMakeLists.txt index 15c21f5..85e8490 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt
@@ -8,8 +8,7 @@ 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 0000000..ba06bd8 --- /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 b881c09..eb29c0e 100644 --- a/Eigen/src/Core/MatrixBase.h +++ b/Eigen/src/Core/MatrixBase.h
@@ -40,6 +40,7 @@ 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 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 0000000..ce4937d --- /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 0000000..b5571f4 --- /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 2f1aa21..9c65b9c 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt
@@ -149,7 +149,7 @@ 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 0000000..d8b11c9 --- /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)) ); + } + +}