稀疏矩阵向量乘（sparse matrix-vector multiplication，SpMV）作为数值计算中的基础算子， 被广泛应用在科学计算及各类工程应用中，包括但不限于图计算、机器学习、流体力学和金融数学等. SpMV将一个矩阵A和向量x相乘得到目标向量y（y=Ax）. 其中矩阵A是稀疏矩阵，x和y是稠密向量. 系统中的很多问题都会被抽象成SpMV问题，譬如迭代法求解稀疏线性方程组^[1]等. 随着矩阵维度的不断增大，SpMV的计算成本随之增高，成为很多系统中的性能瓶颈. SpMV的计算效率亟需提高. 另一方面，硬件架构在持续演进，多核/众核架构^[2-3]和向量处理能力（如Intel AVX512^[4]和ARM NEON^[5]）等在提升应用性能方面展示了巨大的潜力. 然而，稀疏矩阵的不规则性使SpMV无法有效利用这些先进的硬件特性.

SpMV的向量优化面临诸多问题，例如：大量的短向量（单行内非零元较少）导致向量计算不饱和； 向量寄存器把值写回数组y时，会因为数据不整齐造成写冲突问题等. 现有SpMV向量优化方法主要分为3种：逐行向量化（如Intel MKL^[6] 和 CSR5^[7]等）、多行向量化（如CVR^[8]等）和零元填充（如SELL-C-σ^[9]等）. 实验分析表明这3种向量优化方法，无法在特征各异的稀疏矩阵上皆取得最优性能. 逐行向量化方法，适合处理长向量（单行内非零元较多）且访存跨度不大的稀疏矩阵；多行向量化方法，同时从多行加载元素，适合处理相邻行的非零元数量和分布模式较为接近的稀疏矩阵；零元填充作为前2种向量化方法的补充，难点在于填充比例和向量优化效果中间的权衡. 综合来看，这3种方法各有优势，可用于应对具备特定非零元分布特征的稀疏矩阵，但尚无法较好处理非零元分布较为复杂和多样化的稀疏矩阵.

另一方面，向量处理架构已经成为当前提升数据处理性能的主要方向之一. 因设计目标和服务场景各异，不同向量硬件间虽有较大的向量指令重合，但也通常会引入一些特定的向量特性和指令. 这为研究通用的SpMV向量优化方法带来了巨大的挑战.

故而，把极不规则的稀疏矩阵映射到极度规则的向量硬件上进行计算，需首先解决稀疏矩阵的规则化问题，以及向量硬件计算能力的通用化表示问题. 同时，为达成最优的SpMV性能，需在提升向量化性能的同时，不影响多线程并行执行性能.

基于此，本文提出一种混合向量化的SpMV优化机制（hybrid vectorization-optimized mechanism of SpMV，HVMS），以应对复杂和非零元分布较为多样化的稀疏矩阵. HVMS首先分析向量硬件，抽象其基本操作，例如常用的规约求和操作、前向求和操作、零元填充操作和掩码运算操作，以及由基础操作组合而成的较高级向量操作，如前缀求和等. 其后，HVMS设计多条规则，指导稀疏矩阵进行规则化转换，弱化其不规则的非零元分布特征为“规则（regular）”或“弱不规则（less irregular）”. 此处需要指出的是，HVMS的每种规则对应一种指定特征的非零元分布模式，并可由对应的向量操作完成计算. 由此，SpMV的向量优化问题转化为一个较规则的稀疏矩阵到一系列基础向量操作的映射问题.

最后，本文基于Intel Xeon，在30个稀疏矩阵数据集上测试HVMS的性能和有效性. 实验结果表明，相比3种代表性方案CVR，SELL-C-σ，Intel MKL，HVMS可分别获得1.60倍、1.72倍和1.93倍的平均加速比. 同时，本文分析HVMS的可扩展性，随着线程数量从1增加至64，HVMS相比单线程Intel MKL的平均加速比可达20.35倍，分别高于其他3种方案10.93倍、11.91倍和12.76倍.

1.   研究背景和研究动机

1.1   基于CSR格式的SpMV算法

基于CSR （compressed sparse row）格式的SpMV算法是基于CSR格式的SpMV代码实现的. 如算法1所示. 图1是稀疏矩阵A的CSR格式存储的示例描述，开辟3个数组来存储稀疏矩阵： csrVal，csrColIdx，csrRowPtr. 其中数组csrVal存储非零元素值；数组csrColIdx存储非零元素的列索引；数组csrRowPtr存储每行第1个非零元素在数组csrVal中的偏移量. 本文使用CSR格式作为基本存储格式.

图  1  CSR格式示意图

Figure  1.  Illustration of CSR format

下载: 全尺寸图片 幻灯片

算法1. 基于CSR格式的串行SpMV算法.

输入：稀疏矩阵A（m×n），稠密向量x；

输出：稠密向量y.

① for i = 0 ：m do

② 　sum = 0；

③ 　for j = csrRowPtr[i] ：csrRowPtr[i+1] do

④ 　　sum += csrVal[j]×x[csrColIdx[j]]；

⑤ 　end for

⑥ 　y[i] = sum；

⑦ end for

算法1逐行遍历稀疏矩阵A，每次计算矩阵的1行与数组x的点积（行③~⑤），并将结果写回数组y.

1.2   向量优化技术

