用于浮点阈值操作的 SIMD

SIMD for float threshold operation

本文关键字：SIMD 操作 阈值 用于      更新时间：2023-10-16

我想使一些矢量计算更快，我相信用于浮点比较和操作的 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) ;
}
}

主要是，如果该值高于thresholdright它会将left值降低drop。

大小约为 128-256(不是那么大)，但计算被大量调用。

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

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

Clang已经按照Soonts建议的手动方式自动矢量化了。在指针上使用__restrict，这样它就不需要用于某些数组之间重叠的回退版本。 它仍然会自动矢量化，但它使函数膨胀。

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

使用该选项，GCC 也将(v)pblendvpd与-march=nehalem或-march=znver1一起使用。在Godbolt上看到它

此外，您的 C 函数已损坏。th和drop是标量双精度，但您将它们声明为const double *

AVX512F可以让您进行!(right[i] >= thresh)比较并使用生成的掩码进行合并掩码减法。

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

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

将AVX512VL用于 256 位矢量(以防程序的其余部分无法有效地使用 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));

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

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

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

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

根据比较结果生成0.0或常量非常有效：只需一条andps指令。 (0.0的位模式为全零，SIMD 比较生成的全 1 位或全 0 位的向量。 所以 AND 保留旧值或将其归零。

我们也可以加-drop而不是减去drop. 这需要对输入进行额外的否定，但使用 AVX 允许vaddpd的内存源操作数。 不过，GCC 选择使用索引寻址模式，因此这实际上无助于减少英特尔 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从内存中获取一个常量(到 XOR 符号位)，这是此函数需要的唯一静态常量，因此对于循环运行次数不多的情况，值得考虑权衡。 (除非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

你去吧(未经测试)，我试图在评论中解释他们的作用。

void func_sse41( const double* left, const double* right, double* res,
const size_t size, double th, double drop )
{
// Verify the size is even.
// If it's not, you'll need extra code at the end to process last value the old way.
assert( 0 == ( size % 2 ) );
// Load scalar values into 2 registers.
const __m128d threshold = _mm_set1_pd( th );
const __m128d dropVec = _mm_set1_pd( drop );
for( size_t i = 0; i < size; i += 2 )
{
// Load 4 double values into registers, 2 from right, 2 from left
const __m128d r = _mm_loadu_pd( right + i );
const __m128d l = _mm_loadu_pd( left + i );
// Compare ( r >= threshold ) for 2 values at once
const __m128d comp = _mm_cmpge_pd( r, threshold );
// Compute ( left[ i ] - drop ), for 2 values at once
const __m128d dropped = _mm_sub_pd( l, dropVec );
// Select either left or ( left - drop ) based on the comparison.
// This is the only instruction here that requires SSE 4.1.
const __m128d result = _mm_blendv_pd( l, dropped, comp );
// Store the 2 result values
_mm_storeu_pd( res, result );
}
}

如果 CPU 没有 SSE 4.1，代码将崩溃并显示"无效指令"运行时错误。为获得最佳结果，请使用 CPU ID 进行检测以正常失败。我认为现在在 2019 年假设它受到支持是相当合理的，英特尔在 2008 年这样做，AMD 在 2011 年这样做，蒸汽调查显示"96.3%"。如果您想支持较旧的 CPU，可以使用其他 3 条指令(_mm_and_pd、_mm_andnot_pd_mm_or_pd)模拟_mm_blendv_pd。

如果您可以保证数据对齐，则用_mm_load_pd替换负载会稍微快一些，_mm_cmpge_pd编译成CMPPD https://www.felixcloutier.com/x86/cmppd 可以直接从RAM获取其中一个参数。

潜在地，您可以通过编写AVX版本获得进一步的2倍改进。但我希望即使是 SSE 版本也比你的代码更快，它每次迭代处理 2 个值，并且循环内部没有条件。如果你运气不好，AVX会变慢，许多CPU需要一些时间来打开他们的AVX单元，需要数千个周期。在通电之前，AVX 代码运行速度非常慢。

您可以使用 GCC 和 Clang 的向量扩展来实现三元选择函数(参见 https://stackoverflow.com/a/48538557/2542702)。

#include <stddef.h>
#include <inttypes.h>
#if defined(__clang__)
typedef  double double4 __attribute__ ((ext_vector_type(4)));
typedef int64_t   long4 __attribute__ ((ext_vector_type(4)));
#else
typedef  double double4 __attribute__ ((vector_size (sizeof(double)*4)));
typedef int64_t   long4 __attribute__ ((vector_size (sizeof(int64_t)*4)));
#endif
double4 select(long4 s, double4 a, double4 b) {
double4 c;
#if defined(__GNUC__) && !defined(__INTEL_COMPILER) && !defined(__clang__)
c = s ? a : b;
#else
for(int i=0; i<4; i++) c[i] = s[i] ? a[i] : b[i];
#endif
return c;
}
void func(double* left, double* right, double* res, size_t size, double th, double drop) {
size_t i;
for (i = 0; i<(size&-4); i+=4) {
double4 leftv = *(double4*)&left[i];
double4 rightv = *(double4*)&right[i];
*(double4*)&res[i] = select(rightv >= th, leftv, leftv - drop);
}
for(;i<size; i++) res[i] = right[i] >= th ? left[i] : (left[i] - drop);
}

https://godbolt.org/z/h4OKMl

相关文章：