当今社会各行各业都在广泛使用计算机和互联网，产生的巨量数据在接入云或互联网后被存入数据中心，其存储的数据可靠性、一致性、安全性和实时性至关重要，其中可靠性是数据存储的关键. 传统的RAID5在k(k>1)个数据块上保存数据，任意(k–1)个数据块都有完整的数据，只允许单个数据块发生故障. 然而随着k的数据量增加，多个数据块发生故障的概率随之增长，RAID5不再适用. 纠删码（erasure code，EC）能够在多个数据块发生故障时恢复出原始数据，提高了系统容错性^[1]. 在允许相同数量级的数据块发生故障的算法中，纠删码具有低存储开销优势^[2]，被广泛应用于各种数据中心的存储系统，包括磁盘阵列系统^[3]、点对点存储系统^[4]、分布式存储系统^[5-6]、云存储系统^[7]等.

虽然纠删码具有高可靠性和低存储开销的优势，但其在编码和解码过程会引发大量密集型计算，这不仅花费过长的运算时间，也对处理器的运算性能和访存带宽带来巨大挑战^[8-9]. 纠删码的运算通常是在有限域内进行，常用的有限域是伽罗华域Galois Field，GF(2^w). 伽罗华域内的运算包括加法和乘法，其中加法运算等同于异或操作，乘法运算等同于多项式运算. 乘法运算效率是伽罗华域运算乃至纠删码运行性能的一个重要的瓶颈点^[10,12]，现有的伽罗华域乘法运算操作主要分2类：基于异或操作和基于查表操作.

基于异或操作方法将伽罗华域的乘法展开为2个或多个异或操作，该方法由Blomer等人^[12]首次提出. Detchart等人^[10]在文献[12]方法基础上使用多项式环变换对异或操作的展开过程进行了优化，减少了展开后的总操作数. Luo等人^[13]分析了异或操作的执行顺序对执行效率的影响，提出数据词引导和数据包引导的调度算法以提升异或操作的缓存命中率，进而提高纠删码的运算性能. Uezato等人^[11]对冗余的异或操作进行合并和压缩，改进了缓存命中率，减少了不必要的访存操作. Zhou等人^[2]对以上各类方法进行了总结，将异或操作的合理调度、缓存命中的改善和数据级并行相结合，旨在最大程度上优化纠删码的运算性能. 基于异或操作的方法能够将复杂的伽罗华域乘法转化为便于计算机实现的简单异或操作，但依旧存在以下关键挑战：

1）优化效率依赖于复杂的冗余合并算法和调度优化算法，优化难度大；

2）异或操作之间存在强数据依赖关系，不利于数据级并行化处理；

3）伽罗华域GF(2^w)中位宽为w位的数据需要展开成维数为w×w的矩阵，增加了数据运算的空间复杂度.

基于查表操作的方法避免了异或操作所面临的挑战，提供了log查表法和完全乘法查表法2种方法来实现伽罗华域的乘法运算. log查表法利用乘法的对数表达将伽罗华域乘法转化为对数查表和指数查表，只需3次查表操作就可等价实现伽罗华域乘法运算. 完全乘法查表法则是将所有可能的伽罗华域乘法运算结果保存到1个数据表中，使用时只需1次查表操作. 此外，基于查表操作的方法可以最大程度地利用数据级并行化来优化有限域乘法运算效率，提升纠删码运算性能^[14-16]. 基于查表操作的方法具有实现简单、运行效率高的优点，但也存在以下挑战.

1）查表过程引入的访存开销大. 查表避免了过渡消耗处理器计算资源，但强依赖于表自身的容量大小和处理器的访存带宽.

2）按字节为单位的细粒度低效数据处理：当有限域的数据位宽小于64 b，即GF(2^w)中的w<64时，在数据存取和数据运算过程对现代处理器的硬件优势发挥不足.

3）利用体系结构相关的单指令多数据（single instruction multiple data，SIMD）指令加速优化设计研究较少.

针对以上挑战，本文基于查表操作的伽罗华域乘法提出了3种优化思路.

1）降低运算和查表开销. 设计4 b分割法，组合采用查表方法和异或操作. 减少了查表空间的同时降低了运算数据之间的依赖程度，为数据并行化提供基础.

2）提升数据处理粒度. 优化了GF(2⁸)查表方法的数据处理粒度，将循环内数据访问和数据运算的粒度从字节提升到长字，避免了重复查表操作.

3）支持数据级并行处理. 利用SIMD指令向量化内存存取、与/异或、移位和查表操作，充分发挥现代处理器算力优势.

本文基于Intel存储加速库（Intel storage acceleration library，ISA-L）^[17]实现了上述优化，并在申威平台和x86平台进行了代码实现. 基于不同数据规模分析和验证了优化效果.

1.   相关背景

1.1   EC算法原理

EC算法的基本原理：对k个同样大小的数据块，添加m个校验块，保证(k+m)个数据中任意丢失n(n≤m)个数据块或校验块时都能恢复到原始数据. 其中，由k个数据块计算得出m个校验块的过程称为数据编码，由(k+m–n)个数据块和校验块恢复丢失数据的过程称为数据解码. EC解码和解码过程如图1. 相较“1份数据存储多个副本”的方式，EC算法在保证高可靠性的同时，减少了数据存储空间. 相较RAID5，EC允许用户按照自身需求设定冗余信息的个数，提高了容错性.

图  1  EC的编码与解码

Figure  1.  Encoding and decoding of erasure codes

下载: 全尺寸图片 幻灯片

EC算法编码过程的运算示意图如图2，k个数据矩阵D左乘编码矩阵M后生成数据和校验矩阵E，表示为

图  2  编码过程的运算示意图

Figure  2.  Matrix representation of encoding

下载: 全尺寸图片 幻灯片