SIMD（single instruction multiple data）^[10]即单指令多数据流，即一条指令处理多个数据. 自问世以来，SIMD技术得到了蓬勃的发展，Intel的SSE系列^[11]、AVX系列^[11]、ARM 的 NEON^[5]系列以及AMD的3D Now! ^[12]指令集是其中的代表. 从AVX2（256 b）到AVX512（512 b），通过扩大向量指令的宽度，Intel 处理器在数据级并行（data level parallelism，DLP）方面取得了不错的性能，AVX512指令可同时处理8个双精度浮点类型或者16个单精度浮点类型的数据，表1列举了部分常用的向量指令操作.

表  1  部分常用的向量指令

Table  1.  Some Common Vector Instructions

指令操作	指令行为说明
load(&addr)	将数据（例如双精度浮点类型的数据）从内存地址addr处加载到寄存器.
fmadd_pd(Reg1, Reg2, Reg3)	执行向量乘加操作，即目标寄存器的值=Reg1×Reg2 + Reg3.
gather(idx, &array)	使用idx寄存器提供的索引值从数组array中gather多个数据（例如8个double数据），并将其加载入寄存器中.
scatter(&array, idx, Reg)	使用idx寄存器提供的索引值从Reg中scatter多个数据（例如8个double数据），并将其存储到内存数组array.
reduction_add(Reg)	返回将寄存器Reg内所有值相加后的结果.
mask_reduction_add(Mask, Reg)	利用给定的掩码值（Mask）将寄存器Reg内相应位置的值相加并返回.

下载: 导出CSV 

| 显示表格

因设计目标和服务场景各异，不同的向量硬件如Intel AVX，AMD 3D Now!，ARM NEON，GPGPU^[13]等的向量特性和指令虽有重合，但也有较大区别. 对向量硬件进行抽象，支持通用SpMV向量优化方法，是很重要的.

为便于表述，本文约定向量寄存器宽度为ω，且ω是指向量寄存器可存储的非零元数量，亦即如果寄存器宽度为512 b. 若稀疏矩阵中元素为整型（int）或单精度浮点型（single float），则ω= (512 b/32 b) =16；若稀疏矩阵中元素为双精度浮点型（double float），则ω= (512 b/64 b) = 8. 鉴于ω的值对表述向量优化机制并无影响，本文为便于表述和图表展示，设定 ω=4.

1.3   研究动机

为研究通用的SpMV向量优化方法，本文首先分析现有SpMV向量优化方法的技术路线和性能表现. 图2描绘了3种常用的SpMV向量化策略及其代表方案：逐行向量化、多行向量化和零元填充.

图  2  稀疏矩阵示例以及Intel MKL，CVR，SELL-C-σ在计算时的对应数据排布

Figure  2.  A example of sparse matrix and the corresponding data layouts of Intel MKL, CVR and SELL-C-σ in the calculation

下载: 全尺寸图片 幻灯片

1）基于CSR格式的逐行向量化. 代表性工作有Intel MKL库中的csrSpMV算子实现. 该算子对CSR进行自动向量优化，每次对单行内的ω个非零元进行自动向量化；不足ω个则执行标量计算.

2）同时计算多行的向量化实现. 代表性工作有 CVR等. CVR按序同时读入ω行，从每行各取一个元素加载入ω个SIMD通道中.

3）基于零元填充的向量化方法. SELL-C-σ将矩阵按照行内非零元数量进行排序，并采取零元填充的方式来达到数据对齐的效果，目的是确保同时计算的ω行也可同时写回.

如上3种SpMV优化方法各具特点且分别代表了一类SpMV优化思路. 本文选择4个稀疏矩阵数据集，评测以上3种方案的性能，如图3所示. Intel MKL，CVR，SELL-C-σ 这3个方案分别在不同的稀疏矩阵数据集上取得了最优性能，亦即没有一种方案可以在全部4个稀疏矩阵上都取得最好性能.

图  3  Intel MKL，CVR，SELL-C-σ在4个数据集上的性能结果

Figure  3.  Performance results of Intel MKL, CVR and SELL-C-σ on four datasets

下载: 全尺寸图片 幻灯片

为分析其中原因，本文描绘了4个矩阵的非零元分布密度图，如图4所示，不同颜色代表不同密度，黄色为密度最高部分. 4个稀疏矩阵的非零元分布差异极大. 因稀疏矩阵非零元分布的不规则性，以及不同稀疏矩阵非零元分布的差异性，使得SpMV很难借助单一优化手段在差异较大的稀疏矩阵上都取得较好性能. 一般来讲，单一优化方法如CVR和SELL-C-σ等仅对稀疏矩阵符合特定规则或矩阵中的一些行符合特定规则时较为有效.

图  4  4个典型矩阵数据集的密度分布

Figure  4.  Density distribution of four typical matrix datasets

下载: 全尺寸图片 幻灯片

基于如上分析，本文将研究一种SpMV通用向量优化机制以应对多样化的稀疏矩阵非零元分布模式和差异化的向量硬件特性等，把SpMV映射到规则的向量硬件上进行计算.

2.   HVMS方案设计

为使用极度规则的向量化计算方法来提升极不规则的稀疏矩阵运算问题（如SpMV）的性能，本文提出一种混合向量化的SpMV优化机制HVMS，图5刻画了HVMS的设计思想，主要分为3部分：

图  5  HVMS的设计思想

Figure  5.  The design philosophy of HVMS

下载: 全尺寸图片 幻灯片

1）对硬件的向量计算能力进行抽象建模，转化为规约求和操作（reduction-based operation）、前向求和操作（forward-sum operation）、零元填充操作（zero-padding operation）和掩码运算操作（mask-based operation）等基本操作. 鉴于不同的向量化硬件，如AVX512，ARM NEON，GPGPU等支持的向量化操作有较大重合但也各有差异，HVMS针对不同硬件建模后的向量化操作集合是不同的. 本文主要以AVX512为例进行建模，但本文方法可以扩展到其他向量硬件.

