Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 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 | |
Rasmus Munk Larsen | 66c640f | 2023-08-22 14:16:58 -0700 | [diff] [blame] | 14 | // IWYU pragma: private |
C. Antonio Sanchez | dcfd770 | 2021-10-26 11:41:24 -0700 | [diff] [blame] | 15 | #include "./InternalHeaderCheck.h" |
| 16 | |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 17 | namespace Eigen { |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 18 | |
| 19 | namespace internal { |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 20 | template <int n> |
| 21 | struct decrement_size { |
| 22 | enum { ret = n == Dynamic ? n : n - 1 }; |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 23 | }; |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 24 | } // namespace internal |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 25 | |
| 26 | /** Computes the elementary reflector H such that: |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 27 | * \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 | */ |
| 42 | template <typename Derived> |
| 43 | EIGEN_DEVICE_FUNC void MatrixBase<Derived>::makeHouseholderInPlace(Scalar& tau, RealScalar& beta) { |
| 44 | VectorBlock<Derived, internal::decrement_size<Base::SizeAtCompileTime>::ret> essentialPart(derived(), 1, size() - 1); |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 45 | makeHouseholder(essentialPart, tau, beta); |
| 46 | } |
| 47 | |
| 48 | /** Computes the elementary reflector H such that: |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 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> |
| 65 | EIGEN_DEVICE_FUNC void MatrixBase<Derived>::makeHouseholder(EssentialPart& essential, Scalar& tau, |
| 66 | RealScalar& beta) const { |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 67 | using numext::conj; |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 68 | using numext::sqrt; |
| 69 | |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 70 | EIGEN_STATIC_ASSERT_VECTOR_ONLY(EssentialPart) |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 71 | VectorBlock<const Derived, EssentialPart::SizeAtCompileTime> tail(derived(), 1, size() - 1); |
| 72 | |
| 73 | RealScalar tailSqNorm = size() == 1 ? RealScalar(0) : tail.squaredNorm(); |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 74 | Scalar c0 = coeff(0); |
| 75 | const RealScalar tol = (std::numeric_limits<RealScalar>::min)(); |
| 76 | |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 77 | if (tailSqNorm <= tol && numext::abs2(numext::imag(c0)) <= tol) { |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 78 | tau = RealScalar(0); |
| 79 | beta = numext::real(c0); |
| 80 | essential.setZero(); |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 81 | } else { |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 82 | beta = sqrt(numext::abs2(c0) + tailSqNorm); |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 83 | if (numext::real(c0) >= RealScalar(0)) beta = -beta; |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 84 | essential = tail / (c0 - beta); |
| 85 | tau = conj((beta - c0) / beta); |
| 86 | } |
| 87 | } |
| 88 | |
| 89 | /** Apply the elementary reflector H given by |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 90 | * \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 | */ |
| 104 | template <typename Derived> |
| 105 | template <typename EssentialPart> |
| 106 | EIGEN_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()); |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 114 | 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 Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 122 | * \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 | */ |
| 136 | template <typename Derived> |
| 137 | template <typename EssentialPart> |
| 138 | EIGEN_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); |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 146 | 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 Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 153 | } // end namespace Eigen |
Googler | 45874d8 | 2019-08-21 12:06:47 -0700 | [diff] [blame] | 154 | |
Rasmus Munk Larsen | 2434cfd | 2023-12-06 12:02:41 -0800 | [diff] [blame] | 155 | #endif // EIGEN_HOUSEHOLDER_H |