其中M为(k+m)×k维的编码矩阵，D为k×len维的数据矩阵，E为数据块+校验块组成的(k+m)×len维矩阵. 编码矩阵M的前k行组成一个(k×k)的单位阵，其后m行取自范德蒙矩阵或柯西矩阵. 编码矩阵M具备如下特性：任意k行线性不相关，即任意k行组成的方阵是可逆的. 假设有m个数据块/校验块丢失，形成了新的数据+校验矩阵${\hat {\boldsymbol E}}$，将M的对应行去掉后新组成矩阵${\boldsymbol{\hat M}}$，此时编码过程可表示为${\hat {\boldsymbol M}}{\boldsymbol D} = {\hat {\boldsymbol E}}$. 编码矩阵的特性决定了矩阵${\boldsymbol{\hat M}}$是可逆的k阶方阵，原始数据块D可以由${ ( {{\boldsymbol{\hat M}}} )^{ - 1}}{\hat {\boldsymbol E}}$计算得到，丢失的数据得以恢复. 因此，EC算法解码过程可表示为

本质上，EC算法的编码与解码过程均为矩阵相乘，唯一不同是编码使用了编码矩阵，而解码使用了编码矩阵的逆矩阵.

在实际应用中，数据块通常以8位、16位、32位、64位几种位宽存储在计算机中，数据块直接左乘编码矩阵很容易导致校验数据的溢出. 为了将校验数据控制在特定数值范围，EC算法采用了伽罗华域（Galois fields，GF）的代数运算. GF的加法运算（记为$\oplus$）等同于异或操作，GF的乘法运算（记为$\otimes$）等同于与操作. GF(2^w)表示位宽为w位的伽罗华域，数值范围为[0, 2^w–1]，w通常取值为8位、16位、32位、64位等. GF(2^w)的数值可看作是一个多项式，数值的每个对应位数是多项式的系数，比如GF(2⁸)内的数值a=10011101，数据块矩阵元素a的多项式形式可写为

GF(2^w)内数据块矩阵元素对应的多项式a和编码矩阵元素对应的多项式b的加法表示为

GF(2^w)内2个多项式a和b的乘法可理解为2个多项式的乘法，写为

式（4）中的c_j可写为

显然，多项式的乘法结果超出[0, 2^w–1]，这时将乘法结果对本原多项式p(X)取模，即可将计算结果重新限制在预设的数值范围内^[18]. 本原多项式是不能被GF(2^w)中任何多项式整除的质多项式，取模的过程就是逐步削去高阶多项式的过程.

由GF(2^w)的加法和乘法运算规则可知，2个数据的加法运算只需按位执行1次异或操作即可实现. 而乘法运算按照式（4）执行的步骤如下：

1）使用移位和与操作按位提取对应操作数位值. 对于w位的数据提取需要(2w–1)次操作，这样2个操作数共需(4w–2)次提取操作.

2）执行多项式乘法. 计算c_k需(k+1)次$\otimes$操作和k次$\oplus$操作. 这样遍历c₀至c_2(w–1)，共需要(2w–1)²次操作.

3）c₀至c_2(w–1)拼接成一个(2w–1)位的数据，共需(4w–3)次操作.

4）执行(w–1)次，削去$c_wX^w$至$c_{2(w-1)}X^{2(w-1)}$高阶多项式，若c_k为1，$c_kX^k$削去的过程本质上是将$p(X)$移位后与$c(X)$执行异或操作，从$c_wX^w$至$c_{2(w-1)}X^{2(w-1)}$，共(2w–2)次操作.

总结式（3）中的GF(2^w)乘法，共需(4w²+6w–6)次操作. 显然伽罗华乘法的算法复杂度远大于伽罗华加法的复杂度.

1.2   EC函数及应用

EC算法的核心计算就是伽罗华域乘法运算，实现方法有基于异或操作和基于查表操作2种. 基于异或操作被Jerasure^[19]广泛使用，其支持GF(2⁴)至GF(2¹²⁸)之间的广域运算，但基于异或操作计算方法强依赖异或操作的冗余合并算法和调度优化；基于查表操作被ISA-L^[17]库广泛使用，ISA-L是一套用于数据恢复、数据安全、数据压缩和数据加密的开源库. 为了通用性考虑，ISA-L对数据位宽处理较为固定，运算过程全部以字节单位为粒度，即仅支持GF(2⁸)运算，但其在运算吞吐量方面远优于Jerasure^[20].

为了支持实际应用场景中的整体数据写入和部分数据修改，ISA-L为纠删码设计了2类接口：编解码生成函数和编解码更新函数.

编解码生成函数支持整体数据的写入，可将所有数据块一次性集中编码生成校验块，或在数据修复时根据校验块生成数据块. 考虑到编解码生成函数的核心为伽罗华域的矩阵乘，ISA-L为编码和解码设计了ENCODE接口函数. ENCODE的运算规则为$\boldsymbol M \otimes \boldsymbol D = \boldsymbol E$，当2维矩阵M退化为1维矩阵时，运算规则变为点乘函数（记为DOTP）. 当2维矩阵M进一步退化为常数时，运算规则变为基础乘函数（记为MUL）.

编解码更新函数（记为UPDATE）支持对数据块的修改. UPDATE的功能是更新内存数据E. 继续以图2为例，将式（2）可展开为

其中M_j表示编码矩阵M的第j行，D_j表示第j个数据块. 若部分数据块D_j修改为${\hat {\boldsymbol D}_j}$，那么校验块E也应修改为$\hat {\boldsymbol E}$. 此时使用编码公式全部重新计算校验块，计算复杂度等于(k+m)×k维矩阵M和k×len维矩阵D在伽罗华域的乘法计算复杂度. 为了降低计算复杂度，编解码更新函数从式（6）演变为

