SIMD 减少 4 个向量而没有 hadd

Posted 2023-02-16

【中文标题】SIMD 减少 4 个向量而没有 hadd

【英文标题】：SIMD reduce 4 vectors without hadd

【发布时间】：2020-03-21 17:28:21

【问题描述】：

我正在尝试优化一些代码，并且我处于有 4 个向量 __m256d 的状态，我想将每个向量的总和存储在另一个 __m256d 中。 所以基本上result = [sum(a), sum(b), sum(c), sum(d)]。我知道有一种方法可以使用 2 个 hadds 混合和置换来做到这一点，但我意识到 hadd 太贵了。

所以我想知道是否有一个内在函数可以更快地做到这一点。

【问题讨论】：

这能回答你的问题吗？ SIMD optimization of a curve computed from the second derivative

Mmh 不，因为最终结果不是像我的情况那样必须汇总的 4 个变量，但它们利用了问题定义。还是我错过了他们如何按照我的要求做的事情？

您可以使用 2 vshufpd + 垂直添加来模拟 hadd，但这与某些 CPU 的成本完全相同。 （尽管在 Zen1 上 IIRC 稍微便宜一点，其中 hadd 出于某种原因每条通道 4 微秒。）转置和求和本质上并不便宜，但是将向量混洗在一起以提供垂直添加是您最好的选择。也许 2x vperm2f128 可能会很好，虽然这在 Zen1 上很慢。

你需要转置你的矩阵，然后一个普通的 SIMD "add" 就可以了。

@PeterCordes 您能否详细说明 2x vperm2f128。因为我尝试使用它以及一些混合和 3 个添加，但我得到了与使用 hadds 相同的性能。也许以您想要的方式使用 vperm2f128 更快。不胜感激，谢谢

【参考方案1】：

三个选项：

1 矩阵转置，然后垂直求和

好：概念上简单，使用普遍有用的算法（矩阵转置），可移植代码

不好：代码大小、延迟、吞吐量

2 有效地使用vhaddpd

好：小代码（适用于 Icache）、英特尔 uArchs 上的良好延迟和吞吐量

不好：需要特定于架构的代码，在某些 uArch 上存在问题

3部分转置、求和、部分转置、求和

good：良好的延迟，良好的吞吐量

不好：不像vhaddpd-code那么小，不像完整的矩阵转置那么容易理解

矩阵转置，垂直和

让你的编译器为你优化这个。使用gcc 向量扩展*，对转置矩阵求和的代码可能如下所示：

#include <stdint.h>

typedef uint64_t v4u64 __attribute__((vector_size(32)));
typedef double v4f64  __attribute__((vector_size(32)));

v4f64 dfoo(v4f64 sv0, v4f64 sv1, v4f64 sv2, v4f64 sv3)

  v4f64 tv[4];
  tv[0] = __builtin_shuffle(sv0, sv1, (v4u64)0,4,2,6);
  tv[1] = __builtin_shuffle(sv0, sv1, (v4u64)1,5,3,7);
  tv[2] = __builtin_shuffle(sv2, sv3, (v4u64)0,4,2,6);
  tv[3] = __builtin_shuffle(sv2, sv3, (v4u64)1,5,3,7);
  v4f64 fv[4];
  fv[0] = __builtin_shuffle(tv[0], tv[2], (v4u64)0,1,4,5);
  fv[1] = __builtin_shuffle(tv[0], tv[2], (v4u64)2,3,6,7);
  fv[2] = __builtin_shuffle(tv[1], tv[3], (v4u64)0,1,4,5);
  fv[3] = __builtin_shuffle(tv[1], tv[3], (v4u64)2,3,6,7);
  return fv[0]+fv[1]+fv[2]+fv[3];

gcc-9.2.1 生成以下程序集：

dfoo:
    vunpcklpd   %ymm3, %ymm2, %ymm5
    vunpcklpd   %ymm1, %ymm0, %ymm4
    vunpckhpd   %ymm1, %ymm0, %ymm0
    vinsertf128 $1, %xmm5, %ymm4, %ymm1
    vperm2f128  $49, %ymm5, %ymm4, %ymm4
    vunpckhpd   %ymm3, %ymm2, %ymm2
    vaddpd  %ymm4, %ymm1, %ymm1
    vinsertf128 $1, %xmm2, %ymm0, %ymm3
    vperm2f128  $49, %ymm2, %ymm0, %ymm0
    vaddpd  %ymm3, %ymm1, %ymm1
    vaddpd  %ymm0, %ymm1, %ymm0
    ret

Agner Fog 的表格说：

vunpck[h/l]pd：1 个周期延迟，每周期 1 个吞吐量，1 个 uOP 端口5。 vinsertf128：3 个周期延迟，每周期 1 个吞吐量，1 个 uOP 端口5。 vperm2f128：3 个周期延迟，每周期 1 个吞吐量，1 个 uOP 端口5。 vaddpd：4 个周期延迟，每个周期 2 个吞吐量，1 个 uOP 端口01。