2）鉴于SpMV计算性能受稀疏矩阵A的非零元分布模式的影响较大，HVMS需着重处理稀疏矩阵A的不规则问题，把极度不规则的稀疏矩阵，通过统计、切块和分类等方式，进行规则化转换（regularization transformation）. 亦即，结合向量硬件的特性，把不规则的稀疏矩阵分解为若干相对规则的部分.

3）把稀疏矩阵的计算（如SpMV），按照不同的方式映射到向量计算单元上，由此达到最优的向量化性能，最大程度利用当前先进体系架构AVX512等强大的向量计算能力.

区别于现有工作，HVMS并不关注单独针对特定稀疏矩阵的非零元分布模式进行优化，而是通过抽象建模的方式，把稀疏矩阵运算按照指定规则映射到指定的向量计算特性上. 故而，其他稀疏矩阵运算的向量化方法，以及新增的向量特性或自定义指令，都可以定义为规则集成到HVMS中.

2.1   向量硬件的计算能力抽象

很多不同硬件的向量处理能力，例如Intel AVX2，AVX512，ARM NEON，GPGPU等，都可通过SIMD的方式提供向量处理能力从而提升性能，但在具体的特性实现上存在较大差异. 体现在硬件特性上，这些差异包括但不限于向量处理位宽，以及对load，store，gather，scatter，reduction，mask等指令及其变种的支持. 例如：mask_reduction指令可以对向量单元里的指定元素进行规约操作；gather指令可以从内存中加载非连续数据到向量寄存器；permute和shuffle指令可以按照指定要求重新排布向量单元内的元素位置等. 这些指令特性可以增加设计稀疏矩阵向量优化方法的灵活性. 对向量计算单元的能力进行抽象和建模，首先需要分析指定硬件的向量指令或特性，将其中有助于实施稀疏矩阵运算的特性抽象为基本操作. 因不同的向量计算单元支持的特性各异，故而抽象后的模型也会有较大区别，例如AVX2仅支持gather， AVX512则同时支持gather和scatter.

HVMS将以Intel AVX512为示例硬件，对稀疏矩阵运算的过程进行建模. 该模型中将覆盖大部分通用向量计算特性，如fmadd，gather等及部分Intel AVX512 的专有特性，如mask_reduction等. 该模型具体分为4个部分：

1）规约求和. 如图6（a）所示，reduction指令对向量寄存器内的所有元素求和，是稀疏矩阵运算中最常用的向量操作之一. 规约的元素数量取决于向量寄存器的宽度ω. 当稀疏矩阵以CSR格式进行运算时，逐行计算并规约后写回到数组y，是常见的向量优化方法. 规约求和一般用于处理矩阵中同一行的非零元. 当行内非零元数量是 ω 的整数倍时，可调用reduction指令，把运算结果写回数组y指定位置.

图  6  4种向量计算基本操作

Figure  6.  Four basic operations of vector computing

下载: 全尺寸图片 幻灯片

2）前向求和. 如图6（b）所示，向量寄存器读入ω个非零元计算后，顺序读入另外ω个非零元继续计算，每次计算的结果自动求和. 该过程通过使用 fmadd指令完成. 前向求和在CSR格式的稀疏矩阵运算中通常会同时读入多行的非零元，计算后把多个结果同时写回数组y.

3）零元填充和掩码运算. 如图6 （c）所示，为使得数据较为整齐，方便进行向量计算，通常会填充部分零元. 同样，图6 （d）的掩码运算在应对分支语句和数据不规则等影响向量化性能的问题上也有较大帮助. 通常，在稀疏矩阵运算中，每行的非零元数量并不是ω的整数倍，故而零元填充和掩码运算便成为了向量优化方法的重要辅助手段.

HVMS是把不规则的稀疏矩阵运算映射到规则的向量计算单元的机制和方法，可以覆盖到各类向量计算单元，例如Intel SSE，AVX2/AVX512，AMD3D Now!，ARM NEON，GPGPU以及自定义指令的向量加速器等. 本文为简化方案描述，仅阐述部分通用基础操作的建模过程和方法. 但除此之外，HVMS可以覆盖更为复杂且特性各异的向量化操作，以及诸如前缀求和、CSR5和CVR等高级向量化方法.

2.2   稀疏矩阵的规则化转换

稀疏矩阵极不规则的非零元分布，使得稀疏矩阵运算的访存过程和计算过程也都极不规则. 其中，不规则的访存过程，亦即非连续的向量存取操作，可以借助向量单元中的gather 和scatter 等操作完成. 但是，不规则的计算过程为稀疏矩阵运算的向量优化带来了额外的阻碍. HVMS对稀疏矩阵进行规则化转换、重新排布稀疏矩阵的非零元、弱化甚至消除稀疏矩阵运算的不规则性，使其更加适合利用规则向量化硬件的计算能力.

如图7所示，对照向量计算硬件建模后的操作类型，以规约求和、掩码运算、前向求和以及零元填充为例，HVMS对稀疏矩阵中满足指定要求的行进行统计和分类. 具体规则为：

图  7  稀疏矩阵的规则化转换

Figure  7.  Regularization transformation of sparse matrix

下载: 全尺寸图片 幻灯片

1）规则1. 对应规约求和操作. 统计非零元数量为向量位宽ω整数倍的行，亦即满足规则nnz_row[i]%ω== 0的行，标记为使用规约求和操作计算.

2）规则2. 对应掩码运算操作. 统计2行可以同时进行运算的情况，亦即满足规则nnz_row [i]+ nnz_row [j] == ω的2行，标记为合并后使用掩码运算.

