c++ - 用于 float 阈值操作的 SIMD

Question

我想让一些 vector 计算更快，我相信用于浮点比较和操作的 SIMD 指令会有所帮助，这是操作:

void func(const double* left, const double* right, double* res, const size_t size, const double th, const double drop) {
        for (size_t i = 0; i < size; ++i) {
            res[i] = right[i] >= th ? left[i] : (left[i] - drop) ;
        }
    }

主要是，如果 left 值高于 drop ，它会通过 right 删除 threshold 值。

大小在 128-256 左右(不是很大)，但计算量很大。

我尝试从循环展开开始，但没有赢得很多性能，但可能需要一些编译指令。

您能否建议对代码进行一些改进以加快计算速度？

Answer 1

Clang 已经自动矢量化了，这几乎就像 Soonts 建议的手动做的那样。 在您的指针上使用 __restrict，因此它不需要用于某些数组之间重叠的回退版本。它仍然自动矢量化，但它使函数膨胀。

不幸的是，gcc 只使用 -ffast-math 自动矢量化。结果证明只需要 -fno-trapping-math :这通常是安全的，特别是如果您不使用 fenv 访问来取消屏蔽任何 FP 异常 ( feenableexcept ) 或查看 MXCSR 粘性 FP 异常标志 ( fetestexcept )。

使用该选项，GCC 也将使用 (v)pblendvpd 和 -march=nehalem 或 -march=znver1 。 See it on Godbolt

此外，您的 C 函数已损坏。 th 和 drop 是双标量，但您将它们声明为 const double *
AVX512F 可以让您进行 !(right[i] >= thresh) 比较并使用结果掩码进行合并掩码减法。

谓词为真的元素将获得 left[i] - drop ，其他元素将保留其 left[i] 值，因为您将 info 合并为 left 值的 vector 。

不幸的是，带有 -march=skylake-avx512 的 GCC 使用一个普通的 vsubpd 然后一个单独的 vmovapd zmm2{k1}, zmm5 来混合，这显然是一个错过的优化。混合目标已经是 SUB 的输入之一。

将 AVX512VL 用于 256 位 vector (以防您的程序的其余部分无法有效地使用 512 位，因此您不会降低涡轮时钟速度):

__m256d left = ...;
__m256d right = ...;
__mmask8 cmp = _mm256_cmp_pd_mask(right, set1(th), _CMP_NGE_UQ);
__m256d res = _mm256_mask_sub_pd (left, cmp, left, set1(drop));

所以(除了加载和存储)它是 AVX512F/VL 的 2 条指令。

如果您不需要版本的特定 NaN 行为，GCC 也可以自动矢量化

它对所有编译器都更有效率，因为您只需要一个 AND，而不是一个可变混合。 因此，仅使用 SSE2 就明显更好，并且即使在大多数 CPU 支持 SSE4.1 blendvpd 时也更好，因为该指令效率不高。

您可以根据比较结果从 0.0 中减去 drop 或 left[i]。

根据比较结果生成 0.0 或常量非常有效:只是一条 andps 指令。 ( 0.0 的位模式是全零，SIMD 比较产生全 1 或全 0 位的 vector 。因此 AND 保留旧值或将其归零。)

我们还可以添加 -drop 而不是减去 drop 。这需要对输入进行额外的否定，但 AVX 允许 vaddpd 的内存源操作数。不过，GCC 选择使用索引寻址模式，因此实际上无助于减少 Intel CPU 上的前端 uop 计数；它会“解压”。但即使使用 -ffast-math ，gcc 也不会自行进行这种优化以允许折叠负载。 (不过，除非我们展开循环，否则不值得进行单独的指针增量。)

void func3(const double *__restrict left, const double *__restrict right, double *__restrict res,
  const size_t size, const double th, const double drop)
{
    for (size_t i = 0; i < size; ++i) {
        double add = right[i] >= th ? 0.0 : -drop;
        res[i] = left[i] + add;
    }
}

GCC 9.1 的内部循环(没有任何 -march 选项并且没有 -ffast-math )来自上面的 Godbolt 链接:

# func3 main loop
# gcc -O3 -march=skylake       (without fast-math)
.L33:
    vcmplepd        ymm2, ymm4, YMMWORD PTR [rsi+rax]
    vandnpd ymm2, ymm2, ymm3
    vaddpd  ymm2, ymm2, YMMWORD PTR [rdi+rax]
    vmovupd YMMWORD PTR [rdx+rax], ymm2
    add     rax, 32
    cmp     r8, rax
    jne     .L33

或者普通的 SSE2 版本有一个内部循环，它与 left - zero_or_drop 而不是 left + zero_or_minus_drop 基本相同，所以除非你能保证编译器 16 字节对齐或者你正在制作 AVX 版本，否则否定 drop 只是额外的开销。

否定 drop 从内存中获取一个常量(对符号位进行异或)，这是该函数需要的唯一静态常量 ，因此对于循环运行次数不多的情况，值得考虑权衡。 (除非 th 或 drop 在内联后也是编译时常量，并且无论如何都会被加载。或者特别是如果 -drop 可以在编译时计算。或者如果你可以让你的程序使用负值 drop 。)

加法和减法的另一个区别是减法不会破坏零的符号。 -0.0 - 0.0 = -0.0 , +0.0 - 0.0 = +0.0 。万一那很重要。

# gcc9.1 -O3
.L26:
    movupd  xmm5, XMMWORD PTR [rsi+rax]
    movapd  xmm2, xmm4                    # duplicate  th
    movupd  xmm6, XMMWORD PTR [rdi+rax]
    cmplepd xmm2, xmm5                    # destroy the copy of th
    andnpd  xmm2, xmm3                    # _mm_andnot_pd
    addpd   xmm2, xmm6                    # _mm_add_pd
    movups  XMMWORD PTR [rdx+rax], xmm2
    add     rax, 16
    cmp     r8, rax
    jne     .L26

GCC 使用未对齐的加载，因此(没有 AVX)它不能将内存源操作数折叠成 cmppd 或 subpd