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>
makeHouseholderInPlace(Scalar & tau,RealScalar & beta)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>
makeHouseholder(EssentialPart & essential,Scalar & tau,RealScalar & beta)65 void 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
79 if(tailSqNorm == RealScalar(0) && numext::imag(c0)==RealScalar(0))
80 {
81 tau = RealScalar(0);
82 beta = numext::real(c0);
83 essential.setZero();
84 }
85 else
86 {
87 beta = sqrt(numext::abs2(c0) + tailSqNorm);
88 if (numext::real(c0)>=RealScalar(0))
89 beta = -beta;
90 essential = tail / (c0 - beta);
91 tau = conj((beta - c0) / beta);
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>
applyHouseholderOnTheLeft(const EssentialPart & essential,const Scalar & tau,Scalar * workspace)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>
applyHouseholderOnTheRight(const EssentialPart & essential,const Scalar & tau,Scalar * workspace)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
172