用于浮动阈值操作的 SIMD
SIMD for float threshold operation
我想更快地进行一些矢量计算，我相信用于浮点数比较和操作的 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) ;
}
}
主要是，如果 right 值高于 threshold，它会将 left 值降低 drop。
大小在128-256左右（不算大），但是计算调用量很大
我尝试从循环展开开始，但没有赢得很多性能，但可能需要一些编译指令。
您能否建议对代码进行一些改进以加快计算速度？
给你（未经测试），我试图在评论中解释他们的作用。
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 年支持，Steam 调查显示“96.3%”。如果你想支持旧的 CPUs，可以用其他 3 条指令模拟 _mm_blendv_pd，_mm_and_pd、_mm_andnot_pd、_mm_or_pd.
如果能保证数据对齐，用_mm_load_pd替换load会稍微快一些，_mm_cmpge_pd编译成CMPPD https://www.felixcloutier.com/x86/cmppd可以直接从RAM中取一个参数.
有可能，您可以通过编写 AVX 版本进一步提高 2 倍。但我希望即使是 SSE 版本也比你的代码更快，它每次迭代处理 2 个值，并且循环内没有条件。如果你运气不好，AVX 会更慢，许多 CPUs 需要一些时间来启动他们的 AVX 单元，需要数千个周期。在通电之前，AVX 代码运行非常缓慢。
Clang 已经 auto-vectorizes 这几乎是 Soonts 建议手动执行的方式。 在你的指针上使用 __restrict 所以它不需要后备版本适用于某些阵列之间的重叠。它仍然是 auto-vectorizes，但它使函数膨胀。
不幸的是，gcc 仅 auto-vectorizes 和 -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) 比较，并将生成的掩码用于 merge-masked 减法。
谓词为真的元素将获得 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 也可以auto-vectorize
它对所有编译器都更有效，因为你只需要一个 AND，而不是 variable-blend。 所以只用 SSE2 就明显更好，而且在大多数情况下也更好CPU 即使支持 SSE4.1 blendvpd，因为该指令效率不高。
您可以根据比较结果将left[i]减去0.0或drop。
根据比较结果生成 0.0 或常量非常有效：只需一条 andps 指令。 （0.0 的 bit-pattern 是 all-zeros，并且 SIMD 比较生成全 1 或全 0 位的向量。因此 AND 保留旧值或将其归零。）
我们也可以加上 -drop 而不是减去 drop。这需要对输入进行额外的取反，但 AVX 允许 memory-source 操作数用于 vaddpd。不过，GCC 选择使用索引寻址模式，因此实际上并不能帮助减少 Intel CPU 上的 front-end uop 计数；它会 "unlaminate"。但即使使用 -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 link：
# 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 也是 compile-time 内联后的常量，并且无论如何都会被加载。或者特别是如果 -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
您可以使用 GCC 和 Clang 的向量扩展来实现三元 select 函数（参见 ）。
#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);
}
我想更快地进行一些矢量计算，我相信用于浮点数比较和操作的 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) ;
}
}
主要是，如果 right 值高于 threshold，它会将 left 值降低 drop。
大小在128-256左右（不算大），但是计算调用量很大
我尝试从循环展开开始，但没有赢得很多性能，但可能需要一些编译指令。
您能否建议对代码进行一些改进以加快计算速度？
给你（未经测试），我试图在评论中解释他们的作用。
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 年支持，Steam 调查显示“96.3%”。如果你想支持旧的 CPUs，可以用其他 3 条指令模拟 _mm_blendv_pd，_mm_and_pd、_mm_andnot_pd、_mm_or_pd.
如果能保证数据对齐，用_mm_load_pd替换load会稍微快一些，_mm_cmpge_pd编译成CMPPD https://www.felixcloutier.com/x86/cmppd可以直接从RAM中取一个参数.
有可能，您可以通过编写 AVX 版本进一步提高 2 倍。但我希望即使是 SSE 版本也比你的代码更快，它每次迭代处理 2 个值，并且循环内没有条件。如果你运气不好，AVX 会更慢，许多 CPUs 需要一些时间来启动他们的 AVX 单元，需要数千个周期。在通电之前，AVX 代码运行非常缓慢。
Clang 已经 auto-vectorizes 这几乎是 Soonts 建议手动执行的方式。 在你的指针上使用 __restrict 所以它不需要后备版本适用于某些阵列之间的重叠。它仍然是 auto-vectorizes，但它使函数膨胀。
不幸的是，gcc 仅 auto-vectorizes 和 -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) 比较，并将生成的掩码用于 merge-masked 减法。
谓词为真的元素将获得 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 也可以auto-vectorize
它对所有编译器都更有效，因为你只需要一个 AND，而不是 variable-blend。 所以只用 SSE2 就明显更好，而且在大多数情况下也更好CPU 即使支持 SSE4.1 blendvpd，因为该指令效率不高。
您可以根据比较结果将left[i]减去0.0或drop。
根据比较结果生成 0.0 或常量非常有效：只需一条 andps 指令。 （0.0 的 bit-pattern 是 all-zeros，并且 SIMD 比较生成全 1 或全 0 位的向量。因此 AND 保留旧值或将其归零。）
我们也可以加上 -drop 而不是减去 drop。这需要对输入进行额外的取反，但 AVX 允许 memory-source 操作数用于 vaddpd。不过，GCC 选择使用索引寻址模式，因此实际上并不能帮助减少 Intel CPU 上的 front-end uop 计数；它会 "unlaminate"。但即使使用 -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 link：
# 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 也是 compile-time 内联后的常量，并且无论如何都会被加载。或者特别是如果 -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
您可以使用 GCC 和 Clang 的向量扩展来实现三元 select 函数（参见 ）。
#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);
}