UPDATE重新计算了部分数据块的修改值（即$\Delta {\boldsymbol D_j}$）所对应的校验块（即${\boldsymbol M_j}\Delta {\boldsymbol D_j}$），并与原校验块相加得到更新后的校验块. 由于伽罗华加的计算复杂度远低于伽罗华乘，UPDATE的计算复杂度仅需考虑${\boldsymbol M_j}$与$\Delta {\boldsymbol D_j}$相乘部分，编码更新函数的计算复杂度等同于维数为(k+m)×1的矩阵M_j和维数为1×len的矩阵$\Delta {\boldsymbol D_j}$在伽罗华域的乘法计算复杂度. 可见，使用编码更新的方法在计算复杂度上远低于全部重新计算校验块的方法.

对式（7）进行抽象，得到UPDATE的运算规则，记为$\left( {M \otimes D} \right) \oplus \boldsymbol E = \hat {\boldsymbol E}$，这里的M是M_j的简写，D是$\Delta {\boldsymbol D_j}$的简写. 具体实现时，E和$\hat {\boldsymbol E}$处于同一段内存空间. 当UPDATE的一维矩阵M退化为常数a时，运算规则退化为乘加函数（记为MAD）.

总结起来，编解码函数和编码更新函数细化为5类：

1）ENCODE表示2维矩阵M与数据矩阵D相乘；

2）DOTP表示1维矩阵M与数据矩阵D相乘；

3）MUL表示2维矩阵M中的常数a与数据矩阵D相乘；

4）UPDATE表示1维矩阵M与数据矩阵D相乘后与校验矩阵E相加；

5）MAD表示1维矩阵M中的常数a与数据矩阵D相乘后与校验矩阵E相加.

5类函数涵盖了纠删码的主要使用场景，属于EC纠删码的核心函数. 表1所列函数的核心算法都是多层循环求解伽罗华乘法与伽罗华加法，在原理上是一致的.

表  1  EC算法的实现函数

Table  1.  The Implementation Functions of Erasure Codes

说明	函数名	运算规则	类别
基础乘	gf_vect_mul	$a \otimes\boldsymbol D = \boldsymbol E$	MUL
点乘	gf_vect_dot_prod	$\boldsymbol M \otimes \boldsymbol D = \boldsymbol E$	DOTP
乘加	gf_vect_mad	$\left( {a \otimes \boldsymbol D} \right) \oplus \boldsymbol E = \hat {\boldsymbol E}$	MAD
编解码	gf_encode_data	$\boldsymbol M \otimes \boldsymbol D = \boldsymbol E$	ENCODE
编码更新	gf_encode_data_update	$\left( {\boldsymbol M \otimes \boldsymbol D} \right) \oplus \boldsymbol E = \hat {\boldsymbol E}$	UPDATE

下载: 导出CSV 

| 显示表格

接下来以编码生成函数ENCODE为例开展优化算法.

编码生成函数的数据块D维数为k×len，编码矩阵M维数为m×k，校验块E维数为m×len. 依据式（1），编码函数的计算过程如算法1所示，属于一个m—len—k的3层循环. 在第2层循环中申请字节类型临时变量temp（第③行），在最内层循环中，循环执行M[l][j]和D[j][i]的GF(2⁸)乘积并累加至temp（第⑤行），最后将temp运算结果赋予E[l][i]（第⑦行），编码计算结束.

算法1. 编码生成函数核心算法.

① for l = 0 : m–1

②　for i = 0 : len–1

③　　temp = 0；

④　　for j = 0 : k–1

⑤　　　temp = temp$\oplus$(M[l][j]$\otimes$D[j][i])；

⑥　　end for

⑦　　E[l][i] = temp；

⑧　end for

⑨ end for

算法1中，GF(2⁸)加法使用异或操作实现，GF(2⁸)乘法涉及到log查表法（logarithm table，LT）和完全乘法查表法（multiplication table，MT）. LT方法将2个字节类型的数据a和b的GF(2⁸)乘法求解转换为查表过程，包括2次对数表查找与1次指数表查找，具体计算过程见表2. LT会提前将1B大小的数据的所有可能的指数结果（表2的gf_exp）与对数结果（表2的gf_log）保存起来，计算乘积时，借助$a \otimes b = 2 \wedge ( {\rm lb} a + {\rm lb} b )$即可完成运算. 考虑到字节数据处理粒度，LT方法的指数表与对数表通常比较小，不会超出处理器的缓存大小，因此其运行效率对处理器缓存大小不敏感，GF(2⁸)乘法计算默认采取该方法. 但LT方法的查表次数较多，查表引入访存开销较大.

表  2  GF(2⁸)乘法的4种计算方法汇总

Table  2.  The Summary of Four Methods for GF(2⁸) Multiplication

方法	伪代码	表大小/B	查表次数	运算次数
LT	c = gf_exp[gf_log[a] + gf_log[b]]	512	3	1
MT	c = gf_mul_table[b×256 + a]	65536	1	2
SH	c = 0;	0	0	24
	c = c$\wedge$((a & 0x1)? t0: 0);
	c = c$\wedge$((a & 0x2)? t1: 0);
	c = c$\wedge$((a & 0x4)? t2: 0);
	c = c$\wedge$((a & 0x8)? t3: 0);
	c = c$\wedge$((a & 0x10)? t4: 0);
	c = c$\wedge$((a & 0x20)? t5: 0);
	c = c$\wedge$((a & 0x40)? t6: 0);
	c = c$\wedge$((a & 0x80)? t7: 0);
SP	low = a & 0x0f;	32	2	4
	high = ($a\gg 4$) & 0x0f;
	c = gf_split[low]$\wedge$gf_split[high+16];
$\wedge$表示异或操作.

下载: 导出CSV 

| 显示表格

