为什么向量长度SIMD代码比普通的C语言代码更慢

Question

为什么我的SIMD向量4长度函数比朴素的向量长度方法慢3倍？

SIMD向量4长度函数：

__extern_always_inline float vec4_len(const float *v) {
    __m128 vec1 = _mm_load_ps(v);
    __m128 xmm1 = _mm_mul_ps(vec1, vec1);
    __m128 xmm2 = _mm_hadd_ps(xmm1, xmm1);
    __m128 xmm3 = _mm_hadd_ps(xmm2, xmm2);
    return sqrtf(_mm_cvtss_f32(xmm3));
}

天真的实现：

sqrtf(V[0] * V[0] + V[1] * V[1] + V[2] * V[2] + V[3] * V[3])

SIMD版本迭代1000000000次需要16110毫秒。相比之下，朴素版本快了约3倍，只需4746毫秒。

#include <math.h>
#include <time.h>
#include <stdint.h>
#include <stdio.h>
#include <x86intrin.h>

static float vec4_len(const float *v) {
    __m128 vec1 = _mm_load_ps(v);
    __m128 xmm1 = _mm_mul_ps(vec1, vec1);
    __m128 xmm2 = _mm_hadd_ps(xmm1, xmm1);
    __m128 xmm3 = _mm_hadd_ps(xmm2, xmm2);
    return sqrtf(_mm_cvtss_f32(xmm3));
}

int main() {
    float A[4] __attribute__((aligned(16))) = {3, 4, 0, 0};

    struct timespec t0 = {};
    clock_gettime(CLOCK_MONOTONIC, &t0);

    double sum_len = 0;
    for (uint64_t k = 0; k < 1000000000; ++k) {
        A[3] = k;
        sum_len += vec4_len(A);
//        sum_len += sqrtf(A[0] * A[0] + A[1] * A[1] + A[2] * A[2] + A[3] * A[3]);
    }
    struct timespec t1 = {};
    clock_gettime(CLOCK_MONOTONIC, &t1);

    fprintf(stdout, "%f\n", sum_len);

    fprintf(stdout, "%ldms\n", (((t1.tv_sec - t0.tv_sec) * 1000000000) + (t1.tv_nsec - t0.tv_nsec)) / 1000000);

    return 0;
}

我在一台搭载Intel(R) Core(TM) i7-8550U处理器的计算机上，使用以下命令运行。首先是使用vec4_len版本，然后是使用普通的C语言。

我使用GCC (Ubuntu 7.4.0-1ubuntu1~18.04.1) 7.4.0进行编译：

Original Answer翻译成"最初的回答"

gcc -Wall -Wextra -O3 -msse -msse3 sse.c -lm && ./a.out

SSE版本输出：

499999999500000128.000000
13458ms

普通C版本输出：

最初的回答：

499999999500000128.000000
4441ms

- user1679740

---

"haddps"很糟糕，虽然3倍听起来还是很多，但反汇编中有什么有趣的东西吗？ - undefined

你是如何进行基准测试的？你确定标量版本没有提升一些工作吗？使用了什么编译器/版本和什么硬件，以及这个内联到什么周围代码中？另外，由于你的值在__m128i中，可以通过使用_mm_sqrt_ss避免sqrtf()的愚蠢的设置NaN时的errno行为，因此您无需使用-fno-math-errno编译。 - undefined

需要编译器版本。需要提供最小可复现代码。Godbolt 链接很好用。 - undefined

你的实际用途是什么？偶尔计算单个向量的范数，还是计算一系列向量的范数？你能重新组织你的数据吗（SoA vs AoS）？ - undefined

看起来编译器在优化方面比你更好，godbolt。你的函数对于任意数据似乎更快（汇编代码更短），但是你测量执行速度的方法有点无效，编译器优化了循环。你需要使用两个函数并告诉编译器副作用。例如使用__attribute__((__noinline__))，就像这个例子中的做法一样。 - undefined

@KamilCuk：但是这个汇编使用了haddps，这个指令很糟糕，尤其对于Sandybridge系列来说。它需要2个洗牌微操作和1个垂直加法微操作，而SnB系列每个时钟周期只能执行1个洗牌微操作（在端口5上）。直到IceLake在另一个端口上添加第二个洗牌单元之前，这种情况下较短的汇编代码并不总是更好的选择。 - undefined

Answer 1