3）规则3. 对应前向求和操作. 统计非零元数量相等的ω行，亦即满足规则nnz_row [i] == nnz_row [j] == nnz_row [k] ==…== nnz_row [l]的ω行，标记为使用前向求和操作计算.

4）规则4. 对应零元填充. 对于适合规约求和操作和前向求和操作但是略带偏差的行. 例如，填充比例在阈值以下时（阈值可由用户指定），可判定为零元填充，同时归入规则1、规则2或规则3.

HVMS对稀疏矩阵进行规则化转换所遵循的规则，取决于向量化硬件建模后的操作类型，亦即由向量化硬件可以提供哪些类型的向量化方法来决定. 故而，随着HVMS覆盖更多向量特性和硬件类型，稀疏矩阵的转换规则也会随之丰富. 值得指出的是，HVMS通常会设定好规则间的优先级，以解决指定行满足不止1条规则时的冲突问题. 同时，本文也将在实验部分分析，当规则间优先级不同可能会带来的性能差异.

2.3   稀疏矩阵运算到向量计算单元的映射

稀疏矩阵A被规则化并分解为多个块（block），且块内各行的非零元分布规则一致. 符合指定规则的块会被分配，由对应的向量操作来执行. 具体如图8所示，0~3行属于单行非零元数量为ω倍数的行，交由规约求和操作计算；4~7行每行的非零元数量相等，故可同时读入多行非零元，交由前向求和操作计算；8~11行经零元填充后交由前向求和操作计算；12行经零元填充后交由规约求和操作计算；13和14行合并后元素数量为ω，可使用掩码进行计算；15行只有2个非零元素，可执行标量计算，无需映射到向量计算单元.

图  8  转换后的稀疏矩阵到基本向量计算操作的映射

Figure  8.  Converted sparse matrix maps to the basic operations of vector computing

下载: 全尺寸图片 幻灯片

该映射关系并不是固定的，例如8~11行在进行零元填充后，同时满足规约求和及前向求和2个操作的要求. 此时，这4行的非零元分布以及执行指令的延迟等因素都可能会影响性能. 通常在此种情况下，HVMS采用默认优先级顺序执行，即规约求和→掩码运算→前向求和→CSR计算.

3.   HVMS的实现及其SpMV计算方法

HVMS把非规则的SpMV运算映射到规则的向量硬件，其方法可覆盖到非零元分布各异的稀疏矩阵，以及不同特性的向量硬件. 本节描述HVMS的实现细节，首先介绍稀疏矩阵的规则化数据转换过程，然后阐述基于HVMS进行SpMV计算的过程.

3.1   稀疏矩阵的数据转换

HVMS的输入矩阵可以是任意格式. 本文为便于表述，选择以CSR格式为例描述稀疏矩阵到HVMS格式的转换. 如图9所示，经HVMS转换后的稀疏矩阵存储在4个辅助数组中：数组value存储非零元的值；数组col_idx存储非零元的列索引；数组row_idx存储对应非零元的行索引；数组offset存储偏移量，亦即每执行一次向量计算基本操作，从内存加载的元素值在value中的偏移量. 除此之外，HVMS引入数组rest，用于统计非零元数量少于ω 的行的数量. 数组rest中共有ω 个元素，且rest[i] = k表示非零元数量为i的行共有k个.

图  9  HVMS的数据转换过程

Figure  9.  Data conversion process of HVMS

下载: 全尺寸图片 幻灯片

整体来看，HVMS的数据转换过程可分为4个步骤，在此之前需要初始化上述4个数组. 算法2是数据转化过程的伪代码实现，接下来本文以图9为例，分别对每一步骤做一个简单的介绍.

1）步骤1. 遍历稀疏矩阵A，找出所有满足规则1的行. 本文定义nnz_row[i] ≥ ω和ω－nnz_row[i]%ω ≤ δ的所有行执行规约求和操作，其中δ的值可由用户指定. 本文默认δ=ω.

图9中，稀疏矩阵A的第1行满足规则1. HVMS将该行的非零元值、行索引以及非零元值的列索引依次存储至数组value，row_idx，col_idx. 数组offset存储该行首个非零元在数组value中的偏移量. 同时，在遍历稀疏矩阵A的过程中，rest的值也可同步统计出来.

2）步骤2. 利用辅助数组rest，找出所有满足规则2的行. 本文定义nnz_row [i]+nnz_row[j] ≥ θ和nnz_row [i]+nnz_row[j] ≤ ω的所有行执行掩码运算操作，其中θ的值可由用户指定. 本文默认θ = ω. 该步骤需要借助数组rest完成.

图9中，矩阵A的第2行和第 9行的非零元数量分别为3和1，相加之后满足上述定义，因此将上述2行的数据信息顺序存储至4个数组中. 此处需要注意的是，数组offset存储的是2行非零元数量相加之后的偏移量，并非单行非零元数量的偏移量. 最后，数组rest也被更新，例如在本例中，非零元数量为3的行已被标记掩码运算操作，因此rest[3]=rest[3]–1. 同理，第1行和第7行亦符合规则2，因此操作同上.

3）步骤3. 遍历矩阵A的剩余行，找出所有满足规则3的行. 本文定义2点来选择满足要求的行：

1）如果nnz_row [i] == nnz_row [j]==…== nnz_row [l]，并且行数相加等于ω，定义ω行执行前向求和操作；

2）如果行数相加不足ω，那么继续遍历. 当第k行满足

行数加1，直至行数相加等于ω. 最后标记ω行执行前向求和操作. 其中ϵ的值可由用户指定，本文默认ϵ值等于0.

