1 引言

目前, 伴随着高性能计算机的快速发展, 大规模、长时程的科学计算已广泛应用于各个领域之中. 线性代数计算是高性能计算非常重要的领域之一. 通常, 这些计算通过调用基本线性代数子程序库(basic linear algebra subprograms, BLAS)来执行, BLAS库是线性代数应用程序的计算内核. 目前主流的 BLAS 库主要包括两大类: 一类是针对特定处理器专门优化的BLAS库, 例如, Intel 的 MKL数学核心库和 AMD 的 ACML库等, 一类是开源BLAS库, 不同平台都可以对开源库进行优化, 例如, GotoBLAS^[1]和 OpenBLAS^{[2, 3]}等.

现在, 越来越多的应用开始部署在以申威处理器为代表的国产高性能计算平台上, 包括科学计算、气象预报、武器装备、金融统计等领域, 均极度依赖底层BLAS接口. 现有的BLAS库函数一般支持普适性强的单精度和双精度两种浮点计算精度, 由于浮点数表示及浮点数计算不可避免地引入舍入误差, 使得函数在实现过程中会造成计算精度的损失. 双精度函数有时不能满足某些应用的精度要求, 因为浮点计算存在误差导致的事故比比皆是, 例如, 在金融交易领域, 温哥华证券交易所推出的一项初始值为1 000的股票指数, 由于舍入误差累积, 当22个月后发现该问题时, 指数为524.811, 而实际的指数是1 098.892. 在武器装备领域, 美军部署的“爱国者”导弹, 由于软件计时系统存在累计误差, 导致拦截伊军“飞毛腿”导弹失败, 造成了28名士兵死亡^[4]. 一系列事件都表明了误差的危害性^[5], 基于这一现实状况, 设计实现高精度、高可靠性的浮点数值算法很有必要.

点积函数是BLAS1中具有代表性的子程序, 是线性代数应用中最为基础的运算之一, 同时, 它也可以应用于更高级别的运算, 例如矩阵乘法. 精确求点积算法在许多领域中有着不同的应用, 具体应用描述见文献[6, 7]. 目前科研工作者对包括点积函数在内的高精度算法有较多的研究成果. 文献[7]基于double-double算法, 提出了扩展混合精度的XBLAS库. 文献[8]基于无误差变换技术, 设计实现了高精度的浮点数求和及点积算法, 该算法相比XBLAS库中点积和求和算法性能有明显的提升. 文献[9]在GPU上实现了BLAS函数库的4倍精度版本. 文献[10]在GPU上以SUM、DOT、SCAL和AXPY函数为例实现了BLAS一级函数的多精度版本. 文献[11]面向飞腾处理器实现并优化了高精度的求和与点积算法. 文献[12]基于无误差变换技术, 实现并优化了QGEMM算法. 上述高精度点积算法, 并没有结合申威平台有特殊优化. 文献[13-15]基于国产申威平台, 在BLAS库函数性能优化方面做了大量工作.

本文面向国产申威1621平台, 在现有BLAS库的基础上实现了高精度点积函数. 同时, 针对该算法实现过程中的访存延迟和数据依赖问题, 运用访存和指令重排等优化策略, 充分利用国产申威平台硬件特性对其手工优化. 在满足精度提升的基础上, 使性能和精度达到尽可能的平衡.

2 相关理论 2.1 IEEE浮点数标准

IEEE-754是二进制浮点数算术标准, 其定义了如表1所示的几类二进制浮点数格式, 该标准被大部分处理器所采用.

目前, 包括申威处理器在内的绝大多数通用高性能处理器, 在硬件上通常支持普适性较强的单精度和双精度, 并没有对更高精度浮点算术提供相应的硬件支持. 因此, 目前一般采用软件算法来控制浮点运算中的舍入误差, 来模拟实现数值精度的提升. 本文以高精度点积函数为例, 在已有数值算法的基础上, 通过设计误差补偿算法, 从而实现更高精度.

2.2 无误差变换

