================================
快速数无向图中的三角形个数的程序
================================


运行方式
========

依照 Rust 语言的传统和推荐做法，使用了 cargo 作为依赖管理工具和构建系统。程序使用了很多 Rust 的预览特性，请使用 nightly 工具链构建。

运行程序的机器的 CPU 必须支持 AVX512 指令集。

程序有 memmap2、rayon 和 rayon-k-way-merge 三个外部依赖。其中 memmmap2 是封装了跨平台的 mmap 调用的很好用的库，rayon 是 Rust 世界中最完善且通用的线程池，rayon-k-way-merge 是一个并行的多路归并算法，因为觉得它可能可以复用到更多程序中，所以我单独写了一个库。

程序运行时必须接收一个参数即输入文件路径。输入文件要求以 CRLF 为换行符，每行两个十进制表示的 32 位无符号整数，用制表符分隔，表示图中的一条边。两个整数分别为边的两点的编号，输入文件中允许存在重边。

编译程序只需要 `cargo build --release` 即可。编译出来的二进制程序会放在 `target/release/triangle-counting-lab` 这里。

一个运行程序的命令行示例：`./target/release/triangle-counting-lab /data/soc-LiveJournal1.txt`


运行结果
========

测试环境如下：

    CPU: 16-Core Intel Xeon Gold 5218 (-MCP-) speed: 2295 MHz
    Kernel: 5.15.0-107-generic x86_64 Up: 23h 18m Mem: 2518.9/15985.7 MiB (15.8%)
    Storage: 40.00 GiB (88.5% used) Procs: 504 Shell: bash 5.0.17 inxi: 3.0.38

程序可以在 1.28 秒内完成对 https://snap.stanford.edu/data/soc-LiveJournal1.html 数据集的读取和处理，能够正确计算出给定无向图中的三角形个数。数据集包含 4,847,571 个节点和 68,993,773 条边，输入大小约为 1.1 GB。

    test-KVM 15:36 (master) ~/tc/target/release
    0 ./triangle-counting-lab /data/soc-LiveJournal1.txt
    compute chunk size ratio: 8
    scan chunk size ratio: 2
    [ScopeTimer] 0.44753418 seconds | scanning and parsing
    [ScopeTimer] 0.13225858 seconds | merging thread results
    [ScopeTimer] 0.10247132 seconds | deduping edge list
    [ScopeTimer] 0.05148247 seconds | adjusting memory layout
    [ScopeTimer] 0.05831781 seconds | building adjacency list
    [ScopeTimer] 0.42185576 seconds | computation
    285730264
    [ScopeTimer] 1.27462609 seconds | totals


算法流程
========

定义：Lower(i) 为编号为 i 的节点的邻居中，编号小于 i 的节点所组成的集合。

1. 读取文件，获得边的列表。
2. 对边的列表排序去重，得到边的集合。
3. 根据边的集合建立邻接表。
4. 对每个节点 i，求 Lower(i)。Lower(i) 中的每个元素 j 对答案的贡献为 Lower(i) 与 Lower(j) 的交集的大小

这是一个 Vertex Ordering to List Triangles[1] 算法的简化版本。

[1]: https://github.com/lecfab/volt


为什么它能跑得比朴素实现快
==========================

~~当然是因为它是 Rust 写的！~~


读入文件的优化
--------------

* 并行读取文件，并且解析数据
* 使用了 mmap 系统调用来完成文件读写

首先明确：C 标准库的 fread 函数在我们的场景下并不适用。因为 fread 函数本质上是 read 系统调用套了个缓冲区，每次调用 fread 时如果缓冲区是空的，fread 函数就会调用 read 系统调用将缓冲区填满；如果缓冲区不空，fread 就直接将缓冲区中的数据作为结果返回。在少量多次顺序读写的场景下，这个策略可以通过减少系统调用次数的方式来降低所有读取操作的平均开销。

而我们操作数据的模式主要是大量的，次数很少的读写，每次可能读取几十几百 MB 的数据。这个时候，fread 函数所带的缓冲区没有任何意义，反而额外增加了内存拷贝的开销。数据要先从内核拷到 FILE 结构体内部的缓冲区，再拷到 fread 函数的参数所指向的内存，比直接使用 read 调用读反而更加耗时。