算法2. 数据转换过程的伪代码实现.

输入：稀疏矩阵A（m×n）；

输出：value[]，col_idx[]，row_idx[]，offset[].

① for（i = 0；i < m；i++）

② 　if（nnz_row[i] ≥ ω &（ω－nnz_row[i]%ω）< δ）

③ 　　update(value，col_idx，row_idx，offset)；

④ 　end if

⑤ end for

⑥ for（i = 0；i < m；i++）

⑦ 　if （（nnz_row[i]+nnz_row[j]）≥θ &

⑧ 　　（nnz_row[i]+nnz_row[j]）≤ω）

⑨ 　　update(value，col_idx，row_idx，offset)；

⑩ 　end if

⑪ end for

⑫ for（i = 0；i < m；i++）

⑬ 　if（nnz_row [i] == nnz_row [k]/

$\left|nnz_{\text{row}}\text{[}{i}\text{]}-nnz_{\text{row}}\text{[}{k}\text{]}\right|/nnz_{\text{row}}\text{[}{i}\text{]}$≤ ϵ）

⑭ 　　row_num++；

⑮ 　end if

⑯ 　if ( row_num == ω )

⑰ 　　update(value，col_idx，row_idx，offset)；

⑱ 　　row_num = 0；

⑲ 　end if

⑳ end for

图9中，矩阵A的第3，4，5，6行满足要求. 因此将4行的数据信息分别存储至4个数组，图9中的示例比较特殊，4行是连续的，不过本文对ω行是否连续没有硬性要求. 最后，rest[3]的值也要减4.

对于稀疏矩阵来说，严格满足规则1~3的行在矩阵中的比例并不算大. 本文提出参数化的方式（δ，θ，ϵ）扩充该比例. 通过追加规则4，即零元填充的方式，使略有偏差的行能够满足要求.

4）步骤4. 剩余的所有行将执行CSR计算. 经过步骤1~3的筛选后，矩阵A的非零元数量已经极少，本文在30个数据集上的结果表明数量平均占比只有0.06%，因此对该部分的计算不再进行任何优化.

3.2   HVMS的SpMV实现

本节以图10为例，介绍HVMS的向量计算流程. 此时，输入的稀疏矩阵已经完成规则化转换. 本文设定 ω=4，亦即图10中的向量单元可同时处理4个双精度浮点数据. 从整体上看，基于HVMS的SpMV计算流程分为3个步骤：初始化、向量计算以及数组写回. 这3个步骤结合图10分别进行介绍.

图  10  HVMS的SpMV实现

Figure  10.  SpMV implementation of HVMS

下载: 全尺寸图片 幻灯片

1）初始化阶段. 向量寄存器Reg₁，Reg₂，Reg₃，Reg₄顺序执行规约求和、掩码运算以及前向求和操作. 其中规约求和以及前向求和包含累加计算，因此Reg₁和Reg₄在初始化阶段需被初始化为0.

2）向量计算. HVMS从数组col_idx中取得索引值，并使用gather指令从数组x中读取多个非连续的值；使用load指令将矩阵A中非零元的值加载入向量寄存器；使用fmadd指令实现乘加操作；使用reduction_add指令将向量寄存器内的数据进行累加，其结果即可写回数组y的对应位置；使用mask_reduction_add指令对向量寄存器中对应数据求和，并将结果写回数组y的对应位置. 其中，前4条指令分别对应规约求和、前向求和以及掩码运算操作，第5条指令对应掩码运算操作.

如图10所示，HVMS首先将步骤1中转换出的数据（对应规则1和规则4）载入Reg₁，使用fmadd指令跟数组x中的对应数据做乘加运算，并借助reduction_add指令将结果值相加. 此后，HVMS将步骤2中转换出的数据（对应规则2）载入Reg₂和Reg₃，与数组x中的对应数据做乘法运算后，借助mask_reduction_add指令将对应的结果值相加. 最后，HVMS把步骤3中转换出的数据（对应规则3）载入Reg₄，并使用fmadd指令做乘加运算.

3）数组写回. 前向求和操作可同时处理ω行数据，从而同时写回ω个y值. 在此种情况下，因ω行未必连续，故HVMS会调用scatter指令完成数据的间接内存写回（以row_idx为其索引值）. 使用其他操作计算的数据，如Reg₁，Reg₂，Reg₃，则在相加后使用store或vstore指令分别写回数组y的对应位置.

4.   实　　验

4.1   实验配置

1）实验平台. 基于Intel Xeon（Gold 6130，Skylake）和Intel Core（i9-10940X，Cascade Lake）开展实验对比. 其中Intel Xeon平台包含16个处理器核，最多可开启64个线程，支持AVX512向量指令集，内存为126 GB；Intel Core平台包含14个处理器核，最多可开启28个线程，支持AVX512向量指令集，内存为62 GB. 表2列出本文实验平台的详细参数.

表  2  实验平台的配置信息

Table  2.  Configuration Information of Experimental Platforms

参数	Intel Xeon	Intel Core
处理器	Gold 6130 (Skylake)	i9-10940X (Cascade Lake)
处理器个数	2	1
核心数	16	14
每核线程数	2	2
L1 cache大小/KB	32(I)+32 (D)	32 (I)+32 (D)
L2 cache大小/MB	1	1
L3 cache大小/MB	22	19
内存大小/GB	126	62
操作系统	Ubuntu18.04.5	Ubuntu18.04.3
注：I表示指令缓存，D表示数据缓存.

下载: 导出CSV 

| 显示表格

