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 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 */
107 template<typename Derived>
108 template<typename EssentialPart>
applyHouseholderOnTheLeft(const EssentialPart & essential,const Scalar & tau,Scalar * workspace)109 void 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 */
144 template<typename Derived>
145 template<typename EssentialPart>
applyHouseholderOnTheRight(const EssentialPart & essential,const Scalar & tau,Scalar * workspace)146 void 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