面向飛騰處理器的高精度求和與點乘算法實現(xiàn)和優(yōu)化*

黃 春，姜 浩，谷同祥，齊 進(jìn)，劉文超

(1.國防科技大學(xué)計算機(jī)學(xué)院，湖南 長沙 410073；2.北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京 100088)

1 引言

近些年來，隨著高性能計算機(jī)的不斷發(fā)展，數(shù)值計算對于科學(xué)研究和社會應(yīng)用越來越重要。規(guī)模巨大和運算能力強(qiáng)勁的計算機(jī)系統(tǒng)，被廣泛應(yīng)用于生物醫(yī)療、航空航天、人工智能、信息安全和氣象預(yù)報等領(lǐng)域。高性能計算作為關(guān)鍵領(lǐng)域的核心技術(shù)，已經(jīng)成為代表國家綜合實力的科技要素之一。高性能計算技術(shù)的發(fā)展包括2個主要方面：硬件與軟件。當(dāng)前，計算機(jī)體系結(jié)構(gòu)不斷演進(jìn)，新的計算硬件不斷出現(xiàn)和優(yōu)化。近年來國產(chǎn)處理器迅猛發(fā)展，在一定程度上實現(xiàn)了自主可控。飛騰處理器在自主研發(fā)微內(nèi)核的基礎(chǔ)上，實現(xiàn)了對 ARMv8 架構(gòu)的兼容，其在較低的功耗下獲得了較好的性能。

當(dāng)前軟件的發(fā)展和應(yīng)用相對于國產(chǎn)硬件的發(fā)展還有很多需要提升的方面。高精度數(shù)學(xué)庫的設(shè)計和實現(xiàn)，是解決大規(guī)模數(shù)值計算的可靠性和提升處理器性能的有效途徑之一。線性代數(shù)數(shù)學(xué)庫是高性能計算系統(tǒng)中最為基礎(chǔ)的軟件組件，因此研究線性代數(shù)數(shù)學(xué)庫的高精度實現(xiàn)對發(fā)展高性能計算技術(shù)具有重要意義。經(jīng)過近幾十年的發(fā)展，目前比較廣泛地用于科學(xué)計算的開源BLAS(Basic Linear Algebra Subroutine)庫有GotoBLAS[1]、OpenBLAS[2]和ATLAS[3]等。其中，OpenBLAS是以GotoBLAS為基礎(chǔ)，由中國科學(xué)院計算技術(shù)研究所張云泉和中國科學(xué)院軟件研究所張先秩等人開發(fā)、維護(hù)的。OpenBLAS在多個平臺上生成的代碼性能比較優(yōu)秀，因此本文選用OpenBLAS作為軟件平臺，實現(xiàn)高精度算法。

求和與點乘運算是數(shù)值計算中特別基礎(chǔ)的組成部分，在大規(guī)模的線性代數(shù)計算中被頻繁調(diào)用[4]。在現(xiàn)代科學(xué)計算中，需要處理的數(shù)據(jù)規(guī)模很龐大，而遞歸求和的舍入誤差會隨著數(shù)據(jù)規(guī)模n的增大而逐漸增加[5]，這就使得誤差控制變得尤為重要。2005年，基于Knuth[6]、 Dekker[7]提出的補(bǔ)償加法和乘法，日本早稻田大學(xué)的Oishi、Ogita和德國漢堡大學(xué)的Rump[8]提出了無誤差變換(Error-Free Transformation)的概念，并設(shè)計了高精度的浮點數(shù)求和[8,9]及點乘[8]運算。2008 年，Yamanaka 等[10]設(shè)計了簡單的高精度點乘運算的并行算法。2010年，香港科技大學(xué)的Lu 等[11]在GPU硬件加速平臺上基于QD(Quad-Double library)設(shè)計了GQD(GPU-based Quad-Double library)高精度函數(shù)庫。同年，日本筑波大學(xué)的Mukunoki等[12]在GPU上實現(xiàn)了BLAS 函數(shù)庫4倍精度版本。日本理化研究所的Maho[13]基于QD library[14]及其它高精度軟件庫開發(fā)了高精度計算軟件庫Mpack，其針對BLAS 和LAPACK(Linear Algebra PACKage)實現(xiàn)了高精度的MBLAS(Modular Basic Linear Algebra Subroutine)和MLAPACK(Multiple precision Linear Algebra PACKage)。國防科技大學(xué)的Su等[15 - 17]在數(shù)學(xué)庫的性能和精度方面做了大量工作。

