MPI和OpenMP实现蒙特卡罗算法

MPI和OpenMP实现蒙特卡罗算法

一、蒙特卡洛算法介绍

1. 基本思想

   当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。

2. 数学应用：

   通常蒙特·卡罗方法通过构造符合一定规则的随机数来解决数学上的各种问题。对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题，蒙特·卡罗方法是一种有效的求出数值解的方法。一般蒙特·卡罗方法在数学中最常见的应用就是蒙特·卡罗积分。

3. 案例：

   通过在正方形内随机撒点，落在圆内的点 / 落在正方形内的点，就约等于圆的面积 / 正方形的面积 = $\frac{\pi }{4}$，随着撒点的数量增多，$\pi$的数值越来越接近真实值
   $$\frac{N_{圆}}{N_{正方形}}\approx \frac{S_{圆}}{S_{正方形}}=\frac{\pi }{4}(N为点的数量，S为面积)$$

二、MPI实现蒙特卡罗算法计算$\pi$值

1. 首先要配置要MPI多节点集群环境，没有配置好的可以参考我的前几篇博客

2. 原理

   1. 生成0-1内的随机数$x,y$记为$(x,y)$，总共生成的点数为$N$
   2. 计算每一个$(x,y)$是否满足$x^2+y^2\le 1$，将符合条件的点数记为$C$
   3. $\pi$的近似值可以看作$\frac{4C}{N}$

3. 运行

   1. 编写profile文件，与上一节一样

   node01:2
   node02:2
   

   2. 编写源代码文件mpimc.cpp，两个虚拟机都需要

   3. 编译mpimc.cpp文件，两个虚拟机都需要编译

   mpic++ mpimc.cpp -o mpimc
   

   4. mpi多节点运行，-n参数推荐为2，只有main节点运行

   mpiexec -f profile -n 2 ./mpimc
   

4. 运行结果

5. 代码

#include <iostream>
#include "mpi.h"

using namespace std;
const int maxn = 1e6; //随机点数量

MPI_Status status;
int main(int argc, char *argv[]) {
    int pCnt, pId, slaves, source, dest;
    int cnt, n, offset;
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &pId);
    MPI_Comm_size(MPI_COMM_WORLD, &pCnt);

    slaves = pCnt - 1;
    offset = maxn / slaves;
    if(pId == 0) {
        for(dest = 1; dest <= slaves; dest++) {
            MPI_Send(&offset, 1, MPI_INT, dest, 1, MPI_COMM_WORLD);
            MPI_Send(&n, 1, MPI_INT, dest, 1, MPI_COMM_WORLD);
        }
        cnt = 0;
        for(int i = 1; i <= slaves; i++) {
            source = i;
            MPI_Recv(&offset, 1, MPI_INT, source, 2, MPI_COMM_WORLD, &status);
            MPI_Recv(&n, 1, MPI_INT, source, 2, MPI_COMM_WORLD, &status);
            cnt += n;
        }

        double pi = 4.0 * cnt / maxn;
        cout << "PI is " << pi << endl;
    }

    if(pId > 0) {
        source = 0;
        MPI_Recv(&offset, 1, MPI_INT, source, 1, MPI_COMM_WORLD, &status);
        MPI_Recv(&n, 1, MPI_INT, source, 1, MPI_COMM_WORLD, &status);
        n = 0;
        srand(time(NULL));
        for(int i = 0; i < offset; i++) {
            // 生成0-1内的随机数
            double x = (double)rand() / RAND_MAX;
            double y = (double)rand() / RAND_MAX;
            if(x * x + y * y <= 1)
                n++;
        }
        MPI_Send(&offset, 1, MPI_INT, 0, 2, MPI_COMM_WORLD);
        MPI_Send(&n, 1, MPI_INT, 0, 2, MPI_COMM_WORLD);
    }

    MPI_Finalize();
    return 0;
}

三、OpenMP实现蒙特卡罗算法计算$\pi$值

1. 不需要MPI集群环境，只需要一个节点就行
2. 原理（同上）
3. 编译运行

g++ ompmc.cpp -fopenmp -o ompmc
./ompmc

4. 运行结果

5. 代码

#include <iostream>
#include <pthread.h>
#include <omp.h>
using namespace std;
const int maxn = 1e6; //随机点数量


double dis(double x, double y) {
    return x * x + y * y;
}

int main(int argc, char *argv[]) {
    int cnt = 0;
    double x, y;
    int i;

    omp_set_num_threads(omp_get_num_procs());
    srand48(time(NULL));

    #pragma omp parallel for private(i)
    for(i = 0; i < maxn; i++) {
        // 生成0-1内的随机数
        x = (double)rand() / RAND_MAX;
        y = (double)rand() / RAND_MAX;
        if(dis(x, y) <= 1) cnt++;
    }
    double pi = 4.0 * cnt / maxn;
    cout << "Pi is " << pi << endl;
    return 0;
}