blob: c4fa771e222fc630b94733068df379d052959852 [file] [log] [blame]
// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2009 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/.
#ifndef EIGEN_GENERAL_BLOCK_PANEL_H
#define EIGEN_GENERAL_BLOCK_PANEL_H
// IWYU pragma: private
#include "../InternalHeaderCheck.h"
namespace Eigen {
namespace internal {
enum GEBPPacketSizeType { GEBPPacketFull = 0, GEBPPacketHalf, GEBPPacketQuarter };
template <typename LhsScalar_, typename RhsScalar_, bool ConjLhs_ = false, bool ConjRhs_ = false,
int Arch = Architecture::Target, int PacketSize_ = GEBPPacketFull>
class gebp_traits;
/** \internal \returns b if a<=0, and returns a otherwise. */
inline std::ptrdiff_t manage_caching_sizes_helper(std::ptrdiff_t a, std::ptrdiff_t b) { return a <= 0 ? b : a; }
#if defined(EIGEN_DEFAULT_L1_CACHE_SIZE)
#define EIGEN_SET_DEFAULT_L1_CACHE_SIZE(val) EIGEN_DEFAULT_L1_CACHE_SIZE
#else
#define EIGEN_SET_DEFAULT_L1_CACHE_SIZE(val) val
#endif // defined(EIGEN_DEFAULT_L1_CACHE_SIZE)
#if defined(EIGEN_DEFAULT_L2_CACHE_SIZE)
#define EIGEN_SET_DEFAULT_L2_CACHE_SIZE(val) EIGEN_DEFAULT_L2_CACHE_SIZE
#else
#define EIGEN_SET_DEFAULT_L2_CACHE_SIZE(val) val
#endif // defined(EIGEN_DEFAULT_L2_CACHE_SIZE)
#if defined(EIGEN_DEFAULT_L3_CACHE_SIZE)
#define EIGEN_SET_DEFAULT_L3_CACHE_SIZE(val) EIGEN_DEFAULT_L3_CACHE_SIZE
#else
#define EIGEN_SET_DEFAULT_L3_CACHE_SIZE(val) val
#endif // defined(EIGEN_DEFAULT_L3_CACHE_SIZE)
#if EIGEN_ARCH_i386_OR_x86_64
const std::ptrdiff_t defaultL1CacheSize = EIGEN_SET_DEFAULT_L1_CACHE_SIZE(32 * 1024);
const std::ptrdiff_t defaultL2CacheSize = EIGEN_SET_DEFAULT_L2_CACHE_SIZE(256 * 1024);
const std::ptrdiff_t defaultL3CacheSize = EIGEN_SET_DEFAULT_L3_CACHE_SIZE(2 * 1024 * 1024);
#elif EIGEN_ARCH_PPC
const std::ptrdiff_t defaultL1CacheSize = EIGEN_SET_DEFAULT_L1_CACHE_SIZE(64 * 1024);
#ifdef _ARCH_PWR10
const std::ptrdiff_t defaultL2CacheSize = EIGEN_SET_DEFAULT_L2_CACHE_SIZE(2 * 1024 * 1024);
const std::ptrdiff_t defaultL3CacheSize = EIGEN_SET_DEFAULT_L3_CACHE_SIZE(8 * 1024 * 1024);
#else
const std::ptrdiff_t defaultL2CacheSize = EIGEN_SET_DEFAULT_L2_CACHE_SIZE(512 * 1024);
const std::ptrdiff_t defaultL3CacheSize = EIGEN_SET_DEFAULT_L3_CACHE_SIZE(4 * 1024 * 1024);
#endif
#else
const std::ptrdiff_t defaultL1CacheSize = EIGEN_SET_DEFAULT_L1_CACHE_SIZE(16 * 1024);
const std::ptrdiff_t defaultL2CacheSize = EIGEN_SET_DEFAULT_L2_CACHE_SIZE(512 * 1024);
const std::ptrdiff_t defaultL3CacheSize = EIGEN_SET_DEFAULT_L3_CACHE_SIZE(512 * 1024);
#endif
#undef EIGEN_SET_DEFAULT_L1_CACHE_SIZE
#undef EIGEN_SET_DEFAULT_L2_CACHE_SIZE
#undef EIGEN_SET_DEFAULT_L3_CACHE_SIZE
/** \internal */
struct CacheSizes {
CacheSizes() : m_l1(-1), m_l2(-1), m_l3(-1) {
int l1CacheSize, l2CacheSize, l3CacheSize;
queryCacheSizes(l1CacheSize, l2CacheSize, l3CacheSize);
m_l1 = manage_caching_sizes_helper(l1CacheSize, defaultL1CacheSize);
m_l2 = manage_caching_sizes_helper(l2CacheSize, defaultL2CacheSize);
m_l3 = manage_caching_sizes_helper(l3CacheSize, defaultL3CacheSize);
}
std::ptrdiff_t m_l1;
std::ptrdiff_t m_l2;
std::ptrdiff_t m_l3;
};
/** \internal */
inline void manage_caching_sizes(Action action, std::ptrdiff_t* l1, std::ptrdiff_t* l2, std::ptrdiff_t* l3) {
static CacheSizes m_cacheSizes;
if (action == SetAction) {
// set the cpu cache size and cache all block sizes from a global cache size in byte
eigen_internal_assert(l1 != 0 && l2 != 0);
m_cacheSizes.m_l1 = *l1;
m_cacheSizes.m_l2 = *l2;
m_cacheSizes.m_l3 = *l3;
} else if (action == GetAction) {
eigen_internal_assert(l1 != 0 && l2 != 0);
*l1 = m_cacheSizes.m_l1;
*l2 = m_cacheSizes.m_l2;
*l3 = m_cacheSizes.m_l3;
} else {
eigen_internal_assert(false);
}
}
/* Helper for computeProductBlockingSizes.
*
* Given a m x k times k x n matrix product of scalar types \c LhsScalar and \c RhsScalar,
* this function computes the blocking size parameters along the respective dimensions
* for matrix products and related algorithms. The blocking sizes depends on various
* parameters:
* - the L1 and L2 cache sizes,
* - the register level blocking sizes defined by gebp_traits,
* - the number of scalars that fit into a packet (when vectorization is enabled).
*
* \sa setCpuCacheSizes */
template <typename LhsScalar, typename RhsScalar, int KcFactor, typename Index>
void evaluateProductBlockingSizesHeuristic(Index& k, Index& m, Index& n, Index num_threads = 1) {
typedef gebp_traits<LhsScalar, RhsScalar> Traits;
// Explanations:
// Let's recall that the product algorithms form mc x kc vertical panels A' on the lhs and
// kc x nc blocks B' on the rhs. B' has to fit into L2/L3 cache. Moreover, A' is processed
// per mr x kc horizontal small panels where mr is the blocking size along the m dimension
// at the register level. This small horizontal panel has to stay within L1 cache.
std::ptrdiff_t l1, l2, l3;
manage_caching_sizes(GetAction, &l1, &l2, &l3);
#ifdef EIGEN_VECTORIZE_AVX512
// We need to find a rationale for that, but without this adjustment,
// performance with AVX512 is pretty bad, like -20% slower.
// One reason is that with increasing packet-size, the blocking size k
// has to become pretty small if we want that 1 lhs panel fit within L1.
// For instance, with the 3pX4 kernel and double, the size of the lhs+rhs panels are:
// k*(3*64 + 4*8) Bytes, with l1=32kBytes, and k%8=0, we have k=144.
// This is quite small for a good reuse of the accumulation registers.
l1 *= 4;
#endif
if (num_threads > 1) {
typedef typename Traits::ResScalar ResScalar;
enum {
kdiv = KcFactor * (Traits::mr * sizeof(LhsScalar) + Traits::nr * sizeof(RhsScalar)),
ksub = Traits::mr * (Traits::nr * sizeof(ResScalar)),
kr = 8,
mr = Traits::mr,
nr = Traits::nr
};
// Increasing k gives us more time to prefetch the content of the "C"
// registers. However once the latency is hidden there is no point in
// increasing the value of k, so we'll cap it at 320 (value determined
// experimentally).
// To avoid that k vanishes, we make k_cache at least as big as kr
const Index k_cache = numext::maxi<Index>(kr, (numext::mini<Index>)((l1 - ksub) / kdiv, 320));
if (k_cache < k) {
k = k_cache - (k_cache % kr);
eigen_internal_assert(k > 0);
}
const Index n_cache = (l2 - l1) / (nr * sizeof(RhsScalar) * k);
const Index n_per_thread = numext::div_ceil(n, num_threads);
if (n_cache <= n_per_thread) {
// Don't exceed the capacity of the l2 cache.
eigen_internal_assert(n_cache >= static_cast<Index>(nr));
n = n_cache - (n_cache % nr);
eigen_internal_assert(n > 0);
} else {
n = (numext::mini<Index>)(n, (n_per_thread + nr - 1) - ((n_per_thread + nr - 1) % nr));
}
if (l3 > l2) {
// l3 is shared between all cores, so we'll give each thread its own chunk of l3.
const Index m_cache = (l3 - l2) / (sizeof(LhsScalar) * k * num_threads);
const Index m_per_thread = numext::div_ceil(m, num_threads);
if (m_cache < m_per_thread && m_cache >= static_cast<Index>(mr)) {
m = m_cache - (m_cache % mr);
eigen_internal_assert(m > 0);
} else {
m = (numext::mini<Index>)(m, (m_per_thread + mr - 1) - ((m_per_thread + mr - 1) % mr));
}
}
} else {
// In unit tests we do not want to use extra large matrices,
// so we reduce the cache size to check the blocking strategy is not flawed
#ifdef EIGEN_DEBUG_SMALL_PRODUCT_BLOCKS
l1 = 9 * 1024;
l2 = 32 * 1024;
l3 = 512 * 1024;
#endif
// Early return for small problems because the computation below are time consuming for small problems.
// Perhaps it would make more sense to consider k*n*m??
// Note that for very tiny problem, this function should be bypassed anyway
// because we use the coefficient-based implementation for them.
if ((numext::maxi)(k, (numext::maxi)(m, n)) < 48) return;
typedef typename Traits::ResScalar ResScalar;
enum {
k_peeling = 8,
k_div = KcFactor * (Traits::mr * sizeof(LhsScalar) + Traits::nr * sizeof(RhsScalar)),
k_sub = Traits::mr * (Traits::nr * sizeof(ResScalar))
};
// ---- 1st level of blocking on L1, yields kc ----
// Blocking on the third dimension (i.e., k) is chosen so that an horizontal panel
// of size mr x kc of the lhs plus a vertical panel of kc x nr of the rhs both fits within L1 cache.
// We also include a register-level block of the result (mx x nr).
// (In an ideal world only the lhs panel would stay in L1)
// Moreover, kc has to be a multiple of 8 to be compatible with loop peeling, leading to a maximum blocking size of:
const Index max_kc = numext::maxi<Index>(((l1 - k_sub) / k_div) & (~(k_peeling - 1)), 1);
const Index old_k = k;
if (k > max_kc) {
// We are really blocking on the third dimension:
// -> reduce blocking size to make sure the last block is as large as possible
// while keeping the same number of sweeps over the result.
k = (k % max_kc) == 0 ? max_kc
: max_kc - k_peeling * ((max_kc - 1 - (k % max_kc)) / (k_peeling * (k / max_kc + 1)));
eigen_internal_assert(((old_k / k) == (old_k / max_kc)) && "the number of sweeps has to remain the same");
}
// ---- 2nd level of blocking on max(L2,L3), yields nc ----
// TODO find a reliable way to get the actual amount of cache per core to use for 2nd level blocking, that is:
// actual_l2 = max(l2, l3/nb_core_sharing_l3)
// The number below is quite conservative: it is better to underestimate the cache size rather than overestimating it)
// For instance, it corresponds to 6MB of L3 shared among 4 cores.
#ifdef EIGEN_DEBUG_SMALL_PRODUCT_BLOCKS
const Index actual_l2 = l3;
#else
const Index actual_l2 = 1572864; // == 1.5 MB
#endif
// Here, nc is chosen such that a block of kc x nc of the rhs fit within half of L2.
// The second half is implicitly reserved to access the result and lhs coefficients.
// When k<max_kc, then nc can arbitrarily growth. In practice, it seems to be fruitful
// to limit this growth: we bound nc to growth by a factor x1.5.
// However, if the entire lhs block fit within L1, then we are not going to block on the rows at all,
// and it becomes fruitful to keep the packed rhs blocks in L1 if there is enough remaining space.
Index max_nc;
const Index lhs_bytes = m * k * sizeof(LhsScalar);
const Index remaining_l1 = l1 - k_sub - lhs_bytes;
if (remaining_l1 >= Index(Traits::nr * sizeof(RhsScalar)) * k) {
// L1 blocking
max_nc = remaining_l1 / (k * sizeof(RhsScalar));
} else {
// L2 blocking
max_nc = (3 * actual_l2) / (2 * 2 * max_kc * sizeof(RhsScalar));
}
// WARNING Below, we assume that Traits::nr is a power of two.
Index nc = numext::mini<Index>(actual_l2 / (2 * k * sizeof(RhsScalar)), max_nc) & (~(Traits::nr - 1));
if (n > nc) {
// We are really blocking over the columns:
// -> reduce blocking size to make sure the last block is as large as possible
// while keeping the same number of sweeps over the packed lhs.
// Here we allow one more sweep if this gives us a perfect match, thus the commented "-1"
n = (n % nc) == 0 ? nc : (nc - Traits::nr * ((nc /*-1*/ - (n % nc)) / (Traits::nr * (n / nc + 1))));
} else if (old_k == k) {
// So far, no blocking at all, i.e., kc==k, and nc==n.
// In this case, let's perform a blocking over the rows such that the packed lhs data is kept in cache L1/L2
// TODO: part of this blocking strategy is now implemented within the kernel itself, so the L1-based heuristic
// here should be obsolete.
Index problem_size = k * n * sizeof(LhsScalar);
Index actual_lm = actual_l2;
Index max_mc = m;
if (problem_size <= 1024) {
// problem is small enough to keep in L1
// Let's choose m such that lhs's block fit in 1/3 of L1
actual_lm = l1;
} else if (l3 != 0 && problem_size <= 32768) {
// we have both L2 and L3, and problem is small enough to be kept in L2
// Let's choose m such that lhs's block fit in 1/3 of L2
actual_lm = l2;
max_mc = (numext::mini<Index>)(576, max_mc);
}
Index mc = (numext::mini<Index>)(actual_lm / (3 * k * sizeof(LhsScalar)), max_mc);
if (mc > Traits::mr)
mc -= mc % Traits::mr;
else if (mc == 0)
return;
m = (m % mc) == 0 ? mc : (mc - Traits::mr * ((mc /*-1*/ - (m % mc)) / (Traits::mr * (m / mc + 1))));
}
}
}
template <typename Index>
inline bool useSpecificBlockingSizes(Index& k, Index& m, Index& n) {
#ifdef EIGEN_TEST_SPECIFIC_BLOCKING_SIZES
if (EIGEN_TEST_SPECIFIC_BLOCKING_SIZES) {
k = numext::mini<Index>(k, EIGEN_TEST_SPECIFIC_BLOCKING_SIZE_K);
m = numext::mini<Index>(m, EIGEN_TEST_SPECIFIC_BLOCKING_SIZE_M);
n = numext::mini<Index>(n, EIGEN_TEST_SPECIFIC_BLOCKING_SIZE_N);
return true;
}
#else
EIGEN_UNUSED_VARIABLE(k)
EIGEN_UNUSED_VARIABLE(m)
EIGEN_UNUSED_VARIABLE(n)
#endif
return false;
}
/** \brief Computes the blocking parameters for a m x k times k x n matrix product
*
* \param[in,out] k Input: the third dimension of the product. Output: the blocking size along the same dimension.
* \param[in,out] m Input: the number of rows of the left hand side. Output: the blocking size along the same dimension.
* \param[in,out] n Input: the number of columns of the right hand side. Output: the blocking size along the same
* dimension.
*
* Given a m x k times k x n matrix product of scalar types \c LhsScalar and \c RhsScalar,
* this function computes the blocking size parameters along the respective dimensions
* for matrix products and related algorithms.
*
* The blocking size parameters may be evaluated:
* - either by a heuristic based on cache sizes;
* - or using fixed prescribed values (for testing purposes).
*
* \sa setCpuCacheSizes */
template <typename LhsScalar, typename RhsScalar, int KcFactor, typename Index>
void computeProductBlockingSizes(Index& k, Index& m, Index& n, Index num_threads = 1) {
if (!useSpecificBlockingSizes(k, m, n)) {
evaluateProductBlockingSizesHeuristic<LhsScalar, RhsScalar, KcFactor, Index>(k, m, n, num_threads);
}
}
template <typename LhsScalar, typename RhsScalar, typename Index>
inline void computeProductBlockingSizes(Index& k, Index& m, Index& n, Index num_threads = 1) {
computeProductBlockingSizes<LhsScalar, RhsScalar, 1, Index>(k, m, n, num_threads);
}
template <typename RhsPacket, typename RhsPacketx4, int registers_taken>
struct RhsPanelHelper {
private:
static constexpr int remaining_registers =
(std::max)(int(EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS) - registers_taken, 0);
public:
typedef std::conditional_t<remaining_registers >= 4, RhsPacketx4, RhsPacket> type;
};
template <typename Packet>
struct QuadPacket {
Packet B_0, B1, B2, B3;
const Packet& get(const FixedInt<0>&) const { return B_0; }
const Packet& get(const FixedInt<1>&) const { return B1; }
const Packet& get(const FixedInt<2>&) const { return B2; }
const Packet& get(const FixedInt<3>&) const { return B3; }
};
template <int N, typename T1, typename T2, typename T3>
struct packet_conditional {
typedef T3 type;
};
template <typename T1, typename T2, typename T3>
struct packet_conditional<GEBPPacketFull, T1, T2, T3> {
typedef T1 type;
};
template <typename T1, typename T2, typename T3>
struct packet_conditional<GEBPPacketHalf, T1, T2, T3> {
typedef T2 type;
};
#define PACKET_DECL_COND_POSTFIX(postfix, name, packet_size) \
typedef typename packet_conditional< \
packet_size, typename packet_traits<name##Scalar>::type, typename packet_traits<name##Scalar>::half, \
typename unpacket_traits<typename packet_traits<name##Scalar>::half>::half>::type name##Packet##postfix
#define PACKET_DECL_COND(name, packet_size) \
typedef typename packet_conditional< \
packet_size, typename packet_traits<name##Scalar>::type, typename packet_traits<name##Scalar>::half, \
typename unpacket_traits<typename packet_traits<name##Scalar>::half>::half>::type name##Packet
#define PACKET_DECL_COND_SCALAR_POSTFIX(postfix, packet_size) \
typedef typename packet_conditional< \
packet_size, typename packet_traits<Scalar>::type, typename packet_traits<Scalar>::half, \
typename unpacket_traits<typename packet_traits<Scalar>::half>::half>::type ScalarPacket##postfix
#define PACKET_DECL_COND_SCALAR(packet_size) \
typedef typename packet_conditional< \
packet_size, typename packet_traits<Scalar>::type, typename packet_traits<Scalar>::half, \
typename unpacket_traits<typename packet_traits<Scalar>::half>::half>::type ScalarPacket
/* Vectorization logic
* real*real: unpack rhs to constant packets, ...
*
* cd*cd : unpack rhs to (b_r,b_r), (b_i,b_i), mul to get (a_r b_r,a_i b_r) (a_r b_i,a_i b_i),
* storing each res packet into two packets (2x2),
* at the end combine them: swap the second and addsub them
* cf*cf : same but with 2x4 blocks
* cplx*real : unpack rhs to constant packets, ...
* real*cplx : load lhs as (a0,a0,a1,a1), and mul as usual
*/
template <typename LhsScalar_, typename RhsScalar_, bool ConjLhs_, bool ConjRhs_, int Arch, int PacketSize_>
class gebp_traits {
public:
typedef LhsScalar_ LhsScalar;
typedef RhsScalar_ RhsScalar;
typedef typename ScalarBinaryOpTraits<LhsScalar, RhsScalar>::ReturnType ResScalar;
PACKET_DECL_COND_POSTFIX(_, Lhs, PacketSize_);
PACKET_DECL_COND_POSTFIX(_, Rhs, PacketSize_);
PACKET_DECL_COND_POSTFIX(_, Res, PacketSize_);
enum {
ConjLhs = ConjLhs_,
ConjRhs = ConjRhs_,
Vectorizable = unpacket_traits<LhsPacket_>::vectorizable && unpacket_traits<RhsPacket_>::vectorizable,
LhsPacketSize = Vectorizable ? unpacket_traits<LhsPacket_>::size : 1,
RhsPacketSize = Vectorizable ? unpacket_traits<RhsPacket_>::size : 1,
ResPacketSize = Vectorizable ? unpacket_traits<ResPacket_>::size : 1,
NumberOfRegisters = EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS,
// register block size along the N direction must be 1 or 4
nr = 4,
// register block size along the M direction (currently, this one cannot be modified)
default_mr = (plain_enum_min(16, NumberOfRegisters) / 2 / nr) * LhsPacketSize,
#if defined(EIGEN_HAS_SINGLE_INSTRUCTION_MADD) && !defined(EIGEN_VECTORIZE_ALTIVEC) && \
!defined(EIGEN_VECTORIZE_VSX) && ((!EIGEN_COMP_MSVC) || (EIGEN_COMP_MSVC >= 1914))
// we assume 16 registers or more
// See bug 992, if the scalar type is not vectorizable but that EIGEN_HAS_SINGLE_INSTRUCTION_MADD is defined,
// then using 3*LhsPacketSize triggers non-implemented paths in syrk.
// Bug 1515: MSVC prior to v19.14 yields to register spilling.
mr = Vectorizable ? 3 * LhsPacketSize : default_mr,
#else
mr = default_mr,
#endif
LhsProgress = LhsPacketSize,
RhsProgress = 1
};
typedef std::conditional_t<Vectorizable, LhsPacket_, LhsScalar> LhsPacket;
typedef std::conditional_t<Vectorizable, RhsPacket_, RhsScalar> RhsPacket;
typedef std::conditional_t<Vectorizable, ResPacket_, ResScalar> ResPacket;
typedef LhsPacket LhsPacket4Packing;
typedef QuadPacket<RhsPacket> RhsPacketx4;
typedef ResPacket AccPacket;
EIGEN_STRONG_INLINE void initAcc(AccPacket& p) { p = pset1<ResPacket>(ResScalar(0)); }
template <typename RhsPacketType>
EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacketType& dest) const {
dest = pset1<RhsPacketType>(*b);
}
EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacketx4& dest) const {
pbroadcast4(b, dest.B_0, dest.B1, dest.B2, dest.B3);
}
template <typename RhsPacketType>
EIGEN_STRONG_INLINE void updateRhs(const RhsScalar* b, RhsPacketType& dest) const {
loadRhs(b, dest);
}
EIGEN_STRONG_INLINE void updateRhs(const RhsScalar*, RhsPacketx4&) const {}
EIGEN_STRONG_INLINE void loadRhsQuad(const RhsScalar* b, RhsPacket& dest) const { dest = ploadquad<RhsPacket>(b); }
template <typename LhsPacketType>
EIGEN_STRONG_INLINE void loadLhs(const LhsScalar* a, LhsPacketType& dest) const {
dest = pload<LhsPacketType>(a);
}
template <typename LhsPacketType>
EIGEN_STRONG_INLINE void loadLhsUnaligned(const LhsScalar* a, LhsPacketType& dest) const {
dest = ploadu<LhsPacketType>(a);
}
template <typename LhsPacketType, typename RhsPacketType, typename AccPacketType, typename LaneIdType>
EIGEN_STRONG_INLINE void madd(const LhsPacketType& a, const RhsPacketType& b, AccPacketType& c, RhsPacketType& tmp,
const LaneIdType&) const {
conj_helper<LhsPacketType, RhsPacketType, ConjLhs, ConjRhs> cj;
// It would be a lot cleaner to call pmadd all the time. Unfortunately if we
// let gcc allocate the register in which to store the result of the pmul
// (in the case where there is no FMA) gcc fails to figure out how to avoid
// spilling register.
#ifdef EIGEN_HAS_SINGLE_INSTRUCTION_MADD
EIGEN_UNUSED_VARIABLE(tmp);
c = cj.pmadd(a, b, c);
#else
tmp = b;
tmp = cj.pmul(a, tmp);
c = padd(c, tmp);
#endif
}
template <typename LhsPacketType, typename AccPacketType, typename LaneIdType>
EIGEN_STRONG_INLINE void madd(const LhsPacketType& a, const RhsPacketx4& b, AccPacketType& c, RhsPacket& tmp,
const LaneIdType& lane) const {
madd(a, b.get(lane), c, tmp, lane);
}
EIGEN_STRONG_INLINE void acc(const AccPacket& c, const ResPacket& alpha, ResPacket& r) const {
r = pmadd(c, alpha, r);
}
template <typename ResPacketHalf>
EIGEN_STRONG_INLINE void acc(const ResPacketHalf& c, const ResPacketHalf& alpha, ResPacketHalf& r) const {
r = pmadd(c, alpha, r);
}
};
template <typename RealScalar, bool ConjLhs_, int Arch, int PacketSize_>
class gebp_traits<std::complex<RealScalar>, RealScalar, ConjLhs_, false, Arch, PacketSize_> {
public:
typedef std::complex<RealScalar> LhsScalar;
typedef RealScalar RhsScalar;
typedef typename ScalarBinaryOpTraits<LhsScalar, RhsScalar>::ReturnType ResScalar;
PACKET_DECL_COND_POSTFIX(_, Lhs, PacketSize_);
PACKET_DECL_COND_POSTFIX(_, Rhs, PacketSize_);
PACKET_DECL_COND_POSTFIX(_, Res, PacketSize_);
enum {
ConjLhs = ConjLhs_,
ConjRhs = false,
Vectorizable = unpacket_traits<LhsPacket_>::vectorizable && unpacket_traits<RhsPacket_>::vectorizable,
LhsPacketSize = Vectorizable ? unpacket_traits<LhsPacket_>::size : 1,
RhsPacketSize = Vectorizable ? unpacket_traits<RhsPacket_>::size : 1,
ResPacketSize = Vectorizable ? unpacket_traits<ResPacket_>::size : 1,
NumberOfRegisters = EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS,
nr = 4,
#if defined(EIGEN_HAS_SINGLE_INSTRUCTION_MADD) && !defined(EIGEN_VECTORIZE_ALTIVEC) && !defined(EIGEN_VECTORIZE_VSX)
// we assume 16 registers
mr = 3 * LhsPacketSize,
#else
mr = (plain_enum_min(16, NumberOfRegisters) / 2 / nr) * LhsPacketSize,
#endif
LhsProgress = LhsPacketSize,
RhsProgress = 1
};
typedef std::conditional_t<Vectorizable, LhsPacket_, LhsScalar> LhsPacket;
typedef std::conditional_t<Vectorizable, RhsPacket_, RhsScalar> RhsPacket;
typedef std::conditional_t<Vectorizable, ResPacket_, ResScalar> ResPacket;
typedef LhsPacket LhsPacket4Packing;
typedef QuadPacket<RhsPacket> RhsPacketx4;
typedef ResPacket AccPacket;
EIGEN_STRONG_INLINE void initAcc(AccPacket& p) { p = pset1<ResPacket>(ResScalar(0)); }
template <typename RhsPacketType>
EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacketType& dest) const {
dest = pset1<RhsPacketType>(*b);
}
EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacketx4& dest) const {
pbroadcast4(b, dest.B_0, dest.B1, dest.B2, dest.B3);
}
template <typename RhsPacketType>
EIGEN_STRONG_INLINE void updateRhs(const RhsScalar* b, RhsPacketType& dest) const {
loadRhs(b, dest);
}
EIGEN_STRONG_INLINE void updateRhs(const RhsScalar*, RhsPacketx4&) const {}
EIGEN_STRONG_INLINE void loadRhsQuad(const RhsScalar* b, RhsPacket& dest) const {
loadRhsQuad_impl(b, dest, std::conditional_t<RhsPacketSize == 16, true_type, false_type>());
}
EIGEN_STRONG_INLINE void loadRhsQuad_impl(const RhsScalar* b, RhsPacket& dest, const true_type&) const {
// FIXME we can do better!
// what we want here is a ploadheight
RhsScalar tmp[4] = {b[0], b[0], b[1], b[1]};
dest = ploadquad<RhsPacket>(tmp);
}
EIGEN_STRONG_INLINE void loadRhsQuad_impl(const RhsScalar* b, RhsPacket& dest, const false_type&) const {
eigen_internal_assert(RhsPacketSize <= 8);
dest = pset1<RhsPacket>(*b);
}
EIGEN_STRONG_INLINE void loadLhs(const LhsScalar* a, LhsPacket& dest) const { dest = pload<LhsPacket>(a); }
template <typename LhsPacketType>
EIGEN_STRONG_INLINE void loadLhsUnaligned(const LhsScalar* a, LhsPacketType& dest) const {
dest = ploadu<LhsPacketType>(a);
}
template <typename LhsPacketType, typename RhsPacketType, typename AccPacketType, typename LaneIdType>
EIGEN_STRONG_INLINE void madd(const LhsPacketType& a, const RhsPacketType& b, AccPacketType& c, RhsPacketType& tmp,
const LaneIdType&) const {
madd_impl(a, b, c, tmp, std::conditional_t<Vectorizable, true_type, false_type>());
}
template <typename LhsPacketType, typename RhsPacketType, typename AccPacketType>
EIGEN_STRONG_INLINE void madd_impl(const LhsPacketType& a, const RhsPacketType& b, AccPacketType& c,
RhsPacketType& tmp, const true_type&) const {
#ifdef EIGEN_HAS_SINGLE_INSTRUCTION_MADD
EIGEN_UNUSED_VARIABLE(tmp);
c.v = pmadd(a.v, b, c.v);
#else
tmp = b;
tmp = pmul(a.v, tmp);
c.v = padd(c.v, tmp);
#endif
}
EIGEN_STRONG_INLINE void madd_impl(const LhsScalar& a, const RhsScalar& b, ResScalar& c, RhsScalar& /*tmp*/,
const false_type&) const {
c += a * b;
}
template <typename LhsPacketType, typename AccPacketType, typename LaneIdType>
EIGEN_STRONG_INLINE void madd(const LhsPacketType& a, const RhsPacketx4& b, AccPacketType& c, RhsPacket& tmp,
const LaneIdType& lane) const {
madd(a, b.get(lane), c, tmp, lane);
}
template <typename ResPacketType, typename AccPacketType>
EIGEN_STRONG_INLINE void acc(const AccPacketType& c, const ResPacketType& alpha, ResPacketType& r) const {
conj_helper<ResPacketType, ResPacketType, ConjLhs, false> cj;
r = cj.pmadd(c, alpha, r);
}
protected:
};
template <typename Packet>
struct DoublePacket {
Packet first;
Packet second;
};
template <typename Packet>
DoublePacket<Packet> padd(const DoublePacket<Packet>& a, const DoublePacket<Packet>& b) {
DoublePacket<Packet> res;
res.first = padd(a.first, b.first);
res.second = padd(a.second, b.second);
return res;
}
// note that for DoublePacket<RealPacket> the "4" in "downto4"
// corresponds to the number of complexes, so it means "8"
// it terms of real coefficients.
template <typename Packet>
const DoublePacket<Packet>& predux_half_dowto4(const DoublePacket<Packet>& a,
std::enable_if_t<unpacket_traits<Packet>::size <= 8>* = 0) {
return a;
}
template <typename Packet>
DoublePacket<typename unpacket_traits<Packet>::half> predux_half_dowto4(
const DoublePacket<Packet>& a, std::enable_if_t<unpacket_traits<Packet>::size == 16>* = 0) {
// yes, that's pretty hackish :(
DoublePacket<typename unpacket_traits<Packet>::half> res;
typedef std::complex<typename unpacket_traits<Packet>::type> Cplx;
typedef typename packet_traits<Cplx>::type CplxPacket;
res.first = predux_half_dowto4(CplxPacket(a.first)).v;
res.second = predux_half_dowto4(CplxPacket(a.second)).v;
return res;
}
// same here, "quad" actually means "8" in terms of real coefficients
template <typename Scalar, typename RealPacket>
void loadQuadToDoublePacket(const Scalar* b, DoublePacket<RealPacket>& dest,
std::enable_if_t<unpacket_traits<RealPacket>::size <= 8>* = 0) {
dest.first = pset1<RealPacket>(numext::real(*b));
dest.second = pset1<RealPacket>(numext::imag(*b));
}
template <typename Scalar, typename RealPacket>
void loadQuadToDoublePacket(const Scalar* b, DoublePacket<RealPacket>& dest,
std::enable_if_t<unpacket_traits<RealPacket>::size == 16>* = 0) {
// yes, that's pretty hackish too :(
typedef typename NumTraits<Scalar>::Real RealScalar;
RealScalar r[4] = {numext::real(b[0]), numext::real(b[0]), numext::real(b[1]), numext::real(b[1])};
RealScalar i[4] = {numext::imag(b[0]), numext::imag(b[0]), numext::imag(b[1]), numext::imag(b[1])};
dest.first = ploadquad<RealPacket>(r);
dest.second = ploadquad<RealPacket>(i);
}
template <typename Packet>
struct unpacket_traits<DoublePacket<Packet> > {
typedef DoublePacket<typename unpacket_traits<Packet>::half> half;
enum { size = 2 * unpacket_traits<Packet>::size };
};
// template<typename Packet>
// DoublePacket<Packet> pmadd(const DoublePacket<Packet> &a, const DoublePacket<Packet> &b)
// {
// DoublePacket<Packet> res;
// res.first = padd(a.first, b.first);
// res.second = padd(a.second,b.second);
// return res;
// }
template <typename RealScalar, bool ConjLhs_, bool ConjRhs_, int Arch, int PacketSize_>
class gebp_traits<std::complex<RealScalar>, std::complex<RealScalar>, ConjLhs_, ConjRhs_, Arch, PacketSize_> {
public:
typedef std::complex<RealScalar> Scalar;
typedef std::complex<RealScalar> LhsScalar;
typedef std::complex<RealScalar> RhsScalar;
typedef std::complex<RealScalar> ResScalar;
PACKET_DECL_COND_POSTFIX(_, Lhs, PacketSize_);
PACKET_DECL_COND_POSTFIX(_, Rhs, PacketSize_);
PACKET_DECL_COND_POSTFIX(_, Res, PacketSize_);
PACKET_DECL_COND(Real, PacketSize_);
PACKET_DECL_COND_SCALAR(PacketSize_);
enum {
ConjLhs = ConjLhs_,
ConjRhs = ConjRhs_,
Vectorizable = unpacket_traits<RealPacket>::vectorizable && unpacket_traits<ScalarPacket>::vectorizable,
ResPacketSize = Vectorizable ? unpacket_traits<ResPacket_>::size : 1,
LhsPacketSize = Vectorizable ? unpacket_traits<LhsPacket_>::size : 1,
RhsPacketSize = Vectorizable ? unpacket_traits<RhsScalar>::size : 1,
RealPacketSize = Vectorizable ? unpacket_traits<RealPacket>::size : 1,
// FIXME: should depend on NumberOfRegisters
nr = 4,
mr = ResPacketSize,
LhsProgress = ResPacketSize,
RhsProgress = 1
};
typedef DoublePacket<RealPacket> DoublePacketType;
typedef std::conditional_t<Vectorizable, ScalarPacket, Scalar> LhsPacket4Packing;
typedef std::conditional_t<Vectorizable, RealPacket, Scalar> LhsPacket;
typedef std::conditional_t<Vectorizable, DoublePacketType, Scalar> RhsPacket;
typedef std::conditional_t<Vectorizable, ScalarPacket, Scalar> ResPacket;
typedef std::conditional_t<Vectorizable, DoublePacketType, Scalar> AccPacket;
// this actually holds 8 packets!
typedef QuadPacket<RhsPacket> RhsPacketx4;
EIGEN_STRONG_INLINE void initAcc(Scalar& p) { p = Scalar(0); }
EIGEN_STRONG_INLINE void initAcc(DoublePacketType& p) {
p.first = pset1<RealPacket>(RealScalar(0));
p.second = pset1<RealPacket>(RealScalar(0));
}
// Scalar path
EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, ScalarPacket& dest) const { dest = pset1<ScalarPacket>(*b); }
// Vectorized path
template <typename RealPacketType>
EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, DoublePacket<RealPacketType>& dest) const {
dest.first = pset1<RealPacketType>(numext::real(*b));
dest.second = pset1<RealPacketType>(numext::imag(*b));
}
EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacketx4& dest) const {
loadRhs(b, dest.B_0);
loadRhs(b + 1, dest.B1);
loadRhs(b + 2, dest.B2);
loadRhs(b + 3, dest.B3);
}
// Scalar path
EIGEN_STRONG_INLINE void updateRhs(const RhsScalar* b, ScalarPacket& dest) const { loadRhs(b, dest); }
// Vectorized path
template <typename RealPacketType>
EIGEN_STRONG_INLINE void updateRhs(const RhsScalar* b, DoublePacket<RealPacketType>& dest) const {
loadRhs(b, dest);
}
EIGEN_STRONG_INLINE void updateRhs(const RhsScalar*, RhsPacketx4&) const {}
EIGEN_STRONG_INLINE void loadRhsQuad(const RhsScalar* b, ResPacket& dest) const { loadRhs(b, dest); }
EIGEN_STRONG_INLINE void loadRhsQuad(const RhsScalar* b, DoublePacketType& dest) const {
loadQuadToDoublePacket(b, dest);
}
// nothing special here
EIGEN_STRONG_INLINE void loadLhs(const LhsScalar* a, LhsPacket& dest) const {
dest = pload<LhsPacket>((const typename unpacket_traits<LhsPacket>::type*)(a));
}
template <typename LhsPacketType>
EIGEN_STRONG_INLINE void loadLhsUnaligned(const LhsScalar* a, LhsPacketType& dest) const {
dest = ploadu<LhsPacketType>((const typename unpacket_traits<LhsPacketType>::type*)(a));
}
template <typename LhsPacketType, typename RhsPacketType, typename ResPacketType, typename TmpType,
typename LaneIdType>
EIGEN_STRONG_INLINE std::enable_if_t<!is_same<RhsPacketType, RhsPacketx4>::value> madd(const LhsPacketType& a,
const RhsPacketType& b,
DoublePacket<ResPacketType>& c,
TmpType& /*tmp*/,
const LaneIdType&) const {
c.first = padd(pmul(a, b.first), c.first);
c.second = padd(pmul(a, b.second), c.second);
}
template <typename LaneIdType>
EIGEN_STRONG_INLINE void madd(const LhsPacket& a, const RhsPacket& b, ResPacket& c, RhsPacket& /*tmp*/,
const LaneIdType&) const {
c = cj.pmadd(a, b, c);
}
template <typename LhsPacketType, typename AccPacketType, typename LaneIdType>
EIGEN_STRONG_INLINE void madd(const LhsPacketType& a, const RhsPacketx4& b, AccPacketType& c, RhsPacket& tmp,
const LaneIdType& lane) const {
madd(a, b.get(lane), c, tmp, lane);
}
EIGEN_STRONG_INLINE void acc(const Scalar& c, const Scalar& alpha, Scalar& r) const { r += alpha * c; }
template <typename RealPacketType, typename ResPacketType>
EIGEN_STRONG_INLINE void acc(const DoublePacket<RealPacketType>& c, const ResPacketType& alpha,
ResPacketType& r) const {
// assemble c
ResPacketType tmp;
if ((!ConjLhs) && (!ConjRhs)) {
tmp = pcplxflip(pconj(ResPacketType(c.second)));
tmp = padd(ResPacketType(c.first), tmp);
} else if ((!ConjLhs) && (ConjRhs)) {
tmp = pconj(pcplxflip(ResPacketType(c.second)));
tmp = padd(ResPacketType(c.first), tmp);
} else if ((ConjLhs) && (!ConjRhs)) {
tmp = pcplxflip(ResPacketType(c.second));
tmp = padd(pconj(ResPacketType(c.first)), tmp);
} else if ((ConjLhs) && (ConjRhs)) {
tmp = pcplxflip(ResPacketType(c.second));
tmp = psub(pconj(ResPacketType(c.first)), tmp);
}
r = pmadd(tmp, alpha, r);
}
protected:
conj_helper<LhsScalar, RhsScalar, ConjLhs, ConjRhs> cj;
};
template <typename RealScalar, bool ConjRhs_, int Arch, int PacketSize_>
class gebp_traits<RealScalar, std::complex<RealScalar>, false, ConjRhs_, Arch, PacketSize_> {
public:
typedef std::complex<RealScalar> Scalar;
typedef RealScalar LhsScalar;
typedef Scalar RhsScalar;
typedef Scalar ResScalar;
PACKET_DECL_COND_POSTFIX(_, Lhs, PacketSize_);
PACKET_DECL_COND_POSTFIX(_, Rhs, PacketSize_);
PACKET_DECL_COND_POSTFIX(_, Res, PacketSize_);
PACKET_DECL_COND_POSTFIX(_, Real, PacketSize_);
PACKET_DECL_COND_SCALAR_POSTFIX(_, PacketSize_);
#undef PACKET_DECL_COND_SCALAR_POSTFIX
#undef PACKET_DECL_COND_POSTFIX
#undef PACKET_DECL_COND_SCALAR
#undef PACKET_DECL_COND
enum {
ConjLhs = false,
ConjRhs = ConjRhs_,
Vectorizable = unpacket_traits<RealPacket_>::vectorizable && unpacket_traits<ScalarPacket_>::vectorizable,
LhsPacketSize = Vectorizable ? unpacket_traits<LhsPacket_>::size : 1,
RhsPacketSize = Vectorizable ? unpacket_traits<RhsPacket_>::size : 1,
ResPacketSize = Vectorizable ? unpacket_traits<ResPacket_>::size : 1,
NumberOfRegisters = EIGEN_ARCH_DEFAULT_NUMBER_OF_REGISTERS,
// FIXME: should depend on NumberOfRegisters
nr = 4,
mr = (plain_enum_min(16, NumberOfRegisters) / 2 / nr) * ResPacketSize,
LhsProgress = ResPacketSize,
RhsProgress = 1
};
typedef std::conditional_t<Vectorizable, LhsPacket_, LhsScalar> LhsPacket;
typedef std::conditional_t<Vectorizable, RhsPacket_, RhsScalar> RhsPacket;
typedef std::conditional_t<Vectorizable, ResPacket_, ResScalar> ResPacket;
typedef LhsPacket LhsPacket4Packing;
typedef QuadPacket<RhsPacket> RhsPacketx4;
typedef ResPacket AccPacket;
EIGEN_STRONG_INLINE void initAcc(AccPacket& p) { p = pset1<ResPacket>(ResScalar(0)); }
template <typename RhsPacketType>
EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacketType& dest) const {
dest = pset1<RhsPacketType>(*b);
}
EIGEN_STRONG_INLINE void loadRhs(const RhsScalar* b, RhsPacketx4& dest) const {
pbroadcast4(b, dest.B_0, dest.B1, dest.B2, dest.B3);
}
template <typename RhsPacketType>
EIGEN_STRONG_INLINE void updateRhs(const RhsScalar* b, RhsPacketType& dest) const {
loadRhs(b, dest);
}
EIGEN_STRONG_INLINE void updateRhs(const RhsScalar*, RhsPacketx4&) const {}
EIGEN_STRONG_INLINE void loadLhs(const LhsScalar* a, LhsPacket& dest) const { dest = ploaddup<LhsPacket>(a); }
EIGEN_STRONG_INLINE void loadRhsQuad(const RhsScalar* b, RhsPacket& dest) const { dest = ploadquad<RhsPacket>(b); }
template <typename LhsPacketType>
EIGEN_STRONG_INLINE void loadLhsUnaligned(const LhsScalar* a, LhsPacketType& dest) const {
dest = ploaddup<LhsPacketType>(a);
}
template <typename LhsPacketType, typename RhsPacketType, typename AccPacketType, typename LaneIdType>
EIGEN_STRONG_INLINE void madd(const LhsPacketType& a, const RhsPacketType& b, AccPacketType& c, RhsPacketType& tmp,
const LaneIdType&) const {
madd_impl(a, b, c, tmp, std::conditional_t<Vectorizable, true_type, false_type>());
}
template <typename LhsPacketType, typename RhsPacketType, typename AccPacketType>
EIGEN_STRONG_INLINE void madd_impl(const LhsPacketType& a, const RhsPacketType& b, AccPacketType& c,
RhsPacketType& tmp, const true_type&) const {
#ifdef EIGEN_HAS_SINGLE_INSTRUCTION_MADD
EIGEN_UNUSED_VARIABLE(tmp);
c.v = pmadd(a, b.v, c.v);
#else
tmp = b;
tmp.v = pmul(a, tmp.v);
c = padd(c, tmp);
#endif
}
EIGEN_STRONG_INLINE void madd_impl(const LhsScalar& a, const RhsScalar& b, ResScalar& c, RhsScalar& /*tmp*/,
const false_type&) const {
c += a * b;
}
template <typename LhsPacketType, typename AccPacketType, typename LaneIdType>
EIGEN_STRONG_INLINE void madd(const LhsPacketType& a, const RhsPacketx4& b, AccPacketType& c, RhsPacket& tmp,
const LaneIdType& lane) const {
madd(a, b.get(lane), c, tmp, lane);
}
template <typename ResPacketType, typename AccPacketType>
EIGEN_STRONG_INLINE void acc(const AccPacketType& c, const ResPacketType& alpha, ResPacketType& r) const {
conj_helper<ResPacketType, ResPacketType, false, ConjRhs> cj;
r = cj.pmadd(alpha, c, r);
}
protected:
};
/* optimized General packed Block * packed Panel product kernel
*
* Mixing type logic: C += A * B
* | A | B | comments
* |real |cplx | no vectorization yet, would require to pack A with duplication
* |cplx |real | easy vectorization
*/
template <typename LhsScalar, typename RhsScalar, typename Index, typename DataMapper, int mr, int nr,
bool ConjugateLhs, bool ConjugateRhs>
struct gebp_kernel {
typedef gebp_traits<LhsScalar, RhsScalar, ConjugateLhs, ConjugateRhs, Architecture::Target> Traits;
typedef gebp_traits<LhsScalar, RhsScalar, ConjugateLhs, ConjugateRhs, Architecture::Target, GEBPPacketHalf>
HalfTraits;
typedef gebp_traits<LhsScalar, RhsScalar, ConjugateLhs, ConjugateRhs, Architecture::Target, GEBPPacketQuarter>
QuarterTraits;
typedef typename Traits::ResScalar ResScalar;
typedef typename Traits::LhsPacket LhsPacket;
typedef typename Traits::RhsPacket RhsPacket;
typedef typename Traits::ResPacket ResPacket;
typedef typename Traits::AccPacket AccPacket;
typedef typename Traits::RhsPacketx4 RhsPacketx4;
typedef typename RhsPanelHelper<RhsPacket, RhsPacketx4, 15>::type RhsPanel15;
typedef typename RhsPanelHelper<RhsPacket, RhsPacketx4, 27>::type RhsPanel27;
typedef gebp_traits<RhsScalar, LhsScalar, ConjugateRhs, ConjugateLhs, Architecture::Target> SwappedTraits;
typedef typename SwappedTraits::ResScalar SResScalar;
typedef typename SwappedTraits::LhsPacket SLhsPacket;
typedef typename SwappedTraits::RhsPacket SRhsPacket;
typedef typename SwappedTraits::ResPacket SResPacket;
typedef typename SwappedTraits::AccPacket SAccPacket;
typedef typename HalfTraits::LhsPacket LhsPacketHalf;
typedef typename HalfTraits::RhsPacket RhsPacketHalf;
typedef typename HalfTraits::ResPacket ResPacketHalf;
typedef typename HalfTraits::AccPacket AccPacketHalf;
typedef typename QuarterTraits::LhsPacket LhsPacketQuarter;
typedef typename QuarterTraits::RhsPacket RhsPacketQuarter;
typedef typename QuarterTraits::ResPacket ResPacketQuarter;
typedef typename QuarterTraits::AccPacket AccPacketQuarter;
typedef typename DataMapper::LinearMapper LinearMapper;
enum {
Vectorizable = Traits::Vectorizable,
LhsProgress = Traits::LhsProgress,
LhsProgressHalf = HalfTraits::LhsProgress,
LhsProgressQuarter = QuarterTraits::LhsProgress,
RhsProgress = Traits::RhsProgress,
RhsProgressHalf = HalfTraits::RhsProgress,
RhsProgressQuarter = QuarterTraits::RhsProgress,
ResPacketSize = Traits::ResPacketSize
};
EIGEN_DONT_INLINE void operator()(const DataMapper& res, const LhsScalar* blockA, const RhsScalar* blockB, Index rows,
Index depth, Index cols, ResScalar alpha, Index strideA = -1, Index strideB = -1,
Index offsetA = 0, Index offsetB = 0);
};
template <typename LhsScalar, typename RhsScalar, typename Index, typename DataMapper, int mr, int nr,
bool ConjugateLhs, bool ConjugateRhs,
int SwappedLhsProgress =
gebp_traits<RhsScalar, LhsScalar, ConjugateRhs, ConjugateLhs, Architecture::Target>::LhsProgress>
struct last_row_process_16_packets {
typedef gebp_traits<LhsScalar, RhsScalar, ConjugateLhs, ConjugateRhs, Architecture::Target> Traits;
typedef gebp_traits<RhsScalar, LhsScalar, ConjugateRhs, ConjugateLhs, Architecture::Target> SwappedTraits;
typedef typename Traits::ResScalar ResScalar;
typedef typename SwappedTraits::LhsPacket SLhsPacket;
typedef typename SwappedTraits::RhsPacket SRhsPacket;
typedef typename SwappedTraits::ResPacket SResPacket;
typedef typename SwappedTraits::AccPacket SAccPacket;
EIGEN_STRONG_INLINE void operator()(const DataMapper& res, SwappedTraits& straits, const LhsScalar* blA,
const RhsScalar* blB, Index depth, const Index endk, Index i, Index j2,
ResScalar alpha, SAccPacket& C0) {
EIGEN_UNUSED_VARIABLE(res);
EIGEN_UNUSED_VARIABLE(straits);
EIGEN_UNUSED_VARIABLE(blA);
EIGEN_UNUSED_VARIABLE(blB);
EIGEN_UNUSED_VARIABLE(depth);
EIGEN_UNUSED_VARIABLE(endk);
EIGEN_UNUSED_VARIABLE(i);
EIGEN_UNUSED_VARIABLE(j2);
EIGEN_UNUSED_VARIABLE(alpha);
EIGEN_UNUSED_VARIABLE(C0);
}
};
template <typename LhsScalar, typename RhsScalar, typename Index, typename DataMapper, int mr, int nr,
bool ConjugateLhs, bool ConjugateRhs>
struct last_row_process_16_packets<LhsScalar, RhsScalar, Index, DataMapper, mr, nr, ConjugateLhs, ConjugateRhs, 16> {
typedef gebp_traits<LhsScalar, RhsScalar, ConjugateLhs, ConjugateRhs, Architecture::Target> Traits;
typedef gebp_traits<RhsScalar, LhsScalar, ConjugateRhs, ConjugateLhs, Architecture::Target> SwappedTraits;
typedef typename Traits::ResScalar ResScalar;
typedef typename SwappedTraits::LhsPacket SLhsPacket;
typedef typename SwappedTraits::RhsPacket SRhsPacket;
typedef typename SwappedTraits::ResPacket SResPacket;
typedef typename SwappedTraits::AccPacket SAccPacket;
EIGEN_STRONG_INLINE void operator()(const DataMapper& res, SwappedTraits& straits, const LhsScalar* blA,
const RhsScalar* blB, Index depth, const Index endk, Index i, Index j2,
ResScalar alpha, SAccPacket& C0) {
typedef typename unpacket_traits<typename unpacket_traits<SResPacket>::half>::half SResPacketQuarter;
typedef typename unpacket_traits<typename unpacket_traits<SLhsPacket>::half>::half SLhsPacketQuarter;
typedef typename unpacket_traits<typename unpacket_traits<SRhsPacket>::half>::half SRhsPacketQuarter;
typedef typename unpacket_traits<typename unpacket_traits<SAccPacket>::half>::half SAccPacketQuarter;
SResPacketQuarter R = res.template gatherPacket<SResPacketQuarter>(i, j2);
SResPacketQuarter alphav = pset1<SResPacketQuarter>(alpha);
if (depth - endk > 0) {
// We have to handle the last row(s) of the rhs, which
// correspond to a half-packet
SAccPacketQuarter c0 = predux_half_dowto4(predux_half_dowto4(C0));
for (Index kk = endk; kk < depth; kk++) {
SLhsPacketQuarter a0;
SRhsPacketQuarter b0;
straits.loadLhsUnaligned(blB, a0);
straits.loadRhs(blA, b0);
straits.madd(a0, b0, c0, b0, fix<0>);
blB += SwappedTraits::LhsProgress / 4;
blA += 1;
}
straits.acc(c0, alphav, R);
} else {
straits.acc(predux_half_dowto4(predux_half_dowto4(C0)), alphav, R);
}
res.scatterPacket(i, j2, R);
}
};
template <int nr, Index LhsProgress, Index RhsProgress, typename LhsScalar, typename RhsScalar, typename ResScalar,
typename AccPacket, typename LhsPacket, typename RhsPacket, typename ResPacket, typename GEBPTraits,
typename LinearMapper, typename DataMapper>
struct lhs_process_one_packet {
typedef typename GEBPTraits::RhsPacketx4 RhsPacketx4;
EIGEN_STRONG_INLINE void peeled_kc_onestep(Index K, const LhsScalar* blA, const RhsScalar* blB, GEBPTraits traits,
LhsPacket* A0, RhsPacketx4* rhs_panel, RhsPacket* T0, AccPacket* C0,
AccPacket* C1, AccPacket* C2, AccPacket* C3) {
EIGEN_ASM_COMMENT("begin step of gebp micro kernel 1X4");
EIGEN_ASM_COMMENT("Note: these asm comments work around bug 935!");
traits.loadLhs(&blA[(0 + 1 * K) * LhsProgress], *A0);
traits.loadRhs(&blB[(0 + 4 * K) * RhsProgress], *rhs_panel);
traits.madd(*A0, *rhs_panel, *C0, *T0, fix<0>);
traits.madd(*A0, *rhs_panel, *C1, *T0, fix<1>);
traits.madd(*A0, *rhs_panel, *C2, *T0, fix<2>);
traits.madd(*A0, *rhs_panel, *C3, *T0, fix<3>);
#if EIGEN_GNUC_STRICT_AT_LEAST(6, 0, 0) && defined(EIGEN_VECTORIZE_SSE) && !(EIGEN_COMP_LCC)
__asm__("" : "+x,m"(*A0));
#endif
EIGEN_ASM_COMMENT("end step of gebp micro kernel 1X4");
}
EIGEN_STRONG_INLINE void operator()(const DataMapper& res, const LhsScalar* blockA, const RhsScalar* blockB,
ResScalar alpha, Index peelStart, Index peelEnd, Index strideA, Index strideB,
Index offsetA, Index offsetB, int prefetch_res_offset, Index peeled_kc, Index pk,
Index cols, Index depth, Index packet_cols4) {
GEBPTraits traits;
Index packet_cols8 = nr >= 8 ? (cols / 8) * 8 : 0;
// loops on each largest micro horizontal panel of lhs
// (LhsProgress x depth)
for (Index i = peelStart; i < peelEnd; i += LhsProgress) {
#if EIGEN_ARCH_ARM64
EIGEN_IF_CONSTEXPR(nr >= 8) {
for (Index j2 = 0; j2 < packet_cols8; j2 += 8) {
const LhsScalar* blA = &blockA[i * strideA + offsetA * (LhsProgress)];
prefetch(&blA[0]);
// gets res block as register
AccPacket C0, C1, C2, C3, C4, C5, C6, C7;
traits.initAcc(C0);
traits.initAcc(C1);
traits.initAcc(C2);
traits.initAcc(C3);
traits.initAcc(C4);
traits.initAcc(C5);
traits.initAcc(C6);
traits.initAcc(C7);
LinearMapper r0 = res.getLinearMapper(i, j2 + 0);
LinearMapper r1 = res.getLinearMapper(i, j2 + 1);
LinearMapper r2 = res.getLinearMapper(i, j2 + 2);
LinearMapper r3 = res.getLinearMapper(i, j2 + 3);
LinearMapper r4 = res.getLinearMapper(i, j2 + 4);
LinearMapper r5 = res.getLinearMapper(i, j2 + 5);
LinearMapper r6 = res.getLinearMapper(i, j2 + 6);
LinearMapper r7 = res.getLinearMapper(i, j2 + 7);
r0.prefetch(prefetch_res_offset);
r1.prefetch(prefetch_res_offset);
r2.prefetch(prefetch_res_offset);
r3.prefetch(prefetch_res_offset);
r4.prefetch(prefetch_res_offset);
r5.prefetch(prefetch_res_offset);
r6.prefetch(prefetch_res_offset);
r7.prefetch(prefetch_res_offset);
const RhsScalar* blB = &blockB[j2 * strideB + offsetB * 8];
prefetch(&blB[0]);
LhsPacket A0;
for (Index k = 0; k < peeled_kc; k += pk) {
RhsPacketx4 rhs_panel;
RhsPacket T0;
#define EIGEN_GEBGP_ONESTEP(K) \
do { \
EIGEN_ASM_COMMENT("begin step of gebp micro kernel 1pX8"); \
traits.loadLhs(&blA[(0 + 1 * K) * LhsProgress], A0); \
traits.loadRhs(&blB[(0 + 8 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C0, T0, fix<0>); \
traits.updateRhs(&blB[(1 + 8 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C1, T0, fix<1>); \
traits.updateRhs(&blB[(2 + 8 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C2, T0, fix<2>); \
traits.updateRhs(&blB[(3 + 8 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C3, T0, fix<3>); \
traits.loadRhs(&blB[(4 + 8 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C4, T0, fix<0>); \
traits.updateRhs(&blB[(5 + 8 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C5, T0, fix<1>); \
traits.updateRhs(&blB[(6 + 8 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C6, T0, fix<2>); \
traits.updateRhs(&blB[(7 + 8 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C7, T0, fix<3>); \
EIGEN_ASM_COMMENT("end step of gebp micro kernel 1pX8"); \
} while (false)
EIGEN_ASM_COMMENT("begin gebp micro kernel 1pX8");
EIGEN_GEBGP_ONESTEP(0);
EIGEN_GEBGP_ONESTEP(1);
EIGEN_GEBGP_ONESTEP(2);
EIGEN_GEBGP_ONESTEP(3);
EIGEN_GEBGP_ONESTEP(4);
EIGEN_GEBGP_ONESTEP(5);
EIGEN_GEBGP_ONESTEP(6);
EIGEN_GEBGP_ONESTEP(7);
blB += pk * 8 * RhsProgress;
blA += pk * (1 * LhsProgress);
EIGEN_ASM_COMMENT("end gebp micro kernel 1pX8");
}
// process remaining peeled loop
for (Index k = peeled_kc; k < depth; k++) {
RhsPacketx4 rhs_panel;
RhsPacket T0;
EIGEN_GEBGP_ONESTEP(0);
blB += 8 * RhsProgress;
blA += 1 * LhsProgress;
}
#undef EIGEN_GEBGP_ONESTEP
ResPacket R0, R1;
ResPacket alphav = pset1<ResPacket>(alpha);
R0 = r0.template loadPacket<ResPacket>(0);
R1 = r1.template loadPacket<ResPacket>(0);
traits.acc(C0, alphav, R0);
traits.acc(C1, alphav, R1);
r0.storePacket(0, R0);
r1.storePacket(0, R1);
R0 = r2.template loadPacket<ResPacket>(0);
R1 = r3.template loadPacket<ResPacket>(0);
traits.acc(C2, alphav, R0);
traits.acc(C3, alphav, R1);
r2.storePacket(0, R0);
r3.storePacket(0, R1);
R0 = r4.template loadPacket<ResPacket>(0);
R1 = r5.template loadPacket<ResPacket>(0);
traits.acc(C4, alphav, R0);
traits.acc(C5, alphav, R1);
r4.storePacket(0, R0);
r5.storePacket(0, R1);
R0 = r6.template loadPacket<ResPacket>(0);
R1 = r7.template loadPacket<ResPacket>(0);
traits.acc(C6, alphav, R0);
traits.acc(C7, alphav, R1);
r6.storePacket(0, R0);
r7.storePacket(0, R1);
}
}
#endif
// loops on each largest micro vertical panel of rhs (depth * nr)
for (Index j2 = packet_cols8; j2 < packet_cols4; j2 += 4) {
// We select a LhsProgress x nr micro block of res
// which is entirely stored into 1 x nr registers.
const LhsScalar* blA = &blockA[i * strideA + offsetA * (LhsProgress)];
prefetch(&blA[0]);
// gets res block as register
AccPacket C0, C1, C2, C3;
traits.initAcc(C0);
traits.initAcc(C1);
traits.initAcc(C2);
traits.initAcc(C3);
// To improve instruction pipelining, let's double the accumulation registers:
// even k will accumulate in C*, while odd k will accumulate in D*.
// This trick is crutial to get good performance with FMA, otherwise it is
// actually faster to perform separated MUL+ADD because of a naturally
// better instruction-level parallelism.
AccPacket D0, D1, D2, D3;
traits.initAcc(D0);
traits.initAcc(D1);
traits.initAcc(D2);
traits.initAcc(D3);
LinearMapper r0 = res.getLinearMapper(i, j2 + 0);
LinearMapper r1 = res.getLinearMapper(i, j2 + 1);
LinearMapper r2 = res.getLinearMapper(i, j2 + 2);
LinearMapper r3 = res.getLinearMapper(i, j2 + 3);
r0.prefetch(prefetch_res_offset);
r1.prefetch(prefetch_res_offset);
r2.prefetch(prefetch_res_offset);
r3.prefetch(prefetch_res_offset);
// performs "inner" products
const RhsScalar* blB = &blockB[j2 * strideB + offsetB * 4];
prefetch(&blB[0]);
LhsPacket A0, A1;
for (Index k = 0; k < peeled_kc; k += pk) {
EIGEN_ASM_COMMENT("begin gebp micro kernel 1/half/quarterX4");
RhsPacketx4 rhs_panel;
RhsPacket T0;
internal::prefetch(blB + (48 + 0));
peeled_kc_onestep(0, blA, blB, traits, &A0, &rhs_panel, &T0, &C0, &C1, &C2, &C3);
peeled_kc_onestep(1, blA, blB, traits, &A1, &rhs_panel, &T0, &D0, &D1, &D2, &D3);
peeled_kc_onestep(2, blA, blB, traits, &A0, &rhs_panel, &T0, &C0, &C1, &C2, &C3);
peeled_kc_onestep(3, blA, blB, traits, &A1, &rhs_panel, &T0, &D0, &D1, &D2, &D3);
internal::prefetch(blB + (48 + 16));
peeled_kc_onestep(4, blA, blB, traits, &A0, &rhs_panel, &T0, &C0, &C1, &C2, &C3);
peeled_kc_onestep(5, blA, blB, traits, &A1, &rhs_panel, &T0, &D0, &D1, &D2, &D3);
peeled_kc_onestep(6, blA, blB, traits, &A0, &rhs_panel, &T0, &C0, &C1, &C2, &C3);
peeled_kc_onestep(7, blA, blB, traits, &A1, &rhs_panel, &T0, &D0, &D1, &D2, &D3);
blB += pk * 4 * RhsProgress;
blA += pk * LhsProgress;
EIGEN_ASM_COMMENT("end gebp micro kernel 1/half/quarterX4");
}
C0 = padd(C0, D0);
C1 = padd(C1, D1);
C2 = padd(C2, D2);
C3 = padd(C3, D3);
// process remaining peeled loop
for (Index k = peeled_kc; k < depth; k++) {
RhsPacketx4 rhs_panel;
RhsPacket T0;
peeled_kc_onestep(0, blA, blB, traits, &A0, &rhs_panel, &T0, &C0, &C1, &C2, &C3);
blB += 4 * RhsProgress;
blA += LhsProgress;
}
ResPacket R0, R1;
ResPacket alphav = pset1<ResPacket>(alpha);
R0 = r0.template loadPacket<ResPacket>(0);
R1 = r1.template loadPacket<ResPacket>(0);
traits.acc(C0, alphav, R0);
traits.acc(C1, alphav, R1);
r0.storePacket(0, R0);
r1.storePacket(0, R1);
R0 = r2.template loadPacket<ResPacket>(0);
R1 = r3.template loadPacket<ResPacket>(0);
traits.acc(C2, alphav, R0);
traits.acc(C3, alphav, R1);
r2.storePacket(0, R0);
r3.storePacket(0, R1);
}
// Deal with remaining columns of the rhs
for (Index j2 = packet_cols4; j2 < cols; j2++) {
// One column at a time
const LhsScalar* blA = &blockA[i * strideA + offsetA * (LhsProgress)];
prefetch(&blA[0]);
// gets res block as register
AccPacket C0;
traits.initAcc(C0);
LinearMapper r0 = res.getLinearMapper(i, j2);
// performs "inner" products
const RhsScalar* blB = &blockB[j2 * strideB + offsetB];
LhsPacket A0;
for (Index k = 0; k < peeled_kc; k += pk) {
EIGEN_ASM_COMMENT("begin gebp micro kernel 1/half/quarterX1");
RhsPacket B_0;
#define EIGEN_GEBGP_ONESTEP(K) \
do { \
EIGEN_ASM_COMMENT("begin step of gebp micro kernel 1/half/quarterX1"); \
EIGEN_ASM_COMMENT("Note: these asm comments work around bug 935!"); \
/* FIXME: why unaligned???? */ \
traits.loadLhsUnaligned(&blA[(0 + 1 * K) * LhsProgress], A0); \
traits.loadRhs(&blB[(0 + K) * RhsProgress], B_0); \
traits.madd(A0, B_0, C0, B_0, fix<0>); \
EIGEN_ASM_COMMENT("end step of gebp micro kernel 1/half/quarterX1"); \
} while (false);
EIGEN_GEBGP_ONESTEP(0);
EIGEN_GEBGP_ONESTEP(1);
EIGEN_GEBGP_ONESTEP(2);
EIGEN_GEBGP_ONESTEP(3);
EIGEN_GEBGP_ONESTEP(4);
EIGEN_GEBGP_ONESTEP(5);
EIGEN_GEBGP_ONESTEP(6);
EIGEN_GEBGP_ONESTEP(7);
blB += pk * RhsProgress;
blA += pk * LhsProgress;
EIGEN_ASM_COMMENT("end gebp micro kernel 1/half/quarterX1");
}
// process remaining peeled loop
for (Index k = peeled_kc; k < depth; k++) {
RhsPacket B_0;
EIGEN_GEBGP_ONESTEP(0);
blB += RhsProgress;
blA += LhsProgress;
}
#undef EIGEN_GEBGP_ONESTEP
ResPacket R0;
ResPacket alphav = pset1<ResPacket>(alpha);
R0 = r0.template loadPacket<ResPacket>(0);
traits.acc(C0, alphav, R0);
r0.storePacket(0, R0);
}
}
}
};
template <int nr, Index LhsProgress, Index RhsProgress, typename LhsScalar, typename RhsScalar, typename ResScalar,
typename AccPacket, typename LhsPacket, typename RhsPacket, typename ResPacket, typename GEBPTraits,
typename LinearMapper, typename DataMapper>
struct lhs_process_fraction_of_packet
: lhs_process_one_packet<nr, LhsProgress, RhsProgress, LhsScalar, RhsScalar, ResScalar, AccPacket, LhsPacket,
RhsPacket, ResPacket, GEBPTraits, LinearMapper, DataMapper> {
EIGEN_STRONG_INLINE void peeled_kc_onestep(Index K, const LhsScalar* blA, const RhsScalar* blB, GEBPTraits traits,
LhsPacket* A0, RhsPacket* B_0, RhsPacket* B1, RhsPacket* B2, RhsPacket* B3,
AccPacket* C0, AccPacket* C1, AccPacket* C2, AccPacket* C3) {
EIGEN_ASM_COMMENT("begin step of gebp micro kernel 1X4");
EIGEN_ASM_COMMENT("Note: these asm comments work around bug 935!");
traits.loadLhsUnaligned(&blA[(0 + 1 * K) * (LhsProgress)], *A0);
traits.broadcastRhs(&blB[(0 + 4 * K) * RhsProgress], *B_0, *B1, *B2, *B3);
traits.madd(*A0, *B_0, *C0, *B_0);
traits.madd(*A0, *B1, *C1, *B1);
traits.madd(*A0, *B2, *C2, *B2);
traits.madd(*A0, *B3, *C3, *B3);
EIGEN_ASM_COMMENT("end step of gebp micro kernel 1X4");
}
};
template <typename LhsScalar, typename RhsScalar, typename Index, typename DataMapper, int mr, int nr,
bool ConjugateLhs, bool ConjugateRhs>
EIGEN_DONT_INLINE void gebp_kernel<LhsScalar, RhsScalar, Index, DataMapper, mr, nr, ConjugateLhs,
ConjugateRhs>::operator()(const DataMapper& res, const LhsScalar* blockA,
const RhsScalar* blockB, Index rows, Index depth,
Index cols, ResScalar alpha, Index strideA, Index strideB,
Index offsetA, Index offsetB) {
Traits traits;
SwappedTraits straits;
if (strideA == -1) strideA = depth;
if (strideB == -1) strideB = depth;
conj_helper<LhsScalar, RhsScalar, ConjugateLhs, ConjugateRhs> cj;
Index packet_cols4 = nr >= 4 ? (cols / 4) * 4 : 0;
Index packet_cols8 = nr >= 8 ? (cols / 8) * 8 : 0;
const Index peeled_mc3 = mr >= 3 * Traits::LhsProgress ? (rows / (3 * LhsProgress)) * (3 * LhsProgress) : 0;
const Index peeled_mc2 =
mr >= 2 * Traits::LhsProgress ? peeled_mc3 + ((rows - peeled_mc3) / (2 * LhsProgress)) * (2 * LhsProgress) : 0;
const Index peeled_mc1 =
mr >= 1 * Traits::LhsProgress ? peeled_mc2 + ((rows - peeled_mc2) / (1 * LhsProgress)) * (1 * LhsProgress) : 0;
const Index peeled_mc_half =
mr >= LhsProgressHalf ? peeled_mc1 + ((rows - peeled_mc1) / (LhsProgressHalf)) * (LhsProgressHalf) : 0;
const Index peeled_mc_quarter =
mr >= LhsProgressQuarter
? peeled_mc_half + ((rows - peeled_mc_half) / (LhsProgressQuarter)) * (LhsProgressQuarter)
: 0;
enum { pk = 8 }; // NOTE Such a large peeling factor is important for large matrices (~ +5% when >1000 on Haswell)
const Index peeled_kc = depth & ~(pk - 1);
const int prefetch_res_offset = 32 / sizeof(ResScalar);
// const Index depth2 = depth & ~1;
//---------- Process 3 * LhsProgress rows at once ----------
// This corresponds to 3*LhsProgress x nr register blocks.
// Usually, make sense only with FMA
if (mr >= 3 * Traits::LhsProgress) {
// Here, the general idea is to loop on each largest micro horizontal panel of the lhs (3*Traits::LhsProgress x
// depth) and on each largest micro vertical panel of the rhs (depth * nr). Blocking sizes, i.e., 'depth' has been
// computed so that the micro horizontal panel of the lhs fit in L1. However, if depth is too small, we can extend
// the number of rows of these horizontal panels. This actual number of rows is computed as follow:
const Index l1 = defaultL1CacheSize; // in Bytes, TODO, l1 should be passed to this function.
// The max(1, ...) here is needed because we may be using blocking params larger than what our known l1 cache size
// suggests we should be using: either because our known l1 cache size is inaccurate (e.g. on Android, we can only
// guess), or because we are testing specific blocking sizes.
const Index actual_panel_rows =
(3 * LhsProgress) * std::max<Index>(1, ((l1 - sizeof(ResScalar) * mr * nr - depth * nr * sizeof(RhsScalar)) /
(depth * sizeof(LhsScalar) * 3 * LhsProgress)));
for (Index i1 = 0; i1 < peeled_mc3; i1 += actual_panel_rows) {
const Index actual_panel_end = (std::min)(i1 + actual_panel_rows, peeled_mc3);
#if EIGEN_ARCH_ARM64
EIGEN_IF_CONSTEXPR(nr >= 8) {
for (Index j2 = 0; j2 < packet_cols8; j2 += 8) {
for (Index i = i1; i < actual_panel_end; i += 3 * LhsProgress) {
const LhsScalar* blA = &blockA[i * strideA + offsetA * (3 * LhsProgress)];
prefetch(&blA[0]);
// gets res block as register
AccPacket C0, C1, C2, C3, C4, C5, C6, C7, C8, C9, C10, C11, C12, C13, C14, C15, C16, C17, C18, C19, C20,
C21, C22, C23;
traits.initAcc(C0);
traits.initAcc(C1);
traits.initAcc(C2);
traits.initAcc(C3);
traits.initAcc(C4);
traits.initAcc(C5);
traits.initAcc(C6);
traits.initAcc(C7);
traits.initAcc(C8);
traits.initAcc(C9);
traits.initAcc(C10);
traits.initAcc(C11);
traits.initAcc(C12);
traits.initAcc(C13);
traits.initAcc(C14);
traits.initAcc(C15);
traits.initAcc(C16);
traits.initAcc(C17);
traits.initAcc(C18);
traits.initAcc(C19);
traits.initAcc(C20);
traits.initAcc(C21);
traits.initAcc(C22);
traits.initAcc(C23);
LinearMapper r0 = res.getLinearMapper(i, j2 + 0);
LinearMapper r1 = res.getLinearMapper(i, j2 + 1);
LinearMapper r2 = res.getLinearMapper(i, j2 + 2);
LinearMapper r3 = res.getLinearMapper(i, j2 + 3);
LinearMapper r4 = res.getLinearMapper(i, j2 + 4);
LinearMapper r5 = res.getLinearMapper(i, j2 + 5);
LinearMapper r6 = res.getLinearMapper(i, j2 + 6);
LinearMapper r7 = res.getLinearMapper(i, j2 + 7);
r0.prefetch(0);
r1.prefetch(0);
r2.prefetch(0);
r3.prefetch(0);
r4.prefetch(0);
r5.prefetch(0);
r6.prefetch(0);
r7.prefetch(0);
// performs "inner" products
const RhsScalar* blB = &blockB[j2 * strideB + offsetB * 8];
prefetch(&blB[0]);
LhsPacket A0, A1;
for (Index k = 0; k < peeled_kc; k += pk) {
EIGEN_ASM_COMMENT("begin gebp micro kernel 3pX8");
// 27 registers are taken (24 for acc, 3 for lhs).
RhsPanel27 rhs_panel;
RhsPacket T0;
LhsPacket A2;
#if EIGEN_ARCH_ARM64 && defined(EIGEN_VECTORIZE_NEON) && EIGEN_GNUC_STRICT_LESS_THAN(9, 0, 0)
// see http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1633
// without this workaround A0, A1, and A2 are loaded in the same register,
// which is not good for pipelining
#define EIGEN_GEBP_3Px8_REGISTER_ALLOC_WORKAROUND __asm__("" : "+w,m"(A0), "+w,m"(A1), "+w,m"(A2));
#else
#define EIGEN_GEBP_3Px8_REGISTER_ALLOC_WORKAROUND
#endif
#define EIGEN_GEBP_ONESTEP(K) \
do { \
EIGEN_ASM_COMMENT("begin step of gebp micro kernel 3pX8"); \
traits.loadLhs(&blA[(0 + 3 * K) * LhsProgress], A0); \
traits.loadLhs(&blA[(1 + 3 * K) * LhsProgress], A1); \
traits.loadLhs(&blA[(2 + 3 * K) * LhsProgress], A2); \
EIGEN_GEBP_3Px8_REGISTER_ALLOC_WORKAROUND traits.loadRhs(blB + (0 + 8 * K) * Traits::RhsProgress, rhs_panel); \
traits.madd(A0, rhs_panel, C0, T0, fix<0>); \
traits.madd(A1, rhs_panel, C8, T0, fix<0>); \
traits.madd(A2, rhs_panel, C16, T0, fix<0>); \
traits.updateRhs(blB + (1 + 8 * K) * Traits::RhsProgress, rhs_panel); \
traits.madd(A0, rhs_panel, C1, T0, fix<1>); \
traits.madd(A1, rhs_panel, C9, T0, fix<1>); \
traits.madd(A2, rhs_panel, C17, T0, fix<1>); \
traits.updateRhs(blB + (2 + 8 * K) * Traits::RhsProgress, rhs_panel); \
traits.madd(A0, rhs_panel, C2, T0, fix<2>); \
traits.madd(A1, rhs_panel, C10, T0, fix<2>); \
traits.madd(A2, rhs_panel, C18, T0, fix<2>); \
traits.updateRhs(blB + (3 + 8 * K) * Traits::RhsProgress, rhs_panel); \
traits.madd(A0, rhs_panel, C3, T0, fix<3>); \
traits.madd(A1, rhs_panel, C11, T0, fix<3>); \
traits.madd(A2, rhs_panel, C19, T0, fix<3>); \
traits.loadRhs(blB + (4 + 8 * K) * Traits::RhsProgress, rhs_panel); \
traits.madd(A0, rhs_panel, C4, T0, fix<0>); \
traits.madd(A1, rhs_panel, C12, T0, fix<0>); \
traits.madd(A2, rhs_panel, C20, T0, fix<0>); \
traits.updateRhs(blB + (5 + 8 * K) * Traits::RhsProgress, rhs_panel); \
traits.madd(A0, rhs_panel, C5, T0, fix<1>); \
traits.madd(A1, rhs_panel, C13, T0, fix<1>); \
traits.madd(A2, rhs_panel, C21, T0, fix<1>); \
traits.updateRhs(blB + (6 + 8 * K) * Traits::RhsProgress, rhs_panel); \
traits.madd(A0, rhs_panel, C6, T0, fix<2>); \
traits.madd(A1, rhs_panel, C14, T0, fix<2>); \
traits.madd(A2, rhs_panel, C22, T0, fix<2>); \
traits.updateRhs(blB + (7 + 8 * K) * Traits::RhsProgress, rhs_panel); \
traits.madd(A0, rhs_panel, C7, T0, fix<3>); \
traits.madd(A1, rhs_panel, C15, T0, fix<3>); \
traits.madd(A2, rhs_panel, C23, T0, fix<3>); \
EIGEN_ASM_COMMENT("end step of gebp micro kernel 3pX8"); \
} while (false)
EIGEN_GEBP_ONESTEP(0);
EIGEN_GEBP_ONESTEP(1);
EIGEN_GEBP_ONESTEP(2);
EIGEN_GEBP_ONESTEP(3);
EIGEN_GEBP_ONESTEP(4);
EIGEN_GEBP_ONESTEP(5);
EIGEN_GEBP_ONESTEP(6);
EIGEN_GEBP_ONESTEP(7);
blB += pk * 8 * RhsProgress;
blA += pk * 3 * Traits::LhsProgress;
EIGEN_ASM_COMMENT("end gebp micro kernel 3pX8");
}
// process remaining peeled loop
for (Index k = peeled_kc; k < depth; k++) {
RhsPanel27 rhs_panel;
RhsPacket T0;
LhsPacket A2;
EIGEN_GEBP_ONESTEP(0);
blB += 8 * RhsProgress;
blA += 3 * Traits::LhsProgress;
}
#undef EIGEN_GEBP_ONESTEP
ResPacket R0, R1, R2;
ResPacket alphav = pset1<ResPacket>(alpha);
R0 = r0.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r0.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r0.template loadPacket<ResPacket>(2 * Traits::ResPacketSize);
traits.acc(C0, alphav, R0);
traits.acc(C8, alphav, R1);
traits.acc(C16, alphav, R2);
r0.storePacket(0 * Traits::ResPacketSize, R0);
r0.storePacket(1 * Traits::ResPacketSize, R1);
r0.storePacket(2 * Traits::ResPacketSize, R2);
R0 = r1.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r1.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r1.template loadPacket<ResPacket>(2 * Traits::ResPacketSize);
traits.acc(C1, alphav, R0);
traits.acc(C9, alphav, R1);
traits.acc(C17, alphav, R2);
r1.storePacket(0 * Traits::ResPacketSize, R0);
r1.storePacket(1 * Traits::ResPacketSize, R1);
r1.storePacket(2 * Traits::ResPacketSize, R2);
R0 = r2.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r2.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r2.template loadPacket<ResPacket>(2 * Traits::ResPacketSize);
traits.acc(C2, alphav, R0);
traits.acc(C10, alphav, R1);
traits.acc(C18, alphav, R2);
r2.storePacket(0 * Traits::ResPacketSize, R0);
r2.storePacket(1 * Traits::ResPacketSize, R1);
r2.storePacket(2 * Traits::ResPacketSize, R2);
R0 = r3.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r3.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r3.template loadPacket<ResPacket>(2 * Traits::ResPacketSize);
traits.acc(C3, alphav, R0);
traits.acc(C11, alphav, R1);
traits.acc(C19, alphav, R2);
r3.storePacket(0 * Traits::ResPacketSize, R0);
r3.storePacket(1 * Traits::ResPacketSize, R1);
r3.storePacket(2 * Traits::ResPacketSize, R2);
R0 = r4.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r4.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r4.template loadPacket<ResPacket>(2 * Traits::ResPacketSize);
traits.acc(C4, alphav, R0);
traits.acc(C12, alphav, R1);
traits.acc(C20, alphav, R2);
r4.storePacket(0 * Traits::ResPacketSize, R0);
r4.storePacket(1 * Traits::ResPacketSize, R1);
r4.storePacket(2 * Traits::ResPacketSize, R2);
R0 = r5.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r5.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r5.template loadPacket<ResPacket>(2 * Traits::ResPacketSize);
traits.acc(C5, alphav, R0);
traits.acc(C13, alphav, R1);
traits.acc(C21, alphav, R2);
r5.storePacket(0 * Traits::ResPacketSize, R0);
r5.storePacket(1 * Traits::ResPacketSize, R1);
r5.storePacket(2 * Traits::ResPacketSize, R2);
R0 = r6.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r6.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r6.template loadPacket<ResPacket>(2 * Traits::ResPacketSize);
traits.acc(C6, alphav, R0);
traits.acc(C14, alphav, R1);
traits.acc(C22, alphav, R2);
r6.storePacket(0 * Traits::ResPacketSize, R0);
r6.storePacket(1 * Traits::ResPacketSize, R1);
r6.storePacket(2 * Traits::ResPacketSize, R2);
R0 = r7.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r7.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r7.template loadPacket<ResPacket>(2 * Traits::ResPacketSize);
traits.acc(C7, alphav, R0);
traits.acc(C15, alphav, R1);
traits.acc(C23, alphav, R2);
r7.storePacket(0 * Traits::ResPacketSize, R0);
r7.storePacket(1 * Traits::ResPacketSize, R1);
r7.storePacket(2 * Traits::ResPacketSize, R2);
}
}
}
#endif
for (Index j2 = packet_cols8; j2 < packet_cols4; j2 += 4) {
for (Index i = i1; i < actual_panel_end; i += 3 * LhsProgress) {
// We selected a 3*Traits::LhsProgress x nr micro block of res which is entirely
// stored into 3 x nr registers.
const LhsScalar* blA = &blockA[i * strideA + offsetA * (3 * LhsProgress)];
prefetch(&blA[0]);
// gets res block as register
AccPacket C0, C1, C2, C3, C4, C5, C6, C7, C8, C9, C10, C11;
traits.initAcc(C0);
traits.initAcc(C1);
traits.initAcc(C2);
traits.initAcc(C3);
traits.initAcc(C4);
traits.initAcc(C5);
traits.initAcc(C6);
traits.initAcc(C7);
traits.initAcc(C8);
traits.initAcc(C9);
traits.initAcc(C10);
traits.initAcc(C11);
LinearMapper r0 = res.getLinearMapper(i, j2 + 0);
LinearMapper r1 = res.getLinearMapper(i, j2 + 1);
LinearMapper r2 = res.getLinearMapper(i, j2 + 2);
LinearMapper r3 = res.getLinearMapper(i, j2 + 3);
r0.prefetch(0);
r1.prefetch(0);
r2.prefetch(0);
r3.prefetch(0);
// performs "inner" products
const RhsScalar* blB = &blockB[j2 * strideB + offsetB * 4];
prefetch(&blB[0]);
LhsPacket A0, A1;
for (Index k = 0; k < peeled_kc; k += pk) {
EIGEN_ASM_COMMENT("begin gebp micro kernel 3pX4");
// 15 registers are taken (12 for acc, 3 for lhs).
RhsPanel15 rhs_panel;
RhsPacket T0;
LhsPacket A2;
#if EIGEN_ARCH_ARM64 && defined(EIGEN_VECTORIZE_NEON) && EIGEN_GNUC_STRICT_LESS_THAN(9, 0, 0)
// see http://eigen.tuxfamily.org/bz/show_bug.cgi?id=1633
// without this workaround A0, A1, and A2 are loaded in the same register,
// which is not good for pipelining
#define EIGEN_GEBP_3PX4_REGISTER_ALLOC_WORKAROUND __asm__("" : "+w,m"(A0), "+w,m"(A1), "+w,m"(A2));
#else
#define EIGEN_GEBP_3PX4_REGISTER_ALLOC_WORKAROUND
#endif
#define EIGEN_GEBP_ONESTEP(K) \
do { \
EIGEN_ASM_COMMENT("begin step of gebp micro kernel 3pX4"); \
EIGEN_ASM_COMMENT("Note: these asm comments work around bug 935!"); \
internal::prefetch(blA + (3 * K + 16) * LhsProgress); \
if (EIGEN_ARCH_ARM || EIGEN_ARCH_MIPS) { \
internal::prefetch(blB + (4 * K + 16) * RhsProgress); \
} /* Bug 953 */ \
traits.loadLhs(&blA[(0 + 3 * K) * LhsProgress], A0); \
traits.loadLhs(&blA[(1 + 3 * K) * LhsProgress], A1); \
traits.loadLhs(&blA[(2 + 3 * K) * LhsProgress], A2); \
EIGEN_GEBP_3PX4_REGISTER_ALLOC_WORKAROUND \
traits.loadRhs(blB + (0 + 4 * K) * Traits::RhsProgress, rhs_panel); \
traits.madd(A0, rhs_panel, C0, T0, fix<0>); \
traits.madd(A1, rhs_panel, C4, T0, fix<0>); \
traits.madd(A2, rhs_panel, C8, T0, fix<0>); \
traits.updateRhs(blB + (1 + 4 * K) * Traits::RhsProgress, rhs_panel); \
traits.madd(A0, rhs_panel, C1, T0, fix<1>); \
traits.madd(A1, rhs_panel, C5, T0, fix<1>); \
traits.madd(A2, rhs_panel, C9, T0, fix<1>); \
traits.updateRhs(blB + (2 + 4 * K) * Traits::RhsProgress, rhs_panel); \
traits.madd(A0, rhs_panel, C2, T0, fix<2>); \
traits.madd(A1, rhs_panel, C6, T0, fix<2>); \
traits.madd(A2, rhs_panel, C10, T0, fix<2>); \
traits.updateRhs(blB + (3 + 4 * K) * Traits::RhsProgress, rhs_panel); \
traits.madd(A0, rhs_panel, C3, T0, fix<3>); \
traits.madd(A1, rhs_panel, C7, T0, fix<3>); \
traits.madd(A2, rhs_panel, C11, T0, fix<3>); \
EIGEN_ASM_COMMENT("end step of gebp micro kernel 3pX4"); \
} while (false)
internal::prefetch(blB);
EIGEN_GEBP_ONESTEP(0);
EIGEN_GEBP_ONESTEP(1);
EIGEN_GEBP_ONESTEP(2);
EIGEN_GEBP_ONESTEP(3);
EIGEN_GEBP_ONESTEP(4);
EIGEN_GEBP_ONESTEP(5);
EIGEN_GEBP_ONESTEP(6);
EIGEN_GEBP_ONESTEP(7);
blB += pk * 4 * RhsProgress;
blA += pk * 3 * Traits::LhsProgress;
EIGEN_ASM_COMMENT("end gebp micro kernel 3pX4");
}
// process remaining peeled loop
for (Index k = peeled_kc; k < depth; k++) {
RhsPanel15 rhs_panel;
RhsPacket T0;
LhsPacket A2;
EIGEN_GEBP_ONESTEP(0);
blB += 4 * RhsProgress;
blA += 3 * Traits::LhsProgress;
}
#undef EIGEN_GEBP_ONESTEP
ResPacket R0, R1, R2;
ResPacket alphav = pset1<ResPacket>(alpha);
R0 = r0.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r0.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r0.template loadPacket<ResPacket>(2 * Traits::ResPacketSize);
traits.acc(C0, alphav, R0);
traits.acc(C4, alphav, R1);
traits.acc(C8, alphav, R2);
r0.storePacket(0 * Traits::ResPacketSize, R0);
r0.storePacket(1 * Traits::ResPacketSize, R1);
r0.storePacket(2 * Traits::ResPacketSize, R2);
R0 = r1.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r1.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r1.template loadPacket<ResPacket>(2 * Traits::ResPacketSize);
traits.acc(C1, alphav, R0);
traits.acc(C5, alphav, R1);
traits.acc(C9, alphav, R2);
r1.storePacket(0 * Traits::ResPacketSize, R0);
r1.storePacket(1 * Traits::ResPacketSize, R1);
r1.storePacket(2 * Traits::ResPacketSize, R2);
R0 = r2.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r2.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r2.template loadPacket<ResPacket>(2 * Traits::ResPacketSize);
traits.acc(C2, alphav, R0);
traits.acc(C6, alphav, R1);
traits.acc(C10, alphav, R2);
r2.storePacket(0 * Traits::ResPacketSize, R0);
r2.storePacket(1 * Traits::ResPacketSize, R1);
r2.storePacket(2 * Traits::ResPacketSize, R2);
R0 = r3.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r3.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r3.template loadPacket<ResPacket>(2 * Traits::ResPacketSize);
traits.acc(C3, alphav, R0);
traits.acc(C7, alphav, R1);
traits.acc(C11, alphav, R2);
r3.storePacket(0 * Traits::ResPacketSize, R0);
r3.storePacket(1 * Traits::ResPacketSize, R1);
r3.storePacket(2 * Traits::ResPacketSize, R2);
}
}
// Deal with remaining columns of the rhs
for (Index j2 = packet_cols4; j2 < cols; j2++) {
for (Index i = i1; i < actual_panel_end; i += 3 * LhsProgress) {
// One column at a time
const LhsScalar* blA = &blockA[i * strideA + offsetA * (3 * Traits::LhsProgress)];
prefetch(&blA[0]);
// gets res block as register
AccPacket C0, C4, C8;
traits.initAcc(C0);
traits.initAcc(C4);
traits.initAcc(C8);
LinearMapper r0 = res.getLinearMapper(i, j2);
r0.prefetch(0);
// performs "inner" products
const RhsScalar* blB = &blockB[j2 * strideB + offsetB];
LhsPacket A0, A1, A2;
for (Index k = 0; k < peeled_kc; k += pk) {
EIGEN_ASM_COMMENT("begin gebp micro kernel 3pX1");
RhsPacket B_0;
#define EIGEN_GEBGP_ONESTEP(K) \
do { \
EIGEN_ASM_COMMENT("begin step of gebp micro kernel 3pX1"); \
EIGEN_ASM_COMMENT("Note: these asm comments work around bug 935!"); \
traits.loadLhs(&blA[(0 + 3 * K) * LhsProgress], A0); \
traits.loadLhs(&blA[(1 + 3 * K) * LhsProgress], A1); \
traits.loadLhs(&blA[(2 + 3 * K) * LhsProgress], A2); \
traits.loadRhs(&blB[(0 + K) * RhsProgress], B_0); \
traits.madd(A0, B_0, C0, B_0, fix<0>); \
traits.madd(A1, B_0, C4, B_0, fix<0>); \
traits.madd(A2, B_0, C8, B_0, fix<0>); \
EIGEN_ASM_COMMENT("end step of gebp micro kernel 3pX1"); \
} while (false)
EIGEN_GEBGP_ONESTEP(0);
EIGEN_GEBGP_ONESTEP(1);
EIGEN_GEBGP_ONESTEP(2);
EIGEN_GEBGP_ONESTEP(3);
EIGEN_GEBGP_ONESTEP(4);
EIGEN_GEBGP_ONESTEP(5);
EIGEN_GEBGP_ONESTEP(6);
EIGEN_GEBGP_ONESTEP(7);
blB += int(pk) * int(RhsProgress);
blA += int(pk) * 3 * int(Traits::LhsProgress);
EIGEN_ASM_COMMENT("end gebp micro kernel 3pX1");
}
// process remaining peeled loop
for (Index k = peeled_kc; k < depth; k++) {
RhsPacket B_0;
EIGEN_GEBGP_ONESTEP(0);
blB += RhsProgress;
blA += 3 * Traits::LhsProgress;
}
#undef EIGEN_GEBGP_ONESTEP
ResPacket R0, R1, R2;
ResPacket alphav = pset1<ResPacket>(alpha);
R0 = r0.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r0.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r0.template loadPacket<ResPacket>(2 * Traits::ResPacketSize);
traits.acc(C0, alphav, R0);
traits.acc(C4, alphav, R1);
traits.acc(C8, alphav, R2);
r0.storePacket(0 * Traits::ResPacketSize, R0);
r0.storePacket(1 * Traits::ResPacketSize, R1);
r0.storePacket(2 * Traits::ResPacketSize, R2);
}
}
}
}
//---------- Process 2 * LhsProgress rows at once ----------
if (mr >= 2 * Traits::LhsProgress) {
const Index l1 = defaultL1CacheSize; // in Bytes, TODO, l1 should be passed to this function.
// The max(1, ...) here is needed because we may be using blocking params larger than what our known l1 cache size
// suggests we should be using: either because our known l1 cache size is inaccurate (e.g. on Android, we can only
// guess), or because we are testing specific blocking sizes.
Index actual_panel_rows =
(2 * LhsProgress) * std::max<Index>(1, ((l1 - sizeof(ResScalar) * mr * nr - depth * nr * sizeof(RhsScalar)) /
(depth * sizeof(LhsScalar) * 2 * LhsProgress)));
for (Index i1 = peeled_mc3; i1 < peeled_mc2; i1 += actual_panel_rows) {
Index actual_panel_end = (std::min)(i1 + actual_panel_rows, peeled_mc2);
#if EIGEN_ARCH_ARM64
EIGEN_IF_CONSTEXPR(nr >= 8) {
for (Index j2 = 0; j2 < packet_cols8; j2 += 8) {
for (Index i = i1; i < actual_panel_end; i += 2 * LhsProgress) {
const LhsScalar* blA = &blockA[i * strideA + offsetA * (2 * Traits::LhsProgress)];
prefetch(&blA[0]);
AccPacket C0, C1, C2, C3, C4, C5, C6, C7, C8, C9, C10, C11, C12, C13, C14, C15;
traits.initAcc(C0);
traits.initAcc(C1);
traits.initAcc(C2);
traits.initAcc(C3);
traits.initAcc(C4);
traits.initAcc(C5);
traits.initAcc(C6);
traits.initAcc(C7);
traits.initAcc(C8);
traits.initAcc(C9);
traits.initAcc(C10);
traits.initAcc(C11);
traits.initAcc(C12);
traits.initAcc(C13);
traits.initAcc(C14);
traits.initAcc(C15);
LinearMapper r0 = res.getLinearMapper(i, j2 + 0);
LinearMapper r1 = res.getLinearMapper(i, j2 + 1);
LinearMapper r2 = res.getLinearMapper(i, j2 + 2);
LinearMapper r3 = res.getLinearMapper(i, j2 + 3);
LinearMapper r4 = res.getLinearMapper(i, j2 + 4);
LinearMapper r5 = res.getLinearMapper(i, j2 + 5);
LinearMapper r6 = res.getLinearMapper(i, j2 + 6);
LinearMapper r7 = res.getLinearMapper(i, j2 + 7);
r0.prefetch(prefetch_res_offset);
r1.prefetch(prefetch_res_offset);
r2.prefetch(prefetch_res_offset);
r3.prefetch(prefetch_res_offset);
r4.prefetch(prefetch_res_offset);
r5.prefetch(prefetch_res_offset);
r6.prefetch(prefetch_res_offset);
r7.prefetch(prefetch_res_offset);
const RhsScalar* blB = &blockB[j2 * strideB + offsetB * 8];
prefetch(&blB[0]);
LhsPacket A0, A1;
for (Index k = 0; k < peeled_kc; k += pk) {
RhsPacketx4 rhs_panel;
RhsPacket T0;
// NOTE: the begin/end asm comments below work around bug 935!
// but they are not enough for gcc>=6 without FMA (bug 1637)
#if EIGEN_GNUC_STRICT_AT_LEAST(6, 0, 0) && defined(EIGEN_VECTORIZE_SSE)
#define EIGEN_GEBP_2Px8_SPILLING_WORKAROUND __asm__("" : [a0] "+x,m"(A0), [a1] "+x,m"(A1));
#else
#define EIGEN_GEBP_2Px8_SPILLING_WORKAROUND
#endif
#define EIGEN_GEBGP_ONESTEP(K) \
do { \
EIGEN_ASM_COMMENT("begin step of gebp micro kernel 2pX8"); \
traits.loadLhs(&blA[(0 + 2 * K) * LhsProgress], A0); \
traits.loadLhs(&blA[(1 + 2 * K) * LhsProgress], A1); \
traits.loadRhs(&blB[(0 + 8 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C0, T0, fix<0>); \
traits.madd(A1, rhs_panel, C8, T0, fix<0>); \
traits.updateRhs(&blB[(1 + 8 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C1, T0, fix<1>); \
traits.madd(A1, rhs_panel, C9, T0, fix<1>); \
traits.updateRhs(&blB[(2 + 8 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C2, T0, fix<2>); \
traits.madd(A1, rhs_panel, C10, T0, fix<2>); \
traits.updateRhs(&blB[(3 + 8 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C3, T0, fix<3>); \
traits.madd(A1, rhs_panel, C11, T0, fix<3>); \
traits.loadRhs(&blB[(4 + 8 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C4, T0, fix<0>); \
traits.madd(A1, rhs_panel, C12, T0, fix<0>); \
traits.updateRhs(&blB[(5 + 8 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C5, T0, fix<1>); \
traits.madd(A1, rhs_panel, C13, T0, fix<1>); \
traits.updateRhs(&blB[(6 + 8 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C6, T0, fix<2>); \
traits.madd(A1, rhs_panel, C14, T0, fix<2>); \
traits.updateRhs(&blB[(7 + 8 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C7, T0, fix<3>); \
traits.madd(A1, rhs_panel, C15, T0, fix<3>); \
EIGEN_GEBP_2Px8_SPILLING_WORKAROUND EIGEN_ASM_COMMENT("end step of gebp micro kernel 2pX8"); \
} while (false)
EIGEN_ASM_COMMENT("begin gebp micro kernel 2pX8");
EIGEN_GEBGP_ONESTEP(0);
EIGEN_GEBGP_ONESTEP(1);
EIGEN_GEBGP_ONESTEP(2);
EIGEN_GEBGP_ONESTEP(3);
EIGEN_GEBGP_ONESTEP(4);
EIGEN_GEBGP_ONESTEP(5);
EIGEN_GEBGP_ONESTEP(6);
EIGEN_GEBGP_ONESTEP(7);
blB += pk * 8 * RhsProgress;
blA += pk * (2 * Traits::LhsProgress);
EIGEN_ASM_COMMENT("end gebp micro kernel 2pX8");
}
// process remaining peeled loop
for (Index k = peeled_kc; k < depth; k++) {
RhsPacketx4 rhs_panel;
RhsPacket T0;
EIGEN_GEBGP_ONESTEP(0);
blB += 8 * RhsProgress;
blA += 2 * Traits::LhsProgress;
}
#undef EIGEN_GEBGP_ONESTEP
ResPacket R0, R1, R2, R3;
ResPacket alphav = pset1<ResPacket>(alpha);
R0 = r0.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r0.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r1.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R3 = r1.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
traits.acc(C0, alphav, R0);
traits.acc(C8, alphav, R1);
traits.acc(C1, alphav, R2);
traits.acc(C9, alphav, R3);
r0.storePacket(0 * Traits::ResPacketSize, R0);
r0.storePacket(1 * Traits::ResPacketSize, R1);
r1.storePacket(0 * Traits::ResPacketSize, R2);
r1.storePacket(1 * Traits::ResPacketSize, R3);
R0 = r2.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r2.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r3.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R3 = r3.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
traits.acc(C2, alphav, R0);
traits.acc(C10, alphav, R1);
traits.acc(C3, alphav, R2);
traits.acc(C11, alphav, R3);
r2.storePacket(0 * Traits::ResPacketSize, R0);
r2.storePacket(1 * Traits::ResPacketSize, R1);
r3.storePacket(0 * Traits::ResPacketSize, R2);
r3.storePacket(1 * Traits::ResPacketSize, R3);
R0 = r4.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r4.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r5.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R3 = r5.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
traits.acc(C4, alphav, R0);
traits.acc(C12, alphav, R1);
traits.acc(C5, alphav, R2);
traits.acc(C13, alphav, R3);
r4.storePacket(0 * Traits::ResPacketSize, R0);
r4.storePacket(1 * Traits::ResPacketSize, R1);
r5.storePacket(0 * Traits::ResPacketSize, R2);
r5.storePacket(1 * Traits::ResPacketSize, R3);
R0 = r6.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r6.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r7.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R3 = r7.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
traits.acc(C6, alphav, R0);
traits.acc(C14, alphav, R1);
traits.acc(C7, alphav, R2);
traits.acc(C15, alphav, R3);
r6.storePacket(0 * Traits::ResPacketSize, R0);
r6.storePacket(1 * Traits::ResPacketSize, R1);
r7.storePacket(0 * Traits::ResPacketSize, R2);
r7.storePacket(1 * Traits::ResPacketSize, R3);
}
}
}
#endif
for (Index j2 = packet_cols8; j2 < packet_cols4; j2 += 4) {
for (Index i = i1; i < actual_panel_end; i += 2 * LhsProgress) {
// We selected a 2*Traits::LhsProgress x nr micro block of res which is entirely
// stored into 2 x nr registers.
const LhsScalar* blA = &blockA[i * strideA + offsetA * (2 * Traits::LhsProgress)];
prefetch(&blA[0]);
// gets res block as register
AccPacket C0, C1, C2, C3, C4, C5, C6, C7;
traits.initAcc(C0);
traits.initAcc(C1);
traits.initAcc(C2);
traits.initAcc(C3);
traits.initAcc(C4);
traits.initAcc(C5);
traits.initAcc(C6);
traits.initAcc(C7);
LinearMapper r0 = res.getLinearMapper(i, j2 + 0);
LinearMapper r1 = res.getLinearMapper(i, j2 + 1);
LinearMapper r2 = res.getLinearMapper(i, j2 + 2);
LinearMapper r3 = res.getLinearMapper(i, j2 + 3);
r0.prefetch(prefetch_res_offset);
r1.prefetch(prefetch_res_offset);
r2.prefetch(prefetch_res_offset);
r3.prefetch(prefetch_res_offset);
// performs "inner" products
const RhsScalar* blB = &blockB[j2 * strideB + offsetB * 4];
prefetch(&blB[0]);
LhsPacket A0, A1;
for (Index k = 0; k < peeled_kc; k += pk) {
EIGEN_ASM_COMMENT("begin gebp micro kernel 2pX4");
RhsPacketx4 rhs_panel;
RhsPacket T0;
// NOTE: the begin/end asm comments below work around bug 935!
// but they are not enough for gcc>=6 without FMA (bug 1637)
#if EIGEN_GNUC_STRICT_AT_LEAST(6, 0, 0) && defined(EIGEN_VECTORIZE_SSE) && !(EIGEN_COMP_LCC)
#define EIGEN_GEBP_2PX4_SPILLING_WORKAROUND __asm__("" : [a0] "+x,m"(A0), [a1] "+x,m"(A1));
#else
#define EIGEN_GEBP_2PX4_SPILLING_WORKAROUND
#endif
#define EIGEN_GEBGP_ONESTEP(K) \
do { \
EIGEN_ASM_COMMENT("begin step of gebp micro kernel 2pX4"); \
traits.loadLhs(&blA[(0 + 2 * K) * LhsProgress], A0); \
traits.loadLhs(&blA[(1 + 2 * K) * LhsProgress], A1); \
traits.loadRhs(&blB[(0 + 4 * K) * RhsProgress], rhs_panel); \
traits.madd(A0, rhs_panel, C0, T0, fix<0>); \
traits.madd(A1, rhs_panel, C4, T0, fix<0>); \
traits.madd(A0, rhs_panel, C1, T0, fix<1>); \
traits.madd(A1, rhs_panel, C5, T0, fix<1>); \
traits.madd(A0, rhs_panel, C2, T0, fix<2>); \
traits.madd(A1, rhs_panel, C6, T0, fix<2>); \
traits.madd(A0, rhs_panel, C3, T0, fix<3>); \
traits.madd(A1, rhs_panel, C7, T0, fix<3>); \
EIGEN_GEBP_2PX4_SPILLING_WORKAROUND \
EIGEN_ASM_COMMENT("end step of gebp micro kernel 2pX4"); \
} while (false)
internal::prefetch(blB + (48 + 0));
EIGEN_GEBGP_ONESTEP(0);
EIGEN_GEBGP_ONESTEP(1);
EIGEN_GEBGP_ONESTEP(2);
EIGEN_GEBGP_ONESTEP(3);
internal::prefetch(blB + (48 + 16));
EIGEN_GEBGP_ONESTEP(4);
EIGEN_GEBGP_ONESTEP(5);
EIGEN_GEBGP_ONESTEP(6);
EIGEN_GEBGP_ONESTEP(7);
blB += pk * 4 * RhsProgress;
blA += pk * (2 * Traits::LhsProgress);
EIGEN_ASM_COMMENT("end gebp micro kernel 2pX4");
}
// process remaining peeled loop
for (Index k = peeled_kc; k < depth; k++) {
RhsPacketx4 rhs_panel;
RhsPacket T0;
EIGEN_GEBGP_ONESTEP(0);
blB += 4 * RhsProgress;
blA += 2 * Traits::LhsProgress;
}
#undef EIGEN_GEBGP_ONESTEP
ResPacket R0, R1, R2, R3;
ResPacket alphav = pset1<ResPacket>(alpha);
R0 = r0.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r0.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r1.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R3 = r1.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
traits.acc(C0, alphav, R0);
traits.acc(C4, alphav, R1);
traits.acc(C1, alphav, R2);
traits.acc(C5, alphav, R3);
r0.storePacket(0 * Traits::ResPacketSize, R0);
r0.storePacket(1 * Traits::ResPacketSize, R1);
r1.storePacket(0 * Traits::ResPacketSize, R2);
r1.storePacket(1 * Traits::ResPacketSize, R3);
R0 = r2.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R1 = r2.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
R2 = r3.template loadPacket<ResPacket>(0 * Traits::ResPacketSize);
R3 = r3.template loadPacket<ResPacket>(1 * Traits::ResPacketSize);
traits.acc(C2, alphav, R0);
traits.acc(C6, alphav, R1);
traits.acc(C3, alphav, R2);
traits.acc(C7, alphav, R3);
r2.storePacket(0 * Traits::ResPacketSize, R0);
r2.storePacket(1 * Traits::ResPacketSize, R1);
r3.storePacket(0 * Traits::ResPacketSize, R2);
r3.storePacket(1 * Traits::ResPacketSize, R3);
}
}
// Deal with remaining columns of the rhs
for (Index j2 = packet_cols4; j2 < cols; j2++) {
for (Index i = i1; i < actual_panel_end; i += 2 * LhsProgress) {
// One column at a time
const LhsScalar* blA = &blockA[i * strideA + offsetA * (2 * Traits::LhsProgress)];
prefetch(&blA[0]);
// gets res block as register
AccPacket C0, C4;
traits.initAcc(C0);
traits.initAcc(C4);
LinearMapper r0 = res.getLinearMapper(i, j2);
r0.prefetch(prefetch_res_offset);
// performs "inner" products
const RhsScalar* blB = &blockB[j2 * strideB + offsetB];
LhsPacket A0, A1;
for (Index k = 0; k < peeled_kc; k += pk) {
EIGEN_ASM_COMMENT("begin gebp micro kernel 2pX1");
RhsPacket B_0, B1;
#define EIGEN_GEBGP_ONESTEP(K) \
do { \
EIGEN_ASM_COMMENT("begin step of gebp micro kernel 2pX1"); \
EIGEN_ASM_COMMENT("Note: these asm comments work around bug 935!"); \
traits.loadLhs(&blA[(0 + 2 * K) * LhsProgress], A0); \
traits.loadLhs(&blA[(1 + 2 * K) * LhsProgress], A1); \
traits.loadRhs(&blB[(0 + K) * RhsProgress], B_0); \
traits.madd(A0, B_0, C0, B1, fix<0>); \
traits.madd(A1, B_0, C4, B_0, fix<0>); \
EIGEN_ASM_COMMENT("end step of gebp micro kernel 2pX1"); \
} while (false)
EIGEN_GEBGP_ONESTEP(0);
EIGEN_GEBGP_ONESTEP(1);
EIGEN_GEBGP_ONESTEP(2);
EIGEN_GEBGP_ONESTEP(3);
EIGEN_GEBGP_ONESTEP(4);
EIGEN_GEBGP_ONESTEP(5);
EIGEN_GEBGP_ONESTEP(6);
EIGEN_GEBGP_ONESTEP(7);
blB += int(