基于MPI的高精度歸約函數(shù)設(shè)計與實現(xiàn)*

何 康，黃 春，姜 浩，谷同祥，齊 進，劉 杰,3,4

(1.國防科技大學(xué)計算機學(xué)院,湖南 長沙 410073;2北京應(yīng)用物理與計算數(shù)學(xué)研究所,北京 100000;3.國防科技大學(xué)并行與分布處理重點實驗室,湖南 長沙 410073;4.國防科技大學(xué)復(fù)雜系統(tǒng)軟件工程湖南省重點實驗室,湖南 長沙 410073)

1 引言

近年來，高性能計算HPC(High Performance Computing)在國內(nèi)外取得了高速發(fā)展，在科學(xué)研究、工程技術(shù)和軍事模擬等各個方面有著越來越廣泛的應(yīng)用。

并行計算(Parallel Computing)是以高性能計算機為平臺，應(yīng)用于科學(xué)與工程領(lǐng)域，使用多個中央處理單元或多臺計算機以協(xié)同工作方式解決大規(guī)模運算問題的計算模式[1]。并行計算可以加快計算速度，在更短的時間內(nèi)解決相同的問題或者在相同的時間內(nèi)解決更多的問題。隨著多核處理器和云計算系統(tǒng)的廣泛應(yīng)用，并行已成為有效利用資源的首要手段。目前，國內(nèi)外在高性能計算系統(tǒng)中最廣泛使用的并行編程接口是MPI(Message-Passing Interface)。

MPI是一種基于信息傳遞的并行編程技術(shù)，它定義了一組具有可移植性的編程接口，已成為國際上的一種并行程序標(biāo)準(zhǔn)[1]。MPICH(a high performance portable MPI implementation)是一種最重要的MPI實現(xiàn)。MPICH的開發(fā)與MPI規(guī)范的制定是同步進行的，每當(dāng)MPI推出新版本，就會有相應(yīng)的MPICH的實現(xiàn)版本，所以MPICH最能反映MPI的變化與發(fā)展。MPI_REDUCE是MPI中的歸約操作函數(shù)，該函數(shù)對通信子 (Communicator) 內(nèi)所有進程上的數(shù)據(jù)進行歸約操作，并將計算結(jié)果保存至根進程中， 是在并行計算中經(jīng)常使用的通信函數(shù)。

隨著信息化社會的飛速發(fā)展，人們對于信息處理的要求變得越來越高，計算的大規(guī)模、大尺度、長時程和高維數(shù)的特點變得越來越明顯。浮點計算的舍入誤差的累積效應(yīng)，往往會導(dǎo)致不可信的計算結(jié)果，甚至使最終的結(jié)果失效。設(shè)計高精度的算法，是提高數(shù)值計算結(jié)果準(zhǔn)確性和穩(wěn)定性的有效途徑之一。

基于上述分析，本文基于MPICH提出了一種高精度的歸約函數(shù)MPI_ACCU_REDUCE，采用無誤差變換技術(shù)對數(shù)值計算的舍入誤差進行有效控制。該函數(shù)提供了3種高精度的歸約運算操作，提供更加豐富的計算的同時，能更進一步提高計算結(jié)果的準(zhǔn)確性。

2 基本理論

目前，絕大部分的計算機都支持IEEE-754(1985)[2]標(biāo)準(zhǔn)，該標(biāo)準(zhǔn)定義了二進制32位單精度(single)、64位雙精度(double)2種類型的浮點算術(shù)系統(tǒng)。浮點算術(shù)系統(tǒng)的采用使得舍入誤差不可避免，在這種超大規(guī)模的科學(xué)計算中，由于舍入誤差具有累積性，每次計算產(chǎn)生的極小誤差在累積起來之后，就會使計算結(jié)果失去有效性和準(zhǔn)確性。所以，控制舍入誤差累積，提升數(shù)值算法精度成為了研究的重點。

2.1 舍入誤差

