蒙特卡洛方法计算圆周率的三种实现-MPI openmp pthread

蒙特卡洛方法实现计算圆周率的方法比较简单，其思想是假设我们向一个正方形的标靶上随机投掷飞镖，靶心在正中央，标靶的长和宽都是2 英尺。同时假设有一个圆与标靶内切。圆的半径是1英尺，面积是π平方英尺。如果击中点在标靶上是均匀分布的（我们总会击中正方形），那么飞镖击中圆的数量近似满足等式

飞镖落在圆内的次数/飞镖落在标靶内的总次数=π/4

因为环包含的面积与正方形面积的比值是π/4。

因为环所包含的面积与正方形面积的比值是π/4。

我们可以用这个公式和随机数产生器来估计π的值。

伪代码如下：

number_in_circle=;

for(toss=;toss<number_of_tosses;toss++){

     x=random double between - and ;

     y=random double between - and ;

     distance_squared=x*x+y*y;

    if(distance_squared<=) number_in_cycle++;

}

pi_estimate=*number_in_cycle/((double)number_in_tosses);

这种采用了随机（随机投掷）的方法称为蒙特卡洛（Monte Carlo）方法。

编写了一个采用蒙特卡洛方法的MPI，Pthread，OpenMP程序估计π的值。

使用MPI编写时，进程0读取总的投掷次数，并把它们广播给各个进程。使用MPI_Reduce求出局部变量number_in_cycle的全局总和，并让进程0打印它。

使用Pthread编写时，有主线程读入总的投掷数，然后输出估算值。

使用OpenMP编写时，在开启任何线程前读取总的投掷次数。使用reduction子句计算飞镖集中环内的次数。在合并所有的线程后，打印结果。

这三个程序中，投掷次数可能都非常大，可能总的投掷次数和击中环内的次数都得用long int型来表示。

下面是MPI程序代码

 #include<stdio.h>

 #include<stdlib.h>

 #include<math.h>

 #include<time.h>

 #include<mpi.h>

 void read_num(long long int *num_point,int my_rank,MPI_Comm comm);

 void compute_pi(long long int num_point,long long int* num_in_cycle,long long int* local_num_point,int comm_sz,long long int *total_num_in_cycle,MPI_Comm comm,int my_rank);

 int main(int argc,char** argv){

     long long int num_in_cycle,num_point,total_num_in_cycle,local_num_point;

     int my_rank,comm_sz;

     MPI_Comm comm;

     MPI_Init(NULL,NULL);//初始化

     comm=MPI_COMM_WORLD;

     MPI_Comm_size(comm,&comm_sz);//得到进程总数

     MPI_Comm_rank(comm,&my_rank);//得到进程编号

     read_num(&num_point,my_rank,comm);//读取输入数据

     compute_pi(num_point,&num_in_cycle,&local_num_point,comm_sz,&total_num_in_cycle,comm,my_rank);

     MPI_Finalize();

     return ;

 }

 void read_num(long long int* num_point,int my_rank,MPI_Comm comm){

     if(my_rank==){

         printf("please input num in sqaure \n");

         scanf("%lld",num_point);

     }

     /*

     广播函数

     int MPI_Bcast(

         void*    data_p //in/out

         int    count  //in

         MPI_Datatype    datatype //in

         int    source_proc  //in

         MPI_Comm    comm  //in

     )

     */

     MPI_Bcast(num_point,,MPI_LONG_LONG,,comm);

 }

 void compute_pi(long long int num_point,long long int* num_in_cycle,long long int* local_num_point,int comm_sz,long long int *total_num_in_cycle,MPI_Comm comm,int my_rank){

     *num_in_cycle=;

     *local_num_point=num_point/comm_sz;

     double x,y,distance_squared;

     srand(time(NULL));

     for(long long int i=;i< *local_num_point;i++){

         x=(double)rand()/(double)RAND_MAX;

         x=x*-;

         y=(double)rand()/(double)RAND_MAX;

         y=y*-;

         distance_squared=x*x+y*y;

         if(distance_squared<=)

         *num_in_cycle=*num_in_cycle+;

     }

     /*

     全局函数

     MPI_Reduce(

         void*    input_data_p    //in

         void*    output_data_p    //out

         int    count        //in

         MPI_Datatype    datatype     //in

         MPI_Op    oprtator    //in

         int    dest_process    //in

         MPI_Comm    comm    //in

     )

     */

       MPI_Reduce(num_in_cycle,total_num_in_cycle,,MPI_LONG_LONG,MPI_SUM,,comm);

     if(my_rank==){

         double pi=(double)*total_num_in_cycle/(double)num_point*;

         printf("the estimate value of pi is %lf\n",pi);

     }

 }

