7.6 举办了 PAC2023 的开幕赛，以轻松娱乐为主，我还没有体验过在短时间内完成优化工作（不过今天的工作量倒也不大）。这场娱乐赛给了我别样的体验，也让我发现了我有好的好多不足，以及一些以前没有注意到的坏习惯。
赛题介绍
赛题是针对一个朴素实现的康威生命游戏的优化。由于代码写得很朴素，所以比较轻松地就可以达到比较好看的加速比。
代码的核心部分如下：
for (int iter = 0; iter < max_iter; ++iter) {
printf("Iter %d...\n", iter);
for (int i = 1; i < Ny - 1; ++i) {
for (int j = 1; j < Nx - 1; ++j) {
int cnt = 0;
for (int k = 0; k < 8; ++k) {
int ii = i + dy[k];
int jj = j + dx[k];
if (a[ii][jj] == 1) {
cnt++;
}
}
if (a[i][j] == 1) {
if (cnt == 2 || cnt == 3) {
tmp[i][j] = 1;
} else {
tmp[i][j] = 0;
}
} else { // a[i][j] == 0
if (cnt == 3) {
tmp[i][j] = 1;
} else {
tmp[i][j] = 0;
}
}
}
}
for (int i = 0; i < Ny; ++i) {
for (int j = 0; j < Nx; ++j) {
a[i][j] = tmp[i][j];
}
}
代码分析
可以看到生命游戏一共要进行 1000 次迭代，每次都遍历一遍整个地图，对每个点都计算一下邻居的个数，再根据生命游戏的 B3/S23 规则决定下一次迭代时的状态。
B3/S23 规则是对生命游戏规则的概括缩写
这个朴素写法非常的暴力，不仅是直接访问每个点周围的八个点，而且还使用了 tmp 数组暂存元素，最后将 tmp 复制到地图数组 a 中，开销不可谓不大。
在本地的 64 核心集群上运行 baseline ，运行时间为 741335.164000 ms
在北京超级云计算中心上运行 baseline ，运行时间为 634306.056000 ms
优化策略
算法微调 && 访存优化
首先我想到的是对算法进行微调，将之前的对每个点求邻居个数改为，将每行与上下两行进行叠加，这样每个元素就是三个点和。这就把对一个「九宫格」的计算改成了两次对三个元素的加和。这样并不能减少计算量，但是能减少访存上的开销。
这样做的目的还有简化代码逻辑，方便编译器生成 SIMD 指令
在一番取舍之后没有使用前缀和算法来，这是因为担心在改变数据类型后，数据类型不够宽，无法存放前缀和。
其实这是后话，后面把
int改成了uint8_t
核心代码改成：
for (int iter = 0; iter < max_iter; ++iter) {
printf("Iter %d...\n", iter);
memcpy(sum, a, Nx * Ny * sizeof(int));
for (int i = 1; i < Ny - 1; ++i) {
for (int j = 0; j < Nx - 1; ++j) {
sum[i][j] += a[i + 1][j];
}
}
for (int i = 1; i < Ny - 1; ++i) {
for (int j = 0; j < Nx; ++j) {
sum[i][j] += a[i - 1][j];
}
}
for (int i = 1; i < Ny - 1; ++i) {
for (int j = 1; j < Nx - 1; ++j) {
int neighbors = sum[i][j] + sum[i][j - 1] + sum[i][j + 1] - a[i][j];
if (neighbors == 3) {
tmp[i][j] = 1;
} else if (neighbors < 2 || neighbors > 3) {
tmp[i][j] = 0;
} else {
tmp[i][j] = a[i][j];
}
}
}
memcpy(a, tmp, Nx * Ny * sizeof(int));
在本地的 64 核心集群上运行时间为 293731.276000 ms
在北京超级云计算中心上运行时间为 260813.238000 ms ，加速比约为 2.4x 。
openmp && 合并循环
之后就是最喜欢的「openmp，启动！」环节了。openmp 相对其他的并行框架来说非常好用，无脑推荐。
在敲代码时又突然发现叠加的两个循环其实可以合并，这下傻逼了。
所以核心代码变成了：
printf("Iter %d...\n", iter);
memcpy(sum, a, Nx * Ny * sizeof(int));
#pragma omp parallel for
for (int i = 1; i < Ny - 1; ++i) {
for (int j = 0; j < Nx - 1; ++j) {
sum[i][j] += a[i + 1][j];
sum[i][j] += a[i - 1][j];
}
}
#pragma omp parallel for
...
在本地的 64 核心集群上运行时间为 18918.941000 ms
在北京超级云计算中心上运行时间为 18864.766000 ms ，加速比约为 33x 。
O3 && SIMD
之前使用的编译参数都是 -O0 现在改成 -O3 ，大量的整数运算操作可以向量化。
在本地的 64 核心集群上运行时间为 10292.194000 ms
在北京超级云计算中心上运行时间为 15419.145000 ms ，加速比约为 41x 。
这里就可以看到本地集群和云计算中心集群机器性能的差别，本地集群单核性能弱于计算中心（体现在 baseline 本地跑了更久），但是计算中心的核心少于本地，因此使用 openmp 的情况下本地加速更明显。
改变数据类型
由于在九宫格中的计算结果在 0-9 之间，因此不需要很宽的数据类型来存储，而且也不需要考虑负数，因此我使用了 uint8_t 类型来替换 int ，这样操作之后加速比又提高了不少。
在本地的 64 核心集群上运行时间为 3004.615000 ms
在北京超级云计算中心上运行时间为 3858.775000 ms ，加速比约为 164x 。
去掉不必要的拷贝
最后再看看自己的代码，发现很多细节还可以处理。例如数据的拷贝，为了计算与相邻两行的和，我新增了一个数组，由于这个数组可以存储中间过程，我可以去掉 tmp 数组，直接对 a 数组进行操作，这样减少了 memcpy 的开销。
小插曲：baseline 对边缘的处理非常差，导致第一次迭代(iter=0)时需要将
a数组的边缘全部置为 0
还有就是在一开始将 a 数组拷贝到 sum 数组，这个其实也可以带到第一次的循环中进行，这对访存没有什么副作用，而且可以并行执行。
最后的代码：
for (int iter = 0; iter < max_iter; ++iter) {
#pragma omp parallel for proc_bind(spread)
for (int i = 0; i < Ny - 1; ++i) {
for (int j = 0; j < Nx - 1; ++j) {
sum[i][j] = 0;
sum[i][j] += a[i + 1][j];
sum[i][j] += a[i][j];
sum[i][j] += a[i - 1][j];
}
}
#pragma omp parallel for proc_bind(spread)
for (int i = 1; i < Ny - 1; ++i) {
for (int j = 1; j < Nx - 1; ++j) {
int neighbors = sum[i][j] + sum[i][j - 1] + sum[i][j + 1] - a[i][j];
if (neighbors == 3) {
a[i][j] = 1;
} else if (neighbors < 2 || neighbors > 3) {
a[i][j] = 0;
}
}
}
if (iter == 0) {
memset(a[0], 0, sizeof(a[0]));
memset(a[Ny - 1], 0, sizeof(a[Ny - 1]));
for (int i = 0; i < Ny; ++i) {
a[i][Nx - 1] = 0;
a[i][0] = 0;
}
}
在本地的 64 核心集群上运行时间为 247.119000 ms
本地由于单核弱而核数多，跑到 3000x 没什么鸟用，不具备参考价值😭
在北京超级云计算中心上运行时间为 318.769000 ms ，加速比约为 1990x 。
总结
之前打高性能比赛都是组队参加，总是需要大佬的指导和安排任务，这次个人赛暴露了很多不足，完全没有发挥出机器的全部性能。
不过往好处想，整体的优化思路是短时间内就想出来了的，都怪愚蠢的 baseline ，我花了很多时间在解决正确性的问题。
这次比赛没有机会尝试使用 MPI ，有时间的话用 MPI 继续做做这道题。