對于如何有效地控制浮點運算中的舍入誤差，最有效的辦法就是提高浮點運算的工作精度。1991年，Goldberg[3]闡述了浮點數(shù)系統(tǒng)中舍入誤差、有效精度等問題對于計算機科研人員的重要性。2008年，IEEE組織考慮到舍入誤差累積的影響，對IEEE-754(1985)標(biāo)準(zhǔn)進行擴展，增加了四精度(quadruple，128 bit)浮點算術(shù)和十進制浮點算術(shù)(decimal arithmetic)等，形成了新的算術(shù)標(biāo)準(zhǔn)，簡稱IEEE-754(2008)[4]，下文簡稱IEEE-754。根據(jù)實現(xiàn)層次的不同，高精度浮點運算的實現(xiàn)可以分為軟件和硬件2個層次[5]。軟件方法主要是從算法層面實現(xiàn)高精度運算，其靈活性要高于硬件方法。

一個標(biāo)準(zhǔn)的浮點計算模型[6]如式(1)所示：

aopb=fl(a°b)=

(a°b)(1+ε1)=(a°b)/(1+ε2),?a∈R

(1)

其中op∈ {加，減，乘，除}，°∈ {+,-,×,÷}，且|ε1|，|ε2|≤u。u為基本算術(shù)運算所使用的機器工作精度，又稱為單位舍入單元(unit round- off)。在IEEE-754浮點標(biāo)準(zhǔn)的單精度中μ近似等于10-8，雙精度中μ近似等于10-16。

該模型給出了浮點數(shù)基本運算的誤差界如式(2)所示：

|a°b-fl(a°b)|≤u|a°b|,

|a°b-fl(a°b)|≤u|fl(a°b)|

(2)

該過程就是由于計算機字長有限而導(dǎo)致計算產(chǎn)生舍入誤差的基本過程。此模型僅在沒有下溢情況時才成立。從模型中可以看出，n個浮點數(shù)的基本運算的向后誤差界限會隨著n的增加不斷增大。

為了進行誤差分析，本文引入2個誤差分析符號θn和γn，設(shè)n為正整數(shù)且nu<1，則有以下結(jié)論：

若εi≤u,ρi=±1,對i=1:n,且nu<1，有:

(3)

其中|θn|≤γn=nu/(1-nu)。

2.2 無誤差變換技術(shù)

無誤差變換技術(shù)(Error-Free Transformation)是設(shè)計補償模式的高精度數(shù)值算法的基本思想。無誤差變換的思想是在二十世紀(jì)六七十年代由Kahan[7]和Dekker[8]提出的。

無誤差變換的思想如下所示：

設(shè)a,b是2個浮點數(shù)a,b∈F，且fl(a°b)∈F?？芍獙τ诨镜倪\算，浮點數(shù)的誤差仍是一個浮點數(shù)，所以可以得到：

x=fl(a±b)?a±b=x+y,y∈F

(4)

x=fl(a·b)?a·b=x+y,y∈F

(5)

使用補償?shù)姆椒▽τ嬎愕慕Y(jié)果進行改進，即使用一個巧妙設(shè)計的修正項來改善結(jié)果，這就是從浮點數(shù)(a,b)到浮點數(shù)(x,y)的無誤差變換。

2.3 基于無誤差變換的補償算法

算法1[8]FastTwoSum

輸入：a,b。

輸出：x,y。

步驟1x=a+b；

步驟2y=b-(x-a)

FastTwoSum是由Dekker[8]于1971年提出的，算法需要滿足|a|≥|b|的條件，共計3個浮點運算量。

算法2[9]TwoSum

輸入：a,b。

輸出：x,y。

步驟1x=a+b；

步驟2z=x-a；

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

TwoSum算法是由Knuth[9]提出的，需要6個浮點運算量。TwoSum不需要先驗條件，且在下溢發(fā)生時仍然有效。

算法3[8]Split

輸入：a。

輸出：x,y。

步驟1c=factor×a；%factor=2s+1

步驟2x=c-(c-a)；

步驟3y=a-x。

算法4[8]TwoProd

輸入：a,b。

輸出：x,y。

步驟1x=a×b；