本文利用無誤差變換技術(shù)，結(jié)合飛騰處理器的硬件特點，實現(xiàn)并優(yōu)化了高精度求和算法與點乘算法。應(yīng)用本文設(shè)計的高精度算法，科研工作者可以進(jìn)行更高精度的數(shù)值計算和模擬，對數(shù)值計算的舍入誤差進(jìn)行有效控制，從而使得計算結(jié)果的可靠性和穩(wěn)定性得到進(jìn)一步的保障。

2 理論背景簡介

2.1 無誤差變換

2005 年，Oishi等[8]提出了無誤差變換的概念。設(shè)°∈{+,-,×}，2個浮點數(shù)(a,b)∈F且有x=fl(a°b)∈F，F(xiàn)為浮點數(shù)集合，在沒有上下溢出，舍入模式為就近舍入時，存在:

(a°b)=x+y

(1)

其中，x表示計算結(jié)果最好的浮點數(shù)近似，y∈F表示精確的舍入誤差結(jié)果。滿足上述表達(dá)式的過程即為一對浮點數(shù)加、減和乘法運算的無誤差變換。無誤差變換的基本原理是：在每次浮點操作的過程中，通過給定的計算方法計算出該次浮點操作產(chǎn)生的舍入誤差，然后記錄舍入誤差，最終將累積的舍入誤差補(bǔ)償?shù)礁↑c計算的結(jié)果中，屬于誤差補(bǔ)償算法。

2個浮點數(shù)加法的無誤差變換[6]如算法1所示。

算法1TwoSum

輸入：a,b。

輸出：x,y。

步驟1x=a+b；

步驟2z=x-a；

步驟3y=(a-(x-z))+(b-z)。

對于2個浮點數(shù)乘法的無誤差變換，最為常用的是TwoProduct算法[7]，考慮到TwoProduct算法涉及浮點數(shù)分割操作，導(dǎo)致運算量比較大。并且，由于飛騰處理器的硬件可以直接支持混合乘加指令，結(jié)合硬件特點，本文使用基于混合乘加的乘法無誤差變換TwoProductFMA算法[18]。TwoProductFMA算法如算法2所示。

算法2TwoProductFMA

輸入：a,b。

輸出：x,y。

步驟1x=a×b；

步驟2y=FMA(a,b,-x)。

2.2 基于無誤差變換的高精度算法

基于無誤差加法變換的高精度求和AccurateSum算法[8]如算法3所示。其中，步驟1是基于TwoSum算法循環(huán)對向量a的元素進(jìn)行無誤差求和變換；步驟2是對步驟1求和產(chǎn)生的舍入誤差進(jìn)行累加；步驟3將舍入誤差補(bǔ)償?shù)阶罱K結(jié)果。i和j表示向量a的下標(biāo)索引。

算法3AccurateSum

輸入：a。

輸出：res。

步驟1 fori=2:n

[ai,ai-1]=TwoSum(ai,ai-1)；

步驟2 forj=1:n-1

temp=sum(aj)；

步驟3res=temp+an。

AccurateSum的相對誤差界[8]為：

(2)

其中，eps=2-t，s為計算機(jī)中計算的結(jié)果，本文實驗平臺對應(yīng)的t=64，γn=n·eps/(1-n·eps)，cond項表示求和向量中元素pi對應(yīng)的條件數(shù)。

