当前位置：首页 > article >正文

(done) openMP学习 (Day13: 线程私有数据和如何支持库(Pi again)，蒙特卡洛计算 Pi，线性同余法)

article 2025/2/10 13:06:22

url: https://dazuozcy.github.io/posts/introdution-to-openmp-intel/#23-%E5%8F%AF%E6%80%95%E7%9A%84%E4%B8%9C%E8%A5%BF%E5%86%85%E5%AD%98%E6%A8%A1%E5%9E%8Batomicsflushpairwise%E5%90%8C%E6%AD%A5%20

视频：https://www.bilibili.com/video/BV1SW411s7ST?spm_id_from=333.788.videopod.episodes&vd_source=7a1a0bc74158c6993c7355c5490fc600&p=26

先来回顾 threadprivate

threadprivate: makes global data private to a thread.

例子：

int counter = 0;
#pragma omp threadprivate(counter) //每个线程都有一个counter变量，并且被初始化为0

int increment_counter()
{
    counter++;
    return (counter);
}

一个例子：蒙特卡洛法计算pi值

#include <random>

static long num_trails = 10000;

void calc_pi_montacalo()
{
    long i;
    long Ncircle = 0;
    double pi, x, y;
    double r = 1.0;

    std::random_device rd;
    std::mt19937 rng(rd());
    std::uniform_real_distribution<> dist(-1.0, 1.0);
    
    #pragma omp parallel for private(x, y) reduction(+:Ncircle)
    for(i = 0; i < num_trails; i++) {
        x = dist(rng);
        y = dist(rng);
        if ((x * x + y * y) <= r * r) {
            Ncircle++;
        }
    }

    pi = 4.0 * ((double)Ncircle / (double)num_trails);
    printf("\n%d trails, pi is %f\n", num_trails, pi);
}

蒙特卡洛算法中要使用随机数生成器。这里出现了一个问题：随机数生成器也许不支持并行，从而导致上述代码无法并行。

线性同余法随机数生成器

在这里插入图片描述

源码：

/* simple single thread version */
static long MULTIPLIER = 1366;
static long ADDEND = 150889;
static long PMOD = 714025;
long random_last = 0;

double random()
{
    long random_next;
    random_next = (MULTIPLIER* random_last + ADDEND) % PMOD;
    random_last = random_next;

    return ((double)random_next / (double)PMOD);
}

经过测试：并行时，同样参数生成的伪随机数并不同
在这里插入图片描述

原因是，上述代码中实际上存在数据竞争：

/* simple single thread version */
static long MULTIPLIER = 1366;
static long ADDEND = 150889;
static long PMOD = 714025;
long random_last = 0;

double random()
{
    long random_next;
    random_next = (MULTIPLIER* random_last + ADDEND) % PMOD;
    random_last = random_next;

    return ((double)random_next / (double)PMOD);
}