误差补偿的思想由Dekker^[16]和Kahan^[17]首先进行了相关研究, 2005年, Ogita等人^[8]正式提出了基于误差补偿思想的无误差变换的概念, 其概念如定义1.

定义1. 设 $ \circ \in \left\{ {{{ + , - , }} \times } \right\}$ , a和b为两个浮点数, $a, b \in F$ , 且有 $x = fl\left( {a \circ b} \right) \in F$ , 在没有上下溢出, 且舍入模式是就近舍入时, 存在:

$ \left( {a \circ b} \right) = x + y $

(1)

其中, x表示计算结果的最佳浮点数近似, y表示舍入误差. 把浮点数(a, b)转化为浮点数(x, y)的过程就是无误差变换.

对于两个浮点数加法的无误差变换, 目前已知的最优的是TwoSum算法^[18], 如算法1所示, 该算法在下溢发生时仍然有效.

算法1. TwoSum算法
输入: 浮点数a, b 输出: 浮点数x, y
1) x=a+b2) z=x−a3) y=(a−(x−z))+(b−z)

对于两个浮点数乘法的无误差变换, 最著名的是TwoProd算法^[16]. 由于申威处理器支持混合乘减指令, 可以简化TwoProd算法, 从而有效提升算法的效率, 简化后的算法TwoProductFMA^[8]如算法2所示.

算法2. TwoProductFMA算法
输入: 浮点数a, b 输出: 浮点数x, y
1) x=a×b2) y=FMSx(a, b, x)

这里的FMSx为申威平台上浮点乘减指(a×b−c), 其所实现的功能是: 对于浮点数Fa和Fb相乘后减去Fc, 只对最终的结果进行舍入. 本文选择其作为高精度点积函数基础实现算法之一.

基于上述无误差变换算法的高精度点积算法, 如算法3^[8]所示. 该算法利用TwoProductFMA算法计算出每次浮点数乘积的舍入误差, 利用 TwoSum 算法记录迭代求和过程中的舍入误差. 结合循环求点积运算, 将记录的舍入误差累加并与原点积计算结果求和, 通过误差补偿的方式修正了原点积算法的数值结果, 从而有效提升数值结果的精度. 该算法的数值结果精度与应用双倍工作精度下的dot算法的数值结果精度相同, 就像以四精度计算一样. 文献[8]给出了这一算法的相关理论分析证明, 论证了该算法的可行性.

算法3. AccurateDot算法
输入: 浮点数x, y输出: 浮点数res
1) [p, s0]=TwoProductFMA(x1, y1)2) for i=2:n3)　[h, s1]=TwoProductFMA(xi, yi)4)　[p, s2]=Twosum(p, h)5)　 s0=fl(s+(s0+s1))6) res=fl(p+s0)

3 高精度点积在申威1621上的实现与优化 3.1 高精度点积算法的实现

高精度点积算法实现的主要实现过程如图1所示, 图中的x_i和y_i表示输入的向量元素. 由于访问一次存储器的时间消耗, 远大于访问一次寄存器的时间消耗, 因此, 对于运算的中间误差结果, 这里使用寄存器来进行保存. 每次迭代使用s0寄存器对误差结果进行保存, 最后, 将浮点最佳近似值p与误差和进行求和运算.

该算法只使用与数据相同的工作精度, 不需要对输入数据进行排序, 并且不需要特殊的计算机体系结构, 具备库函数实现的基本条件. 本文以移植到申威平台的Gotoblas库为基础, 在原有单精度和双精度的基础上, 新增高精度的点积函数接口, 用户直接调用新增的高精度函数接口名就可以使用. 对于BLAS库实现框架来说, Kernel核心计算代码处使用汇编代码, 其汇编核心代码如下所示(其中S1–S2指令对应算法1, S3–S8对应算法2), 而其他的框架部分使用C语言来实现.

…

S1: fmuld a0, b0, t0 //对应算法1

S2: fmsd a0, b0, t0, s0

…