为了避免LT方法引入的多次查表问题，MT方法将2个字节类型变量所有可能的相乘结果保存起来，计算乘积时直接从乘法表（表2的gf_mul_table）中查找即可. 相比于LT的3次查表，MT查表次数仅为1次，具体计算过程见表2. MT方法虽然在查表次数上有所减少，但乘法表占用空间较大，通常为64 KB，可能超出L1d缓存的大小（如Intel Xeon E5-2620的L1d Cache为32 KB，SW3231的L1d Cache为32 KB）. 实际运行时，MT方法的查表操作会从L1d缓存访问跨越对L2缓存访问，访问时间成倍增加，导致查表效率下降.

2.   EC算法的优化研究

分析算法1的运算耗时集中在2个部分：1）伽罗华域的乘法运算；2）访问矩阵M，D，E自身开销. 对于伽罗华域的乘法运算，LT方法的对数表和指数表占用空间较小，但查表次数较多；MT方法的查表次数较少，但乘法表占用空间大. 对于M，D，E矩阵的访问，算法1在存取矩阵数据时以字节为单位，属于一种细粒度的数据操作. 但事实上，64位的现代处理器支持以长字为单位的数据粒度操作. 另外，对矩阵D采取按行存储和按列读取的方式，属于离散访存，不利于缓存命中.

这样，针对伽罗华域乘法，在查表时间开销与表存储空间占用的均衡就变得很关键. 假如矩阵访存使用64位内存访问指令一次性读取矩阵D的8列数据，将访存粒度从8 b提升到64 b，则可以极大提升循环内数据处理效率. 此外，本文基于申威指令集提出了SIMD数据并行访问优化，以提升EC函数的运算效率.

2.1   伽罗华域乘法优化

已有的LT和MT伽罗华乘计算方法存在查表次数多和表占用空间大的问题，在性能上依据有提升空间. 首先考虑到的是增加运算量来减少查表操作，本文尝试使用移位分解法（shift，SH），采用移位和异或操作按位进行运算以完成所有字节的乘法运算. 然而以位粒度为单位的运算会显著增加运算开销，额外引入的运算开销会抵消查表减少所带来的收益. 因此，本文结合查表法与移位分解法的优点，提出4 b分割法，将运算单位粒度从1 b提升至4 b. 在减少运算开销的同时控制查表开销.

2.1.1   SH移位分解法

SH方法将多项式以位为单位进行分割，然后分别对a的每个位与多项式b进行运算，最后将所有位的运算结果进行累加.

具体来说，SH方法将2个字节类型的数据a和b的GF(2⁸)乘法运算分解为：1）a的每个位与b做GF(2⁸)乘积；2）所有乘积结果做GF(2⁸)加. SH方法可表达为

其中$a_i^{\rm bit}$表示a中第i个位左移i位. 当${a_i}$为1时，$a_i^{\rm bit} \otimes b$为$\left( {1 \ll i} \right) \otimes b$，记为${t_i}$=$\left( {1 \ll i} \right) \otimes b$，当${a_i}$为0时，$a_i^{\rm bit} \otimes b$为0，因此式（8）可写为

SH方法预先计算出所有${t_i},{\text{ }}i = 0,1,…,7$，并且在运算前从内存中取出，这样可以减少SH方法的运算次数. 对于式（9），由于a按位做分割，SH方法的循环次数为8，每次循环中，共有1次与操作，1次三目运算（也叫条件选择操作），1次异或操作. 完成a与b的乘法计算，共需要执行24次运算，具体运算过程见表2.

SH方法无需查表操作，但运算次数显著增加. 假设每次运算所耗平均时间为N_ca，L1d缓存访问时间为N_L1d，内存访问时间为N_m. 现代处理器中，N_ca通常为0.5~1个时钟周期，N_L1d为3~4个时钟周期，N_m从几十到上百个时钟周期不等^[21-22]. 所有查表操作的平均耗时记为N_av，那么N_av为

其中ratio为缓存命中率，N_av的值介于N_c与N_m之间.

考虑到现代计算机的指令执行多以乱序执行，不存在数据依赖的指令可依赖超标量引擎并行发射^[21-22]，LT方法的执行流程见图3（a），可总结为：1）同时执行2个gf_log查表，耗时N_av；2）执行指令加法，耗时N_ca；3）执行gf_exp查表，耗时N_av. LT方法所耗总时间为(N_av+N_ca+N_av) = (2N_av+N_ca).

图  3  GF(2⁸)乘法的执行流程

Figure  3.  Execution flow of GF(2⁸) multiplication

下载: 全尺寸图片 幻灯片

SH方法的执行流程如图3（b），总结为：1）并行地从内存中读取所有${t_i},{\text{ }}i = 0,1,…,7$，同时执行所有的与操作，耗时N_av；2）同时执行所有的条件选择操作，耗时N_ca；3）逐个执行异或操作共7次，耗时7N_ca. SH方法所耗总时间为(N_av+N_ca+7N_ca) = (N_av+8N_ca).

对比LT方法与SH方法所耗时间. 当N_av<7N_ca时，(2N_av+N_ca)<( N_av+8N_ca)，意味着SH方法耗时比LT方法耗时长. 可见缓存命中率ratio较高时，SH方法运行效率低于LT方法. 显然一味的避免表查找并不能提升伽罗华域乘法的运算效率. 另外，以位为运算单位的方式会导致运算数量激增，额外引入的运算开销会抵消查表减少所带来的收益.

为了进一步减少伽罗华域乘法的运算耗时，本文结合MT方法与SH方法，提出了4 b分割法（4-bit splitting，SP），在减少运算开销的同时控制查表开销.

2.1.2   SP分割法

SP方法在SH方法基础上，将数据块中每个多项式a分割为2部分，同时利用查表法将SH方法的运算转换为查表操作.

具体来说，SP方法将数据块中每个多项式a分割为位宽为4 b的high与位宽为4 b的low这2部分，预先计算high，low与编码矩阵中每一位相乘的所有可能的值存入gf_split表，在计算乘积时只需开展2次查表即可得到a与b的运算结果，具体计算过程见表2.

