Lost's Blog

PKU HPCGame 2nd 2025 Writeup

 
  • Jan 27, 2025
  • 0 views

这是 PKU HPCGame 2025 的本人题解与思路

本人代码

官方题解与记录

前言

作为一个普通非信科人，HPC 不算太离题但也和本职专业差的不是一点半点的，所以几乎所有题目都是靠搜索来搞定的，剩余部分打几行注释写个开头 Copilot 就自动补全了，我的社工能力得到了极大锻炼，脚本小子能力得到了极大成长……呃不对，那个“雷方块”真的是我思考了一整天带一个晚上，又写了一整天，在完赛前一个小时成功赶出了第七个 bitset 优化的手动向量化版本的程序，成功拿到了仅仅一半的分数……哦，还有“着火的森林”必须自己写。

本题解大部分是在赛后写的，因为 HPC 充满了各种搜索、尝试、看不懂的迷惑时间，而且也不知道对不对，每题也不一定能满分。所以和 CTF 不太一样，除了前三题，就是签到、真正的社工搜索题“小北问答”和一道简单编译题，其它部分我没有选择一边做题一边写。

在本文中我不想贴太多代码，如果你想看就去上面的我的代码仓库里看看吧。

各题思路

A. 签到

把全部内容看完可一位一位得知答案 1898。

B. 小北问答

  1. 鸡兔同笼：某厂的 CPU 采用了大小核架构，大核有超线程，小核没有超线程。已知物理核心数为 12，逻辑核心数为 16，大核数量为 4，小核数量为 8。
  2. 编程语言：C 语言中，假设有函数 void f(const void **p);，我们有 void **q;，请问不使用强制类型转换，直接调用 f(q) 是否符合 C 的规范？否。
  3. CPU Architecture：ARM 架构的 sve2 指令集具有可变向量长度，且无需重新编译代码，就可以在不同向量长度的硬件上运行。sve2 指令集的最小硬件向量长度是 128，最大硬件向量长度是 2048。（百度搜索得到知乎文章）
  4. MISC：fp4 是一种新的数字格式，近期发布的许多硬件都增加了对 fp4 的支持。SE2M1（一位符号，两位 exponent，一位 mantissa）条件下，fp4 能精确表示的最大数字是 3，能精确表示的最小的正数是 0.5。（维基百科底下有一张表直接写出了所有可能）
  5. 储存：ZNS（Zoned Namespaces）SSD 是一种新型储存设备，关于传统 SSD 与 ZNS（Zoned Namespaces）SSD 的行为差异，以下哪些说法是正确的？（多选） A. 当写入一个已有数据的位置时，传统 SSD 会直接原地覆盖，而 ZNS SSD 必须先执行 Zone Reset 操作 B. 传统 SSD 的 FTL 会维护逻辑地址到物理地址的映射表，而 ZNS SSD 可以显著简化或消除这个映射过程 C. 当可用空间不足时，传统 SSD 会自动触发垃圾回收，而 ZNS SSD 需要主机端主动管理并执行显式擦除 D. 传统 SSD 一般支持任意位置的随机读取，而 ZNS SSD 只支持顺序读取 E. 传统 SSD 通常需要较大比例的预留空间 (Over-Provisioning)，而 ZNS SSD 可以将这部分空间暴露给用户使用 BCE（搜索后参考文章、文章和文章可知，D 选项应该是必须顺序写入，可以任意读取。） 第一次尝试错误，至于 A 选项，我使用了 GPT-4o，直接给了题干得到了回答：

    传统 SSD 并不会直接原地覆盖已有数据，因为 NAND 闪存的特性决定了写入数据之前必须擦除现有数据。传统 SSD 通过 FTL（Flash Translation Layer）来管理写入和擦除操作，通常会将新的数据写入空闲块，旧数据标记为无效，稍后由垃圾回收机制处理。

  6. OpenMPI：OpenMPI 是一个开源的消息传递接口 (MPI) 实现，在高性能计算领域被广泛使用。截至2025年1月18日，OpenMPI 发布的最新稳定版本为 5.0.6，在此版本的 OpenMPI 中内置使用的 PRRTE 的版本为 3.0.7。大家可以了解一下 PRRTE 的作用，OpenMPI 4 到 5 的架构变化，还挺有趣的。（官网下载页面和文档）
  7. RDMA：RDMA 是一种高性能网络通信技术，它允许计算机直接访问远程内存，从而大大降低了通信延迟和 CPU 开销。目前，主流的 RDMA 实现包括 InfiniBand、RoCE、RoCEv2 和 iWARP。下图中从左到右的四列展示了四种 RDMA 实现的架构图，请你说出按照从左到右的顺序，说出下图中的四列分别对应了什么 RDMA 的实现架构 DABC。（使用 Google 搜图搜到原图片就行，比如这篇文章后半段有这个图）
  8. HPCKit：HPCKit 是针对鲲鹏平台深度优化的HPC基础软件，请选择以下组件的具体作用。A. BiSheng B. HMPI C. KML D. KBLAS E. EXAGEAR，选项：1 高性能数学计算加速库、2 基础线性代数过程库、3 高性能通信库、4 X86到ARM的二进制指令动态翻译软件、5 编译器套件 53124（参考一下文档，然后根据常识猜一下缩写意思，最后用排除法）
  9. CXL：假设有以下条件：每次批处理需要传输的数据量为 1GB。GPU 每秒钟可以完成 10 次这样的批处理。传统架构下，CPU 到 GPU 的 PCIe 传输延迟为 50μs，传输带宽为 10GB/s。CXL 架构下，传输延迟降至 10μs，且数据访问可直接完成，无需显式传输。假设总训练任务包含 10000 次批处理。比较传统架构和 CXL 架构下完成任务所需的总时间，计算加速比（传统架构时间 / CXL架构时间），保留两位有效数字。\(\frac{50\times10^{-6} + 1/10 + 1/10}{10\times10^{-6} + 1/10} =\) 2.0
  10. 量子计算：初始状态为 \(\ket{0}\) 的量子比特，经过一次 Hadamard 门操作后，测量得到 \(\ket{0}\) 的概率是 0.5？经过两次 Hadamard 门操作后，测量得到的 \(\ket{0}\) 概率是 1？（这下专业对口了，用矩阵算一下就知道了）

