blob: 96b1daf5fa939cb84794e656e971e235a87b0ce6 [file] [log] [blame]
Googler45874d82019-08-21 12:06:47 -07001// 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
Rasmus Munk Larsen66c640f2023-08-22 14:16:58 -070014// IWYU pragma: private
C. Antonio Sanchezdcfd7702021-10-26 11:41:24 -070015#include "./InternalHeaderCheck.h"
16
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080017namespace Eigen {
Googler45874d82019-08-21 12:06:47 -070018
19namespace internal {
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080020template <int n>
21struct decrement_size {
22 enum { ret = n == Dynamic ? n : n - 1 };
Googler45874d82019-08-21 12:06:47 -070023};
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080024} // namespace internal
Googler45874d82019-08-21 12:06:47 -070025
26/** Computes the elementary reflector H such that:
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080027 * \f$ H *this = [ beta 0 ... 0]^T \f$
28 * where the transformation H is:
29 * \f$ H = I - tau v v^*\f$
30 * and the vector v is:
31 * \f$ v^T = [1 essential^T] \f$
32 *
33 * The essential part of the vector \c v is stored in *this.
34 *
35 * On output:
36 * \param tau the scaling factor of the Householder transformation
37 * \param beta the result of H * \c *this
38 *
39 * \sa MatrixBase::makeHouseholder(), MatrixBase::applyHouseholderOnTheLeft(),
40 * MatrixBase::applyHouseholderOnTheRight()
41 */
42template <typename Derived>
43EIGEN_DEVICE_FUNC void MatrixBase<Derived>::makeHouseholderInPlace(Scalar& tau, RealScalar& beta) {
44 VectorBlock<Derived, internal::decrement_size<Base::SizeAtCompileTime>::ret> essentialPart(derived(), 1, size() - 1);
Googler45874d82019-08-21 12:06:47 -070045 makeHouseholder(essentialPart, tau, beta);
46}
47
48/** Computes the elementary reflector H such that:
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080049 * \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 */
63template <typename Derived>
64template <typename EssentialPart>
65EIGEN_DEVICE_FUNC void MatrixBase<Derived>::makeHouseholder(EssentialPart& essential, Scalar& tau,
66 RealScalar& beta) const {
Googler45874d82019-08-21 12:06:47 -070067 using numext::conj;
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080068 using numext::sqrt;
69
Googler45874d82019-08-21 12:06:47 -070070 EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart)
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080071 VectorBlock<const Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size() - 1);
72
73 RealScalar tailSqNorm = size() == 1 ? RealScalar(0) : tail.squaredNorm();
Googler45874d82019-08-21 12:06:47 -070074 Scalar c0 = coeff(0);
75 const RealScalar tol = (std::numeric_limits<RealScalar>::min)();
76
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080077 if (tailSqNorm <= tol && numext::abs2(numext::imag(c0)) <= tol) {
Googler45874d82019-08-21 12:06:47 -070078 tau = RealScalar(0);
79 beta = numext::real(c0);
80 essential.setZero();
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080081 } else {
Googler45874d82019-08-21 12:06:47 -070082 beta = sqrt(numext::abs2(c0) + tailSqNorm);
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080083 if (numext::real(c0) >= RealScalar(0)) beta = -beta;
Googler45874d82019-08-21 12:06:47 -070084 essential = tail / (c0 - beta);
85 tau = conj((beta - c0) / beta);
86 }
87}
88
89/** Apply the elementary reflector H given by
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -080090 * \f$ H = I - tau v v^*\f$
91 * with
92 * \f$ v^T = [1 essential^T] \f$
93 * from the left to a vector or matrix.
94 *
95 * On input:
96 * \param essential the essential part of the vector \c v
97 * \param tau the scaling factor of the Householder transformation
98 * \param workspace a pointer to working space with at least
99 * this->cols() entries
100 *
101 * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(),
102 * MatrixBase::applyHouseholderOnTheRight()
103 */
104template <typename Derived>
105template <typename EssentialPart>
106EIGEN_DEVICE_FUNC void MatrixBase<Derived>::applyHouseholderOnTheLeft(const EssentialPart& essential, const Scalar& tau,
107 Scalar* workspace) {
108 if (rows() == 1) {
109 *this *= Scalar(1) - tau;
110 } else if (!numext::is_exactly_zero(tau)) {
111 Map<typename internal::plain_row_type<PlainObject>::type> tmp(workspace, cols());
112 Block<Derived, EssentialPart::SizeAtCompileTime, Derived::ColsAtCompileTime> bottom(derived(), 1, 0, rows() - 1,
113 cols());
Googler45874d82019-08-21 12:06:47 -0700114 tmp.noalias() = essential.adjoint() * bottom;
115 tmp += this->row(0);
116 this->row(0) -= tau * tmp;
117 bottom.noalias() -= tau * essential * tmp;
118 }
119}
120
121/** Apply the elementary reflector H given by
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -0800122 * \f$ H = I - tau v v^*\f$
123 * with
124 * \f$ v^T = [1 essential^T] \f$
125 * from the right to a vector or matrix.
126 *
127 * On input:
128 * \param essential the essential part of the vector \c v
129 * \param tau the scaling factor of the Householder transformation
130 * \param workspace a pointer to working space with at least
131 * this->rows() entries
132 *
133 * \sa MatrixBase::makeHouseholder(), MatrixBase::makeHouseholderInPlace(),
134 * MatrixBase::applyHouseholderOnTheLeft()
135 */
136template <typename Derived>
137template <typename EssentialPart>
138EIGEN_DEVICE_FUNC void MatrixBase<Derived>::applyHouseholderOnTheRight(const EssentialPart& essential,
139 const Scalar& tau, Scalar* workspace) {
140 if (cols() == 1) {
141 *this *= Scalar(1) - tau;
142 } else if (!numext::is_exactly_zero(tau)) {
143 Map<typename internal::plain_col_type<PlainObject>::type> tmp(workspace, rows());
144 Block<Derived, Derived::RowsAtCompileTime, EssentialPart::SizeAtCompileTime> right(derived(), 0, 1, rows(),
145 cols() - 1);
Googler45874d82019-08-21 12:06:47 -0700146 tmp.noalias() = right * essential;
147 tmp += this->col(0);
148 this->col(0) -= tau * tmp;
149 right.noalias() -= tau * tmp * essential.adjoint();
150 }
151}
152
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -0800153} // end namespace Eigen
Googler45874d82019-08-21 12:06:47 -0700154
Rasmus Munk Larsen2434cfd2023-12-06 12:02:41 -0800155#endif // EIGEN_HOUSEHOLDER_H