步驟2[a1,a2]=Split(a)；

步驟3[b1,b2]=Split(b)；

步驟4y=a2×b2-(((x-a1×b1)-a2×b1)-a1×b2)。

TwoProd算法是由Dekker[8]提出的，該算法首先通過Split算法將輸入的參數(shù)分成2部分再進行計算，需要17個浮點計算量。

當(dāng)數(shù)值計算需要近似2倍工作精度時，double-double 數(shù)據(jù)格式是最有效、最常用的選擇。下面介紹double-double數(shù)據(jù)格式的數(shù)值算法,首先介紹double-double格式數(shù)的加法算法add_dd_dd，算法的輸入為2個double-double格式的數(shù)據(jù)a,b,其中ah和bh分別代表a和b的高位，al和bl分別代表a和b的低位，算法輸出為一個double-double格式的數(shù)據(jù)r，rh和rl分別代表r的高位和低位。

算法5[10]add_dd_dd

輸入：a=(ah,al),b=(bh,bl)。

輸出：r=(rh,rl)。

步驟1[sh,sl]=TwoSum(ah,bh)；

步驟2[th,tl]=TwoSum(al,bl)；

步驟3sl=sl+th；

步驟4[th,sl]=FastTwoSum(sh,sl)；

步驟5tl=tl+sl；

步驟6[rh,rl]=FastTwoSum(th,tl)。

接下來介紹double-double格式數(shù)的乘法算法prod_dd_dd。與算法add_dd_dd類似，prod_dd_dd的輸入也為2個double-double格式的數(shù)據(jù)。

算法6[10]prod_dd_dd

輸入：a=(ah,al),b=(bh,bl)。

輸出：r=(rh,rl)。

步驟1[th,tl]=TwoProd(ah,bh)；

步驟2tl=ah×bl+al×bh+tl；

步驟3[rh,rl]=FastTwoSum(th,tl)。

3 高精度歸約函數(shù)

求和和求積運算是科學(xué)工程計算的基礎(chǔ),隨著工程計算的規(guī)模越來越大，提高基本運算的準(zhǔn)確性對于大規(guī)模工程運算具有非常重要的意義。本文以無誤差變換技術(shù)為基礎(chǔ)，提出了高精度的歸約函數(shù)MPI_ACCU_REDUCE，其包括求和、求積和求L2范數(shù)3種高精度歸約運算。

3.1 MPI_REDUCE(歸約操作)工作原理

MPI_REDUCE是MPI中的歸約操作，對通信子(communicator)內(nèi)所有進程上的數(shù)據(jù)進行歸約操作(比如求和、求極大值和邏輯與等)，這個歸約操作即可以是MPI定義的操作,也可以是用戶自定義的操作[12]。

MPI_REDUCE函數(shù)定義為：

intMPI_REDUCE(void*sendbuf,void*recvbuf,intcount,MPI_Datatypedatatype,MPI_Opop,introot,MPI_Commcomm)

函數(shù)接口中的參數(shù)定義如表1所示。

Table 1 Parameter definition of MPI_REDUCE

MPI_REDUCE將組內(nèi)每個進程輸入緩沖區(qū)中的數(shù)據(jù)按op操作組合起來，并將其結(jié)果返回到序列號為root的進程的輸出緩沖區(qū)中。輸入緩沖區(qū)由參數(shù)sendbuf、count和datatype定義，輸出緩沖區(qū)由參數(shù)recvbuf、count和datatype定義。兩者的元素數(shù)目和類型都相同。所有組成員都用同樣的參數(shù)count、datatype、op、root和comm來調(diào)用此例程,因此所有進程都提供長度相同、元素類型相同的輸入和輸出緩沖區(qū)。每個進程可能提供一個元素或一系列元素，組合操作針對每個元素進行。

3.2 MPI用戶自定義操作MPI_Op_create