2）数据集. 本文从稀疏矩阵集合^[14]中选取30个有代表性的稀疏矩阵，该矩阵覆盖电路仿真、社交网络图、高性能计算等多个领域，涵盖不同的矩阵规模和非零元分布模式等. 详细信息如表3所示.

表  3  评测数据集

Table  3.  Datasets Used for Evaluation

数据集	维度	非零元个数	平均每行非零元个数
cage15	5.1×106×5.1×106	99×106	19
caidaRouterLevel	192×103×192×103	1.2×106	6
c-big	345×103×345×103	2.3×106	6
circuit5M	5.5×106×5.5×106	59×106	10
citationCiteseer	268×103×268×103	2.3×106	8
coAuthorsDBLP	299×103×299×103	1.9×106	6
com-DBLP	317×103×317×103	2.0×106	6
com-Orkut	3×106×3×106	234×106	76
com-Youtube	1.1×106×1.1×106	5.9×106	5
dblp-2010	326×103×326×103	1.6×106	4
flickr	820×103×820×103	9.8×106	11
wiki-topcats	1×106×1×106	28×106	15
higgs-twitter	456×103×456×103	14×106	32
in-2004	1×106×1×106	16×106	12
language	399×103×399×103	1.2×106	3
loc-Gowalla	196×103×196×103	1.9×106	9
mip1	66×103×66×103	10×106	155
mouse_gene	45×103×45×103	28×106	642
NotreDame_actors	392×103×127×103	1.5×106	3
nxp1	414×103×414×103	2.6×106	6
preferential	100×103×100×103	1×106	9
rail4284	4×103×1×106	11×106	2633
road_central	14×106×14×106	33×106	2
pwtk	217×103×217×103	11×106	52
soc-Pokec	1.6×106×1.6×106	30×106	18
soc-sign-epinions	131×103×131×103	0.8×106	6
sx-superuser	194×103×194×103	0.9×106	4
webbase-1M	1×106×1×106	3.1×106	3
ldoor	952×103×952×103	42×106	44
CurlCurl_2	806×103×806×103	8×106	11

下载: 导出CSV 

| 显示表格

3）实验对比. 本文选取3个现有的SpMV工作进行实验对比：CVR，SELL-C-σ，Intel MKL（部分现有工作存在编译错误或运行时间过长等问题，本文不予对比. 例如：SpV8^[15]运行时间过长；CSR5^[7]无法在AVX512上成功运行等）. 选取的SpMV具体情况为：

① Intel MKL^[6]是由Intel开发的基于x86平台的基础数学库，提供了绝大部分的线性/非线性算子函数，在众多领域有较广泛的应用.

② CVR^[8]是较具代表性的多行同时计算（前向求和）向量优化方法，经常用于SpMV实验对比.

③ SELL-C-σ^[9]采用零元填充进行向量优化. 本文按照其方案建议，选择C=8（C是向量寄存器的位宽），同时调整σ的值，使用最优结果进行对比.

本文所有实验使用双精度浮点类型，HVMS及对比工作均使用“icc -O3 -xCORE-AVX512 -fopenmp”进行编译，程序在运行过程中开启64（28）个线程. 同时，为准确获取各方案的SpMV执行时间，本文迭代运行10000次并使用平均时间.

4.2   SpMV性能结果分析

首先，本文对比Intel MKL，CVR，SELL-C-σ，HVMS在2个平台上的SpMV运行时间. 图11描绘了30个数据集的实验对比结果. 其中，Intel MKL 的时间被标准化为1. 纵坐标表示CVR，SELL-C-σ，HVMS相对于Intel MKL的加速比. 从图11可知，在Intel Xeon平台上，HVMS相比Intel MKL的加速比最高可达5.16倍，在Intel Core平台上则为3.16倍. 同时，在30个数据集上，HVMS在Intel Xeon平台上分别获得了1.60倍（CVR）、1.72倍（SELL-C-σ）和1.93倍（Intel MKL）的平均加速比；在Intel Core平台上分别获得了1.02倍（CVR）、1.42倍（SELL-C-σ）和1.39倍（Intel MKL）的平均加速比.

图  11  4种方案在2个平台上的实验性能对比

Figure  11.  Experimental performance comparison of four schemes on two platforms

下载: 全尺寸图片 幻灯片

由图11可知，HVMS在Intel Xeon平台上的性能结果优于Intel Core平台，该结果与2个服务器核的CPU架构、内存系统配置、算法本身特性等多方面因素相关. 整体来看，HVMS在高性能计算场景下具备更好的适用性，该优势跟SpMV的应用领域重合度较高，例如深度学习、数值模拟等. 并且在通用场景下，HVMS可取得不低于其他方案的性能.

具体分析来看，HVMS在大部分数据集上的性能都是最优的，且在soc-sign-epinions和preferential上表现较为突出. 这2个数据集的特点是：行与行之间的非零元数量相差很大，非零元分布极不规则. HVMS机制首先对稀疏矩阵进行规则化转换，并使用不同的向量操作处理不同特征的行，从而缓解非规则性对向量优化的影响. 而在com-DBLP，in-2004，ldoor，mip1等矩阵上，HVMS的性能表现较为普通. 分析这些矩阵可知，它们的非零元分布相对比较规则，使用单一规则即可获得较好性能.

另一方面，HVMS在对稀疏矩阵进行规则化转换的过程中，也会影响一部分行的计算顺序，从而影响稀疏矩阵A和数组x内元素的访问顺序，并一定程度上影响访存局部性. 虽然在绝大部分稀疏矩阵上，HVMS的规则化转换是有益的. 但是对于部分稀疏矩阵，例如CurlCurl_2 和road_central上，也会因影响访存局部性而使性能有轻微下降.

