稀疏矩阵向量乘法（sparse matrix-vector multiplication，SpMV）^[1]是数值计算中的核心运算之一，在众多科学与工程领域中发挥着关键作用. 它不仅是求解线性方程组和优化问题的基础，而且在多个高性能计算（high performance computing，HPC）应用中具有广泛的应用场景. SpMV是许多科学计算和工程应用中的基本组成部分，尤其在求解大型稀疏线性系统^[2]时至关重要. 典型的应用包括流体力学、结构分析、电磁场仿真以及地质学中的地震波模拟等. 这些领域中的问题通常会被转化为大规模的稀疏矩阵形式，SpMV用于快速求算这些问题的解. 此外，SpMV在机器学习扮演着重要角色，尤其是在处理大规模稀疏数据和高维特征空间时. 随着数据量的不断增长，许多机器学习任务涉及到稀疏矩阵的运算，如协同过滤^[3]和支持向量机^[4]. 随着硬件架构的不断发展，多核/众核处理器、GPU，以及现代向量处理单元（如Intel AVX512和ARM NEON）为SpMV提供了加速的潜力，优化SpMV成为提升其性能的关键研究方向.

传统的SpMV实现通常需要开发人员手动进行优化，这包括选择合适的稀疏矩阵存储格式^[5]、调整计算策略^[6]、优化内存访问模式^[7]以及手动编写并行化代码等，而编译器优化能够消除手动调优的繁琐工作，同时避免由于开发人员的经验限制而带来的优化不足. 而传统编译器（如LLVM^[8]和GCC^[9]）的优化策略通常是针对稠密数据结构设计的，依赖静态分析技术进行循环展开、指令级并行和数据局部性优化等. 然而，由于稀疏矩阵的非零元素分布是高度动态的，传统的编译器优化方法无法在编译期对其进行有效的静态分析，从而导致编译器缺乏对不同稀疏模式的适配能力.

近年来，张量编译器和代码生成器（如TACO^[10]和TVM^[11]）的出现，使得开发人员通过几行领域特定语言（domain specific language，DSL）就能实现高性能算子的生成，然而它们通常只能针对特定模式进行优化，缺乏自适应调优^[12]的机制，无法根据不同类型的稀疏矩阵动态选择最优的数据存储格式和计算策略. 例如，某些稀疏矩阵可能更适合使用压缩稀疏列（compressed sparse column，CSC）格式或压缩稀疏行（compressed sparse row，CSR）^[13]格式，而另一些稀疏矩阵则可能需要块压缩稀疏行（block compressed sparse row，BCSR）格式或混合格式，然而现有的编译器通常无法自动识别并选择合适的存储格式，从而限制了计算性能的提升. 因此，要有效提升SpMV操作的性能，必须开发能够感知稀疏矩阵结构，并根据硬件架构和稀疏模式动态优化计算策略的编译器框架，这将能够最大化地发挥现代硬件的计算潜力.

针对上述传统张量编译器在处理SpMV的挑战，本文设计了SparseMode，用于高效地利用SpMV代码生成稀疏编译器. 本文的主要贡献有3点：

1）设计了新型的SpMV计算领域描述语言SpMV-DSL，简洁且高效地表达稀疏矩阵的结构和稀疏计算操作，能够自适应不同的稀疏模式，方便与编译器进行集成；

2）提出了基于稀疏模式感知的稀疏计算编译器，能够根据SpMV-DSL设定的存储格式和输入的稀疏矩阵非零元素分布，从而选择最优的计算策略和内存布局进行优化；

3）提出了生成高性能SpMV代码的流程，通过稀疏模式分析和SpMV-DSL指定的并行和向量化策略，生成高效的并行化SpMV算子代码，充分利用SIMD指令加速计算，同时减少冗余操作和内存访问.

1.   研究背景和研究动机

1.1   SpMV

SpMV的计算表达式为$\boldsymbol y=\boldsymbol {Ax}$，其中$\boldsymbol A$为稀疏矩阵，$\boldsymbol x$为输入的稠密向量，而$\boldsymbol y$为输出的结果，同样为稠密向量. 稀疏矩阵$\boldsymbol A$中存在大量零元，如果采用稠密矩阵的存储形式会造成大量内存浪费，通常$\boldsymbol A$采用压缩存储格式如CSR，如图1所示. CSR格式采用3个数组表示，$A.indptr$表示矩阵每一行的非零元在$A.idx$和$A.data$的起始位置，$A.idx$和$A.data$分别表示非零元列的位置以及存储的值，2个数组长度相同且一一对应.

图  1  CSR格式示意图

Figure  1.  Illustration of CSR format

下载: 全尺寸图片 幻灯片

