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

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

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

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

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

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

伪代码如下:

number_in_circle=0;
for(toss=0;toss<number_of_tosses;toss++){
     x=random double between -1 and 1;
     y=random double between -1 and 1;
     distance_squared=x*x+y*y;
    if(distance_squared<=1) number_in_cycle++;    
}
pi_estimate=4*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程序代码

 1 #include<stdio.h>
 2 #include<stdlib.h>
 3 #include<math.h>
 4 #include<time.h>
 5 #include<mpi.h>
 6 void read_num(long long int *num_point,int my_rank,MPI_Comm comm);
 7 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);
 8 int main(int argc,char** argv){
 9     long long int num_in_cycle,num_point,total_num_in_cycle,local_num_point;
10     int my_rank,comm_sz;
11     MPI_Comm comm;
12     MPI_Init(NULL,NULL);//初始化
13     comm=MPI_COMM_WORLD;
14     MPI_Comm_size(comm,&comm_sz);//得到进程总数
15     MPI_Comm_rank(comm,&my_rank);//得到进程编号
16     read_num(&num_point,my_rank,comm);//读取输入数据
17     compute_pi(num_point,&num_in_cycle,&local_num_point,comm_sz,&total_num_in_cycle,comm,my_rank);
18     MPI_Finalize();
19     return 0;
20 }
21 void read_num(long long int* num_point,int my_rank,MPI_Comm comm){
22     if(my_rank==0){
23         printf("please input num in sqaure \n");
24         scanf("%lld",num_point);
25     }
26     /*
27     广播函数
28     int MPI_Bcast(
29         void*    data_p //in/out
30         int    count  //in
31         MPI_Datatype    datatype //in
32         int    source_proc  //in
33         MPI_Comm    comm  //in    
34     )
35     */
36     MPI_Bcast(num_point,1,MPI_LONG_LONG,0,comm);
37 
38 }
39 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){
40     *num_in_cycle=0;
41     *local_num_point=num_point/comm_sz;
42     double x,y,distance_squared; 
43     srand(time(NULL));
44     for(long long int i=0;i< *local_num_point;i++){     
45         x=(double)rand()/(double)RAND_MAX;
46         x=x*2-1;
47         y=(double)rand()/(double)RAND_MAX;
48         y=y*2-1;
49         distance_squared=x*x+y*y;
50         if(distance_squared<=1)
51         *num_in_cycle=*num_in_cycle+1;
52     }
53     /*
54     全局函数
55     MPI_Reduce(
56         void*    input_data_p    //in
57         void*    output_data_p    //out
58         int    count        //in
59         MPI_Datatype    datatype     //in
60         MPI_Op    oprtator    //in
61         int    dest_process    //in
62         MPI_Comm    comm    //in
63     )    
64     */
65       MPI_Reduce(num_in_cycle,total_num_in_cycle,1,MPI_LONG_LONG,MPI_SUM,0,comm);
66     if(my_rank==0){
67         double pi=(double)*total_num_in_cycle/(double)num_point*4;
68         printf("the estimate value of pi is %lf\n",pi);
69     }
70 }

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

执行 mpiexec -n 4 mpi_mete 

输入数据

下面是Pthread程序代码

 

 1 #include<stdlib.h>
 2 #include<stdio.h>
 3 #include<math.h>
 4 #include<pthread.h>
 5 int thread_count;
 6 long long int num_in_circle,n;
 7 pthread_mutex_t mutex;
 8 void* compute_pi(void* rank);
 9 int main(int argc,char* argv[]){
10     long    thread;
11     pthread_t* thread_handles;
12     thread_count=strtol(argv[1],NULL,10);
13     printf("please input the number of point\n");
14     scanf("%lld",&n);
15     thread_handles=(pthread_t*)malloc(thread_count*sizeof(pthread_t));
16     pthread_mutex_init(&mutex,NULL);
17     for(thread=0;thread<thread_count;thread++){
18         //创建线程
19         /*
20         int pthread_create(
21             pthread_t*    thread_p  //out
22             const    pthread_attr_t*    attr_p
23             void*    (*start_routine)(void*)    //in
24             void*    arg_p    //in        
25         )
26         第一个参数是一个指针,指向对应的pthread_t对象。注意,pthread_t对象不是pthread_create函数分配的,必须在调用pthread_create函数前就为pthread_t
27         对象分配内存空间。第二个参数不用,所以只是函数调用时把NULL传递参数。第三个参数表示该线程将要运行的函数;最后一个参数也是一个指针,指向传给函数start_routine的参数
28         */
29         pthread_create(&thread_handles[thread],NULL,compute_pi,(void*)thread);    
30     }
31     //停止线程
32     /*
33     int pthread_join(
34         pthread_t    thread    /in
35         void**    ret_val_p    /out    
36 37     第二个参数可以接收任意由pthread_t对象所关联的那个线程产生的返回值。
38     */
39     for(thread=0;thread<thread_count;thread++){
40         pthread_join(thread_handles[thread],NULL);    
41     }
42     pthread_mutex_destroy(&mutex);
43     double pi=4*(double)num_in_circle/(double)n;
44     printf("the esitimate value of pi is %lf\n",pi);
45 }
46 void* compute_pi(void* rank){
47     long long int local_n;
48     local_n=n/thread_count;
49     double x,y,distance_squared;
50     for(long long int i=0;i<local_n;i++){
51         x=(double)rand()/(double)RAND_MAX;
52         y=(double)rand()/(double)RAND_MAX;
53         distance_squared=x*x+y*y;
54         if(distance_squared<=1){
55             pthread_mutex_lock(&mutex);
56             num_in_circle++;
57             pthread_mutex_unlock(&mutex);    
58         }
59     }
60     return NULL;
61 }    

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

运行 ./pth_mete 4 

输入数据结果如下

 下面是OpenMP代码

 1 #include<stdlib.h>
 2 #include<stdio.h>
 3 #include<time.h>
 4 #include<omp.h>
 5 int main(int argc,char** argv){
 6     long long int num_in_cycle=0;
 7     long long int num_point;
 8     int thread_count;
 9     thread_count=strtol(argv[1],NULL,10);
10     printf("please input the number of point\n");
11     scanf("%lld",&num_point);
12     srand(time(NULL));
13     double x,y,distance_point;
14     long long int i;
15     #pragma omp parallel for num_threads(thread_count) default(none) \
16         reduction(+:num_in_cycle) shared(num_point) private(i,x,y,distance_point)
17     for( i=0;i<num_point;i++){
18         x=(double)rand()/(double)RAND_MAX;
19         y=(double)rand()/(double)RAND_MAX;
20         distance_point=x*x+y*y;
21         if(distance_point<=1){
22             num_in_cycle++;
23         }
24     }
25     double estimate_pi=(double)num_in_cycle/num_point*4;
26     printf("the estimate value of pi is %lf\n",estimate_pi);
27     return 0;
28 }

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

执行 ./omp_mete 4 

结果如下

 

posted @ 2014-11-29 20:56  海上的风  阅读(11487)  评论(0编辑  收藏  举报