Eigen/src/Core/GeneralProduct.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2006-2008 Benoit Jacob <jacob.benoit.1@gmail.com>
// Copyright (C) 2008-2011 Gael Guennebaud <gael.guennebaud@inria.fr>
//
// 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_GENERAL_PRODUCT_H
#define EIGEN_GENERAL_PRODUCT_H

// IWYU pragma: private
#include "./InternalHeaderCheck.h"

namespace Eigen {

enum { Large = 2, Small = 3 };

// Runtime-size threshold for falling back from the generic matrix-matrix product
// implementation (heavy) to the lightweight coeff-based product one. See
// generic_product_impl<Lhs,Rhs,DenseShape,DenseShape,GemmProduct> in
// products/GeneralMatrixMatrix.h for more details.
// The crossover belongs to the kernel the GEMM path would select, not to the
// build: a scalar type that falls back to the generic gebp kernel keeps this
// value even in an SME build. The SME kernels have their own crossovers
// -- sme_gemm_to_coeffbased_threshold and
// EIGEN_SME_FIXED_SIZE_GEMM_TO_COEFFBASED_THRESHOLD below.
#ifndef EIGEN_GEMM_TO_COEFFBASED_THRESHOLD
// This default value has been obtained on a Haswell architecture.
#define EIGEN_GEMM_TO_COEFFBASED_THRESHOLD 20
#endif

// Fixed-size products can reach the GEMM product path even when the
// coeff-based evaluator is still faster. Keep this threshold separate so
// runtime-size dispatch preserves the historical heuristic above. This default
// was tuned on the same Haswell system as the runtime threshold. It does not
// follow a user-defined EIGEN_GEMM_TO_COEFFBASED_THRESHOLD, which still applies
// at run time, so a fixed-size product's crossover is the larger of the two.
#ifndef EIGEN_FIXED_SIZE_GEMM_TO_COEFFBASED_THRESHOLD
#define EIGEN_FIXED_SIZE_GEMM_TO_COEFFBASED_THRESHOLD 40
#endif

#ifdef EIGEN_VECTORIZE_SME
// Read only for SME-kernel scalar pairs; fitted on Apple M4 Pro to today's NEON small-block bounds (re-measure: #3194).
#ifndef EIGEN_SME_FIXED_SIZE_GEMM_TO_COEFFBASED_THRESHOLD
#define EIGEN_SME_FIXED_SIZE_GEMM_TO_COEFFBASED_THRESHOLD 61
#endif
#endif

// The dimension-sum bound in sme_gemm_to_coeffbased_threshold below grows with
// the depth, so it cannot reach a small output over a long depth: a 2x2 result
// at k=1024 sums to 1028 and takes the SME path with nearly the whole ZA grid
// predicated off. The crossover for that shape family is an output *area* rather
// than a dimension sum, and it tracks the scalar width.
//
// Measured on Apple M4 at SVL=512, every m x n with m,n in [2,12] against the
// coeff-based path: the largest area at which *no* shape loses is 27, 14, 15 and
// 6 for float, complex<float>, double and complex<double>. 96/sizeof (24, 12,
// 12, 6) sits just inside all four, and holds out to k=32768 -- the worst cell
// there is float 4x6 at 1.03x, i.e. neutral, while 2x2 is 13x and 3x3 2.0x.
// A flat constant does not work: the value that is safe for complex<double>
// admits only 2x2 for float, and float's would cost complex<double> up to 2x.
#ifndef EIGEN_SME_GEMM_TO_COEFFBASED_OUTPUT_AREA_THRESHOLD
#define EIGEN_SME_GEMM_TO_COEFFBASED_OUTPUT_AREA_THRESHOLD(Scalar) (96 / int(sizeof(Scalar)))
#endif

namespace internal {

#ifdef EIGEN_VECTORIZE_SME
// Defined in products/GeneralBlockPanelKernel.h, which Core includes after this
// header; only the declaration is needed to form the dependent type below.
template <typename LhsScalar, typename RhsScalar>
struct sme_has_gebp_kernel;
// Dimension sum below which the coeff-based product beats the SME GEMM kernel.
// mr x nr is sized in ZA tiles, so a product too small to fill the grid wastes
// most of the outer products the kernel issues, and the block holds fewer
// scalars as the scalar widens: the crossover sits well above the Haswell-tuned
// EIGEN_GEMM_TO_COEFFBASED_THRESHOLD, and it moves with the scalar type.
//
// First cube (m=n=k) at which the SME kernel wins, from two sweeps that each
// built the same source twice with only this threshold changed, so that the real
// dispatch selected the path in both:
//
//   scalar            M4   M4 Pro   M4 Pro, result 64-byte aligned
//   complex<double>   14   14       15
//   double            16   18       16
//   complex<float>    20   19       19
//   float             26   22       23
//
// The entries below are 3n for the largest of each row, so neither host loses at
// the crossover it gets. Over the sizes where the two disagree the paths are
// within 1.2x on both, while at n=14 a too-low bound costs 1.7x to 2.9x.
//
// Defining EIGEN_SME_GEMM_TO_COEFFBASED_THRESHOLD replaces every entry with one
// value, which is what such a sweep wants. The primary template answers with the
// generic threshold: sme_has_gebp_kernel gates the dispatch, so a scalar pair
// without an SME kernel never reads the table.
// With the NEON small-block path (arch/SME/GeneralBlockPanelKernel.h) the SME
// kernels take the generic crossover: below it the small products run the NEON
// packers and kernel, so the SME build matches a NEON build there.
#ifdef EIGEN_SME_GEMM_TO_COEFFBASED_THRESHOLD
template <typename Scalar>
struct sme_gemm_to_coeffbased_threshold : std::integral_constant<int, EIGEN_SME_GEMM_TO_COEFFBASED_THRESHOLD> {};
#else
template <typename Scalar>
struct sme_gemm_to_coeffbased_threshold : std::integral_constant<int, EIGEN_GEMM_TO_COEFFBASED_THRESHOLD> {};
#endif
#endif

template <int Rows, int Cols, int Depth>
struct product_type_selector;

template <int Size, int MaxSize>
struct product_size_category {
  enum {
#ifndef EIGEN_GPU_COMPILE_PHASE
    is_large = MaxSize == Dynamic || Size >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD ||
               (Size == Dynamic && MaxSize >= EIGEN_CACHEFRIENDLY_PRODUCT_THRESHOLD),
#else
    is_large = 0,
#endif
    value = is_large    ? Large
            : Size == 1 ? 1
                        : Small
  };
};

template <typename Lhs, typename Rhs>
struct product_type {
  using Lhs_ = remove_all_t<Lhs>;
  using Rhs_ = remove_all_t<Rhs>;
  enum {
    MaxRows = traits<Lhs_>::MaxRowsAtCompileTime,
    Rows = traits<Lhs_>::RowsAtCompileTime,
    MaxCols = traits<Rhs_>::MaxColsAtCompileTime,
    Cols = traits<Rhs_>::ColsAtCompileTime,
    MaxDepth = min_size_prefer_fixed(traits<Lhs_>::MaxColsAtCompileTime, traits<Rhs_>::MaxRowsAtCompileTime),
    Depth = min_size_prefer_fixed(traits<Lhs_>::ColsAtCompileTime, traits<Rhs_>::RowsAtCompileTime)
  };