算法1所示为基于CSR格式的SpMV实现. SpMV优化的难度在于对$\boldsymbol x$是间接访存，可能由于访存跨度过大，导致无法高效地利用缓存的局部性特性.

算法1. CSR格式的SpMV实现.

输入：CSR $\boldsymbol A[M,N]$，稠密向量$\boldsymbol x$；

输出：稠密向量$\boldsymbol y$.

① for $i= 1,…,M$ do

②　$sum= 0$；

③　for $j= A.indptr\left[i\right],…, A.indptr[i+1]$ do

④　　$sum+= A.data[j] \times \boldsymbol x[ A.idx[j\left]\right]$；

⑤　end for

⑥　$y\left[i\right]=sum；$

⑦ end if

1.2   张量编译器

随着深度学习模型的快速发展，张量编译器如TACO和TVM已经成为提升计算效率、加速深度学习模型部署的重要工具. 这些编译器通过实现自动化的高性能算子代码生成，极大地简化了算子优化过程. 张量编译器通常定义了一系列特定领域语言，用以精确描述算子的计算过程和优化逻辑. 开发人员只需编写几行DSL代码，便可生成高效的底层执行代码，从而避免了繁琐且容易出错的手动算子优化.

通过张量编译器，开发者能够在不深入底层硬件实现的情况下，针对特定硬件架构（如GPU、TPU^[14]等）进行优化. 这种方法不仅提升了算子执行效率，也降低了手动优化过程中的开发成本. 张量编译器利用数学模型、编译优化技术及自动调优算法，能够在多种硬件架构之间实现最佳的算子映射和调度策略，从而大幅提升计算性能.

1.3   研究动机

尽管张量编译器表现优异，但它们的设计主要针对的是稠密计算模式. 稠密矩阵的计算模式相对固定，数据访问和计算过程的并行性容易预测和优化. 针对这些特性，张量编译器能够有效地生成高度优化的代码，从而显著提升计算性能. 然而，在面对稀疏计算时，现有的张量编译器面临一定的挑战.

稀疏矩阵的计算模式与稠密矩阵的大不相同. 由于稀疏矩阵中大部分元素为零，非零元素通常分布不均匀，导致其计算模式高度不规则. 对于这种稀疏模式，现有的张量编译器往往无法很好地适应其复杂性. 尤其是在SpMV这类问题中，存在着大量的跳跃访问. 因此，传统的稠密计算优化策略无法直接迁移到稀疏计算任务中. 更为重要的是，当前的张量编译器通常缺乏对稀疏模式的深入分析与优化能力. 尽管一些张量编译器能够进行简单的内存优化和调度，但它们在稀疏矩阵的特殊模式下，往往无法自动识别出不同稀疏性带来的计算和存储优化机会. 缺乏对稀疏结构的细致分析使得编译器无法充分发挥稀疏矩阵的潜力，导致计算效率未能达到预期.

针对张量编译器在SpMV计算中的挑战，本文提出了高效的SpMV代码生成稀疏编译器SparseMode，其由SpMV-DSL和编译优化框架组成，其结构如图2所示. SpMV-DSL负责描述SpMV的计算逻辑，并提供包括格式选择、向量化和并行化调度指令等功能. 在编译阶段，编译框架根据SpMV-DSL中实现的计算和调度描述，以及稀疏矩阵的非零元素分布，执行稀疏模式分析等操作，并根据稀疏模式分析的结果进一步实现数据分块和调度优化，最终生成高性能的SpMV内核代码.

图  2  SparseMode编译优化框架

Figure  2.  Compiler optimization framework of SparseMode

下载: 全尺寸图片 幻灯片

2.   SpMV的领域特定语言设计

2.1   方案设计

SpMV的计算表达式为$\boldsymbol y=\boldsymbol {Ax}$. 受到TACO张量代数表达式的启发^[10]，为了准确描述SpMV的计算逻辑和稀疏矩阵$\boldsymbol A$的存储格式，本文设计了针对SpMV的领域特定语言SpMV-DSL，具体包括计算、格式和调度等描述形式.

1）计算描述

该部分主要描述SpMV的计算逻辑，根据计算过程中矩阵和向量的访存关系设计了如下的表达式：

2）格式描述

式（1）仅描述了计算中的访存关系，例如在计算$y\left[i\right]$时，需要遍历稀疏矩阵$\boldsymbol A$的第$i$行数据与向量$\boldsymbol x$，并进行规约求和. 然而，在实际运算过程中稀疏矩阵$\boldsymbol A$通常采用的是压缩存储的格式，如CSR或CSC格式，因此还需要额外的格式描述：

A.set_format(i,“dense”, j,“sparse”)，

A.set_format(i,“sparse”, j,“dense”).

