OpenMP实现的声波方程传播代码无法在2线程以上扩展性能的问题排查

问题背景 #

你在实现基于OpenMP的声波方程传播代码时，发现2线程能达到接近2倍的加速比，但进一步增加线程数（比如到4、6甚至更多），加速比始终停留在2倍左右。你的运行环境和代码细节整理如下：

系统与编译信息

• 编译器：GCC 12.2.0，编译标志：-O3 -funroll-loops -ftree-vectorize -fopenmp
• CPU：AMD Ryzen 5 5600G（6核12线程）
• OpenMP版本：4.5

关键变量参数

• 模型尺寸：扩展后nz=2101、nx=2601（原模型nz=501、nx=1001，每边扩展800边界）
• 缓存相关：CACHE_SZ=16（对应64字节缓存行，16个float），nbzcache=2101/16≈131
• 其他：NDERIV=3，nbx=2601

核心代码片段

#pragma omp parallel default(none)                                             \
    shared(p1, p2, wavemovie, isrc, srcindex, ishot, nbz, idxs, deriv, ricker, \
               dx, velextend, nbzcache, nbx, i0_extended_model, i0_model,      \
               gamma_x, gamma_z, nx, nxx, nzz, nz, nt, ndt, nshots, nthreads)
{
    int iconv, ii, ix, it = 0;
    int mm;
    int imodel;
    int iz, niz;
    float gama, mgamma, invpgamma;
    float laplacian;
    float tempfield[CACHE_SZ];
    nthreads = omp_get_num_threads();
#pragma omp single
    {
        printf("Number of threads used:%d\n", nthreads);
    }
    while (ishot < nshots) {
#pragma omp single
        {
            it = 0;
            srcindex = isrc + idxs * ishot;
        }
#pragma omp for
        for (int iswap = 0; iswap < nxx * nzz; iswap++) {
            p1[iswap] = 0.0;
            p2[iswap] = 0.0;
        }

        //  Forward field time update
        for (it = 0; it < nt; it++) {
#pragma omp for schedule(static)
            for (mm = 0; mm < nbzcache * nbx; mm++) {
                ix = (int)mm / nbzcache;
                niz = mm % nbzcache;
                iz = niz * CACHE_SZ;
                // Loop to write in cachelines
                for (ii = 0; ii < CACHE_SZ; ii++) {
                    imodel = i0_extended_model + ix * nzz + iz + ii;
                    gama = gamma_x[ix + NDERIV - 1] +
                           gamma_z[iz + NDERIV - 1 + ii];
                    mgamma = -(1 - gama) / (1 + gama);
                    invpgamma = 0.5 * (1 - mgamma);
                    laplacian = 2 * deriv[0] * p2[imodel];
                    for (iconv = 1; iconv < NDERIV; iconv++) {
                        laplacian += deriv[iconv] * p2[imodel - iconv];
                        laplacian += deriv[iconv] * p2[imodel + iconv];
                        laplacian += deriv[iconv] * p2[imodel - nzz * iconv];
                        laplacian += deriv[iconv] * p2[imodel + nzz * iconv];
                    }
                    tempfield[ii] = invpgamma * velextend[imodel] * laplacian;
                    tempfield[ii] += invpgamma * 2.0 * p2[imodel];
                    tempfield[ii] += mgamma * p1[imodel];
                }
                for (ii = 0; ii < CACHE_SZ; ii++) {
                    p1[imodel - (CACHE_SZ - 1) + ii] = tempfield[ii];
                }
            }
#pragma omp single
            {
                p1[srcindex] +=
                    -pow(dx, 2.0) * ricker[it] * velextend[srcindex];
                if (it % ndt == 0) {
                    for (ii = 0; ii < nx; ii++) {
                        fwrite(p1 + i0_model + ii * nzz, sizeof(float), nz,
                               wavemovie);
                    }
                }
                swap(&p1, &p2);
            }
        }
#pragma omp single
        {
            ishot++;
            printf("ishot:%d\n", ishot);
            printf("nshots:%d\n", nshots);
        }
    }
}

性能无法扩展的核心原因分析 #

结合你的代码和硬件环境，我认为主要瓶颈来自串行区域的开销占比过高，其次是缓存局部性与循环划分的优化不足，具体如下：