  static constexpr int ProductType =
      product_type_selector<product_size_category<Rows, MaxRows>::value, product_size_category<Cols, MaxCols>::value,
                            product_size_category<Depth, MaxDepth>::value>::value;

  // Same rule as the runtime dispatch in products/GeneralMatrixMatrix.h: the
  // threshold follows the kernel the GEMM path would select, so a scalar pair
  // without an SME kernel keeps the generic crossover in an SME build.
  static constexpr int FixedSizeThreshold =
#ifdef EIGEN_VECTORIZE_SME
      sme_has_gebp_kernel<typename traits<Lhs_>::Scalar, typename traits<Rhs_>::Scalar>::value
          ? EIGEN_SME_FIXED_SIZE_GEMM_TO_COEFFBASED_THRESHOLD
          :
#endif
          EIGEN_FIXED_SIZE_GEMM_TO_COEFFBASED_THRESHOLD;

  static constexpr bool FixedSizeCoeffBasedProduct = ProductType == GemmProduct && Rows != Dynamic && Cols != Dynamic &&
                                                     Depth != Dynamic &&
                                                     (int(Rows) + int(Cols) + int(Depth) < FixedSizeThreshold);

  static constexpr int value = FixedSizeCoeffBasedProduct ? CoeffBasedProductMode : ProductType;
#ifdef EIGEN_DEBUG_PRODUCT
  static void debug() {
    const int rows_select = product_size_category<Rows, MaxRows>::value;
    const int cols_select = product_size_category<Cols, MaxCols>::value;
    const int depth_select = product_size_category<Depth, MaxDepth>::value;
    EIGEN_DEBUG_VAR(Rows);
    EIGEN_DEBUG_VAR(Cols);
    EIGEN_DEBUG_VAR(Depth);
    EIGEN_DEBUG_VAR(rows_select);
    EIGEN_DEBUG_VAR(cols_select);
    EIGEN_DEBUG_VAR(depth_select);
    EIGEN_DEBUG_VAR(ProductType);
    EIGEN_DEBUG_VAR(FixedSizeCoeffBasedProduct);
    EIGEN_DEBUG_VAR(value);
  }
#endif
};

/* The following allows to select the kind of product at compile time
 * based on the three dimensions of the product.
 * This is a compile time mapping from {1,Small,Large}^3 -> {product types} */
// FIXME: the current compile-time product-type mapping may not be optimal.
template <int M, int N>
struct product_type_selector<M, N, 1> : std::integral_constant<int, OuterProduct> {};
template <int M>
struct product_type_selector<M, 1, 1> : std::integral_constant<int, LazyCoeffBasedProductMode> {};
template <int N>
struct product_type_selector<1, N, 1> : std::integral_constant<int, LazyCoeffBasedProductMode> {};
template <int Depth>
struct product_type_selector<1, 1, Depth> : std::integral_constant<int, InnerProduct> {};
template <>
struct product_type_selector<1, 1, 1> : std::integral_constant<int, InnerProduct> {};
template <>
struct product_type_selector<Small, 1, Small> : std::integral_constant<int, CoeffBasedProductMode> {};
template <>
struct product_type_selector<1, Small, Small> : std::integral_constant<int, CoeffBasedProductMode> {};
template <>
struct product_type_selector<Small, Small, Small> : std::integral_constant<int, CoeffBasedProductMode> {};
template <>
struct product_type_selector<Small, Small, 1> : std::integral_constant<int, LazyCoeffBasedProductMode> {};
template <>
struct product_type_selector<Small, Large, 1> : std::integral_constant<int, LazyCoeffBasedProductMode> {};
template <>
struct product_type_selector<Large, Small, 1> : std::integral_constant<int, LazyCoeffBasedProductMode> {};
template <>
struct product_type_selector<1, Large, Small> : std::integral_constant<int, CoeffBasedProductMode> {};
template <>
struct product_type_selector<1, Large, Large> : std::integral_constant<int, GemvProduct> {};
template <>
struct product_type_selector<1, Small, Large> : std::integral_constant<int, CoeffBasedProductMode> {};
template <>
struct product_type_selector<Large, 1, Small> : std::integral_constant<int, CoeffBasedProductMode> {};
template <>
struct product_type_selector<Large, 1, Large> : std::integral_constant<int, GemvProduct> {};
template <>
struct product_type_selector<Small, 1, Large> : std::integral_constant<int, CoeffBasedProductMode> {};
template <>
struct product_type_selector<Small, Small, Large> : std::integral_constant<int, GemmProduct> {};
template <>
struct product_type_selector<Large, Small, Large> : std::integral_constant<int, GemmProduct> {};
template <>
struct product_type_selector<Small, Large, Large> : std::integral_constant<int, GemmProduct> {};
template <>
struct product_type_selector<Large, Large, Large> : std::integral_constant<int, GemmProduct> {};
template <>
struct product_type_selector<Large, Small, Small> : std::integral_constant<int, CoeffBasedProductMode> {};
template <>
struct product_type_selector<Small, Large, Small> : std::integral_constant<int, CoeffBasedProductMode> {};
template <>
struct product_type_selector<Large, Large, Small> : std::integral_constant<int, GemmProduct> {};

}  // end namespace internal

/***********************************************************************
 *  Implementation of Inner Vector Vector Product
 ***********************************************************************/

// FIXME: consider returning a Scalar instead of a 1x1 matrix for inner products.
// Pro: more natural for the user.
// Con: in a meta-unrolled algorithm a matrix-matrix product may reduce to a
// row-vector times column-vector product. To handle this, we could specialize
// Block<MatrixType,1,1> with operator=(Scalar x).

/***********************************************************************
 *  Implementation of Outer Vector Vector Product
 ***********************************************************************/

/***********************************************************************
 *  Implementation of General Matrix Vector Product
 ***********************************************************************/

/*  According to the shape/flags of the matrix we have to distinguish 3 different cases:
 *   1 - the matrix is col-major, BLAS compatible and M is large => call fast BLAS-like colmajor routine
 *   2 - the matrix is row-major, BLAS compatible and N is large => call fast BLAS-like rowmajor routine
 *   3 - all other cases are handled using a simple loop along the outer-storage direction.
 *  Therefore we need a lower level meta selector.
 *  Furthermore, if the matrix is the rhs, then the product has to be transposed.
 */
namespace internal {

template <int Side, int StorageOrder, bool BlasCompatible>
struct gemv_dense_selector;

}  // end namespace internal

namespace internal {

template <typename Scalar, int Size, int MaxSize, bool Cond>
struct gemv_static_vector_if;

template <typename Scalar, int Size, int MaxSize>
struct gemv_static_vector_if<Scalar, Size, MaxSize, false> {
  EIGEN_DEVICE_FUNC constexpr Scalar* data() {
    eigen_internal_assert(false && "should never be called");
    return 0;
  }
};

template <typename Scalar, int Size>
struct gemv_static_vector_if<Scalar, Size, Dynamic, true> {
  EIGEN_DEVICE_FUNC constexpr Scalar* data() { return 0; }
};

template <typename Scalar, int Size, int MaxSize>
struct gemv_static_vector_if<Scalar, Size, MaxSize, true> {
#if EIGEN_MAX_STATIC_ALIGN_BYTES != 0
  internal::plain_array<Scalar, internal::min_size_prefer_fixed(Size, MaxSize), 0, AlignedMax> m_data;
  constexpr Scalar* data() { return m_data.array; }
#else
  // Some architectures cannot align on the stack,
  // => let's manually enforce alignment by allocating more data and return the address of the first aligned element.
  internal::plain_array<Scalar, internal::min_size_prefer_fixed(Size, MaxSize) + EIGEN_MAX_ALIGN_BYTES, 0> m_data;
  constexpr Scalar* data() {
    return reinterpret_cast<Scalar*>((std::uintptr_t(m_data.array) & ~(std::size_t(EIGEN_MAX_ALIGN_BYTES - 1))) +
                                     EIGEN_MAX_ALIGN_BYTES);
  }
#endif
};

template <typename ResScalar>
using gemv_mapped_destination =
    Map<Matrix<ResScalar, Dynamic, 1>, plain_enum_min(AlignedMax, internal::packet_traits<ResScalar>::size)>;

template <typename Dest>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE gemv_mapped_destination<typename Dest::Scalar> gemv_construct_mapped_destination(
    Dest& dest, typename Dest::Scalar* actual_dest_ptr) {
#ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
  constexpr int Size = Dest::SizeAtCompileTime;
  Index size = dest.size();
  EIGEN_DENSE_STORAGE_CTOR_PLUGIN
#endif
  return gemv_mapped_destination<typename Dest::Scalar>(actual_dest_ptr, dest.size());
}

// Prepares the temporary destination shared by general, triangular, and selfadjoint GEMV kernels.
template <bool EvalToDest, typename Dest>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void gemv_prepare_destination(Dest& dest,
                                                                    typename Dest::Scalar* actual_dest_ptr) {
  EIGEN_IF_CONSTEXPR (!EvalToDest) gemv_construct_mapped_destination(dest, actual_dest_ptr) = dest;
}

template <bool EvalToDest, typename Dest>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void gemv_prepare_destination(Dest& dest, typename Dest::Scalar* actual_dest_ptr,
                                                                    bool initialize_to_zero) {
  if (initialize_to_zero) {
    gemv_construct_mapped_destination(dest, actual_dest_ptr).setZero();
  } else {
    gemv_prepare_destination<EvalToDest>(dest, actual_dest_ptr);
  }
}

template <bool EvalToDest, typename Dest>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void gemv_copy_destination(Dest& dest, typename Dest::Scalar* actual_dest_ptr) {
  EIGEN_IF_CONSTEXPR (!EvalToDest) {
    dest = gemv_mapped_destination<typename Dest::Scalar>(actual_dest_ptr, dest.size());
  }
}

// Adds complex-by-real scalar adaptation to the shared destination handling when necessary.
template <typename RhsScalar, typename ResScalar, bool EvalToDestAtCompileTime, bool ComplexByReal>
class gemv_destination_policy {
 public:
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE explicit gemv_destination_policy(const ResScalar& alpha) : alpha_(alpha) {}

  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool eval_to_dest() const {
    return EvalToDestAtCompileTime && alpha_is_compatible();
  }

  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE RhsScalar compatible_alpha() const {
    return alpha_is_compatible() ? get_factor<ResScalar, RhsScalar>::run(alpha_) : RhsScalar(1);
  }