考虑到high与low均为4 b，那么high，low与字节类型的变量b相乘的结果表均为16 B. 这里将该表称为gf_split表，它的前16 B相乘结果，总大小为32 B. SP方法计算时使用1次移位和2次与操作提取high和low，然后对gf_split查表2次，执行1次异或操作后得到运算结果.

SP方法中，变量b的不同取值得到不同gf_split表，这里将编码矩阵中每个字节看作b，在进行编码操作前，预先计算出所有gf_split表，使用时查表完成. 因此，SP方法在对编码矩阵中不同变量b进行运算时，需频繁读取gf_split表，伴随有内存数据的换入换出. 此外，SP方法会增加存储开销，存储编码矩阵的内存空间增加至32倍：编码矩阵的维数为(k+m)×k，前k行组成单位矩阵，只有后m行保存在内存中，因此SP方法为编码矩阵提供的内存空间大小为(32×k×m)B. 考虑到纠删码的运算性能，k和m的取值一般较小，工程实现中(k, m)通常取值有(3, 2)，(6, 3)，(10, 4)，(17, 3)^[9]，因此SP方法为编码矩阵提供的内存空间大小约为100~2 048 B，存储开销在可接受范围.

SP方法对应实现为算法2，实施步骤如下.

1）将D[j][i]一次性加载到临时变量src中（第⑤行）.

2）将指针P指向M[l][j]所对应的gf_split查找表（第⑥行）.

3）执行移位、与操作（第⑦⑧行）.

4）对P查表，并执行与操作（第⑨行）.

5）异或操作得到伽罗华和（第⑩行）.

算法2. 编码函数的SP方法实现.

① for l = 0 : m–1

②　for i = 0 : len–1

③　　temp = 0；

④　　for j = 0 : k–1

⑤　　　Load D[j][i] to src；

⑥　　　Point P to gf_spit of M[l][j]；

⑦　　　low = src & 0x0f；

⑧　　　high = ($src \gg 4$) & 0x0f；

⑨　　　s = P[low]$\wedge$P[high+16]；

⑩　　　temp$\wedge$= s；

⑪　　end for

⑫　　Store temp to E[l][i]；

⑬　end for

⑭ end for

SP方法的执行流程为见图3（c），运算开销总结为：1）同时计算low和high，计算low耗时N_ca，计算high耗时2N_ca，总耗时记为2N_ca；2）并行执行2个查表操作，耗时记为N_av；3）执行异或操作，耗时记N_ca. SP方法所耗总时间为(2N_ca+N_av+N_ca) = (N_av+3N_ca). 可见，SP方法耗时小于SH方法耗时，考虑到N_av > 2N_ca，SP方法耗时也小于LT方法耗时.

2.2   利用CPU资源优化

SP方法平衡了运算和查表访问开销，但其依旧存在一个问题：继承了GF(2⁸)的细粒度计算，即运算仍然按字节为处理粒度开展. 事实上，各个字节之间的读取过程并不总是存在数据依赖关系，对a和b的数据装入和运算粒度可以提升至字或者长字，这样既利用了现代处理器高位宽的寄存器和指令特点，同时减少了不必要的内存访问.

2.2.1   数据访问粒度扩展

分析算法2，发现有3点不足.

1）数组D在内存中按行存储，相应地，数据存储时也按行存储，但变量D[j][i]的值是按列加载到src中. 引发离散访存，对缓存局部性利用不足.

2）移位、与操作都以字节为单位开展，没有充分发挥现代处理器高位宽的寄存器和指令特点.

3）在最内层循环中，指针P被不断访问，所指向的内存被频繁地从缓存换入换出，引发缓存颠簸.

针对这3点不足，本文在算法2基础上进一步扩增了数据访问粒度，从原来的按字节访问，扩展为按长字操作，形成优化算法3：基于一次内存访问，多次数据利用的思路，减少了访存次数. 具体实施步骤如下.

1）访存粒度从字节提升为长字，一次性读取64 b位宽的数据（D[j][8×i:8×i+7]）到src. 对应的最外层的循环次数减少为原来的1/8.

2）移位、与的运算处理粒度从字节提升为长字.

3）将64 b位宽的数据low对P的查表按字节展开成8个8 b位宽的数据对P的查表；将1个字节型变量temp扩展为8个字节型temp变量（记为temp₀至temp₇）用于保存临时计算结果. 分别对8个temp变量执行查表计算，最后将8个变量结果依次保存到E[l][8×i:8×i+7]中.

算法3. 数据访问粒度扩展优化算法.

① mask = 0x0f0f0f0f0f0f0f0f；

② for l = 0 : m–1

③　for i = 0 : (len–1)/8

④　　temp₀ = 0，temp₁ = 0，temp₂ = 0，temp₃ = 0；

⑤　　temp₄ = 0，temp₅ = 0，temp₆ = 0，temp₇ = 0；

⑥　　for j = 0 : k–1

⑦　　　Load D[j][8×i : 8×i+7] to src；

⑧　　　Point P to gf_spit of M[l][j]；

⑨　　　low = src & mask；

⑩　　　high = ($src \gg 4$) & mask；

⑪　　　temp_0~7$\wedge$= P[($low \gg 4$(8×(0~7))) & 0x0f] 　　　$\wedge$P[(($high \gg 4$(8×(0~7))) & 0x0f)+16]；

⑫　　end for

⑬　　Store temp_0~7 to E[l][8×i:8×i+7]；

⑭　end for

⑮ end for

2.2.2   数据级并行优化

从算法2到算法3的优化，扩展了GF(2⁸)在单循环内的数据访问粒度，提升了内存访问效率，另外变量temp₀至temp₇之间的独立性为数据并行优化提供了条件. 但是算法3还是以标量的代码实现为主，对于支持超长字的处理器，提出用SIMD对数据开展进一步并行化.