一般来说，Linux 系统下的文件 IO 都是有内核维护的缓冲机制的，这个缓冲机制主要就是 page cache。page cache 是内核中缓存磁盘文件的数据结构，单位为内存页（通常大小为 4 KB）。每个内存页与磁盘文件上的一段数据直接关联——即读写磁盘时会先读写 page cache，等因为被修改而与磁盘数据不一致的页面（即 “脏页”）占比达到一定程度时，再统一写回到磁盘。

我们最终选用的 IO 方式为 mmap。mmap 系统调用主要作用于页表，它会指定一块内存区域代表被映射的文件，读写这块内存时通过缺页中断来真正完成对文件的访问。

与 mmap 竞争的读入方式为基于 read 系统调用的读写方式，两者都和 page cache 紧密关联。

read 系统调用的读入流程如下：

1. 发起系统调用
2. （如果数据没在 page cache 中）从磁盘将数据读入 page cache 中
3. 将 page cache 中的数据拷贝到调用者提供的缓冲区中
4. 系统调用完成

使用 mmap 方法时读入流程如下：

1. 程序试图访问 mmap 所映射的内存
2. （如果之前没有访问过这个内存页）触发缺页中断
    1. （如果 page cache 中没有被访问的内存页所关联的磁盘文件的内容）从磁盘将数据读入 page cache 中
    2. 设置页表，将内存页映射到用户空间
    3. 缺页中断返回
3. 读取完成

可以发现，read 调用的优点是执行路径简单、可以在单次调用中读取大量数据，缺点是存在额外的数据拷贝；mmap 调用的优点是没有额外的数据拷贝、编程模型简单（访存即是读取文件），缺点是每次访问新页时都需要陷入缺页中断。对比两者优缺点之后，我们可以很粗略地认为，如果拷贝一个内存页的代价大于一次缺页中断的代价，那么 read 调用占优；否则，mmap 调用点优。

对于 read 来说，内核会通过 copy_to_user 函数来将数据拷贝到用户空间，copy_to_user 最后调用的还是 memcpy，内核的 memcpy 的实现是 4 个字节为一组拷贝的；对于 mmap 来说，我只找到了一个十年前的讨论[1]，里面提到了一次 page fault 需要 1000 cycles 的开销。

[1]: https://plus.google.com/app/basic/stream/z12fybdg3uqeh5nhu04cj1yx1py4xrsb1ps0k

懒得管了，脑测了一下觉得两者耗时应该会是差不多的，所以决定用写起来更方便的 mmap 做法来写。


读入文件的可能的进一步优化
--------------------------

在我们对 mmap 和 read 两种方案做评估时，我们认为在 page cache 未命中时，两者读盘的开销是相同的。实际上，虽然读盘开销不可变更，但是我们可以将磁盘读写放到后台进行，让数据处理（分行、将字符串解析成整数）和磁盘读写交叠进行，从而掩盖掉部分磁盘读写的开销。

以 read 调用为例，记磁盘 IO 的部分为 L，数据处理的部分为 A，那么现有的先读取再处理的方案是这么做：

    LLLLLLLLLLLLAAAA
    时间 --->

如果我们交替读取和数据处理，那么它执行时应该会是这样的结果：

    LLLALLLALLLALLLA
    时间 --->

而 read 调用实际上是有预读机制的，简单说就是在 read 调用会在读到所要求的字节数量后立刻返回，但是还是会向前多读几个页面，这几个页面便是预读页面，保存在 page cache 中。所以实际上，我们程序执行应该是这样的结果：

    LLLALLALLALLA      <-- 用户态的程序
       L  L  L         <-- 后台预读
    时间 --->

可以看出，这里数据处理和磁盘读写形成了交叠的结构，从而掩盖了部分磁盘读取的开销，让整个过程加速了一些。

在 Linux 下，可以用 readahead 系统调用控制 read 的预读行为；可以用 readaround 系统调用控制 mmap 的预读行为。


并行读取任务划分设计
--------------------

最简单的并行读取便是将文件大小分割成 N 块，每一块由一个线程来读来处理。遗憾地是这种方法在我们的问题里是行不通的，因为如果分块边界没有落在换行符上，就会让同一行被两个线程分别处理，产生两个错误的结果。

