Openmp 多线程代码在使用多个线程时给出不同的答案

Posted 2023-02-22

【中文标题】Openmp 多线程代码在使用多个线程时给出不同的答案

【英文标题】：Openmp multithreaded code giving different answer when using multiple threads

【发布时间】：2020-10-21 14:38:46

【问题描述】：

我有以下带有openmp 多线程的蒙特卡罗代码

int main()

  std::uniform_int_distribution<int> dir(0, 1);
  std::uniform_int_distribution<int> xyz(0, 2);
  std::uniform_real_distribution<double> angle(0,360);
  mt19937 gen0;

  auto M = 20;
  long double sum = 0;
  auto num_trials = 10000;
  // omp_set_num_threads(12);

  #pragma omp parallel
  
    double loc_sum = 0.0;
    #pragma omp for
    for(int i = 0; i < num_trials; ++i)
    
      double x = 0;
      double y = 0;
      double z = 0;
      double r = 0; 
      auto N = 0;
      
      while(r < M)
      
        auto d = dir(gen);
        auto p = xyz(gen);
        if(p == 0)
        
          x += (d == 1) ? -1 : 1;
        
        else if(p == 1)
        
          y += (d == 1) ? -1 : 1;
        
        else
        
          z += (d == 1) ? -1 : 1;
        

        r = std::sqrt(x * x + y * y + z * z);
        ++N;
      

      loc_sum += N;
    

    #pragma omp critical
      sum += loc_sum;
  

变量sum 在串行和多线程中执行是完全不同的。由于对随机均匀分布的调用，我预计会略有不同，但我观察到的差异太显着了，不可能是由于随机性，我怀疑我的多线程代码中存在错误。

此代码中是否存在影响sum 的竞争条件或数据竞争？

【问题讨论】：

我很确定你不能在没有锁的情况下调用这些生成器（dir 和xyz）。您还在使用 PRNG (gen) 而不锁定。如果您使用#pragma omp critical 标记生成d 和p 的行，代码是否有效？

@Darhuuk 啊，我认为你是对的。用关键部分包裹d,p，它可以工作。我原以为生成器会以原子方式产生值。实际上，我不想在d,p 周围放置一个关键部分。这里最好的解决方案是在并行区域内创建生成器吗？

是的，在这种情况下，您最好的选择是创建生成器的副本。那些应该是相当轻量级的。您可能不想复制 PRNG 对象，因为那时您将在每个线程中拥有完全相同的状态。所以我会用proper initialization在for循环中创建它。

以正确的答案扩展了我的 cmets。

【参考方案1】：

问题是您调用生成器（dir 和 xyz）时没有锁定它们。您还使用了没有锁定的 PRNG (gen)。

这些调用都不是原子的，因为默认情况下实现它们会使单线程代码比需要的慢。

使用#pragma omp critical 标记生成d 和p 的行应该可以解决问题。

如果您不想要临界区，则需要在每个线程中分别使用 dir、xyz 和 gen 对象。可以简单地复制生成器（dir 和 xyz）。 PRNG (gen) 应该在每个线程中正确初始化，否则您最终会得到每个线程中具有完全相同状态的 PNG。例如：

std::random_device rd; /* Outside the parallel section. */

// Code below once per thread.
/* Initialization of the PRNG calls std::random_device::operator() which 
 * needs a lock around it when called in parallel. */
std::mt199937 gen;    
#pragma omp critical

  gen.seed(rd());


// for-loop starts here.

【讨论】：

rd() 调用中是否存在竞争条件？

除非明确提到有问题的函数是原子的，否则我会假设会有竞争条件。标准库中的函数极不可能是原子的，因为这会使它比不使其成为原子的要慢。 IE。通过使其具有原子性，您可以惩罚从单个线程中使用该函数的每个人。相反，如果您编写非原子函数，无论如何都可以通过在调用周围加锁来使它们成为原子函数。现在，您只需在真正需要时支付锁定/原子调用的代价。

此外，围绕生成器播种的关键部分不应对程序的运行时产生任何明显影响，因为它每个线程只运行一次。