以申威SIMD指令为分析对象，对算法2开展数据级并行优化. 申威SIMD支持多种向量数据类型，包括：长度为32 B的整数向量（32×8 b）、长度为16半字的整数向量（16×16 b）、长度为8 b的整数向量（8×32 b）等.

然而要对算法2进行数据并行优化，需要解决2个关键问题：C1. 如何支持位宽为256 b的数据在GF(2⁸)的运算；C2. 如何提升查表效率，减少对数组P的访问次数.

对于C1问题，预先装入表P（即gf_tbl表）的内容到4个宽度为8×32 b的向量寄存器，记为Va₁，Vb₁，Va₂，Vb₂. 这4个向量寄存器的每一个32 b数据都由表P中的一个8 b数据扩展得到. 如图4，Va₁对应表P的第0~7个8 b数据，Vb₁对应表P的第8~15个8 b数据；Va₂对应表P的第16~23个8 b数据，Vb₂对应表P的第24~31个8 b数据. 那么算法2中P[low]是针对Va₁, Vb₁的操作，P[high+16]是针对Va₂, Vb₂的操作.

图  4  Va₁，Vb₁，Va₂，Vb₂的数据结构

Figure  4.  The data structure of Va₁, Vb₁, Va₂, Vb₂

下载: 全尺寸图片 幻灯片

对于C2问题，申威SIMD指令将P一次性载入到向量寄存器中，然后使用字向量混洗（word vector shuffle，VSHFW）指令对向量寄存器中的数据进行重排以替代多字节的查表操作. 然而，将P直接载入到向量寄存器时，会引入一个新的重复载入问题：算法2的3层循环次序为m—len—k，表P指向的内存M[l][j]的循环次序是m—k，那么在算法2的最内层循环里，表P重复len次指向M[l][j]，表P重复len次载入到向量寄存器. 为了避免该问题，我们调整了循环之间的顺序.

结合内存存取、与/异或、移位和查表指令的SIMD实现，申威平台上实现EC算法的SIMD优化版本见算法4.

算法4. 算法2的向量优化版本.

① mask = 0x0f0f0f0f;

② for l = 0 : m–1

③　for j = 0 : k–1

④　　Set Va₁，Vb₁，Va₂，Vb₂ by gf_tbl of M[l][j];

⑤　　for i = 0 : (len–1)/4

⑥　　　Load D[j][4×i : 4×i+3] to src₃₂;

⑦　　　low₃₂ = src₃₂ & mask;

⑧　　　high₃₂ = ($src_{32} \gg 4$) & mask;

⑨　　　Vc = VSHFW (Va₁, Vb₁, low₃₂);

⑩　　　Vd = VSHFW (Va₂, Vb₂, high₃₂);

⑪　　　Vd = Vd$\wedge$Vc;

⑫　　　Extract temp from Vd;

⑬　　　E[l][4×i : 4×i+3] $\wedge$= temp;

⑭　　end for

⑮　end for

⑯ end for

算法4的核心是3层循环m—len—k，最内层循环的计算步骤为：

1）将32 b数据D[j][4×i:4×i+3]加载到src₃₂中（第⑥行）.

2）对src₃₂执行移位、与操作得到low₃₂与high₃₂（第⑦⑧行）.

3）使用查表指令VSHFW得到Vc和Vd（第⑨⑩行）.

4）Vc和Vd异或，结果存入Vd（第⑪行）.

5）将Vd的有效部分抽取为32 b数据temp.

6）temp和E[l][4×i:4×i+3]异或，结果存入E[l][4×i:4×i+3].

3.   实验结果

为了验证本文所提出的SP方法和利用CPU资源优化的加速效果，我们在x86平台和申威平台做了性能加速比测试与分析. 基于Intel ISA-L开源库实现了上述优化，并对其纠删码中最核心的5类函数进行验证.

3.1   实验设计

3.1.1   实验环境

所选择的x86平台和申威平台测试环境配置如表3所示.

表  3  实验环境

Table  3.  Experimental Environment

环境细则	申威平台	x86平台
CPU	SW3231 @ 2.40 GHz 64核	Intel Xeon E5-2620 @ 2.10 GHz 32核
一级缓存/KB	32	32
二级缓存/KB	512	256
三级缓存/KB	65 536	20 480
操作系统	UnionTech OS Server v20 Enterprise	CentOS Linux release 8.1.1911
编译器	GCC 8.3.0	GCC 8.3.1
ISA-L库	V2.30.0	V2.30.0

下载: 导出CSV 

| 显示表格

3.1.2   测试函数

测试函数选用ISA-L中EC纠删码所用到的MUL，DOTP，MAD，ENCODE，UPDATE共5类函数. 选取ISA-L的EC算法中自带的perf测试集来调用上述5种函数功能.

在数据规模选取上（矩阵E和D占用的总空间量可描述为数据规模），考虑到SW3231和E5-2620的实际缓存容量大小，我们通过动态调整矩阵E和D的数据规模，确保数据规模覆盖缓存范围内和缓存范围外的不同情况（范围选取8~256 MB之间）. M的维数为4×10.

3.1.3   实验方法

本文基于吞吐量来评估优化效果. 吞吐量作为衡量纠删码算法的常用标准^[17]，其表示在单位时间内处理的数据规模能力. 为了减少吞吐量计算带来的误差，所有测试数据均采取重复运行10遍取算术平均值得到.

当前开源ISA-L库中EC纠删码运算默认采取LT方法实现，因此我们选取LT方法为对比基准. 不同方法的优化效果用加速比来衡量，加速比的计算方式为不同优化方法的吞吐量与LT方法吞吐量的比值.

3.2   不同伽罗华域乘法运算性能对比

