问题描述
对于为蒙特卡罗模拟生成随机数的以下代码,我需要接收每次运行的确切总和,但这不会发生,尽管我已经修复了种子。如果有人能指出这段代码的问题,我将不胜感激
#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;
}
解决方法
TL;DR 问题有两方面:
- 浮点加法不具有关联性;
- 您正在为每个线程生成不同的随机数。
我需要收到每次运行的确切总和,但这不会 发生了,虽然我已经修复了种子。我会很感激,如果 任何人都可以指出这段代码的问题
首先,您在 rank=omp_get_thread_num();
上有一个 竞争条件,变量 rank
在所有线程之间共享,以解决您可以声明变量 {{1} } 在并行区域内,因此,对每个线程来说都是私有的。
rank
在您的代码中,您不应期望 #pragma omp parallel num_threads(6)
{
int rank=omp_get_thread_num();
...
}
的值对于不同的线程数会相同。为什么?
-
因为您要并行添加
sum
doubles
从 What Every Computer Scientist Should Know about Floating Point Arithmetic 可以读到:
另一个灰色区域涉及括号的解释。 由于舍入误差,代数的结合律不一定适用于浮点数。例如, 表达式 (x+y)+z 与 x+(y+z) 的答案完全不同,当 x = 1e30,y = -1e30 和 z = 1(前一种情况为 1,前一种情况为 0 后者)。
因此,您由此得出结论,浮点加法不是 关联,以及为什么对于不同数量的线程,您可能具有不同的
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); }
值。 -
您正在每个线程生成不同的随机值:
sum
因此,对于不同数量的线程,变量
for (int i = 0; i < 6; i++) { r[i].split(6,i); }
也会得到不同的结果。
正如 jérôme-richard 在评论中指出的那样:
请注意,像 Kahan summation 这样更精确的算法可以 显着减少舍入问题,同时仍然相对 快。