MPI中已經(jīng)定義好了一些操作，它們?yōu)楹瘮?shù)MPI_REDUCE和其他的相關(guān)函數(shù)提供調(diào)用。這些操作對應(yīng)相應(yīng)的op。例如：MPI_SUM求和操作，MPI_PROD求積操作等。MPI中也提供了一種用戶自定義操作的方式：通過MPI_Op_create()函數(shù)將用戶自定義的操作和自定義的操作符綁定在一起，實現(xiàn)類似的調(diào)用。

MPI_Op_create函數(shù)定義如下：

intMPI_Op_create(MPI_User_function *function,intcommute,MPI_Op *op)

其中，function為用戶自定義的函數(shù),必須具備4個參數(shù):invec、inoutvec、len和datatype。其中invec和inoutvec分別表示將要被歸約的數(shù)據(jù)所在的緩沖區(qū)的首地址,len表示將要歸約的元素個數(shù),datatype
表示歸約對象的數(shù)據(jù)類型。

3.3 高精度歸約函數(shù)MPI_ACCU_REDUCE

雖然MPI中已經(jīng)定義好了一些簡單的操作，然而在大規(guī)模計算中，這些操作運算結(jié)果的精度無法得到有效的保障?；诖耍疚奶岢隽司哂懈呔鹊臍w約函數(shù)MPI_ACCU_REDUCE，其包含求和、求積和求L2范數(shù)3種高精度的歸約運算，提高了歸約計算的精度。

MPI_ACCU_REDUCE函數(shù)定義為:

doubleMPI_ACCU_REDUCE(void *sendbuf,void *recvbuf,intcount,intoptype,introot,MPI_Commcomm)

函數(shù)接口中的參數(shù)定義如表2所示。

Table 2 Parameter definition of MPI_ACCU_REDUCE

用戶在調(diào)用MPI_ACCU_REDUCE進行高精度歸約時，根據(jù)計算需求輸入相應(yīng)的參數(shù)，MPI_ACCU_REDUCE函數(shù)會根據(jù)不同的歸約操作符調(diào)用不同的高精度運算操作，并將計算結(jié)果發(fā)送到根進程的接收消息緩沖區(qū)中。

3.3.1 高精度求和運算MPI_DDSUM

本文在第2節(jié)中介紹了基于無誤差變換技術(shù)實現(xiàn)的double-double格式數(shù)據(jù)的加法算法add_dd_dd。MPI_DDSUM操作便是以算法add_dd_dd為基礎(chǔ)實現(xiàn)的。

MPI_DDSUM的流程圖如圖1所示。

Figure 1 Flow chart of MPI_DDSUM圖1 MPI_DDSUM流程圖

MPI_DDSUM可以實現(xiàn)對一組double-double數(shù)據(jù)的高精度求和，通過算法add_dd_dd實現(xiàn)了自定義函數(shù)ddsum，使用用戶自定義歸約操作函數(shù)MPI_Op_create將ddsum函數(shù)和歸約操作符DDSUM聯(lián)系起來，這樣定義的操作DDSUM可以像MPI預(yù)定義的歸約操作一樣應(yīng)用于各種MPI的歸約函數(shù)中。

MPI_DDSUM同樣可實現(xiàn)一組double數(shù)據(jù)的求和。用戶可以通過MPI_ACCU_REDUCE靜態(tài)庫提供的getDoubleDoubleNum函數(shù)將輸入的double格式的數(shù)據(jù)轉(zhuǎn)換成double-double數(shù)據(jù)。

MPI_DDSUM算法的核心實現(xiàn)如下所示：

…

MPI_Comm_rank(MPI_COMM_WORLD,&rank);

MPI_Comm_size(MPI_COMM_WORLD,&size);

MPI_Type_contiguous(2,MPI_DOUBLE,&ctype);

MPI_Type_commit(&ctype);

MPI_Op_create((MPI_User_function*)ddsum,1,&DDSUM);

MPI_REDUCE(in,inout,count,ctype,DDSUM,root,comm);

MPI_Op_free(&DDSUM);

其中自定義函數(shù)ddsum()的核心實現(xiàn)為：

for(i=0;i< *len;i++)

{

temp=add_dd_dd(*inout,*in);

*inout=temp;

in++;

inout++;

}