C. 不简单的编译

依我之见这题应该只是想让我们选择正确的编译器、最快的编译选项，代码大概是不需要改的。

注意到 CPU 为 Intel Xeon 8358，那我会优先试试 Intel 的编译器，按照题给的文件搭好环境后，修改 CMakeLists.txt：

# Set the Intel compilers
set(CMAKE_C_COMPILER icx)
set(CMAKE_CXX_COMPILER icpx)

cmake_minimum_required(VERSION 3.20)

project(FilterProject LANGUAGES C CXX)

# SIMD & O3
set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -fast")
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fast")

add_executable(program main.cpp filter.F90)
target_link_libraries(program stdc++)
set_source_files_properties(filter.F90 PROPERTIES LANGUAGE C)

注意到我使用的编译选项是 -fast，这个嘛，直接在 Google 上搜到中科大超算中心文档，里面有介绍。其实手动试试，最大影响就是 -O3 优化和 -xHost 开启向量化。手动指定 AVX512 似乎不太快，不过这是我最初的想法，不然我也不会发现 -xHost 选项了。

评测机好像比集群上测试慢一些，最终获得 94.12 分，后来才知道应该再用 Fortran 改写就能满了，看起来是向量化不完全的样子。

D. 最长公共子序列

嘶，我错估了本题的难度，这 100 分的题其实还挺难的。最长公共子序列（LCS）是一个非常非常经典的问题，所以我从始至终的大想法都是去开源社区搜一下并行版本的代码，然后抄过来，可惜，我高估开源的水平了（

如果你搜一下就可以发现，并行 LCS 同样是 2023 年南京大学操作系统课程的一个作业题，但是那题不仅搜不到高性能的、有缓存优化的答案，也推荐的是 pthread 库而不是 openmp，所以最多提供一下其中一个思路。

第一个思路，也是我最终采用的思路，就是沿着反对角线去计算。这是一个非常棒的想法，LCS 的数据依赖是上面、左边和左上角三个格子，直接并行需要要求按着顺序计算，我试过，超慢。沿着反对角线就没有冲突了，可以简单地 #pragma omp parallel for 加在内循环上来并行。具体而言我拼尽全力进行了搜索后，确定思路参考了这个博客和这个问题，当两个字符串长度一样的时候，我可以抄这个代码，不过他没处理长度不一样的，所以我又融合了这个代码的串行版本，这边的串行版本就是沿着反对角线算的，直接并行就可以了。

这题有一个关键点就是，你需要尽量地省内存，而且需要按反对角线方式存数据来实现缓存优化，所以第一个思路符合这点的就只有上面两个项目的而已，可恶，没得抄。

顺带一提，我还有别的一些乱七八糟的想法，比如说：我一看这题就 100 分，又在第四题的位置，前面题正好是练习编译选项，后面题正好是 MPI，按照上一届的固有经验，这题用个 openmp 改个一两行就行了……

哈哈，我真对这题只有 100 分表示怀疑.jpg

第二个思路我提一下，在乱搜过程中我确实看到了在南大的 PLCS 题上有人使用分块的方法进行处理了，不过我想了一下这写起来太麻烦了，也没找到能直接抄的代码，分块的前几个周期会让别的线程空等，边界不好处理，所以我懒了……

嘛还有第三个思路，就是 LCS 可能有别的算法，诶，它确实有，在 DNA 分析领域因为碱基只有四五种，所以可以对更短的那个字符串进行一遍扫描化成一个新的矩阵，处理每种字符的位置之类的信息。具体可以看看这个项目，包括那篇论文我也看了一遍。本题我确实试了一下这个办法，毕竟看生成的代码字符只有 65536 种，但是并没有太快，还爆内存了.jpg

