1

对于以下为蒙特卡罗模拟生成随机数的代码，我需要接收每次运行的确切总和，但这不会发生，尽管我已经修复了种子。如果有人能指出此代码的问题，我将不胜感激

#include <cmath>
#include <random>
#include <iostream>
#include <chrono>
#include <cfloat>
#include <iomanip>
#include <cstdlib>
#include <omp.h>
#include <trng/yarn2.hpp>
#include <trng/mt19937_64.hpp>
#include <trng/uniform01_dist.hpp>

using namespace std;
using namespace chrono;
const double landa = 1; 
const double exact_solution = landa / (pow(landa, 2) + 1); 
double function(double x) {
    return cos(x) / landa;
}

int main() {
    int rank; 
    const int N = 1000000;
    double sum = 0.0;
    trng::yarn2 r[6];
    for (int i = 0; i <6; i++)
    { 
        r[i].seed(0); 
    }
    
    for (int i = 0; i < 6; i++)
    {
          r[i].split(6,i);
    }
    
     trng::uniform01_dist<double> u;
     auto start = high_resolution_clock::now();
     #pragma omp parallel  num_threads(6)  
     {
         rank=omp_get_thread_num();
         #pragma omp for reduction (+: sum)
         for (int i = 0; i<N; ++i) {
            //double x = distribution(g);
      
            double x= u(r[rank]); 
            x = (-1.0 / landa) * log(1.0 - x); 
            sum = sum+function(x);
         }
    }
     double app   = sum / static_cast<double> (N);
     auto end = high_resolution_clock::now();
        
    auto diff=duration_cast<milliseconds>(end-start);
    
    cout << "Approximation is: " <<setprecision(17) << app << "\t"<<"Time: "<< setprecision(17) << diff.count()<<" Error: "<<(app-exact_solution)<< endl; 

    return 0;
}
4

1 回答 1

1

TL;DR问题有两个方面：

  1. 浮点加法不是关联的；
  2. 您正在为每个线程生成不同的随机数。

我需要收到每次运行的确切金额，但这不会发生，尽管我已经修复了种子。如果有人能指出此代码的问题，我将不胜感激

首先，您有一个竞争条件，rank=omp_get_thread_num();该变量rank在所有线程之间共享，以修复您可以rank在并行区域内声明该变量，从而使其对每个线程都是私有的。

 #pragma omp parallel  num_threads(6)  
 {
     int rank=omp_get_thread_num();
     ...
  }

在您的代码中，您不应期望sum不同数量的线程的值相同。为什么 ？

  1. 因为您正在doubles并行添加

        double sum = 0.0;
    ...
    #pragma omp for reduction (+: sum)
    for (int i = 0; i<N; ++i) {
        //double x = distribution(g);

        double x= u(r[rank]); 
        x = (-1.0 / landa) * log(1.0 - x); 
        sum = sum+function(x);
    }

    并且从每个计算机科学家应该知道的关于浮点算术的内容中可以阅读：

    另一个灰色地带涉及括号的解释。由于舍入误差，代数的结合定律不一定适用于浮点数。例如，当 x = 1e30、y = -1e30 和 z = 1 时，表达式 (x+y)+z 与 x+(y+z) 的答案完全不同（前者为 1，后者为 0 ）。

    因此，从中您得出结论，浮点加法不是关联的，以及为什么对于不同数量的线程您可能具有不同sum值的原因。

  2. 您正在为每个线程生成不同的随机值：

      for (int i = 0; i < 6; i++)
  {
      r[i].split(6,i);
  }

    因此，对于不同数量的线程，变量sum 也会得到不同的结果。

正如jérôme-richard在评论中指出的那样：

请注意，更精确的算法（如Kahan 求和）可以显着减少舍入问题，同时仍然相对较快。

于 2021-03-23T07:43:08.000 回答