基于C语言实现超级计算机原理与操作

♻️ 资源大小6.75MB➡️资源下载https://download.csdn.net/download/s1t16/87453219超级计算机原理与操作一、实验题目nbody 问题输入文件为 nbody.txt输入每一行为一个 body每列分别是质量x 轴位置y 轴位置z 轴位置x 轴速度y 轴速度z 轴速度输出文件格式与输入一致为 20 轮迭代后的结果相关参数设置为 dT0.005G1迭代次数为 20tsp 问题任选一个数据集大规模输入文件为 tsp.txt每一行为每个城市的坐标小规模输入文件为 tsp2.txt文件内是城市之间的距离矩阵补充说明可以自行调整数据规模例如删减至 10 个城市提交时需要将输入文件也一起打包输出文件为经过城市的顺序从 0 开始编号本次实验将求解 nbody 问题实现一个串行版本和 MPI、OpenMP、pthread 三种版本的解法。二、实验内容nbody 问题给定一系列物体的初始位置坐标、初始速度求出在物体间引力作用下一段时间后它们的位置坐标、速度。一个 nbody 问题的求解器通过模拟粒子的行为求解 nbody 问题。输入物体数量 n、物体的质量 masses、物体的初始位置、物体的初始速度、经过时间 dT 输出dT 后物体的位置、速度由牛顿的万有引力定律可以计算出两个物体之间的引力对任意一个物体它和其余 n-1 个物体间的引力之和即它所受到的合力据此由牛顿第二定律可得到该物体的加速度由加速度可得到它下一时刻的速度由物体的起始位置和速度可得它下一时刻的速度。物体间作用力的计算公式原子受到的作用力之和串行版本伪代码如下///////////////////////////////////////////////////////////////////////////// // Declare variables // SCALE 1024 // DIM 3 double masses[SCALE]; double pos[SCALE][DIM]; double vel[SCALE][DIM]; double force_qk[DIM]; double forces[SCALE][DIM]; procedure serial solver Get input data get time start for each timestep do // Initialize forces to 0 for each particle q do forces[q] 0 end // Compute total forces on q for each particle q do for each particle kq do Compute force_qk(force between q and k) forces[q] force_qk forces[k] - force_qk end end // Compute position and velocity of q for each particle q do Compute vel[q] Compute pos[q] end end get time stop Print positions and velocities Print time stop-start end procedure解释get time 部分通过 GET_TIME 函数获取当前时间并返回到变量中。# define GET_TIME(now) { \ struct timeval t; \ gettimeofday(t, NULL); \ now t.tv_sec t.tv_usec/1000000.0; \ }get input 读取输入文件 nbody.txt 的信息并存储在 masses、pos、vel 中。FILE* data fopen(nbody.txt, rb); if(dataNULL) fprintf(stderr, Can not open data file.\n); for(i0; iSCALE; i) { fscanf(data, %lf %lf %lf %lf %lf %lf %lf, masses[i], pos[i][X], pos[i][ , pos[i][Z], vel[i][X], vel[i][Y], vel[i][Z]); } fclose(data);计算物体所受合力的具体过程如下for(q0; qSCALE; q) { for(kq1; kSCALE; k) { x_diff pos[q][X] - pos[k][X]; y_diff pos[q][Y] - pos[k][Y]; z_diff pos[q][Z] - pos[k][Z]; dist sqrt(x_diff*x_diff y_diff*y_diff z_diff*z_diff); dist_cubed dist*dist*dist; force_qk[X] G*masses[q]*masses[k]/dist_cubed*x_diff; force_qk[Y] G*masses[q]*masses[k]/dist_cubed*y_diff; force_qk[Z] G*masses[q]*masses[k]/dist_cubed*z_diff; forces[q][X] - force_qk[X]; forces[q][Y] - force_qk[Y]; forces[q][Z] - force_qk[Z]; forces[k][X] force_qk[X]; forces[k][Y] force_qk[Y]; forces[k][Z] force_qk[Z]; } }计算物体的位置、速度的具体过程如下for(q0; qSCALE; q) { vel[q][X] dT/masses[q]*forces[q][X]; vel[q][Y] dT/masses[q]*forces[q][Y]; vel[q][Z] dT/masses[q]*forces[q][Z]; pos[q][X] vel[q][X]*dT; pos[q][Y] vel[q][Y]*dT; pos[q][Z] vel[q][Z]*dT; }OpenMP 版本伪代码如下///////////////////////////////////////////////////////////////////////////// // Declare shared variables // SCALE 1024 // DIM 3 double masses[SCALE]; double pos[SCALE][DIM]; double vel[SCALE][DIM]; double forces[SCALE][DIM]; procedure openmp solver Get thread_count double loc_forces[thread_count][SCALE][DIM]; Get input data get time start for each timestep do // Initialize loc_forces to 0 # pragma omp for for each particle q do for each thread do loc_forces[thread][q] 0 end end // Compute local forces on q # pragma omp for for each particle q do for each particle kq do Compute force_qk(force between q and k) loc_forces[q] force_qk loc_forces[k] - force_qk end end // Gather all the local forces # pragma omp for for each particle q do forces[q] 0 for each thread do forces[q] loc_forces[thread][q] end end // Compute position and velocity of q # pragma omp for for each particle q do Compute vel[q] Compute pos[q] end end get time stop Print positions and velocities Print time stop-start end procedure解释为了防止竞争用 loc_forces 存储各个线程的计算结果最后再汇总到 forces 中。为了减少不必要的空间浪费根据输入的线程数申请 loc_forces 的空间。if (argc ! 2) Usage(argv[0]); Get_args(argv, thread_count); //获取线程数 double loc_forces[thread_count][SCALE][DIM];获取输入的过程和并行版本相同。初始化 loc_forces 的过程如下# pragma omp for for(q0; qSCALE; q){ for(thread0; threadthread_count; thread) loc_forces[thread][q][X] loc_forces[thread][q][Y] loc_forces[th read][q][Z] 0; }对于计算 loc_forces 的过程中使用到的变量 q, k, x_diff, y_diff, z_diff, dist, dist_cubed 等若将它们用作共享变量中将会导致竞争冲突使得计算结果出错因此需要把它们的声明放在各个线程的内部。# pragma omp parallel num_threads(thread_count) { int q, k; int my_rank omp_get_thread_num(); double x_diff, y_diff, z_diff, dist, dist_cubed; double force_qk[DIM]; # pragma omp for for(q0; qSCALE; q){ for(kq1; kSCALE; k){ x_diff pos[q][X] - pos[k][X]; y_diff pos[q][Y] - pos[k][Y]; z_diff pos[q][Z] - pos[k][Z]; dist sqrt(x_diff*x_diff y_diff*y_diff z_diff*z_diff); dist_cubed dist*dist*dist; force_qk[X] G*masses[q]*masses[k]/dist_cubed*x_diff; force_qk[Y] G*masses[q]*masses[k]/dist_cubed*y_diff; force_qk[Z] G*masses[q]*masses[k]/dist_cubed*z_diff; loc_forces[my_rank][q][X] - force_qk[X]; loc_forces[my_rank][q][Y] - force_qk[Y]; loc_forces[my_rank][q][Z] - force_qk[Z]; loc_forces[my_rank][k][X] force_qk[X]; loc_forces[my_rank][k][Y] force_qk[Y]; loc_forces[my_rank][k][Z] force_qk[Z]; } } }对每个物体 q将各个线程计算出的对 q 的作用力汇总# pragma omp for for(q0; qSCALE; q){ forces[q][X] forces[q][Y] forces[q][Z] 0; for(thread0; threadthread_count; thread){ forces[q][X] loc_forces[thread][q][X]; forces[q][Y] loc_forces[thread][q][Y]; forces[q][Z] loc_forces[thread][q][Z]; } }计算物体的位置、速度的具体过程如下# pragma omp for for(q0; qSCALE; q){ vel[q][X] dT/masses[q]*forces[q][X]; vel[q][Y] dT/masses[q]*forces[q][Y]; vel[q][Z] dT/masses[q]*forces[q][Z]; pos[q][X] vel[q][X]*dT; pos[q][Y] vel[q][Y]*dT; pos[q][Z] vel[q][Z]*dT; }Pthread 版本伪代码如下//////////////////////////////////////////////////////////////// // Declare shared variables // SCALE 1024 // DIM 3 int thread_count; double masses[SCALE]; double pos[SCALE][DIM]; double vel[SCALE][DIM]; double force_qk[DIM]; double forces[SCALE][DIM]; double loc_forces[MAX_THREADS][SCALE][DIM]; procedure pthread solver Get thread_count Malloc space of thread_handles Get input data get time start for each timestep do // Initialize forces to 0 for each thread do pthread_create thread_initForces end for each thread do pthread_join end // Compute local forces on q for each thread do pthread_create thread_compForce end for each thread do pthread_join end // Gather all the local forces for each thread do pthread_create thread_reducForce end for each thread do pthread_join end // Compute position and velocity of q for each thread do pthread_create thread_compPosAndVel end for each thread do pthread_join end end get time stop Print positions and velocities Print time stop-start end procedure function thread_initForces input: my_rank, loc_forces output: set loc_forces[my_rank] to 0 for each particle q do loc_forces[my_rank][q] 0 end end function function thread_compForce input: my_rank, masses, pos output: loc_forces my_n SCALE/thread_count my_first_q my_rank*my_n my_last_q my_first_qmy_n for each particle q in [my_first_q, my_last_q] do for each particle kq do Compute force_qk (force between q and k) loc_forces[q] force_qk loc_forces[k] - force_qk end end end function function thread_reducForce input: my_rank, loc_forces output: forces my_n SCALE/thread_count my_first_q my_rank*my_n my_last_q my_first_qmy_n for each particle q in [my_first_q, my_last_q] do forces[q] 0 for each thread do forces[q] loc_forces[thread][q] end end end function function thread_compPosAndVel input: my_rank, forces, pos, vel output: updated pos and vel my_n SCALE/thread_count my_first_q my_rank*my_n my_last_q my_first_qmy_n for each particle q in [my_first_q, my_last_q] do Compute vel[q] Compute pos[q] end end function解释各个线程函数声明如下void* thread_initForces(void* args_p);void* thread_compForce(void* args_p);void* thread_reducForce(void* args_p);void* thread_compPosAndVel(void* args_p);对于每一个线程它们执行时的输入除了 my_rank 都是共享的因此只需要传递参数 thread对应线程的 my_rank即可。对线程函数 thread_initForces任务分配方式是每个线程处理自己被分配的 loc_forces 的部分。对线程函数 thread_compForce、thread_reducForce、thread_compPosAndVel任务分配方式是将物体按块划分到各个线程去处理每个线程处理自己范围内的物体my_first_q, my_last_q]。每次调用线程函数后都要使用 pthread_join 设置路障因为各个步骤间有先后关系必须等待所有线程都得出运算结果才能进行下一步的运算。其余部分和 OpenMP 版本大致相同。MPI 版本由于进程间不共享内存空间因此 MPI 版本的解法和 OpenMP、Pthread 大不相同。每个进程都需要保存物体的质量 masses。每个进程处理 loc_n SCALE/comm_sz 个物体使用 loc_pos 保存进程所处理的物体的位置使用 loc_vel 保存进程所处理的物体的速度使用 loc_forces 存储进程所处理的物体受到的合力使用 tmp_pos 和 tmp_forces 暂时存储目标物体的质量和受到力。为了高效地使用 MPI 的库函数使用一维数组来存储物体的位置、速度信息在前述解法中都是使用二维数组来存储多维度的位置、速度信息例如 double pos[SCALE][DIM]物体 q 在 X 轴的坐标为 pos[q][X]在 MPI 版本中将使用一维数组来存储例如double* pos;pos malloc(SCALE*DIM*sizeof(double));物体 q 在 X 轴的坐标为 pos[q*DIMX]使用这种方法的存储的好处是可以使用MPI_Send、MPI_Recv 、MPI_Scatter等MPI库函数因为这些函数只支持一维数组的操作。进程间呈环形传递信息。各个进程一开始时只持有自己所处理的物体的位置信息在计算完进程内的物体间作用力后将自己所持有的物体的位置信息及其所受的力传递给它的邻居进程同时从它的另一个邻居进程获取它所持有的物体的位置信息及其所受的力再计算自己所处理的物体和接收到的物体的物体间作用力再将接收的物体的位置信息及其所受的力传递给它的邻居进程……最后对每个进程它所处理的物体和其他进程所处理的物体之间的力都被计算了此时 loc_forces 存储的便是自己所处理的物体所受的合力。伪代码如下////////////////////////////////////////////////////////// procedure mpi solver Get comm_sz Get my_rank source (my_rank1)%comm_sz dest (my_rank-1comm_sz)%comm_sz Get input data MPI_Scatter pos to every process as loc_pos MPI_Scatter vel to every process as loc_vel MPI_Bcast masses get time start for each timestep do copy loc_pos into tmp_pos loc_forces tmp_forces 0 // Compute forces due to interactions among local particles for each local particle q do for each local particle kq do Compute force_qk (force between q and k) loc_forces force_qk tmp_forces - force_qk end end for phase 1 to comm_sz-1 do Send tmp_pos to dest Send tmp_forces to dest Recv tmp_pos from source Recv tmp_forces from source owner (my_rankphase)%comm_szCompute forces due to interactions among my particles and owner’s particle s end// Recv my particles pos and forces Send tmp_pos to dest Send tmp_forces to dest Recv tmp_pos from source Recv tmp_forces from source // Add tmp_forces to loc_forces for each local particle q do loc_forces tmp_forces end for each local particle do Compute loc_vel[q] Compute loc_pos[q] end end get time stop if my_rank 0 do Malloc space of pos, vel end MPI_Gather loc_pos to pos MPI_Gather loc_vel to vel if my_rank 0 do Print pos, vel Print time stop-start end end procedure三、实验结果串行版本运行结果如下经分析与参考输出存在的误差 10-7。运行 7 次输出结果相同运行时间不同每次的运行时间和平均运行时间如下OpenMP 版本运行结果如下经分析与参考输出存在的误差 10-7与串行版本的输出完全相同。分别以 1、2、4、8、16、32 个线程运行 OpenMP 版本的求解器每种情况运行 3 次并计算平均运行时间如下所示Pthread 版本运行结果如下经分析与参考输出存在的误差 10-7与串行版本的输出完全相同。分别以 1、2、4、8、16、32 个线程运行 Pthread 版本的求解器每种情况运行 3 次并计算平均运行时间如下所示MPI 版本运行结果如下经分析与参考输出存在的误差 10-7与串行版本的输出完全相同。分别以 2、4、8、16、32 个线程运行 MPI 版本的求解器每种情况运行 3 次并计算平均运行时间如下所示上述结果的运行时间汇总如下线程/进程数串行版本OpenMPPthreadMPI16.987495×10-18.578487×10-18.256597×10-1/2/6.429719×10-16.252744×10-15.457426×10-14/3.966370×10-13.812043×10-13.362603×10-18/2.151856×10-12.636734×10-11.810127×10-116/1.715460×10-12.186649×10-15.137097×10-132/1.294960×10-12.379883×10-1/对各个并行程序若只以 1 个线程/进程运行运行时间将会比普通的串行程序长OpenMP 版本的运行时间随着线程数增加而减短Pthread 版本的运行时间随着线程数增加而减短但当线程数从 16 增加到 32 时运行时间反而增加 MPI 版本的运行时间随着线程数增加而减短但当进程数大于 8 时运行时间反而增加且是十分显著的增加。只以 1 个线程/进程运行并行程序并行程序比串行程序多了申请空间、判断语句等使得运行时间比串行程序要慢在本次实验中OpenMP 版本是在确定线程数之后才申请 loc_forces 的空间而 Pthread 直接以最大线程数申请空间且 Pthread 版本需要一次次的创建线程函数再调度运行这些原因可能造成 Pthread 版本比 OpenMP 版本要慢此外与 Pthread 相比OpenMP 更多的部分是有编译器来执行的编译器所使用的算法可能比我用 Pthread 写的更高效这也是 OpenMP 更快的一个可能原因。MPI 之间的数据共享需要通过消息传递因为 MPI 同步的程序属于不同的进程甚至不同的主机上的不同进程相反由于 OpenMP 和 Pthread 共享内存由于数据不共享每个进程都需要存储大量的相同数据比如 masses我认为这是导致 MPI 版本速度受限的一大原因当进程数为 16 时运行时间大大增加估计是因为分配较大的内存空间拖慢速度此外进程间的通信速度也是影响 MPI 版本的求解器的运行时间的一大原因。