自定義函數(shù)ddsum的主體是算法add_dd_dd。該算法每進行一次加法計算都要進行一次歸一化處理，即FastTwoSum操作。歸一化處理的目的是保證double-double數(shù)的高位和低位嚴(yán)格滿足一定的關(guān)系，本文對ddsum函數(shù)進行改進，提出了統(tǒng)一歸一化處理的算法CompDDsum。該算法在最后統(tǒng)一進行歸一化處理，然后補償回原結(jié)果。

接下來對統(tǒng)一歸一化處理的double-double數(shù)據(jù)加法算法CompDDsum進行介紹，算法的輸入是一組double-double格式的數(shù)據(jù)xi(i=1,…,n)，xi.hi和xi.lo分別代表數(shù)據(jù)的高位和低位。

算法7CompDDsum

輸入：一組double-double格式的數(shù)據(jù)xi(i=1,…,n),xi=(xi.hi,xi.lo)。

輸出：res。

步驟1 fori=1:n

步驟2xi+1=TwoSum(xi.hi,xx+1.hi)；

步驟3ri+1=ri+xi.lo+xi+1.lo；

步驟4 end

步驟5temp_res=rn+xn.lo；

步驟6[h,l]=FastTwoSum(xn.hi.temp_res)；

步驟7res=h+l。

在自定義函數(shù)CompDDsum的基礎(chǔ)上實現(xiàn)了更加高效的MPI_CompDDsum操作。比起原始的MPI_DDSUM操作，MPI_CompDDsum在計算的最后統(tǒng)一進行歸一化處理，降低了計算成本的同時，幾乎沒有降低計算精度。

其中自定義函數(shù)CompDDsum的核心實現(xiàn)為：

for(i=0;i< *len;i++)

{

temp=two_sum(inout→hi,in→hi);

r[i] +=inout.lo+in.lo;

inout→hi=temp.hi;

in++;

inout++;

}

r[*len-1] +=inout[*len-1].lo;

inout[*len-1].lo=r[*len-1];

3.3.2 高精度求積運算MPI_DDPROD

本小節(jié)在雙精度乘法算法prod_dd_dd的基礎(chǔ)上實現(xiàn)了高精度求積操作MPI_DDPROD，并比較了普通乘法算法與高精度乘法算法prod_dd_dd的誤差界。

算法8Prod

輸入：一組double格式的數(shù)據(jù)ai(i=1,…,n)。

輸出：res。

步驟1x1=a1；

步驟2 fori=2:n

步驟3xi=xi-1×a1；

步驟4 end

步驟5res=xn。

普通的乘法運算需要n-1個浮點運算量，我們對其誤差界進行分析，其中res代表算法的輸出結(jié)果，a1a2…an為輸入數(shù)據(jù)的精確乘積，eps代表機器精度，該算法誤差界為：

|a1a2…an-res|≤γn-1|res|≤

基于算法prod_dd_dd提出了計算一組double-double數(shù)據(jù)乘積的高精度算法DDProd,算法的輸入是一組double-double格式的數(shù)據(jù) ，ai(i=1,…,n),ai.hi和ai.lo分別代表數(shù)據(jù)的高位和低位。

算法9DDProd

輸入：一組double-double格式的數(shù)據(jù)ai(i,…,n),ai=(ai.hi,ai.lo)。

輸出：res。

步驟1 fori=2:n

步驟2ai+1=prod_dd_dd(ai,ai+1)；

步驟3 end

步驟4res=an.hi+an.lo。

算法DDprod需要 25n-24 個浮點計算量。

假設(shè)在IEEE-754 標(biāo)準(zhǔn)的雙精度格式下，此時機器精度eps=2-53，若輸入數(shù)據(jù)長度n滿足n<249，則可以獲得一個完整準(zhǔn)確的舍入結(jié)果，即算法DDprod會比算法Prod具有更高的精度。

MPI_DDPROD操作通過算法DDprod實現(xiàn)了用戶自定義函數(shù)ddprod，通MPI_Op_create函數(shù)將ddprod函數(shù)和DDPROD操作符聯(lián)系起來，實現(xiàn)了對數(shù)據(jù)的高精度求積操作。

