1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2010 Benoit Jacob <[email protected]>
5 // Copyright (C) 2009 Gael Guennebaud <[email protected]>
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
makeHouseholderInPlace(Scalar & tau,RealScalar & beta)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
makeHouseholder(EssentialPart & essential,Scalar & tau,RealScalar & beta)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
applyHouseholderOnTheLeft(const EssentialPart & essential,const Scalar & tau,Scalar * workspace)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
applyHouseholderOnTheRight(const EssentialPart & essential,const Scalar & tau,Scalar * workspace)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
177