VectorwiseOp.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2019 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2006-2008 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_PARTIAL_REDUX_H
#define EIGEN_PARTIAL_REDUX_H
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {

 
template <typename MatrixType, typename MemberOp, int Direction>
class PartialReduxExpr;
 
namespace internal {
 
template <typename MatrixType, typename MemberOp, int Direction>
struct traits<PartialReduxExpr<MatrixType, MemberOp, Direction>> : traits<MatrixType> {
  using Scalar = typename MemberOp::result_type;
  using StorageKind = typename traits<MatrixType>::StorageKind;
  using XprKind = typename traits<MatrixType>::XprKind;
  using InputScalar = typename MatrixType::Scalar;
  enum {
    RowsAtCompileTime = Direction == Vertical ? 1 : MatrixType::RowsAtCompileTime,
    ColsAtCompileTime = Direction == Horizontal ? 1 : MatrixType::ColsAtCompileTime,
    MaxRowsAtCompileTime = Direction == Vertical ? 1 : MatrixType::MaxRowsAtCompileTime,
    MaxColsAtCompileTime = Direction == Horizontal ? 1 : MatrixType::MaxColsAtCompileTime,
    Flags = RowsAtCompileTime == 1 ? RowMajorBit : 0,
    TraversalSize = Direction == Vertical ? MatrixType::RowsAtCompileTime : MatrixType::ColsAtCompileTime
  };
};
}  // namespace internal
 
template <typename MatrixType, typename MemberOp, int Direction>
class PartialReduxExpr : public internal::dense_xpr_base<PartialReduxExpr<MatrixType, MemberOp, Direction>>::type,
                         internal::no_assignment_operator {
 public:
  using Base = typename internal::dense_xpr_base<PartialReduxExpr>::type;
  EIGEN_DENSE_PUBLIC_INTERFACE(PartialReduxExpr)
 
  EIGEN_DEVICE_FUNC explicit PartialReduxExpr(const MatrixType& mat, const MemberOp& func = MemberOp())
      : m_matrix(mat), m_functor(func) {}
 
  EIGEN_DEVICE_FUNC constexpr Index rows() const noexcept { return (Direction == Vertical ? 1 : m_matrix.rows()); }
  EIGEN_DEVICE_FUNC constexpr Index cols() const noexcept { return (Direction == Horizontal ? 1 : m_matrix.cols()); }
 
  EIGEN_DEVICE_FUNC typename MatrixType::Nested nestedExpression() const { return m_matrix; }
 
  EIGEN_DEVICE_FUNC const MemberOp& functor() const { return m_functor; }
 
 protected:
  typename MatrixType::Nested m_matrix;
  const MemberOp m_functor;
};
 
template <typename A, typename B>
struct partial_redux_dummy_func;
 
#define EIGEN_MAKE_PARTIAL_REDUX_FUNCTOR(MEMBER, COST, VECTORIZABLE, BINARYOP)              \
  template <typename ResultType, typename Scalar>                                           \
  struct member_##MEMBER {                                                                  \
    typedef ResultType result_type;                                                         \
    typedef BINARYOP<Scalar, Scalar> BinaryOp;                                              \
    template <int Size>                                                                     \
    struct Cost : std::integral_constant<int, COST> {};                                     \
    enum { Vectorizable = VECTORIZABLE };                                                   \
    template <typename XprType>                                                             \
    EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ResultType operator()(const XprType& mat) const { \
      return mat.MEMBER();                                                                  \
    }                                                                                       \
    BinaryOp binaryFunc() const { return BinaryOp(); }                                      \
  }
 
#define EIGEN_MEMBER_FUNCTOR(MEMBER, COST) EIGEN_MAKE_PARTIAL_REDUX_FUNCTOR(MEMBER, COST, 0, partial_redux_dummy_func)
 
