src/Householder/Householder.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
0005 // Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_HOUSEHOLDER_H
0012 #define EIGEN_HOUSEHOLDER_H
0013
0014 namespace Eigen {
0015
0016 namespace internal {
0017 template<int n> struct decrement_size
0018 {
0019   enum {
0020     ret = n==Dynamic ? n : n-1
0021   };
0022 };
0023 }
0024
0025 /** Computes the elementary reflector H such that:
0026   * \f$ H *this = [ beta 0 ... 0]^T \f$
0027   * where the transformation H is:
0028   * \f$ H = I - tau v v^*\f$
0029   * and the vector v is:
0030   * \f$ v^T = [1 essential^T] \f$
0031   *
0032   * The essential part of the vector \c v is stored in *this.
0033   *
0034   * On output:
0035   * \param tau the scaling factor of the Householder transformation
0036   * \param beta the result of H * \c *this
0037   *
0038   * \sa MatrixBase::makeHouseholder(), MatrixBase::applyHouseholderOnTheLeft(),
0039   *     MatrixBase::applyHouseholderOnTheRight()
0040   */
0041 template<typename Derived>
0042 EIGEN_DEVICE_FUNC
0043 void MatrixBase<Derived>::makeHouseholderInPlace(Scalar& tau, RealScalar& beta)
0044 {
0045   VectorBlock<Derived, internal::decrement_size<Base::SizeAtCompileTime>::ret> essentialPart(derived(), 1, size()-1);
0046   makeHouseholder(essentialPart, tau, beta);
0047 }
0048
0049 /** Computes the elementary reflector H such that:
0050   * \f$ H *this = [ beta 0 ... 0]^T \f$
0051   * where the transformation H is:
0052   * \f$ H = I - tau v v^*\f$
0053   * and the vector v is:
0054   * \f$ v^T = [1 essential^T] \f$
0055   *
0056   * On output:
0057   * \param essential the essential part of the vector \c v
0058   * \param tau the scaling factor of the Householder transformation
0059   * \param beta the result of H * \c *this
0060   *
0061   * \sa MatrixBase::makeHouseholderInPlace(), MatrixBase::applyHouseholderOnTheLeft(),
0062   *     MatrixBase::applyHouseholderOnTheRight()
0063   */
0064 template<typename Derived>
0065 template<typename EssentialPart>
0066 EIGEN_DEVICE_FUNC
0067 void MatrixBase<Derived>::makeHouseholder(
0068   EssentialPart& essential,
0069   Scalar& tau,
0070   RealScalar& beta) const
0071 {
0072   using std::sqrt;
0073   using numext::conj;
0074
0075   EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart)
0076   VectorBlock<const Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1);
0077
0078   RealScalar tailSqNorm = size()==1 ? RealScalar(0) : tail.squaredNorm();
0079   Scalar c0 = coeff(0);
0080   const RealScalar tol = (std::numeric_limits<RealScalar>::min)();
0081
0082   if(tailSqNorm <= tol && numext::abs2(numext::imag(c0))<=tol)
0083   {
0084     tau = RealScalar(0);
0085     beta = numext::real(c0);
0086     essential.setZero();
0087   }
0088   else
0089   {
0090     beta = sqrt(numext::abs2(c0) + tailSqNorm);
0091     if (numext::real(c0)>=RealScalar(0))
0092       beta = -beta;
0093     essential = tail / (c0 - beta);
0094     tau = conj((beta - c0) / beta);
0095   }
0096 }
0097
0098 /** Apply the elementary reflector H given by
0099   * \f$ H = I - tau v v^*\f$
0100   * with
0101   * \f$ v^T = [1 essential^T] \f$
0102   * from the left to a vector or matrix.
0103   *
0104   * On input:
0105   * \param essential the essential part of the vector \c v
0106   * \param tau the scaling factor of the Householder transformation
0107   * \param workspace a pointer to working space with at least
0108   *                  this->cols() entries
0109   *
0110   * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(),
0111   *     MatrixBase::applyHouseholderOnTheRight()
0112   */
0113 template<typename Derived>
0114 template<typename EssentialPart>
0115 EIGEN_DEVICE_FUNC
0116 void MatrixBase<Derived>::applyHouseholderOnTheLeft(
0117   const EssentialPart& essential,
0118   const Scalar& tau,
0119   Scalar* workspace)
0120 {
0121   if(rows() == 1)
0122   {
0123     *this *= Scalar(1)-tau;
0124   }
0125   else if(tau!=Scalar(0))
0126   {
0127     Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace,cols());
0128     Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows()-1, cols());
0129     tmp.noalias() = essential.adjoint() * bottom;
0130     tmp += this->row(0);
0131     this->row(0) -= tau * tmp;
0132     bottom.noalias() -= tau * essential * tmp;
0133   }
0134 }
0135
0136 /** Apply the elementary reflector H given by
0137   * \f$ H = I - tau v v^*\f$
0138   * with
0139   * \f$ v^T = [1 essential^T] \f$
0140   * from the right to a vector or matrix.
0141   *
0142   * On input:
0143   * \param essential the essential part of the vector \c v
0144   * \param tau the scaling factor of the Householder transformation
0145   * \param workspace a pointer to working space with at least
0146   *                  this->rows() entries
0147   *
0148   * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(),
0149   *     MatrixBase::applyHouseholderOnTheLeft()
0150   */
0151 template<typename Derived>
0152 template<typename EssentialPart>
0153 EIGEN_DEVICE_FUNC
0154 void MatrixBase<Derived>::applyHouseholderOnTheRight(
0155   const EssentialPart& essential,
0156   const Scalar& tau,
0157   Scalar* workspace)
0158 {
0159   if(cols() == 1)
0160   {
0161     *this *= Scalar(1)-tau;
0162   }
0163   else if(tau!=Scalar(0))
0164   {
0165     Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace,rows());
0166     Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols()-1);
0167     tmp.noalias() = right * essential;
0168     tmp += this->col(0);
0169     this->col(0) -= tau * tmp;
0170     right.noalias() -= tau * tmp * essential.adjoint();
0171   }
0172 }
0173
0174 } // end namespace Eigen
0175
0176 #endif // EIGEN_HOUSEHOLDER_H