4.3   预处理开销分析

为提升SpMV的性能，不少研究工作提出比CSR更为复杂的数据存储格式，故而格式转换过程（也被称之为预处理过程）亦备受关注，它对程序的整体性能会产生一定的影响. 本节利用指标${A}_{\mathrm{p}\mathrm{r}\mathrm{e}}$来分析预处理开销所带来的性能损失，其值表示平摊预处理开销所需的SpMV迭代次数. 计算公式为：

其中${T}_{\mathrm{p}\mathrm{r}\mathrm{e}}$表示预处理时间，${T}_{\rm SpMV}^{\rm MKL}$表示Intel MKL的SpMV算子执行时间，${T}_{\rm SpMV}^{\mathrm{n}\mathrm{e}\mathrm{w}}$表示对比方案的SpMV执行时间（本文中的对比方案为SELL-C-σ，CVR，HVMS）. 本文指定Intel MKL的SpMV执行时间作为衡量其他方案的基准数据. 故而，若Intel MKL优于对比方案，此时${A}_{\mathrm{p}\mathrm{r}\mathrm{e}}$的值被置为∞，表示对比方案的执行时间大于Intel MKL的执行时间. ${A}_{\mathrm{p}\mathrm{r}\mathrm{e}}$的值越小，表示平摊预处理开销所需的SpMV迭代次数越少. 因篇幅受限，本文仅列举Intel Xeon平台上，CVR，SELL-C-σ，HVMS在30个数据集上的${A}_{\mathrm{p}\mathrm{r}\mathrm{e}}$值，如表4所示.

表  4  3种方案平摊格式转换开销所需的SpMV迭代次数

Table  4.  SpMV Iterations that Need to Amortize the Format Conversion Overhead on Three Schemes

数据集	CVR	SELL-C-σ	HVMS （本文）
cage15	∞	∞	∞
caidaRouterLevel	∞	161.47	4.16
c-big	1.73	105.80	3.29
circuit5M	2.66	∞	3.56
citationCiteseer	∞	172.75	16.08
coAuthorsDBLP	∞	88.11	4.01
com-DBLP	∞	119.02	6.60
com-Orkut	∞	474.72	7.35
com-Youtube	2.06	96.67	4.09
dblp-2010	∞	133.69	7.50
flickr	2.53	21.63	4.09
wiki-topcats	369.29	∞	78.67
higgs-twitter	51.17	71.15	46.99
in-2004	19.28	∞	56.19
language	8.14	∞	3.73
loc-Gowalla	3.50	11015.26	1.98
mip1	32.85	∞	115.37
mouse_gene	60.95	∞	133.37
NotreDame_actors	∞	89.13	4.62
nxp1	∞	∞	28.12
preferential	1.22	215.17	1.93
rail4284	12.25	∞	24.19
road_central	28.39	∞	∞
pwtk	153.19	186.86	∞
soc-Pokec	∞	84.20	19.49
soc-sign-epinions	1.66	137.77	2.38
sx-superuser	11.00	193.27	2.95
webbase-1M	4.61	98.82	6.49
ldoor	∞	∞	∞
CurlCurl_2	∞	118.43	∞
注：加粗数字表示3种方案中的最少迭代次数，亦即最优方案.

下载: 导出CSV 

| 显示表格

由表4可知，SELL-C-σ的预处理开销较大，需更多的迭代次数平摊开销，如在2个平台上的最大迭代次数分别为11016和315，平均迭代次数分别为715和102. 而HVMS在所有数据集上的平摊迭代次数均在可接受范围之内，如在2个平台上的最大迭代次数分别为134和221，平均迭代次数分别为24和28，CVR亦能获得不错的效果，如在2个平台上的最大迭代次数分别为370和135，平均迭代次数分别为43和17.

4.4   整体性能分析

真实应用场景下，SpMV过程通常会被迭代计算多次，譬如迭代法求解稀疏线性方程组，该过程将会被计算成千上万次以达到结果值收敛的目的. 而预处理过程仅需执行1次，故而为做更准确的性能评测，本文考虑迭代场景下的程序执行，包含1次预处理过程以及多次SpMV计算过程. 本文将Intel MKL的SpMV算子作为基准测试程序，与其相比，本文其他方案的性能加速比计算公式为：

其中n是迭代计算的次数，$n{T}_{\rm SpMV}^{\rm MKL}$表示n次迭代后Intel MKL的SpMV总执行时间，${T}_{\mathrm{p}\mathrm{r}\mathrm{e}}$表示预处理时间，$n{T}_{\rm SpMV}^{\mathrm{n}\mathrm{e}\mathrm{w}}$表示n次迭代后对比方案的SpMV总执行时间.

图12展示了随着迭代次数的增加，相比于Intel MKL，其他3种方案在30个数据集上的平均加速比. 可以看出，在2个平台上HVMS均能达到最优加速比，故而其在真实应用场景下具有较好的适用性.

图  12  更改迭代次数后4种方案相比于Intel MKL的平均加速比

Figure  12.  Average speedups of four shemes compared with Intel MKL after the iteration numbers change

下载: 全尺寸图片 幻灯片

4.5   可扩展性分析