我们从2个方面评估不同伽罗华域乘法性能：1）计算基于SH，MT，LT，SP这4种方法时EC算法执行的吞吐量；2）分析上述4种方法所引入的存储开销.

3.2.1   不同伽罗华域乘法的吞吐量对比

图5和图6为申威平台和x86平台上基于SH，MT，LT，SP这4种方法时EC算法执行的吞吐量对比情况.

图  5  申威平台GF(2⁸)乘法的4种方法吞吐量对比

Figure  5.  Throughput performance comparison of four methods for GF(2⁸) multiplication on the Sunway platform

下载: 全尺寸图片 幻灯片

图  6  x86平台GF(2⁸)乘法的4种方法吞吐量对比

Figure  6.  Throughput performance comparison of four methods for GF(2⁸) multiplication on the x86 platform

下载: 全尺寸图片 幻灯片

在申威平台，MT方法的平均吞吐量是LT方法的1.19倍，SH方法的平均吞吐量是LT方法的0.64倍，SP方法的平均吞吐量是LT方法的2.25倍；在x86平台，MT方法的平均吞吐量是LT方法的1.00倍，SH方法的平均吞吐量是LT方法的0.32倍，SP方法的平均吞吐量是LT方法的1.85倍.

从结果看出，SP方法的吞吐量最高，MT方法和LT方法次之，SH方法最低.

具体看来，申威平台上SH方法在MUL，MAD和UPDATE上的吞吐量分别达到LT方法的86.24%，81.29%，96.26%，但DOTP和ENCODE函数只达到LT方法的28.11%和26.79%. 这是由于MUL，MAD，UPDATE函数采取2层循环运算，而DOTP，ENCODE函数是3层循环运算. SH方法每一次执行最内层循环时，都需要根据不同的编码元素M[l][j]重新从内存中读取计算表（式（10）中的t_i），影响对缓存数据的刷新频率与命中率.

进一步分析SH，MT，LT，SP在运算时对内存操作的影响，统计其缓存失效率. 在申威平台上，4种GF(2⁸)乘法的缓存失效率结果如图7，SH方法执行MUL，MAD，UPDATE函数的缓存失效率是LT方法的0.97倍、1.00倍和0.98倍，执行DOTP和ENCODE函数时缓存失效率则为1.40倍和1.27倍. SH方法在执行DOTP和ENCODE函数时性能差，主要归结于缓存失效率较高.

图  7  申威平台GF(2⁸)乘法的4种方法缓存失效率对比

Figure  7.  Cache miss ratios comparison of four methods for GF(2⁸) multiplication on the Sunway platform

下载: 全尺寸图片 幻灯片

同样地，x86平台上4种GF(2⁸)乘法的缓存失效率测试结果见图8，同样验证了SH方法在执行DOTP和ENCODE函数时存在缓存未命中率较高的现象.

图  8  x86平台GF(2⁸)乘法的4种方法缓存失效率对比

Figure  8.  Cache miss ratios comparison of four methods for GF(2⁸) multiplication on the x86 platform

下载: 全尺寸图片 幻灯片

3.2.2   不同伽罗华域乘法的存储开销对比

本文的存储开销是指纠删码在运算过程中对内存资源的占用，因此以程序执行过程中的最大内存消耗作为量化指标，以评估不同伽罗华域乘法方法的存储开销.

表4和表5为申威平台和x86平台上SH，MT，LT，SP这4种不同计算方法的存储开销对比情况.

表  4  申威平台GF(2⁸)乘法4种方法的存储开销对比

Table  4.  Memory Overhead Comparison of Four Methods for GF(2⁸) Multiplication on the Sunway Platform KB

函数	SH	MT	LT	SP
MUL	136 752	136 824	136 824	136 824
DOTP	44 424	44 424	44 424	44 436
MAD	44 304	44 368	44 304	44 304
ENCODE	43 888	43 928	43 936	44 004
UPDATE	193 932	193 960	193 968	194 024

下载: 导出CSV 

| 显示表格

表  5  x86平台GF(2⁸)乘法4种方法的存储开销对比

Table  5.  Memory Overhead Comparison of Four Methods for GF(2⁸) Multiplication on the x86 Platform KB

函数	SH	MT	LT	SP
MUL	132 876	132 848	132 880	132 880
DOTP	44 384	44 432	44 408	44 436
MAD	44 376	44 376	44 376	44 376
ENCODE	43 948	44 004	43 980	44 032
UPDATE	193 992	19 4016	193 992	19 4016

下载: 导出CSV 

| 显示表格

在申威平台，SH方法的平均存储开销与LT方法的相差–0.036%，MT方法与LT方法相差0.024%，SP方法与LT方法相差0.042%；在x86平台，SH方法的平均存储开销与LT方法的相差–0.026%，MT方法与LT方法相差0.019%，SP方法与LT方法相差0.039%.

从结果看，4种方法的存储开销差异很小，平均存储开销最高的SP方法与LT方法相差不超过0.05%. 综上，SP方法的存储开销略高于LT方法，但在可控范围内. 综合考虑SP方法带来的性能提升，SP方法优于LT方法，适用于纠删码的运算.

3.3   SP方法性能分析

图9为申威平台在不同数据规模情况下，基于算法2的SP方法的性能加速比，其中MUL，DOTP，MAD，ENCODE，UPDATE分别表示5类函数的加速比. 平均值为5类函数的加速比的算术平均值. 整体看来，在不同数据规模下，算法2的加速比在1.55~3.02倍之间，平均加速比为2.10倍. 其中，算法2对于MAD和MUL的加速效果最显著. 分析原因，MAD表示1维矩阵M中的常数a与数据矩阵D相乘后与校验矩阵E相加，MUL表示2维矩阵M中的常数a与数据矩阵D相乘，两者功能存在重合，造成加速比趋同的现象.

图  9  申威平台采用SP方法在不同数据规模的加速比