最终我使用第一个想法抄了两个开源项目代码得到了 73 分，在一堆 dalao 的满分中显得非常菜！后来群里讨论我才知道这玩意好像没有自动向量化，要手写才能快，分块当然更是比较正解的了。

E. 着火的森林

中规中矩的 MPI 练习题，看来是固定剧目了，不过我犯了好几个错导致浪费不少时间。

本题使用 64 个单核进程去做计算，也就是尽量并行，那“魔法事件阶段”肯定是不能并行的，而“状态转移阶段”可以把格子切成 64 份进行计算，问题是两两之间的边界上的数据怎么处理。

切块是简单的，横着切成 \((n/64) \times n\) 就行了。我一开始的想法是让一个主线程来统一数据，在每次进行“状态转移阶段”开始前把“森林”的状态发送给其它节点，算完后在回合结束时回收所有数据。诶，正好我搜索到了 MPI_Bcast 这个广播函数，那算一下内存，我觉得把全部数据同步到所有节点是非常可行的！然后我吭哧吭哧写出串行版本，“验证正确”后吭哧吭哧改成 MPI 花了几小时后，成功地超时了！

这题非常不好的一点就是评测超时的报错很奇怪，导致我迷惑了好久，当然我后来才发现我不仅超时也算错了。

我想了想，同步到所有节点传递的数据量可能有些大了，每个时间步迭代就广播数据非常非常慢，那还得是让所有进程去同时读数据、同时做计算，只好在“状态转移阶段”前把自己负责的“森林”的部分的上下两个边界分别发送给前后两个节点就行了。

核心代码就这一点，我这里让所有节点在一开始都保存了整个“森林”，当然后续非负责部分是用不到的，状态也是错的，这只是为了代码好改方便而已：

// 覆盖 all_forest
for (int x = x_start; x < x_end; x++) {
    for (int y = 0; y < n; y++) {
        all_forest[idx(x, y, n)] = new_forest[idx(x % my_n, y, n)];
    }
}