对于需要设置稀疏格式的矩阵，set_format可以帮助矩阵设定每个维度的稀疏和稠密属性，例如(i,“dense”, j,“sparse”)表示维度$i$是稠密的和维度$j$是稀疏的，能够适配于CSR格式，同理(i,“sparse”, j,“dense”)适配于CSC格式.

3）调度描述

SpMV-DSL提供了向量化执行和并行优化的调度原语，用以实现高性能SpMV内核的生成. 对于并行执行，设计的调度为：

y.set_parallel(i, 2).

该调度表示在维度$i$的方向上进行并行化处理，且线程数设置为2. 而对于向量化执行，考虑到不同的SIMD架构所支持的向量化指令不同，用户可以设置底层向量化执行的操作：

y.set_simd_lane(2)；

y.set_vectorize(“gather”,“vmul”,“reduce”,“write”)；

y.set_vectorize(“gather”,“vmul”,“vadd”,“scatter”).

其中，set_simd_lane用于设置SIMD的数据通道数，对于AVX512指令集，底层在计算SpMV时可以采用2种向量化方式：行向量计算和列向量计算. 如图3所示，行向量计算对稀疏矩阵的某一行进行向量化计算，采用（“gather”,“vmul”,“reduce”,“write”）的运算模式，首先通过“gather”指令将指定位置的稠密向量$\boldsymbol x$加载到向量寄存器中，然后通过“vmul”实现向量乘法，并进行“reduce”规约求和，最终“write”写回到结果$\boldsymbol y$的对应位置. 而列向量计算对稀疏矩阵的多个列进行向量化计算，采用（“gather”,“vmul”,“vadd”,“scatter”）的运算模式，“gather”，“vmul”与上述行向量计算相同，列向量采用Element-wise的计算，因此采用“vadd”的向量加操作和“scatter”将结果写回到$\boldsymbol y$.

图  3  向量化计算模式

Figure  3.  Pattern of vectorized computation

下载: 全尺寸图片 幻灯片

2.2   与编译优化框架的关联

在本文设计的SpMV编译优化框架中，SpMV-DSL提供的格式和调度描述（如设置线程数、向量化操作等）用于指导编译框架进行合理的代码生成，而输入的稀疏数据，有助于编译框架进行稀疏模式分析，进而优化稀疏矩阵的存储和SpMV的向量化执行.

3.   SpMV的代码生成与优化

3.1   方案设计

SparseMode提供了针对SpMV算子的稀疏数据感知编译优化框架，其由稀疏模式分析、调度优化和代码生成等模块组成. 框架结构及各模块之间的关系如图2所示.

首先，根据SpMV-DSL的表达式描述，前端解析器会生成一个循环迭代树. 该迭代树的每个节点表示一个循环的变量，对迭代树进行先序遍历即可得到循环的顺序. 叶子节点表示相应的计算逻辑，根据设置的向量化调度原语，提取适合向量化的计算形式. 在节点j的叶子节点形成图中所示的向量化计算逻辑.

其次，针对多线程并行计算的需求，调度优化模块会根据SpMV-DSL中预定的线程数量以及输入的稀疏数据非零元的分布，对数据块进行合理划分，并优化调度策略，以确保每个线程分配到均衡的计算负载，从而提升计算资源的利用率和整体执行效率（图2中所示的稀疏矩阵分块和并行的线程数）.

最后，代码生成模块基于稀疏模式分析和调度优化的结果，输出高效的SpMV内核代码（图2中所示的2个内核代码，每个内核代码根据数据块的排布采用不同的向量化执行），并针对稀疏矩阵的存储方式进行优化，使得计算过程更加高效.

3.2   稀疏模式分析

为了生成高性能的SpMV代码，需要根据输入的稀疏矩阵非零元分布，划分多个可进行向量化运算的数据块. 假设SIMD的数据通道数为$d$，对于稀疏矩阵$\boldsymbol A$，可以将其划分为多个可向量化数据块$VP=$ $\{ [({r}_{1}^{1},{c}_{1}^{1}),…,({r}_{1}^{d},{c}_{1}^{d})],[({r}_{2}^{1},{c}_{2}^{1}),…,({r}_{2}^{d},{c}_{2}^{d})],…,[({r}_{V}^{1}, {c}_{V}^{1}),…,({r}_{V}^{d},$ ${c}_{V}^{d})]\}$，其中$\left[\right({r}_{i}^{1},{c}_{i}^{1}),…,({r}_{i}^{d},{c}_{i}^{d}\left)\right]$表示第$i$个向量化数据块中每个元素的行列位置. 如果${r}_{i}^{1}={r}_{i}^{2}=…={r}_{i}^{d}$，则表明该数据块采用行向量计算模式，用${B}_{i}^{r}$表示；如果$\forall j,k, {r}_{i}^{j}\ne {r}_{i}^{k}$，则表明该数据块采用列向量计算模式，用${B}_{i}^{c}$表示.