本文对4种SpMV方案在30个稀疏矩阵上的扩展性进行对比分析. 由实验结果可知，随着线程数目的增加，HVMS获得了较好的性能加速比，具备较好的可扩展性. 在64个线程上，较之其他3种方案的平均加速比（Intel MKL为10.93倍，CVR为11.91倍，SELL-C-σ为12.76倍）， HVMS可获得平均20.35倍的性能提升. 在28个线程上，较之其他3种方案的平均加速比（Intel MKL为5.17倍，CVR为6.92倍，SELL-C-σ为5.35倍），HVMS可获得平均7.43倍的性能提升. 此外，因篇幅受限，本文选取4个有代表性的稀疏矩阵，图13刻画了随着线程数目的增加，4个矩阵在不同优化方案下的性能变化趋势. 由图可知，HVMS的可扩展性优于其他方案，最高可分别获得41倍（Intel Xeon平台）和18倍（Intel Core平台）的性能提升.

图  13  4种方案在4种稀疏矩阵上的可扩展性对比

Figure  13.  Scaling comparison of four schemes on four sparse matrices

下载: 全尺寸图片 幻灯片

4.6   规则优先级对HVMS性能的影响

HVMS中默认的规则优先级为：规约求和→掩码运算→前向求和→零元填充/掩码运算→标量计算. 本文更改规则的优先级，并分析其对性能的影响.

通过对前3种基本操作排列组合，可得到6种不同的优先级组合. 图14展示了在2种平台上、不同的优先级组合下，HVMS相比Intel MKL的性能提升. “增量加速比”指6种组合中最优结果相比默认优先级的性能提升. 其中，几乎所有稀疏矩阵都在不同的规则优先级下获得了不同程度的性能增长. 但大多数矩阵调整规则优先级后的性能增长较小，仅部分矩阵，例如在Intel Xeon平台上，c-big和 preferential获得了较大幅度增长，分别为0.83和1.19. 原因在于这2个稀疏矩阵的行内非零元数量偏少，且具有相同非零元个数的行数较多，故而更适合基于前向求和的操作进行计算. HVMS倾向使用默认的规则优先级来获得较好性能，而非自动调优，以免引入额外的成本影响性能优化效果.

图  14  更改HVMS中规则的优先级对性能的影响

Figure  14.  Effect of changing priorities of rules in HVMS on performace

下载: 全尺寸图片 幻灯片

5.   相关工作

针对SpMV的优化，已有数十年的研究历史. 研究人员致力于在各个平台、各个领域，通过改进现有的存储格式、挖掘硬件特性、设计自动调优框架等等技术手段，从多个维度提升SpMV的性能.

ALBUS^[16]设计了一种简单的逐行向量计算机制，将单行内的非零元依次加载入向量寄存器计算；CSR2^[17]利用零元填充和切块的方式进一步提升向量计算的效率；VHCC^[18]通过拼接多行来提高向量计算的效率，同时对矩阵做2维划分来提升数据局部性； CSR5^[7]通过拼接多行以及分段求和的方法设计向量优化机制，并且还提供了在多个平台上的具体优化实现；CVR^[8]将非零行依次加载入SIMD通道内计算，同时处理SIMD 通道宽度的行，可大大提升向量计算的效率；BCSR^[19]通过切分稠密子块，从而将稀疏矩阵计算转换成整齐的数据计算模式；SELL-C-σ^[9]可看做是对ELL格式的改进，矩阵重排和行划分能够大量减少零元填充的数量；BCCOO^[20]采用分块策略，利用位压缩行索引存储来缓解内存压力；Merge-SpMV^[21]设计了一种更优的负载均衡机制，并且无需数据预处理即可获得不错的性能加速；LAV^[22]针对符合幂律分布的矩阵，设计了行/列划分的机制来提升数据访问局部性；CFS^[23]利用图着色算法缓解对称矩阵的写冲突问题，可大大减轻多线程情况下的内存访问压力；Regu2D^[24]利用矩阵划分、规则重排以及索引压缩的方式提升向量计算的效率；SpV8^[15]聚合规则子矩阵，目的是让更多的数据进行向量计算.

SpMV的优化平台也非常的广泛，文献[25]是在神威平台上的优化工作，文献[26]是在POWER8系统上的优化工作，文献[27]是在FPGA平台上的优化工作，文献[28]是在众核平台上的优化工作，文献[29–31]是在GPU平台上的优化工作.

近几年，由于海量数据的出现，数据类型以及数据分布多种多样，因此自动调优和性能建模的相关工作成为热门的研究内容，代表性工作有文献[32-40].

6.   总　　结

本文深入分析了向量优化SpMV计算的难点以及所带来的挑战. 向量化一直以来都是加速程序性能的关键技术手段之一，然而由于稀疏矩阵数据的不规则性，规则化的向量计算难以为SpMV带来成正比的性能提升. 尤其随着越来越丰富的矩阵数据类型的产生，单个矩阵内的非零元分布以及矩阵之间的非零元分布特征千差万别，并且向量化硬件的多样性也给通用SpMV的向量优化带来了不小的挑战. 基于此，本文抽象了4种向量计算的基本操作：规约求和操作、掩码运算操作、前向求和操作以及零元填充，并且提出了基于混合向量化的SpMV优化机制——HVMS. 该机制根据抽象的向量基本操作，设计相应的规则指导稀疏矩阵的规则化数据转换，目的是为了将不规则的稀疏矩阵运算映射到规则的向量运算上进行计算，最终利用规则的向量指令加速SpMV运算，充分发挥向量计算的效率. 本文以AVX512为例，完成基于HVMS的SpMV优化实现，在30个常用矩阵数据集上进行实验设计与分析，并且相比CVR，SELL-C-σ，Intel MKL，HVMS获得了不错的性能加速比以及更优的线程可扩展性.

作者贡献声明：颜志远提出方案、完成实验并撰写论文；解壁伟提出实验方案并指导论文写作；包云岗提供研究思路、给予指导意见.