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

♻️ 资源

大小：6.75MB

➡️资源下载：https://download.csdn.net/download/s1t16/87453219

超级计算机原理与操作

一、实验题目

nbody 问题

输入文件为 nbody.txt，输入每一行为一个 body，每列分别是质量，x 轴位置，y 轴位置，z 轴位置，x 轴速度，y 轴速度，z 轴速度输出文件格式与输入一致，为 20 轮迭代后的结果相关参数设置为 dT=0.005，G=1，迭代次数为 20

• tsp 问题（任选一个数据集）大规模输入文件为 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 k>q 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(data==NULL) fprintf(stderr, "Can not open data file.\n"); for(i=0; i<SCALE; 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(q=0; q<SCALE; q++) { for(k=q+1; k<SCALE; 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(q=0; q<SCALE; 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 k>q 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(q=0; q<SCALE; q++){ for(thread=0; thread<thread_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(q=0; q<SCALE; q++){ for(k=q+1; k<SCALE; 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(q=0; q<SCALE; q++){ forces[q][X] = forces[q][Y] = forces[q][Z] = 0; for(thread=0; thread<thread_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(q=0; q<SCALE; 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_q+my_n for each particle q in [my_first_q, my_last_q] do for each particle k>q 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_q+my_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_q+my_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*DIM+X]，使用这种方法的存储的好处是可以使用MPI_Send、MPI_Recv 、MPI_Scatter等MPI库函数，因为这些函数只支持一维数组的操作。

进程间呈环形传递信息。

各个进程一开始时只持有自己所处理的物体的位置信息，在计算完进程内的物体间作用力后，将自己所持有的物体的位置信息及其所受的力传递给它的邻居进程，同时从它的另一个邻居进程获取它所持有的物体的位置信息及其所受的力，再计算自己所处理的物体和接收到的物体的物体间作用力，再将接收的物体的位置信息及其所受的力传递给它的邻居进程……最后对每个进程，它所处理的物体和其他进程所处理的物体之间的力都被计算了，此时 loc_forces 存储的便是自己所处理的物体所受的合力。

伪代码如下：

////////////////////////////////////////////////////////// procedure mpi solver Get comm_sz Get my_rank source = (my_rank+1)%comm_sz dest = (my_rank-1+comm_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 k>q 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_rank+phase)%comm_sz

Compute 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 次并计算平均运行时间，如下所示：

上述结果的运行时间汇总如下：

线程/进程数	串行版本	OpenMP	Pthread	MPI
1	6.987495×10-1	8.578487×10-1	8.256597×10-1	/
2	/	6.429719×10-1	6.252744×10-1	5.457426×10-1
4	/	3.966370×10-1	3.812043×10-1	3.362603×10-1
8	/	2.151856×10-1	2.636734×10-1	1.810127×10-1
16	/	1.715460×10-1	2.186649×10-1	5.137097×10-1
32	/	1.294960×10-1	2.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 版本的求解器的运行时间的一大原因。