高精度的MPI_DDPROD運算具有廣泛的應(yīng)用，可用來計算三角形矩陣的行列式和求浮點數(shù)的冪等。

3.3.3 高精度求L2范數(shù)操作MPI_NORM

算法10CommonNorm

輸入：一組double格式的數(shù)據(jù)xi(i=1,…,n)。

輸出：res。

步驟1 fori=1:n

步驟2acc=acc+xi*xi；

步驟3 end

步驟4res=sqrt(acc)。

接下來介紹帶有補償方案的高精度的求L2范數(shù)算法ComNorm()。其中S和P均為double-double格式的數(shù)據(jù)，sh和ph分別代表s和p的高位，sl和pl分別代表s和p的低位，最終輸出的結(jié)果res為double格式數(shù)據(jù)。

算法11ComNorm

輸入：一組double格式的數(shù)據(jù)xi(i=1,…,n)。

輸出：res。

步驟1S=[sh,sl]=[0,0]；

步驟2 fori=1:n

步驟3[ph,pl]=TwoProd(xi,xi)；

步驟4[sh,sl]=add_dd_dd(sh,sl,ph,pl);

步驟5 end

步驟6res=sqrt(sh+sl)。

同理，本文通過MPI_Op_create函數(shù)實現(xiàn)了用戶自定義的歸約操作MPI_NORM，實現(xiàn)了高精度的求L2范數(shù)函數(shù)，豐富了MPI的歸約操作。

4 數(shù)值實驗結(jié)果與分析

本文中的所有數(shù)值實驗都是在 IEEE-754(2008)標(biāo)準(zhǔn)雙精度下進行的，計算使用數(shù)據(jù)均為病態(tài)浮點數(shù)。其中3種高精度的歸約操作均在MPICH下使用C語言實現(xiàn)，數(shù)值圖表則是使用Matlab生成的。選用多精度浮點運算庫MPFR作為比較的基準(zhǔn)。

實驗均在Ubuntu 16.04系統(tǒng)中進行，gcc版本為4.7，MPICH的版本為使用MPI-3標(biāo)準(zhǔn)的MPICH 3.3.2。

4.1 MPI_DDSUM算法的數(shù)值測試實驗

在測試MPI_DDSUM時，本文選擇多精度浮點運算庫MPFR中的加法來作為判斷精度是否提升的標(biāo)準(zhǔn)。通過比較MPI_DDSUM和MPI_SUM在不同病態(tài)數(shù)據(jù)量n情況下的相對誤差，判斷計算結(jié)果的準(zhǔn)確性。相對誤差的計算方式為|res-sum|/|sum|，其中res代表算法的輸出結(jié)果，sum為精確的加法和，選取MPFR加法的計算結(jié)果作為精確的加法和sum。

ReproBLAS的求和用例中提供了一種產(chǎn)生正弦波數(shù)據(jù)的方式，生成的數(shù)據(jù)在進行加法運算時具有顯著的病態(tài)性，本文使用正弦波數(shù)據(jù)作為測試數(shù)據(jù)。其數(shù)據(jù)生成方式為：

sin(2 *M_PI* (rank/((double)size)-0.5))

其中，rank為進程號，size為進程總數(shù)，M_PI是C語言中標(biāo)準(zhǔn)庫定義的宏。

由圖2可以看出，MPI_SUM在病態(tài)數(shù)據(jù)量n=103時，其與MPFR加法求和結(jié)果的相對誤差已經(jīng)大于1，即此時MPI_SUM的結(jié)果已經(jīng)失去了準(zhǔn)確性。而隨著病態(tài)數(shù)據(jù)量n的增大，MPI_DDSUM算法的相對誤差穩(wěn)定在10-15～10-10，較小的相對誤差表明MPI_DDSUM的計算結(jié)果具有更好的準(zhǔn)確性。由此可以得出，相比MPI_SUM求和，MPI_DDSUM求和運算提高了計算結(jié)果的準(zhǔn)確性。