// 通信，交换边界
// 第一步，向后传递
if (rank < size - 1) {
    MPI_Send(new_forest.data() + (my_n - 1) * n, n, MPI_BYTE, rank + 1,
             0, MPI_COMM_WORLD);
}
if (rank > 0) {
    MPI_Recv(all_forest.data() + (rank * my_n - 1) * n, n, MPI_BYTE,
             rank - 1, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
}

// 第二步，向前传递
if (rank > 0) {
    MPI_Send(new_forest.data(), n, MPI_BYTE, rank - 1, 1,
             MPI_COMM_WORLD);
}
if (rank < size - 1) {
    MPI_Recv(all_forest.data() + (rank + 1) * my_n * n, n, MPI_BYTE,
             rank + 1, 1, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
}

MPI_Barrier(MPI_COMM_WORLD);

吭哧吭哧又改了一两小时，提交，诶，还就那个 Wrong Answer，和之前一模一样……嗯？我开始怀疑人生了。不过这次真的是 WA 了就是喽，只不过评测反馈完全看不出来和之前超时的区别。

一个巧合就是，这题给的两个样例数据都无法测出我写错了……我还是用没并行的版本去交了一下才意识到可能真的弄错了。

当时原话：我是傻逼，树被雷劈才着火，我把灰烬被雷劈也变成火了！草草草，查了我整整一天！！！

修改完轻松满分。

F. 雷方块

这题是我最花费时间，做得最累的，改的版本最多的一题。我十分怀疑这题没有给出 baseline 凭什么只有 200 分啊喂（当然换个思路幸好这题只有两百分，因为我没拿满，倒是很多 dalao 满了）！从本身算法，到优化复杂度，到最后的卡常数操作，步步惊心，全都是坑，真的恶心至极，体验差到爆炸。甚至为了这题在最后一天前面的半夜里躺着满脑子都在思考，完全睡不着，根据手环记录，当晚躺了八九小时，实际只睡了三四小时……

那么这题的记录我会写的详细一些，让人感受一下带痛苦面具的滋味。

问题的算法

这题本身就是一个改版的开关灯问题（Light Out Problem），域上的线性代数算法我是知道的，因为以前搓过小规模的玩具，但是这题的规模很大，是 $n_2 \times n_1$ 的矩形格子，内存考虑下就需要稀疏矩阵，还不让装一些第三方库……

搜一搜，就能发现很多关于类似问题的讨论，我的第一个想法是采用复杂度更低的别的算法，比如这个回答提到了一个非常精妙的算法，或者你可以在这个回答中看到一种更简单的说法，当然二者原理是一样的。这个算法的原理是仅仅第一行有缺陷的状态会通过一种朴素的操作传递到最后一行变成最后一行有缺陷的状态，而这朴素的操作就是对当前行变换以消除上一行的缺陷。也就是说，仅第一行的缺陷和仅最后一行的缺陷存在一种线性对应关系，能写成一个矩阵进行求解，再考虑对初始状态进行朴素操作留下最后一行的缺陷，此时就可以求解整个问题了，这样时间复杂度大概是 \(\order{\max (n_1^3, n_1^2 n_2)}\)，当然我没有具体实现一下所以也不清楚这对不对。

我觉得这个想法非常的 amazing，想了几乎一个晚上，直到我彻底否定这个算法对于这题的可能性。本题不行，因为它有空的块！这非常难处理，因为从上往下关灯的操作，遇上空的块就会出现各种奇怪的情况，理论上甚至有能卡死我的构造，更加地费脑子，所以我最终还是放弃了这个算法思路。

那我们还是回到最朴素的算法吧，就是把所有格子的 $(n_2 n_1) \times (n_2 n_1)$ 维邻接矩阵 $A$ 写出来，把每个格子上的操作和从初始到目标的所需操作看做是 $(n_2 n_1) \times 1$ 维数向量 \(\vb*{x}\) 和 \(\vb*{b}\)，解一个数域 \(\mathrm{GF}(3)\) 上的线性方程 \(A \vb*{x} = \vb*{b}\)。

线性方程算法

本题测试点最大 $n_2 = 2 n_1 = 1024$，如果你进行一个基本的数量级估计，那么你就会惊奇地发现这题的内存开销大到离谱，所以稀疏矩阵是必须的。

很好，问题变成的快速地求稀疏矩阵线性方程的办法，众所周知，LU 分解是优于朴素高斯消元法的，所以我的先去试着找找能简单抄写的开源库……诶，真有，一个叫做 SpaSM 库进入了我的视野，它是我找到的唯一一个在模 $p$ 上解稀疏矩阵方程的库，使用的确实是 LU 分解算法，符合要求，就是不太好抄，所以这题到这里被我搁置了几天。

嘶，做到最后几天，看了眼发现排名掉得有点多，感觉这题不做甚至就二等奖难保了，所以不得不捡起来接着做。

那怎么办呢，开源库直接抄过来！我大概花了五到六个小时去把我所需要的函数合在一个文件里，又改写了一点，取消了它去调用别的第三方库的密集矩阵解法的引用，终于在三千行左右搞定了一切，拿到了 10 分基本分……

对，没错，超时了，我试了下，大概正好比最低要求多一分钟左右。

仔细研究了一下可以优化的地方，我发现这个库会去试图提取主元，然后将剩下来的部分判断，如果变成了密集矩阵或者长宽差距极大的矩阵，就交给第三方库去处理。但是我这里显然没这个精力和时间把另外两个第三方库抄进来了，而且那玩意看起来更加的复杂，不是能容易解决的样子。

那，怎么办呢，真的只能高斯消元了吗……想法很多，我又搜了搜，这种对称稀疏矩阵有巧妙的办法去优化求解，比如这里提到了 Cholesky 分解，只能对正定对称矩阵操作什么的……但是我完全不想看了，又去做别的题了。

高斯消元法优化

这题无论矩阵稀疏与否，如果采用朴素的高斯消元法，那么时间复杂度应该是 \(\order{n_1^3 n_2^3}\)，对于测评数据来说肯定是会超时的。

我最后只剩这一题和那个完全看不懂的 J 题了，排名又掉了不少，看了眼榜，连续几个人分差非常小，思来想去，还是这题有点思路，再努努力吧。

这时候只剩最后 24 小时了，我实在是睡不着，在各种想法中反复确认，对着题目提示上的一句话“复杂度为 \(\order{n_1^3 n_2}\) 的朴素高斯消元算法确能拿到满分”想了一宿，在迷迷糊糊快睡着的时候，我突然悟到了……好像邻接矩阵上每个点和别的点有关联最多向下跨格子上的一行也就是 $n_1$ 个位置，也就是说向下消元只需要从当前主元在被换上来前的那个位置开始，再往下扫 $n_1$ 行就行了，这样复杂度就对上了！

结果这下更睡不着了，那晚上全在想这对不对，向上回代怎么办之类的……当然你知道的，没草稿纸全心算我只能想个大概，得一遍遍反复确认反复思考，当时就是觉得回代也就是最多往上 $2 n_1$ 行就行，不过确实想对了就是喽。

稀疏矩阵实现优化

顺带一提，在上面的想法正式确认可行前，我还手写了普通的高斯消元法，分别使用了 unordered_map 和普通的两个数组来存储每一行的非零元信息，当然这俩都挺慢，属于没分的范畴内。

此时我已经确认了使用普通的两个数组，一个记录非零元列号，一个记录非零元的值，这样是比那个常数极大的字典更快的。所以我起床后的第一个实现就是这样修改出来的，可惜，还是很慢，不过是接近有分的了。

又卡了，明明算法复杂度到位了，那这次还好，想了一两个小时就有思路了，感觉还是缓存优化的问题，用密集矩阵就可以连续访问和向量化了，但是稀疏的这种需要反复访问和重建数据，会很慢。那怎么办呢，还记得我说过最多跨 $2n_1$ 行非零吗，那想一想就知道，每一行从第一个非零元到最后一个非零元最多有 $2n_1 + 1$ 个元素，所以我们只要记录每一行第一个非零元的位置和接下来 $2n_1 + 1$ 个元素就行了！这样甚至很巧，对于高斯消元这种需要找第一个非零元位置的算法，有更多优化的空间了。

噼里啪啦火急火燎地写出了代码，又花了很久去 check，毕竟写错了是很正常的事情，这里算第一个非零元的偏移就很容易错。最终在 openmp 的加持下，于下午六点左右拿到了 78 分，还有三倍左右的优化空间才能满分。

此时如果我接着手动向量化说不定能发现对齐的问题，能拿到更多分，但是这里我的思路就又变了，我想的是正好差四倍，那位运算可以 2 bit 表示一个位置，正好也是四倍优化啊，这“绝对”是正解吧。

位运算优化与手动向量化

很好，还有不到六小时，都花了这么长时间了，我再写一个版本不过分吧~

完全没用过，所以现场学了下 bitset 的用法，然后我又走入了歧途，我直接用一整个 bitset 去存一行的信息。这样有个问题，当计算模 3 的加减法的时候，我需要去把奇数位和偶数位分离出来进行计算，这显然不太好自动向量化。

写完一测发现速度反而更慢了，还剩不到三小时，嘶，再来一个优化！我们可以把奇数位偶数位分开，做成两个长度 $2n_1 + 1$ 的 bitset，一个叫做非零判断位，一个叫做一二判断位。对于有零和没有零的加法我们分别处理：两个数中有零的话，直接异或就行了；两个数都不是零，那么新的非零判断位等于两个数的一二判断位异或后取反，新的一二判断位等于两个数的一二判断位取或后取反。这样加法的逻辑就搞定了，减法很简单，只要把减数取逆当加法即可，取逆也很方便，将非零判断位和一二判断位异或后当成新的一二判断位即可。

做完整套逻辑分析后，我又噼里啪啦整出来一版，测试一遍通过，不过速度上也只是刚好和之前数组版本的一样而已。那为啥呢，接着网上搜索一番发现，哎，bitset 没有直接暴露数据接口，向量化可能存在问题……

不过也不是没有办法，这个问题底下的回答给了我思路，我试了一下，借着 copilot 强大的自动写码能力，成功改了一个 SSE 的版本出来，稍微检查修改后，跑了一下发现对了，速度上终于达到了新高。然后没时间了，我只能止步于此，稍稍试了下发现 AVX512 确实是更快一点的，交了完事。

最终核心代码如下：

#define MAX_N 1025
typedef std::bitset<MAX_N> ARR;
inline void bit_add(ARR &a1, ARR &a2, const ARR &b1, const ARR &b2) {
    // ARR nonzero = a1 & b1;  // 两个都是非零
    // a1 = (nonzero & ~(a2 ^ b2)) | (~nonzero & (a1 ^ b1));
    // a2 = (nonzero & ~(a2 | b2)) | (~nonzero & (a2 ^ b2));

    char *a1_ptr = (char *)&a1;
    char *a2_ptr = (char *)&a2;
    const char *b1_ptr = (const char *)&b1;
    const char *b2_ptr = (const char *)&b2;

    __m512i *a1_ptr_avx = (__m512i *)a1_ptr;
    __m512i *a2_ptr_avx = (__m512i *)a2_ptr;
    const __m512i *b1_ptr_avx = (const __m512i *)b1_ptr;
    const __m512i *b2_ptr_avx = (const __m512i *)b2_ptr;

    #pragma unroll
    for (int i = 0; i < MAX_N / 512; i++) {
        __m512i a1_val = _mm512_loadu_si512(&a1_ptr_avx[i]);
        __m512i a2_val = _mm512_loadu_si512(&a2_ptr_avx[i]);
        __m512i b1_val = _mm512_loadu_si512(&b1_ptr_avx[i]);
        __m512i b2_val = _mm512_loadu_si512(&b2_ptr_avx[i]);

        __m512i nonzero = _mm512_and_si512(a1_val, b1_val);
        __m512i a2_b2 = _mm512_xor_si512(a2_val, b2_val);
        __m512i a1_result = _mm512_or_si512(_mm512_andnot_si512(a2_b2, nonzero), _mm512_andnot_si512(nonzero, _mm512_xor_si512(a1_val, b1_val)));
        __m512i a2_result = _mm512_or_si512(_mm512_andnot_si512(_mm512_or_si512(a2_val, b2_val), nonzero), _mm512_andnot_si512(nonzero, a2_b2));
        _mm512_storeu_si512(&a1_ptr_avx[i], a1_result);
        _mm512_storeu_si512(&a2_ptr_avx[i], a2_result);
    }
    // 最后一个 1025
    auto nnz = a1[MAX_N - 1] & b1[MAX_N - 1];
    a1[MAX_N - 1] = (nnz & ~(a2[MAX_N - 1] ^ b2[MAX_N - 1])) | (~nnz & (a1[MAX_N - 1] ^ b1[MAX_N - 1]));
    a2[MAX_N - 1] = (nnz & ~(a2[MAX_N - 1] | b2[MAX_N - 1])) | (~nnz & (a2[MAX_N - 1] ^ b2[MAX_N - 1]));

}

inline int vec_sub(ARR &a1, ARR &a2, const ARR &b1, const ARR &b2,
                   bool factor) {
    // 注意：密集 bits，a1 a2 会被改
    if (factor) {
        // a - b
        // 按照 b1 反转 b2，注意 b1 也是 inv_b1
        bit_add(a1, a2, b1, b1 ^ b2);
    } else {
        // a + b
        bit_add(a1, a2, b1, b2);
    }

    // 找到第一个非零元素
    int nnz_idx = a1._Find_first();

    // 重新赋值
    a1 >>= nnz_idx;
    a2 >>= nnz_idx;
    return nnz_idx;
}

void solve(SparseMatrixRow &A, std::vector<u8> &B, std::vector<u8> &X, int n1) {
    int n = A.n;
    std::vector<int> &nnz_idx = A.nnz_col;
    std::vector<ARR> &rows_first = A.first;
    std::vector<ARR> &rows_second = A.second;

    // 高斯消元 mod 3
    for (int i = 0; i < n; i++) {
        // 找到第 i 列的主元
        int j;
        int j_end = std::min(i + 1 + n1, n);
        for (j = i; j < j_end; j++) {
            // 直接看非零列号
            if (nnz_idx[j] != i) continue;
            if (i == j) break;

            // 交换行
            std::swap(nnz_idx[i], nnz_idx[j]);
            std::swap(rows_first[i], rows_first[j]);
            std::swap(rows_second[i], rows_second[j]);
            std::swap(B[i], B[j]);
            break;
        }

        ARR &pivot_row_first = rows_first[i];
        ARR &pivot_row_second = rows_second[i];
        u8 pivot = pivot_row_second[0] + 1;

        int k_end = std::min(j + 2 + n1, n);
        // 消元
        #pragma omp parallel for
        for (int k = j + 1; k < k_end; k++) {
            // 直接看非零列号
            if (nnz_idx[k] != i) continue;

            // 两个向量相减
            // 注意 nnz_idx 是需要加上偏移的
            bool factor = get_factor(rows_second[k][0] + 1, pivot);
            nnz_idx[k] += vec_sub(rows_first[k], rows_second[k],
                                  pivot_row_first, pivot_row_second, factor);
            B[k] = sub(B[k], B[i], factor);
        }
    }

    // 回代
    for (int i = n - 1; i > 0; i--) {
        u8 pivot = rows_second[i][0] + 1;
        int j_end = std::max(i - 2 * n1, 0);

        // #pragma omp parallel for schedule(static)
        for (int j = i - 1; j >= j_end; j--) {
            // 需要算一下偏移
            int col_i_idx = i - j;
            if (rows_first[j][col_i_idx] == 0) continue;

            u8 element = rows_second[j][col_i_idx] + 1;
            bool factor = get_factor(element, pivot);
            B[j] = sub(B[j], B[i], factor);
        }

    }

    // 解
    #pragma omp parallel for schedule(static)
    for (int i = 0; i < n; i++) {
        u8 pivot = rows_second[i][0] + 1;
        if (pivot == 1) {
            X[i] = B[i];
        } else {
            X[i] = (3 - B[i]) % 3;
        }
    }
}

你知道吗，我做了这么多优化，花了这么长时间，最后只拿到了仅仅一半的分数，94 分而已，而很多 dalao 都是满分，人与人的差距真是大啊……还不清楚具体少了啥，不过可以肯定的是，我手写向量化忘记对齐了，我一直以为速度是差不多的来着。

唉，打 HPC 害人不浅啊……

G. TopK

一开始没想到怎么做，试了下分块后 partialsortperm 的并行，拿不到多少分，所以真的按提示去写了一个基于二进制的基数筛，筛掉数字太小的再调函数，嘿，更慢了！

然后我先去做别的题了，直到群里看到有人问装第三方包的事情，诶，为什么要这样呢……感谢这位老哥哈，这下验证了我之前搜到了东西是有用的，就是 DataStructures 库的 Nlargest 函数，可以看看这个贴子甚至写出了更符合本题要求的实现。

直接替换之前分块并行的函数，随手 196 分，差一点满不了就算了，反正我不会 Julia。

当时原话：学了半天基数排序，看了一两天整数浮点数结构，搓了一个还没原生分块多线程跑得快……甚至远不如调库！！！

H. Posit GEMM

这我真的不会了，看了眼代码，copilot 自动改了个 cuda 的版本但是无法编译，看起来很费时间也没头绪，就放弃了。

I. 刀锐奶化

有点可惜，不是很了解浮点误差的相关知识，丢了起码一倍分数。不过除了浮点优化，分块近似的那部分我确实想不到了，我一直以为这题的正解应该是波动光学习题，利用远场衍射的一些性质去近似……

看了眼 baseline，外层循环可以直接交给 cuda，很容易配合着 copilot 写出了一份代码，跑出来大概耗时一分钟，相当快了。这时候仔细看看代码，sqrt 和 sin 以及 cos 函数肯定是最花时间的地方，那我去翻了翻文档，真的找到了一个 sincospi 函数同时计算 sin 和 cos，经测试，确实省了一些时间，大概到了 50 秒左右。

接着仔细研究了一下，计算距离时我们内部循环扫描的是镜子，那镜子到源的距离对于各个传感器是一致的，可以算好后缓存查表，这一步带来了一点点优化，时间来到了 40 秒左右。

最后核心代码就是这样的：

__global__ void compute(d3_t* __restrict__ mir, d3_t* __restrict__ sen, d_t * __restrict__ d_src_mirn_norm, d_t* __restrict__ data, int64_t mirn,
                        int64_t senn) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i >= senn) return;

    d3_t sen_i = sen[i];

    d_t a = 0;
    d_t b = 0;
    d_t tmp_sin = 0;
    d_t tmp_cos = 0;
    for (int64_t j = 0; j < mirn; j++) {
        // d_t l = norm(mir[j] - src) + norm(mir[j] - sen_i);
        d_t l = d_src_mirn_norm[j] + norm(mir[j] - sen_i);
        // d_t tmp = 4000 * l;
        // a += cospi(tmp);
        // b += sinpi(tmp);
        sincospi(4000 * l, &tmp_sin, &tmp_cos);
        a += tmp_cos;
        b += tmp_sin;
    }
    // data[i] = sqrt(a * a + b * b);
    data[i] = hypot(a, b);
}