  template <typename Dest>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void prepare(Dest& dest, ResScalar* actual_dest_ptr) const {
    gemv_prepare_destination<EvalToDestAtCompileTime>(dest, actual_dest_ptr, !alpha_is_compatible());
  }

  template <typename Dest>
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void copy_back(Dest& dest, ResScalar* actual_dest_ptr) const {
    if (!alpha_is_compatible()) {
      dest.matrix() += alpha_ * gemv_mapped_destination<ResScalar>(actual_dest_ptr, dest.size());
    } else {
      gemv_copy_destination<EvalToDestAtCompileTime>(dest, actual_dest_ptr);
    }
  }

 private:
  EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE bool alpha_is_compatible() const {
    return !ComplexByReal || numext::is_exactly_zero(numext::imag(alpha_));
  }

  ResScalar alpha_;
};

// Copies a non-contiguous GEMV right-hand side into the aligned buffer allocated by the caller.
template <bool DirectlyUseRhs, typename ActualRhsType>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void gemv_prepare_rhs(
    const ActualRhsType& rhs, typename remove_all_t<ActualRhsType>::Scalar* actual_rhs_ptr) {
  using ActualRhsTypeCleaned = remove_all_t<ActualRhsType>;
  EIGEN_IF_CONSTEXPR (!DirectlyUseRhs) {
#ifdef EIGEN_DENSE_STORAGE_CTOR_PLUGIN
    constexpr int Size = ActualRhsTypeCleaned::SizeAtCompileTime;
    Index size = rhs.size();
    EIGEN_DENSE_STORAGE_CTOR_PLUGIN
#endif
    // The caller obtains this temporary from ei_declare_aligned_stack_constructed_variable.
    Map<typename ActualRhsTypeCleaned::PlainObject, AlignedMax>(actual_rhs_ptr, rhs.size()) = rhs;
  }
}

// The vector is on the left => transposition
template <int StorageOrder, bool BlasCompatible>
struct gemv_dense_selector<OnTheLeft, StorageOrder, BlasCompatible> {
  template <typename Lhs, typename Rhs, typename Dest>
  static void run(const Lhs& lhs, const Rhs& rhs, Dest& dest, const typename Dest::Scalar& alpha) {
    Transpose<Dest> destT(dest);
    enum { OtherStorageOrder = StorageOrder == RowMajor ? ColMajor : RowMajor };
    gemv_dense_selector<OnTheRight, OtherStorageOrder, BlasCompatible>::run(rhs.transpose(), lhs.transpose(), destT,
                                                                            alpha);
  }
};

template <>
struct gemv_dense_selector<OnTheRight, ColMajor, true> {
  template <typename Lhs, typename Rhs, typename Dest>
  static inline void run(const Lhs& lhs, const Rhs& rhs, Dest& dest, const typename Dest::Scalar& alpha) {
    using LhsScalar = typename Lhs::Scalar;
    using RhsScalar = typename Rhs::Scalar;
    using ResScalar = typename Dest::Scalar;

    using LhsBlasTraits = internal::blas_traits<Lhs>;
    using ActualLhsType = typename LhsBlasTraits::DirectLinearAccessType;
    using RhsBlasTraits = internal::blas_traits<Rhs>;
    using ActualRhsType = typename RhsBlasTraits::DirectLinearAccessType;

    ActualLhsType actualLhs = LhsBlasTraits::extract(lhs);
    ActualRhsType actualRhs = RhsBlasTraits::extract(rhs);

    ResScalar actualAlpha = combine_scalar_factors(alpha, lhs, rhs);

    // make sure Dest is a compile-time vector type (bug 1166)
    using ActualDest = std::conditional_t<Dest::IsVectorAtCompileTime, Dest, typename Dest::ColXpr>;

    enum {
      // FIXME: find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1
      // on the other hand it is good for the cache to pack the vector anyways...
      EvalToDestAtCompileTime = (ActualDest::InnerStrideAtCompileTime == 1),
      ComplexByReal = (NumTraits<LhsScalar>::IsComplex) && (!NumTraits<RhsScalar>::IsComplex),
      MightCannotUseDest = ((!EvalToDestAtCompileTime) || ComplexByReal) && (ActualDest::MaxSizeAtCompileTime != 0)
    };

    using LhsMapper = const_blas_data_mapper<LhsScalar, Index, ColMajor>;
    using RhsMapper = const_blas_data_mapper<RhsScalar, Index, RowMajor>;
    EIGEN_IF_CONSTEXPR (!MightCannotUseDest) {
      // shortcut if we are sure to be able to use dest directly,
      // this eases the compiler to generate cleaner and more optimized code for most common cases
      general_matrix_vector_product<
          Index, LhsScalar, LhsMapper, ColMajor, LhsBlasTraits::NeedToConjugate, RhsScalar, RhsMapper,
          RhsBlasTraits::NeedToConjugate>::run(actualLhs.rows(), actualLhs.cols(),
                                               LhsMapper(actualLhs.data(), actualLhs.outerStride()),
                                               RhsMapper(actualRhs.data(), actualRhs.innerStride()), dest.data(), 1,
                                               get_factor<ResScalar, RhsScalar>::run(actualAlpha));
    } else {
      gemv_static_vector_if<ResScalar, ActualDest::SizeAtCompileTime, ActualDest::MaxSizeAtCompileTime,
                            MightCannotUseDest>
          static_dest;

      gemv_destination_policy<RhsScalar, ResScalar, EvalToDestAtCompileTime, ComplexByReal> destPolicy(actualAlpha);

      ei_declare_aligned_stack_constructed_variable(ResScalar, actualDestPtr, dest.size(),
                                                    destPolicy.eval_to_dest() ? dest.data() : static_dest.data());

      destPolicy.prepare(dest, actualDestPtr);

      general_matrix_vector_product<Index, LhsScalar, LhsMapper, ColMajor, LhsBlasTraits::NeedToConjugate, RhsScalar,
                                    RhsMapper, RhsBlasTraits::NeedToConjugate>::run(actualLhs.rows(), actualLhs.cols(),
                                                                                    LhsMapper(actualLhs.data(),
                                                                                              actualLhs.outerStride()),
                                                                                    RhsMapper(actualRhs.data(),
                                                                                              actualRhs.innerStride()),
                                                                                    actualDestPtr, 1,
                                                                                    destPolicy.compatible_alpha());

      destPolicy.copy_back(dest, actualDestPtr);
    }
  }
};

template <>
struct gemv_dense_selector<OnTheRight, RowMajor, true> {
  template <typename Lhs, typename Rhs, typename Dest>
  static void run(const Lhs& lhs, const Rhs& rhs, Dest& dest, const typename Dest::Scalar& alpha) {
    using LhsScalar = typename Lhs::Scalar;
    using RhsScalar = typename Rhs::Scalar;
    using ResScalar = typename Dest::Scalar;

    using LhsBlasTraits = internal::blas_traits<Lhs>;
    using ActualLhsType = typename LhsBlasTraits::DirectLinearAccessType;
    using RhsBlasTraits = internal::blas_traits<Rhs>;
    using ActualRhsType = typename RhsBlasTraits::DirectLinearAccessType;
    using ActualRhsTypeCleaned = internal::remove_all_t<ActualRhsType>;

    std::add_const_t<ActualLhsType> actualLhs = LhsBlasTraits::extract(lhs);
    std::add_const_t<ActualRhsType> actualRhs = RhsBlasTraits::extract(rhs);

    ResScalar actualAlpha = combine_scalar_factors(alpha, lhs, rhs);

    enum {
      // FIXME: find a way to allow an inner stride on the result if packet_traits<Scalar>::size==1
      // on the other hand it is good for the cache to pack the vector anyways...
      DirectlyUseRhs =
          ActualRhsTypeCleaned::InnerStrideAtCompileTime == 1 || ActualRhsTypeCleaned::MaxSizeAtCompileTime == 0
    };

    gemv_static_vector_if<RhsScalar, ActualRhsTypeCleaned::SizeAtCompileTime,
                          ActualRhsTypeCleaned::MaxSizeAtCompileTime, !DirectlyUseRhs>
        static_rhs;

    ei_declare_aligned_stack_constructed_variable(
        RhsScalar, actualRhsPtr, actualRhs.size(),
        DirectlyUseRhs ? const_cast<RhsScalar*>(actualRhs.data()) : static_rhs.data());

    gemv_prepare_rhs<DirectlyUseRhs>(actualRhs, actualRhsPtr);

    using LhsMapper = const_blas_data_mapper<LhsScalar, Index, RowMajor>;
    using RhsMapper = const_blas_data_mapper<RhsScalar, Index, ColMajor>;
    general_matrix_vector_product<Index, LhsScalar, LhsMapper, RowMajor, LhsBlasTraits::NeedToConjugate, RhsScalar,
                                  RhsMapper, RhsBlasTraits::NeedToConjugate>::
        run(actualLhs.rows(), actualLhs.cols(), LhsMapper(actualLhs.data(), actualLhs.outerStride()),
            RhsMapper(actualRhsPtr, 1), dest.data(),
            dest.col(0).innerStride(),  // NOTE  if dest is not a vector at compile-time, then dest.innerStride() might
                                        // be wrong. (bug 1166)
            actualAlpha);
  }
};

template <>
struct gemv_dense_selector<OnTheRight, ColMajor, false> {
  template <typename Lhs, typename Rhs, typename Dest>
  static void run(const Lhs& lhs, const Rhs& rhs, Dest& dest, const typename Dest::Scalar& alpha) {
    EIGEN_STATIC_ASSERT((!nested_eval<Lhs, 1>::Evaluate),
                        EIGEN_INTERNAL_COMPILATION_ERROR_OR_YOU_MADE_A_PROGRAMMING_MISTAKE);
    // TODO: if rhs is large enough it might be beneficial to make sure that dest is sequentially stored in memory,
    // otherwise use a temp
    typename nested_eval<Rhs, 1>::type actual_rhs(rhs);
    const Index size = rhs.rows();
    for (Index k = 0; k < size; ++k) dest += (alpha * actual_rhs.coeff(k)) * lhs.col(k);
  }
};

template <>
struct gemv_dense_selector<OnTheRight, RowMajor, false> {
  template <typename Lhs, typename Rhs, typename Dest>
  static void run(const Lhs& lhs, const Rhs& rhs, Dest& dest, const typename Dest::Scalar& alpha) {
    EIGEN_STATIC_ASSERT((!nested_eval<Lhs, 1>::Evaluate),
                        EIGEN_INTERNAL_COMPILATION_ERROR_OR_YOU_MADE_A_PROGRAMMING_MISTAKE);
    typename nested_eval<Rhs, Lhs::RowsAtCompileTime>::type actual_rhs(rhs);
    const Index rows = dest.rows();
    for (Index i = 0; i < rows; ++i)
      dest.coeffRef(i) += alpha * (lhs.row(i).cwiseProduct(actual_rhs.transpose())).sum();
  }
};

}  // end namespace internal

/***************************************************************************
 * Implementation of matrix base methods
 ***************************************************************************/

/** \returns the matrix product of \c *this and \a other.
 *
 * \note If instead of the matrix product you want the coefficient-wise product, see Cwise::operator*().
 *
 * \sa lazyProduct(), operator*=(const MatrixBase&), Cwise::operator*()
 */
template <typename Derived>
template <typename OtherDerived>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Product<Derived, OtherDerived> MatrixBase<Derived>::operator*(
    const MatrixBase<OtherDerived>& other) const {
  // A note regarding the function declaration: In MSVC, this function will sometimes
  // not be inlined since DenseStorage is an unwindable object for dynamic
  // matrices and product types are holding a member to store the result.
  // Thus it does not help tagging this function with EIGEN_STRONG_INLINE.
  enum {
    ProductIsValid = Derived::ColsAtCompileTime == Dynamic || OtherDerived::RowsAtCompileTime == Dynamic ||
                     int(Derived::ColsAtCompileTime) == int(OtherDerived::RowsAtCompileTime),
    AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime,
    SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived, OtherDerived)
  };
  // note to the lost user:
  //    * for a dot product use: v1.dot(v2)
  //    * for a coeff-wise product use: v1.cwiseProduct(v2)
  EIGEN_STATIC_ASSERT(
      ProductIsValid || !(AreVectors && SameSizes),
      INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
  EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
                      INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
  EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)
#ifdef EIGEN_DEBUG_PRODUCT
  internal::product_type<Derived, OtherDerived>::debug();
#endif