Figure 2 Relative error comparison between MPI_SUM and MPI_DDSUM under different n圖2 不同病態(tài)數(shù)據(jù)量n的情況下MPI_SUM與MPI_DDSUM相對誤差對比

4.2 MPI_NORM算法的數(shù)值測試實驗

本小節(jié)選擇高精度的求L2范數(shù)算法MPI_NORM與常規(guī)的求L2范數(shù)算法CommonNorm進行比較，使用多精度浮點運算庫MPFR實現(xiàn)精確的求L2范數(shù)的算法MPFRNorm并作為比較的標(biāo)準(zhǔn)。通過比較在不同病態(tài)數(shù)據(jù)量n下CommonNorm和MPI_NORM的相對誤差，判斷其結(jié)果的準(zhǔn)確性。相對誤差的計算方式為|res-norm|/|norm|，其中，res代表算法的輸出結(jié)果，norm為精確的L2范數(shù)和，本文選取MPFRNorm算法的計算結(jié)果作為精確的范數(shù)和norm。

Graillat等[13]提出了一種生成多種類型隨機浮點數(shù)的方法，其大致思想為針對輸入的指數(shù)值，分別生成了值域上均勻分布的指數(shù)值和有效值，然后根據(jù)這個指數(shù)值和有效值產(chǎn)生浮點數(shù)值。

Graillat等[13]提供的方法可以生成多種不同特點的浮點數(shù)據(jù)，本文選擇范數(shù)逐漸向上溢出的向量和一組值極小的向量這2種類型的數(shù)據(jù)分別進行測試。

先使用一組值極小的向量進行測試，所得結(jié)果如圖3所示。

Figure 3 Relative error comparison when testing with extremely small vectors圖3 使用值極小的向量進行測試時的相對誤差比較

再使用范數(shù)逐漸向上溢出的向量進行測試，此時若求得的相對誤差大于1，則使其等于1，所得結(jié)果如圖4所示。

Figure 4 Relative error comparison when testing with vector for which the norm gradually underflows圖4 使用范數(shù)逐漸上溢的向量進行測試時的相對誤差比較

由圖3可知，當(dāng)使用值極小的一組向量進行測試時，此時MPI_NORM算法的相對誤差小于CommonNorm的相對誤差，且兩者的相對誤差都小于10-12，表明此時2種算法的結(jié)果均具有準(zhǔn)確性。隨著病態(tài)數(shù)據(jù)量n的增大，MPI_NORM和CommonNorm的相對誤差都在增大，由圖3可知，CommonNorm算法相對誤差上升的速度大于MPI_NORM算法的。

如圖4所示，當(dāng)使用范數(shù)值逐漸上溢的向量進行測試時，由于此時必定發(fā)生上溢，數(shù)據(jù)極度病態(tài)，CommonNorm算法的相對誤差始終大于或等于1，表明此時該算法的結(jié)果已經(jīng)失效。而隨著n的增大,MPI_NORM算法的相對誤差緩慢上升，處于10-15～10-10，表明此時MPI_NORM計算的結(jié)果仍保持準(zhǔn)確性。由此可以得出，相比于常規(guī)的CommonNorm算法，MPI_NORM算法提高了計算精度。

4.3 性能比較分析

本小節(jié)將對高精度歸約函數(shù)MPI_ACCU_REDUCE的性能進行測試。在不同進程規(guī)模的情況下，分別測試MPI_ACCU_REDUCE中的加法操作MPI_DDSUM和乘法操作MPI_DDPROD的運行時間，并與MPI_REDUCE的加法和乘法操作的運行時間進行比較。以MPI_REDUCE中的MPI_SUM和MPI_PROD操作的計算時間作為基準(zhǔn)，分別求得加法和乘法計算時間開銷的比值，結(jié)果如圖5所示。

Figure 5 Calculation time ratio of the summation and quadrature algorithms under different process numbers圖5 不同進程數(shù)下的求和和求積算法的計算時間比

