▶ 《并行程序设计导论》第六章中讨论了 n 体问题，分别使用了 MPI，Pthreads，OpenMP 来进行实现，这里是 MPI 的代码，分为基本算法和简化算法（引力计算量为基本算法的一半，但是消息传递较为复杂）
● 基本算法
1 1 //mpi_nbody_basic.c，MPI 基本算法 2 2 #include <stdio.h> 3 3 #include <stdlib.h> 4 4 #include <string.h> 5 5 #include <math.h> 6 6 #include <mpi.h> 7 7 8 8 #define OUTPUT // 要求输出结果 9 9 #define DEBUG // debug 模式，每个函数给出更多的输出 10 10 #define DIM 2 // 二维系统 11 11 #define X 0 // X 坐标 12 12 #define Y 1 // Y 坐标 13 13 typedef double vect_t[DIM]; // 向量数据类型 14 14 const double G = 6.673e-11; // 万有引力常量 15 15 int my_rank, comm_sz; // 进程编号和总进程数 16 16 MPI_Datatype vect_mpi_t; // 使用的派生数据类型 17 17 vect_t *vel = NULL; // 全局颗粒速度，用于 0 号进程的输出 18 18 19 19 void Usage(char* prog_name)// 输入说明 20 20 { 21 21 fprintf(stderr, "usage: mpiexec -n <nProcesses> %s\n", prog_name); 22 22 fprintf(stderr, "<nParticle> <nTimestep> <sizeTimestep> <outputFrequency> <g|i>\n"); 23 23 fprintf(stderr, " 'g': inite condition by random\n"); 24 24 fprintf(stderr, " 'i': inite condition from stdin\n"); 25 25 exit(0); 26 26 } 27 27 28 28 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)// 获取参数信息 29 29 { // 所有进程均调用该函数，因为有集合通信，但只有 0 号进程处理参数 30 30 if (my_rank == 0) 31 31 { 32 32 if (argc != 6) 33 33 Usage(argv[0]); 34 34 *n_p = strtol(argv[1], NULL, 10); 35 35 *n_steps_p = strtol(argv[2], NULL, 10); 36 36 *delta_t_p = strtod(argv[3], NULL); 37 37 *output_freq_p = strtol(argv[4], NULL, 10); 38 38 *g_i_p = argv[5][0]; 39 39 if (*n_p <= 0 || *n_p % comm_sz || *n_steps_p < 0 || *delta_t_p <= 0 || *g_i_p != 'g' && *g_i_p != 'i')// 不合要求的输入情况 40 40 { 41 41 printf("haha\n"); 42 42 if (my_rank == 0) 43 43 Usage(argv[0]); 44 44 MPI_Finalize(); 45 45 exit(0); 46 46 } 47 47 } 48 48 MPI_Bcast(n_p, 1, MPI_INT, 0, MPI_COMM_WORLD); 49 49 MPI_Bcast(n_steps_p, 1, MPI_INT, 0, MPI_COMM_WORLD); 50 50 MPI_Bcast(delta_t_p, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); 51 51 MPI_Bcast(output_freq_p, 1, MPI_INT, 0, MPI_COMM_WORLD); 52 52 MPI_Bcast(g_i_p, 1, MPI_CHAR, 0, MPI_COMM_WORLD); 53 53 # ifdef DEBUG// 确认各进程中的参数情况 54 54 printf("Get_args, rank%2d n %d n_steps %d delta_t %e output_freq %d g_i %c\n", 55 55 my_rank, *n_p, *n_steps_p, *delta_t_p, *output_freq_p, *g_i_p); 56 56 fflush(stdout); 57 57 # endif 58 58 } 59 59 60 60 void Gen_init_cond(double masses[], vect_t pos[], vect_t loc_vel[], int n, int loc_n)// 自动生成初始条件，所有进程均调用该函数，因为有集合通信 61 61 { // 生成的颗粒位于原点和 X 正半轴，速度大小相等，方向平行于 Y 轴，交错向上下 62 62 const double mass = 5.0e24, gap = 1.0e5, speed = 3.0e4;// 使用了地球的质量和公转速度 63 63 if (my_rank == 0) 64 64 { 65 65 // srand(2);// 使用随机方向和速度大小，下同 66 66 for (int i = 0; i < n; i++) 67 67 { 68 68 masses[i] = mass; 69 69 pos[i][X] = i * gap; 70 70 pos[i][Y] = 0.0; 71 71 vel[i][X] = 0.0; 72 72 // vel[i][Y] = speed * (2 * rand() / (double)RAND_MAX) - 1); 73 73 vel[i][Y] = (i % 2) ? -speed : speed; 74 74 } 75 75 } 76 76 // 同步质量，位置信息，分发速度信息 77 77 MPI_Bcast(masses, n, MPI_DOUBLE, 0, MPI_COMM_WORLD); 78 78 MPI_Bcast(pos, n, vect_mpi_t, 0, MPI_COMM_WORLD); 79 79 MPI_Scatter(vel, loc_n, vect_mpi_t, loc_vel, loc_n, vect_mpi_t, 0, MPI_COMM_WORLD); 80 80 # ifdef DEBUG// 确认各进程中第一个颗粒的初始条件 81 81 printf("Gen_init_cond, rank%2d %10.3e %10.3e %10.3e %10.3e %10.3e\n", 82 82 my_rank, masses[0], pos[0][X], pos[0][Y], loc_vel[0][X], loc_vel[0][Y]); 83 83 fflush(stdout); 84 84 # endif 85 85 } 86 86 87 87 void Get_init_cond(double masses[], vect_t pos[], vect_t loc_vel[], int n, int loc_n)// 手工输入初始条件，类似函数 Gen_init_cond() 88 88 { 89 89 if (my_rank == 0) 90 90 { 91 91 printf("For each particle, enter (in order): mass x-coord y-coord x-velocity y-velocity\n"); 92 92 for (int i = 0; i < n; i++) 93 93 { 94 94 scanf_s("%lf", &masses[i]); 95 95 scanf_s("%lf", &pos[i][X]); 96 96 scanf_s("%lf", &pos[i][Y]); 97 97 scanf_s("%lf", &vel[i][X]); 98 98 scanf_s("%lf", &vel[i][Y]); 99 99 } 100100 } 101101 MPI_Bcast(masses, n, MPI_DOUBLE, 0, MPI_COMM_WORLD); 102102 MPI_Bcast(pos, n, vect_mpi_t, 0, MPI_COMM_WORLD); 103103 MPI_Scatter(vel, loc_n, vect_mpi_t, loc_vel, loc_n, vect_mpi_t, 0, MPI_COMM_WORLD); 104104 # ifdef DEBUG 105105 printf("Get_init_cond, rank%2d %10.3e %10.3e %10.3e %10.3e %10.3e\n", 106106 my_rank, masses[0], pos[0][X], pos[0][Y], loc_vel[0][X], loc_vel[0][Y]); 107107 fflush(stdout); 108108 # endif 109109 } 110110 111111 void Output_state(double time, double masses[], vect_t pos[], vect_t loc_vel[], int n, int loc_n)// 输出当前状态 112112 { 113113 MPI_Gather(loc_vel, loc_n, vect_mpi_t, vel, loc_n, vect_mpi_t, 0, MPI_COMM_WORLD);// 从各进程聚集速度信息用于输出 114114 if (my_rank == 0) 115115 { 116116 printf("Output_state, time = %.2f\n", time); 117117 for (int i = 0; i < n; i++) 118118 printf(" %2d %10.3e %10.3e %10.3e %10.3e %10.3e\n", i, masses[i], pos[i][X], pos[i][Y], vel[i][X], vel[i][Y]); 119119 printf("\n"); 120120 fflush(stdout); 121121 } 122122 } 123123 124124 void Compute_force(int loc_part, double masses[], vect_t loc_forces[], vect_t pos[], int n, int loc_n)// 计算颗粒 part 受到的万有引力 125125 { 126126 const int part = my_rank * loc_n + loc_part;// 将局部颗粒编号转化为全局颗粒编号 127127 int k; 128128 vect_t f_part_k; 129129 double len, fact; 130130 for (loc_forces[loc_part][X] = loc_forces[loc_part][Y] = 0.0, k = 0; k < n; k++) 131131 { 132132 if (k != part) 133133 { 134134 f_part_k[X] = pos[part][X] - pos[k][X]; 135135 f_part_k[Y] = pos[part][Y] - pos[k][Y]; 136136 len = sqrt(f_part_k[X] * f_part_k[X] + f_part_k[Y] * f_part_k[Y]); 137137 fact = -G * masses[part] * masses[k] / (len * len * len); 138138 f_part_k[X] *= fact; 139139 f_part_k[Y] *= fact; 140140 loc_forces[loc_part][X] += f_part_k[X]; 141141 loc_forces[loc_part][Y] += f_part_k[Y]; 142142 # ifdef DEBUG// 确认计算结果 143143 printf("Compute_force, rank%2d k%2d> %10.3e %10.3e %10.3e %10.3e\n", my_rank, k, len, fact, f_part_k[X], f_part_k[Y]); 144144 fflush(stdout);// 函数 printf() 中引用 vel 会导致非 0 号进程中断退出，吃了大亏 145145 # endif 146146 } 147147 } 148148 } 149149 150150 void Update_part(int loc_part, double masses[], vect_t loc_forces[], vect_t loc_pos[], vect_t loc_vel[], int n, int loc_n, double delta_t)// 更新颗粒位置 151151 { 152152 const int part = my_rank * loc_n + loc_part; 153153 const double fact = delta_t / masses[part]; 154154 # ifdef DEBUG// 输出在更新前和更新后的数据 155155 printf("Update_part before, part%2d %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n", 156156 part, loc_pos[loc_part][X], loc_pos[loc_part][Y], loc_vel[loc_part][X], loc_vel[loc_part][Y], loc_forces[loc_part][X], loc_forces[loc_part][Y]); 157157 fflush(stdout); 158158 # endif 159159 loc_pos[loc_part][X] += delta_t * loc_vel[loc_part][X]; 160160 loc_pos[loc_part][Y] += delta_t * loc_vel[loc_part][Y]; 161161 loc_vel[loc_part][X] += fact * loc_forces[loc_part][X]; 162162 loc_vel[loc_part][Y] += fact * loc_forces[loc_part][Y]; 163163 # ifdef DEBUG 164164 printf("Update_part after, part%2d %10.3e %10.3e %10.3e %10.3e\n", 165165 part, loc_pos[loc_part][X], loc_pos[loc_part][Y], loc_vel[loc_part][X], loc_vel[loc_part][Y]); 166166 fflush(stdout); 167167 # endif 168168 } 169169 170170 int main(int argc, char* argv[]) 171171 { 172172 int n, loc_n, loc_part; // 颗粒数，每进程颗粒数，当前颗粒（循环变量） 173173 int n_steps, step; // 计算时间步数，当前时间片（循环变量） 174174 double delta_t; // 计算时间步长 175175 int output_freq; // 数据输出频率 176176 double *masses; // 颗粒质量，每个进程都有，一经初始化和同步就不再改变 177177 vect_t *pos, *loc_pos; // 颗粒位置，每个时间片计算完成后需要同步 178178 vect_t *loc_vel; // 颗粒速度，由各进程分开保存，不到输出时不用同步 179179 vect_t *loc_forces; // 各进程的颗粒所受引力 180180 char g_i; // 初始条件选项，g 为自动生成，i 为手工输入 181181 double start, finish; // 计时器 182182 183183 MPI_Init(&argc, &argv); 184184 MPI_Comm_size(MPI_COMM_WORLD, &comm_sz); 185185 MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); 186186 MPI_Type_contiguous(DIM, MPI_DOUBLE, &vect_mpi_t);// 提交需要的派生类型 187187 MPI_Type_commit(&vect_mpi_t); 188188 189189 // 获取参数，初始化数组 190190 Get_args(argc, argv, &n, &n_steps, &delta_t, &output_freq, &g_i); 191191 loc_n = n / comm_sz; // 要求 n % comm_sz == 0 192192 masses = (double*)malloc(n * sizeof(double)); 193193 pos = (vect_t*)malloc(n * sizeof(vect_t)); 194194 loc_pos = pos + my_rank * loc_n; 195195 loc_vel = (vect_t*)malloc(loc_n * sizeof(vect_t)); 196196 loc_forces = (vect_t*)malloc(loc_n * sizeof(vect_t)); 197197 if (my_rank == 0) 198198 vel = (vect_t*)malloc(n * sizeof(vect_t)); 199199 if (g_i == 'g') 200200 Gen_init_cond(masses, pos, loc_vel, n, loc_n); 201201 else 202202 Get_init_cond(masses, pos, loc_vel, n, loc_n); 203203 204204 // 开始计算并计时 205205 if (my_rank == 0) 206206 start = MPI_Wtime(); 207207 # ifdef OUTPUT// 输出出初始态 208208 Output_state(0.0, masses, pos, loc_vel, n, loc_n); 209209 # endif 210210 for (step = 1; step <= n_steps; step++) 211211 { 212212 // 计算每颗粒受力，更新颗粒状态，然后同步颗粒位置 213213 for (loc_part = 0; loc_part < loc_n; Compute_force(loc_part++, masses, loc_forces, pos, n, loc_n)); 214214 for (loc_part = 0; loc_part < loc_n; Update_part(loc_part++, masses, loc_forces, loc_pos, loc_vel, n, loc_n, delta_t)); 215215 MPI_Allgather(MPI_IN_PLACE, loc_n, vect_mpi_t, pos, loc_n, vect_mpi_t, MPI_COMM_WORLD); 216216 # ifdef OUTPUT// 每隔一个输出时间间隔就输出一次 217217 if (step % output_freq == 0) 218218 Output_state(step * delta_t, masses, pos, loc_vel, n, loc_n); 219219 # endif 220220 } 221221 // 报告计时，释放资源 222222 if (my_rank == 0) 223223 { 224224 finish = MPI_Wtime(); 225225 printf("Elapsed time = %e ms\n", (finish - start) * 1000); 226226 free(vel); 227227 } 228228 MPI_Type_free(&vect_mpi_t); 229229 free(masses); 230230 free(pos); 231231 free(loc_vel); 232232 free(loc_forces); 233233 MPI_Finalize(); 234234 return 0; 235235 }
● 输出结果。8 进程 16 体，3 秒，时间步长 1 秒，舍去 debug 输出 1.264884e+00 ms；8 进程 1024 体，3600 秒，时间步长 1 秒，舍去 output 和 debug 输出 1.328894e+04 ms
1D:\Code\MPI\MPIProjectTemp\x64\Debug>mpiexec -n 8 MPIProjectTemp.exe 16 3 1 1 g 2Output_state, time = 0.00 3 0 5.000e+24 0.000e+00 0.000e+00 0.000e+00 3.000e+04 4 1 5.000e+24 1.000e+05 0.000e+00 0.000e+00 -3.000e+04 5 2 5.000e+24 2.000e+05 0.000e+00 0.000e+00 3.000e+04 6 3 5.000e+24 3.000e+05 0.000e+00 0.000e+00 -3.000e+04 7 4 5.000e+24 4.000e+05 0.000e+00 0.000e+00 3.000e+04 8 5 5.000e+24 5.000e+05 0.000e+00 0.000e+00 -3.000e+04 9 6 5.000e+24 6.000e+05 0.000e+00 0.000e+00 3.000e+04 10 7 5.000e+24 7.000e+05 0.000e+00 0.000e+00 -3.000e+04 11 8 5.000e+24 8.000e+05 0.000e+00 0.000e+00 3.000e+04 12 9 5.000e+24 9.000e+05 0.000e+00 0.000e+00 -3.000e+04 13 10 5.000e+24 1.000e+06 0.000e+00 0.000e+00 3.000e+04 14 11 5.000e+24 1.100e+06 0.000e+00 0.000e+00 -3.000e+04 15 12 5.000e+24 1.200e+06 0.000e+00 0.000e+00 3.000e+04 16 13 5.000e+24 1.300e+06 0.000e+00 0.000e+00 -3.000e+04 17 14 5.000e+24 1.400e+06 0.000e+00 0.000e+00 3.000e+04 18 15 5.000e+24 1.500e+06 0.000e+00 0.000e+00 -3.000e+04 19 20Output_state, time = 1.00 21 0 5.000e+24 0.000e+00 3.000e+04 5.273e+04 3.000e+04 22 1 5.000e+24 1.000e+05 -3.000e+04 1.922e+04 -3.000e+04 23 2 5.000e+24 2.000e+05 3.000e+04 1.071e+04 3.000e+04 24 3 5.000e+24 3.000e+05 -3.000e+04 6.802e+03 -3.000e+04 25 4 5.000e+24 4.000e+05 3.000e+04 4.485e+03 3.000e+04 26 5 5.000e+24 5.000e+05 -3.000e+04 2.875e+03 -3.000e+04 27 6 5.000e+24 6.000e+05 3.000e+04 1.614e+03 3.000e+04 28 7 5.000e+24 7.000e+05 -3.000e+04 5.213e+02 -3.000e+04 29 8 5.000e+24 8.000e+05 3.000e+04 -5.213e+02 3.000e+04 30 9 5.000e+24 9.000e+05 -3.000e+04 -1.614e+03 -3.000e+04 31 10 5.000e+24 1.000e+06 3.000e+04 -2.875e+03 3.000e+04 32 11 5.000e+24 1.100e+06 -3.000e+04 -4.485e+03 -3.000e+04 33 12 5.000e+24 1.200e+06 3.000e+04 -6.802e+03 3.000e+04 34 13 5.000e+24 1.300e+06 -3.000e+04 -1.071e+04 -3.000e+04 35 14 5.000e+24 1.400e+06 3.000e+04 -1.922e+04 3.000e+04 36 15 5.000e+24 1.500e+06 -3.000e+04 -5.273e+04 -3.000e+04 37 38Output_state, time = 2.00 39 0 5.000e+24 5.273e+04 6.000e+04 9.288e+04 1.641e+04 40 1 5.000e+24 1.192e+05 -6.000e+04 3.818e+04 -3.791e+03 41 2 5.000e+24 2.107e+05 6.000e+04 2.116e+04 3.791e+03 42 3 5.000e+24 3.068e+05 -6.000e+04 1.356e+04 -3.101e+03 43 4 5.000e+24 4.045e+05 6.000e+04 8.930e+03 3.101e+03 44 5 5.000e+24 5.029e+05 -6.000e+04 5.739e+03 -2.959e+03 45 6 5.000e+24 6.016e+05 6.000e+04 3.218e+03 2.959e+03 46 7 5.000e+24 7.005e+05 -6.000e+04 1.043e+03 -2.929e+03 47 8 5.000e+24 7.995e+05 6.000e+04 -1.043e+03 2.929e+03 48 9 5.000e+24 8.984e+05 -6.000e+04 -3.218e+03 -2.959e+03 49 10 5.000e+24 9.971e+05 6.000e+04 -5.739e+03 2.959e+03 50 11 5.000e+24 1.096e+06 -6.000e+04 -8.930e+03 -3.101e+03 51 12 5.000e+24 1.193e+06 6.000e+04 -1.356e+04 3.101e+03 52 13 5.000e+24 1.289e+06 -6.000e+04 -2.116e+04 -3.791e+03 53 14 5.000e+24 1.381e+06 6.000e+04 -3.818e+04 3.791e+03 54 15 5.000e+24 1.447e+06 -6.000e+04 -9.288e+04 -1.641e+04 55 56Output_state, time = 3.00 57 0 5.000e+24 1.456e+05 7.641e+04 1.273e+05 -1.575e+03 58 1 5.000e+24 1.574e+05 -6.379e+04 5.867e+04 2.528e+04 59 2 5.000e+24 2.319e+05 6.379e+04 2.685e+04 -2.069e+04 60 3 5.000e+24 3.204e+05 -6.310e+04 1.884e+04 2.228e+04 61 4 5.000e+24 4.134e+05 6.310e+04 1.235e+04 -2.151e+04 62 5 5.000e+24 5.086e+05 -6.296e+04 8.111e+03 2.179e+04 63 6 5.000e+24 6.048e+05 6.296e+04 4.537e+03 -2.163e+04 64 7 5.000e+24 7.016e+05 -6.293e+04 1.493e+03 2.169e+04 65 8 5.000e+24 7.984e+05 6.293e+04 -1.493e+03 -2.169e+04 66 9 5.000e+24 8.952e+05 -6.296e+04 -4.537e+03 2.163e+04 67 10 5.000e+24 9.914e+05 6.296e+04 -8.111e+03 -2.179e+04 68 11 5.000e+24 1.087e+06 -6.310e+04 -1.235e+04 2.151e+04 69 12 5.000e+24 1.180e+06 6.310e+04 -1.884e+04 -2.228e+04 70 13 5.000e+24 1.268e+06 -6.379e+04 -2.685e+04 2.069e+04 71 14 5.000e+24 1.343e+06 6.379e+04 -5.867e+04 -2.528e+04 72 15 5.000e+24 1.354e+06 -7.641e+04 -1.273e+05 1.575e+03 73 74Elapsed time = 1.264884e+00 ms
● 简化算法
1 1 // mpi_nbody_red.c，MPI 简化算法，与 mppi_nbody_basic.c 相同的地方去掉了注释 2 2 #include <stdio.h> 3 3 #include <stdlib.h> 4 4 #include <string.h> 5 5 #include <math.h> 6 6 #include <mpi.h> 7 7 8 8 #define OUTPUT 9 9 #define DEBUG 10 10 #define DIM 2 11 11 #define X 0 12 12 #define Y 1 13 13 typedef double vect_t[DIM]; 14 14 const double G = 6.673e-11; 15 15 int my_rank, comm_sz; 16 16 MPI_Datatype vect_mpi_t; 17 17 MPI_Datatype cyclic_mpi_t; // 环形传输的数据类型 18 18 vect_t *vel = NULL; 19 19 vect_t *pos = NULL; // 位置信息变成了全局变量 20 20 21 21 void Usage(char* prog_name) 22 22 { 23 23 fprintf(stderr, "usage: mpiexec -n <nProcesses> %s\n", prog_name); 24 24 fprintf(stderr, "<nParticle> <nTimestep> <sizeTimestep> <outputFrequency> <g|i>\n"); 25 25 fprintf(stderr, " 'g': inite condition by random\n"); 26 26 fprintf(stderr, " 'i': inite condition from stdin\n"); 27 27 exit(0); 28 28 } 29 29 30 30 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) 31 31 { 32 32 if (my_rank == 0) 33 33 { 34 34 if (argc != 6) 35 35 Usage(argv[0]); 36 36 *n_p = strtol(argv[1], NULL, 10); 37 37 *n_steps_p = strtol(argv[2], NULL, 10); 38 38 *delta_t_p = strtod(argv[3], NULL); 39 39 *output_freq_p = strtol(argv[4], NULL, 10); 40 40 *g_i_p = argv[5][0]; 41 41 if (*n_p <= 0 || *n_p % comm_sz || *n_steps_p < 0 || *delta_t_p <= 0 || *g_i_p != 'g' && *g_i_p != 'i') 42 42 { 43 43 printf("haha\n"); 44 44 if (my_rank == 0) 45 45 Usage(argv[0]); 46 46 MPI_Finalize(); 47 47 exit(0); 48 48 } 49 49 } 50 50 MPI_Bcast(n_p, 1, MPI_INT, 0, MPI_COMM_WORLD); 51 51 MPI_Bcast(n_steps_p, 1, MPI_INT, 0, MPI_COMM_WORLD); 52 52 MPI_Bcast(delta_t_p, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); 53 53 MPI_Bcast(output_freq_p, 1, MPI_INT, 0, MPI_COMM_WORLD); 54 54 MPI_Bcast(g_i_p, 1, MPI_CHAR, 0, MPI_COMM_WORLD); 55 55 # ifdef DEBUG 56 56 printf("Get_args, rank%2d n %d n_steps %d delta_t %e output_freq %d g_i %c\n", 57 57 my_rank, *n_p, *n_steps_p, *delta_t_p, *output_freq_p, *g_i_p); 58 58 fflush(stdout); 59 59 # endif 60 60 } 61 61 62 62 void Build_cyclic_mpi_type(int loc_n)// 生成大小为 loc_n 的循环分配数据类型 63 63 { 64 64 MPI_Datatype temp_mpi_t; 65 65 MPI_Aint lb, extent; 66 66 MPI_Type_vector(loc_n, 1, comm_sz, vect_mpi_t, &temp_mpi_t);// 将跨度为 comm_sz 的 loc_n 个 “连续 2 个 double” 封装为 temp_mpi_t 67 67 MPI_Type_get_extent(vect_mpi_t, &lb, &extent); 68 68 MPI_Type_create_resized(temp_mpi_t, lb, extent, &cyclic_mpi_t); 69 69 MPI_Type_commit(&cyclic_mpi_t); 70 70 } 71 71 72 72 void Gen_init_cond(double masses[], vect_t loc_pos[], vect_t loc_vel[], int n, int loc_n) 73 73 { 74 74 const double mass = 5.0e24, gap = 1.0e5, speed = 3.0e4; 75 75 if (my_rank == 0) 76 76 { 77 77 // srand(2); 78 78 for (int i = 0; i < n; i++) 79 79 { 80 80 masses[i] = mass; 81 81 pos[i][X] = i * gap; 82 82 pos[i][Y] = 0.0; 83 83 vel[i][X] = 0.0; 84 84 // vel[i][Y] = speed * (2 * rand() / (double)RAND_MAX) - 1); 85 85 vel[i][Y] = (i % 2) ? -speed : speed; 86 86 } 87 87 } 88 88 MPI_Bcast(masses, n, MPI_DOUBLE, 0, MPI_COMM_WORLD); 89 89 MPI_Scatter(pos, 1, cyclic_mpi_t, loc_pos, loc_n, vect_mpi_t, 0, MPI_COMM_WORLD);// loc_pos 和 loc_vel 需要分别分发 90 90 MPI_Scatter(vel, 1, cyclic_mpi_t, loc_vel, loc_n, vect_mpi_t, 0, MPI_COMM_WORLD); 91 91 # ifdef DEBUG 92 92 printf("Gen_init_cond, rank%2d %10.3e %10.3e %10.3e %10.3e %10.3e\n", 93 93 my_rank, masses[0], loc_pos[0][X], loc_pos[0][Y], loc_vel[0][X], loc_vel[0][Y]); 94 94 fflush(stdout); 95 95 # endif 96 96 } 97 97 98 98 void Get_init_cond(double masses[], vect_t loc_pos[], vect_t loc_vel[], int n, int loc_n) 99 99 { 100100 if (my_rank == 0) 101101 { 102102 printf("For each particle, enter (in order): mass x-coord y-coord x-velocity y-velocity\n"); 103103 for (int i = 0; i < n; i++) 104104 { 105105 scanf_s("%lf", &masses[i]); 106106 scanf_s("%lf", &pos[i][X]); 107107 scanf_s("%lf", &pos[i][Y]); 108108 scanf_s("%lf", &vel[i][X]); 109109 scanf_s("%lf", &vel[i][Y]); 110110 } 111111 } 112112 MPI_Bcast(masses, n, MPI_DOUBLE, 0, MPI_COMM_WORLD); 113113 MPI_Scatter(pos, 1, cyclic_mpi_t, loc_pos, loc_n, vect_mpi_t, 0, MPI_COMM_WORLD);// loc_pos 和 loc_vel 需要分别分发 114114 MPI_Scatter(vel, 1, cyclic_mpi_t, loc_vel, loc_n, vect_mpi_t, 0, MPI_COMM_WORLD); 115115 # ifdef DEBUG 116116 printf("Get_init_cond, rank%2d %10.3e %10.3e %10.3e %10.3e %10.3e\n", 117117 my_rank, masses[0], loc_pos[0][X], loc_pos[0][Y], loc_vel[0][X], loc_vel[0][Y]); 118118 fflush(stdout); 119119 # endif 120120 } 121121 122122 void Output_state(double time, double masses[], vect_t loc_pos[], vect_t loc_vel[], int n, int loc_n) 123123 { 124124 MPI_Gather(loc_pos, loc_n, vect_mpi_t, pos, 1, cyclic_mpi_t, 0, MPI_COMM_WORLD);// loc_pos 和 loc_vel 需要分别聚集 125125 MPI_Gather(loc_vel, loc_n, vect_mpi_t, vel, 1, cyclic_mpi_t, 0, MPI_COMM_WORLD); 126126 if (my_rank == 0) 127127 { 128128 printf("Output_state, time = %.2f\n", time); 129129 for (int i = 0; i < n; i++) 130130 printf(" %2d %10.3e %10.3e %10.3e %10.3e %10.3e\n", i, masses[i], pos[i][X], pos[i][Y], vel[i][X], vel[i][Y]); 131131 printf("\n"); 132132 fflush(stdout); 133133 } 134134 } 135135 136136 int First_index(int gbl1, int rk1, int rk2, int proc_count) 137137 { 138138 return gbl1 + (rk2 - rk1) + (rk1 < rk2 ? 0 : proc_count); 139139 } 140140 141141 int Local_to_global(int loc_part, int proc_rk, int proc_count)// 进程局部编号转全局编号 142142 { 143143 return loc_part * proc_count + proc_rk; 144144 } 145145 146146 int Global_to_local(int gbl_part, int proc_rk, int proc_count)// 全局编号转进程局部编号 147147 { 148148 return (gbl_part - proc_rk) / proc_count; 149149 } 150150 151151 void Compute_force_pair(double m1, double m2, vect_t pos1, vect_t pos2, vect_t force1, vect_t force2)// 计算两个颗粒之间的引力 152152 { 153153 const double mg = -G * m1 * m2; 154154 vect_t f_part_k; 155155 double len, fact; 156156 157157 f_part_k[X] = pos1[X] - pos2[X]; 158158 f_part_k[Y] = pos1[Y] - pos2[Y]; 159159 len = sqrt(f_part_k[X] * f_part_k[X] + f_part_k[Y] * f_part_k[Y]); 160160 fact = mg / (len * len * len); 161161 f_part_k[X] *= fact; 162162 f_part_k[Y] *= fact; 163163 force1[X] += f_part_k[X]; 164164 force1[Y] += f_part_k[Y]; 165165 force2[X] -= f_part_k[X]; 166166 force2[Y] -= f_part_k[Y]; 167167 } 168168 169169 void Compute_proc_forces(double masses[], vect_t tmp_data[], vect_t loc_forces[], vect_t pos1[], int loc_n1, int rk1, int loc_n2, int rk2, int n, int p) 170170 { // 计算 pos1 与 tmp_data 中满足下标条件的颗粒之间的引力 171171 //masses 颗粒质量表 172172 //tmp_data 环形传输数据 173173 //loc_forces 本进程颗粒受力数据 174174 //pos1 本进程颗粒位置数据 175175 //loc_n1 pos1 中颗粒数量 176176 //rk1 pos1 中 process owning particles 177177 //loc_n2 temp_data 中颗粒数量 178178 int gbl_part1, gbl_part2, loc_part1, loc_part2; //rk2 tmp_data 中 process owning contributed particles 179179 for (gbl_part1 = rk1, loc_part1 = 0; loc_part1 < loc_n1; loc_part1++, gbl_part1 += p) //n 总颗粒数 180180 { //p 参与计算的进程数 181181 for (gbl_part2 = First_index(gbl_part1, rk1, rk2, p), loc_part2 = Global_to_local(gbl_part2, rk2, p); loc_part2 < loc_n2; loc_part2++, gbl_part2 += p) 182182 { 183183 # ifdef DEBUG 184184 printf("Compute_proc_forces before, rank%2d> part%2d %10.3e %10.3e part%2d %10.3e %10.3e\n", 185185 my_rank, gbl_part1, loc_forces[loc_part1][X], loc_forces[loc_part1][Y], gbl_part2, tmp_data[loc_n2 + loc_part2][X], tmp_data[loc_n2 + loc_part2][Y]); 186186 # endif 187187 Compute_force_pair(masses[gbl_part1], masses[gbl_part2], pos1[loc_part1], tmp_data[loc_part2], loc_forces[loc_part1], tmp_data[loc_n2 + loc_part2]); 188188 # ifdef DEBUG 189189 printf("Compute_proc_forces before, rank%2d> part%2d %10.3e %10.3e part%2d %10.3e %10.3e\n", 190190 my_rank, gbl_part1, loc_forces[loc_part1][X], loc_forces[loc_part1][Y], gbl_part2, tmp_data[loc_n2 + loc_part2][X], tmp_data[loc_n2 + loc_part2][Y]); 191191 # endif 192192 } 193193 } 194194 } 195195 196196 void Compute_forces(double masses[], vect_t tmp_data[], vect_t loc_forces[], vect_t loc_pos[], int n, int loc_n)// 计算本线程各颗粒受到的引力 197197 { // masses 颗粒质量表 198198 const int src = (my_rank + 1) % comm_sz, dest = (my_rank - 1 + comm_sz) % comm_sz; // tmp_data 环形传输数据 199199 int i, other_proc, loc_part; // loc_forces 本进程颗粒受力数据 200200 // loc_pos 本进程颗粒位置数据 201201 memcpy(tmp_data, loc_pos, loc_n * sizeof(vect_t)); // 将本进程分到的位置数据放入 tmp_data // n 总颗粒数 202202 memset(tmp_data + loc_n, 0, loc_n * sizeof(vect_t));// 初始化 tmp_data 的引力数据 // loc_n 本进程分到的颗粒数 203203 memset(loc_forces, 0, loc_n * sizeof(vect_t)); // 初始化本进程引力数据 204204 205205 Compute_proc_forces(masses, tmp_data, loc_forces, loc_pos, loc_n, my_rank, loc_n, my_rank, n, comm_sz);// 计算本进程的颗粒间的引力作用 206206 for (i = 1; i < comm_sz; i++)// 计算本进程的颗粒与收到的新颗粒之间的引力作用 207207 { 208208 other_proc = (my_rank + i) % comm_sz;// 每次交换信息的对象不同 209209 MPI_Sendrecv_replace(tmp_data, 2 * loc_n, vect_mpi_t, dest, 0, src, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); 210210 Compute_proc_forces(masses, tmp_data, loc_forces, loc_pos, loc_n, my_rank, loc_n, other_proc, n, comm_sz); 211211 } 212212 MPI_Sendrecv_replace(tmp_data, 2 * loc_n, vect_mpi_t, dest, 0, src, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE); 213213 for (loc_part = 0; loc_part < loc_n; loc_part++) 214214 { 215215 loc_forces[loc_part][X] += tmp_data[loc_n + loc_part][X]; 216216 loc_forces[loc_part][Y] += tmp_data[loc_n + loc_part][Y]; 217217 } 218218 } 219219 220220 void Update_part(int loc_part, double masses[], vect_t loc_forces[], vect_t loc_pos[], vect_t loc_vel[], int n, int loc_n, double delta_t) 221221 { 222222 const int part = my_rank * loc_n + loc_part; 223223 const double fact = delta_t / masses[part]; 224224 # ifdef DEBUG// 输出在更新前和更新后的数据 225225 printf("Update_part before, part%2d %10.3e %10.3e %10.3e %10.3e %10.3e %10.3e\n", 226226 part, loc_pos[loc_part][X], loc_pos[loc_part][Y], loc_vel[loc_part][X], loc_vel[loc_part][Y], loc_forces[loc_part][X], loc_forces[loc_part][Y]); 227227 fflush(stdout); 228228 # endif 229229 loc_pos[loc_part][X] += delta_t * loc_vel[loc_part][X]; 230230 loc_pos[loc_part][Y] += delta_t * loc_vel[loc_part][Y]; 231231 loc_vel[loc_part][X] += fact * loc_forces[loc_part][X]; 232232 loc_vel[loc_part][Y] += fact * loc_forces[loc_part][Y]; 233233 # ifdef DEBUG 234234 printf("Update_part after, part%2d %10.3e %10.3e %10.3e %10.3e\n", 235235 part, loc_pos[loc_part][X], loc_pos[loc_part][Y], loc_vel[loc_part][X], loc_vel[loc_part][Y]); 236236 fflush(stdout); 237237 # endif 238238 } 239239 240240 241241 int main(int argc, char* argv[]) 242242 { 243243 int n, loc_n, loc_part; 244244 int n_steps, step; 245245 double delta_t; 246246 int output_freq; 247247 double *masses; 248248 vect_t* loc_pos; // *pos 变成了全局变量 249249 vect_t* tmp_data; // 用于环形交换的数据 250250 vect_t* loc_vel; 251251 vect_t* loc_forces; 252252 char g_i; 253253 double start, finish; 254254 255255 MPI_Init(&argc, &argv); 256256 MPI_Comm_size(MPI_COMM_WORLD, &comm_sz); 257257 MPI_Comm_rank(MPI_COMM_WORLD, &my_rank); 258258 MPI_Type_contiguous(DIM, MPI_DOUBLE, &vect_mpi_t); 259259 MPI_Type_commit(&vect_mpi_t); 260260 261261 Get_args(argc, argv, &n, &n_steps, &delta_t, &output_freq, &g_i); 262262 loc_n = n / comm_sz; 263263 Build_cyclic_mpi_type(loc_n); 264264 masses = (double*)malloc(n * sizeof(double)); 265265 tmp_data = (vect_t*)malloc(2 * loc_n * sizeof(vect_t)); // 前半段 tmp_pos 后半段 temp_force 266266 loc_pos = (vect_t*)malloc(loc_n * sizeof(vect_t)); 267267 loc_vel = (vect_t*)malloc(loc_n * sizeof(vect_t)); 268268 loc_forces = (vect_t*)malloc(loc_n * sizeof(vect_t)); 269269 if (my_rank == 0) 270270 { 271271 pos = (vect_t*)malloc(n * sizeof(vect_t)); 272272 vel = (vect_t*)malloc(n * sizeof(vect_t)); 273273 } 274274 if (g_i == 'g') 275275 Gen_init_cond(masses, loc_pos, loc_vel, n, loc_n); 276276 else 277277 Get_init_cond(masses, loc_pos, loc_vel, n, loc_n); 278278 279279 if(my_rank == 0) 280280 start = MPI_Wtime(); 281281 # ifdef OUTPUT 282282 Output_state(0.0, masses, loc_pos, loc_vel, n, loc_n); 283283 # endif 284284 for (step = 1; step <= n_steps; step++) 285285 { 286286 Compute_forces(masses, tmp_data, loc_forces, loc_pos, n, loc_n); 287287 for (loc_part = 0; loc_part < loc_n; Update_part(loc_part++, masses, loc_forces, loc_pos, loc_vel, n, loc_n, delta_t)); 288288 # ifdef OUTPUT 289289 if (step % output_freq == 0) 290290 Output_state(step * delta_t, masses, loc_pos, loc_vel, n, loc_n); 291291 # endif 292292 } 293293 294294 if (my_rank == 0) 295295 { 296296 finish = MPI_Wtime(); 297297 printf("Elapsed time = %e ms\n", (finish - start) * 1000); 298298 free(pos); 299299 free(vel); 300300 } 301301 MPI_Type_free(&vect_mpi_t); 302302 MPI_Type_free(&cyclic_mpi_t); 303303 free(masses); 304304 free(tmp_data); 305305 free(loc_forces); 306306 free(loc_pos); 307307 free(loc_vel); 308308 MPI_Finalize(); 309309 return 0; 310310 }
● 输出结果，与基本算法类似。8 进程 16 体，3 秒，时间步长 1 秒，舍去 debug 输出 1.476754e+00 ms；8 进程 1024 体，3600 秒，时间步长 1 秒，舍去 output 和 debug 输出 1.329172e+04 ms，在较大数据规模上稍有优势