由于CPU的缓存具有局部性特征，即使在向量化执行过程中，同样要考虑数据读写的局部性. 对于稀疏矩阵$\boldsymbol A$中非零元的读取，因为在编译过程中能获取$\boldsymbol A$的非零元分布，可以提前将需要进行向量化计算的非零元存放在相邻位置. 然而，对于稠密向量$\boldsymbol x$和$\boldsymbol y$的读写是根据$\boldsymbol A$中的行列位置确定的，因此在挖掘$\boldsymbol A$中的向量化数据块时，要考虑对于$\boldsymbol x$和$\boldsymbol y$的访存跨度.

定义1. 稠密向量的访存跨度. 给定一个稠密向量$\boldsymbol x$，以及其访存位置的顺序为$P=[{p}_{1},{p}_{2},…,{p}_{T}]$，则该向量$\boldsymbol x$的访存跨度：

具体而言，对于第$i$个向量化数据块$\left[\right({r}_{i}^{1},{c}_{i}^{1}),…, ({r}_{i}^{d},{c}_{i}^{d}\left)\right]$，若该数据块采用行向量计算模式，则其对于稠密向量$\boldsymbol x$的跨度为$\displaystyle\sum _{t=2}^{d}|{c}_{i}^{t}-{c}_{i}^{t-1}|$，最后计算结果为标量，写回向量$\boldsymbol y$还需要一次访存，因此${S}_{{B}_{i}^{r}}=1+\displaystyle\sum _{t=2}^{d}|{c}_{i}^{t}-{c}_{i}^{t-1}|$；若该数据块采用列向量计算模式，除了对稠密向量$\boldsymbol x$的读取外，最后还要将计算结果写回向量$\boldsymbol y$，因此还需要加上$\boldsymbol y$的跨度，即${S}_{ {B}_{i}^{c}}=\displaystyle\sum _{t=2}^{d}|{r}_{i}^{t}-{r}_{i}^{t-1}|+|{c}_{i}^{t}-{c}_{i}^{t-1}|$. 综上，稀疏模式分析的优化目标为划分数据块$VP=\{{B}_{1},…, {B}_{V}\}$，${B}_{i}=\left\{\right({r}_{i}^{1},{c}_{i}^{1}),…,({r}_{i}^{d},{c}_{i}^{d}\left)\right\}$，最小化以下目标函数：

然而，直接求解最优的向量化数据块的时间复杂度是指数级的，且由于稀疏矩阵的规模通常较大，这在编译过程中会导致稀疏分析的开销过大. 稀疏矩阵中的非零元素分布通常是不规则的，且其稀疏性会因矩阵的规模、非零元素的分布和计算平台的硬件架构不同而有所变化. 为了避免这一瓶颈，本文提出了近似算法，旨在尽可能接近最优的向量化数据块选择. 通过对稀疏矩阵存储格式（如压缩稀疏行格式）的特性进行分析，该算法能够高效地选择适合SIMD计算的向量化数据块，具体如算法2~4所示.

算法2. 行向量块选择函数.

输入：CSR格式 $\boldsymbol A[M,N]$，SIMD通道数$d$，初始行$t$，访问标记$vis$，计数器$cnt$；

输出：行向量块$rcand$及其访存跨度${S}_{ 1}.$

① ${S}_{ 1}$←INF；

② if $cnt\left[t\right]\ge d$ then

③　for $i= A.indptr\left[t\right]\;{\rm{to}}\; A.indptr[t+1]$ do

④　　if $vis$.find ([t, A.idx[i]]) == false then

⑤　　　$rcand$.append ([$t, A.idx\left[i\right]]$)；

⑥　　end if

⑦　　if size ($rcand$) == $d$ then

⑧　　　${{S}}_{1}=\displaystyle\sum _{{j}=2}^{{d}}\left|{r}{c}{a}{n}{d}\right[{j},2]-{r}{c}{a}{n}{d}[{j}-1,2\left]\right|$；

⑨　　　break；

⑩　　end if

⑪　end for

⑫ end if

算法3. 列向量块选择函数.

输入：CSR格式$\boldsymbol A[M,N]$，SIMD通道数$d$，初始行$t$，访问标记$vis$，计数器$cnt；$

输出：行向量块$ccand$，访存跨度${S}_{ 2}$.

① ${S}_{ 2}\leftarrow INF,L\leftarrow t$, $R\leftarrow t$；

② while $L\ge 0$ and $cnt\left[L\right] > 0$

③　$L-=1$；