namespace internal {
 
EIGEN_MEMBER_FUNCTOR(norm, (Size + 5) * NumTraits<Scalar>::MulCost + (Size - 1) * NumTraits<Scalar>::AddCost);
// These multi-pass reductions must not inherit the cheaper one-pass norm cost,
// which could suppress materialization of an enclosing expression.
EIGEN_MEMBER_FUNCTOR(stableNorm, Size* NumTraits<Scalar>::ReadCost + (2 * Size + 5) * NumTraits<Scalar>::MulCost +
                                     (2 * Size - 1) * NumTraits<Scalar>::AddCost);
EIGEN_MEMBER_FUNCTOR(blueNorm, Size* NumTraits<Scalar>::ReadCost + (2 * Size + 8) * NumTraits<Scalar>::MulCost +
                                   (3 * Size - 1) * NumTraits<Scalar>::AddCost);
EIGEN_MEMBER_FUNCTOR(
    hypotNorm,
    ((NumTraits<Scalar>::IsComplex ? 2 * Size : Size) - 1) *
        functor_traits<
            scalar_hypot_op<typename stable_norm_accumulator<typename NumTraits<Scalar>::Real>::type>>::Cost);
EIGEN_MEMBER_FUNCTOR(all, (Size - 1) * NumTraits<Scalar>::AddCost);
EIGEN_MEMBER_FUNCTOR(any, (Size - 1) * NumTraits<Scalar>::AddCost);
EIGEN_MEMBER_FUNCTOR(count, (Size - 1) * NumTraits<Scalar>::AddCost);
 
EIGEN_MAKE_PARTIAL_REDUX_FUNCTOR(sum, (Size - 1) * NumTraits<Scalar>::AddCost, 1, internal::scalar_sum_op);
EIGEN_MAKE_PARTIAL_REDUX_FUNCTOR(minCoeff, (Size - 1) * NumTraits<Scalar>::AddCost, 1, internal::scalar_min_op);
EIGEN_MAKE_PARTIAL_REDUX_FUNCTOR(maxCoeff, (Size - 1) * NumTraits<Scalar>::AddCost, 1, internal::scalar_max_op);
EIGEN_MAKE_PARTIAL_REDUX_FUNCTOR(prod, (Size - 1) * NumTraits<Scalar>::MulCost, 1, internal::scalar_product_op);
 
template <typename ResultType, typename Scalar>
struct member_squaredNorm {
  using result_type = ResultType;
  static constexpr bool Vectorizable = false;
  template <int Size>
  struct Cost : std::integral_constant<int, Size * functor_traits<scalar_abs2_op<Scalar>>::Cost +
                                                (Size - 1) * NumTraits<ResultType>::AddCost> {};
  template <typename XprType>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE ResultType operator()(const XprType& mat) const {
    return mat.matrix().squaredNorm();
  }
};
 
template <typename ExpressionType, int Direction,
          bool ComplexInnerReduction =
              NumTraits<typename ExpressionType::Scalar>::IsComplex &&
              complex_array_access<typename ExpressionType::Scalar>::value &&
              (Direction == (ExpressionType::IsRowMajor ? Horizontal : Vertical)) &&
              ((int(evaluator<ExpressionType>::Flags) & (DirectAccessBit | PacketAccessBit)) != 0)>
struct partial_squared_norm {
  using Scalar = typename ExpressionType::Scalar;
  using RealScalar = typename ExpressionType::RealScalar;
  using Type = PartialReduxExpr<const CwiseUnaryOp<scalar_abs2_op<Scalar>, const ExpressionType>,
                                member_sum<RealScalar, RealScalar>, Direction>;
  EIGEN_DEVICE_FUNC static Type run(const ExpressionType& matrix) { return Type(matrix.cwiseAbs2()); }
};
 
// Keep scalars without array-oriented access and scalar-only expressions on abs2(): realView()
// copies/evaluates their coefficients twice. Packet expressions still evaluate scalar tails
// and reductions shorter than one packet twice.
template <typename ExpressionType, int Direction>
struct partial_squared_norm<ExpressionType, Direction, true> {
  using Type =
      PartialReduxExpr<const ExpressionType,
                       member_squaredNorm<typename ExpressionType::RealScalar, typename ExpressionType::Scalar>,
                       Direction>;
  EIGEN_DEVICE_FUNC static Type run(const ExpressionType& matrix) { return Type(matrix); }
};
 
template <int p, typename ResultType, typename Scalar>
struct member_lpnorm {
  using result_type = ResultType;
  enum { Vectorizable = 0 };
  template <int Size>
  struct Cost
      : std::integral_constant<int, (Size + 5) * NumTraits<Scalar>::MulCost + (Size - 1) * NumTraits<Scalar>::AddCost> {
  };
  EIGEN_DEVICE_FUNC member_lpnorm() = default;
  template <typename XprType>
  EIGEN_DEVICE_FUNC inline ResultType operator()(const XprType& mat) const {
    return mat.template lpNorm<p>();
  }
};
 
template <typename BinaryOpT, typename Scalar>
struct member_redux {
  using BinaryOp = BinaryOpT;
  using result_type = typename result_of<BinaryOp(const Scalar&, const Scalar&)>::type;
 
