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