④ end while

⑤ while $R < M$ and $cnt\left[R\right] > 0$

⑥　$R+=1$；

⑦ end while

⑧ if $R-L\ge d$ then

⑨　从$[L,R]$中搜索最优的范围$[l,r](l\le t\le r, r-l=$ 　$d)$，使得每行中第一个符合要求的非零元 　组成$ccand$的访存跨度${S}_{ 2}$最小；

⑩ else

⑪　将$[L,R]$中每行第一个未被访问到的非零元 　加入到$ccand$；

⑫　在$[L,R]$之外，找到与第$L$和第$R$行距最近的 　$d+L-R$个行，并且每个行的$cnt > 0$；

⑬　将找到的$d+L-R$个行中第一个未被访问到 　的非零元加入到$ccand$；

⑭　计算$ccand$的访问跨度${S}_{ 2}$；

⑮ end if

算法4. 稀疏模式分析函数.

输入：CSR格式$\boldsymbol A[M,N]$，SIMD数据通道数$d$；

输出：向量化数据块$VP$.

① 初始化VP←$\varnothing$；

② 初始化集合$vis\leftarrow \varnothing$，用于标记已经访问过的 非零元位置；

③ 初始化数组$cnt$，使得$cnt\left[i\right]= A.indptr[i+ 1]-$ $A.indptr\left[i\right]$；

④ do while

⑤　$t\leftarrow \mathop{\rm{arg\,max}}\limits_{i}cnt\left[i\right]$,${S}_{ 1}\leftarrow INF$,${S}_{ 2}\leftarrow INF$；

⑥　$rcand= \left[\right]$,$ccand= \left[\right]$；

⑦　$rcand,{S}_{ 1}\leftarrow$算法2；

⑧　$ccand,{S}_{ 2}\leftarrow$算法3；

⑨　if ${S}_{ 1}$ < ${S}_{ 2}$ then

⑩　　$VP$.append ($rcand$),$cnt\left[rcand\left[:, 1\right]\right]-= 1；$

⑪　　$vis$.append ($rcand[:]$)；

⑫　else if ${S}_{ 1}$ $\ge$ ${S}_{ 2}$ && ${S}_{ 1}\ne INF$ then

⑬　　$VP$.append ($ccand$),$cnt\left[ccand\right[:, 1\left]\right] -=$1；

⑭　　$vis$.append ($ccand[:]$)；

⑮　end if

⑯ end while

⑰ 对剩余的非零元按照行/列向量填充零元，并 加入到$VP$；

⑱ end while

算法2和算法3分别介绍了如何在稀疏矩阵$\boldsymbol A$的第$t$行中分别找到行向量块和列向量块，而算法4为稀疏模式分析的方法，在执行过程中会调用算法2和算法3，并判断选择哪种向量化执行方式. 具体而言，算法2的思路比较直接，如第$t$行存在大于$d$数量的非零元，则将该行中邻近的$d$个非零元加入到行向量块，该算法的时间复杂度为O(nnz)，nnz为第$t$行的非零元数目.

由于列向量需要考虑读和写的访存，算法3相较于算法2更为复杂，如图4所示，算法3首先围绕着第$t$行搜索范围最大的连续行$[L,R]$，如果$R-L\ge d$，则直接在该范围内找到能最小化${S}_{ 2}$的连续$d$个行，并将它们的第一个符合要求的非零元加入到列向量块中（行⑨），否则先将$[L,R]$行符合要求的非零元添加到列向量块，剩余的$d+L-R$个非零元则从$d+L-R$个与$[L,R]$行距离最近的行中被找到（行⑪~⑭）. 该算法的时间复杂度为O($N$).

图  4  列向量块选择

Figure  4.  Selection of column vector block

下载: 全尺寸图片 幻灯片

为了防止非零元被计算多次，在挖掘向量块时，算法2和算法3都采用了辅助数据结构$vis$，用于标记已经添加的非零元. 算法4根据算法2和算法3输出的访存跨步${S}_{ 1}$和${S}_{ 2}$进行比较，判断当前的$d$个非零元为行向量块还是列向量块. 在算法4的行⑰中，如果当前找不到向量块，但还存在未处理的非零元，则需要将它们进行零元填充，以适应于向量化处理. 对于未处理的非零元的分布，最多存在$d-1$个行和$d-1$个列（如果有大于等于$d$的行或列，则会被加入到向量块中），本文以图5为例说明零元的填充方法.

图  5  向量块零元填充方法

Figure  5.  A method of zero filling for vector block

下载: 全尺寸图片 幻灯片