  enum { Vectorizable = functor_traits<BinaryOp>::PacketAccess };
  template <int Size>
  struct Cost : std::integral_constant<int, (Size - 1) * functor_traits<BinaryOp>::Cost> {};
  EIGEN_DEVICE_FUNC explicit member_redux(const BinaryOp func) : m_functor(func) {}
  template <typename Derived>
  EIGEN_DEVICE_FUNC inline result_type operator()(const DenseBase<Derived>& mat) const {
    return mat.redux(m_functor);
  }
  const BinaryOp& binaryFunc() const { return m_functor; }
  const BinaryOp m_functor;
};
 
template <typename Scalar>
struct scalar_replace_zero_with_one_op {
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Scalar operator()(const Scalar& x) const {
    return numext::is_exactly_zero(x) ? Scalar(1) : x;
  }
  template <typename Packet>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE Packet packetOp(const Packet& x) const {
    return pselect(pcmp_eq(x, pzero(x)), pset1<Packet>(Scalar(1)), x);
  }
};
template <typename Scalar>
struct functor_traits<scalar_replace_zero_with_one_op<Scalar>> {
  enum { Cost = 1, PacketAccess = packet_traits<Scalar>::HasCmp };
};
 
}  // namespace internal

template <typename ExpressionType, int Direction>
class VectorwiseOp {
 public:
  using Scalar = typename ExpressionType::Scalar;
  using RealScalar = typename ExpressionType::RealScalar;
  using Index = Eigen::Index;  
  using ExpressionTypeNested = typename internal::ref_selector<ExpressionType>::non_const_type;
  using ExpressionTypeNestedCleaned = internal::remove_all_t<ExpressionTypeNested>;
 
  template <template <typename OutScalar, typename InputScalar> class Functor, typename ReturnScalar = Scalar>
  struct ReturnType {
    using Type = PartialReduxExpr<ExpressionType, Functor<ReturnScalar, Scalar>, Direction>;
  };
 
  template <typename BinaryOp>
  struct ReduxReturnType {
    using Type = PartialReduxExpr<ExpressionType, internal::member_redux<BinaryOp, Scalar>, Direction>;
  };
 
  enum { isVertical = (Direction == Vertical) ? 1 : 0, isHorizontal = (Direction == Horizontal) ? 1 : 0 };
 
 protected:
  template <typename OtherDerived>
  struct ExtendedType {
    using Type = Replicate<OtherDerived, isVertical ? 1 : ExpressionType::RowsAtCompileTime,
                           isHorizontal ? 1 : ExpressionType::ColsAtCompileTime>;
  };

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC typename ExtendedType<OtherDerived>::Type extendedTo(const DenseBase<OtherDerived>& other) const {
    EIGEN_STATIC_ASSERT(internal::check_implication(isVertical, OtherDerived::MaxColsAtCompileTime == 1),
                        YOU_PASSED_A_ROW_VECTOR_BUT_A_COLUMN_VECTOR_WAS_EXPECTED)
    EIGEN_STATIC_ASSERT(internal::check_implication(isHorizontal, OtherDerived::MaxRowsAtCompileTime == 1),
                        YOU_PASSED_A_COLUMN_VECTOR_BUT_A_ROW_VECTOR_WAS_EXPECTED)
    return typename ExtendedType<OtherDerived>::Type(other.derived(), isVertical ? 1 : m_matrix.rows(),
                                                     isHorizontal ? 1 : m_matrix.cols());
  }
 
  template <typename OtherDerived>
  struct OppositeExtendedType {
    using Type = Replicate<OtherDerived, isHorizontal ? 1 : ExpressionType::RowsAtCompileTime,
                           isVertical ? 1 : ExpressionType::ColsAtCompileTime>;
  };

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC typename OppositeExtendedType<OtherDerived>::Type extendedToOpposite(
      const DenseBase<OtherDerived>& other) const {
    EIGEN_STATIC_ASSERT(internal::check_implication(isHorizontal, OtherDerived::MaxColsAtCompileTime == 1),
                        YOU_PASSED_A_ROW_VECTOR_BUT_A_COLUMN_VECTOR_WAS_EXPECTED)
    EIGEN_STATIC_ASSERT(internal::check_implication(isVertical, OtherDerived::MaxRowsAtCompileTime == 1),
                        YOU_PASSED_A_COLUMN_VECTOR_BUT_A_ROW_VECTOR_WAS_EXPECTED)
    return typename OppositeExtendedType<OtherDerived>::Type(other.derived(), isHorizontal ? 1 : m_matrix.rows(),
                                                             isVertical ? 1 : m_matrix.cols());
  }
 