总之有

4 [解包] + 2 [插入] + 2 [置换] = 8 个端口 5 uOP。 3 [添加] = 3 port01 uOP。

端口 5 的吞吐量将成为瓶颈。 大约 18 个周期的延迟非常糟糕。 代码大小约为 60 字节。

水平总和

使用vhadd 的代码（明智地）不容易通过 gcc 向量扩展获得，因此代码需要特定于 Intel 的内在函数：

v4f64 dfoo_hadd(v4f64 sv0, v4f64 sv1, v4f64 sv2, v4f64 sv3)

  v4f64 hv[2];
  hv[0] = __builtin_ia32_haddpd256(sv0, sv1); //[00+01, 10+11, 02+03, 12+13]
  hv[1] = __builtin_ia32_haddpd256(sv2, sv3); //[20+21, 30+31, 22+23, 32+33]
  v4f64 fv[2];
  fv[0] = __builtin_shuffle(hv[0], hv[1], (v4u64)0, 1, 4, 5); //[00+01, 10+11, 20+21, 30+31]
  fv[1] = __builtin_shuffle(hv[0], hv[1], (v4u64)2, 3, 6, 7); //[02+03, 12+13, 22+23, 32+33]
  return fv[0] + fv[1]; //[00+01+02+03, 10+11+12+13, 20+21+22+23, 30+31+32+33]

这会生成以下程序集：

dfoo_hadd:
    vhaddpd %ymm3, %ymm2, %ymm2
    vhaddpd %ymm1, %ymm0, %ymm0
    vinsertf128 $1, %xmm2, %ymm0, %ymm1
    vperm2f128  $49, %ymm2, %ymm0, %ymm0
    vaddpd  %ymm0, %ymm1, %ymm0
    ret

根据 Agner Fog 的指令表，

vhaddpd：6 个周期延迟，每周期 0.5 个吞吐量，3 uOPS port01 + 2*port5。

总之有

4 [hadd] + 2 [插入/置换] = 6 uOPs port5。 3 [hadd/add] = 3 uOPs port01。

吞吐量也受到端口5的限制，这比转置代码的吞吐量更大。 延迟应该约为 16 个周期，也比转置代码快。 代码大小约为 25 字节。

部分转置、总和、部分转置、总和

实现@PeterCordes 评论：

v4f64 dfoo_PC(v4f64 sv0, v4f64 sv1, v4f64 sv2, v4f64 sv3)

  v4f64 tv[4];
  v4f64 av[2];
  tv[0] = __builtin_shuffle(sv0, sv1, (v4u64)0,4,2,6);//[00, 10, 02, 12]
  tv[1] = __builtin_shuffle(sv0, sv1, (v4u64)1,5,3,7);//[01, 11, 03, 13]
  av[0] = tv[0] + tv[1];//[00+01, 10+11, 02+03, 12+13]
  tv[2] = __builtin_shuffle(sv2, sv3, (v4u64)0,4,2,6);//[20, 30, 22, 32]
  tv[3] = __builtin_shuffle(sv2, sv3, (v4u64)1,5,3,7);//[21, 31, 23, 33]
  av[1] = tv[2] + tv[3];//[20+21, 30+31, 22+23, 32+33]
  v4f64 fv[2];
  fv[0] = __builtin_shuffle(av[0], av[1], (v4u64)0,1,4,5);//[00+01, 10+11, 20+21, 30+31]
  fv[1] = __builtin_shuffle(av[0], av[1], (v4u64)2,3,6,7);//[02+03, 12+13, 22+23, 32+33]
  return fv[0]+fv[1];//[00+01+02+03, 10+11+12+13, 20+21+22+23, 30+31+32+33]

这会生成：

dfoo_PC:
    vunpcklpd   %ymm1, %ymm0, %ymm4
    vunpckhpd   %ymm1, %ymm0, %ymm1
    vunpcklpd   %ymm3, %ymm2, %ymm0
    vunpckhpd   %ymm3, %ymm2, %ymm2
    vaddpd  %ymm1, %ymm4, %ymm1
    vaddpd  %ymm2, %ymm0, %ymm2
    vinsertf128 $1, %xmm2, %ymm1, %ymm0
    vperm2f128  $49, %ymm2, %ymm1, %ymm1
    vaddpd  %ymm1, %ymm0, %ymm0
    ret

总之有

4 [解包] + 2 [插入/置换] = 6 个端口 5 uOP。 3 [添加] = 3 port01 uOP。

这会得到与hadd-code 相同数量的 port5 uOP。代码仍然是端口 5 的瓶颈，延迟约为 16 个周期。 代码大小约为 41 字节。