该填充方法与算法2和算法3相似，只是不需要再考虑能否将SIMD的数据通道填满非零元，该方法首先选取非零元数量最多的行，根据算法2和算法3的相同思路，计算能获取的行向量（图4中{a,b,c}）和列向量（图5中{a,e,f}）的访存跨度，根据它们的大小选择对行向量或者列向量进行零元填充. 此时选择了{a,e,f}填充为列向量块，但填充零元的位置需要结合向量化运算时“gather”和“scatter”的访存，选择跨度最小的位置填充，如图中（4,1）的位置. 接下来，还有{b,c}要进行填充，经分析，在（2,3）和（2,4）的位置进行填充最优，并将{b,0,0,c}作为行向量块.

3.3   调度优化

在现代CPU架构中，通常采用多核（multi-core）架构，因此在优化过程中还需要考虑如何支持多线程并行处理的需求. 在经过稀疏模式分析后，可以得到多个可以向量化处理的数据块$VP$. 调度优化目标是将向量数据块$VP$按照设定的线程数进行划分，使得每个线程尽可能获得均匀的计算负载. 该优化方法如算法5所示，首先对$VP$中的数据块进行升序排列，由于每个数据块的长度为$d$，对于任意2个数据块，可以先比较第1个元素的行列位置，如果相同就比较第2个，依次类推. 算法5先按线程数对$VP$进行平均划分，然而在多线程运行过程中，考虑到不同线程分配的向量块可能存在相同的行，在写回到结果向量$\boldsymbol y$过程中，可能存在线程冲突的问题，如图6所示. 为了避免不同的线程同时写向量$\boldsymbol y$中相同的列，算法5先遍历矩阵的每个行$r$，再根据每个${VP}_{i}$是否存在行位置等于$r$的数据块，将这些${VP}_{i}$和冲突的数据块记录下来，从记录的${VP}_{i}$选择最小的数据块数量$t$，将有冲突的数据块由$t$进行处理（行④~⑨）.

图  6  避免线程冲突

Figure  6.  Avoidation of thread conflict

下载: 全尺寸图片 幻灯片

算法5. 稀疏模式分析函数.

输入：向量化数据块$VP$，线程数$T$，稀疏矩阵行维度$M$；

输出：$T$个互不相交的数据块$\{{VP}_{1},{VP}_{2},…,{VP}_{{T}}\}$，且${VP}_{1}\cup …\cup {VP}_{{T}}$=$VP$.

① 对$VP$按每个数据块的行列位置升序排列；

② $avg\leftarrow \left|VP\right|/T$；

③ 按照每个${VP}_{i}(1 \le i\le T)$的数据块数量为$avg$进 行均匀划分；

④ for r = 1, 2, …, M do

⑤　$\mathrm{\omega }\leftarrow$存在数据块的行位置等于r的${VP}_{i}(1 \le$ 　$i\le T)$；

⑥　$t\leftarrow {\rm{arg\,min}}\left|f\right|, \forall f\in \omega$；

⑦　$\forall f\in \omega ,f\ne t$，将$f$中所有行位置等于r的数 　据块去除，并加入到$t$中；

⑧ end for

此外，在向量化执行的过程中，同一线程在计算不同的向量块时可以进行融合. 线程在计算{i, j, k, l}和{m, n, o, p}时，由于它们计算结果写入向量$\boldsymbol y$的相同位置，可以使用同一个向量寄存器对它们的结果进行累加，最后一起写入$\boldsymbol y$，从而减少访存次数，过程描述如图6所示.

3.4   代码生成

根据3.1~3.3节分析，最后代码生成模块会产生向量化的多线程SpMV内核，以及将稀疏矩阵转化为向量数据块VP并存储到磁盘. 在运行时，主线程先将转化后的VP加载到内存后创建多个工作线程，每个线程调用相应的SpMV代码，并对其负责的${VP}_{i}$进行独立计算.

4.   实　　验

4.1   实验配置

1）实验平台. 基于 Intel Xeon、 Intel Core和鲲鹏920开展实验对比. 其中 Intel Xeon 平台包含16个处理器核，每个核可开启2个线程，最多有32个线程，支持 AVX 512向量指令集，内存为 512 GB；Intel Core 平台包含 14 个处理器核，每个核可开启2个线程，最多可开启28 个线程，支持 AVX 512 向量指令集，内存为 64 GB. 鲲鹏920包含64个处理器核，每个核有1个线程，最多可开启64个线程，支持ARM NEON向量指令集，内存为256 GB，表1列出本文实验平台的详细参数.

表  1  实验平台配置参数

Table  1.  Configuration Parameters of Experimental Platforms

配置/参数	Intel Xeon	Intel Core	鲲鹏920
处理器	Silver 4110	i5-13600KF	KP920
核心数	16	14	64
每核心线程数	2	2	1
L1 Cache	512 KB	1.2 MB	8 MB
L2 Cache	16 MB	20 MB	64 MB
L3 Cache	22 MB	24 MB	256 MB
内存大小	512 GB	64 GB	256 GB

