OpenMP中atomic互斥与reduction的区别及曼德博计算疑问

核心问题根源 #

你对reduction的理解是正确的：reduction确实是通过为每个线程创建私有副本、并行阶段仅操作私有副本、结束后合并结果的方式消除竞态，本身不需要额外原子操作，性能也优于原子操作方案。你修改后结果不一致的原因是reduction生成的私有副本仅在parallel for所属的代码块作用域内可见，你在testpoint函数中操作的仍然是全局共享的numoutside变量，没有用到reduction私有副本，去掉atomic后自然会出现多线程写竞态，导致结果错误。

OpenMP reduction的工作逻辑 #

当你添加reduction(+:numoutside)子句后，OpenMP的执行流程是：

• 并行区域开始时，为每个参与计算的线程创建一个numoutside的私有副本，初始化为加法的单位值0
• 并行循环执行过程中，线程所有对numoutside的操作，都是操作自己的私有副本，不会访问全局的共享变量
• 并行循环结束后，OpenMP将所有线程的私有副本相加，结果赋值给原来的全局numoutside变量

正确修改方式 #

你需要将计数逻辑从testpoint函数中移到parallel for的循环体内，保证操作的是reduction的私有副本，修改方案如下：

第一步：修改testpoint函数返回值 #

// 返回1表示当前点在曼德博集合外，0表示在集合内
int testpoint(struct d_complex c){
   struct d_complex z;
   int iter;
   double temp;

   z=c;
   for (iter=0; iter<MAXITER; iter++){
     temp = (z.r*z.r)-(z.i*z.i)+c.r;
     z.i = z.r*z.i*2+c.i;
     z.r = temp;
     if ((z.r*z.r+z.i*z.i)>4.0) {
       return 1;
     }
   }
   return 0;
}

第二步：调整主函数循环逻辑 #

#pragma omp parallel for default(shared) private(c,j) firstprivate(eps) reduction(+:numoutside)
for(i = 0; i<NPOINTS; i++){
    for(j=0; j < NPOINTS; j++){
         c.r = -2.0+2.5*(double)(i)/(double)(NPOINTS)+eps;
         c.i = 1.125*(double)(j)/(double)(NPOINTS)+eps;
         // 这里操作的是当前线程的reduction私有副本，无竞态无需atomic
         numoutside += testpoint(c);
    }
}

修改后即可得到正确结果，且性能比原atomic版本更高，避免了每次计数都要做原子同步的开销。

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