 public:
  EIGEN_DEVICE_FUNC explicit inline VectorwiseOp(ExpressionType& matrix) : m_matrix(matrix) {}

  EIGEN_DEVICE_FUNC inline const ExpressionType& _expression() const { return m_matrix; }
 
#ifdef EIGEN_PARSED_BY_DOXYGEN
  random_access_iterator_type iterator;
  random_access_iterator_type const_iterator;
#else
  using iterator = internal::subvector_stl_iterator<ExpressionType, DirectionType(Direction)>;
  using const_iterator = internal::subvector_stl_iterator<const ExpressionType, DirectionType(Direction)>;
  using reverse_iterator = internal::subvector_stl_reverse_iterator<ExpressionType, DirectionType(Direction)>;
  using const_reverse_iterator =
      internal::subvector_stl_reverse_iterator<const ExpressionType, DirectionType(Direction)>;
#endif

  iterator begin() { return iterator(m_matrix, 0); }
  const_iterator begin() const { return const_iterator(m_matrix, 0); }
  const_iterator cbegin() const { return const_iterator(m_matrix, 0); }

  reverse_iterator rbegin() {
    return reverse_iterator(m_matrix, m_matrix.template subVectors<DirectionType(Direction)>() - 1);
  }

  const_reverse_iterator rbegin() const {
    return const_reverse_iterator(m_matrix, m_matrix.template subVectors<DirectionType(Direction)>() - 1);
  }

  const_reverse_iterator crbegin() const {
    return const_reverse_iterator(m_matrix, m_matrix.template subVectors<DirectionType(Direction)>() - 1);
  }

  iterator end() { return iterator(m_matrix, m_matrix.template subVectors<DirectionType(Direction)>()); }
  const_iterator end() const {
    return const_iterator(m_matrix, m_matrix.template subVectors<DirectionType(Direction)>());
  }

  const_iterator cend() const {
    return const_iterator(m_matrix, m_matrix.template subVectors<DirectionType(Direction)>());
  }

  reverse_iterator rend() { return reverse_iterator(m_matrix, -1); }
  const_reverse_iterator rend() const { return const_reverse_iterator(m_matrix, -1); }
  const_reverse_iterator crend() const { return const_reverse_iterator(m_matrix, -1); }

  template <typename BinaryOp>
  EIGEN_DEVICE_FUNC const typename ReduxReturnType<BinaryOp>::Type redux(const BinaryOp& func = BinaryOp()) const {
    eigen_assert(redux_length() > 0 && "you are using an empty matrix");
    return typename ReduxReturnType<BinaryOp>::Type(_expression(), internal::member_redux<BinaryOp, Scalar>(func));
  }
 
  using MinCoeffReturnType = typename ReturnType<Eigen::internal::member_minCoeff>::Type;
  using MaxCoeffReturnType = typename ReturnType<Eigen::internal::member_maxCoeff>::Type;
  using SquaredNormReturnType = typename internal::partial_squared_norm<ExpressionTypeNestedCleaned, Direction>::Type;
  using NormReturnType = CwiseUnaryOp<internal::scalar_sqrt_op<RealScalar>, const SquaredNormReturnType>;
  using BlueNormReturnType = typename ReturnType<Eigen::internal::member_blueNorm, RealScalar>::Type;
  using StableNormReturnType = typename ReturnType<Eigen::internal::member_stableNorm, RealScalar>::Type;
  using HypotNormReturnType = typename ReturnType<Eigen::internal::member_hypotNorm, RealScalar>::Type;
  using SumReturnType = typename ReturnType<Eigen::internal::member_sum>::Type;
  using MeanReturnType = EIGEN_EXPR_BINARYOP_SCALAR_RETURN_TYPE(SumReturnType, Scalar, internal::scalar_quotient_op);
  using AllReturnType = typename ReturnType<Eigen::internal::member_all, bool>::Type;
  using AnyReturnType = typename ReturnType<Eigen::internal::member_any, bool>::Type;
  using CountReturnType = PartialReduxExpr<ExpressionType, internal::member_count<Index, Scalar>, Direction>;
  using ProdReturnType = typename ReturnType<Eigen::internal::member_prod>::Type;
  using ConstReverseReturnType = Reverse<const ExpressionType, Direction>;
  using ReverseReturnType = Reverse<ExpressionType, Direction>;
 
