如何使用 C 中的 SSE 内在函数计算单向量点积

Question

我试图将两个向量相乘，其中一个向量的每个元素都乘以另一个向量的相同索引中的元素。然后我想对结果向量的所有元素求和以获得一个数字。例如，向量 {1,2,3,4} 和 {5,6,7,8} 的计算如下所示：

1*5 + 2*6 + 3*7 + 4*8

本质上，我取两个向量的点积。我知道有一个 SSE 命令可以执行此操作，但该命令没有与之关联的内在函数。此时，我不想在 C 代码中编写内联汇编，因此我只想使用内部函数。这似乎是一个常见的计算，所以我对自己在谷歌上找不到答案感到惊讶。

注意：我正在针对支持 SSE 4.2 的特定微架构进行优化。

Answer 1

如果您要对较长向量进行点积，请在内循环中使用乘法和常规 _mm_add_ps（或 FMA）。保存水平总和直至结束。

---

但是，如果您只对一对 SIMD 向量进行点积：

GCC（至少版本 4.3）包含带有 SSE4.1 级别内在函数的 ，包括 single 和 double - 精度点积：

_mm_dp_ps (__m128 __X, __m128 __Y, const int __M);
_mm_dp_pd (__m128d __X, __m128d __Y, const int __M);

在 Intel 主流 CPU（不是 Atom/Silvermont）上，这些比使用多条指令手动执行要快一些。

但在 AMD（包括 Ryzen）上，dpps 速度明显慢一些。 （请参阅Agner Fog 的指令表）

---

作为旧处理器的后备方案，您可以使用此算法来创建点积向量 a 和 b：

__m128 r1 = _mm_mul_ps(a, b);

然后使用 在 x86 上进行水平浮点向量和的最快方法（请参阅此处的评论版本，以及为什么它更快。）

__m128 shuf   = _mm_shuffle_ps(r1, r1, _MM_SHUFFLE(2, 3, 0, 1));
__m128 sums   = _mm_add_ps(r1, shuf);
shuf          = _mm_movehl_ps(shuf, sums);
sums          = _mm_add_ss(sums, shuf);
float result =  _mm_cvtss_f32(sums);

较慢的替代方案每个 hadd 需要 2 次 shuffle，这很容易成为 shuffle 吞吐量的瓶颈，尤其是在 Intel CPU 上。

r2 = _mm_hadd_ps(r1, r1);
r3 = _mm_hadd_ps(r2, r2);
_mm_store_ss(&result, r3);

Answer 2

我想说最快的 SSE 方法是：

static inline float CalcDotProductSse(__m128 x, __m128 y) {
    __m128 mulRes, shufReg, sumsReg;
    mulRes = _mm_mul_ps(x, y);

    // Calculates the sum of SSE Register - https://stackoverflow.com/a/35270026/195787
    shufReg = _mm_movehdup_ps(mulRes);        // Broadcast elements 3,1 to 2,0
    sumsReg = _mm_add_ps(mulRes, shufReg);
    shufReg = _mm_movehl_ps(shufReg, sumsReg); // High Half -> Low Half
    sumsReg = _mm_add_ss(sumsReg, shufReg);
    return  _mm_cvtss_f32(sumsReg); // Result in the lower part of the SSE Register
}

我遵循 - 在 x86 上进行水平浮点向量和的最快方法。

Answer 3

我写了这个并用 gcc -O3 -S -ftree-vectorize -ftree-vectorizer-verbose=2 sse.c 编译它

void f(int * __restrict__ a, int * __restrict__ b, int * __restrict__ c, int * __restrict__ d,
       int * __restrict__ e, int * __restrict__ f, int * __restrict__ g, int * __restrict__ h,
       int * __restrict__ o)
{
    int i;

    for (i = 0; i < 8; ++i)
        o[i] = a[i]*e[i] + b[i]*f[i] + c[i]*g[i] + d[i]*h[i];
}

并且 GCC 4.3.0 自动矢量化它：

sse.c:5: note: LOOP VECTORIZED.
sse.c:2: note: vectorized 1 loops in function.

但是，它只会做如果我使用具有足够迭代次数的循环，否则详细的输出将表明矢量化无利可图或循环太小。如果没有 __restrict__ 关键字，它必须生成单独的非向量化版本来处理输出 o 可能指向输入之一的情况。

我将粘贴这些指令作为示例，但由于矢量化的一部分展开了循环，因此可读性不太好。

Answer 4

英特尔此处有一篇文章涉及点积实现。