基于無誤差加法變換和無誤差乘法FMA變換的高精度點乘AccurateDot算法[8]如算法4所示。其中，步驟1使用了向量r存放舍入誤差和計算結(jié)果，其長度2n是向量a長度的2倍；步驟2的循環(huán)中，第1個操作是使用TwoProductFMA算法計算向量元素相乘的結(jié)果和舍入誤差，第2個操作是使用TwoSum算法將相乘得到的結(jié)果進(jìn)行求和，得到求和的計算結(jié)果和舍入誤差，結(jié)合循環(huán)計算，實現(xiàn)向量點乘的計算。步驟3和步驟4將點乘的計算結(jié)果和每次計算過程中產(chǎn)生的舍入誤差相加，通過補(bǔ)償?shù)姆绞綄崿F(xiàn)高精度算法。

算法4AccurateDot

輸入：a,b。

輸出：res。

步驟1[p,r1]=TwoProductFMA(a1,b1)；

步驟2 fori=2:n{

[h,ri]=TwoProductFMA(ai,bi)

[p,rn+i-1]=TwoSum(p,h)}；

步驟3r2n=p；

步驟4res=Sum(r)。

計算向量x,y的點乘時，AccurateDot的相對誤差界[8]為：

(3)

其中，eps=2-t，t=64，γn=n·eps/(1-n·eps)，cond項表示點乘對應(yīng)的條件數(shù)。

3 高精度算法的實現(xiàn)和性能優(yōu)化

3.1 高精度算法的實現(xiàn)

由于高精度點乘和求和算法實現(xiàn)的思想類似，限于篇幅，以高精度求和算法的實現(xiàn)為例進(jìn)行描述。高精度求和算法以加法的無誤差變換為基礎(chǔ)，對求和過程產(chǎn)生的浮點舍入誤差進(jìn)行循環(huán)補(bǔ)償，從而實現(xiàn)精度的提高。AccurateSum算法的主要實現(xiàn)過程如圖1所示。

Figure 1 Realization process of high precision summation圖1 高精度求和實現(xiàn)的過程

圖1中的ai表示求和輸入的向量元素，ei表示每次求和產(chǎn)生的舍入誤差，使用Sum_i對求和結(jié)果進(jìn)行記錄，最終將誤差值與Sum_n相加，實現(xiàn)求和舍入誤差的補(bǔ)償。可以看到，AccurateSum算法以TwoSum算法作為核心，計算求和的值和舍入誤差值。最終實現(xiàn)高精度求和的核心代碼如下所示：

ld {v2.2d,v3.2d,v4.2d.v5.2d,[X]}

add X,X,#64

