HouseholderSequence.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2010 Benoit Jacob <jacob.benoit.1@gmail.com>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
// SPDX-License-Identifier: MPL-2.0
 
#ifndef EIGEN_HOUSEHOLDER_SEQUENCE_H
#define EIGEN_HOUSEHOLDER_SEQUENCE_H
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {

 
namespace internal {
 
template <typename VectorsType, typename CoeffsType, int Side>
struct traits<HouseholderSequence<VectorsType, CoeffsType, Side> > {
  using Scalar = typename VectorsType::Scalar;
  using StorageIndex = typename VectorsType::StorageIndex;
  using StorageKind = typename VectorsType::StorageKind;
  enum {
    RowsAtCompileTime =
        Side == OnTheLeft ? traits<VectorsType>::RowsAtCompileTime : traits<VectorsType>::ColsAtCompileTime,
    ColsAtCompileTime = RowsAtCompileTime,
    MaxRowsAtCompileTime =
        Side == OnTheLeft ? traits<VectorsType>::MaxRowsAtCompileTime : traits<VectorsType>::MaxColsAtCompileTime,
    MaxColsAtCompileTime = MaxRowsAtCompileTime,
    Flags = 0
  };
};
 
struct HouseholderSequenceShape {};
 
template <typename VectorsType, typename CoeffsType, int Side>
struct evaluator_traits<HouseholderSequence<VectorsType, CoeffsType, Side> >
    : public evaluator_traits_base<HouseholderSequence<VectorsType, CoeffsType, Side> > {
  using Shape = HouseholderSequenceShape;
};
 
template <typename VectorsType, typename CoeffsType, int Side>
struct hseq_side_dependent_impl {
  using EssentialVectorType = Block<const VectorsType, Dynamic, 1>;
  using HouseholderSequenceType = HouseholderSequence<VectorsType, CoeffsType, OnTheLeft>;
  static EIGEN_DEVICE_FUNC inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k) {
    Index start = k + 1 + h.m_shift;
    return Block<const VectorsType, Dynamic, 1>(h.m_vectors, start, k, h.rows() - start, 1);
  }
};
 
template <typename VectorsType, typename CoeffsType>
struct hseq_side_dependent_impl<VectorsType, CoeffsType, OnTheRight> {
  using EssentialVectorType = Transpose<Block<const VectorsType, 1, Dynamic>>;
  using HouseholderSequenceType = HouseholderSequence<VectorsType, CoeffsType, OnTheRight>;
  static inline const EssentialVectorType essentialVector(const HouseholderSequenceType& h, Index k) {
    Index start = k + 1 + h.m_shift;
    return Block<const VectorsType, 1, Dynamic>(h.m_vectors, k, start, 1, h.rows() - start).transpose();
  }
};
 
template <typename OtherScalarType, typename MatrixType>
struct matrix_type_times_scalar_type {
  using ResultScalar = typename ScalarBinaryOpTraits<OtherScalarType, typename MatrixType::Scalar>::ReturnType;
  using Type =
      Matrix<ResultScalar, MatrixType::RowsAtCompileTime, MatrixType::ColsAtCompileTime,
             (MatrixType::MaxRowsAtCompileTime == 1 && MatrixType::MaxColsAtCompileTime != 1) ? RowMajor : ColMajor,
             MatrixType::MaxRowsAtCompileTime, MatrixType::MaxColsAtCompileTime>;
};
 
// A permutation has no scalar type; its product with a Householder sequence takes the sequence's.
template <typename Scalar, typename PermutationType>
struct permutation_type_as_matrix_type {
  using Type = Matrix<Scalar, PermutationType::RowsAtCompileTime, PermutationType::ColsAtCompileTime, 0,
                      PermutationType::MaxRowsAtCompileTime, PermutationType::MaxColsAtCompileTime>;
};
 
// A diagonal or triangular operand is assigned to the result directly when the scalar types agree. The triangular
// assignment kernel does not convert scalars, so a promoted result evaluates the operand and casts it, as the dense
// product casts its operand.
template <typename ResultType, typename OtherDerived>
std::enable_if_t<std::is_same<typename ResultType::Scalar, typename OtherDerived::Scalar>::value>
assign_householder_operand(ResultType& res, const EigenBase<OtherDerived>& other) {
  res = other.derived();
}
template <typename ResultType, typename OtherDerived>
std::enable_if_t<!std::is_same<typename ResultType::Scalar, typename OtherDerived::Scalar>::value>
assign_householder_operand(ResultType& res, const EigenBase<OtherDerived>& other) {
  res = other.derived().toDenseMatrix().template cast<typename ResultType::Scalar>();
}
 
}  // end namespace internal
 