Figure  9.  Performance speedup with different data scalability by using SP on the Sunway platform

下载: 全尺寸图片 幻灯片

图10为x86平台在不同数据规模情况下，基于算法2的SP方法加速比. 整体看来，在不同数据规模下，算法2的加速比在1.63~2.16倍之间，平均加速比为1.85倍. 其中，算法2对于DOTP和UPDATE的加速效果最显著，同样地，DOTP和UPDATE的功能存在重合，原因与申威平台类似.

图  10  x86平台采用SP方法在不同数据规模的加速比

Figure  10.  Performance speedup with different data scalability by using SP splitting on the x86 platform

下载: 全尺寸图片 幻灯片

3.4   数据访问粒度扩展优化实验

图11为申威平台在不同数据规模情况下，基于算法3的数据访问粒度扩展优化的加速比. 算法3对于DOTP运算加速效果最显著，对于MUL加速效果相对较弱. 整体看来，在不同数据规模下，算法3的加速比在2.63~4.32倍之间，平均加速比为3.28倍. 对比算法2中SP方法的加速效果，算法3提出的数据访问粒度扩展优化能够在算法2的基础上进一步地提升加速比. 有趣的是，随着数据规模继续增大，算法3在申威平台的加速比呈下降趋势，当数据规模达到128MB时，下降幅度有所减缓. 分析原因，当数据规模增加，访问内存数据块变大时，查表操作成为制约运算性能的一个重要因素. 此时运算矩阵数据主要是从内存中读取，缓存失效数量趋于稳定. 总结看来，数据粒度扩展优化在数据规模较小的情况可以充分利用缓存数据，效果更好.

图  11  申威平台采用数据粒度扩展优化在不同数据规模加速比

Figure  11.  Performance speedup with different data scalability by using data granularity amplification on the Sunway platform

下载: 全尺寸图片 幻灯片

图12为x86平台在不同数据规模情况下，基于算法3的数据访问粒度扩展优化的加速比. 算法3对于UPDATE运算加速效果最为显著，对于ENCODE加速效果最弱. 整体看来在不同数据规模下，算法3的加速比在1.93~2.69倍之间，平均加速比为2.36倍. 相比于申威平台，x86平台上加速比随着数据规模的变化幅度趋于平稳，主要原因是随着数据规模增加，申威平台的访存瓶颈凸显，而x86平台上尚未达到访存瓶颈.

图  12  x86平台采用数据粒度扩展优化在不同数据规模加速比

Figure  12.  Performance speedup with different data scalability by using data granularity amplification on the x86 platform

下载: 全尺寸图片 幻灯片

3.5   数据级并行优化实验

在申威平台基于算法4实现的SIMD并行优化加速效果见图13，对于DOTP运算加速效果最显著，对于MUL加速效果相对较弱. 整体看来，在不同数据规模下，算法4的加速比在2.66~4.19倍之间,平均加速比为3.25倍. 对比算法3的加速效果，算法4基于的SIMD并行优化效果整体上与算法3效果相当. 算法4实现利用了申威SIMD指令并行特性，基于256 b位宽的访存指令并行访问，缓解了算法3中的访存带宽瓶颈问题.

图  13  申威平台采用SIMD并行优化在不同数据规模的加速比

Figure  13.  Performance speedup with different data scalability by using SIMD on the Sunway platform

下载: 全尺寸图片 幻灯片

总结图11和图13发现，数据规模小于64MB时，申威平台数据访问粒度扩展优化的加速比大于数据级并行优化的加速比；当数据规模大于64MB时，数据级并行优化加速效果更加显著. 算法3和算法4优化可以针对不用的应用场景的数据规模灵活切换算法. 算法4引入的基于SIMD指令的数据并行优化，其在实现上会依赖于特定的体系结构SIMD指令支持情况，在向其他平台的可扩展性方面存在一定的限制，但是算法3的方法则无此限制.

4.   结　　论

纠删码的编码与解码具有计算和访存密集的特点，主要瓶颈点聚焦于多层循环包裹的伽罗华域乘法运算. 基于查表操作的LT和MT方法实现伽罗华域乘法运算，被业界广泛使用，但是其执行效率受限于查表的构造和数据运算次数，此外以字节为单位的伽罗华域乘法操作在数据并行度利用方面也存在不足.

本文首先针对LT算法的不足，提出了SP方法，将GF(2⁸)的数据块元素对应的多项式a分割为2部分，分别与编码矩阵对应多项式b相乘并存入查询表中，然后在实际GF(2⁸)运算时基于查表和异或操作完成. 此外，在SP方法基础上，利用现代64位处理器数据位宽特性，提出了数据访问粒度扩展优化和数据级并行方法. 本文基于Intel ISA-L开源库实现了上述优化，在申威平台和x86平台开展实验分析. 对比分析了不同优化方法相比于基准LT方法的吞吐量加速比. 实验表明，本文提出的优化方法在吞吐量方面优于LT与MT方法， 在不同的数据规模下，数据访问粒度扩展优化和数据级并行优化均有加速效果. 申威平台，数据粒度扩展优化相比LT平均加速比达3.28倍，SIMD向量化相比LT平均加速比达3.25倍；x86平台，数据粒度扩展优化相比LT平均加速比达2.36倍. 表明了本文提出的优化方法的有效性.

此外，SP方法和数据访问粒度扩展优化的实现方法与架构无关，具备较强的通用性和普适性，可直接应用于其他不同硬件平台. 而基于申威指令集的SIMD并行方案为其他指令集架构的数据级并行优化提供了参考意义.

作者贡献声明：谢汶兵提出了算法思路和设计了论文框架；关睿雪完成论文初稿撰写、完成实验环境部署与实验测试；张艺鸣完成实验数据处理与实验分析；李佳梅完成实验数据整理与实验分析；王俊提出指导意见，参与论文修改.