为什么向量长度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])

#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;
}

我用GCC(Ubuntu7.4.0-1Ubuntu1~18.04.1)7.4.0:

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

SSE版本输出：

499999999500000128.000000
13458ms

纯C版本输出：

499999999500000128.000000
4441ms

Answer 1

最明显的问题是使用效率低下的点积（使用haddps，其成本为2倍shuffle uops+1x add uop)，而不是shuffle+add。关于在_mm_mul_ps之后要做什么，请参见在x86上做水平浮点向量和的最快方法，这并不太费劲。但这仍然不是x86能够非常有效地完成的事情。

但无论如何，真正的问题是您的基准循环。

A[3]=K；然后使用_mm_load_ps(A)创建一个存储转发延迟，如果它是简单编译而不是向量洗牌的话。如果加载只从单个存储指令加载数据，而没有其他数据，则存储+重新加载可以以大约5个延迟周期有效地转发。否则，它必须对整个存储缓冲区进行较慢的扫描来组装字节。这会给存储转发增加大约10个延迟周期。

参见Agner Fog的microarch指南和/或优化手册。

• Intel开发人员手册中的“存储缓冲区转发”是什么意思？
• 在内存访问不对齐的情况下，存储区如何进行加载转发？
• 现代x86实现是否可以从一个以上的先前存储进行存储转发？
• 编译器为什么要生成这个程序集？引用英特尔的优化手册RE：存储转发。（在这个问题中，旧的gcc版本分别存储了一个8字节结构的两个dword半部分，然后用qword加载/存储。Super braindead复制该结构。）

另外，通过让编译器将v[0]*v[0]+v[1]*v[1]+v[2]*v[2]的计算从循环中移出，您对SSE的影响更大。

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

我不知道为什么海合会认为这是个好主意。对于将向量加载分成8字节的CPU（如Pentium M或Bobcat）可能会有所帮助，以避免存储转发停滞。但对于“通用”的现代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与scalar具有相同的每3个循环1的吞吐量，因此通过一次水平执行1个向量，而不是并行执行4个向量的4个长度，您只能获得CPU sqrt吞吐量能力的1/4。