并行化牛顿运动方程求解循环后执行时间增加的问题排查

以下是求解牛顿运动方程、计算粒子下一时刻位置与速度的串行代码片段：

• deltat_s和dtSqr为常量
• X[]、Z[]、V_X[]、V_Z[]为存储位置与速度的双精度数组
• a_x[]、a_z[]存储加速度分量
• ComputeForces()用于计算力并更新a_x[]和a_z[]，第一个循环需更新V_X[]、V_Z[]供其使用

for(i=0; i<N_TOTAL; i++)
  {
    X[i] = X[i] + deltat_s * V_X[i] + 0.5 * dtSqr * a_x[i];
    V_X[i] = V_X[i] + 0.5 * deltat_s * a_x[i];
    Z[i] = Z[i] + deltat_s * V_Z[i] + 0.5 * dtSqr * a_z[i];
    V_Z[i] = V_Z[i] + 0.5 * deltat_s * a_z[i];
  }

ComputeForces ();

for(i=0; i<N_TOTAL; i++)
  {
    V_X[i] = V_X[i] + 0.5 * deltat_s * a_x[i];
    V_Z[i] = V_Z[i] + 0.5 * deltat_s * a_z[i];
  }

我使用OpenMP对代码进行了并行化：

#pragma omp parallel for ordered schedule (guided, 4) default(shared) private(i)
for(i=0; i<N_TOTAL; i++)
      {
        X[i] = X[i] + deltat_s * V_X[i] + 0.5 * dtSqr * a_x[i];
        V_X[i] = V_X[i] + 0.5 * deltat_s * a_x[i];
        Z[i] = Z[i] + deltat_s * V_Z[i] + 0.5 * dtSqr * a_z[i];
        V_Z[i] = V_Z[i] + 0.5 * deltat_s * a_z[i];
      }

    ComputeForces ();

    #pragma omp parallel for schedule (guided, 16) default(shared) private(i)
    for(i=0; i<N_TOTAL; i++)
      {
        V_X[i] = V_X[i] + 0.5 * deltat_s * a_x[i];
        V_Z[i] = V_Z[i] + 0.5 * deltat_s * a_z[i];
      }

但并行化后代码执行时间反而比串行版本更长：

• 单次执行：串行耗时3.546000e-02，并行耗时4.632966e-02
• 30000步迭代后：串行耗时1.102543e+03，并行耗时1.363923e+03

我怀疑V_X[]和V_Z[]存在伪共享，尝试调整分块大小和调度方式后仍未解决，寻求优化方案。

硬件环境：Intel(R) Xeon(R) CPU E5-2630 v4 @ 2.20GHz（40核），N_TOTAL=1000

1. 聚焦计算密集型部分并行，避免小循环并行开销 #

当前N_TOTAL=1000，两个循环的单迭代计算量极小，OpenMP线程创建、调度的开销远超过并行带来的收益。40核调度1000个迭代，每个线程仅处理25个左右元素，调度成本占比过高。

优化动作：

• 取消两个小循环的并行化，仅对ComputeForces()内部的计算逻辑并行（如果该函数是计算密集型，比如包含粒子间成对力计算）。
• 若ComputeForces()已实现并行，两个小循环直接串行执行即可，其计算量远低于线程调度开销。

2. 解决伪共享：调整数据布局 #

当前采用的数组结构（SoA）中，V_X[i]、V_Z[i]分散在不同数组，双精度元素占8字节，64字节缓存行可容纳8个同数组元素，多线程修改同一缓存行内不同元素会导致缓存行频繁失效。

优化动作：

• 切换为结构体数组（AoS），将单个粒子的所有数据打包，确保每个粒子的数据落在同一缓存行：

typedef struct {
    double x, z;
    double vx, vz;
    double ax, az;
} Particle;

Particle particles[N_TOTAL];

线程处理连续粒子时，只会加载对应粒子的缓存行，避免跨线程缓存冲突。

• 若坚持使用SoA，可在数组元素间填充空白，确保每个线程处理的元素落在独立缓存行：

#define PAD 8
double V_X[N_TOTAL * PAD];
double V_Z[N_TOTAL * PAD];
// 访问时用V_X[i*PAD]、V_Z[i*PAD]，间隔64字节避免伪共享

3. 优化并行调度策略 #

40核线程数过多，对1000个元素的循环而言，线程切换成本极高。

优化动作：

• 设置线程数为物理核心数（E5-2630 v4为10核，可通过omp_set_num_threads(10)或环境变量OMP_NUM_THREADS=10配置），避免超线程带来的额外开销。
• 对循环使用static调度，替代guided调度。当前每个迭代计算量完全相同，静态调度能最小化调度开销：

#pragma omp parallel for schedule(static) private(i)

4. 开启编译器最高级别优化 #

编译时启用-O3优化和OpenMP支持，编译器会自动进行循环展开、向量优化等操作，提升串行与并行代码的执行效率：

gcc -O3 -fopenmp your_code.c -o your_program

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