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>
42EIGEN_DEVICE_FUNC
43void 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 */
64template<typename Derived>
65template<typename EssentialPart>
66EIGEN_DEVICE_FUNC
67void 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 */
113template<typename Derived>
114template<typename EssentialPart>
115EIGEN_DEVICE_FUNC
116void 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 */
151template<typename Derived>
152template<typename EssentialPart>
153EIGEN_DEVICE_FUNC
154void 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
177

source code of qtmultimedia/src/3rdparty/eigen/Eigen/src/Householder/Householder.h