PRFM FLDL1KEEP,[X,#1024]

fadd v6.2d,v2.2d,v3.2d

fsub v14.2d,v6.2d,v2.2d

fsub v15.2d,v6.2d,v14.2d

fsub v16.2d,v2.2d,v15.2d

fsub v17.2d,v3.2d,v14.2d

fsub v10.2d,v16.2d,v17.2d

本文基于OpenBLAS軟件實現(xiàn)高精度算法，通過新增函數(shù)接口的方式將實現(xiàn)的算法嵌入到庫中，包含相應(yīng)頭文件即可直接調(diào)用。對于C語言和匯編實現(xiàn)的高精度算法分別添加了新增函數(shù)接口，如表1所示。

Table 1 New function interface表1 新增函數(shù)接口

表1中，函數(shù)接口名稱中的“A”代表高精度實現(xiàn)，“C”和“S”分別表示使用C語言和匯編語言實現(xiàn)的算法。

3.2 性能優(yōu)化

本文主要從寄存器分配、預(yù)取操作和循環(huán)優(yōu)化等角度進(jìn)行優(yōu)化，采用手寫匯編的方式實現(xiàn)核心部分，獲得了較高的程序執(zhí)行效率和較好的系統(tǒng)性能。飛騰處理器是基于ARM架構(gòu)設(shè)計的，是load-store型的體系結(jié)構(gòu)處理器,訪問寄存器數(shù)據(jù)的速度要比訪問存儲器數(shù)據(jù)的快得多。為提高效率通常一次性進(jìn)行多個字節(jié)的加載或存儲，但同時也需要考慮寄存器數(shù)量的限制。為了盡可能地有效使用寄存器，對高精度算法進(jìn)行了如下分析：高精度求和的核心為TwoSum算法，高精度點乘的核心為TwoSum算法和TwoProductFMA算法。TwoProductFMA算法和TwoSum算法的有向圖如圖2所示，其中變量a,b,x,y,z均為浮點數(shù)。

Figure 2 Directed graph of TwoProductFMA algorithm and TwoSum algorithm圖2 TwoProductFMA算法和TwoSum算法的有向圖

從圖2可以看到，在不考慮寄存器復(fù)用的情況下，運行1次TwoProductFMA 算法需要4個寄存器，運行1次TwoSum算法需要8個寄存器。設(shè)Fn表示計算核的分塊大小，對于計算核的4,8和16分塊設(shè)計，不考慮寄存器復(fù)用時，所需要的寄存器數(shù)目N可以使用以下公式估計：對于高精度求和，N=Fn/4×8+Fn/4×6+2，而對于高精度點乘，N=2Fn/4×4+2Fn/4×6+3。具體結(jié)果如表2所示。

由表2可以知道，對于高精度求和，在不考慮寄存器復(fù)用的情況下，4分塊設(shè)計會導(dǎo)致每次循環(huán)中有大量寄存器空閑；而16分塊的設(shè)計則明顯超出寄存器限制數(shù)目太多，即使考慮存放中間變量的寄存器進(jìn)行復(fù)用，還是會有較大的寄存器資源競爭。8分塊的設(shè)計，在理論上則較為合適，不會導(dǎo)致寄存器過多閑置，也能盡量避免寄存器資源限制導(dǎo)致的競爭。

Table 2 Number of registers required表2 寄存器需求數(shù)量

對于高精度點乘，在不考慮寄存器復(fù)用的情況下，16分塊的設(shè)計明顯會導(dǎo)致劇烈的寄存器資源競爭；4分塊設(shè)計對寄存器利用不足；而8分塊設(shè)計則超過限制數(shù)目，但考慮到實際計算過程中，串行計算部分的存在，存放中間變量的寄存器可以適當(dāng)復(fù)用。因此，4分塊的設(shè)計和8分塊的設(shè)計都需要考慮，最終以實際測試效果作為選擇依據(jù)。

將計算核做了分塊后，程序設(shè)計的流程圖如圖3所示。

Figure 3 Flow chart of calculation kernel after partitioning圖3 計算核分塊設(shè)計后的流程圖

考慮到每次循環(huán)中做了展開之后，循環(huán)計數(shù)值不是展開次數(shù)的倍數(shù)的情況，設(shè)計了F1計算核，每次取1個字?jǐn)?shù)據(jù)進(jìn)行計算，用于對小于展開次數(shù)的剩余數(shù)據(jù)進(jìn)行計算。F8表示每次取8個字?jǐn)?shù)據(jù)進(jìn)行計算。

為了提高程序的性能，需要有效地命中Cache。對于核心數(shù)據(jù)，插入預(yù)取指令。在高精度算法的匯編實現(xiàn)中，使用PRFM預(yù)取指令對計算核心部分要用到的數(shù)據(jù)進(jìn)行預(yù)取。對不同預(yù)取參數(shù)進(jìn)行測試，得到的結(jié)果如表3所示。

Table 3 Test of time overhead to prefetch parameters表3 預(yù)取參數(shù)的時間開銷測試

