c语言 - 为什么我没有使用从工作代码复制的并行化 (OpenMP) 算法获得相同的 PI 估计?

Question

下面的代码是YouTube视频的直接翻译，该视频使用OpenMP和Monte Carlo估算PI。即使使用相同的输入，我也不会在这里获得它们的输出。事实上，似乎大约一半的价值是我得到的。

int main() {
int num; // number of iterations
printf("Enter number of iterations you want the loop to run for:   ");
scanf_s("%d", &num);
double x, y, z, pi;
long long int i;
int count = 0;
int num_thread;
printf("Enter number of threads you want to run to parallelize the process:t");
scanf_s("%d", &num_thread);
printf("n");
#pragma omp parallel firstprivate(x,y,z,i) shared(count) num_threads(num_thread)
{
srand((int)time(NULL) ^ omp_get_thread_num());
for (i = 0; i < num; i++) {
x = (double)rand() / (double)RAND_MAX;
y = (double)rand() / (double)RAND_MAX;
z = pow(((x * x) + (y * y)), .5);
if (z <= 1) {
count++;
}
}
} // END PRAGMA
pi = ((double)count / (double)(num * num_thread)) * 4;
printf("The value of pi obtained is %fn", pi);
return 0;
}

我还直接使用了橡树岭国家实验室网站(https://www.olcf.ornl.gov/tutorials/monte-carlo-pi/)的类似算法：

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

int main(int argc, char* argv[])
{
int niter = 1000000;            //number of iterations per FOR loop
double x,y;                     //x,y value for the random coordinate
int i;                          //loop counter
int count=0;                //Count holds all the number of how many good coordinates
double z;                       //Used to check if x^2+y^2<=1
double pi;                      //holds approx value of pi
int numthreads = 16;

#pragma omp parallel firstprivate(x, y, z, i) shared(count) num_threads(numthreads)
{
srandom((int)time(NULL) ^ omp_get_thread_num());    //Give random() a seed value
for (i=0; i<niter; ++i)              //main loop
{
x = (double)random()/RAND_MAX;      //gets a random x coordinate
y = (double)random()/RAND_MAX;      //gets a random y coordinate
z = sqrt((x*x)+(y*y));          //Checks to see if number is inside unit circle
if (z<=1)
{
++count;            //if it is, consider it a valid random point
}
}
//print the value of each thread/rank
}
pi = ((double)count/(double)(niter*numthreads))*4.0;
printf("Pi: %fn", pi);

return 0;
}

而且我有确切的问题，所以我认为这不是代码，而是我的机器。

我正在运行VS Studio 22，Windows 11，具有16核i9-12900kf和32 gb内存。

编辑：我忘了提到我确实更改了第二种算法以使用 srand() 和 rand() 代替。

Answer 1

代码中有许多错误：

1. 正如@JeromeRichard和@JohnBollingerrandsrandrandom所指出的，不是线程安全的，您应该使用线程安全解决方案。

2. 第++count;行存在争用条件(不同的线程读取和写入共享变量)。您应该使用缩减来避免它。

3. 该代码假定您使用numthreads线程，但 OpenMP 不保证您实际获得了请求的所有线程。我认为如果您因此获得了 PI/2，问题应该是请求的线程数和获得的线程数之间的差异。如果在循环之前使用#pragma omp parallel for...，则不需要对线程数进行任何假设(即在这种情况下，计算 PI 的公式不包含线程数)。

4. 一个小评论是，您不需要使用耗时的pow功能。

把它放在一起，你的代码应该是这样的：

#pragma omp parallel for reduction(+:count) num_threads(num_thread)    
for (long long int i = 0; i < num; i++) {
const double x = threadsafe_random_number_between_0_1();
const double y = threadsafe_random_number_between_0_1();
const double z = x * x + y * y;
if (z <= 1) {
count++;
}
}
double pi = ((double) count / (double) num ) * 4.0;

Answer 2

一个假设，但我可能是错的：你随时间随机初始化，所以它可能发生在不同的线程使用相同的时间，这可能会导致生成相同的随机数，所以结果会非常糟糕，因为你多次得到相同的值。这是蒙特卡洛方法的一个问题，其中 2 个相同的点将产生错误的结果。