解决方案也很简单：首先根据文件大小平分成 N 块的做法求出大致的块边界，接着将每个块边界向右移动直到遇到换行符为止。这样便是一个合理合法的分块方案。

关于 N 的大小如何选取的问题：假设 CPU 的数量是 ncpu，一个简单的想法是取 N = ncpu。然而这相当坏，因为并不能保证数据大小相近的块计算量也相近。一般我们会采用的做法是取 N = M * ncpu，总共划分出 M * ncpu 个任务。负载均衡通过线程池的任务窃取机制来实现。

M 的要求是，既不能大到增加数据合并和任务调度的开销，又不能小到保证负载均衡的效果。这里我们按照经验取了 M = 2。实际最佳值是多少需要消耗大量核时依靠枚举或者优选法试出来。



SIMD 加速的分行算法
-------------------

考虑一行最多可能出现多少个字符：32 位无符号整数不会超过 10 个字符，中间用制表符隔开，行末是 \r\n 换行符，也就是说单行不会超过 10 + 1 + 10 + 2 = 23 个字符。

avx512 指令集能够同时处理 512 位的向量，也就是 64 个 8 位整数，即 64 个字符。

64 > 23，所以，我们有这个重要结论：将行首及其后 64 个字节读入一个 512 位向量寄存器后，要么行尾也在同一个寄存器里，要么这是一个无关紧要的注释行。

这个结论提示我们分割出一行甚至不用读两次，可以大大减少我们需要考虑的边界情况数量，也彻底消除了分行过程的循环。

分行算法的过程如下：

1. 将行首及其后 64 个字节读入 u8x64 的向量 data 中
2. 将 data 和 '\n' 比较，得到一个 64 位的结果 result，如果 data 某个位置上的值是 '\n'，则结果对应的位会是 1；否则为 0
3. 数 result 的前导零数量，将其作为本次处理的行的长度


SIMD 加速的字符串转整数算法
---------------------------

考虑将 512 位向量看作 u32x16 的向量。32 位无符号整数不会超过 10 个字符，10 < 16。

也就是说，如果我们让 32 位无符号整数的十进制字符串中的每一个字符占用 4 byte，这个字符串仍然可以放进一个 512 位向量中。这启发我们直接将十进制字符串中的每个字符当作 32 位整数，放进 512 位向量中后乘上对应的系数（指 1、10、100、1000 等等）后，求归约和即可得到结果。

算法过程如下：

1. 将字符串读入 u8x16 的向量中
2. 将向量中每个 u8 减去 '0'
3. 调整向量中字符串所在的位置，由头对头调整为尾对尾，方便后续乘上系数
4. 将 u8x16 转换成 u32x16
5. 将 u32x16 的向量乘上 [1, 10, 100, 1000, 10000 ... ] 的系数
6. 求 u32x16 的归约和作为结果

相比逐位乘系数求和的朴素算法，SIMD 加速的方法会比原方法快约 1/3。


并行多路归并算法
----------------

每个磁盘读取任务完成后，会对自己读取到的所有数据排序。在所有磁盘读取任务完成后，只要对所有结果做多路归并就可以得到一个完整的排好序的边表。

为什么需要对所有的边排序呢？

1. 我们需要处理的是无向图，而输入文件给的是有向图。需要将边去重以避免出现错误的结果
2. 排序后可以让源点相同的边排在相邻的位置，在后面建立邻接表时可以获得更好的 cache 亲和性
3. 我们需要知道节点编号的最大值，提前分配内存以减少数据扩容的开销。排序后正好可以从数组末尾获取这一信息

普通的串行多路归并算法是一个迭代的过程，它会给每一路的最小值建一个小根堆，每次从堆顶取出下一个值。这个算法因为用到了共用的数据结构（堆），并且所有求出的最小值需要按照顺序写入内存，导致它根本不适合做并行化。

如果要获得很高的并行度，就势必要抛弃共享的数据结构，转而思考数据并行的可能性。所以，我们需要考察原问题是否存在一些特质，使它可以被划分成若干个之间不需要相互依赖的子任务，且子任务与原问题性质相似但是规模较小。

