在SIMD循环中打开MP:SIMD兼容功能?

是,启用了优化( -O3 -march=native )现代编译器可以通过函数指针可靠地内联 如果 满足这些条件:

• 函数指针具有编译时常量值
• 它指向的是一个函数,编译器可以看到它的定义

这听起来很简单,但 如果此代码在UNIX/Linux上的共享lib中使用(使用 -fPIC ) 那么符号插入规则意味着 float halve(float x) { return x * 0.5f; } ^一 即使在同一个翻译单元中也不能内联。见 Sorry state of dynamic libraries on Linux 。

使用 inline 关键字,即使在构建共享库时也允许内联 像 static 如果编译器决定在每个调用站点内联,它就允许编译器根本不发出函数的独立定义。

使用 内联的 在 halve , square 和 apply_simd . (因为 应用模拟人生 需要与通过 减半 作为函数arg。独立定义 应用模拟人生 是无用的,因为它不能内联未知函数。)如果它们位于 .cpp 而不是 .h 你也可以做它们 静止的 或者相反,否则就把它们做成 内联的 .

一次完成尽可能多的工作

我怀疑你想写一些非常低效的东西,比如:

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(或add)吞吐量(每个时钟2个256位矢量)是它们的存储带宽(每个时钟到l1d 1个256位矢量)的两倍。 计算强度很重要。不要通过编写多个执行单独的琐碎操作的循环来给代码添加gimp。

在任何循环展开时,或者如果数据不适合L1D,则为SIMD的吞吐量 x[i] = 0.25f * y[i]*y[i] 将与单独的这些操作相同。

我检查了g++8.2和clang++6.0的ASM输出。 on the Godbolt compiler explorer . 即使有 __restrict 为了说明x和y不重叠,编译器仍然做了两个单独的循环。

将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的一个bug),这目前不可用。

但只使用一个函数指针,它实际上并不能帮助捕获忘记使用 内联的 关键字或其他方式破坏了编译器的内联能力。它只意味着函数地址必须是动态链接时间常量(即,直到动态库的运行时绑定时才知道),因此它不会将您从符号插入中保存下来。AT 编译 时间与 -FPIC ,编译器仍然不知道它所能看到的全局函数的版本是否是链接时实际解析的版本,或者 LD_PRELOAD 或者主可执行文件中的符号将覆盖它。所以它只发出从got加载函数指针的代码,并在循环中调用它。西姆德当然不可能。

但是,它通过传递函数指针的方式来阻止你用脚射自己,这种方式并不总是内联的。可能和 constexpr 不过,在模板中使用之前,您仍然可以将它们作为参数传递。 所以,如果不是GCC的bug阻止了你将它与lambda一起使用的话,你可能会想使用它。

C++ 17允许不带捕获功能的自动存储LAMBDA作为函数对象。(以前的标准要求外部或内部( 静止的 )作为模板参数传递的函数的链接。)

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有一个bug( 83258 )用这种方式使用lambda。 见 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; 而不是做 减半 首先,看看我们得到了什么样的坏代码。这就是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 哈斯韦尔/天空湖。所以这个循环实际上是7个Uops,每次迭代最多运行7/4个周期。

根据英特尔的优化手册,对于宽向量,接近每时钟2读1写的限制显然是一个持续运行的问题。(他们说,Skylake每时钟可以支持82个字节,而不是一个时钟中存储的96个以上的峰值。)如果数据不知道是对齐的,这就特别不明智了,而且GCC8已经转向了一种未知对齐数据的乐观策略:使用未对齐的加载/存储,并让硬件处理CA当没有32字节对齐时。GCC7和早期版本将指针对准主循环之前,并且只加载一次向量。

x / 2. 进入之内 x * 0.5f ,避免晋升到 double .

不用乘法而不用除法 -ffast-math 因为 0.5 完全可以表示为 float 与分母不是2的幂的分数不同。

但请注意 0.5 * x 做 不 优化到 0.5f * x ;GCC和Clang实际上扩展到 双重的 然后回来。我不确定这是否是一个遗漏的优化vs. x/2。 或者,如果存在真正的语义差异,使得它无法将双常量优化为一个浮点,而此时它可以精确表示为 浮动 。