$Loop:

…

S3: faddd t0, t1, t2 //对应算法2

S4: fsubd t2, t0, t3

S5: fsubd t2, t3, a5

S6: fsubd t1,t3, t3

S7: fsubd t0, a5, a5

S8: faddd a5, t3, s2

…

S9: faddd s1, s2,s2

S10: faddd s0, s2, s0

S11: fcpys t2, t2, t0

S12: subl I, 1, I

S13: bgt I, $loop

指令说明: fmuld是双精度浮点乘指令, fmsd为双精度浮点乘减指令, 指令功能是浮点寄存器a0与b0相乘再减去t0, 最后结果写入s0, 只对最终结果s0进行舍入. faddd是双精度浮点加指令, fsubd是双精度浮点减指令. fcpys是符号拷贝指令, 其中当源操作数相同时, 可以实现浮点寄存器间的传送功能. bgt是条件转移指令, 当浮点数大于“0”时转移.

3.2 性能优化

性能是评价底层数学库函数的重要衡量标准, 在兼顾高精度、高可靠性的基础上, 需要尽可能地平衡好性能指标. 本文在对高精度点积算法深入分析的基础上, 结合申威1621处理器的结构特点, 从循环、访存、指令重排等多个方面展开性能优化.

循环展开是常用的优化技术之一, 通过循环展开可以减少判断指令的数量和循环变量改变的计算次数, 从而提高处理器流水线的执行性能. 循环展开最重要的就是确定合适的循环展开因子. 对于高精度点积算法而言, 若展开因子过低, 不能充分利用寄存器数量, 若展开因子过大, 则对寄存器压力过大. 申威平台提供了32个浮点寄存器, 其中31个可供自由使用. 本文经过大量实验测试, 以及考虑后续需要对相关指令重排的基础上, 发现当展开因子为4时, 既不会对后续的优化工作产生较大寄存器压力, 也具有良好的优化效果. 故本文设计了循环展开4次的高精度点积核心计算分支. 在循环展开具体实现中, 还使用了指令替换的优化操作, 例如用位与操作来代替模余操作等.

对于包括申威处理器在内的大多数通用处理器来说, 各存储层次访问延迟为主存>缓存>寄存器^[19]. 一般来说, 访问主存有较大的延迟, 为了满足高性能需求, 这里需要对访存优化, 具体可以分为两类解决方案: 一是尽量地减少访存操作次数, 对高精度点积算法在实现过程中, 由于需要对误差数据重复使用, 这里采用将误差数据存储到寄存器中, 来减少访存次数. 此外, 申威1621提供LDx cache控制指令, 支持将内存中的数据预取到cache中, 申威1621处理器的 cache line大小为 128 B, 一个缓存行可以存放16个双精度浮点数据. 对于高精度点积核心数据处, 插入预取指令, 如果所需数据能直接在cache中存取, 则可大幅度减少处理器直接访问主存储器的次数. 二是隐藏访存延迟. 在高精度点积核心数据计算前, 为了避免访存读取造成的流水线停顿, 这里提出如图2所示的数据取用交替策略, 其主要思想是: 在第1次循环计算前, 提前将所需数据加载到寄存器中, 然后, 在做第1次计算的同时, 读取第2次循环所需要的值到寄存器中, 这样, 在经循环展开和数据取用交替策略下, 每个数据 load 和 use 的间隔为4次循环, 完全隐藏了迭代时的读取开销, 避免了“访存一计算”的数据依赖问题, 在 cache 命中的情况下流水线不会有任何延迟.

