首页 > 代码库 > 蒙特卡洛方法计算圆周率的三种实现-MPI openmp pthread
蒙特卡洛方法计算圆周率的三种实现-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/out30 int count //in31 MPI_Datatype datatype //in32 int source_proc //in33 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 //in57 void* output_data_p //out58 int count //in59 MPI_Datatype datatype //in60 MPI_Op oprtator //in61 int dest_process //in62 MPI_Comm comm //in63 ) 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 //out22 const pthread_attr_t* attr_p23 void* (*start_routine)(void*) //in24 void* arg_p //in 25 )26 第一个参数是一个指针，指向对应的pthread_t对象。注意，pthread_t对象不是pthread_create函数分配的，必须在调用pthread_create函数前就为pthread_t27 对象分配内存空间。第二个参数不用，所以只是函数调用时把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 /in35 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
结果如下
蒙特卡洛方法计算圆周率的三种实现-MPI openmp pthread
声明：以上内容来自用户投稿及互联网公开渠道收集整理发布，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任，若内容有误或涉及侵权可进行投诉： 投诉/举报 工作人员会在5个工作日内联系你，一经查实，本站将立刻删除涉嫌侵权内容。