上面的代码拿到了 50 分，理论上没有任何精度损失，没有做任何近似。另外无关紧要地提一下，hypot 函数并没有速度优势，但是 norm3d 会比直接写平方和开根号慢不少。

然后我就开始突发奇想，开始无限的“有点”错误尝试了，降低精度是很容易想到的，这题题面说了允许误差，所以降低到单精度浮点数是必须试一试的，然后试一试就炸了……好吧，此时我认为是这道题不允许我使用更低精度的浮点数进行计算。

接着就是想到经典的那个快速平方根算法，我抄了个双精度的版本，可惜啊，速度反而更慢了，而且精度也爆掉了，因为迭代次数完全不够。

这时候我就觉得这题是道物理题了，于是把输入数据的三维图样画出来，把输出结果的衍射图样也画出来看一看，诶，我发现镜子有四重旋转对称性，而传感器和光源的位置又非常巧，把光路整出了镜面对称性，也就是说样例输出的上半和下半部分的对称的，太棒了，这就可以省掉一半计算量。可惜，测评的样例貌似没有这个对称性，直接 wrong answer 炸了。

嘶，不对，那就看看距离是不是什么特别的数得了，诶，全在 602 左右，那难道是取某些特殊值附近就可以近似了吗……试了好几个小时，反正取各种值误差都很大，看起来是不太行的。

