Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 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 | void MatrixBase<Derived>::makeHouseholderInPlace(Scalar& tau, RealScalar& beta) |
| 43 | { |
| 44 | VectorBlock<Derived, internal::decrement_size<Base::SizeAtCompileTime>::ret> essentialPart(derived(), 1, size()-1); |
| 45 | makeHouseholder(essentialPart, tau, beta); |
| 46 | } |
| 47 | |
| 48 | /** Computes the elementary reflector H such that: |
| 49 | * \f$ H *this = [ beta 0 ... 0]^T \f$ |
| 50 | * where the transformation H is: |
| 51 | * \f$ H = I - tau v v^*\f$ |
| 52 | * and the vector v is: |
| 53 | * \f$ v^T = [1 essential^T] \f$ |
| 54 | * |
| 55 | * On output: |
| 56 | * \param essential the essential part of the vector \c v |
| 57 | * \param tau the scaling factor of the Householder transformation |
| 58 | * \param beta the result of H * \c *this |
| 59 | * |
| 60 | * \sa MatrixBase::makeHouseholderInPlace(), MatrixBase::applyHouseholderOnTheLeft(), |
| 61 | * MatrixBase::applyHouseholderOnTheRight() |
| 62 | */ |
| 63 | template<typename Derived> |
| 64 | template<typename EssentialPart> |
| 65 | void MatrixBase<Derived>::makeHouseholder( |
| 66 | EssentialPart& essential, |
| 67 | Scalar& tau, |
| 68 | RealScalar& beta) const |
| 69 | { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 70 | using std::sqrt; |
| 71 | using numext::conj; |
| 72 | |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 73 | EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart) |
| 74 | VectorBlock<const Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1); |
| 75 | |
| 76 | RealScalar tailSqNorm = size()==1 ? RealScalar(0) : tail.squaredNorm(); |
| 77 | Scalar c0 = coeff(0); |
| 78 | |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 79 | if(tailSqNorm == RealScalar(0) && numext::imag(c0)==RealScalar(0)) |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 80 | { |
| 81 | tau = RealScalar(0); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 82 | beta = numext::real(c0); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 83 | essential.setZero(); |
| 84 | } |
| 85 | else |
| 86 | { |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 87 | beta = sqrt(numext::abs2(c0) + tailSqNorm); |
| 88 | if (numext::real(c0)>=RealScalar(0)) |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 89 | beta = -beta; |
| 90 | essential = tail / (c0 - beta); |
Carlos Hernandez | 7faaa9f | 2014-08-05 17:53:32 -0700 | [diff] [blame] | 91 | tau = conj((beta - c0) / beta); |
Narayan Kamath | c981c48 | 2012-11-02 10:59:05 +0000 | [diff] [blame] | 92 | } |
| 93 | } |
| 94 | |
| 95 | /** Apply the elementary reflector H given by |
| 96 | * \f$ H = I - tau v v^*\f$ |
| 97 | * with |
| 98 | * \f$ v^T = [1 essential^T] \f$ |
| 99 | * from the left to a vector or matrix. |
| 100 | * |
| 101 | * On input: |
| 102 | * \param essential the essential part of the vector \c v |
| 103 | * \param tau the scaling factor of the Householder transformation |
| 104 | * \param workspace a pointer to working space with at least |
| 105 | * this->cols() * essential.size() entries |
| 106 | * |
| 107 | * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(), |
| 108 | * MatrixBase::applyHouseholderOnTheRight() |
| 109 | */ |
| 110 | template<typename Derived> |
| 111 | template<typename EssentialPart> |
| 112 | void MatrixBase<Derived>::applyHouseholderOnTheLeft( |
| 113 | const EssentialPart& essential, |
| 114 | const Scalar& tau, |
| 115 | Scalar* workspace) |
| 116 | { |
| 117 | if(rows() == 1) |
| 118 | { |
| 119 | *this *= Scalar(1)-tau; |
| 120 | } |
| 121 | else |
| 122 | { |
| 123 | Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace,cols()); |
| 124 | Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows()-1, cols()); |
| 125 | tmp.noalias() = essential.adjoint() * bottom; |
| 126 | tmp += this->row(0); |
| 127 | this->row(0) -= tau * tmp; |
| 128 | bottom.noalias() -= tau * essential * tmp; |
| 129 | } |
| 130 | } |
| 131 | |
| 132 | /** Apply the elementary reflector H given by |
| 133 | * \f$ H = I - tau v v^*\f$ |
| 134 | * with |
| 135 | * \f$ v^T = [1 essential^T] \f$ |
| 136 | * from the right to a vector or matrix. |
| 137 | * |
| 138 | * On input: |
| 139 | * \param essential the essential part of the vector \c v |
| 140 | * \param tau the scaling factor of the Householder transformation |
| 141 | * \param workspace a pointer to working space with at least |
| 142 | * this->cols() * essential.size() entries |
| 143 | * |
| 144 | * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(), |
| 145 | * MatrixBase::applyHouseholderOnTheLeft() |
| 146 | */ |
| 147 | template<typename Derived> |
| 148 | template<typename EssentialPart> |
| 149 | void MatrixBase<Derived>::applyHouseholderOnTheRight( |
| 150 | const EssentialPart& essential, |
| 151 | const Scalar& tau, |
| 152 | Scalar* workspace) |
| 153 | { |
| 154 | if(cols() == 1) |
| 155 | { |
| 156 | *this *= Scalar(1)-tau; |
| 157 | } |
| 158 | else |
| 159 | { |
| 160 | Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace,rows()); |
| 161 | Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols()-1); |
| 162 | tmp.noalias() = right * essential.conjugate(); |
| 163 | tmp += this->col(0); |
| 164 | this->col(0) -= tau * tmp; |
| 165 | right.noalias() -= tau * tmp * essential.transpose(); |
| 166 | } |
| 167 | } |
| 168 | |
| 169 | } // end namespace Eigen |
| 170 | |
| 171 | #endif // EIGEN_HOUSEHOLDER_H |