如果您想提高吞吐量，则必须将工作从端口 5 转移出去。不幸的是，几乎所有的 permute/insert/shuffle 指令都需要端口 5，并且车道交叉指令（此处需要）具有至少 3 个周期的延迟。 几乎 有帮助的一条有趣指令是vblendpd，它具有 3 个/周期的吞吐量、1 个周期的延迟，并且可以在端口 015 上执行，但使用它来替换其中一个置换/插入/随机播放需要向量的 128 位通道的 64 位移位，由 vpsrldq/vpslldq 实现，你猜对了，它需要一个端口 5 uOP（所以这个 有助于 32 位向量位float，因为vpsllq/vpsrlq 不需要端口5）。这里没有免费的午餐。

* gcc 向量扩展快速说明：

代码使用 gcc 向量扩展，允许在向量上使用基本运算符（+-*/=><>><< 等），按元素操作。它们还包括一些 __builtin_* 函数，特别是 __builtin_shuffle()，它具有 3 操作数形式，其中前两个是相同类型 T 的两个（相同长度 N）向量，它们（逻辑上）连接到T 类型的双倍长度 (2N) 向量，第三个是整数类型 (IT) 的向量，其宽度和长度 (N) 与原始向量的类型相同。结果是一个与原始向量具有相同类型 T 和宽度 N 的向量，其元素由整数类型向量中的索引选择。

最初，我的回答是关于uint64_t，留在这里作为上下文：

 #include <stdint.h>

typedef uint64_t v4u64 __attribute__((vector_size(32)));

v4u64 foo(v4u64 sv0, v4u64 sv1, v4u64 sv2, v4u64 sv3)

  v4u64 tv[4];
  tv[0] = __builtin_shuffle(sv0, sv1, (v4u64)0,4,2,6);
  tv[1] = __builtin_shuffle(sv0, sv1, (v4u64)1,5,3,7);
  tv[2] = __builtin_shuffle(sv2, sv3, (v4u64)0,4,2,6);
  tv[3] = __builtin_shuffle(sv2, sv3, (v4u64)1,5,3,7);
  v4u64 fv[4];
  fv[0] = __builtin_shuffle(tv[0], tv[2], (v4u64)0,1,4,5);
  fv[1] = __builtin_shuffle(tv[0], tv[2], (v4u64)2,3,6,7);
  fv[2] = __builtin_shuffle(tv[1], tv[3], (v4u64)0,1,4,5);
  fv[3] = __builtin_shuffle(tv[1], tv[3], (v4u64)2,3,6,7);
  return fv[0]+fv[1]+fv[2]+fv[3];

gcc-9.2.1 在 skylake-avx2 上生成的翻译可能如下所示：

foo:
    vpunpcklqdq %ymm3, %ymm2, %ymm5
    vpunpcklqdq %ymm1, %ymm0, %ymm4
    vpunpckhqdq %ymm3, %ymm2, %ymm2
    vpunpckhqdq %ymm1, %ymm0, %ymm0
    vperm2i128  $32, %ymm2, %ymm0, %ymm3
    vperm2i128  $32, %ymm5, %ymm4, %ymm1
    vperm2i128  $49, %ymm2, %ymm0, %ymm0
    vperm2i128  $49, %ymm5, %ymm4, %ymm4
    vpaddq  %ymm4, %ymm1, %ymm1
    vpaddq  %ymm0, %ymm3, %ymm0
    vpaddq  %ymm0, %ymm1, %ymm0
    ret

请注意，该程序集几乎有一行与 gcc 向量扩展的行对应。

根据 Agner Fog 的 Skylake 指令表，

vpunpck[h/l]qdq：1 个周期延迟，每周期 1 个吞吐量，端口 5。 vperm2i128：3 个周期延迟，每个周期 1 个吞吐量，端口 5。 vpaddq：1 个周期延迟，每周期 3 个吞吐量，端口 015。

因此转置需要 10 个周期（解包 4 个，吞吐量 4 个 + 置换延迟 2 个）。三个加法中，只有两个可以并行执行，所以需要 2 个周期，总共 12 个周期。

【讨论】：

__m256d 是双精度，而不是 int64_t。 FP 版本的所有功能都具有相同的性能，除了 FP add 的吞吐量为 0.5c，SKL 的延迟为 4c。

您的性能分析应该更清楚地说明“需要 12 个周期”是延迟，而不是吞吐量。如果您为 independent 输入背靠背地执行此操作，则可以每 8 个时钟进行一次迭代（4 个向量 -> 1 个），仅在 shuffle 吞吐量上成为瓶颈。剩余大量前端带宽用于加载/存储结果。

顺便说一句，我们有没有办法可以hadd 第一步，减少到只有 2 个向量？或者这是否会使洗牌复杂化，以至于你需要 AVX512 vpermt2pd 或其他东西来在车道上排列正确的元素以进行最终的垂直添加，或者像 AVX2 vpermpd 在 vperm2f128 之后进行额外的车道交叉洗牌？

我忘了提到我使用双打，所以vpunpcklqdq 不会很好。可能会使用vunpcklpd，但吞吐量更差。第一个代码块代表什么？因为我无法理解排列是如何准确执行的

@PeterCordes 好点。添加了一个在改组之间添加的版本，与hadd port5 利用率相匹配，谢谢。