  return Product<Derived, OtherDerived>(derived(), other.derived());
}

/** \returns an expression of the matrix product of \c *this and \a other without implicit evaluation.
 *
 * The returned product will behave like any other expressions: the coefficients of the product will be
 * computed once at a time as requested. This might be useful in some extremely rare cases when only
 * a small and no coherent fraction of the result's coefficients have to be computed.
 *
 * \warning This version of the matrix product can be much much slower. So use it only if you know
 * what you are doing and that you measured a true speed improvement.
 *
 * \sa operator*(const MatrixBase&)
 */
template <typename Derived>
template <typename OtherDerived>
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Product<Derived, OtherDerived, LazyProduct>
MatrixBase<Derived>::lazyProduct(const MatrixBase<OtherDerived>& other) const {
  enum {
    ProductIsValid = Derived::ColsAtCompileTime == Dynamic || OtherDerived::RowsAtCompileTime == Dynamic ||
                     int(Derived::ColsAtCompileTime) == int(OtherDerived::RowsAtCompileTime),
    AreVectors = Derived::IsVectorAtCompileTime && OtherDerived::IsVectorAtCompileTime,
    SameSizes = EIGEN_PREDICATE_SAME_MATRIX_SIZE(Derived, OtherDerived)
  };
  // note to the lost user:
  //    * for a dot product use: v1.dot(v2)
  //    * for a coeff-wise product use: v1.cwiseProduct(v2)
  EIGEN_STATIC_ASSERT(
      ProductIsValid || !(AreVectors && SameSizes),
      INVALID_VECTOR_VECTOR_PRODUCT__IF_YOU_WANTED_A_DOT_OR_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTIONS)
  EIGEN_STATIC_ASSERT(ProductIsValid || !(SameSizes && !AreVectors),
                      INVALID_MATRIX_PRODUCT__IF_YOU_WANTED_A_COEFF_WISE_PRODUCT_YOU_MUST_USE_THE_EXPLICIT_FUNCTION)
  EIGEN_STATIC_ASSERT(ProductIsValid || SameSizes, INVALID_MATRIX_PRODUCT)

  return Product<Derived, OtherDerived, LazyProduct>(derived(), other.derived());
}

}  // end namespace Eigen

#endif  // EIGEN_GENERAL_PRODUCT_H