1. 时间步循环内的串行瓶颈（最关键）

你在每个时间步it的末尾都加了一个#pragma omp single块，里面包含三个操作：源点注入、写wavemovie文件、交换p1/p2指针。其中：

• 写wavemovie的循环是完全串行的：for (ii = 0; ii < nx; ii++) { fwrite(...) }，这个循环的规模是nx=1001次迭代，每次写nz=501个float。如果nt很大（比如几千上万步），这个串行操作的总开销会被放大，当线程数越多，所有线程等待单个线程完成串行任务的时间占比就越高，直接限制了加速比的上限。
• single块本身会触发线程同步，增加额外的等待开销，进一步压缩并行部分的收益。

2. 缓存局部性与循环划分的问题

你的mm循环采用了schedule(static)划分，但循环的迭代方式可能没有最大化缓存利用率：

• mm循环的迭代逻辑是先遍历完一个ix对应的所有nbzcache个niz，再进入下一个ix。这种划分下，不同线程会处理连续的ix范围，但ix*nzz是跨越大块内存的跨步（nzz=2101，每个跨步是2101个float≈8KB），导致p2的邻域访问（imodel ± nzz*iconv）会频繁触发缓存失效，尤其是当线程数增加时，每个核心的L1/L2缓存被多个线程共享，缓存命中率进一步下降，内存访问延迟成为新的瓶颈。
• 你用tempfield按缓存行写入p1的逻辑是对的，但代码中p1[imodel - (CACHE_SZ - 1) + ii]的写法有点绕（其实等价于i0_extended_model + ix*nzz + iz + ii），虽然功能没问题，但可能会让编译器的优化器难以进一步优化内存访问模式。

3. 超线程的边际效益递减

AMD Ryzen 5 5600G是6核12线程的CPU，超线程（SMT）的核心是让同一物理核心的两个线程共享L1/L2缓存。当线程数超过6（物理核心数）时，新增的超线程线程会加剧缓存竞争，性能提升非常有限甚至下降，这也会让你觉得“线程数增加后加速比不涨”。

针对性的优化建议 #

1. 并行化串行瓶颈区域

• 并行化wavemovie写入循环：这个循环的迭代之间没有依赖，完全可以并行化。修改成：

  if (it % ndt == 0) {
      #pragma omp for
      for (ii = 0; ii < nx; ii++) {
          fwrite(p1 + i0_model + ii * nzz, sizeof(float), nz, wavemovie);
      }
  }
  

  注：GCC的fwrite是线程安全的，若担心IO锁开销，也可以让每个线程写入独立临时文件，最后合并。
• 简化ishot循环的single块：把ishot的递增用原子操作替代，打印逻辑只让主线程执行，减少同步开销：

  // 替换原有ishot的single块
  #pragma omp atomic
  ishot++;
  if (omp_get_thread_num() == 0) {
      printf("ishot:%d\n", ishot);
      printf("nshots:%d\n", nshots);
  }
  

2. 优化循环划分与缓存局部性

• 调整mm循环的调度策略：用schedule(static, nbzcache)让每个线程连续处理一个完整的ix对应的所有迭代，提升ix维度的数据局部性：

  #pragma omp for schedule(static, nbzcache)
  for (mm = 0; mm < nbzcache * nbx; mm++) {
      // ... 原有逻辑
  }
  

• 重构循环嵌套顺序：显式展开ix和niz循环，用collapse(2)让OpenMP更均衡地分配任务：

  #pragma omp for schedule(static) collapse(2)
  for (int ix = 0; ix < nbx; ix++) {
      for (int niz = 0; niz < nbzcache; niz++) {
          int iz = niz * CACHE_SZ;
          // 原有ii循环的逻辑
      }
  }
  

3. 绑定线程到物理核心

通过设置OpenMP环境变量，让线程固定绑定到物理核心，避免调度带来的缓存失效：

export OMP_PROC_BIND=true
export OMP_PLACES=cores

4. 用性能工具定位热点

建议用perf工具分析程序的性能瓶颈：

# 记录程序的性能数据
perf record -g ./your_program
# 查看性能报告
perf report

通过报告你可以直观看到串行区域的时间占比、缓存失效占比等关键指标，精准定位优化点。

备注：内容来源于stack exchange，提问作者Joan Guastalla