MPI并行化多体问题源码

//mpi_nbody_basic.c，MPI 基本算法
#include 
#include 
#include 
#include 
#include 

#define OUTPUT                          // 要求输出结果
#define DEBUG                           // debug 模式，每个函数给出更多的输出
#define DIM     2                       // 二维系统
#define X       0                       // X 坐标
#define Y       1                       // Y 坐标
typedef double vect_t[DIM];             // 向量数据类型
const double G = 6.673e-11;             // 万有引力常量
int my_rank comm_sz;                   // 进程编号和总进程数
MPI_Datatype vect_mpi_t;                // 使用的派生数据类型
vect_t *vel = NULL;                     // 全局颗粒速度，用于 0 号进程的输出

void Usage（char* prog_name）// 输入说明
{
	fprintf（stderr “usage: mpiexec -n  %s\n“ prog_name）;
	fprintf（stderr “    \n“）;
	fprintf（stderr “    ‘g‘: inite condition by random\n“）;
	fprintf（stderr “    ‘i‘: inite condition from stdin\n“）;
	exit（0）;
}

void Get_args（int argc char* argv[] int* n_p int* n_steps_p double* delta_t_p int* output_freq_p char* g_i_p）// 获取参数信息
{                                                                     // 所有进程均调用该函数，因为有集合通信，但只有 0 号进程处理参数
	if （my_rank == 0）
	{
		if （argc != 6）
			Usage（argv[0]）;
		*n_p = strtol（argv[1] NULL 10）;
		*n_steps_p = strtol（argv[2] NULL 10）;
		*delta_t_p = strtod（argv[3] NULL）;
		*output_freq_p = strtol（argv[4] NULL 10）;
		*g_i_p = argv[5][0];
		if （*n_p <= 0 || *n_p % comm_sz || *n_steps_p < 0 || *delta_t_p <= 0 || *g_i_p != ‘g‘ && *g_i_p != ‘i‘）// 不合要求的输入情况
		{
			printf（“haha\n“）;
			if （my_rank == 0）
				Usage（argv[0]）;
			MPI_Finalize（）;
			exit（0）;
		}
	}
	MPI_Bcast（n_p 1 MPI_INT 0 MPI_COMM_WORLD）;
	MPI_Bcast（n_steps_p 1 MPI_INT 0 MPI_COMM_WORLD）;
	MPI_Bcast（delta_t_p 1 MPI_DOUBLE 0 MPI_COMM_WORLD）;
	MPI_Bcast（output_freq_p 1 MPI_INT 0 MPI_COMM_WORLD）;
	MPI_Bcast（g_i_p 1 MPI_CHAR 0 MPI_COMM_WORLD）;
#   ifdef DEBUG// 确认各进程中的参数情况
	printf（“Get_args rank%2d n %d n_steps %d delta_t %e output_freq %d g_i %c\n“
		my_rank *n_p *n_steps_p *delta_t_p *output_freq_p *g_i_p）;
	fflush（stdout）;
#   endif
}

void Gen_init_cond（double masses[] vect_t pos[] vect_t loc_vel[] int n int loc_n）// 自动生成初始条件，所有进程均调用该函数，因为有集合通信
{                                                          // 生成的颗粒位于原点和 X 正半轴，速度大小相等，方向平行于 Y 轴，交错向上下
	const double mass = 5.0e24 gap = 1.0e5 speed = 3.0e4;// 使用了地球的质量和公转速度
	if （my_rank == 0）
	{
		// srand（2）;// 使用随机方向和速度大小，下同
		for （int i = 0; i < n; i++）
		{
			masses[i] = mass;
			pos[i][X] = i * gap;
			pos[i][Y] = 0.0;
			vel[i][X] = 0.0;
			// vel[i][Y] = speed * （2 * rand（） / （double）RAND_MAX） - 1）;
			vel[i][Y] = （i % 2） ? -speed : speed;
		}
	}
	// 同步质量，位置信息，分发速度信息
	MPI_Bcast（masses n MPI_DOUBLE 0 MPI_COMM_WORLD）;
	MPI_Bcast（pos n vect_mpi_t 0 MPI_COMM_WORLD）;
	MPI_Scatter（vel loc_n vect_mpi_t loc_vel loc_n vect_mpi_t 0 MPI_COMM_W