或者，镜子的坐标值非常密集，我取一半的量进行计算，最后把光强乘上二就行了。这其实是我觉得最棒的一个想法，实际上最后的图案确实没有什么差别，可惜还是爆了精度。

于是止步于此，我怎么也想不到，单精度浮点误差爆掉了就是因为距离全在 602 附近，乘上波长倒数后整数部分太大导致的。所以一个正确的浮点优化就是把整数部分去掉，再转单精度进行计算……甚至可以选用别的库的快速三角函数查表法……另外，正解应该还包括把传感器分块进行一个近似，因为传感器正好再一个平面上，所以距离成等差数列，可以简化……

J. HPL-MxP

这题就很容易看出来怎么做了，想法有了后就是查文档的事情，可惜，copilot 在行列谁优先上坑了我一下，环境变量也是个小问题，卡了一小会儿。

看完代码跑一遍发现耗时很长，看一下最重要的就是 BLAS 函数，自然而然想去调用更好用的 BLAS 或者 MKL 库。在“小北问答”那题中提到过，鲲鹏平台有个 KBLAS 库，太棒了，查一下文档找到迁移指南，主要优化这个 gemm 就行了。

不过实际操作起来还是有问题，无论编译器选什么，都会报错，说是找不到库，-lkblas 编译选项无效……嘶，检查一下，LD_LIBRARY_PATH 确实是有的，目录也包括了我想要的 KBLAS 库，也有动态库在里面，那为什么找不到呢……草，环境变量实在是太坑了，编译器其实读取的是 LIBRARY_PATH！所以脚本应该这么写：