由圖5可知，當(dāng)進程數(shù)比較小時，MPI_ACCU_REDUCE中加法操作MPI_DDSUM的計算時間是MPI_REDUCE中的加法操作MPI_SUM計算時間的103～104倍左右，乘法操作MPI_DDPROD的計算時間是MPI_PROD的104～105倍；然而隨著進程數(shù)的增加，加法和乘法時間開銷的比率均逐漸下降，最終穩(wěn)定在10左右。第3節(jié)中對不同算法所需的浮點計算量進行了分析，比起普通的求積和求和操作，高精度的DDSUM和DDPROD操作需要更多的浮點計算量，高精度的求和和求積操作所需的浮點計算量比普通的求和求積操作多10倍左右。算法帶來高精度的同時也降低了計算性能，所以本文算法目前更加適用于一些對精度要求更高的場合，同時精度和速度的差異也是在將來的工作中需要改進的地方。

MPI_ACCU_REDUCE性能較低是由于該函數(shù)中的高精度運算操作需要更多的浮點運算量，同時還需要調(diào)用MPI_Op_create函數(shù)新建操作符和數(shù)據(jù)類型，所以相對于MPI_REDUCE，MPI_ACCU_REDUCE會花費更多的時間。

本文第3節(jié)對MPI_DDSUM的核心實現(xiàn)進行了改進，提出了統(tǒng)一歸一化處理的CompDDsum。以MPI_SUM計算時間為基準(zhǔn)，對核心實現(xiàn)為CompDDsum的MPI_CompDDsum和核心實現(xiàn)為ddsum的MPI_DDSUM進行性能比較測試，分別計算兩者在相同進程數(shù)下計算時間與MPI_SUM計算時間的比率，結(jié)果如圖6所示。

Figure 6 Comparison of calculation time between MPI_DDSUM and MPI_CompDDsum under different process numbers圖6 不同進程數(shù)下MPI_DDSUM與MPI_CompDDsum計算時間對比

由圖6可知，當(dāng)進程數(shù)比較小時，MPI_DDSUM和MPI_CompDDsum的計算時間均是MPI_SUM計算時間的103～104倍左右；隨著進程數(shù)的增加，MPI_DDSUM與MPI_SUM計算時間的比率逐漸穩(wěn)定在10倍左右，而MPI_CompDDsum與MPI_SUM計算時間的比率逐漸穩(wěn)定在7倍左右。所以，只需要改變歸一化處理方式，在最后統(tǒng)一進行歸一化處理，便可以在幾乎不降低精度的情況下，使計算速度有明顯的提高，這同時也表明本文算法還有一定的改進空間。

精度提升的同時，帶來了性能的下降，所以本文高精度歸約操作更適用于一些對計算速度要求較低，而對計算精度有更高要求的場景。這同樣表明，在接下來的工作中，應(yīng)該想辦法對高精度的算法進行優(yōu)化，使其在提高計算精度的同時，性能方面也得到很好的保障。

5 結(jié)束語

隨著科學(xué)計算的大規(guī)模、高維數(shù)、大尺度和長時程的特性變得越來越明顯，高精度的計算方式在未來的并行計算領(lǐng)域變得越來越重要。本文基于無誤差變換技術(shù)的補償算法，改進了MPICH的歸約函數(shù)MPI_REDUCE，實現(xiàn)了高精度的歸約函數(shù)MPI_ACCU_REDUCE，提出了3種高精度的歸約計算操作，包括求和、求積和計算L2范數(shù)。數(shù)值實驗結(jié)果表明，高精度歸約函數(shù)MPI_ACCU_REDUCE有效提高了歸約計算的精度，保證了計算結(jié)果的準(zhǔn)確性。

高精度算法雖然帶來了計算精度的提升，然而需要更多的浮點計算量，這使得算法需要更多的計算成本。這就給我們帶來了一個新的挑戰(zhàn)——如何在計算精度和計算速度之間達(dá)到均衡，在不增加計算成本的情況下實現(xiàn)更加優(yōu)秀的計算精度，這也是未來工作的主要內(nèi)容。