template <typename VectorsType, typename CoeffsType, int Side>
class HouseholderSequence : public EigenBase<HouseholderSequence<VectorsType, CoeffsType, Side> > {
  using EssentialVectorType =
      typename internal::hseq_side_dependent_impl<VectorsType, CoeffsType, Side>::EssentialVectorType;
 
 public:
  enum {
    RowsAtCompileTime = internal::traits<HouseholderSequence>::RowsAtCompileTime,
    ColsAtCompileTime = internal::traits<HouseholderSequence>::ColsAtCompileTime,
    MaxRowsAtCompileTime = internal::traits<HouseholderSequence>::MaxRowsAtCompileTime,
    MaxColsAtCompileTime = internal::traits<HouseholderSequence>::MaxColsAtCompileTime
  };
  using Scalar = typename internal::traits<HouseholderSequence>::Scalar;
 
  using ConjugateReturnType = HouseholderSequence<
      std::conditional_t<NumTraits<Scalar>::IsComplex,
                         internal::remove_all_t<typename VectorsType::ConjugateReturnType>, VectorsType>,
      std::conditional_t<NumTraits<Scalar>::IsComplex, internal::remove_all_t<typename CoeffsType::ConjugateReturnType>,
                         CoeffsType>,
      Side>;
 
  using AdjointReturnType = HouseholderSequence<
      VectorsType,
      std::conditional_t<NumTraits<Scalar>::IsComplex, internal::remove_all_t<typename CoeffsType::ConjugateReturnType>,
                         CoeffsType>,
      Side>;
 
  using TransposeReturnType = HouseholderSequence<
      std::conditional_t<NumTraits<Scalar>::IsComplex,
                         internal::remove_all_t<typename VectorsType::ConjugateReturnType>, VectorsType>,
      CoeffsType, Side>;
 
  using ConstHouseholderSequence =
      HouseholderSequence<std::add_const_t<VectorsType>, std::add_const_t<CoeffsType>, Side>;

  EIGEN_DEVICE_FUNC HouseholderSequence(const VectorsType& v, const CoeffsType& h)
      : m_vectors(v), m_coeffs(h), m_reverse(false), m_length(v.diagonalSize()), m_shift(0) {}

  EIGEN_DEVICE_FUNC HouseholderSequence(const HouseholderSequence& other)
      : m_vectors(other.m_vectors),
        m_coeffs(other.m_coeffs),
        m_reverse(other.m_reverse),
        m_length(other.m_length),
        m_shift(other.m_shift) {}

  EIGEN_DEVICE_FUNC constexpr Index rows() const noexcept {
    return Side == OnTheLeft ? m_vectors.rows() : m_vectors.cols();
  }

  EIGEN_DEVICE_FUNC constexpr Index cols() const noexcept { return rows(); }

  EIGEN_DEVICE_FUNC const EssentialVectorType essentialVector(Index k) const {
    eigen_assert(k >= 0 && k < m_length);
    return internal::hseq_side_dependent_impl<VectorsType, CoeffsType, Side>::essentialVector(*this, k);
  }

  TransposeReturnType transpose() const {
    return TransposeReturnType(m_vectors.conjugate(), m_coeffs)
        .setReverseFlag(!m_reverse)
        .setLength(m_length)
        .setShift(m_shift);
  }