考虑到我们是在做和有序的序列相关的问题，按照经验，我们可以先取一个结果序列的中间值 pivot 来看看它会不会对我们的问题有所帮助。观察选取了 pivot 后的结果序列，发现 pivot 左侧的所有元素均小于 pivot，而右侧均大于等于 pivot。如果能够把原始的待归并的 K 个数组，每个数组根据 “是否小于 pivot” 拆成两个部分，那么拆分之后获得的 2K 个数组要么所有元素大于 pivot，要么所有元素小于 pivot。我们可以发现，pivot 左右的内存操作互不干扰，这是一个适合并行的结构。可以让所有元素小于 pivot 的数组归并后写入结果数组中 pivot 左侧的部分，让所有元素大于等于 pivot 的数组归并后写入结果数组中 pivot 右侧的部分，这样两个子任务的结构也与我们正在处理的问题相似了。

这是个递归的过程，递归的终止条件很好设置：K 个数组中只有 1 个长度不为零，此时直接使用内存拷贝将数组内容拷贝到结果数组中即可；或者 K 个数组中的元素数量的总和小于某个阈值，继续划分任务的话并行度增加获得的收益无法抵消划分过程及任务调度产生的损耗。

最后是 pivot 选取的过程。最佳的方案肯定是每次都准确地选中结果序列的中点，这样可以稳定地使每次递归之后任务量变成原来的一半。这是可以做到的，只要二分 pivot 的值即可。设所有数组中的元素数量之和为 N，数组中所有数的最大值为 M，则二分的次数为 log(M)，二分过程中每次判定需要用二分过程的中点值 mid 在所有数组中二分查找，求出所有 K 个数组中小于 mid 的元素的数量，复杂度为 K log(N)，总复杂度为 K log(M) log(N)。对于划分这种简单的任务来说，K log(M) log(N) 的开销还是有些不能接受。

我用的 pivot 选择方案为选择待归并的 K 个数组中最大的那个数组的中点。复杂度为 O(K)，至少能够保证每次划分过后最大的那个数组将会变成均匀的两半。

这个算法的性能相当不错，测试的时候并行的 K 路归并算法仅比 K 个并行的 memcpy 操作慢了 40% 。细心做一些剪枝和边界条件的优化，有望使这个算法吃满内存带宽，达到 memcpy 同等水平。


并行去重
--------

相比归并，去重操作的并行化实现起来非常简单。只需要将数据分块，块内去重，对块边缘的元素重复情况做一些处理后合并就可以了。没有什么值得一提的优化点。


边表从 AOS 到 SOA 的转换
------------------------

简而言之，就是把这样的东西：

    struct Edge {
        unsigned src, dst;
    } edges[100];

转换成这样的东西：

    unsigned src[100];
    unsigned dst[100];

前者的排布就是 AOS（Array Of Struct）；后者的排布就是 SOA（Struct Of Array）。之前由于需要对边做排序和去重，所以我们采用了 AOS 的排布。之后因为关注点由边变成了各个节点的邻居，源点不再重要了，为了提高 cache 的亲和性以及让 SIMD 指令一次能处理更多元素，我们将其转换成 SOA 的模式，使一个源点的所有邻居在内存上连续排布。

这个过程是很好并行的，也很适合 SIMD。可惜我太懒了，看到它并行后跑得挺快的就没管它了。


并行建立邻接表
--------------

这个问题抽象出来就是：给出长为 N 的数组 src[]，对每个 i，求出 lbound[i] 和 ubound[i]，使 src[] 中下标为 lbound[i] 到 ubound[i] 的所有数都等于 i。很显而易见的二分查找问题。

这个显然是可以并行的，算得上是并行计算 hello world 级别的问题，如果要做数据并行的话，也很适合 SIMD。相信受过并行算法训练的同学心里已经有 100 种方法来把它用多核实现，这儿就懒得讲了。

如果你真不会并行二分：<https://github.com/rarocat/recruitment-2024-spring/>


SIMD 加速的 bitset
------------------

回顾我们的算法：

    对每个节点 i，求 Lower(i)。Lower(i) 中的每个元素 j 对答案的贡献为 Lower(i) 与 Lower(j) 的交集的大小

