实验三：稀疏矩阵-矩阵乘

负责助教：张齐颢 zqh23@mails.tsinghua.edu.cn

在本实验中，你将通过实现 GPU 加速的稀疏矩阵-矩阵乘法（SpMM）进一步熟悉 CUDA 编程以及 GPU 体系结构。

实验任务

计算公式

在 SpMM 中，其计算公式为 $C = A*B$，其中 $A$ 为稀疏矩阵，其形状为 $M * M$，其中包含非零元 $nnz$ 个；$B$，$C$为稠密矩阵，$B，C$ 的形状均为 $M * K$。可以用稠密矩阵乘法来理解，区别仅仅是 $A$ 中含有大量的零元。

存储方式

CSR格式：稀疏矩阵以CSR格式存储。以如下稀疏矩阵为例，可以用三个数组来存储它（均从0开始编址）。

PTR = [  0  2  4  7  8 ]
IDX = [  0  1  1  3  2  3  4  5 ]   
VAL = [ 10 20 30 40 50 60 70 80 ]

PTR 数组指出了哪些元素属于某一行。在上例中，PTR[0]=0，PTR[1]=2，其意思是在 IDX 和 VAL 数组中的第 [0,2) 个元素属于第 0 行。同样，PTR[1]=2，PTR[2]=4，其意思是在 IDX 和 VAL 数组中的第 [2,4) 个元素属于第 1 行。

IDX 数组指出了在该行上的具体位置，在上例中，我们通过读取 PTR 数组已经知道了 IDX 数组中的第 [0,2) 个元素属于第 0 行，通过读取 IDX[0]=0 IDX[1]=1，我们可以知道，第 0 行所拥有的两个元素在第 0 列和第 1 列。

VAL 数组指出了其具体数值，通过读取 VAL[0]=10 VAL[1]=20，可以知道，第 0 行的两个元素的数值是 10 和 20。

稠密矩阵以行主序方式存储。

项目介绍

.
|-- CMakeLists.txt
|-- cmake
|-- data
|-- include
|   |-- args.hxx
|   |-- data.h
|   |-- dbg.h
|   |-- spmm_base.h
|   |-- spmm_cusparse.h
|   |-- spmm_opt.h
|   |-- spmm_ref.h
|   |-- util.h
|   `-- valid.h
|-- script
|   `-- run_all.sh
|-- src
|   |-- data.cu
|   |-- spmm_cusparse.cu
|   |-- spmm_opt.cu
|   |-- spmm_ref.cu
|   |-- util.cu
|   `-- valid.cu
|-- test
|   |-- CMakeLists.txt
|   |-- main.cpp
|   `-- test_spmm.cu
`-- third_party
    `-- googletest

其中 spmm_base.h 是SpMM实现的基类，spmm_ref.* 是效率很低的 GPU 参考实现，spmm_cusparse.* 是利用 NVIDIA 的稀疏计算库的 GPU 实现，spmm_opt.* 是你需要实现的地方（请只修改 spmm_opt.h, spmm_opt.cu）。你需要实现的是 preprocess 和 run 函数。

在 test_spmm.cu 中，使用 Googletest 实现了多个测试项，其中有验证正确性的以及测试性能的。请在先过了正确性验证之后再运行性能测试。

实验步骤

spack load gcc@10.2.0
spack load cuda
spack load cmake@3.24.4
cp -R /home/course/hpc/assignments/2024/PA3/ ~
cd ~/PA3/
mkdir build
cd build
cmake ..
make -j4
# 运行单个数据点
srun -N 1 --gres=gpu:1 ./test/unit_tests --dataset <datasetname> --len 32 --datadir ~/PA3/data/ # 所有的数据集在 ~/PA3/data/ 中
# 运行全部 GPU 数据点
srun -N 1 --gres=gpu:1 ~/PA3/script/run_all.sh # 在 PA3/script 目录下
# 改变环境变量，仅仅运行单个测试，例如验证正确性（Validation）
GTEST_FILTER="SpMMTest.validation" # 运行全部 GPU 数据点
srun -N 1 --gres=gpu:1 ./test/unit_tests --dataset <datasetname> --len 32 --datadir ~/PA3/data/

其中 dataset 包含许多真实图数据（稀疏矩阵），稀疏矩阵的形状（$M$，在代码中是 kNumV）和非零元（$nnz$，在代码中是 kNumE）也被其决定。例如对于图数据 a，我们有两个文件，a.config 在同一行内存储了 $M$ 和 $nnz$，a.graph 第一行存储了 PTR 数组，第二行存储了 IDX 数组，VAL 数组在运行时随机初始化。

--len 决定了 $B$ 和 $C$ 的 $K$。数据在 PA3/data 中。可以自己造一个小的数据集来 debug，如下：

srun -N 1 --gres=gpu:1 ./test/unit_tests --dataset toy_graph --len 32 --datadir /path/to/your/data

在测试时会测试 len = 32, 256 的情况。

评分

正确性

正确性得分共占 $60\%$，包含两部分：基础分（占 $10 \%$）和加速比得分（占 $50 \%$）。 通过了 validation 的测试, 则得到基础分;

为了避免直接提交下发的参考实现得分，我们还设置了加速比得分。GPU 的加速比得分为需要超过 cusparse 的性能;

只要在 ≥20 个测试中超过 cusparse 的性能且结果正确, 就可以获得 全部 加速比分;

性能

性能得分共占 $30\%$, 针对 GPU 测试 13 个数据集 (script/run_all.sh中指定), 两种 $K$ 的长度;

