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

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

1 引言

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

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

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

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

2 理論背景簡(jiǎn)介

2.1 無(wú)誤差變換

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

(a°b)=x+y

(1)

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

2個(gè)浮點(diǎn)數(shù)加法的無(wú)誤差變換[6]如算法1所示。

算法1TwoSum

輸入：a,b。

輸出：x,y。

步驟1x=a+b；

步驟2z=x-a；

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

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

算法2TwoProductFMA

輸入：a,b。

輸出：x,y。

步驟1x=a×b；

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

2.2 基于無(wú)誤差變換的高精度算法

基于無(wú)誤差加法變換的高精度求和AccurateSum算法[8]如算法3所示。其中，步驟1是基于TwoSum算法循環(huán)對(duì)向量a的元素進(jìn)行無(wú)誤差求和變換；步驟2是對(duì)步驟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的相對(duì)誤差界[8]為：

(2)

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

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

計(jì)算向量x,y的點(diǎn)乘時(shí)，AccurateDot的相對(duì)誤差界[8]為：

(3)

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

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

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

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

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

圖1中的ai表示求和輸入的向量元素，ei表示每次求和產(chǎn)生的舍入誤差，使用Sum_i對(duì)求和結(jié)果進(jìn)行記錄，最終將誤差值與Sum_n相加，實(shí)現(xiàn)求和舍入誤差的補(bǔ)償。可以看到，AccurateSum算法以TwoSum算法作為核心，計(jì)算求和的值和舍入誤差值。最終實(shí)現(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軟件實(shí)現(xiàn)高精度算法，通過(guò)新增函數(shù)接口的方式將實(shí)現(xiàn)的算法嵌入到庫(kù)中，包含相應(yīng)頭文件即可直接調(diào)用。對(duì)于C語(yǔ)言和匯編實(shí)現(xiàn)的高精度算法分別添加了新增函數(shù)接口，如表1所示。

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

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

3.2 性能優(yōu)化

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

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

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

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

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

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

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

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

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

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

Table 3 Test of time overhead to prefetch parameters表3 預(yù)取參數(shù)的時(shí)間開(kāi)銷(xiāo)測(cè)試

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

4 數(shù)值實(shí)驗(yàn)與分析

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

Table 4 Configuration of test platform表4 測(cè)試平臺(tái)配置表

4.1 高精度求和算法的可靠性實(shí)驗(yàn)

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

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

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

4.2 高精度點(diǎn)乘算法的可靠性實(shí)驗(yàn)

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

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

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

4.3 性能比較與分析

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

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

Table 5 Comparison table of algorithm and description symbol表5 算法和描述符號(hào)對(duì)照表

Table 6 Analysis of floating-point operations in summation algorithms表6 求和算法浮點(diǎn)運(yùn)算量分析表

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

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

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

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

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

Table 7 Analysis of floating-point operations in dot product algorithms表7 點(diǎn)乘算法運(yùn)算量分析表

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

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

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

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

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

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

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

綜上所述，數(shù)值實(shí)驗(yàn)表明這類(lèi)非計(jì)算密集型算法的瓶頸主要在帶寬，因?yàn)楦呔扔?jì)算增加的額外計(jì)算量對(duì)于算法整體性能影響相對(duì)較小。

5 結(jié)束語(yǔ)

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