#!/usr/bin/bash

module load bisheng/compiler4.1.0/bishengmodule
module load bisheng/kml2.5.0/kml
# module load gcc/compiler12.3.1/gccmodule
# module load gcc/kml2.5.0/kml

# echo $LD_LIBRARY_PATH
export LIBRARY_PATH=$LD_LIBRARY_PATH

make -j
./hpl-ai

KBLAS 看起来和 OpenBLAS 差不多，Copilot 直接帮我写了个大概，然后跑出来精度爆炸 nan 了，或者就直接段错误……嘶，看了眼代码，该不会是行列优先搞反了吧，换了一下后直接一遍跑过，看了眼精度对了算力够了，交上去直接满分。

K. 画板

明明是送分题却卡了一会儿，我简单说说吧。

提示建议使用 pthread 库进行多线程并行，我的第一个想法是直接开 10 个线程去分别计算，每个负责一个前缀。写完试了一下，没个半小时完全跑不出来，于是开始对代码进行分区计时，看看哪里比较慢。

首先查了一下，看了一眼，这个随机好像挺慢的，这种程序不需要随机数的性质有多好，所以能多快就多快。基于以前的经验，我把它换成了 MT19937 算法，确实快了那么一点，看起来影响并不大。

第二步，我看了下，每次校验都要转成字符串进行字串比较，这实在是太糟糕了，能用字节数组就应该直接用才对，所以我花了一点时间，把前缀的字符串转成了前缀的字节数组，再改了算法最终的输出，直接对比字节，这样应该会快点。不过这题小提升是看不太出来的，没测具体速度快了多少。