下载: 导出CSV 

| 显示表格

2）数据集. 本文从稀疏矩阵集合中选取10个有代表性的稀疏矩阵，覆盖如最优化、流体力学、线性规划等多个领域，其规模和非零元素分布模式各异. 详细信息如表2所示.

表  2  稀疏矩阵数据集

Table  2.  Sparse Matrices Datasets

数据集	维度	非零元数量
com-Orkut	3 072 441	234 370 166
dblp-2010	326 186	1 615 400
circuit5M	5 558 326	59 524 291
road_central	14 081 816	33 866 826
soc-Pokec	1 632 803	30 622 564
nv2	1 453 908	37 475 646
StocF-1465	1 465 137	21 005 389
flickr	820 878	9 837 214
mc2depi	525 825	2 100 225
webbase-1M	1 000 005	3 105 536

下载: 导出CSV 

| 显示表格

3）对比方法. 本文将目前最先进的张量编译器TACO和TVM，以及Intel的MKL数学库作为对比对象，并对生成的代码统一使用GCC编译器，“-O3 -fopenmp”的编译参数生成可执行程序. 在鲲鹏920的实验中，由于ARM NEON不支持 “gather”和“scatter”指令，在稀疏模式分析过程中，通常只能在每一行中寻找连续的非零元素，或对行进行填充以优化其结构. 考虑到实验过程中的数据波动，为了更准确地评估SpMV在不同数据集上的性能，本文采用每个数据集运行100次后的平均时间作为性能指标. 同时，还对单精度和双精度的SpMV进行了测试，进一步验证这些方法在不同精度下的表现.

4.2   SpMV性能分析

图7所示为不同方法在多线程场景下（Intel Xeon 64线程，Core 28线程，KP920 64线程）的性能比较. 本文使用SparseMode相对于TACO、TVM和MKL的加速比作为实验的可视化结果. 红色的柱状图表示SparseMode对TVM的加速比，绿色的柱状图表示SparseMode对 TACO的加速比，蓝色的柱状图表示SparseMode对MKL的加速比. 在Intel平台上，SparseMode 与TVM 的比较中，表现出较为显著的加速比. 在多数数据集中，SparseMode对TVM的加速比都在1.5以上（如 nv2 和 StocF-1465），表明 SparseMode 相较于 TVM 具有显著的性能提升. SparseMode 相较于TACO的加速比也显示出良好的性能提升，尤其是在某些较大或更稀疏的矩阵数据集上. 例如，在webbase-1M和StocF-1465数据集上，SparseMode相较于TACO取得了更高的加速比. 而在flikcr数据集上，SparseMode的性能较差，主要是flickr相较于其他数据集稀疏度较低. 此外，我们对比了Intel上的MKL函数库，针对大多数稀疏矩阵，SparseMode对MKL表现出较显著的加速，最高可达2.52.

图  7  SpMV性能分析

Figure  7.  Performance analysis of SpMV

下载: 全尺寸图片 幻灯片

对于单精度和双精度运算，SparseMode的性能表现也有所不同，主要原因在于SparseMode在计算过程中将非零元素存储在向量寄存器中. 对于AVX 512指令集来说，每个向量寄存器能够存储的双精度数据的数量少于单精度数据. 因此，双精度运算的性能略低于单精度运算. SparseMode对TACO，TVM和MKL的加速效果在不同数据集上有所不同，这可能与数据集的稀疏性、矩阵规模以及计算模式等因素相关. 对于一些稀疏矩阵特别适合的场景（稀疏矩阵维度与非零元数量相近，如 webbase-1M 和 mc2depi），SparseMode能够显著提高计算效率.

在鲲鹏920平台上，SparseMode的性能要低于其在Intel平台的表现，主要是因为测试的鲲鹏920平台的ARM NEON指令集没有“gather”和“scatter”指令，导致在稀疏模式分析过程中只能在每一行中找连续的非零元素或对行进行填充，因此会造成大量的零元素填充，增大计算开销. 此外，ARM NEON指令集中每个向量寄存器的长度要小于AVX 512指令集，因此单条指令能处理的数据量较小.

4.3   扩展性分析

为了探究SparseMode的可扩展能力，本文在设置的不同线程数量条件下，分析SparseMode与TVM，TACO和MKL的性能对比. 本文选取了在4.2节中SparseMode性能表现最优的mc2depi和webbase-1M的数据集，并使用单精度运算进行测试实验结果，如图8所示.

图  8  SpMV扩展性分析

Figure  8.  Scalability analysis of SpMV

