OpenMP、MPI、CUDA 并行蒙特卡洛方法计算圆周率

用 OpenMP、MPI、CUDA 实现蒙特卡洛算法（Monte Carlo method）计算圆周率 $\pi$，并探索算法参数与运算性能以及精度之间的关系。

硬件环境：

• Intel(R) Core(TM) i7-7700HQ CPU @ 2.80GHz
• RAM: 16GB
• GeForce GTX 1050 Ti

可直接编译的完整代码见文末附录。

串行实现

编程

蒙特卡洛方法，又称“随机抽样”或“统计试验”方法，它的原理是通过大量的重复随机试验来对复杂数学系统的仿真。

计算圆周率是蒙特卡洛方法的一个常见应用：先画一个正方形及其内切四分之一圆，然后随机地向正方形内投点，然后统计落在四分之一圆内的点的数量（红色）占点的总数量的比例。当总数量很大时，该比例将接近 $\dfrac{\pi}{4}$。这就是蒙特卡洛方法计算圆周率的原理。

我们先用串行方法实现该算法，作为对照实验。

长整型变量total为点的总数量，长整型变量count保存落在四分之一圆内的点的数量。我们可以使用(double)rand() / RAND_MAX来生成 $[0, 1)$ 内的随机数。串行蒙特卡洛方法的 C 语言代码如下（下面仅为关键部分，完整代码见1_mc_pi_serial.c）：

srand(time(NULL));      // 随机初始化seed
long long total = 1e6;  // 默认100万个样本
long long count = 0;
double x, y;
for(int i = 0; i < total; i++) {
    x = (double)rand() / RAND_MAX;
    y = (double)rand() / RAND_MAX;
    if(x*x + y*y <= 1) {
        count++;
    }
}
double pi = 4 * (double)count / total;

最后，依次打印出total的值、count的值、蒙特卡洛估计的 $\pi$ 值，以及估计与真值之间的误差。注意，$\pi$ 的真值可由acos(-1)给出。

printf("[+] total = %lld\n", total);
printf("[+] count = %lld\n", count);
printf("[+] pi    = %f\n", pi);
printf("[+] loss  = %e\n", acos(-1) - pi);

运行结果

用以下命令编译、运行并计时。注意argv[1]为样本总数，这里取 $1\times 10^8$。

gcc 1_mc_pi_serial.c -o 1_mc_pi_serial && time ./1_mc_pi_serial 100000000

结果如下。以下结果仅供观赏，后面我们还会做详细的结果比较和分析。

[+] total = 100000000
[+] count = 78540789
[+] pi    = 3.141632
[+] loss  = -3.890641e-05

real    0m2.069s
user    0m2.031s
sys     0m0.016s

本文地址：https://www.jeddd.com/article/parallel-monte-carlo-method-estimate-pi.html

OpenMP 并行实现

编程

rand()函数不是线程安全的，若直接对串行程序的 for 循环用 OpenMP 并行化，则结果反而会比串行的更慢。作为替代，我们使用rand_r()，并在每个线程内初始化用于线程内部的随机数种子。

OpenMP 并行实现的蒙特卡洛方法如下（下面仅为关键部分，完整代码见2_mc_pi_openmp.c）：

long long total = 1e6;  // 默认100万个样本
long long count = 0;
double x, y;
#pragma omp parallel num_threads(tn)
{
    unsigned seed = time(NULL);
    
    #pragma omp for private(x, y) reduction(+:count)
    for(int i = 0; i < total; i++) {
        x = (double)rand_r(&seed) / RAND_MAX;
        y = (double)rand_r(&seed) / RAND_MAX;
        if(x*x + y*y <= 1) {
            count++;
        }
    }
}
double pi = 4 * (double)count / total;

以上代码中，并没有直接使用#pragma omp parallel for，而是首先在#pragma omp parallel中用线程编号初始化了随机数种子，其次在并行化 for 循环。for 循环中的 x 和 y 应该是私有变量，因此使用private。最后，我们要的count是各个线程之和，因此使用reduction(+:count)进行规约。

运行结果

编译、运行并计时，采用 4 个线程进行 $1\times 10^8$ 次试验。别忘了使用-fopenmp参数：

gcc -fopenmp 2_mc_pi_openmp.c -o 2_mc_pi_openmp && time ./2_mc_pi_openmp 100000000 4

结果：

[+] total = 100000000
[+] count = 78537369
[+] pi    = 3.141495
[+] loss  = 9.789359e-05

num_threads = 4

real    0m0.471s
user    0m1.844s
sys     0m0.031s

速度比串行实现快了许多。

MPI 并行实现

编程