最明显的问题是使用效率低下的点积（使用haddps，需要2个乘法和1个加法操作），而不是使用洗牌+加法。请参见在x86上进行浮点向量水平求和的最快方法，了解在_mm_mul_ps之后做什么可以减少成本。但是，这仅是x86不能非常高效地执行的东西。

但是，真正的问题在于您的基准测试循环。 A[3] = k;然后使用_mm_load_ps(A)会创建存储-转发延迟，如果编译器朴素地编译而不是矢量洗牌。如果加载仅从单个存储指令加载数据且没有加载该指令外的任何数据，则可以有效地转发存储+重新加载，并具有约5个周期的延迟。否则，它必须对整个存储缓冲区进行较慢的扫描以组装字节。这将增加大约10个周期的延迟到存储转发中。

我不确定这对吞吐量有多大影响，但足以阻止乱序执行重叠足够的循环迭代以隐藏延迟，并仅在sqrtss洗牌吞吐量上成为瓶颈。

（您的Coffee Lake CPU具有每3个周期的1个sqrtss吞吐量，因此令人惊讶的是，SQRT吞吐量不是瓶颈。¹ 相反，它将是洗牌吞吐量或其他因素。）

请参见Agner Fog的微架构指南和/或优化手册。

• Intel开发人员手册中的“存储缓冲区转发”是什么意思？
• 在非对齐内存访问的情况下，存储到加载转发如何进行？
• 现代x86实现是否可以从多个先前的存储中进行存储转发？
• 为什么编译器会生成这个汇编代码？引用了英特尔的优化手册，其中提到了存储转发。（在该问题中，旧版本的gcc将8字节结构的2个双字单独存储，然后使用qword加载/存储复制结构。非常低效）

---

此外，您让编译器在循环之外执行V [0] * V [0] + V [1] * V [1] + V [2] * V [2]的计算，这使得SSE更加不利。

该表达式的该部分与循环无关，因此编译器每次循环迭代只需要进行(float)k平方、加和标量平方根运算（并将结果转换为double以添加到累加器中）。

（@StaceyGirl的删除答案指出了这一点；查看内部循环代码是撰写本答案的很好的开端。）

---

向量版本中A[3]=k存在的额外低效率

GCC9.1的内部循环如Kamil的Godbolt链接所示，看起来很糟糕，并且似乎包括一个循环传递的存储/加载操作，将新的A [3]合并到8字节的A [2..3]对中，进一步限制了CPU重叠多个迭代的能力。

我不确定gcc为什么认为这是一个好主意。这可能有助于像Pentium M或Bobcat这样将向量加载拆分成8字节一半的CPU避免存储转发停顿。但这对于“通用”的现代x86-64 CPU来说不是一个明智的调优。

.L18:
        pxor    xmm4, xmm4
        mov     rdx, QWORD PTR [rsp+8]     ; reload A[2..3]
        cvtsi2ss        xmm4, rbx
        mov     edx, edx                   ; truncate RDX to 32-bit
        movd    eax, xmm4                  ; float bit-pattern of (float)k
        sal     rax, 32
        or      rdx, rax                   ; merge the float bit-pattern into A[3]
        mov     QWORD PTR [rsp+8], rdx     ; store A[2..3] again

        movaps  xmm0, XMMWORD PTR [rsp]    ; vector load: store-forwarding stall
        mulps   xmm0, xmm0
        haddps  xmm0, xmm0
        haddps  xmm0, xmm0
        ucomiss xmm3, xmm0
        movaps  xmm1, xmm0
        sqrtss  xmm1, xmm1
        ja      .L21             ; call sqrtf to set errno if needed; flags set by ucomiss.
.L17:

        add     rbx, 1
        cvtss2sd        xmm1, xmm1
        addsd   xmm2, xmm1            ; total += (double)sqrtf
        cmp     rbx, 1000000000
        jne     .L18                ; }while(k<1000000000);

在标量版本中不存在这种疯狂的情况。

无论如何，gcc成功避免了完整的uint64_t -> float转换的低效率（在AVX512之前的x86硬件上没有）。它可能能够证明使用有符号64位 -> 浮点数转换总是有效的，因为高位不能被设置。

---

注脚1：但sqrtps与标量具有相同的每3周期吞吐量，因此通过水平地一次处理1个向量而不是同时处理4个向量的4个长度，您只能获得CPU平方根吞吐量能力的1/4。

- Peter Cordes