其次, 高精度点积核心计算过程中, 核心计算指令之间存在较多的数据依赖关系, 数据相关指令往往容易造成流水线的停顿. 这里对核心计算指令进行重排序, 尽可能减少由数据相关造成的流水线停顿. 首先, 对高精度点积算法循环展开前核心计算指令依赖关系分析, 其中步骤S1和S2, S3和S4, S4和S5, S5和S7, S7和S8, S8和S9, S9和S10, S12和S13之间分别存在数据相关. 为了避免当前指令的目的寄存器作为下一条指令的源寄存器造成的流水线中断^[14], 这里不去改变原有指令的先后依赖顺序, 而是在具有数据相关的指令间, 尽可能多地插入无关指令来避免流水线的空转. 经过循环展开后, 在新的循环体中对数据无关指令和数据相关指令进行重排序, 上面采用了数据取用交替策略之后, 访存和计算之间已经没有直接数据依赖关系. 对于访存指令, 这里只需要在基本块内合理地排布, 使访存指令作为数据无关指令穿插在计算指令之中即可, 重排以后的指令段如表2所示.

4 实验结果与分析

申威1621处理器, 其作为国产高性能多核平台, 性能指标优异, 本文以该平台作为实验的测试平台, 该平台具体硬件参数如表3所示.

4.1 精度测试

本文中的数值实验都是IEEE-754标准双精度下进行的, 为了验证本文算法的有效性, 精度测试数据使用文献[8]中的病态浮点数. 该文献给出了高精度点积算法的数值误差理论分析, 其相对误差界如式(2)所示. 这里通过比较在不同病态数据下高精度点积和双精度点积的相对误差, 判断其结果的准确性.

$ \frac{\left|res-{x}^{{{{\rm{T}}}}}y\right|}{\left|{x}^{{\rm{T}}}y\right|}\leqslant eps+\frac{1}{2}\cdot{c}^{2}\cdot cond({x}^{{\rm{T}}}y) $

(2)

其中, eps表示相对误差舍入单位, 对于IEEE-754标准双精度来说eps=2⁻⁵³, c=n∙eps/(1– n∙eps), cond(x^Ty)为点积函数的条件数.

从图3测试结果可以看出, 对于原始双精度点积算法, 随着条件数的增大, 当条件数增大到10¹⁵时, 其相对误差已经大于1, 此时的算法已经失去了准确性. 而本文所实现的高精度点积算法, 在条件数为10¹⁵时, 高精度点积算法的相对误差稳定在10⁻¹⁷–10⁻¹⁵, 较小的相对误差表明高精度点积算法具有更好的准确性, 高精度点积算法直到条件数增大到10³⁰时, 该算法才完全失去准确性. 由此可以得出, 相比双精度点积, 高精度点积提高了计算精度.

4.2 性能测试

性能测试的数据集这里采用特定区间内均匀分布的随机浮点数, 采用的性能测试方法是: 计算被测试函数运行的时间, 这里实验采取运行200次取平均值的形式, 使用加速比来表示性能优化效果, 计算公式如式(3)所示:

$ r = \frac{{Ts}}{{Tp}} $

(3)

其中, Ts表示优化前函数的执行时间, Tp表示经过循环展开、访存优化、指令重排后程序的执行时间.

从图4测试结果可以看出, 当数据规模为10³–10⁷时, 其性能平均加速比为1.61, 最高加速比为2.04. 对于不同的数据规模, 高精度点积算法优化后相比优化前性能均有提升, 这表明优化后的高精度点积算法相比未优化前, 在保证精度提升的基础上, 性能尽可能的平衡. 这里需要说明的是该高精度点积算法与双精度点积相比, 优化后效率低于双精度点积, 因此该算法更适用于对计算精度有更高要求的场景.

5 总结与展望

本文在国产申威1621平台上, 基于开源Gotoblas算法库, 实现并优化了高精度点积函数. 数值结果显示, 文中高精度点积算法的计算结果精度, 近似达到了双精度点积的两倍, 有效控制了浮点计算中舍入误差的累积. 与此同时, 使用一系列与平台架构相关的优化方法对该算法充分优化, 使精度和性能得到了尽可能的平衡. 本文的工作可以为国产自主可控处理器提供软件支持, 但底层数学库的完善与优化是一个长期的系统性工程, 本文没有涉及的BLAS库中的其他函数的高精度实现与优化将是下一步工作的方向.