c++ - Open MP:SIMD循环中的SIMD兼容功能？

Question

通常，我可能会编写一个类似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){ ... } 

改变什么？还是仅当我仅根据本机函数/运算符明确写下该操作时，才能期望从SIMD中受益？

Answer 1

是的，启用优化(-O3 -march=native)，如果满足以下条件，现代编译器可以可靠地内联函数指针:

函数指针具有一个编译时常量值

它指向一个函数，编译器可以看到该函数的定义

听起来很容易确定，但是如果在Unix / Linux上的共享库(使用 -fPIC编译)中使用此代码，则使用 ，然后使用符号插入规则意味着 float halve(float x) { return x * 0.5f; } 1即使在同一翻译单元内也不能内联。参见 Sorry state of dynamic libraries on Linux。
即使构建共享库，也可以使用 inline关键字进行内联；像 static一样，如果编译器决定在每个调用站点内联，则它根本不发出该函数的独立定义。
在 inline， halve和 square上使用 apply_simd。 (因为 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的FMA / mul(或增加)吞吐量(每个时钟2x 256位 vector )是其存储带宽的两倍(每个时钟1x 256位 vector 至L1d)。 计算强度很重要。不要通过编写执行单独的琐碎操作的多个循环来简化代码
如果展开任何循环，或者如果数据不适合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位 vector ，乘以2。 (非索引寻址模式可以通过展开隐藏两个指针增量来实现。傻傻的编译器。:/)

Lambda或函数指针作为模板参数
使用本地lambda作为模板参数(在函数内部定义)，编译器肯定可以始终内联。但是(由于gcc错误)目前尚无法使用。
但是仅使用一个函数指针，它实际上并不能帮助您捕获忘记使用 inline关键字或破坏编译器内联功能的情况。这仅意味着函数地址必须是动态链接时间常数(即直到动态库的运行时绑定(bind)才知道)，因此它不会使您免于插入符号的麻烦。在使用 -fPIC进行编译时，编译器仍然不知道它可以看到的全局函数的版本是否是链接时实际解析的版本，或者 LD_PRELOAD或主可执行文件中的符号是否会覆盖它。因此，它仅发出从GOT加载函数指针的代码，并在循环中调用它。 SIMD当然是不可能的。
但是，通过以并非总是内联的方式传递函数指针，确实可以阻止您脚踏实地。不过，也许可以使用 constexpr将它们作为args传递出去，然后再在模板中使用。 因此，如果不是因为gcc错误导致您无法将其与lambda结合使用，则可能要使用它。
C++ 17允许传递没有捕获的自动存储lambda作为函数对象。 (先前的标准要求以外部或内部( static)链接的形式作为模板参数传递的函数。)

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在使用lambda时存在一个错误(83258)。 有关详细信息，请参见 Can I use the result of a C++17 captureless lambda constexpr conversion operator as a function pointer template non-type argument?。
不过，我们可以在模板中使用常规功能。

// `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未使用索引寻址模式(Haswell / Skylake上的 stops it from micro-fusing)，则两次重载相同的 vector 以保存指令将是一个好主意。因此，此循环实际上发出的是7微分，每次迭代的运行最佳为7/4个周期。
根据英特尔的优化手册，展开时，对于宽 vector ，每个时钟限制接近2读取+ 1写入限制显然是持续运行的问题。 (他们说Skylake可能每个时钟维持82个字节，而不是一个时钟中存储96个加载+的峰值。)如果不知道数据是否对齐，并且gcc8转向针对未知的乐观策略，这是不明智的，对齐数据:使用未对齐的加载/存储，并让硬件处理没有32字节对齐的情况。 gcc7和更早版本在主循环之前对齐指针，并且仅将 vector 加载一次。

脚注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常量优化为float。