MPI 采用消息传递方式实现多进程并行。与串行实现或 OpenMP 不同，在编写 MPI 程序时需要思考进程之间的通信时机与作用，因此要稍稍复杂一些。（但也不会太复杂，因为蒙特卡洛方法计算圆周率本身是一个较为简单的应用。）

先初始化 MPI：

int world_size, my_rank;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &world_size);
MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);

在秩为 0 的进程（我们把它称为 0 号进程）中初始化total表示样本总数，默认为 100 万，当然也可以从命令行参数中获取该值。然后，用total除以进程数量world_size获得每个进程的样本数量，保存在 0 号进程的local_total中。注意用户需要考虑是否能整除的问题。

0 号进程负责把local_total广播（MPI_Bcast()）给所有其他进程，然后每个进程内部只需做local_total次实验，并使用local_count统计本进程的结果。最后，把每个进程的local_count规约（MPI_Reduce）求和，得到最终的结果。计算 $\pi$ 的方法与之前相同，不再赘述。

代码（下面仅为关键部分，完整代码见3_mc_pi_mpi.c）：

long long total;
long long local_total;
long long count;
long long local_count = 0;
double x, y;

if(my_rank == 0) {
    total = 1e6;  // 默认100万个样本
    if(argc >= 2) {
        total = atoi(argv[1]);
    }
    // 在root中计算local_total
    local_total = total / world_size; 
}
// 把local_total广播给所有进程
MPI_Bcast(&local_total, 1, MPI_LONG_LONG, 0, MPI_COMM_WORLD);

