OMP_NUM_THREADS=1 时 #pragma ompatomic 的性能问题

我观察到我正在编写的 openmp 代码出现了意外的（对我来说！）行为。代码结构如下：

#pragma omp parallel for
for(int i=0;i<N;i++){ 
 // lots of calculations that produce 3 integers i1,i2,i3 and 3 doubles d1,d2,d3 
 #pragma omp atomic 
 J1[i1] += d1;
 #pragma omp atomic
 J2[i2] += d2; 
 #pragma omp atomic
 J3[i3] += d3; 
}

我编译了该代码的三个不同版本：

1) 使用 openmp (-fopenmp)

2）没有openmp

3）使用openmp，但没有3个原子操作（只是作为测试，因为原子操作是必要的）

当我使用环境变量 OMP_NUM_THREADS=1 运行版本 1) 时，我观察到版本 2) 的速度显着下降；而版本 3) 的运行速度与版本 2) 一样快。

我想知道这种行为的原因（为什么即使是单线程，原子操作也会减慢代码速度？！）以及是否可以以版本 1）运行速度一样快的方式编译/重写代码版本 2)。

我在问题末尾附上一个显示上述行为的工作示例。我编译了 1)：

g++ -fopenmp -o toy_code toy_code.cpp -std=c++11 -O3

2) with:

g++ -o toy_code_NO_OMP toy_code.cpp -std=c++11 -O3

3) 与：

g++ -fopenmp -o toy_code_NO_ATOMIC toy_code_NO_ATOMIC.cpp -std=c++11 -O3

编译器的版本是gcc版本5.3.1 20160519（Debian 5.3.1-20）。 3个版本的执行时间分别为：

1) 1 分 24 秒

2) 51 秒

3）51秒

预先感谢您的任何建议！

// toy_code.cpp 
#include <stdio.h>
#include <iostream>
#include <stdlib.h>
#include <cmath>
#include <omp.h>
#define Np 1000000
#define N 1000

int main (){
        double* Xp, *Yp, *J,*Jb;
        Xp = new double[Np];
        Yp = new double[Np];  
        J = new double [N*N];
        Jb = new double [N*N];

        for(int i=0;i<N*N;i++){
            J[i]=0.0;
            Jb[i]=0.0;
        }

        for(int i=0;i<Np;i++){
            Xp[i] = rand()*1.0/RAND_MAX - 0.5;
            Yp[i] = rand()*1.0/RAND_MAX - 0.5;
        }

        for(int n=0; n<2000; n++){
        #pragma omp parallel for
        for(int p=0;p<Np;p++){
            double rx = (Xp[p]+0.5)*(N-1);
            double ry = (Yp[p]+0.5)*(N-1);
            int xindex = (int)floor(rx+0.5);
            int yindex = (int)floor(ry+0.5);
            int k;
            k=xindex*N+yindex;

            #pragma omp atomic
            J[k]+=1;
            #pragma omp atomic
            Jb[k]+=1;
         }
         }

        delete[] Xp;
        delete[] Yp;
        delete[] J;
        delete[] Jb;

return 0;
}

如果启用 OpenMP，gcc 必须生成适用于仅在运行时已知的任意数量线程的不同代码。

在这种特殊情况下，看一下输出gcc -S（通过标签稍微缩短）。

没有 OpenMP：

.loc 1 38 0 discriminator 2  # Line 38 is J[k]+=1;
movsd   8(%rsp), %xmm1
cvttsd2si   %xmm0, %edx
cvttsd2si   %xmm1, %eax
movsd   .LC3(%rip), %xmm0
imull   $1000, %eax, %eax
addl    %edx, %eax
cltq
salq    $3, %rax
leaq    0(%r13,%rax), %rdx
.loc 1 40 0 discriminator 2   # Line 40 is Jb[k]+=1;
addq    %r12, %rax
.loc 1 29 0 discriminator 2
cmpq    $8000000, %r15
.loc 1 38 0 discriminator 2
addsd   (%rdx), %xmm0
movsd   %xmm0, (%rdx)
.loc 1 40 0 discriminator 2
movsd   .LC3(%rip), %xmm0
addsd   (%rax), %xmm0
movsd   %xmm0, (%rax)

循环展开使得这变得相当复杂。

With -fopenmp:

movsd   (%rsp), %xmm2
cvttsd2si   %xmm0, %eax
cvttsd2si   %xmm2, %ecx
imull   $1000, %ecx, %ecx
addl    %eax, %ecx
movslq  %ecx, %rcx
salq    $3, %rcx
movq    %rcx, %rsi
addq    16(%rbp), %rsi
movq    (%rsi), %rdx
movsd   8(%rsp), %xmm1
jmp .L4
movq    %rax, %rdx
movq    %rdx, (%rsp)
movq    %rdx, %rax
movsd   (%rsp), %xmm3
addsd   %xmm1, %xmm3
movq    %xmm3, %rdi
lock cmpxchgq   %rdi, (%rsi)
cmpq    %rax, %rdx
jne .L9
.loc 1 40 0
addq    24(%rbp), %rcx
movq    (%rcx), %rdx
jmp .L5
.p2align 4,,10
.p2align 3
movq    %rax, %rdx
movq    %rdx, (%rsp)
movq    %rdx, %rax
movsd   (%rsp), %xmm4
addsd   %xmm1, %xmm4
movq    %xmm4, %rsi
lock cmpxchgq   %rsi, (%rcx)
cmpq    %rax, %rdx
jne .L10
addq    $8, %r12
cmpq    %r12, %rbx
jne .L6

我不会尝试解释或理解这里发生的所有细节，但这对于消息来说不是必需的：编译器必须使用可能更昂贵的不同原子指令，尤其是lock cmpxchgq.

除了这个基本问题之外，OpenMP 可能会以任何可以想象的方式干扰优化器，例如干扰展开。我还看到了一个奇怪的案例，英特尔编译器实际上为 OpenMP 循环生成了更高效的串行代码。

附：认为自己很幸运——情况可能会更糟。如果编译器无法将原子指令映射到硬件指令，则必须使用锁，这会更慢。