打开 MP：SIMD 循环中的 SIMD 兼容功能？
Open MP: SIMD compatible function in SIMD loops?
通常，我可能会编写一个 simd 循环，例如：
float * x = (float *) malloc(10 * sizeof(float));
float * y = (float *) malloc(10 * sizeof(float));
for(int i = 0; i < 10; i++)
y[i] = 10;
#pragma omp simd
for(int i = 0; i < 10; i++)
x[i] = y[i]*y[i];
假设我有两个任务：
float square(float x) {
return x * x;
}
float halve(float x) {
return x / 2.;
}
还有一个 omp 循环原语：
void apply_simd(float * x, float * y, int length, float (*simd_func)(float c)){
#pragma omp simd
for(int i = 0; i < length; i++)
x[i] = simd_func(y[i])
}
这在 SIMD 的参数范围内合法吗？或者编译器是否会生成比我显式内联所有内容效率更低的代码？
是否写作：
float inline square(float x){ ... }
改变什么？或者只有当我明确地写下本机 functions/operators 方面的操作时，我才能期望从 SIMD 中受益？
是的，启用优化 (-O3 -march=native)，现代编译器可以通过函数指针可靠地内联 如果满足这些条件：
- 函数指针有一个编译时常量值
- 它指向一个编译器可以看到其定义的函数
这听起来很容易确保，但是 如果此代码用于 Unix/Linux 上的共享库（使用 -fPIC 编译），则符号插入规则意味着 float halve(float x) { return x * 0.5f; }1 即使在同一个翻译单元中也不能内联。参见 Sorry state of dynamic libraries on Linux。
使用 inline 关键字允许内联，即使在构建共享库时也是如此；像 static 它允许编译器根本不发出函数的独立定义，如果它决定在每个调用站点内联。
在 halve、square 和 apply_simd 上使用 inline。 （因为 apply_simd 需要内联到将 halve 作为函数 arg 传递的调用方。 apply_simd 的独立定义是无用的，因为它不能内联未知函数。）如果它们在 .cpp 而不是 .h 中，你也可以将它们设为 static 或者改为 inline.
一次完成尽可能多的工作
我怀疑你想像这样写一些非常低效的东西：
apply_simd(x, y, length, halve); // copy y to x
apply_simd(x, x, length, square); // then update x in-place
// NEVER DO THIS, make one function that does both things
// with gcc and clang, compiles as written to two separate loops.
只进行复制和乘以 0.5f 的循环通常会成为内存带宽的瓶颈。像 Haswell（或 Skylake）这样的现代 CPU 的吞吐量（每个时钟 2x 256 位向量）是它们存储带宽（每个时钟 1x 256 位向量到 L1d）的两倍 FMA/mul（或增加）。 计算强度很重要。不要通过编写多个循环来执行单独的琐碎操作来削弱您的代码
随着任何循环的展开，或者如果数据不适合 L1d，SIMD x[i] = 0.25f * y[i]*y[i] 的吞吐量将与这些操作中的任何一个单独的吞吐量相同。
我检查了 g++ 8.2 和 clang++ 6.0 on the Godbolt compiler explorer 的 asm 输出。即使 __restrict 告诉它 x 和 y 不重叠，编译器仍然进行了 2 个独立的循环。
将 lambda 作为函数指针传递
我们可以使用 lambda 轻松地将任意操作组合到单个函数中，并将其作为函数指针传递。这解决了上述创建两个单独循环的问题，同时仍然为您提供了将循环包装在一个函数中的所需语法。
如果您的 halve(float) 函数是一些不平凡的东西的占位符，您可以在 lambda 中使用它来将它与其他东西组合起来。例如square(halve(a))
在早期的 C++ 标准中，您需要将 lambda 赋值给函数指针。 (Lambda as function parameter)
// your original function mostly unchanged, but with size_t and inline
inline // allows inlining even with -fPIC
void apply_simd(float * x, const float *y, size_t length, float (*simd_func)(float c)){
#pragma omp simd
for(size_t i = 0; i < length; i++)
x[i] = simd_func(y[i]);
}
C++11 调用者：
// __restrict isn't needed with OpenMP, but you might want to assert non-overlapping for better auto-vectorization with non-OpenMP compilers.
void test_lambda(float *__restrict x, const float *__restrict y, size_t length)
{
float (*funcptr)(float) = [](float a) -> float {
float h=0.5f*a; // halve first allows vmulps with a memory source operand
return h*h; // 0.25 * a * a doesn't optimize to that with clang :/
};
apply_simd(x, y, length, funcptr);
}
在 C++17 中，它更容易，并且只适用于文字匿名 lambda：
void test_lambda17(float *__restrict x, const float *__restrict y, size_t length)
{
apply_simd(x, y, length, [](float a) {
float h = 0.5f*a;
return h * h;
}
);
}
这些都可以使用 gcc 和 clang 高效地编译成这样的内部循环 (Godbolt compiler explorer).
.L4:
vmulps ymm0, ymm1, YMMWORD PTR [rsi+rax]
vmulps ymm0, ymm0, ymm0
vmovups YMMWORD PTR [rdi+rax], ymm0
add rax, 32
cmp rax, rcx
jne .L4
clang 展开了一些，并且可能接近一个 256 位矢量加载+每个时钟存储，有 2 次乘法。 （非索引寻址模式可以通过展开隐藏两个指针增量来实现。愚蠢的编译器。:/)
作为模板参数的 Lambda 或函数指针
使用本地 lambda 作为模板参数（在函数内部定义），编译器绝对可以始终内联。但是（由于 gcc 错误）这目前不可用。
但是只有一个函数指针，它实际上并不能帮助捕获您忘记使用 inline 关键字或以其他方式破坏编译器内联能力的情况。这仅意味着函数地址必须是动态 link 时间常量（即直到动态库的运行时绑定才知道），因此它不会使您免于符号插入。在 编译 时 -fPIC，编译器仍然不知道它能看到的全局函数的版本是否是在 [=168= 时实际解析的版本] 时间，或者如果 LD_PRELOAD 或主可执行文件中的符号将覆盖它。所以它只是发出从 GOT 加载函数指针的代码，并在循环中调用它。 SIMD当然不可能。
不过，它确实可以防止您以不能总是内联的方式传递函数指针而搬起石头砸自己的脚。不过，在使用 constexpr 的情况下，您仍然可以在模板中使用之前将它们作为参数传递。 因此，如果不是因为 gcc 错误阻止您将它与 lambda 一起使用，您可能想使用它。
C++17 允许将没有捕获的自动存储 lambda 作为函数对象传递。 （以前的标准要求外部或内部 (static) link作为模板参数传递的函数的年龄。）
template <float simd_func(float c)>
void apply_template(float *x, const float *y, size_t length){
#pragma omp simd
for(size_t i = 0; i < length; i++)
x[i] = simd_func(y[i]);
}
void test_lambda(float *__restrict x, const float *__restrict y, size_t length)
{
// static // even static doesn't help work around the gcc bug
constexpr auto my_op = [](float a) -> float {
float h=0.5f*a; // halve first allows vmulps with a memory source operand
return h*h; // 0.25 * a * a doesn't optimize to that with clang :/
};
// I don't know what the unary + operator is doing here, but some examples use it
apply_lambda<+my_op>(x, y, length); // clang accepts this, gcc doesn't
}
clang 编译得很好，但 g++ 错误地拒绝了它，即使 -std=gnu++17
不幸的是，gcc 有一个错误（83258) with using lambdas this way. See 了解详情。
不过，我们可以在模板中使用常规函数。
// `inline` is still necessary for it to actually inline with -fPIC (in a shared lib)
inline float my_func(float a) { return 0.25f * a*a;}
void test_template(float *__restrict x, const float *__restrict y, size_t length)
{
apply_lambda<my_func>(x, y, length); // not actually a lambda, just a function
}
然后我们从 g++8.2 -O3 -fopenmp -march=haswell 得到一个像这样的内部循环。请注意，我首先使用 0.25f * a * a; 而不是 halve 来查看我们得到了什么样的错误代码。这就是 g++8.2 所做的。
.L25:
vmulps ymm0, ymm1, YMMWORD PTR [rsi+rax] # ymm0 = 0.25f * y[i+0..7]
vmulps ymm0, ymm0, YMMWORD PTR [rsi+rax] # reload the same vector again
vmovups YMMWORD PTR [rdi+rax], ymm0 # store to x[i+0..7]
add rax, 32
cmp rax, rcx
jne .L25
如果 gcc 没有使用索引寻址模式（stops it from micro-fusing 在 Haswell/Skylake 上，重新加载同一个向量两次以保存指令可能是个好主意。所以这个循环实际上发出 7 微指令，运行 每次迭代最多 7/4 个周期。
根据 Intel 的优化手册，随着展开，接近宽向量每个时钟 2 次读取 + 1 次写入的限制显然是持续运行的问题。 （他们说 Skylake 可能每个时钟维持 82 个字节，而峰值为 96 个加载 + 存储在一个时钟中。）如果不知道数据是否对齐，这尤其不明智，并且 gcc8 已切换到未知的乐观策略 -对齐数据：使用未对齐的 loads/stores 并让硬件处理没有 32 字节对齐的情况。 gcc7 及更早版本在主循环之前对齐指针，并且只加载向量一次。
脚注 1：幸运的是 gcc 和 clang 可以将 x / 2. 优化为 x * 0.5f，避免提升为 double.
不用 -ffast-math 就可以使用乘法而不是除法，因为 0.5f 可以精确表示为 float，这与分母不是 2 的幂的分数不同。
但请注意 0.5 * x 不会 优化为 0.5f * x； gcc 和 clang 实际上确实扩展到 double 并返回。我不确定这是否是与 x / 2. 相比的优化遗漏，或者是否存在真正的语义差异阻止它在可以精确表示为 float 时将 double 常量优化为浮点数。
通常，我可能会编写一个 simd 循环，例如：
float * x = (float *) malloc(10 * sizeof(float));
float * y = (float *) malloc(10 * sizeof(float));
for(int i = 0; i < 10; i++)
y[i] = 10;
#pragma omp simd
for(int i = 0; i < 10; i++)
x[i] = y[i]*y[i];
假设我有两个任务：
float square(float x) {
return x * x;
}
float halve(float x) {
return x / 2.;
}
还有一个 omp 循环原语：
void apply_simd(float * x, float * y, int length, float (*simd_func)(float c)){
#pragma omp simd
for(int i = 0; i < length; i++)
x[i] = simd_func(y[i])
}
这在 SIMD 的参数范围内合法吗？或者编译器是否会生成比我显式内联所有内容效率更低的代码？
是否写作：
float inline square(float x){ ... }
改变什么？或者只有当我明确地写下本机 functions/operators 方面的操作时，我才能期望从 SIMD 中受益？
是的，启用优化 (-O3 -march=native)，现代编译器可以通过函数指针可靠地内联 如果满足这些条件：
- 函数指针有一个编译时常量值
- 它指向一个编译器可以看到其定义的函数
这听起来很容易确保，但是 如果此代码用于 Unix/Linux 上的共享库（使用 -fPIC 编译），则符号插入规则意味着 float halve(float x) { return x * 0.5f; }1 即使在同一个翻译单元中也不能内联。参见 Sorry state of dynamic libraries on Linux。
使用 inline 关键字允许内联，即使在构建共享库时也是如此；像 static 它允许编译器根本不发出函数的独立定义，如果它决定在每个调用站点内联。
在 halve、square 和 apply_simd 上使用 inline。 （因为 apply_simd 需要内联到将 halve 作为函数 arg 传递的调用方。 apply_simd 的独立定义是无用的，因为它不能内联未知函数。）如果它们在 .cpp 而不是 .h 中，你也可以将它们设为 static 或者改为 inline.
一次完成尽可能多的工作
我怀疑你想像这样写一些非常低效的东西：
apply_simd(x, y, length, halve); // copy y to x
apply_simd(x, x, length, square); // then update x in-place
// NEVER DO THIS, make one function that does both things
// with gcc and clang, compiles as written to two separate loops.
只进行复制和乘以 0.5f 的循环通常会成为内存带宽的瓶颈。像 Haswell（或 Skylake）这样的现代 CPU 的吞吐量（每个时钟 2x 256 位向量）是它们存储带宽（每个时钟 1x 256 位向量到 L1d）的两倍 FMA/mul（或增加）。 计算强度很重要。不要通过编写多个循环来执行单独的琐碎操作来削弱您的代码
随着任何循环的展开，或者如果数据不适合 L1d，SIMD x[i] = 0.25f * y[i]*y[i] 的吞吐量将与这些操作中的任何一个单独的吞吐量相同。
我检查了 g++ 8.2 和 clang++ 6.0 on the Godbolt compiler explorer 的 asm 输出。即使 __restrict 告诉它 x 和 y 不重叠，编译器仍然进行了 2 个独立的循环。
将 lambda 作为函数指针传递
我们可以使用 lambda 轻松地将任意操作组合到单个函数中，并将其作为函数指针传递。这解决了上述创建两个单独循环的问题，同时仍然为您提供了将循环包装在一个函数中的所需语法。
如果您的 halve(float) 函数是一些不平凡的东西的占位符，您可以在 lambda 中使用它来将它与其他东西组合起来。例如square(halve(a))
在早期的 C++ 标准中，您需要将 lambda 赋值给函数指针。 (Lambda as function parameter)
// your original function mostly unchanged, but with size_t and inline
inline // allows inlining even with -fPIC
void apply_simd(float * x, const float *y, size_t length, float (*simd_func)(float c)){
#pragma omp simd
for(size_t i = 0; i < length; i++)
x[i] = simd_func(y[i]);
}
C++11 调用者：
// __restrict isn't needed with OpenMP, but you might want to assert non-overlapping for better auto-vectorization with non-OpenMP compilers.
void test_lambda(float *__restrict x, const float *__restrict y, size_t length)
{
float (*funcptr)(float) = [](float a) -> float {
float h=0.5f*a; // halve first allows vmulps with a memory source operand
return h*h; // 0.25 * a * a doesn't optimize to that with clang :/
};
apply_simd(x, y, length, funcptr);
}
在 C++17 中，它更容易，并且只适用于文字匿名 lambda：
void test_lambda17(float *__restrict x, const float *__restrict y, size_t length)
{
apply_simd(x, y, length, [](float a) {
float h = 0.5f*a;
return h * h;
}
);
}
这些都可以使用 gcc 和 clang 高效地编译成这样的内部循环 (Godbolt compiler explorer).
.L4:
vmulps ymm0, ymm1, YMMWORD PTR [rsi+rax]
vmulps ymm0, ymm0, ymm0
vmovups YMMWORD PTR [rdi+rax], ymm0
add rax, 32
cmp rax, rcx
jne .L4
clang 展开了一些，并且可能接近一个 256 位矢量加载+每个时钟存储，有 2 次乘法。 （非索引寻址模式可以通过展开隐藏两个指针增量来实现。愚蠢的编译器。:/)
作为模板参数的 Lambda 或函数指针
使用本地 lambda 作为模板参数（在函数内部定义），编译器绝对可以始终内联。但是（由于 gcc 错误）这目前不可用。
但是只有一个函数指针，它实际上并不能帮助捕获您忘记使用 inline 关键字或以其他方式破坏编译器内联能力的情况。这仅意味着函数地址必须是动态 link 时间常量（即直到动态库的运行时绑定才知道），因此它不会使您免于符号插入。在 编译 时 -fPIC，编译器仍然不知道它能看到的全局函数的版本是否是在 [=168= 时实际解析的版本] 时间，或者如果 LD_PRELOAD 或主可执行文件中的符号将覆盖它。所以它只是发出从 GOT 加载函数指针的代码，并在循环中调用它。 SIMD当然不可能。
不过，它确实可以防止您以不能总是内联的方式传递函数指针而搬起石头砸自己的脚。不过，在使用 constexpr 的情况下，您仍然可以在模板中使用之前将它们作为参数传递。 因此，如果不是因为 gcc 错误阻止您将它与 lambda 一起使用，您可能想使用它。
C++17 允许将没有捕获的自动存储 lambda 作为函数对象传递。 （以前的标准要求外部或内部 (static) link作为模板参数传递的函数的年龄。）
template <float simd_func(float c)>
void apply_template(float *x, const float *y, size_t length){
#pragma omp simd
for(size_t i = 0; i < length; i++)
x[i] = simd_func(y[i]);
}
void test_lambda(float *__restrict x, const float *__restrict y, size_t length)
{
// static // even static doesn't help work around the gcc bug
constexpr auto my_op = [](float a) -> float {
float h=0.5f*a; // halve first allows vmulps with a memory source operand
return h*h; // 0.25 * a * a doesn't optimize to that with clang :/
};
// I don't know what the unary + operator is doing here, but some examples use it
apply_lambda<+my_op>(x, y, length); // clang accepts this, gcc doesn't
}
clang 编译得很好，但 g++ 错误地拒绝了它，即使 -std=gnu++17
不幸的是，gcc 有一个错误（83258) with using lambdas this way. See
不过，我们可以在模板中使用常规函数。
// `inline` is still necessary for it to actually inline with -fPIC (in a shared lib)
inline float my_func(float a) { return 0.25f * a*a;}
void test_template(float *__restrict x, const float *__restrict y, size_t length)
{
apply_lambda<my_func>(x, y, length); // not actually a lambda, just a function
}
然后我们从 g++8.2 -O3 -fopenmp -march=haswell 得到一个像这样的内部循环。请注意，我首先使用 0.25f * a * a; 而不是 halve 来查看我们得到了什么样的错误代码。这就是 g++8.2 所做的。
.L25:
vmulps ymm0, ymm1, YMMWORD PTR [rsi+rax] # ymm0 = 0.25f * y[i+0..7]
vmulps ymm0, ymm0, YMMWORD PTR [rsi+rax] # reload the same vector again
vmovups YMMWORD PTR [rdi+rax], ymm0 # store to x[i+0..7]
add rax, 32
cmp rax, rcx
jne .L25
如果 gcc 没有使用索引寻址模式（stops it from micro-fusing 在 Haswell/Skylake 上，重新加载同一个向量两次以保存指令可能是个好主意。所以这个循环实际上发出 7 微指令，运行 每次迭代最多 7/4 个周期。
根据 Intel 的优化手册，随着展开，接近宽向量每个时钟 2 次读取 + 1 次写入的限制显然是持续运行的问题。 （他们说 Skylake 可能每个时钟维持 82 个字节，而峰值为 96 个加载 + 存储在一个时钟中。）如果不知道数据是否对齐，这尤其不明智，并且 gcc8 已切换到未知的乐观策略 -对齐数据：使用未对齐的 loads/stores 并让硬件处理没有 32 字节对齐的情况。 gcc7 及更早版本在主循环之前对齐指针，并且只加载向量一次。
脚注 1：幸运的是 gcc 和 clang 可以将 x / 2. 优化为 x * 0.5f，避免提升为 double.
不用 -ffast-math 就可以使用乘法而不是除法，因为 0.5f 可以精确表示为 float，这与分母不是 2 的幂的分数不同。
但请注意 0.5 * x 不会 优化为 0.5f * x； gcc 和 clang 实际上确实扩展到 double 并返回。我不确定这是否是与 x / 2. 相比的优化遗漏，或者是否存在真正的语义差异阻止它在可以精确表示为 float 时将 double 常量优化为浮点数。