for(int i = 0; i < local_total; i++) {
    x = (double)rand() / RAND_MAX;
    y = (double)rand() / RAND_MAX;
    if(x*x + y*y <= 1) {
        local_count++;
    }
}
MPI_Reduce(&local_count, &count, 1, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
double pi = 4 * (double)count / total;

最后释放 MPI 资源：

MPI_Finalize();

运行结果

编译、运行并计时，采用 4 个节点（1 台主机的 4 个进程）进行 $1\times 10^8$ 次试验。

mpicc 3_mc_pi_mpi.c -o 3_mc_pi_mpi && time mpiexec -n 4 ./3_mc_pi_mpi 100000000

结果：

[+] total = 100000000
[+] count = 78535076
[+] pi    = 3.141403
[+] loss  = 1.896136e-04

Processes = 4

real    0m0.702s
user    0m2.578s
sys     0m0.109s

CUDA 并行实现

编程

可以将大量的蒙特卡洛试验放到 GPU 上进行，每个 thread 进行一次试验。默认有 128 个 thread，blocks 的数量由试验总数和 threads 数量确定。

核函数trial内只进行一次试验，并将结果保存在布尔数组count_d中。参数中的x_d和y_d是提前生成的随机数数组。

__global__ void trial(int seed, bool count_d[], double x_d[], double y_d[]) {
    long long id = blockIdx.x * blockDim.x + threadIdx.x;
    double x = x_d[id], y = y_d[id];
    if(x*x + y*y <= 1) {
        count_d[id] = true;
    }
    else {
        count_d[id] = false;
    }
}

在主函数中，我们首先创建相关变量，申请内存和显存空间，生成随机数数组，将相关数据拷贝到显存。

dim3 threads(tn);
dim3 blocks((total+tn-1) / tn);
long long real_total = threads.x * blocks.x;

bool* count_h = new bool[real_total];
bool* count_d;
double* x_h = new double[real_total];
double* y_h = new double[real_total];
double* x_d, *y_d;
for(long long i = 0; i < real_total; i++) {
    x_h[i] = (double)rand() / RAND_MAX;
    y_h[i] = (double)rand() / RAND_MAX;
}
cudaMalloc(&count_d, real_total * sizeof(bool));  // 用于保存结果的显存
cudaMalloc(&x_d, real_total * sizeof(double));    // 随机数数组x
cudaMalloc(&y_d, real_total * sizeof(double));    // 随机数数组y
cudaMemcpy(x_d, x_h, real_total * sizeof(double), cudaMemcpyHostToDevice);  // 拷贝随机数数组
cudaMemcpy(y_d, y_h, real_total * sizeof(double), cudaMemcpyHostToDevice);  // 拷贝随机数数组

调用核函数：

trial<<<blocks, threads>>>(seed, count_d, x_d, y_d);

将结果拷贝回内存，并统计每个线程的结果，然后计算 $\pi$，最后打印结果。

cudaMemcpy(count_h, count_d, real_total * sizeof(bool), cudaMemcpyDeviceToHost);

long long count = 0;
for(long long i = 0; i < real_total; i++) {
    if(count_h[i]) {
        count++;
    }
}
double pi = 4 * (double)count / real_total;

printf("[+] total = %lld\n", real_total);  // 实际的total可能与参数不同，取决于是否整除
printf("[+] count = %lld\n", count);
printf("[+] pi    = %f\n", pi);
printf("[+] loss  = %e\n", acos(-1) - pi);

printf("\nBlocks  = %d\n", blocks.x);
printf("Threads = %d\n", threads.x);

由于整除问题，实际进行的测试次数real_total与提供的参数total不一定相等。

结果

注，本程序在 Windows 上运行，因 WSL 无法使用 GPU。

编译、运行并计时，采用每个 block 中 128 个线程进行总共 $1\times 10^8$ 次试验。

nvcc 4_mc_pi_cuda.cu -Xcompiler "/wd 4819" -o 4_mc_pi_cuda && nvprof ./4_mc_pi_cuda 100000000

结果：

[+] total = 100000000
[+] count = 78533817
[+] pi    = 3.141353
[+] loss  = 2.399736e-04

Blocks  = 781250
Threads = 128

本文地址：https://www.jeddd.com/article/parallel-monte-carlo-method-estimate-pi.html

结果分析

无论使用那种程序，其结果应该是相似的，也就是说误差是相似的，与具体实现无关。绝对误差的绝对值与样本数量的关系如下（多次实验取平均值）：

下面比较不同程序的性能。对于每一种参数（线程或进程数量），都测试样本数量分别为 100 万、1000 万、1 亿的三种测试；为了排除偶然因素，每种测试分别进行三次，取平均值。可以使用脚本来进行批量测试。

总结

• 在用蒙特卡洛方法求圆周率的实验中，样本数量（即试验次数）越多，得到结果的误差就越小，精度越高。
• OpenMP 是多线程并行，MPI 是多进程并行，CUDA 则是把任务分为多个线程放到 GPU 上执行。
• 比较 OpenMP 和 MPI 程序性能，发现在线程或进程数量相等的情况下，对于相同的样本数量，OpenMP 程序都比 MPI 略快。这可能是由于线程创建和切换的开销比进程创建和切换的开销更小导致的。
• 表中 CUDA 的时间为核函数的执行时间，不包括随机数生成、内存拷贝。与 OpenMP 和 MPI 直接比较可能是不合适的。单看 CUDA 程序的执行时间，发现核函数的执行时间与每个 block 中的线程数量的关系不大（线程总数是相同的）。
• 关于随机数的生成，可以看到在串行程序和 MPI 程序中，我使用rand()来产生随机数，这是最普通的一种方式；在 OpenMP 程序中，使用的是rand_r()，；在 CUDA 中，是现在 host 中用rand()生成随机数并保存到数组中，再将随机数数组拷贝至 device。不同的随机数产生方式对程序性能会有一定影响。
• 本实验主要用来巩固 OpenMP、MPI、CUDA 编程方法。

附录：完整代码

1_mc_pi_serial.c

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

int main(int argc, char* argv[]) {
    srand(time(NULL));      // 随机初始化seed
    long long total = 1e6;  // 默认100万个样本

    if(argc >= 2) {
        total = atoi(argv[1]);  // 从参数中获取样本数
    }

    long long count = 0;
    double x, y;
    for(long long i = 0; i < total; i++) {
        x = (double)rand() / RAND_MAX;
        y = (double)rand() / RAND_MAX;
        if(x*x + y*y <= 1) {
            count++;
        }
    }
    double pi = 4 * (double)count / total;

    printf("[+] total = %lld\n", total);
    printf("[+] count = %lld\n", count);
    printf("[+] pi    = %f\n", pi);
    printf("[+] loss  = %e\n", acos(-1) - pi);

    return 0;
}

2_mc_pi_openmp.c

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

int main(int argc, char* argv[]) {
    long long total = 1e6;  // 默认100万个样本
    int tn = 2;             // 默认2个线程
    if(argc >= 2) {
        total = atoi(argv[1]);  // 从参数中获取样本数
    }
    if(argc >= 3) {
        tn = atoi(argv[2]);     // 从参数中获得线程数
    }

    long long count = 0;
    double x, y;
    #pragma omp parallel num_threads(tn)
    {
        unsigned seed = time(NULL);
        
        #pragma omp for private(x, y) reduction(+:count)
        for(long long i = 0; i < total; i++) {
            x = (double)rand_r(&seed) / RAND_MAX;
            y = (double)rand_r(&seed) / RAND_MAX;
            if(x*x + y*y <= 1) {
                count++;
            }
        }
    }
    double pi = 4 * (double)count / total;

    printf("[+] total = %lld\n", total);
    printf("[+] count = %lld\n", count);
    printf("[+] pi    = %f\n", pi);
    printf("[+] loss  = %e\n", acos(-1) - pi);
    printf("\nnum_threads = %d\n", tn);

    return 0;
}

3_mc_pi_mpi.c

#include <stdio.h>
#include <stdlib.h> 
#include <time.h>
#include <math.h>
#include <mpi.h>

int main(int argc, char** argv) {
    int world_size, my_rank;
    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &world_size);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);

    long long total;
    long long local_total;
    long long count;
    long long local_count = 0;
    double x, y;
    
    srand(time(NULL));      // 随机初始化seed

    if(my_rank == 0) {
        total = 1e6;  // 默认100万个样本
        if(argc >= 2) {
            total = atoi(argv[1]);
        }
        // 在root中计算local_total
        local_total = total / world_size; 
    }
    // 把local_total广播给所有进程
    MPI_Bcast(&local_total, 1, MPI_LONG_LONG, 0, MPI_COMM_WORLD);

    for(long long i = 0; i < local_total; i++) {
        x = (double)rand() / RAND_MAX;
        y = (double)rand() / RAND_MAX;
        if(x*x + y*y <= 1) {
            local_count++;
        }
    }
    MPI_Reduce(&local_count, &count, 1, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
    double pi = 4 * (double)count / total;

    if(my_rank == 0) {
        printf("[+] total = %lld\n", world_size * local_total);
        printf("[+] count = %lld\n", count);
        printf("[+] pi    = %f\n", pi);
        printf("[+] loss  = %e\n", acos(-1) - pi);
        printf("\nProcesses = %d\n", world_size);
    }

    MPI_Finalize();
    return 0;
}

