• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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