從表3可以看到，使用KEEP和STRM的效果沒有明顯差異，這是因為計算過程中取用的數(shù)據(jù)只用到1次，保持或者用完釋放的效果沒有明顯差異。預(yù)取到L1中顯然要比預(yù)取到L2中的效果好。出于一般性考慮，最終預(yù)取參數(shù)為：“PRFM PLDL1KEEP,[X,#1024]”。

4 數(shù)值實驗與分析

本節(jié)所有的數(shù)值實驗都是在 IEEE-754 標(biāo)準(zhǔn)[19]雙精度下進(jìn)行，計算使用的數(shù)據(jù)都是浮點數(shù)的形式。測試平臺的配置如表4所示。

Table 4 Configuration of test platform表4 測試平臺配置表

4.1 高精度求和算法的可靠性實驗

進(jìn)行求和數(shù)值實驗的測試數(shù)據(jù)來自參考文獻(xiàn)[8]，測試得到的結(jié)果如圖4所示。圖4中CommonSum表示普通的遞歸求和測試結(jié)果，AccurateSum表示高精度求和結(jié)果。

Figure 4 Relative error of summation with different condition numbers圖4 不同條件數(shù)下求和的相對誤差

從圖4的測試結(jié)果可以看到，對病態(tài)數(shù)據(jù)求和時，高精度算法AccurateSum的精度比普通求和算法的精度要高。普通遞歸求和在測試數(shù)據(jù)的條件數(shù)增大到1015時，結(jié)果已經(jīng)完全不可靠，而AccurateSum的結(jié)果仍較為良好，直到條件數(shù)增大到1030時才完全失效。測試結(jié)果表明，相較于普通求和算法，本文實現(xiàn)的高精度求和算法，有效地提高了求和計算精度。

4.2 高精度點乘算法的可靠性實驗

在點乘測試的數(shù)值實驗中，計算使用的測試數(shù)據(jù)來源于文獻(xiàn)[8]。圖5中的CommonDot表示普通的點乘結(jié)果，AccurateDot表示高精度點乘結(jié)果。

Figure 5 Relative error of dot product with different condition numbers圖5 不同條件數(shù)下點乘的相對誤差

從圖5的測試結(jié)果可以看到，高精度算法AccurateDot的精度比普通的點乘算法的精度要高。普通的點乘在測試數(shù)據(jù)的條件數(shù)增大到1015時，結(jié)果已經(jīng)完全不可靠。相較于普通的點乘結(jié)果，AccurateDot有更為良好的表現(xiàn)。測試結(jié)果表明，相較于普通的點乘算法，本文實現(xiàn)的高精度點乘算法有效地提高了點乘計算精度。

4.3 性能比較與分析

為了方便描述，使用不同的描述符號表示不同的算法實現(xiàn)方式，算法和描述符號的對照如表5所示。

(1)對普通求和算法和高精度求和算法的浮點運算量進(jìn)行了分析，結(jié)果如表6所示。

Table 5 Comparison table of algorithm and description symbol表5 算法和描述符號對照表

Table 6 Analysis of floating-point operations in summation algorithms表6 求和算法浮點運算量分析表

從表6可以看到，進(jìn)行普通求和的浮點運算量為1，而進(jìn)行一次高精度求和的浮點運算量為7，從浮點運算量來看，高精度求和算法大致為普通求和算法的7倍。

(2)對高精度求和算法的不同實現(xiàn)進(jìn)行了計算速度的測試，在不同數(shù)據(jù)規(guī)模下，以普通求和的匯編實現(xiàn)作為基準(zhǔn)計算時間開銷的比值，結(jié)果如圖6所示。

Figure 6 Calculation time ratio of summation algorithms with different data scales圖6 不同數(shù)據(jù)規(guī)模下求和算法的計算時間比

從圖6可以看出，在不同數(shù)據(jù)規(guī)模下，SUM-A1的計算速度比SUM-A2和SUM-A3都快。進(jìn)行統(tǒng)計平均之后，SUM-A1、SUM-A2和SUM-A3除以SUM-COM的平均時間比分別為1.76，3.32和7.35。這表明在使用3.3節(jié)所述的一系列性能優(yōu)化策略之后，SUM-A1的計算速度得到了有效提升。結(jié)合處理器的硬件特點，對精度提高和計算速度實現(xiàn)了有效平衡。

(3)對普通點乘算法和高精度點乘算法的浮點運算量進(jìn)行了分析，結(jié)果如表7所示。

Table 7 Analysis of floating-point operations in dot product algorithms表7 點乘算法運算量分析表

從表7可以看出，在使用FMA指令之后，進(jìn)行一次普通點乘的浮點運算量為1，而高精度點乘算法的浮點運算量為11，從浮點運算量來看，高精度點乘算法大致為普通點乘算法的11倍。

(4)對高精度點乘算法的不同實現(xiàn)進(jìn)行了計算速度的測試，在不同數(shù)據(jù)規(guī)模下，以普通點乘的匯編實現(xiàn)作為基準(zhǔn)計算時間開銷的比值，結(jié)果如圖7所示。

Figure 7 Calculation time ratio of dot product algorithms with different data scales圖7 不同數(shù)據(jù)規(guī)模下點乘算法的計算時間比

從圖7可以看出，在不同數(shù)據(jù)規(guī)模下，DOT-A1的計算速度優(yōu)于DOT-A2和DOT-A3的。進(jìn)行統(tǒng)計平均之后，DOT-A1、DOT-A2和DOT-A3除以DOT-COM的平均時間比分別為1.57，4.65和7.70。這表明在使用3.3節(jié)所述的性能優(yōu)化策略之后，DOT-A1的計算速度得到了有效提升。在精度提高和計算速度之間實現(xiàn)了有效的平衡。

Figure 8 The change of calculation time of dot product algorithm with the number of threads圖8 點乘算法計算時間隨線程數(shù)目的變化

(5)對點乘的并行版本進(jìn)行了測試，測試結(jié)果如圖8所示。其中DotCommonPF1對應(yīng)的是普通點乘的并行情況，實現(xiàn)內(nèi)核為F1設(shè)計；DotAccuratePF1為高精度點乘并行計算的情況，使用的內(nèi)核為F1設(shè)計；DotAccuratePF8為高精度點乘并行計算的情況，使用的內(nèi)核為F8設(shè)計。

從圖8可以看到，DotAccuratePF1為了提高精度，增加了計算量，總體計算速度比DotCommonPF1都要慢。在進(jìn)行F8內(nèi)核設(shè)計并進(jìn)行優(yōu)化之后，DotAccuratePF8的計算速度有了改善，甚至在某些線程并行情況下，計算速度比DotCommonPF1略快。從并行測試的情況來看，也再次表明優(yōu)化策略是有效的。但是，由于此類帶寬受限的算法在多線程并行實現(xiàn)過程中涉及歸約操作，往往帶來不可避免的Cache 缺失情況，因此多線程的擴(kuò)展性不是很優(yōu)秀。

綜上所述，數(shù)值實驗表明這類非計算密集型算法的瓶頸主要在帶寬，因為高精度計算增加的額外計算量對于算法整體性能影響相對較小。

5 結(jié)束語

隨著科學(xué)工程計算規(guī)模的增加、計算時程的增長，數(shù)值結(jié)果的精度可靠性越來越受到研究者的重視。對于國防安全等需求，如何在國產(chǎn)處理器上實現(xiàn)自主可控的軟件至關(guān)重要。本文基于國產(chǎn)飛騰處理器，面向開源軟件平臺OpenBLAS，實現(xiàn)并優(yōu)化了高精度求和算法與點乘算法。數(shù)值實驗結(jié)果表明，雖然高精度算法增加了一定的浮點運算量，但計算結(jié)果的可靠性得到了進(jìn)一步的保證，有效控制了浮點運算在大規(guī)模計算中舍入誤差的累積。與此同時，結(jié)合國產(chǎn)飛騰處理器上的硬件特點，使用一系列性能優(yōu)化策略，對計算精度和計算速度實現(xiàn)了有效的平衡。本文的工作可以用于提升大規(guī)模科學(xué)工程計算的穩(wěn)定性和可靠性，以及為在國產(chǎn)自主可控處理器上實現(xiàn)自主可控軟件提供了支持。對于算法實現(xiàn)的程序性能，可以在性能平衡上作進(jìn)一步研究，以及以此為基礎(chǔ)，對更上層的計算進(jìn)行研究，為數(shù)學(xué)庫的發(fā)展提供更多支持。