diff options
author | Benoit Jacob <jacob.benoit.1@gmail.com> | 2009-08-03 16:06:57 +0200 |
---|---|---|
committer | Benoit Jacob <jacob.benoit.1@gmail.com> | 2009-08-03 16:06:57 +0200 |
commit | d10c710b15676935ccc93faf59fe0f37e0f1b7cd (patch) | |
tree | 63d946fd7ed3a4c3acf75d6c6c3ea6820fbaad2c /Eigen | |
parent | 66ee2044ce83f508c78dafebaee12d04d9fc9019 (diff) |
add new Householder module
Diffstat (limited to 'Eigen')
-rw-r--r-- | Eigen/Householder | 24 | ||||
-rw-r--r-- | Eigen/src/Core/MatrixBase.h | 14 | ||||
-rw-r--r-- | Eigen/src/Householder/CMakeLists.txt | 6 | ||||
-rw-r--r-- | Eigen/src/Householder/Householder.h | 88 |
4 files changed, 132 insertions, 0 deletions
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 |