使用OpenMP并行化光线追踪循环后输出存在细微差异的问题

问题描述 #

我尝试使用OpenMP对简单光线追踪程序的内层for循环进行并行化，程序运行看似正常，生成的文件与串行输出几乎一致，但始终存在一处固定位置的细微差异。

并行化的循环代码如下：

// Parallelize the inner loop
#pragma omp parallel for private(x, pixel, xx, yy, raydir, ncalls) reduction(+:ncalls_line)
for (x = 0; x < width; ++x) {
  pixel = &image[y * width + x];
  xx = (2 * ((x + 0.5) * invWidth) - 1) * angle * aspectratio;
  yy = (1 - 2 * ((y + 0.5) * invHeight)) * angle;
  Vec3_new(&raydir, xx, yy, -1);
  Vec3_normalize(&raydir);
  ncalls = trace(pixel, &origin, &raydir, size, spheres, 0);
  ncalls_line += ncalls;
}  // End of parallel loop

使用Linux的cmp命令对比串行与并行生成的二进制PPM文件时，发现差异固定在第4行第17字节处。我曾尝试将ncalls_line的reduction改为使用#pragma omp critical指令，但问题依旧，差异位置完全相同。

完整的render函数源代码：

void render(const unsigned size, const Sphere spheres[size])
{
  //sequential start time
  double start = omp_get_wtime();

  //  unsigned width = 640, height = 480;
  unsigned width = 1024, height = 768;
  Vec3 *image, *pixel, origin, raydir;
  Real invWidth = 1 / (Real)width, invHeight = 1 / (Real)height;
  Real fov = 30, aspectratio = width / (Real)height;
  Real xx, yy, angle = tan(M_PI * 0.5 * fov / 180.0);
  unsigned x, y, i;
  int ncalls;
  double ncalls_line, max_ncalls_line, min_ncalls_line;
  int histo[HIST_N_INTV];
  FILE *F;

  image = malloc(width * height * sizeof(*image));
  Vec3_new0(&origin);
  init_histogram(histo);
  max_ncalls_line = 0;
  min_ncalls_line = INFINITY;

  // Trace rays
  for (y = 0; y < height; ++y) {
    ncalls_line = 0;
    for (x = 0; x < width; ++x) {
      pixel = &image[y*width+x];
      xx = (2 * ((x + 0.5) * invWidth) - 1) * angle * aspectratio;
      yy = (1 - 2 * ((y + 0.5) * invHeight)) * angle;
      Vec3_new(&raydir, xx, yy, -1);
      Vec3_normalize(&raydir);
      ncalls = trace(pixel, &origin, &raydir, size, spheres, 0);
      ncalls_line += ncalls;
    }
    ncalls_line = ncalls_line/width;
    if (ncalls_line < min_ncalls_line)
      min_ncalls_line = ncalls_line;
    if (ncalls_line > max_ncalls_line)
      max_ncalls_line = ncalls_line;
    update_histogram(histo, ncalls_line);
  }

  //ending time before writing to console and to file
  double end = omp_get_wtime();
  double execTime = end - start;
  printf("Sequential execution time: %f", execTime);

  print_histogram(histo);
  printf("Least complex line had %6.3f calls to trace per pixel.\n", min_ncalls_line);
  printf("Most  complex line had %6.3f calls to trace per pixel.\n", max_ncalls_line);

  // Save result to a PPM image (keep these flags if you compile under Windows)
  F = fopen("./sequential.ppm", "wb");
  fprintf(F, "P6\n%u %u\n255\n", width, height);
  for (i = 0; i < width * height; ++i) {
    fprintf(F, "%c%c%c",
        (unsigned char)(min(1.0, image[i].x) * 255),
        (unsigned char)(min(1.0, image[i].y) * 255),
        (unsigned char)(min(1.0, image[i].z) * 255));
  }
  fclose(F);
  free(image);
} 

问题分析与解决办法 #

核心原因 #

差异源于浮点数运算的非确定性：

• OpenMP并行时，不同线程可能使用不同的CPU寄存器或运算单元，浮点数的舍入顺序、精度表现会产生细微差异（比如Vec3_normalize中的平方根、除法运算，或是trace函数里的光线相交计算）。
• 固定位置的差异说明该像素的计算对浮点数精度变化特别敏感——比如刚好处于两种舍入结果的边界，乘以255后整数部分差1，最终体现在PPM文件的字节差异上。

具体排查与修复步骤 #

1. 定位差异像素

   • PPM文件第4行开始是像素数据，第17字节对应第(17-1)/3=5个像素（从0计数），即坐标x=5, y=0（第一行对应y=0）。
   • 单独计算该像素的串行与并行结果，对比Vec3的x/y/z分量，确定是哪个颜色通道出现差异。

2. 强制浮点数运算一致性
   编译时添加以下选项，减少线程间的浮点精度差异：

   • GCC/Clang：添加-ffloat-store（禁止寄存器缓存浮点值，强制写回内存）和-frounding-math（启用严格舍入模式）。
   • Intel ICC：添加-fp-model strict。
   • 确保所有浮点运算统一使用相同精度（比如全程用double，避免float与double混合导致的精度波动）。

3. 确认trace函数的线程安全性

   • 检查trace函数是否存在全局变量、静态变量或共享内存的写入操作。如果有，必须用#pragma omp critical或线程私有变量保护。
   • 从现有代码看，trace的参数均为传入值，且spheres是const常量，大概率线程安全，但仍需确认内部是否有隐藏的共享状态。

4. 验证PPM写入阶段的一致性

   • 现有代码中PPM写入是串行执行的，所以问题不在写入阶段，但可以先打印image数组中差异像素的分量值，对比串行与并行版本，确认差异是出现在计算阶段还是写入阶段。

5. 排查边界计算逻辑
   如果强制浮点一致性后差异仍存在，重点检查该像素对应的光线计算逻辑：比如光线是否刚好擦过球体表面，串行与并行时的相交判断结果是否有差异，这类边界情况容易因浮点精度变化产生不同的计算分支。

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