4_mc_pi_cuda.cu

#include <stdio.h>
#include <stdlib.h> 
#include <time.h>
#include <math.h>
// #include <curand_kernel.h>
#include "error_checks.h"

__global__ void trial(int seed, bool count_d[], double x_d[], double y_d[]) {
    // curandState state;
    long long id = blockIdx.x * blockDim.x + threadIdx.x;
    // curand_init(seed, id, 0, &state);
    // double x = (double)curand_uniform_double(&state);
    // double y = (double)curand_uniform_double(&state);
    double x = x_d[id], y = y_d[id];
    if(x*x + y*y <= 1) {
        count_d[id] = true;
    }
    else {
        count_d[id] = false;
    }
}

int main(int argc, char* argv[]) {
    int seed = time(NULL);
    long long total = 1e6;  // 默认100万个样本
    int tn = 128;             // 默认128个线程

    if(argc >= 2) {
        total = atoi(argv[1]);  // 从参数中获取样本数
    }
    if(argc >= 3) {
        tn = atoi(argv[2]);     // 从参数中获得线程数
    }
    dim3 threads(tn);
    dim3 blocks((total+tn-1) / tn);
    long long real_total = threads.x * blocks.x;

    bool* count_h = new bool[real_total];
    bool* count_d;
    double* x_h = new double[real_total];
    double* y_h = new double[real_total];
    double* x_d, *y_d;
    for(long long i = 0; i < real_total; i++) {
        x_h[i] = (double)rand() / RAND_MAX;
        y_h[i] = (double)rand() / RAND_MAX;
    }
    cudaMalloc(&count_d, real_total * sizeof(bool));  // 用于保存结果的显存
    cudaMalloc(&x_d, real_total * sizeof(double));    // 随机数数组x
    cudaMalloc(&y_d, real_total * sizeof(double));    // 随机数数组y
    cudaMemcpy(x_d, x_h, real_total * sizeof(double), cudaMemcpyHostToDevice);  // 拷贝随机数数组
    cudaMemcpy(y_d, y_h, real_total * sizeof(double), cudaMemcpyHostToDevice);  // 拷贝随机数数组

    trial<<<blocks, threads>>>(seed, count_d, x_d, y_d);
    
    cudaMemcpy(count_h, count_d, real_total * sizeof(bool), cudaMemcpyDeviceToHost);

    long long count = 0;
    for(long long i = 0; i < real_total; i++) {
        if(count_h[i]) {
            count++;
        }
    }
    double pi = 4 * (double)count / real_total;

    printf("[+] total = %lld\n", real_total);  // 实际的total可能与参数不同，取决于是否整除
    printf("[+] count = %lld\n", count);
    printf("[+] pi    = %f\n", pi);
    printf("[+] loss  = %e\n", acos(-1) - pi);

    printf("\nBlocks  = %d\n", blocks.x);
    printf("Threads = %d\n", threads.x);

    return 0;
}

本文地址：https://www.jeddd.com/article/parallel-monte-carlo-method-estimate-pi.html