下载: 全尺寸图片 幻灯片

具体而言，在Intel Xeon平台选取了{1, 4, 8, 16, 32}的线程数，而在Intel Core平台选取了{1, 2, 4, 8, 16}的线程数. SparseMode在单线程情况下的性能低于TVM和MKL，然而在大多数情况下的性能还是优于TVM，TACO和MKL，尤其在较高线程数时，SparseMode表现出了更强的加速性能. 从实验结果上看，SparseMode表现出了较好的扩展性.

5.   相关工作

SpMV^[15]的优化已有多年的研究历史. ALBUS^[16]设计了简单的逐行矢量计算机制，将单行内的非零元素依次加载到矢量寄存器中进行计算；VHCC^[17]通过连接多行来提高向量计算效率，同时将矩阵划分为2个维度，以增强数据的局部性；CVR^[18]将非零行顺序加载到SIMD通道中进行计算，同时对SIMD通道宽度的行进行处理，可以大大提高矢量计算效率；BCSR^[19]通过划分密集子块，将稀疏矩阵计算转化为整洁的数据计算模式；BCCOO^[20]采用分区策略，利用位压缩行索引存储，缓解内存压力；LAV^[21]为遵循幂律分布的矩阵设计了行/列分区机制，提高了数据访问局部性； 有部分研究者探索在GPU上实现高效的SpMV^[22]以及不同稀疏格式的性能分析^[23-24].

现代编译器框架，如MLIR^[25]等在稀疏矩阵运算优化方面已经取得了一定进展，它们通常采用符号化表示、自动调度、存储优化和向量化等方法来优化稀疏计算. 然而，这些框架在处理复杂和不规则的稀疏模式时仍然存在局限性，难以自动感知不同存储模式，并为特定硬件架构生成最优代码.

TACO^[10]是专门为稀疏张量计算设计的编译器框架，主要用于自动生成优化的稀疏矩阵运算代码. 其核心方法包括：稀疏迭代表示（sparse iteration representation）、自动格式转换、代码生成优化、自动向量化. 但它缺乏对不规则稀疏模式的自适应优化，对硬件优化的适应性较差，其SIMD向量化能力不足. TVM^[11]是用于深度学习编译优化的框架，但其计算调度（scheduling）机制也可以用于稀疏计算优化. 但缺乏直接的稀疏模式感知能力、不具备对复杂稀疏模式的自动识别能力以及对不规则稀疏模式支持有限. MLIR^[25]在编译优化领域广泛应用，也开始支持稀疏计算的优化. 但其对硬件优化的适配能力不足，MLIR的Sparse Tensor Dialects仍需要用户手动指定稀疏模式，而非自动分析数据分布情况.

近年来，学术界和工业界针对不同稀疏模式的优化提出了模式分析、存储优化、并行计算和SIMD 加速等多种方法. 稀疏感知集合归属于归一化最小均方算法SM-NLMS^[26]，用对数和log-sum函数取代ZASM-NLMS中的ℓ₁范数惩罚项，应用于稀疏信道估计和回声消除. Kumar^[27]等人设计了自适应信号处理方案来利用系统脉冲响应的稀疏性质的自适应滤波方案. SPARQ^[28]在不同的表示粒度下利用非结构化和动态激活稀疏性，实现了较小的精度降低和实用硬件. DZiG^[29]以稀疏感知的增量处理技术进行递归方式表达计算，能够安全地识别和修剪更新，它在保证BSP语义的同时，保持了稀疏计算的效率. WACO^[30]通过考虑稀疏性模式、格式和调度，设计基于神经网络的成本模型，来预测稀疏张量程序的运行，从而选择合适的稀疏格式和调度策略.

6.   结　　论

本文深入分析了SpMV计算过程中稀疏模式带来的性能优化瓶颈，并探索了矩阵稀疏导致的不规则内存访问及难以应用向量化加速等挑战. 基于此，本文构建了基于稀疏模式感知的编译优化框架SparseMode，用于高效的SpMV代码自动生成. 该方案通过设计SpMV领域描述语言SpMV-DSL来简洁地表达稀疏矩阵结构及计算操作，并能够自适应不同稀疏模式. SparseMode结合稀疏模式分析与编译优化，动态选择最优的计算策略和内存布局，最大化硬件性能. 通过并行化和向量化策略生成高效的SpMV算子代码，充分利用SIMD指令加速计算. 经实验分析论证，该方案相较于张量编译器如TACO和TVM，在提升SpMV计算性能方面具有显著优势.

作者贡献声明：王昊天提出方案和撰写论文；丁岩构建实验并提供论文思路；何贤浩实现代码并完成测试；肖国庆改进实验；阳王东指导论文写作并修改论文.