  template <int p>
  struct LpNormReturnType {
    using Type = PartialReduxExpr<ExpressionType, internal::member_lpnorm<p, RealScalar, Scalar>, Direction>;
  };

  EIGEN_DEVICE_FUNC const MinCoeffReturnType minCoeff() const {
    eigen_assert(redux_length() > 0 && "you are using an empty matrix");
    return MinCoeffReturnType(_expression());
  }

  EIGEN_DEVICE_FUNC const MaxCoeffReturnType maxCoeff() const {
    eigen_assert(redux_length() > 0 && "you are using an empty matrix");
    return MaxCoeffReturnType(_expression());
  }

  EIGEN_DEVICE_FUNC const SquaredNormReturnType squaredNorm() const {
    return internal::partial_squared_norm<ExpressionTypeNestedCleaned, Direction>::run(m_matrix);
  }

  EIGEN_DEVICE_FUNC const NormReturnType norm() const { return NormReturnType(squaredNorm()); }

   *
   * \sa DenseBase::norm() */
  template <int p>
  EIGEN_DEVICE_FUNC const typename LpNormReturnType<p>::Type lpNorm() const {
    return typename LpNormReturnType<p>::Type(_expression());
  }

   * This is a vector with real entries, even if the original matrix has complex entries.
   *
   * \sa DenseBase::blueNorm() */
  EIGEN_DEVICE_FUNC const BlueNormReturnType blueNorm() const { return BlueNormReturnType(_expression()); }

  EIGEN_DEVICE_FUNC const StableNormReturnType stableNorm() const { return StableNormReturnType(_expression()); }

  EIGEN_DEVICE_FUNC const HypotNormReturnType hypotNorm() const { return HypotNormReturnType(_expression()); }

  EIGEN_DEVICE_FUNC const SumReturnType sum() const { return SumReturnType(_expression()); }

  EIGEN_DEVICE_FUNC const MeanReturnType mean() const {
    return sum() / Scalar(Direction == Vertical ? m_matrix.rows() : m_matrix.cols());
  }

  EIGEN_DEVICE_FUNC const AllReturnType all() const { return AllReturnType(_expression()); }

  EIGEN_DEVICE_FUNC const AnyReturnType any() const { return AnyReturnType(_expression()); }

  EIGEN_DEVICE_FUNC const CountReturnType count() const { return CountReturnType(_expression()); }

  EIGEN_DEVICE_FUNC const ProdReturnType prod() const { return ProdReturnType(_expression()); }

  EIGEN_DEVICE_FUNC const ConstReverseReturnType reverse() const { return ConstReverseReturnType(_expression()); }

  EIGEN_DEVICE_FUNC ReverseReturnType reverse() { return ReverseReturnType(_expression()); }
 
  using ReplicateReturnType = Replicate<ExpressionType, (isVertical ? Dynamic : 1), (isHorizontal ? Dynamic : 1)>;
  EIGEN_DEVICE_FUNC const ReplicateReturnType replicate(Index factor) const;

  // NOTE implemented here because of sunstudio's compilation errors
  // isVertical*Factor+isHorizontal instead of (isVertical?Factor:1) to handle CUDA bug with ternary operator
  template <int Factor>
  const Replicate<ExpressionType, isVertical * Factor + isHorizontal,
                  isHorizontal * Factor + isVertical> EIGEN_DEVICE_FUNC
  replicate(Index factor = Factor) const {
    return Replicate<ExpressionType, (isVertical ? Factor : 1), (isHorizontal ? Factor : 1)>(
        _expression(), isVertical ? factor : 1, isHorizontal ? factor : 1);
  }