这样还是过不了，起码交上去无输出，嘶……后来我才知道是测评机器错了罢了，其实现在就满分了……

我接着去优化了算法，标记时间后发现，最耗费时间的就是算法生成 hash 值，但这部分我又不怎么会改。那换个思路，每个计算出来的值不仅仅和一个目标前缀比较，而是和全部的十个目标一起比较呢，具体而言也就是说，用 7 个线程作为生产者产生数据，用 1 个线程作为消费者对比全部的目标前缀。

不过测试时发现算完了答案程序就死锁了，呃……我不太想花脑子了，直接使用 exit(0) 强制退出进程，不优雅但绝对好用！最后，这优化使得程序大概两到三分钟就能算完，这实在是太酷了。

L. EDA

不得不说这题是我花的时间除了签到题以外最少的一道题目了，它实在是太简单了，我很高兴轻松拿到了这大题的分。

这么长的文件，看了眼，做了下本地的 profile，啥也没看出来，反正耗时最长的是树的合并。于是我想优化这里，把函数名扔进了 GitHub 看看有没有人做过这样的优化……诶，我搜到了源文件，诶……原来这题的代码和源项目不一样啊……

比如这份 flute 代码，对着看看改了什么，woc 怎么有人在频繁使用的函数里动态分配数组啊！

随手改成静态数组，MAXD 也改成了原来的 150，加了个最外层循环的 openmp 并行，测了一下结果正确也没报错。时间上，intel 平台差个一秒左右，arm 平台可以正好卡在满分线上。

不过这题测评环境爆炸得有点厉害，重测也不知道结果如何，为了保险我交了一堆上去，满分是大概率的事情。

后记

本人最终得分 1187.12 分，总排名第 5，校内排名第 4，作为第二次参加此比赛，这个成绩我相当满意了。

这次比赛本身还是存在不少问题的，第一天的测评机全炸和容器一个也开不起来就十分让人难受，导致我第一个做的正式题居然是 cuda 题，本地能测试那就是好啊。然后，测评机器的波动还是有的，特别是最后一天感觉最为明显，大家都在疯狂交题导致慢了不少，不过比起上一届比赛来说已经算好了不少，波动没有那么离谱了。另外，题目的质量还是平庸的，没有特别能让我眼前一亮的东西。要么题目非常简单搜一搜就出来了，要么有点思路但做了半天拼尽全力也无法满分，要么就是看都不想看的特难题，似乎从易到难的引导性变得更差了。当然不可否认，我还是学到了不少东西的，以这个角度来看还是很值的。嗯换个角度来说，比赛的评分评奖系统就有些不合适，像是一二等奖的分界线仅仅差了不到十分而已，这就是一个小优化，甚至是一个测评波动的分差，但是前后的奖项和奖金差的确实有点多。对于 HPC 这种算力优化比赛，目前的评分机制还不够完善，不足以看出大家的差距，更加公平有趣、减弱内卷的评分评奖是亟需的，比如甚至可以按分数来分奖金，或者按照每道题做的情况分别评奖分钱之类的。

对我而言，今年明显感觉到压力更大了，去年啥也不会随便抄抄做做就能二等奖，今年不认真一点就榜上无名了。庆幸的是比赛日期是放假日子，算是很有时间的，能陪着大伙一起卷一卷，但是身体和精力上的损耗相当之大，仅次于天天熬夜神志不清的 CTF 抢一血大赛。无论如何，希望这种稀有的比赛能够越办越好，我能够进步一点接着捞点小奖。

Search