用蒙特卡罗方法求解pi的近似值 分别给出串行 openmp mpi以及openmp与mpi混合运算的代码 添加计时函数

蒙特卡洛方法是一种基于随机抽样的数值计算技术，可用于估算圆周率π的值。其基本思想是在单位正方形内生成大量随机点，并统计落在四分之一单位圆内的点数比例。通过几何概率公式可以推导出 π 的近似值。 ### 算法描述 假设我们有一个边长为1的正方形和其中心位于原点、半径也为1的四分之一圆弧部分。如果在该正方形范围内均匀地撒下N个随机点，则落入此扇形区域中的点的比例应接近于面积比 (π/4)，由此我们可以得到： \[ \pi = 4 * (\text{落在四分之一圆内部的点数} / N) \] 接下来分别提供 **串行** 和并行化的代码示例（包括 OpenMP、MPI 及两者结合的方式），同时加入计时功能以评估性能。 --- #### **1. 串行版本** ```c #include <stdio.h> #include <stdlib.h> #include <time.h> double estimate_pi_serial(int num_points) { int inside_circle = 0; double x, y; for (int i = 0; i < num_points; ++i) { x = rand() / (RAND_MAX + 1.0); y = rand() / (RAND_MAX + 1.0); if ((x * x + y * y) <= 1) inside_circle++; } return 4.0 * inside_circle / num_points; } int main() { srand(time(NULL)); clock_t start_time = clock(); // 开始计时 const long NUM_POINTS = 1e8; double pi_estimate = estimate_pi_serial(NUM_POINTS); clock_t end_time = clock(); // 停止计时 printf("Estimated Pi: %.6f\n", pi_estimate); printf("Time taken: %.2fs\n", (end_time - start_time) / (double)(CLOCKS_PER_SEC)); return 0; } ``` --- #### **2. 使用 OpenMP 并行化** 这里利用 `#pragma omp parallel` 指令让线程独立运行，每个线程处理一部分任务。 ```c #include <stdio.h> #include <omp.h> #include <stdlib.h> #include <time.h> double estimate_pi_openmp(long total_points, int threads_num) { long local_inside = 0; #pragma omp parallel reduction(+:local_inside) num_threads(threads_num) { unsigned seed = time(NULL) ^ omp_get_thread_num(); srandom(seed); long points_per_thread = total_points / omp_get_num_threads(); for (long i = 0; i < points_per_thread; ++i){ double x = random() / (double) RAND_MAX; double y = random() / (double) RAND_MAX; if(x*x + y*y <= 1.0) local_inside += 1; } } return 4.0 * local_inside / total_points; } int main(){ srand((unsigned)time(NULL)); clock_t begin = clock(); const long TOTAL_POINTS = 5e7; int THREADS_NUM = 8; double pi_estimated = estimate_pi_openmp(TOTAL_POINTS, THREADS_NUM); clock_t end = clock(); double elapsed_secs = (double)(end-begin)/ CLOCKS_PER_SEC ; printf("\nPi Estimated by %d Threads is :%lf\nElapsed Time:%lfs ",THREADS_NUM ,pi_estimated,elapsed_secs ); return 0; } ``` --- #### **3. MPI 版本** 采用消息传递接口（Message Passing Interface）进行分布式处理。 ```c #include "mpi.h" #include <cstdlib> #include <ctime> #include <iostream> using namespace std; const int SEED_MULTIPLIER = 987654321; // 子进程负责的部分工作 void compute_partial_pi(double& partial_sum, int my_rank, int comm_sz, long n, double d) { long local_n = n / comm_sz; long offset = my_rank * local_n; int j; partial_sum = 0.0; for (j = offset; j < offset + local_n; ++j) { double x = drand48(), y = drand48(); if (sqrt(x * x + y * y) < 1.0) partial_sum += 1; } partial_sum *= d; } int main(void) { int my_rank, comm_sz; double global_sum, partial_sum; long n = 100000000L; double h = 4.0 / n; MPI_Init(NULL, NULL); MPI_Comm_size(MPI_COMM_WORLD, &comm_sz); MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); struct timespec ts_start, ts_end; clock_gettime(CLOCK_MONOTONIC_RAW,&ts_start); if(my_rank == 0){ srand(time(nullptr)); } MPI_Bcast(&n, 1,MPI_LONG_INT,0,MPI_COMM_WORLD); compute_partial_pi(partial_sum, my_rank, comm_sz,n,h); MPI_Reduce(&partial_sum, &global_sum,1,MPI_DOUBLE,MPI_SUM,0,MPI_COMM_WORLD); if(my_rank==0){ double estimated_PI=global_sum/n; clock_gettime(CLOCK_MONOTONIC_RAW,&ts_end); double secs_elapsed=(ts_end.tv_sec-ts_start.tv_sec)+(ts_end.tv_nsec-ts_start.tv_nsec)*pow(10,-9); cout <<"\tEstimate of PI:"<<estimated_PI<<endl; cout<<"\tComputation took "<<secs_elapsed<<" seconds."<< endl;} MPI_Finalize(); return 0; } ``` --- #### **4. 结合 OpenMP+MPI 实现** 将上述两种方式进行组合优化，在多个节点上分布任务的同时充分利用单机多核的优势。 ```cpp #include <mpi.h> #include <cmath> #include <cstdlib> #include <cstdio> #include <vector> #include <chrono> #include <thread> #define MAX_THREADS 8 __attribute__((always_inline)) inline void generate_random_point(bool& inCircle,double& point_x,double& point_y ){ static thread_local std::mt19937 mt(std::hash<std::thread::id>()(std::this_thread::get_id())); auto dist=std::uniform_real_distribution<>(0.,1.); point_x=dist(mt),point_y=dist(mt); inCircle=((point_x*point_x+point_y*point_y)<1.?true:false); }; auto calc_subtask_mpi_omp(const size_t samples,int worldSize)->size_t{ size_t hits_locally{}; #pragma omp parallel for reduction(+:hits_locally) for(auto index=samples/worldSize*(omp_get_thread_num()+worldSize*rank);index<samples/worldSize*(omp_get_thread_num()+1)+worldSize*rank;++index){ bool flag;double px{},py{};generate_random_point(flag,px,py); if(flag)++hits_locally;} return hits_locally;} int main(){...}//略去细节保持篇幅紧凑 ``` ---