  template <typename NFactorType>
  EIGEN_DEVICE_FUNC const
      Replicate<ExpressionType, isVertical * internal::get_fixed_value<NFactorType>::value + isHorizontal,
                isHorizontal * internal::get_fixed_value<NFactorType>::value + isVertical>
      replicate(NFactorType nFactor) const {
    return replicate<internal::get_fixed_value<NFactorType>::value>(internal::get_runtime_value(nFactor));
  }

 
  // The broadcast (compound-)assignments below bind the rhs through
  // `nested_eval<.., Dynamic>` so expressions like `colwise().sum()` are
  // materialized once; otherwise an aliased reduction would be recomputed
  // against partially-updated coefficients (issue #1731).

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC ExpressionType& operator=(const DenseBase<OtherDerived>& other) {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
    EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived)
    // eigen_assert((m_matrix.isNull()) == (other.isNull())); FIXME
    typename internal::nested_eval<OtherDerived, Dynamic>::type other_eval(other.derived());
    return m_matrix = extendedTo(other_eval);
  }

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC ExpressionType& operator+=(const DenseBase<OtherDerived>& other) {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
    EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived)
    typename internal::nested_eval<OtherDerived, Dynamic>::type other_eval(other.derived());
    return m_matrix += extendedTo(other_eval);
  }

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC ExpressionType& operator-=(const DenseBase<OtherDerived>& other) {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
    EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived)
    typename internal::nested_eval<OtherDerived, Dynamic>::type other_eval(other.derived());
    return m_matrix -= extendedTo(other_eval);
  }

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC ExpressionType& operator*=(const DenseBase<OtherDerived>& other) {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
    EIGEN_STATIC_ASSERT_ARRAYXPR(ExpressionType)
    EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived)
    typename internal::nested_eval<OtherDerived, Dynamic>::type other_eval(other.derived());
    m_matrix *= extendedTo(other_eval);
    return m_matrix;
  }

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC ExpressionType& operator/=(const DenseBase<OtherDerived>& other) {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
    EIGEN_STATIC_ASSERT_ARRAYXPR(ExpressionType)
    EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived)
    typename internal::nested_eval<OtherDerived, Dynamic>::type other_eval(other.derived());
    m_matrix /= extendedTo(other_eval);
    return m_matrix;
  }

  template <typename OtherDerived>
  EIGEN_STRONG_INLINE EIGEN_DEVICE_FUNC
      CwiseBinaryOp<internal::scalar_sum_op<Scalar, typename OtherDerived::Scalar>, const ExpressionTypeNestedCleaned,
                    const typename ExtendedType<OtherDerived>::Type>
      operator+(const DenseBase<OtherDerived>& other) const {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
    EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived)
    return m_matrix + extendedTo(other.derived());
  }

  template <typename OtherDerived>
  friend EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
      CwiseBinaryOp<internal::scalar_sum_op<typename OtherDerived::Scalar, Scalar>,
                    const typename ExtendedType<OtherDerived>::Type, const ExpressionTypeNestedCleaned>
      operator+(const DenseBase<OtherDerived>& other, const VectorwiseOp& xpr) {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
    EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived)
    return CwiseBinaryOp<internal::scalar_sum_op<typename OtherDerived::Scalar, Scalar>,
                         const typename ExtendedType<OtherDerived>::Type, const ExpressionTypeNestedCleaned>(
        xpr.extendedTo(other.derived()), xpr._expression());
  }

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC CwiseBinaryOp<internal::scalar_difference_op<Scalar, typename OtherDerived::Scalar>,
                                  const ExpressionTypeNestedCleaned, const typename ExtendedType<OtherDerived>::Type>
  operator-(const DenseBase<OtherDerived>& other) const {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
    EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived)
    return m_matrix - extendedTo(other.derived());
  }

  template <typename OtherDerived>
  friend EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
      CwiseBinaryOp<internal::scalar_difference_op<typename OtherDerived::Scalar, Scalar>,
                    const typename ExtendedType<OtherDerived>::Type, const ExpressionTypeNestedCleaned>
      operator-(const DenseBase<OtherDerived>& other, const VectorwiseOp& xpr) {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
    EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived)
    return CwiseBinaryOp<internal::scalar_difference_op<typename OtherDerived::Scalar, Scalar>,
                         const typename ExtendedType<OtherDerived>::Type, const ExpressionTypeNestedCleaned>(
        xpr.extendedTo(other.derived()), xpr._expression());
  }

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC CwiseBinaryOp<internal::scalar_product_op<Scalar, typename OtherDerived::Scalar>,
                                  const ExpressionTypeNestedCleaned, const typename ExtendedType<OtherDerived>::Type>
  operator*(const DenseBase<OtherDerived>& other) const {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
    EIGEN_STATIC_ASSERT_ARRAYXPR(ExpressionType)
    EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived)
    return m_matrix * extendedTo(other.derived());
  }

  template <typename OtherDerived>
  friend EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
      CwiseBinaryOp<internal::scalar_product_op<typename OtherDerived::Scalar, Scalar>,
                    const typename ExtendedType<OtherDerived>::Type, const ExpressionTypeNestedCleaned>
      operator*(const DenseBase<OtherDerived>& other, const VectorwiseOp& xpr) {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
    EIGEN_STATIC_ASSERT_ARRAYXPR(ExpressionType)
    EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived)
    return CwiseBinaryOp<internal::scalar_product_op<typename OtherDerived::Scalar, Scalar>,
                         const typename ExtendedType<OtherDerived>::Type, const ExpressionTypeNestedCleaned>(
        xpr.extendedTo(other.derived()), xpr._expression());
  }

  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC CwiseBinaryOp<internal::scalar_quotient_op<Scalar, typename OtherDerived::Scalar>,
                                  const ExpressionTypeNestedCleaned, const typename ExtendedType<OtherDerived>::Type>
  operator/(const DenseBase<OtherDerived>& other) const {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
    EIGEN_STATIC_ASSERT_ARRAYXPR(ExpressionType)
    EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived)
    return m_matrix / extendedTo(other.derived());
  }

  template <typename OtherDerived>
  friend EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE
      CwiseBinaryOp<internal::scalar_quotient_op<typename OtherDerived::Scalar, Scalar>,
                    const typename ExtendedType<OtherDerived>::Type, const ExpressionTypeNestedCleaned>
      operator/(const DenseBase<OtherDerived>& other, const VectorwiseOp& xpr) {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
    EIGEN_STATIC_ASSERT_ARRAYXPR(ExpressionType)
    EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived)
    return CwiseBinaryOp<internal::scalar_quotient_op<typename OtherDerived::Scalar, Scalar>,
                         const typename ExtendedType<OtherDerived>::Type, const ExpressionTypeNestedCleaned>(
        xpr.extendedTo(other.derived()), xpr._expression());
  }

  template <typename CustomBinaryOp, typename OtherDerived>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const
      CwiseBinaryOp<CustomBinaryOp, const ExpressionTypeNestedCleaned, const typename ExtendedType<OtherDerived>::Type>
      binaryExpr(const DenseBase<OtherDerived>& other, const CustomBinaryOp& func = CustomBinaryOp()) const {
    EIGEN_STATIC_ASSERT_VECTOR_ONLY(OtherDerived)
    EIGEN_STATIC_ASSERT_SAME_XPR_KIND(ExpressionType, OtherDerived)
    return CwiseBinaryOp<CustomBinaryOp, const ExpressionTypeNestedCleaned,
                         const typename ExtendedType<OtherDerived>::Type>(_expression(), extendedTo(other.derived()),
                                                                          func);
  }

  template <int NaNPropagation = PropagateFast, typename OtherDerived>
  EIGEN_DEVICE_FUNC constexpr EIGEN_STRONG_INLINE const
      CwiseBinaryOp<internal::scalar_min_op<Scalar, typename OtherDerived::Scalar, NaNPropagation>,
                    const ExpressionTypeNestedCleaned, const typename ExtendedType<OtherDerived>::Type>
      cwiseMin(const DenseBase<OtherDerived>& other) const {
    return binaryExpr(other, internal::scalar_min_op<Scalar, typename OtherDerived::Scalar, NaNPropagation>());
  }

  template <int NaNPropagation = PropagateFast>
  EIGEN_DEVICE_FUNC constexpr EIGEN_STRONG_INLINE const
      CwiseBinaryOp<internal::scalar_min_op<Scalar, Scalar, NaNPropagation>, const ExpressionTypeNestedCleaned,
                    const typename ExpressionTypeNestedCleaned::ConstantReturnType>
      cwiseMin(const Scalar& other) const {
    return CwiseBinaryOp<internal::scalar_min_op<Scalar, Scalar, NaNPropagation>, const ExpressionTypeNestedCleaned,
                         const typename ExpressionTypeNestedCleaned::ConstantReturnType>(
        _expression(), ExpressionTypeNestedCleaned::Constant(m_matrix.rows(), m_matrix.cols(), other));
  }

  template <int NaNPropagation = PropagateFast, typename OtherDerived>
  EIGEN_DEVICE_FUNC constexpr EIGEN_STRONG_INLINE const
      CwiseBinaryOp<internal::scalar_max_op<Scalar, typename OtherDerived::Scalar, NaNPropagation>,
                    const ExpressionTypeNestedCleaned, const typename ExtendedType<OtherDerived>::Type>
      cwiseMax(const DenseBase<OtherDerived>& other) const {
    return binaryExpr(other, internal::scalar_max_op<Scalar, typename OtherDerived::Scalar, NaNPropagation>());
  }

  template <int NaNPropagation = PropagateFast>
  EIGEN_DEVICE_FUNC constexpr EIGEN_STRONG_INLINE const
      CwiseBinaryOp<internal::scalar_max_op<Scalar, Scalar, NaNPropagation>, const ExpressionTypeNestedCleaned,
                    const typename ExpressionTypeNestedCleaned::ConstantReturnType>
      cwiseMax(const Scalar& other) const {
    return CwiseBinaryOp<internal::scalar_max_op<Scalar, Scalar, NaNPropagation>, const ExpressionTypeNestedCleaned,
                         const typename ExpressionTypeNestedCleaned::ConstantReturnType>(
        _expression(), ExpressionTypeNestedCleaned::Constant(m_matrix.rows(), m_matrix.cols(), other));
  }
 
  using Normalized_NonzeroNormType =
      CwiseUnaryOp<internal::scalar_replace_zero_with_one_op<Scalar>, const NormReturnType>;
  using NormalizedReturnType = CwiseBinaryOp<internal::scalar_quotient_op<Scalar>, const ExpressionTypeNestedCleaned,
                                             const typename OppositeExtendedType<Normalized_NonzeroNormType>::Type>;

  EIGEN_DEVICE_FUNC NormalizedReturnType normalized() const {
    return m_matrix.cwiseQuotient(extendedToOpposite(Normalized_NonzeroNormType(this->norm())));
  }

  EIGEN_DEVICE_FUNC void normalize() { m_matrix = this->normalized(); }
 
  EIGEN_DEVICE_FUNC inline void reverseInPlace();

 
  using HomogeneousReturnType = Homogeneous<ExpressionType, Direction>;
  EIGEN_DEVICE_FUNC HomogeneousReturnType homogeneous() const;
 
  using CrossReturnType = typename ExpressionType::PlainObject;
  template <typename OtherDerived>
  EIGEN_DEVICE_FUNC const CrossReturnType cross(const MatrixBase<OtherDerived>& other) const;
 
  enum {
    HNormalized_Size = Direction == Vertical ? internal::traits<ExpressionType>::RowsAtCompileTime
                                             : internal::traits<ExpressionType>::ColsAtCompileTime,
    HNormalized_SizeMinusOne = HNormalized_Size == Dynamic ? Dynamic : HNormalized_Size - 1
  };
  using HNormalized_Block = Block<const ExpressionType,
                                  Direction == Vertical ? int(HNormalized_SizeMinusOne)
                                                        : int(internal::traits<ExpressionType>::RowsAtCompileTime),
                                  Direction == Horizontal ? int(HNormalized_SizeMinusOne)
                                                          : int(internal::traits<ExpressionType>::ColsAtCompileTime)>;
  using HNormalized_Factors =
      Block<const ExpressionType, Direction == Vertical ? 1 : int(internal::traits<ExpressionType>::RowsAtCompileTime),
            Direction == Horizontal ? 1 : int(internal::traits<ExpressionType>::ColsAtCompileTime)>;
  using HNormalizedReturnType =
      CwiseBinaryOp<internal::scalar_quotient_op<typename internal::traits<ExpressionType>::Scalar>,
                    const HNormalized_Block,
                    const Replicate<HNormalized_Factors, Direction == Vertical ? HNormalized_SizeMinusOne : 1,
                                    Direction == Horizontal ? HNormalized_SizeMinusOne : 1>>;
 
  EIGEN_DEVICE_FUNC const HNormalizedReturnType hnormalized() const;
 
#ifdef EIGEN_VECTORWISEOP_PLUGIN
#include EIGEN_VECTORWISEOP_PLUGIN
#endif
 
 protected:
  EIGEN_DEVICE_FUNC Index redux_length() const { return Direction == Vertical ? m_matrix.rows() : m_matrix.cols(); }
  ExpressionTypeNested m_matrix;
};
 
// const colwise moved to DenseBase.h due to CUDA compiler bug

template <typename Derived>
EIGEN_DEVICE_FUNC inline typename DenseBase<Derived>::ColwiseReturnType DenseBase<Derived>::colwise() {
  return ColwiseReturnType(derived());
}
 
// const rowwise moved to DenseBase.h due to CUDA compiler bug

template <typename Derived>
EIGEN_DEVICE_FUNC inline typename DenseBase<Derived>::RowwiseReturnType DenseBase<Derived>::rowwise() {
  return RowwiseReturnType(derived());
}
 
}  // end namespace Eigen
 
#endif  // EIGEN_PARTIAL_REDUX_H