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