何 康，黃 春，姜 浩，谷同祥，齊 進(jìn)，劉 杰,3,4

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

1 引言

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

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

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

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

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

2 基本理論

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

2.1 舍入誤差

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

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

aopb=fl(a°b)=

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

(1)

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

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

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

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

(2)

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

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

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

(3)

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

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

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

無(wú)誤差變換的思想如下所示：

設(shè)a,b是2個(gè)浮點(diǎn)數(shù)a,b∈F，且fl(a°b)∈F?？芍獙?duì)于基本的運(yùn)算，浮點(diǎn)數(shù)的誤差仍是一個(gè)浮點(diǎn)數(shù)，所以可以得到：

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

(4)

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

(5)

使用補(bǔ)償?shù)姆椒▽?duì)計(jì)算的結(jié)果進(jìn)行改進(jìn)，即使用一個(gè)巧妙設(shè)計(jì)的修正項(xiàng)來(lái)改善結(jié)果，這就是從浮點(diǎn)數(shù)(a,b)到浮點(diǎn)數(shù)(x,y)的無(wú)誤差變換。

2.3 基于無(wú)誤差變換的補(bǔ)償算法

算法1[8]FastTwoSum

輸入：a,b。

輸出：x,y。

步驟1x=a+b；

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

FastTwoSum是由Dekker[8]于1971年提出的，算法需要滿足|a|≥|b|的條件，共計(jì)3個(gè)浮點(diǎn)運(yùn)算量。

算法2[9]TwoSum

輸入：a,b。

輸出：x,y。

步驟1x=a+b；

步驟2z=x-a；

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

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

算法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]提出的，該算法首先通過(guò)Split算法將輸入的參數(shù)分成2部分再進(jìn)行計(jì)算，需要17個(gè)浮點(diǎn)計(jì)算量。

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

接下來(lái)介紹double-double格式數(shù)的乘法算法prod_dd_dd。與算法add_dd_dd類似，prod_dd_dd的輸入也為2個(gè)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ù)

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

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

MPI_REDUCE是MPI中的歸約操作，對(duì)通信子(communicator)內(nèi)所有進(jìn)程上的數(shù)據(jù)進(jìn)行歸約操作(比如求和、求極大值和邏輯與等)，這個(gè)歸約操作即可以是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)每個(gè)進(jìn)程輸入緩沖區(qū)中的數(shù)據(jù)按op操作組合起來(lái)，并將其結(jié)果返回到序列號(hào)為root的進(jìn)程的輸出緩沖區(qū)中。輸入緩沖區(qū)由參數(shù)sendbuf、count和datatype定義，輸出緩沖區(qū)由參數(shù)recvbuf、count和datatype定義。兩者的元素?cái)?shù)目和類型都相同。所有組成員都用同樣的參數(shù)count、datatype、op、root和comm來(lái)調(diào)用此例程,因此所有進(jìn)程都提供長(zhǎng)度相同、元素類型相同的輸入和輸出緩沖區(qū)。每個(gè)進(jìn)程可能提供一個(gè)元素或一系列元素，組合操作針對(duì)每個(gè)元素進(jìn)行。

3.2 MPI用戶自定義操作MPI_Op_create

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

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

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

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

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

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

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

3.3.1 高精度求和運(yùn)算MPI_DDSUM

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

MPI_DDSUM的流程圖如圖1所示。

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

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

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

MPI_DDSUM算法的核心實(shí)現(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()的核心實(shí)現(xiàn)為：

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

{

temp=add_dd_dd(*inout,*in);

*inout=temp;

in++;

inout++;

}

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

接下來(lái)對(duì)統(tǒng)一歸一化處理的double-double數(shù)據(jù)加法算法CompDDsum進(jìn)行介紹，算法的輸入是一組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ǔ)上實(shí)現(xiàn)了更加高效的MPI_CompDDsum操作。比起原始的MPI_DDSUM操作，MPI_CompDDsum在計(jì)算的最后統(tǒng)一進(jìn)行歸一化處理，降低了計(jì)算成本的同時(shí)，幾乎沒(méi)有降低計(jì)算精度。

其中自定義函數(shù)CompDDsum的核心實(shí)現(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 高精度求積運(yùn)算MPI_DDPROD

本小節(jié)在雙精度乘法算法prod_dd_dd的基礎(chǔ)上實(shí)現(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。

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

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

基于算法prod_dd_dd提出了計(jì)算一組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 個(gè)浮點(diǎn)計(jì)算量。

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

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

高精度的MPI_DDPROD運(yùn)算具有廣泛的應(yīng)用，可用來(lái)計(jì)算三角形矩陣的行列式和求浮點(diǎn)數(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)。

接下來(lái)介紹帶有補(bǔ)償方案的高精度的求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)。

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

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

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

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

4.1 MPI_DDSUM算法的數(shù)值測(cè)試實(shí)驗(yàn)

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

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

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

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

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

4.2 MPI_NORM算法的數(shù)值測(cè)試實(shí)驗(yàn)

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

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

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

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

Figure 3 Relative error comparison when testing with extremely small vectors圖3 使用值極小的向量進(jìn)行測(cè)試時(shí)的相對(duì)誤差比較

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

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

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

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

4.3 性能比較分析

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

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

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

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

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

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

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

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

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

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

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