  ConjugateReturnType conjugate() const {
    return ConjugateReturnType(m_vectors.conjugate(), m_coeffs.conjugate())
        .setReverseFlag(m_reverse)
        .setLength(m_length)
        .setShift(m_shift);
  }

  template <bool Cond>
  EIGEN_DEVICE_FUNC inline std::conditional_t<Cond, ConjugateReturnType, ConstHouseholderSequence> conjugateIf() const {
    using ReturnType = std::conditional_t<Cond, ConjugateReturnType, ConstHouseholderSequence>;
    return ReturnType(m_vectors.template conjugateIf<Cond>(), m_coeffs.template conjugateIf<Cond>());
  }

  AdjointReturnType adjoint() const {
    return AdjointReturnType(m_vectors, m_coeffs.conjugate())
        .setReverseFlag(!m_reverse)
        .setLength(m_length)
        .setShift(m_shift);
  }

  AdjointReturnType inverse() const { return adjoint(); }

  template <typename DestType>
  inline EIGEN_DEVICE_FUNC void evalTo(DestType& dst) const {
    Matrix<Scalar, DestType::RowsAtCompileTime, 1, AutoAlign | ColMajor, DestType::MaxRowsAtCompileTime, 1> workspace(
        rows());
    evalTo(dst, workspace);
  }

  template <typename Dest, typename Workspace>
  EIGEN_DEVICE_FUNC void evalTo(Dest& dst, Workspace& workspace) const {
    workspace.resize(rows());
    Index vecs = m_length;
    if (internal::is_same_dense(dst, m_vectors)) {
      // in-place
      dst.diagonal().setOnes();
      dst.template triangularView<StrictlyUpper>().setZero();
      for (Index k = vecs - 1; k >= 0; --k) {
        Index cornerSize = rows() - k - m_shift;
        if (m_reverse)
          dst.bottomRightCorner(cornerSize, cornerSize)
              .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
        else
          dst.bottomRightCorner(cornerSize, cornerSize)
              .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
 
        // clear the off diagonal vector
        dst.col(k).tail(rows() - k - 1).setZero();
      }
      // clear the remaining columns if needed
      for (Index k = 0; k < cols() - vecs; ++k) dst.col(k).tail(rows() - k - 1).setZero();
    } else if (m_length > BlockSize) {
      dst.setIdentity(rows(), rows());
      applyThisOnTheLeft(dst, workspace, true);
    } else {
      dst.setIdentity(rows(), rows());
      for (Index k = vecs - 1; k >= 0; --k) {
        Index cornerSize = rows() - k - m_shift;
        if (m_reverse)
          dst.bottomRightCorner(cornerSize, cornerSize)
              .applyHouseholderOnTheRight(essentialVector(k), m_coeffs.coeff(k), workspace.data());
        else
          dst.bottomRightCorner(cornerSize, cornerSize)
              .applyHouseholderOnTheLeft(essentialVector(k), m_coeffs.coeff(k), workspace.data());
      }
    }
  }

  template <typename Dest>
  inline void applyThisOnTheRight(Dest& dst) const {
    Matrix<Scalar, 1, Dest::RowsAtCompileTime, RowMajor, 1, Dest::MaxRowsAtCompileTime> workspace(dst.rows());
    applyThisOnTheRight(dst, workspace);
  }

