OpenMP并行化曼德博集合C程序时函数调用失效问题求助

问题根源 #

你的代码出现异常的核心原因是全局变量引发的多线程数据竞争：

• c、z、dc被声明为全局变量，函数f(z)直接依赖全局的c进行计算
• 开启OpenMP并行后，多个线程会同时读写全局的c，导致不同线程的计算数据互相覆盖
• 直接写计算式时，代码使用的是OpenMPprivate子句创建的线程局部c副本，因此不会出现竞争；但调用函数时，函数访问的是全局c而非线程局部副本，最终导致计算错误

解决步骤 #

1. 移除全局变量，改用局部变量 #

删除代码开头的全局变量声明：

// 删掉这三行
double _Complex c;
double _Complex z;
double _Complex dc;

将这些变量改为循环内的局部变量，确保每个线程的计算数据完全独立。

2. 修改函数参数，消除全局依赖 #

修改f函数，让它通过参数接收c，不再依赖全局变量：

// 修改前
complex double f(complex double z){ return z*z*z*z*z + c;}

// 修改后
complex double f(complex double z, complex double c){ 
    return z*z*z*z*z + c;
}

d函数不需要c，保持原逻辑即可。

3. 修正函数调用与OpenMP指令 #

在main函数的循环内，声明局部变量并调用修改后的f函数：

// 内层循环内声明局部变量
double _Complex c = x + I * y;
double _Complex dc = 0;
double _Complex z = 0;

// 调用函数时传入当前线程的c
dc = d(z)*dc +1;
z = f(z, c);

同时修正OpenMP的并行指令，明确私有化局部变量：

#pragma omp parallel for private(i, c, z, dc) shared(w, h, n, r, px, r2, img)

修改后的完整代码 #

#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>

const double pi = 3.141592653589793;

complex double f(complex double z, complex double c){ 
    return z*z*z*z*z + c;
}
complex double d(complex double z) {
    return 5*z*z*z*z; 
}

double cnorm(double _Complex z) {
  return creal(z) * creal(z) + cimag(z) * cimag(z);
}

void hsv2rgb(double h, double s, double v, int *red, int *grn, int *blu) {
  double i, f, p, q, t, r, g, b;
  int ii;
  if (s == 0.0) { r = g = b = v; } else {
    h = 6 * (h - floor(h));
    ii = i = floor(h);
    f = h - i;
    p = v * (1 - s);
    q = v * (1 - (s * f));
    t = v * (1 - (s * (1 - f)));
    switch(ii) {
      case 0: r = v; g = t; b = p; break;
      case 1: r = q; g = v; b = p; break;
      case 2: r = p; g = v; b = t; break;
      case 3: r = p; g = q; b = v; break;
      case 4: r = t; g = p; b = v; break;
      default:r = v; g = p; b = q; break;
    }
  }
  *red = fmin(fmax(255 * r + 0.5, 0), 255);
  *grn = fmin(fmax(255 * g + 0.5, 0), 255);
  *blu = fmin(fmax(255 * b + 0.5, 0), 255);
}

int main()
{
  int aa = 4;
  int w = 800 * aa;
  int h = 800 * aa;
  int n = 1024;
  double r = 2;
  double px = r / (h/2);
  double r2 = 25 * 25;
  unsigned char *img = malloc(3 * w * h);
  int i,j;

  #pragma omp parallel for private(i, c, z, dc) shared(w, h, n, r, px, r2, img)
  for ( j = 0; j < h; ++j)
  {
    double y = (h/2 - (j + 0.5)) / (h/2) * r;
    for (i = 0; i < w; ++i)
    {
      double x =  (i + 0.5 - w/2) / (h/2) * r;
      double _Complex c = x + I * y;
      double _Complex dc = 0;
      double _Complex z = 0;
      int k;
      for (k = 0; k < n; ++k)
      { 
        dc = d(z)*dc +1;
        z = f(z, c);
        
        if (cnorm(z) > r2)
          break;
      }
      
      double hue = 0, sat = 0, val = 1;
      
      if (k < n) 
      {
        double _Complex de = 2 * z * log(cabs(z)) / dc;
        hue = fmod(1 + carg(de) / (2 * pi), 1);
        sat = 0.25;
        val = tanh(cabs(de) / px / aa);
      }
      
      int red, grn, blu;
      hsv2rgb(hue, sat, val, &red, &grn, &blu);
      img[3*(j * w + i)+0] = red;
      img[3*(j * w + i)+1] = grn;
      img[3*(j * w + i)+2] = blu;
    }
  }
  
  printf("P6\n%d %d\n255\n", w, h);
  fwrite(img, 3 * w * h, 1, stdout);
  free(img);
  
  return 0;
}

验证编译 #

使用原编译命令即可正常编译运行：

gcc e.c -lm -Wall -fopenmp
./a.out >ed.ppm

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