使用omp和mpi并行技术实现蒙特卡罗算法实验（MonteCarlo）

蒙特卡洛算法思想: 画一个正方形和一个内切圆然后在图形上面随机撒点，数一下落在圆

内和正方形内的点的数量之比，就是二者面积之比, 也就是π/4。

大数定理告诉我们，随着样本数量的增大，我们用这种方式模拟出来的值应该是越来越趋

近于真实值，样本无穷大的时候收敛于真值。这就证明了蒙特卡罗方法的合理性和有效性。

Mpi基本原理：

1.什么是MPI

Massage Passing Interface:是消息传递函数库的标准规范，由MPI论坛开发。

一种新的库描述，不是一种语言。共有上百个函数调用接口，提供与C和Fortran语言的绑定

MPI是一种标准或规范的代表，而不是特指某一个对它的具体实现

MPI是一种消息传递编程模型，并成为这种编程模型的代表和事实上的标准

2.MPI的特点

MPI有以下的特点：

消息传递式并行程序设计

指用户必须通过显式地发送和接收消息来实现处理机间的数据交换。

在这种并行编程中，每个并行进程均有自己独立的地址空间，相互之间访问不能直接进行，必须通过显式的消息传递来实现。

这种编程方式是大规模并行处理机（MPP）和机群（Cluster）采用的主要编程方式。

并行计算粒度大，特别适合于大规模可扩展并行算法

用户决定问题分解策略、进程间的数据交换策略，在挖掘潜在并行性方面更主动,并行计算粒度大,特别适合于大规模可扩展并行算法

消息传递是当前并行计算领域的一个非常重要的并行程序设计方式

二、MPI的基本函数

MPI调用借口的总数虽然庞大，但根据实际编写MPI的经验，常用的MPI函数是以下6个：

MPI_Init(…);

MPI_Comm_size(…);

MPI_Comm_rank(…);

MPI_Send(…);

MPI_Recv(…);

MPI_Finalize();

三、MPI的通信机制

MPI是一种基于消息传递的编程模型，不同进程间通过消息交换数据。

1.MPI点对点通信类型

所谓点对点的通信就是一个进程跟另一个进程的通信，而下面的聚合通信就是一个进程和多个进程的通信。

标准模式：

该模式下MPI有可能先缓冲该消息，也可能直接发送，可理解为直接送信或通过邮局送信。是最常用的发送方式。

由MPI决定是否缓冲消息

没有足够的系统缓冲区时或出于性能的考虑，MPI可能进行直接拷贝：仅当相应的接收完成后，发送语句才能返回。

这里的系统缓冲区是指由MPI系统管理的缓冲区。而非进程管理的缓冲区。

MPI环境定义有三种缓冲区：应用缓冲区、系统缓冲区、用户向系统注册的通信用缓冲区

MPI缓冲消息：发送语句在相应的接收语句完成前返回。

这时后发送的结束或称发送的完成== 消息已从发送方发出，而不是滞留在发送方的系统缓冲区中。

该模式发送操作的成功与否依赖于接收操作，我们称之为非本地的，即发送操作的成功与否跟本地没关系

mpi实现MonteCarlo

#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;
}

Mpi：

Mpicxx ./xxx.cpp -o ./xxx 编译；

Mpirun -n p ./xxx运行；

运行结果：

Mpi：

Omp基本原理：

OpenMP是一种用于共享内存并行系统的多线程程序设计方案，支持的编程语言包括C、C++和Fortran。OpenMP提供了对并行算法的高层抽象描述，特别适合在多核CPU机器上的并行程序设计。编译器根据程序中添加的pragma指令，自动将程序并行处理，使用OpenMP降低了并行编程的难度和复杂度。当编译器不支持OpenMP时，程序会退化成普通（串行）程序。程序中已有的OpenMP指令不会影响程序的正常编译运行。

OpenMP执行模式

OpenMP采用fork-join的执行模式。开始的时候只存在一个主线程，当需要进行并行计算的时候，派生出若干个分支线程来执行并行任务。当并行代码执行完成之后，分支线程会合，并把控制流程交给单独的主线程。

一个典型的fork-join执行模型的示意图如下：

OpenMP编程模型以线程为基础，通过编译制导指令制导并行化，有三种编程要素可以实现并行化控制，他们分别是编译制导、API函数集和环境变量。

编译制导

编译制导指令以#pragma omp 开始，后边跟具体的功能指令，格式如：#pragma omp 指令[子句[,子句] …]。常用的功能指令如下：

parallel：用在一个结构块之前，表示这段代码将被多个线程并行执行；

for：用于for循环语句之前，表示将循环计算任务分配到多个线程中并行执行，以实现任务分担，必须由编程人员自己保证每次循环之间无数据相关性；

parallel for：parallel和for指令的结合，也是用在for循环语句之前，表示for循环体的代码将被多个线程并行执行，它同时具有并行域的产生和任务分担两个功能；

sections：用在可被并行执行的代码段之前，用于实现多个结构块语句的任务分担，可并行执行的代码段各自用section指令标出（注意区分sections和section）；

parallel sections：parallel和sections两个语句的结合，类似于parallel for；

single：用在并行域内，表示一段只被单个线程执行的代码；

critical：用在一段代码临界区之前，保证每次只有一个OpenMP线程进入；

flush：保证各个OpenMP线程的数据影像的一致性；

barrier：用于并行域内代码的线程同步，线程执行到barrier时要停下等待，直到所有线程都执行到barrier时才继续往下执行；

atomic：用于指定一个数据操作需要原子性地完成；

master：用于指定一段代码由主线程执行；

threadprivate：用于指定一个或多个变量是线程专用，后面会解释线程专有和私有的区别。

相应的OpenMP子句为：

private：指定一个或多个变量在每个线程中都有它自己的私有副本；

firstprivate：指定一个或多个变量在每个线程都有它自己的私有副本，并且私有变量要在进入并行域或任务分担域时，继承主线程中的同名变量的值作为初值；

lastprivate：是用来指定将线程中的一个或多个私有变量的值在并行处理结束后复制到主线程中的同名变量中，负责拷贝的线程是for或sections任务分担中的最后一个线程；

reduction：用来指定一个或多个变量是私有的，并且在并行处理结束后这些变量要执行指定的归约运算，并将结果返回给主线程同名变量；

nowait：指出并发线程可以忽略其他制导指令暗含的路障同步；

num_threads：指定并行域内的线程的数目；

schedule：指定for任务分担中的任务分配调度类型；

shared：指定一个或多个变量为多个线程间的共享变量；

ordered：用来指定for任务分担域内指定代码段需要按照串行循环次序执行；

copyprivate：配合single指令，将指定线程的专有变量广播到并行域内其他线程的同名变量中；

copyin：用来指定一个threadprivate类型的变量需要用主线程同名变量进行初始化；

default：用来指定并行域内的变量的使用方式，缺省是shared。

利用omp_set_num_threads()来设置线程数，

利用#pragma omp parallel sections 声明下面大括号中的语句要并行多线程执行；

利用#pragma omp section 分配线程。

omp实现MonteCarlo

#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;
}

Omp：