• 对于每组测试用例，只有当你获得了正确性基础分后，才能得到性能分。每组测试用例的性能分数相同。
• 每组测试用例有一个性能线，超过性能线的同学将得到满分。
• 未达到性能线的同学，根据测试性能在 未达性能线同学 的排名给出每组测试用例的分数：每组测试用例各自排名，性能排名前 $10 \%$ 的同学得到 $100 \%$ 的分数，排名 $10 \%$ - $20 \%$ 的同学得到 $90 \%$ 的分数，依此类推。对于任何测试用例，获得正确性分数的同学将至少获得 $10 \%$ 的性能分数。

性能线

为避免对于每个数据集进行微调导致不必要的工作量，发布两种性能线：整体性能线和针对每个数据集的性能线，达到整体 或者 针对单个数据集的性能线，即可获得性能分。

• 如果你满足了整体性能线, 则得到当前setting下对应的所有数据点的满分(即使有数据集没达线)
• 如果你满足了某个数据集的性能线, 则得到这个数据集的所有性能分

整体性能线

以吞吐量作为性能指标，其计算方法是 avg(nnz/t)。

对于 K=32, 平均吞吐量达到 5e9 nnz/s

对于 K=256, 平均吞吐量达到 6.5e8 nnz/s

单个数据集性能线

以运行时间作为性能指标，单位为us

dataset	k=32	k=256
arxiv	350	2500
collab	620	4500
citation	8900	70000
ddi	250	1500
protein	8200	130000
ppa	8800	80000
reddit.dgl	17000	160000
products	32000	250000
youtube	2800	16000
amazon_cogdl	44000	400000
yelp	3400	27000
wikikg2	4400	25000
am	2300	13000

得分建议

根据往年的评测结果，最有效果且简单的优化是消除GPU的 warp divergence, 并让同warp的线程访问连续的存储.

实验报告

实验报告占 $10\%$。

实验提交

1. 实验代码：
2. 在截止日期之前将完成后的整个实验框架置于自己 home 目录下的 PA3 目录，如 /home/course/hpc/users/2020000000/PA3。
3. 实验报告：
4. 将 PDF 文件 提交至网络学堂。
5. 包含以下内容：
   1. 介绍你的实现方法，可以包括如何解决 ref 实现中的 warp divergence 的问题，如何利用各级缓存，如何改善数据的局部性，如何解决 load imbalance 的问题等；
   2. 展示不同优化对性能产生的影响，可以以单个数据集为例子详细分析；
   3. 在 $len = 32, 256$ 时的运行时间，及相对于 cuSparse 实现的加速比。

优化 Hint

• GPU 访存：在 SpMM 中，稀疏矩阵的一个元素代表了对于稠密矩阵的一行的访问，所以访存量很大，需要考虑到 GPU 的访存行为（coalesce memory access）来优化
• Warp divergence：稀疏矩阵的元素分布不规则，会导致 reference 实现中同一个 warp 内部的线程工作量差距很大，因为 warp divergence 导致一个 warp 的执行时间由最慢的线程决定。而在 reference 实现中，每个线程负责一整行稀疏矩阵相关计算；每行非零元相差大，所以有严重的 warp divergence 问题。可以通过改变并行维度，让同一个 warp 内线程处理相同的工作量；
• 局部性：稀疏结构带来不规则的数据访问，导致局部性很差，可以通过对图数据做预处理（在 preprocess 函数中处理，不占计算性能的运行时间；可以使用现有工具如 METIS，也可以自己实现），改变图数据结构，增加局部性；
• 负载不均衡：可以自己预处理图结构，来减少它的不规则性，优化计算过程中的 load imbalance 的问题。

GPU 优化技巧

可以通过 GPU profiling 工具 nvprof 对程序进行 profile：

# 得到程序运行的总结，包括整个程序运行过程中，各个 kernel 以及 CUDA API 的执行时间和次数
srun -N 1 --gres=gpu:1 nvprof ./test/unit_tests xxxxxxxxxx
# profile 单个 kernel 的执行情况，通过 --kernels 指定要 profile 的 kernel 的名称；通过 --metrics 指定要 profile 的是什么 metric，如 dram_read_bytes, achieved_occupancy 等，也可以指定为 all 来得到所有的 metric
srun -N 1 --gres=gpu:1 nvprof --kernels "KERNEL1|KERNEL2" --metrics "METRIC1|METRIC2" ./test/unit_tests xxxxxxxxxx

关于可以 profile 得到的性能指标以及 nvprof 更高级的使用方法可以参考 https://docs.nvidia.com/cuda/profiler-users-guide/index.html。

往届实现

图中显示了不同的实现在不同数据集上相对 cuSparse 的加速比（$K = 32$）：同学 A 通过改变线程映射的方式，基本消除了 warp divergence 的问题，其代码量相比于 naive 实现仅有小于 20 行不同；同学 B 除此之外，在 CPU 上预处理了稀疏矩阵，在一些不规则的数据集上，如 arxiv，取得了更多的提升。

注意事项

• 禁止任何欺骗评测程序的手段，包括但不限于直接输出时间、干扰校验程序运行、提前保存结果以在测试时直接输出等。一经发现，将取消本次实验得分。
• 你修改的文件应该仅限于 spmm_opt.h, spmm_opt.cu。即使修改了其他文件（如用于调试等目的），也要确保在 不进行这些修改 的情况下，程序能够正确编译运行。助教将替换所有其他文件为下发的版本后进行评测，以确保评分的正确性和公平性。
• 集群的登陆结点与计算结点配备了不同的 GPU，最终得分以计算结点为准（NVIDIA Tesla P100），程序在登陆结点的性能或正确性不作为评分依据。

严格查重，如果有任何参考借鉴，请在报告中写出来。

---

最后更新: 2024年5月15日

作者: Harry Chen (6.88%), Kezhao Huang (61.47%), Shizhi Tang (0.92%), 翟明书 (0.46%), zqh-wz (30.28%)