1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_HOUSEHOLDER_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_HOUSEHOLDER_H 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<int n> struct decrement_size 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ret = n==Dynamic ? n : n-1 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Computes the elementary reflector H such that: 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ H *this = [ beta 0 ... 0]^T \f$ 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * where the transformation H is: 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ H = I - tau v v^*\f$ 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and the vector v is: 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ v^T = [1 essential^T] \f$ 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The essential part of the vector \c v is stored in *this. 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * On output: 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param tau the scaling factor of the Householder transformation 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param beta the result of H * \c *this 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::makeHouseholder(), MatrixBase::applyHouseholderOnTheLeft(), 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * MatrixBase::applyHouseholderOnTheRight() 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixBase<Derived>::makeHouseholderInPlace(Scalar& tau, RealScalar& beta) 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorBlock<Derived, internal::decrement_size<Base::SizeAtCompileTime>::ret> essentialPart(derived(), 1, size()-1); 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath makeHouseholder(essentialPart, tau, beta); 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Computes the elementary reflector H such that: 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ H *this = [ beta 0 ... 0]^T \f$ 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * where the transformation H is: 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ H = I - tau v v^*\f$ 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and the vector v is: 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ v^T = [1 essential^T] \f$ 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * On output: 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param essential the essential part of the vector \c v 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param tau the scaling factor of the Householder transformation 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param beta the result of H * \c *this 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::makeHouseholderInPlace(), MatrixBase::applyHouseholderOnTheLeft(), 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * MatrixBase::applyHouseholderOnTheRight() 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename EssentialPart> 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixBase<Derived>::makeHouseholder( 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EssentialPart& essential, 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar& tau, 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar& beta) const 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart) 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorBlock<const Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1); 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar tailSqNorm = size()==1 ? RealScalar(0) : tail.squaredNorm(); 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar c0 = coeff(0); 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(tailSqNorm == RealScalar(0) && internal::imag(c0)==RealScalar(0)) 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tau = RealScalar(0); 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath beta = internal::real(c0); 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath essential.setZero(); 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath beta = internal::sqrt(internal::abs2(c0) + tailSqNorm); 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (internal::real(c0)>=RealScalar(0)) 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath beta = -beta; 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath essential = tail / (c0 - beta); 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tau = internal::conj((beta - c0) / beta); 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Apply the elementary reflector H given by 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ H = I - tau v v^*\f$ 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ v^T = [1 essential^T] \f$ 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * from the left to a vector or matrix. 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * On input: 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param essential the essential part of the vector \c v 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param tau the scaling factor of the Householder transformation 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param workspace a pointer to working space with at least 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * this->cols() * essential.size() entries 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(), 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * MatrixBase::applyHouseholderOnTheRight() 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename EssentialPart> 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixBase<Derived>::applyHouseholderOnTheLeft( 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const EssentialPart& essential, 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar& tau, 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* workspace) 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(rows() == 1) 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *this *= Scalar(1)-tau; 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace,cols()); 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows()-1, cols()); 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp.noalias() = essential.adjoint() * bottom; 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp += this->row(0); 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath this->row(0) -= tau * tmp; 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bottom.noalias() -= tau * essential * tmp; 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** Apply the elementary reflector H given by 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ H = I - tau v v^*\f$ 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ v^T = [1 essential^T] \f$ 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * from the right to a vector or matrix. 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * On input: 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param essential the essential part of the vector \c v 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param tau the scaling factor of the Householder transformation 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param workspace a pointer to working space with at least 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * this->cols() * essential.size() entries 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(), 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * MatrixBase::applyHouseholderOnTheLeft() 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Derived> 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename EssentialPart> 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid MatrixBase<Derived>::applyHouseholderOnTheRight( 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const EssentialPart& essential, 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar& tau, 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar* workspace) 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(cols() == 1) 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath *this *= Scalar(1)-tau; 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace,rows()); 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols()-1); 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp.noalias() = right * essential.conjugate(); 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath tmp += this->col(0); 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath this->col(0) -= tau * tmp; 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath right.noalias() -= tau * tmp * essential.transpose(); 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_HOUSEHOLDER_H 169