进行编译 mpicc -std=c99 -o mpi_mete mpi_mete.c

执行 mpiexec -n mpi_mete

输入数据

下面是Pthread程序代码

 #include<stdlib.h>

 #include<stdio.h>

 #include<math.h>

 #include<pthread.h>

 int thread_count;

 long long int num_in_circle,n;

 pthread_mutex_t mutex;

 void* compute_pi(void* rank);

 int main(int argc,char* argv[]){

     long    thread;

     pthread_t* thread_handles;

     thread_count=strtol(argv[],NULL,);

     printf("please input the number of point\n");

     scanf("%lld",&n);

     thread_handles=(pthread_t*)malloc(thread_count*sizeof(pthread_t));

     pthread_mutex_init(&mutex,NULL);

     for(thread=;thread<thread_count;thread++){

         //创建线程

         /*

         int pthread_create(

             pthread_t*    thread_p  //out

             const    pthread_attr_t*    attr_p

             void*    (*start_routine)(void*)    //in

             void*    arg_p    //in

         )

         第一个参数是一个指针，指向对应的pthread_t对象。注意，pthread_t对象不是pthread_create函数分配的，必须在调用pthread_create函数前就为pthread_t

         对象分配内存空间。第二个参数不用，所以只是函数调用时把NULL传递参数。第三个参数表示该线程将要运行的函数；最后一个参数也是一个指针，指向传给函数start_routine的参数

         */

         pthread_create(&thread_handles[thread],NULL,compute_pi,(void*)thread);

     }

     //停止线程

     /*

     int pthread_join(

         pthread_t    thread    /in

         void**    ret_val_p    /out

     ）

     第二个参数可以接收任意由pthread_t对象所关联的那个线程产生的返回值。

     */

     for(thread=;thread<thread_count;thread++){

         pthread_join(thread_handles[thread],NULL);

     }

     pthread_mutex_destroy(&mutex);

     double pi=*(double)num_in_circle/(double)n;

     printf("the esitimate value of pi is %lf\n",pi);

 }

 void* compute_pi(void* rank){

     long long int local_n;

     local_n=n/thread_count;

     double x,y,distance_squared;

     for(long long int i=;i<local_n;i++){

         x=(double)rand()/(double)RAND_MAX;

         y=(double)rand()/(double)RAND_MAX;

         distance_squared=x*x+y*y;

         if(distance_squared<=){

             pthread_mutex_lock(&mutex);

             num_in_circle++;

             pthread_mutex_unlock(&mutex);

         }

     }

     return NULL;

 }

进行编译 gcc -std=c99 -o pth_mete pth_mete.c

运行 ./pth_mete

输入数据结果如下

下面是OpenMP代码

 #include<stdlib.h>

 #include<stdio.h>

 #include<time.h>

 #include<omp.h>

 int main(int argc,char** argv){

     long long int num_in_cycle=;

     long long int num_point;

     int thread_count;

     thread_count=strtol(argv[],NULL,);

     printf("please input the number of point\n");

     scanf("%lld",&num_point);

     srand(time(NULL));

     double x,y,distance_point;

     long long int i;

     #pragma omp parallel for num_threads(thread_count) default(none) \

         reduction(+:num_in_cycle) shared(num_point) private(i,x,y,distance_point)

     for( i=;i<num_point;i++){

         x=(double)rand()/(double)RAND_MAX;

         y=(double)rand()/(double)RAND_MAX;

         distance_point=x*x+y*y;

         if(distance_point<=){

             num_in_cycle++;

         }

     }

     double estimate_pi=(double)num_in_cycle/num_point*;

     printf("the estimate value of pi is %lf\n",estimate_pi);

     return ;

 }

编译代码 gcc -std=c99 -fopenmp -o omp_mete omp_mete.c

执行 ./omp_mete

结果如下

巴特西

蒙特卡洛方法计算圆周率的三种实现-MPI openmp pthread

最新文章

热门文章