按照我们算法的流程，我们求出 Lower(i) 后，会用它和多个 Lower(j) 相交，并对交集的大小求和。我一开始是用 HashSet 来存 Lower(i) 和 Lower(j) 的，但是很快我就发现这超级慢。

去偷看了一眼 volt[1] 的代码后发现，他的做法是用 bitset 表示 Lower(i)，接着用每个 Lower(j) 的元素来访问 bitset，最后对结果求和。这下连交集也不用求了，整个算法变成了很纯的访存数数过程。

[1]: https://github.com/lecfab/volt

bitset 可以用 SIMD 加速吗？一般来说是没有这个必要的，但是我们的场景有其特殊性：Lower(j) 是一些下标，这些下标在内存中的位置是连续的。我们需要用这些下标寻址 bitset，对结果求和。

内存连续——可以读入向量寄存器里；对所有下标做相同操作——可以用向量操作；对结果求和——最后需要一个归约过程。好了，我们发现用 SIMD 来做批量的下标访问其实是有利可图的！

bitset 的 SIMD 向量化实现相当简单，就是把普通的寻址工作用 SIMD 写了一遍而已。但是有个值得考虑的问题：我们能用 SIMD 批量访问 bitset 中的元素，我们能批量设置或者翻转它们吗？

答案是不行。大部分 CPU 是不支持直接对 bit 寻址的。如果我们需要修改某个 bit 的值，我们至少需要把它所在的那个字节读出来，修改 bit 的值，再写回去。如果我们用 SIMD 来批量操作，那么会遇到一种情况：两个需要修改的 bit 刚好在同一个字节中，那么后写回的结果将会掩盖先写回的结果，我们会发现我们丢失了一次写入或者翻转操作。

事先扫描所有下标，通过重新排序来避免写入冲突是可以保证算法的正确性的，但是既然都在扫描下标了，那么在扫描的时候直接访问 bitset，把结果算出来肯定是一个更快的方法。

也就是说 SIMD 加速的批量 bitset 翻转要么有正确性问题，要么打不过串行访问。总之无利可图。


Rust 是好的
===========


零开销的迭代器
--------------

Rust 的迭代器实现是惰性求值的，这就保证了一般迭代器连接时不会需要申请新的内存，被过滤掉的元素也不会参与后续计算。

迭代器写起来很好玩。

如果用迭代器来对标循环的话，那么 rayon 提供的并行迭代器就是类似 C 中 OpenMP 一类的东西。它可以在一定程度上完成对数组的数据并行操作。同时，迭代器提供了 map filter 等方法用于转换和过滤元素，也提供了 collect 和 sum 等方法用于归约，这让我们在使用迭代器时能够更自然地以比较函数式的方式思考问题，从而自然地写出类似 map-reduce 的适合并行的结构。


好用的数组切片
--------------

Rust 里面数组切片能够像只读数组一样使用，切片支持迭代器访问，还能再切出小小切片，非常适合用来做数据划分。同时切片的标准库支持也很丰富，让人感觉到了家的温暖。


无畏并发
--------

首先不管 Rust 的 Send 和 Sync 标记机制，单说一个变量不会有两个可变引用就已经能够避免相当大的一部分数据竞争问题了。堆了这么多屎山竟然连一次并发问题都没有遇到过。

并发是这样的，编译器只要把好像有问题的代码拒掉就可以了，工程师想战斗爽要考虑的问题就多了。


Unsafe
------

Rust 提供了很好用的标准库，标准库里有很多有益的假设和保证：比如，保证所有 String 类型都是合法的 UTF-8 串，保证所有的内存都被妥善初始化，运行时执行下标检查等等。

这些检查能避免很多问题。然而我们现在在做的是即使方向盘掉了也得等跑完再捡起来修的 “我不管总之结果对了就是对了” 的超级高性能优化项目，有时这些标准库的检查甚至会成为性能热点。好在 Rust 提供了 unsafe 机制，标准库也有一些低开销但是标记了 unsafe 的函数，让我们可以直接油门踩死——至于跑着跑着会不会突然爆掉那就是人的水平问题了。


结论
----

Rust 确实比 C++ 好写。