  template <typename Dest, typename Workspace>
  inline void applyThisOnTheRight(Dest& dst, Workspace& workspace) const {
    // Use the block path when the reflectors are long enough for GEMM to outperform GEMV.
    // The threshold on rows() (the reflector length) is higher than for the left-side path because
    // the right-side block application has more overhead from the tmp = mat * V product layout.
    if (m_length >= BlockSize && rows() - m_shift >= 4 * BlockSize) {
      applyBlockOnTheRight(dst);
    } else {
      workspace.resize(dst.rows());
      for (Index k = 0; k < m_length; ++k) {
        Index actual_k = m_reverse ? m_length - k - 1 : k;
        dst.rightCols(rows() - m_shift - actual_k)
            .applyHouseholderOnTheRight(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
      }
    }
  }
 
 private:
  template <typename Dest>
  EIGEN_DONT_INLINE void applyBlockOnTheRight(Dest& dst) const {
    // Make sure we have at least 2 useful blocks, otherwise it is point-less:
    Index blockSize = m_length < Index(2 * BlockSize) ? (m_length + 1) / 2 : Index(BlockSize);
    for (Index i = 0; i < m_length; i += blockSize) {
      // Right-side application processes blocks in opposite order to left-side:
      // forward (m_reverse=false): first block first; reversed: last block first.
      Index end = m_reverse ? m_length - i : (std::min)(m_length, i + blockSize);
      Index k = m_reverse ? (std::max)(Index(0), end - blockSize) : i;
      Index bs = end - k;
      Index start = k + m_shift;
 
      using SubVectorsType = Block<internal::remove_all_t<VectorsType>, Dynamic, Dynamic>;
      SubVectorsType sub_vecs1(m_vectors.const_cast_derived(), Side == OnTheRight ? k : start,
                               Side == OnTheRight ? start : k, Side == OnTheRight ? bs : m_vectors.rows() - start,
                               Side == OnTheRight ? m_vectors.cols() - start : bs);
      std::conditional_t<Side == OnTheRight, Transpose<SubVectorsType>, SubVectorsType&> sub_vecs(sub_vecs1);
 
      Index dstCols = rows() - m_shift - k;
      auto sub_dst = dst.rightCols(dstCols);
      internal::apply_block_householder_on_the_right(sub_dst, sub_vecs, m_coeffs.segment(k, bs), !m_reverse);
    }
  }
 
 public:
  template <typename Dest>
  inline void applyThisOnTheLeft(Dest& dst, bool inputIsIdentity = false) const {
    Matrix<Scalar, 1, Dest::ColsAtCompileTime, RowMajor, 1, Dest::MaxColsAtCompileTime> workspace;
    applyThisOnTheLeft(dst, workspace, inputIsIdentity);
  }

  template <typename Dest, typename Workspace>
  inline void applyThisOnTheLeft(Dest& dst, Workspace& workspace, bool inputIsIdentity = false) const {
    if (inputIsIdentity && m_reverse) inputIsIdentity = false;
    // if the entries are large enough, then apply the reflectors by block
    if (m_length >= BlockSize && dst.cols() > 1) {
      // Make sure we have at least 2 useful blocks, otherwise it is point-less:
      Index blockSize = m_length < Index(2 * BlockSize) ? (m_length + 1) / 2 : Index(BlockSize);
      for (Index i = 0; i < m_length; i += blockSize) {
        Index end = m_reverse ? (std::min)(m_length, i + blockSize) : m_length - i;
        Index k = m_reverse ? i : (std::max)(Index(0), end - blockSize);
        Index bs = end - k;
        Index start = k + m_shift;
 
        using SubVectorsType = Block<internal::remove_all_t<VectorsType>, Dynamic, Dynamic>;
        SubVectorsType sub_vecs1(m_vectors.const_cast_derived(), Side == OnTheRight ? k : start,
                                 Side == OnTheRight ? start : k, Side == OnTheRight ? bs : m_vectors.rows() - start,
                                 Side == OnTheRight ? m_vectors.cols() - start : bs);
        std::conditional_t<Side == OnTheRight, Transpose<SubVectorsType>, SubVectorsType&> sub_vecs(sub_vecs1);
 
        Index dstRows = rows() - m_shift - k;
 
        if (inputIsIdentity) {
          Block<Dest, Dynamic, Dynamic> sub_dst = dst.bottomRightCorner(dstRows, dstRows);
          apply_block_householder_on_the_left(sub_dst, sub_vecs, m_coeffs.segment(k, bs), !m_reverse);
        } else {
          auto sub_dst = dst.bottomRows(dstRows);
          apply_block_householder_on_the_left(sub_dst, sub_vecs, m_coeffs.segment(k, bs), !m_reverse);
        }
      }
    } else {
      workspace.resize(dst.cols());
      for (Index k = 0; k < m_length; ++k) {
        Index actual_k = m_reverse ? k : m_length - k - 1;
        Index dstRows = rows() - m_shift - actual_k;
 
        if (inputIsIdentity) {
          Block<Dest, Dynamic, Dynamic> sub_dst = dst.bottomRightCorner(dstRows, dstRows);
          sub_dst.applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
        } else {
          auto sub_dst = dst.bottomRows(dstRows);
          sub_dst.applyHouseholderOnTheLeft(essentialVector(actual_k), m_coeffs.coeff(actual_k), workspace.data());
        }
      }
    }
  }

  template <typename OtherDerived>
  typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(
      const MatrixBase<OtherDerived>& other) const {
    typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type res(
        other.template cast<typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::ResultScalar>());
    applyThisOnTheLeft(res, internal::is_identity<OtherDerived>::value && res.rows() == res.cols());
    return res;
  }

  template <typename OtherDerived>
  typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(
      const DiagonalBase<OtherDerived>& other) const {
    return applyToUpperTriangularOnTheLeft<
        typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type>(other, true);
  }

  template <typename OtherDerived>
  typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type operator*(
      const TriangularBase<OtherDerived>& other) const {
    constexpr bool kUpperTriangular =
        (int(OtherDerived::Mode) & (int(Upper) | int(Lower) | int(SelfAdjoint))) == int(Upper);
    return applyToUpperTriangularOnTheLeft<
        typename internal::matrix_type_times_scalar_type<Scalar, OtherDerived>::Type>(other, kUpperTriangular);
  }

  template <typename OtherDerived>
  typename internal::permutation_type_as_matrix_type<Scalar, OtherDerived>::Type operator*(
      const PermutationBase<OtherDerived>& other) const {
    using ResultType = typename internal::permutation_type_as_matrix_type<Scalar, OtherDerived>::Type;
    eigen_assert(cols() == other.rows());
    ResultType res = ResultType::Identity(rows(), rows());
    applyThisOnTheLeft(res, true);
    // The permutation product runs in place when its operand is the destination.
    res.noalias() = res * other.derived();
    return res;
  }

  template <typename OtherDerived>
  typename internal::permutation_type_as_matrix_type<Scalar, Inverse<OtherDerived>>::Type operator*(
      const InverseImpl<OtherDerived, PermutationStorage>& other) const {
    using ResultType = typename internal::permutation_type_as_matrix_type<Scalar, Inverse<OtherDerived>>::Type;
    eigen_assert(cols() == other.rows());
    ResultType res = ResultType::Identity(rows(), rows());
    applyThisOnTheLeft(res, true);
    res.noalias() = res * other.derived();
    return res;
  }
 
  template <typename VectorsType_, typename CoeffsType_, int Side_>
  friend struct internal::hseq_side_dependent_impl;
 
 private:
  // The identity path of applyThisOnTheLeft() skips the columns left of the rows each reflection acts on, which is
  // valid for any square input that vanishes below the diagonal: a diagonal or upper triangular matrix as much as
  // the identity. Other inputs take the full path.
  template <typename ResultType, typename OtherDerived>
  ResultType applyToUpperTriangularOnTheLeft(const EigenBase<OtherDerived>& other, bool upperTriangular) const {
    eigen_assert(cols() == other.rows());
    ResultType res;
    internal::assign_householder_operand(res, other);
    applyThisOnTheLeft(res, upperTriangular && res.rows() == res.cols());
    return res;
  }
 
 public:
  EIGEN_DEVICE_FUNC HouseholderSequence& setLength(Index length) {
    m_length = length;
    return *this;
  }

  EIGEN_DEVICE_FUNC HouseholderSequence& setShift(Index shift) {
    m_shift = shift;
    return *this;
  }
 
  EIGEN_DEVICE_FUNC Index length() const {
    return m_length;
  } 
 
  EIGEN_DEVICE_FUNC Index shift() const {
    return m_shift;
  } 
 
  /* Necessary for .adjoint() and .conjugate() */
  template <typename VectorsType2, typename CoeffsType2, int Side2>
  friend class HouseholderSequence;
 
 protected:
  HouseholderSequence& setReverseFlag(bool reverse) {
    m_reverse = reverse;
    return *this;
  }
 
  bool reverseFlag() const { return m_reverse; } 
 
  typename VectorsType::Nested m_vectors;
  typename CoeffsType::Nested m_coeffs;
  bool m_reverse;
  Index m_length;
  Index m_shift;
  enum { BlockSize = 48 };
};

template <typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar, OtherDerived>::Type operator*(
    const MatrixBase<OtherDerived>& other, const HouseholderSequence<VectorsType, CoeffsType, Side>& h) {
  typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar, OtherDerived>::Type res(
      other.template cast<typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar,
                                                                           OtherDerived>::ResultScalar>());
  h.applyThisOnTheRight(res);
  return res;
}

template <typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar, OtherDerived>::Type operator*(
    const DiagonalBase<OtherDerived>& other, const HouseholderSequence<VectorsType, CoeffsType, Side>& h) {
  using ResultType = typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar, OtherDerived>::Type;
  eigen_assert(other.cols() == h.rows());
  ResultType res = ResultType::Identity(h.rows(), h.rows());
  h.applyThisOnTheLeft(res, true);
  // The lazy diagonal product multiplies from the left, D(i,i) * H(i,j), as a noncommutative scalar requires, and
  // reads each coefficient only to overwrite it, so it runs in place.
  res = other.derived() * res;
  return res;
}

template <typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar, OtherDerived>::Type operator*(
    const TriangularBase<OtherDerived>& other, const HouseholderSequence<VectorsType, CoeffsType, Side>& h) {
  using ResultType = typename internal::matrix_type_times_scalar_type<typename VectorsType::Scalar, OtherDerived>::Type;
  ResultType res;
  internal::assign_householder_operand(res, other);
  h.applyThisOnTheRight(res);
  return res;
}

template <typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
typename internal::permutation_type_as_matrix_type<typename VectorsType::Scalar, OtherDerived>::Type operator*(
    const PermutationBase<OtherDerived>& other, const HouseholderSequence<VectorsType, CoeffsType, Side>& h) {
  using ResultType =
      typename internal::permutation_type_as_matrix_type<typename VectorsType::Scalar, OtherDerived>::Type;
  eigen_assert(other.cols() == h.rows());
  ResultType res = ResultType::Identity(h.rows(), h.rows());
  h.applyThisOnTheLeft(res, true);
  // The permutation product runs in place when its operand is the destination.
  res.noalias() = other.derived() * res;
  return res;
}

template <typename OtherDerived, typename VectorsType, typename CoeffsType, int Side>
typename internal::permutation_type_as_matrix_type<typename VectorsType::Scalar, Inverse<OtherDerived>>::Type operator*(
    const InverseImpl<OtherDerived, PermutationStorage>& other,
    const HouseholderSequence<VectorsType, CoeffsType, Side>& h) {
  using ResultType =
      typename internal::permutation_type_as_matrix_type<typename VectorsType::Scalar, Inverse<OtherDerived>>::Type;
  eigen_assert(other.cols() == h.rows());
  ResultType res = ResultType::Identity(h.rows(), h.rows());
  h.applyThisOnTheLeft(res, true);
  res.noalias() = other.derived() * res;
  return res;
}

template <typename VectorsType, typename CoeffsType>
HouseholderSequence<VectorsType, CoeffsType> householderSequence(const VectorsType& v, const CoeffsType& h) {
  return HouseholderSequence<VectorsType, CoeffsType, OnTheLeft>(v, h);
}

template <typename VectorsType, typename CoeffsType>
HouseholderSequence<VectorsType, CoeffsType, OnTheRight> rightHouseholderSequence(const VectorsType& v,
                                                                                  const CoeffsType& h) {
  return HouseholderSequence<VectorsType, CoeffsType, OnTheRight>(v, h);
}
 
}  // end namespace Eigen
 
#endif  // EIGEN_HOUSEHOLDER_SEQUENCE_H