Householder.h
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 
41 template<typename Derived>
42 void MatrixBase<Derived>::makeHouseholderInPlace(Scalar& tau, RealScalar& beta)
43 {
45  makeHouseholder(essentialPart, tau, beta);
46 }
47 
63 template<typename Derived>
64 template<typename EssentialPart>
66  EssentialPart& essential,
67  Scalar& tau,
68  RealScalar& beta) const
69 {
70  EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart)
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 
107 template<typename Derived>
108 template<typename EssentialPart>
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  {
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 
144 template<typename Derived>
145 template<typename EssentialPart>
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  {
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