OpenMP并行C代码运行速度与串行持平甚至更慢的解决方案咨询

• 高频原子操作是首要性能损耗源：当前代码每个满足截断条件的粒子对都要触发4次原子操作（3个力分量更新+1个能量更新），原子操作的锁竞争开销远高于多线程带来的计算收益，2线程场景下锁开销甚至会完全抵消并行收益，导致速度不如串行版本。
• 循环索引i错误声明为共享变量：OpenMP for循环的索引默认应为私有属性，手动将i设置为shared会引发多线程同时修改i的 data race，不仅会导致性能异常，还会直接造成计算结果错误。
• 默认调度策略负载严重不均：N体模拟的外层i循环，每个i对应的内层j循环计算量为nbodies-i次，默认静态调度会将前半段计算量更大的i分给同一个线程，后半段计算量小的i分给另一个线程，线程负载差可达数倍，并行效率极低。
• 存在伪共享问题：多线程同时修改forces数组相邻位置的元素时，会频繁触发缓存行失效，额外增加内存访问开销。

核心优化：消除原子操作，改用线程私有副本累加 #

每个线程单独维护私有力数组和能量副本，计算完成后再统一合并到全局变量，全程避免原子操作，仅在合并阶段加一次临界区锁，锁开销可忽略不计。修改后代码参考如下：

pos = calloc(nbodies, sizeof(*pos));
forces = calloc(nbodies, sizeof(*forces));
//...more...
printf("Calculating......\n");
ene = 0.0;

#pragma omp parallel default(none) shared(pos, forces, cut2, nbodies, ene) private(i,j,k,d,d2,d3,rij)
{
    // 声明线程私有副本
    double (*local_forces)[3] = calloc(nbodies, sizeof(*local_forces));
    double local_ene = 0.0;

    // 改用dynamic调度解决负载不均，块大小16可根据nbodies规模调整为32/64
    #pragma omp for schedule(dynamic, 16)
    for(i=0; i<nbodies; ++i){
        for(j=i+1; j<nbodies; ++j) {
            d2 = 0.0;
            for(k=0; k<3; ++k) {
                rij[k] = pos[i][k] - pos[j][k];
                d2 += rij[k]*rij[k];
            }
            if (d2 <= cut2) {
               d = sqrt(d2);
               d3 = d*d2;
               for(k=0; k<3; ++k) {
                    double f = -rij[k]/d3;
                    // 直接更新私有副本，无需原子操作
                    local_forces[i][k] += f;
                    local_forces[j][k] -= f;
               }
               local_ene += -1.0/d; 
           }
        }
    }

    // 临界区合并所有线程的私有结果到全局变量
    #pragma omp critical
    {
        ene += local_ene;
        for(i=0; i<nbodies; ++i) {
            for(k=0; k<3; ++k) {
                forces[i][k] += local_forces[i][k];
            }
        }
    }
    free(local_forces);
}

配套优化 #

• 调整DevCpp编译参数：开启O2优化，编译指令添加-fopenmp -O2 -march=native，禁止在Debug模式下测试性能，Debug模式会完全禁用OpenMP相关优化。
• 确认计算规模：如果nbodies < 1000，并行本身的开销就会大于收益，只有当粒子数足够大时并行优化才会有明显效果。

内容的提问来源于stack exchange，提问作者D0rkDevelop4r