Compute PI in parallel

最新推荐文章于 2026-06-19 07:39:45 发布

转载最新推荐文章于 2026-06-19 07:39:45 发布 · 1k 阅读

标签

#parallel #parsing #fortran #分布式存储 #float

HPC 专栏收录该内容

1 篇文章

订阅专栏

本文介绍了两种计算圆周率的并行计算方法，一种基于数值积分，另一种基于概率统计。通过OpenMP和MPI实现这两种方法，并展示了如何利用并行处理来提高计算效率。

并行计算圆周率

看到这个题目，俗了，大家都在计算圆周率。不过咱们的目的是看一下并行计算的基本流程。

书上计算PI用的是精确的数值计算方法，我这里再给出一种概率计算方法。

OpenMP和MPI将同时亮相。

计算PI的方法

1.tan(PI/4)=1 => PI=4arctan1。知道arctan1转化为定积分的形式是什么吧。

 
     01#include<stdio.h>
 
     02#include<time.h>
 
     03#define N 1000000
 
     04main(){
 
     05    double local,pi=0.0,w;
 
     06    long i;
 
     07    w=1.0/N;
 
     08    clock_t t1=clock();
 
     09    for(i=0;i<N;i++){
 
     10        local=(i+0.5)*w;
 
     11        pi=pi+4.0/(1.0+local*local);
 
     12    }
 
     13    clock_t t2=clock();
 
     14    printf("PI is %.20f/n",pi*w);
 
     15    printf("Time: %.2f seconds/n",(float)(t2-t1)/CLOCKS_PER_SEC);
 
     16}

orisun@orisun-desktop:~/Program$ ./PI1

PI is 3.14159265358976336202

Time: 0.02 seconds

2.以坐标原点为形心，作半径为1的圆和边长为2的正方形。正方形与圆的面积之比即为PI

 
     01#include<stdio.h>
 
     02#include<stdlib.h>
 
     03#include<time.h>
 
     04#include<math.h>
 
     05#define N 1000000
 
     06main(){
 
     07    long i,sum;
 
     08    double x,y;
 
     09    srand((unsigned)time(NULL));
 
     10    sum=0;
 
     11    clock_t t1=clock();
 
     12    for(i=0;i<N;i++){
 
     13        x=(double)rand()/RAND_MAX;
 
     14        y=(double)rand()/RAND_MAX;
 
     15        if(x*x+y*y<1)
 
     16            sum++;
 
     17    }
 
     18    clock_t t2=clock();
 
     19    printf("PI is %.20f/n",4*(double)sum/N);
 
     20    printf("Time: %.2f/n",(float)(t2-t1)/CLOCKS_PER_SEC);
 
     21}

orisun@orisun-desktop:~$ ./PI0

PI is 3.14301599999999980994

Time: 0.16

对比可以看到方法1在计算精度和速度上都具有绝对的优势。在下面的openMP和MPI计算中我们都采用方法1。

OpenMP

OpenMP[OMP]是一个编译器指令和库函数的集合（已包含在gcc中），它用于为共享存储器计算机创建并行程序。OMP组合了C、C++和Fortran。

 
     01#include<stdio.h>
 
     02#include<time.h>
 
     03#include<omp.h>
 
     04#define N 1000000
 
     05main(){
 
     06    double local,pi=0.0,w;
 
     07    long i;
 
     08    w=1.0/N;
 
     09    clock_t t1=clock();
 
     10#pragma omp parallel for private(local) reduction(+:pi)
 
     11    for(i=0;i<N;i++){
 
     12        local=(i+0.5)*w;
 
     13        pi=pi+4.0/(1.0+local*local);
 
     14    }
 
     15    clock_t t2=clock();
 
     16    printf("PI is %.20f/n",pi*w);
 
     17    printf("Time: %.2f seconds/n",(float)(t2-t1)/CLOCKS_PER_SEC);
 
     18}

orisun@orisun-desktop:~/Program$ ./PI2

PI is 3.14159265358976336202

Time: 0.02 seconds

跟串行计算结果是一模一样。

#pragma omp parallel表示下面的一行代码或代码块要分配到多个执行单元中并行计算。

#pragma omp parallel for用在一个for循环的前面

private(local)默认情况下定义在并行代码之外的变量为各并行的执行单元所共享，使用private限制，表示每个执行单元创建该变量的一个副本

reduction(+:pi)表示并行代码执行完毕后对各个执行单元中的pi进行相加操作

MPICH2

ubuntu下mpich2的安装配置见这位老兄的博客http://hi.baidu.com/mousetwo/blog/item/ddba71a9b0adf4f11f17a2fc.html

MPI（Message Parsing Interface）消息传递接口是用于分布式存储器并行计算机的标准编程环境。MPI的核心构造是消息传递：一个进程将信息打包成消息，并将该消息发送给其他进程。MPI最常用的两个实现是LAM/MPI[LAM]和MPICH[MPI]。

在MPI中执行单元（UE）指的就是进程。

 
     01#include<stdio.h>
 
     02#include<mpi.h>
 
     03#include<math.h>
 
     04 
 
     05int main(int argc,char *argv[]){
 
     06    int my_rank,num_procs;
 
     07    int i,n=0;
 
     08    double sum,width,local,mypi,pi;
 
     09    double start=0.0,stop=0.0;
 
     10    int proc_len;
 
     11    char processor_name[MPI_MAX_PROCESSOR_NAME];
 
     12 
 
     13    MPI_Init(&argc,&argv);          //初始化环境
 
     14    MPI_Comm_size(MPI_COMM_WORLD,&num_procs);   //获取并行的进程数
 
     15    MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);     //当前进程在所有进程中的序号
 
     16    MPI_Get_processor_name(processor_name,&proc_len);   //获取总的处理机数和各个处理机的名称
 
     17 
 
     18    printf("Processor %d of %d on %s/n",my_rank,num_procs,processor_name);
 
     19    if(my_rank==0){
 
     20        printf("please give n=");
 
     21        scanf("%d",&n);
 
     22        start=MPI_Wtime();              //MPI计时
 
     23    }
 
     24    MPI_Bcast(&n,1,MPI_INT,0,MPI_COMM_WORLD);   //把n广播给本通信环境中的所有进程
 
     25    width=1.0/n;
 
     26    sum=0.0;
 
     27    for(i=my_rank;i<n;i+=num_procs){
 
     28        local=width*((double)i+0.5);
 
     29        sum+=4.0/(1.0+local*local);
 
     30    }
 
     31    mypi=width*sum;
 
     32    MPI_Reduce(&mypi,&pi,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD);    //由进程0进行归约，把每个进程计算出来的mypi进行相加（MPI_SUM）,赋给pi
 
     33    if(my_rank==0){
 
     34        printf("PI is %.20f/n",pi);
 
     35        stop=MPI_Wtime();
 
     36        printf("Time: %f/n",stop-start);
 
     37        fflush(stdout);
 
     38    }
 
     39    MPI_Finalize();
 
     40    return 0;
 
     41}