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
14namespace Eigen {
15
16namespace internal {
17template<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  */
41template<typename Derived>
42void 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  */
63template<typename Derived>
64template<typename EssentialPart>
65void MatrixBase<Derived>::makeHouseholder(
66  EssentialPart& essential,
67  Scalar& tau,
68  RealScalar& beta) const
69{
70  EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart)
71  VectorBlock<const Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size()-1);
72
73  RealScalar tailSqNorm = size()==1 ? RealScalar(0) : tail.squaredNorm();
74  Scalar c0 = coeff(0);
75
76  if(tailSqNorm == RealScalar(0) && internal::imag(c0)==RealScalar(0))
77  {
78    tau = RealScalar(0);
79    beta = internal::real(c0);
80    essential.setZero();
81  }
82  else
83  {
84    beta = internal::sqrt(internal::abs2(c0) + tailSqNorm);
85    if (internal::real(c0)>=RealScalar(0))
86      beta = -beta;
87    essential = tail / (c0 - beta);
88    tau = internal::conj((beta - c0) / beta);
89  }
90}
91
92/** Apply the elementary reflector H given by
93  * \f$ H = I - tau v v^*\f$
94  * with
95  * \f$ v^T = [1 essential^T] \f$
96  * from the left to a vector or matrix.
97  *
98  * On input:
99  * \param essential the essential part of the vector \c v
100  * \param tau the scaling factor of the Householder transformation
101  * \param workspace a pointer to working space with at least
102  *                  this->cols() * essential.size() entries
103  *
104  * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(),
105  *     MatrixBase::applyHouseholderOnTheRight()
106  */
107template<typename Derived>
108template<typename EssentialPart>
109void MatrixBase<Derived>::applyHouseholderOnTheLeft(
110  const EssentialPart& essential,
111  const Scalar& tau,
112  Scalar* workspace)
113{
114  if(rows() == 1)
115  {
116    *this *= Scalar(1)-tau;
117  }
118  else
119  {
120    Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace,cols());
121    Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows()-1, cols());
122    tmp.noalias() = essential.adjoint() * bottom;
123    tmp += this->row(0);
124    this->row(0) -= tau * tmp;
125    bottom.noalias() -= tau * essential * tmp;
126  }
127}
128
129/** Apply the elementary reflector H given by
130  * \f$ H = I - tau v v^*\f$
131  * with
132  * \f$ v^T = [1 essential^T] \f$
133  * from the right to a vector or matrix.
134  *
135  * On input:
136  * \param essential the essential part of the vector \c v
137  * \param tau the scaling factor of the Householder transformation
138  * \param workspace a pointer to working space with at least
139  *                  this->cols() * essential.size() entries
140  *
141  * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(),
142  *     MatrixBase::applyHouseholderOnTheLeft()
143  */
144template<typename Derived>
145template<typename EssentialPart>
146void MatrixBase<Derived>::applyHouseholderOnTheRight(
147  const EssentialPart& essential,
148  const Scalar& tau,
149  Scalar* workspace)
150{
151  if(cols() == 1)
152  {
153    *this *= Scalar(1)-tau;
154  }
155  else
156  {
157    Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace,rows());
158    Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(), cols()-1);
159    tmp.noalias() = right * essential.conjugate();
160    tmp += this->col(0);
161    this->col(0) -= tau * tmp;
162    right.noalias() -= tau * tmp * essential.transpose();
163  }
164}
165
166} // end namespace Eigen
167
168#endif // EIGEN_HOUSEHOLDER_H
169