【MPI OpenMP】并行蒙特卡洛方法求圆周率(C语言)

本文记录了使用MPI与OpenMP两种并行计算方法实现蒙特卡洛计算圆周率，题目是专业实验课上老师给的，主要分享一下自己的做法，希望大家不吝赐教（使用的语言是C语言）。

蒙特卡洛方法求圆周率

• 蒙特卡洛方法的原理
  该方法的原理还是比较简单的，网上资料也有很多，我的实验课老师也提供了一个参考资料，如果大家想看看蒙特卡洛计算圆周率的原理，可以自行查看一下：蒙特卡洛方法求圆周率与定积分及python实现。这位博主分享的写法是python的，并且除了计算了圆周率，还计算了定积分，但我只是查看了其中的原理部分。
• 为什么要使用并行算法
  因为蒙特卡洛求圆周率就是模拟向正方形区域撒点嘛，计算落在圆中的点数，撒的点越多，计算出来的圆周率越精确，但是此时会存在一个问题，那就是点太多了之后，算法的运行时间会很长；而恰好模拟撒点这件事是可以并行实现的，也就是说，撒的点之间不影响，所以可以让多个进程或者线程同时进行撒点，最后将所有进/线程撒的点一起统计，得到精度更高的圆周率。
• 我使用的方法
  计算圆周率的撒点有很多方式，我是用的是使用单位圆在第一象限的部分进行模拟（如下图），第一象限中正方形的面积为1，四分之一圆的面积为pi/4，因此可以根据四分之一圆与正方形中的点数比来估计圆周率，如正方形中点数为N，四分之一圆内点数为n，则圆周率 pi = 4*N/n。
  

MPI方法

MPI方法是通过多线程方式实现的并行计算，方式很多，本人使用的是两个虚拟机，使用前需要配置两个虚拟机，并需要写相关配置文件，但是本文主要是分享蒙特卡洛方法，就不讲解如何配置虚拟机了。

• 思路
  模拟撒点的总点数是通过参数获取；进而可以判断总点数是否能整除总进程数， 将能够整除的点数部分均分给每个进程，剩余的不能均分的部分由根进程进行模拟；最后将每个进程模拟的结果统计，输出所估计的圆周率，并输出并行算法的时间。
• 代码

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

long long int get_num(long long int n)
{
    double x,y,d;   //模拟生成的点(x,y)，以及该点到原点的距离d
    long long int num = 0, i;   //落入四分之一圆的点数
    srand(time(NULL));   //随机种子
    for(i=0; i<n; i++){
        x=(double)rand()/(double)RAND_MAX;
        y=(double)rand()/(double)RAND_MAX;
        d=x*x+y*y;
        if(d<=1)
            num += 1;
    }
    return num;
}

int main(int argc,char** argv){
    long long int n;
    int size,rank;
    MPI_Comm comm = MPI_COMM_WORLD;
    MPI_Status status;
    MPI_Init(&argc,&argv); //初始化
    MPI_Comm_size(comm,&size);  //进程总数
    MPI_Comm_rank(comm,&rank);  //当前进程标识
    sscanf(argv[1],"%lld",&n);  //通过参数获取一共模拟多少个点

    clock_t startTime, endTime;   //并行时间
    startTime=clock();
    MPI_Bcast(&n,1,MPI_LONG_LONG,0,comm);   //将总数广播到每个进程
    
    long long int Num, num;   //所有进程的落入四分之一圆区域的总点数与每个进程的落入四分之一圆区域的总点数
    num = get_num(n/size);  //计算每个进程能均分的点数
    MPI_Barrier(comm);
    MPI_Reduce(&num,&Num,1,MPI_LONG_LONG,MPI_SUM,0,comm);   //将每个进程计算的结果求和
    if(rank==0){
        if(num%size != 0)   //根进程计算不能整除的点数
            Num += get_num(n%size);
        double pi=(double)Num/(double)n*4;   //根据开头分析的公式计算圆周率
        endTime=clock();
        printf("pi=%lf\n",pi);
        printf("并行计算时间time=%dms\n",int(endTime-startTime));  //输出并行时间
    }
    
    MPI_Finalize();
    return 0;
}

• 结果
  分别使用1、2、4个进程模拟8000000与10000000个点，结果如下
  

OpenMP

• 思路
  OpenMP的思路就十分简单粗暴了，使用多线程的思想，由于多线程是共享内存，所以无脑地将串行代码中的for循环修改成OpenMP并行代码即可，但是此处需要注意一个问题，由于多线程共享内存，直接将num设置为共享变量，容易导致多个线程共同使num加1的时候丢失数据（比如num为1，某个线程让num变成2，但是还没写进内存，另一个线程也将num读取，给num加1，并写入内存，所以此时会将前一个线程的操作丢失），不过，openmp提供了一个解决方案，在进行并行时，使用reduction声明num，即将每个线程都会对num进行拷贝，在所有线程计算结束后，将所有线程的num加起来
• 代码

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

int main(int argc,char** argv){
    long long int n, num = 0, i;
    double x,y,d;
    int num_thread = atoi(argv[2]);
    
    sscanf(argv[1],"%lld",&n);
    
    clock_t startTime, endTime;
    startTime=clock();
    
    srand(time(NULL));
    omp_set_num_threads(num_thread);
    #pragma omp parallel shared(n) private(i,x,y,d) reduction(+:num)
    {
        #pragma omp for schedule(guided)
        for(i=0; i<n; i++) {     
            x=(double)rand()/(double)RAND_MAX;
            y=(double)rand()/(double)RAND_MAX;
            d=x*x+y*y;
            if(d<=1)
                num += 1;
        }
    }
    double pi=(double)num/(double)n*4;
    endTime=clock();
    printf("pi=%lf\n",pi);
    printf("并行计算时间time=%dms\n",int(endTime-startTime));
    
    return 0;
}

• 结果