DRM:基于迭代歸并策略的GPU并行SpMV存儲(chǔ)格式*

王宇華,何俊飛,張宇琪,徐悅竹,崔環(huán)宇

(1.哈爾濱工程大學(xué)計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院,黑龍江 哈爾濱 150001; 2.電子政務(wù)建模仿真國家工程實(shí)驗(yàn)室,黑龍江 哈爾濱 150001)

1 引言

在科學(xué)計(jì)算和工程應(yīng)用領(lǐng)域中,稀疏矩陣向量乘法SpMV(Sparse Matrix-Vector multiplication)被廣泛應(yīng)用于流體力學(xué)[1,2]、結(jié)構(gòu)力學(xué)[3]、空氣動(dòng)力學(xué)[4]、分子動(dòng)力學(xué)[5]和天體物理學(xué)[6]等領(lǐng)域。特別地,在利用迭代法求解線性系統(tǒng)的過程中,SpMV往往占據(jù)了高達(dá)60%～70%的計(jì)算量[7]。因此,對(duì)SpMV及其優(yōu)化的研究對(duì)于加速科學(xué)計(jì)算具有重要的意義,這也使其成為高性能計(jì)算領(lǐng)域的研究熱點(diǎn)之一。

大多數(shù)SpMV過程形如:y=Ax。其中A表示稀疏矩陣,x是輸入向量,y是輸出向量。矩陣的稀疏性導(dǎo)致SpMV運(yùn)算效率低下,SpMV的性能高度依賴于矩陣的非零元分布,非零元分布決定了其內(nèi)存訪問模式,并且不同矩陣之間有很大的差異。稀疏對(duì)角矩陣是一類特殊的稀疏矩陣,其非零元素按照對(duì)角線的形式密集排列。針對(duì)該特殊的結(jié)構(gòu),Bell等人[8,9]首次將DIA(DIAgonal)格式應(yīng)用到了GPU平臺(tái)上,并通過實(shí)驗(yàn)證明了DIA更適合于稀疏對(duì)角矩陣。DIA按照矩陣中的對(duì)角線來進(jìn)行存儲(chǔ),一個(gè)對(duì)角線利用一個(gè)偏移量進(jìn)行索引,相比于其他格式,這樣不僅能夠減少列索引的存儲(chǔ),同時(shí)還具有良好的時(shí)間空間局部性和訪問向量的連續(xù)性。但是,當(dāng)矩陣中存在散點(diǎn)和長(zhǎng)斷行時(shí),DIA格式會(huì)出現(xiàn)零填充現(xiàn)象[10],在運(yùn)算時(shí)會(huì)造成空間和計(jì)算資源的浪費(fèi),最終導(dǎo)致空間和時(shí)間的利用率下降。

近年來,很多研究人員提出了不同的改進(jìn)格式,如CRSD(Compressed Row Segment Diagonal)[11,12]、DDD-SPLIT(Direct Dense Diagonal SPLIT)[13]、 HDIA(Hacked DIAgonal)[14]、IHYB(Improved HYB format)[15]和DIA-Adaptive (DIAgonal-Adaptive)[16]等。其中,大多數(shù)格式是通過對(duì)矩陣進(jìn)行分塊來解決零填充問題,雖然使SpMV性能有所提升,但由于稀疏矩陣中非零元分布不規(guī)則,塊間非零元個(gè)數(shù)存在一定差異,會(huì)出現(xiàn)塊間負(fù)載不平衡的問題。

負(fù)載不平衡會(huì)引起多個(gè)線程的空閑等待,浪費(fèi)了計(jì)算資源的同時(shí)降低了運(yùn)算的并行性,會(huì)對(duì)最終的運(yùn)算效率產(chǎn)生較大的影響。因此,本文提出了一種新的存儲(chǔ)格式DRM(Divide-Rearrange &Merge)。該格式利用基于固定閾值的矩陣劃分策略和基于迭代歸并的矩陣重構(gòu)策略,對(duì)稀疏矩陣進(jìn)行劃分,并對(duì)劃分結(jié)果進(jìn)行迭代歸并,在解決零填充問題的同時(shí)達(dá)到了負(fù)載平衡的目的。

2 相關(guān)工作

目前,國內(nèi)外很多研究人員對(duì)稀疏矩陣向量乘進(jìn)行了廣泛的研究,提出了一系列優(yōu)化方法,針對(duì)不同問題的優(yōu)化主要分為2類:最小化零填充量和負(fù)載平衡。

2.1 最小化零填充量

Choi等人[17]提出了BCSR(Blocked Compressed Sparse Row)和BELL(Blocked ELLpack)存儲(chǔ)格式,通過對(duì)矩陣進(jìn)行分塊處理,提高了SpMV的性能,但會(huì)出現(xiàn)塊間負(fù)載不平衡的問題。Yuan等人[13]提出了DDD-SPLIT格式,利用對(duì)角線將稀疏矩陣分成若干子塊,壓縮對(duì)角線上連續(xù)的相同元素,進(jìn)而壓縮具有相同值的子塊,節(jié)省了大量的存儲(chǔ)空間。Barbieri等人[14]提出了HDIA格式,通過將原始矩陣分成大小相等的行組,并將行組以DIA格式存儲(chǔ),減少了零元素的填充量和緩解了緩存未命中的問題。但是,該格式可能會(huì)出現(xiàn)塊間元素?cái)?shù)量差異過大的現(xiàn)象。Gao等人[16]提出了DIA-Adaptive格式,將矩陣分為3種類型,并自適應(yīng)地選擇DIA、BRCSD-Ⅰ(Diagonal Compressed Storage based on Row-Blocks-Ⅰ)、BRCSD-Ⅱ格式中的一種對(duì)矩陣進(jìn)行存儲(chǔ),改善了稀疏對(duì)角矩陣存在長(zhǎng)偏移、散點(diǎn)和長(zhǎng)斷行時(shí)的零填充問題,然而這樣可能會(huì)出現(xiàn)負(fù)載不平衡的問題。Sun等人[11,12]提出了CRSD存儲(chǔ)格式,通過設(shè)計(jì)對(duì)角線格式,將矩陣的對(duì)角線分成不同的組,并在行方向上將矩陣劃分成行片段,解決了稀疏對(duì)角矩陣存在散點(diǎn)和長(zhǎng)斷行時(shí)的零填充問題。但是,由于不同對(duì)角線格式中對(duì)角線的數(shù)量不一致,該格式會(huì)出現(xiàn)行片段間負(fù)載不平衡的問題。Yang等人[18,19]提出了一種混合存儲(chǔ)格式HDC(Hybrid of DIA and CSR(Compressed Sparse Row)),通過設(shè)定分割閾值,將矩陣分為對(duì)角區(qū)域和離散區(qū)域,并分別采用DIA格式和CSR格式進(jìn)行存儲(chǔ),有效地減少了矩陣存儲(chǔ)時(shí)的零填充量。但是,并行SpMV運(yùn)算性能受分割閾值影響較大,最佳閾值又難以確定。

2.2 負(fù)載平衡

Liu等人[20]提出了CSR5格式,通過將非零元均勻地劃分為多個(gè)相同大小的2D切片,讓一個(gè)SIMD(Single Instruction Multiple Data)通道處理切片的一列,使稀疏矩陣達(dá)到近乎最佳的負(fù)載平衡。但是,從CSR到CSR5的預(yù)處理時(shí)間需要比SpMV的多10～30倍。Yang等人[21]提出了一種基于非零元分布的最優(yōu)分塊策略和一種重排序算法,并利用混合格式BCE(Blocked stored format mixed CSR and ELL)對(duì)矩陣進(jìn)行壓縮存儲(chǔ),在提升SpMV性能的同時(shí),解決了分區(qū)后的負(fù)載平衡問題。Cui等人[22]提出了一種自動(dòng)獲取最優(yōu)塊劃分的策略和一種基于分塊策略的矩陣重排方法,并采用PBC(Partition-Block-CSR)格式來存儲(chǔ)重排序后的矩陣,分塊策略和矩陣重排保證了分塊之間的標(biāo)準(zhǔn)差最小,既滿足了負(fù)載均衡條件,又避免了零元素填充,有效提高了并行計(jì)算效率,但是該方法存在預(yù)處理時(shí)間過長(zhǎng)的問題。Bian等人[23]提出了CSR2(Compressed Sparse Row 2)存儲(chǔ)格式,它可以根據(jù)矩陣的特點(diǎn)自動(dòng)確定單次執(zhí)行的矩陣分塊寬度,并借助SIMD函數(shù)在平臺(tái)上充分發(fā)揮SIMD向量化的效率。與CSR5相比,CSR2具有更少的存儲(chǔ)空間消耗和更低的格式轉(zhuǎn)換開銷。Gao等人[24]提出了AMF-CSR(Adaptive Multi-row Folding of CSR)存儲(chǔ)格式,通過對(duì)多個(gè)非零行進(jìn)行折疊,并合并非零數(shù)量相近的行,增加了行中非零元的密度,提升了x向量的訪問局部性,改善了線程塊間負(fù)載不均衡的問題。但是,部分閾值根據(jù)經(jīng)驗(yàn)設(shè)置,可能存在最后一塊與其他塊的非零元數(shù)量不平衡的現(xiàn)象。

從以上分析可知,目前,針對(duì)稀疏對(duì)角矩陣的主流存儲(chǔ)格式雖然能夠在一定程度上提升并行SpMV效率,但仍然會(huì)出現(xiàn)零填充過度和GPU線程塊間負(fù)載不平衡的問題。

3 稀疏矩陣存儲(chǔ)格式介紹

3.1 DIA存儲(chǔ)格式

DIA存儲(chǔ)格式是處理稀疏對(duì)角矩陣的一種經(jīng)典格式。對(duì)于具有k條對(duì)角線的n維稀疏矩陣,該格式用2個(gè)數(shù)組來存儲(chǔ)矩陣,一個(gè)是大小為k×n的data數(shù)組,用來存儲(chǔ)每條對(duì)角線的元素;另一個(gè)是大小為k的offsets數(shù)組,用來存儲(chǔ)每條對(duì)角線相對(duì)于主對(duì)角線的偏移量。稀疏矩陣示例和DIA存儲(chǔ)格式如圖1所示,其中Ci和Ri表示稀疏矩陣的列向量和行向量。

DIA格式適合存儲(chǔ)存在密集對(duì)角線的稀疏對(duì)角矩陣。若矩陣中存在散點(diǎn)、長(zhǎng)斷行和與主對(duì)角線偏移較大的對(duì)角線時(shí),會(huì)填充大量的零元素,這在浪費(fèi)存儲(chǔ)空間的同時(shí)還會(huì)影響GPU上SpMV的運(yùn)算效率和時(shí)間性能。

a 稀疏矩陣

offsets=[-3,0,1,3]

3.2 HDIA存儲(chǔ)格式

HDIA存儲(chǔ)格式是在DIA格式的基礎(chǔ)上加入了分塊的思想。通過設(shè)置固定的行塊大小,將稀疏矩陣劃分為大小相等的行塊,并對(duì)每一行塊采用DIA格式進(jìn)行存儲(chǔ)。該格式利用3個(gè)數(shù)組存儲(chǔ)矩陣,hackOffsets數(shù)組用于存儲(chǔ)每個(gè)行塊的對(duì)角線個(gè)數(shù),hdi_offsets數(shù)組用于存儲(chǔ)每個(gè)行塊中對(duì)角線的偏移量,hdi_data數(shù)組用于存儲(chǔ)每個(gè)行塊中的元素。對(duì)于圖1中的矩陣A,當(dāng)行塊大小設(shè)置為2時(shí),HDIA存儲(chǔ)格式如圖2所示,其中bi表示矩陣劃分后的行塊。

hackOffsets=[0,3,6,8]

HDIA存儲(chǔ)格式有效地解決了矩陣中存在長(zhǎng)斷行和散點(diǎn)時(shí)的零填充問題,但由于每個(gè)塊中非零元個(gè)數(shù)不一樣,甚至差異很大,可能會(huì)導(dǎo)致塊間的負(fù)載不平衡,從而影響并行SpMV的性能。

3.3 DIA-Adaptive存儲(chǔ)格式

DIA-Adaptive是一種自適應(yīng)存儲(chǔ)格式。首先通過稀疏矩陣的特征對(duì)矩陣進(jìn)行類別判斷,然后根據(jù)矩陣的類別對(duì)矩陣采用DIA、BRCSD-Ⅰ、BRCSD-Ⅱ中的一種格式進(jìn)行壓縮存儲(chǔ)。其中,BRCSD-Ⅰ格式利用矩陣的偏移量和一種計(jì)算策略來獲得矩陣的有效片點(diǎn)和交叉片點(diǎn),利用交叉片點(diǎn)對(duì)矩陣進(jìn)行分塊。采用2個(gè)數(shù)組存儲(chǔ)矩陣,brcsdⅠ_offsets數(shù)組用于存儲(chǔ)每個(gè)塊的起始行號(hào)和塊內(nèi)對(duì)角線的偏移量,brcsdⅠ_data數(shù)組用于存儲(chǔ)每個(gè)塊的元素。對(duì)于圖1中的矩陣A,BRCSD-Ⅰ存儲(chǔ)格式如圖3所示,其中offsets[i]表示對(duì)矩陣該部分的偏移量進(jìn)行統(tǒng)一存儲(chǔ)。

brcsdⅠ_offsets={0|[0,1,3],3|[-3,0,1]}

BRCSD-Ⅱ格式通過設(shè)置固定的nrows值,將矩陣分成大小相等的行塊,并利用一定的合并策略將具有相同偏移量的塊進(jìn)行偏移量壓縮。采用2個(gè)數(shù)組存儲(chǔ)矩陣,brcsdⅡ_offsets數(shù)組用于存儲(chǔ)相同偏移量的塊數(shù)和對(duì)應(yīng)的偏移量,brcsdⅡ_data數(shù)組用于存儲(chǔ)每個(gè)塊的元素。對(duì)于圖1中的矩陣A,nrows=2時(shí)的BRCSD-Ⅱ存儲(chǔ)格式如圖4所示。

brcsdⅡ_offsets={2|[0,1,3],1|[-3,0]}

BRCSD-Ⅰ存儲(chǔ)格式有效解決了矩陣存在長(zhǎng)偏移時(shí)的零填充問題,但由于每塊行數(shù)不一致,可能會(huì)存在塊間負(fù)載不平衡現(xiàn)象。BRCSD-Ⅱ存儲(chǔ)格式有效解決了矩陣存在散點(diǎn)和長(zhǎng)斷行時(shí)的零填充問題,雖然每塊行數(shù)相等,但仍可能存在塊間負(fù)載不平衡的現(xiàn)象。

4 DRM格式以及并行SpMV實(shí)現(xiàn)

針對(duì)稀疏對(duì)角矩陣存儲(chǔ)時(shí)的零填充問題以及現(xiàn)有分塊格式存在的負(fù)載不平衡問題,本文提出了一種DRM存儲(chǔ)格式。該格式利用基于固定閾值的矩陣劃分策略和基于迭代歸并的矩陣重構(gòu)策略,將輸入矩陣劃分為大小相等的行段,對(duì)行段進(jìn)行迭代歸并,旨在減少矩陣存儲(chǔ)時(shí)的零填充量,實(shí)現(xiàn)子塊間負(fù)載平衡,提升SpMV在GPU上的并行效率。

4.1 基于固定閾值的矩陣劃分策略

DIA存儲(chǔ)格式適合存儲(chǔ)對(duì)角線密集的稀疏矩陣,但是在實(shí)際工程計(jì)算問題中得到的稀疏對(duì)角矩陣往往存在一些散點(diǎn)和長(zhǎng)斷行。若直接采用DIA格式進(jìn)行存儲(chǔ),會(huì)出現(xiàn)零填充問題,當(dāng)矩陣維數(shù)較大時(shí),填充的零元素個(gè)數(shù)甚至?xí)h(yuǎn)超非零元素的個(gè)數(shù),這樣不僅浪費(fèi)內(nèi)存空間,還會(huì)嚴(yán)重影響DIA格式的時(shí)間性能。因此,DRM存儲(chǔ)格式首先對(duì)稀疏對(duì)角矩陣進(jìn)行一次固定閾值(nrows)的劃分,使對(duì)角線按行段進(jìn)行存儲(chǔ),這樣就有效解決了矩陣存儲(chǔ)時(shí)的零填充問題。為了說明劃分的重要性,本文進(jìn)行了深入分析。

定義1假設(shè)稀疏矩陣B分為K塊,每個(gè)塊表示為B0,B1,B2,…,BK-1,若塊之間滿足式(1)和式(2),則稱Bi為稀疏矩陣B的子塊:

Bi∩Bj=?,i,j=0,1,…,K-1,i≠j

(1)

(2)

定義2假設(shè)子塊Bi由若干小塊組成,每個(gè)小塊表示為bi0,bi1,bi2,…,若小塊之間滿足式(3)和式(4),則稱bij為行段:

bij∩bik=?,j≠k

(3)

(4)

定義3對(duì)一個(gè)N維的稀疏對(duì)角矩陣,以nrows為行段大小對(duì)矩陣進(jìn)行劃分,矩陣所分的行段數(shù)如式(5)所示,行段中的操作數(shù)個(gè)數(shù)如式(6)所示:

bl=N/nrows

(5)

opeh=digh×nrows,h=0,1,…,bl-1

(6)

其中,bl為矩陣所分的行段數(shù),opeh為行段中的操作數(shù)個(gè)數(shù),digh為行段中的對(duì)角線個(gè)數(shù)。操作數(shù)是矩陣進(jìn)行SpMV運(yùn)算時(shí)需要參與運(yùn)算的元素個(gè)數(shù)。由于矩陣存儲(chǔ)時(shí)的零填充問題,除了非零元外,一部分零元素也會(huì)參與運(yùn)算,同樣屬于操作數(shù)。

定義4對(duì)定義3中的矩陣,若矩陣包含S條對(duì)角線,DIA格式壓縮矩陣時(shí),需要存儲(chǔ)的操作數(shù)個(gè)數(shù)如式(7)所示:

opeDIA=S×N

(7)

以nrows為行段大小對(duì)矩陣進(jìn)行劃分,劃分后需要存儲(chǔ)的操作數(shù)個(gè)數(shù)如式(8)所示:

h=0,1,…,bl-1,digh≤S

(8)

其中,opeDRM滿足式(9):

nrows∈{1,2,…,N/2}

(9)

由式(9)可知,無論nrows為何值,都滿足opeDRM≤opeDIA。假設(shè)矩陣B在定義4的基礎(chǔ)上存在一個(gè)散點(diǎn),則操作數(shù)個(gè)數(shù)的變化如式(10)和式(11)所示:

opeDIA=(S+1)×N

(10)

S×N+nrows

(11)

由式(10)和式(11)可推出式(12):

opeDRM

(12)

因此,當(dāng)稀疏矩陣中存在散點(diǎn)或長(zhǎng)斷行時(shí),相比DIA存儲(chǔ)格式,對(duì)矩陣的劃分一定能減少矩陣存儲(chǔ)時(shí)的零填充量。對(duì)于存在散點(diǎn)或長(zhǎng)斷行較多的矩陣,通過劃分減少的零元素填充更多。

以圖5的稀疏矩陣B為例,當(dāng)直接使用DIA格式存儲(chǔ)該矩陣時(shí),opeDIA=32,填充了14個(gè)零元素。若先使用劃分策略對(duì)矩陣進(jìn)行劃分,當(dāng)nrows=4時(shí),bl=2,opeDRM=24,需要填充6個(gè)零元素;當(dāng)nrows=2時(shí),bl=4,opeDRM=20,只需要填充2個(gè)零元素。因此,對(duì)稀疏矩陣進(jìn)行一次固定閾值的劃分是必要的,能有效減少矩陣存儲(chǔ)時(shí)的零填充量。nrows=2時(shí)矩陣劃分結(jié)果如圖6所示,本文稱圖6中的矩陣為矩陣C。

Figure 5 Example of sparse matrix B圖5 稀疏矩陣B示例

Figure 6 Partition results of the matrix when nrows is 2圖6 nrows=2時(shí)矩陣的劃分結(jié)果

本文以nrows=2為例進(jìn)行說明,但實(shí)驗(yàn)時(shí)的nrows固定為32。首先,在GPU架構(gòu)中線程調(diào)度的基本單位是warp,一個(gè)warp包含32個(gè)線程,若同一warp內(nèi)的線程訪存地址一致或相鄰,可以減少線程的訪存時(shí)間,從而提升SpMV的性能。其次,本文針對(duì)nrows進(jìn)行了大量的測(cè)試和實(shí)驗(yàn),發(fā)現(xiàn)nrows的大小和SpMV性能成反比,但隨著nrows的減小,預(yù)處理時(shí)間會(huì)逐漸增大。因此,為了達(dá)到一個(gè)較好的性能平衡,將nrows固定為32。

4.2 基于迭代歸并的矩陣重構(gòu)策略

對(duì)稀疏對(duì)角矩陣進(jìn)行一次固定閾值的劃分后,減少了矩陣存儲(chǔ)時(shí)的零填充量,但由于行段間的操作數(shù)個(gè)數(shù)不一致,很容易出現(xiàn)GPU線程塊間負(fù)載不平衡的問題。針對(duì)此問題,本節(jié)提出了一種基于迭代歸并的矩陣重構(gòu)策略,先將矩陣劃分的行段按照操作數(shù)個(gè)數(shù)進(jìn)行重新排列,再根據(jù)排列結(jié)果將對(duì)應(yīng)行段迭代歸并成子塊,實(shí)現(xiàn)了子塊間的方差最小化,滿足了并行SpMV運(yùn)算時(shí)負(fù)載平衡的條件。

定義5對(duì)稀疏對(duì)角矩陣進(jìn)行一次固定閾值的劃分后,行段間操作數(shù)個(gè)數(shù)的方差如式(13)所示:

(13)

其中,ave是行段中操作數(shù)個(gè)數(shù)的平均值,bl是矩陣以nrows為閾值進(jìn)行劃分后的總行段數(shù),var是行段間操作數(shù)個(gè)數(shù)的方差。方差是衡量一組數(shù)據(jù)平衡與否的標(biāo)準(zhǔn),方差越小,說明數(shù)據(jù)間的差異越小。矩陣劃分后,每個(gè)行段中操作數(shù)個(gè)數(shù)的大小是否相似,決定了并行SpMV運(yùn)算時(shí)矩陣的負(fù)載是否平衡。

負(fù)載不平衡會(huì)導(dǎo)致線程發(fā)散,引起計(jì)算資源的浪費(fèi),同時(shí)還會(huì)影響GPU上并行SpMV的計(jì)算速度,降低整個(gè)計(jì)算過程的效率。因此,DRM格式使用基于迭代歸并的矩陣重構(gòu)策略,對(duì)矩陣劃分后的行段進(jìn)行迭代歸并,能夠有效實(shí)現(xiàn)GPU線程塊間計(jì)算量的最大平衡。

定義6對(duì)一個(gè)集合P,假設(shè)P含有l(wèi)個(gè)元素,即P={p0,p1,…,pl-1}。若pmax1、pmin1、pmax2和pmin2滿足式(14)～式(17),則分別稱為集合P中的最大值、最小值、次大值和次小值:

pmax1≥pr,r=0,1,…,l-1

(14)

pmin1≤pr,r=0,1,…,l-1

(15)

pmax1≥pmax2≥ps1,ps1∈(P-{pmax1})

(16)

pmin1≤pmin2≤ps2,ps2∈(P-{pmin1})

(17)

基于迭代歸并的矩陣重構(gòu)策略如下所示:

(1)構(gòu)造操作數(shù)集合P?；诿總€(gè)行段中對(duì)角線的個(gè)數(shù)digh和矩陣劃分閾值nrows,利用式(6)計(jì)算每個(gè)行段中操作數(shù)個(gè)數(shù)opeh,并將其存入集合P。

(2)整理集合P。P中含有bl個(gè)元素,可表示為P={p0,p1,…,pbl-1},循環(huán)執(zhí)行以下步驟:

①對(duì)集合P按照元素遞減的順序重排,重排后的P={pmax1,pmax2,…,pmin2,pmin1}。

②取出集合P中的最大值pmax1和最小值pmin1,如果滿足pmax1>2×pmin1,則執(zhí)行③。否則,結(jié)束循環(huán),執(zhí)行(3)。

③取出集合P中的次小值pmin2,如果滿足pmax1>2×pmin2,則將最小值pmin1和次小值pmin2對(duì)應(yīng)的行段進(jìn)行歸并,并將2個(gè)行段中的操作數(shù)個(gè)數(shù)相加,即psum=pmin1+pmin2,將psum存入集合P中。否則,結(jié)束循環(huán),執(zhí)行(3)。

(3)判斷集合P中元素的總數(shù)card(P)是奇數(shù)還是偶數(shù)。若是奇數(shù),則把最小值pmin1和次小值pmin2對(duì)應(yīng)的行段進(jìn)行歸并,并將2個(gè)行段中的操作數(shù)個(gè)數(shù)之和存入集合P中,再根據(jù)遞減順序進(jìn)行一次重排。若是偶數(shù),直接執(zhí)行(4)。

(4)判斷card(P)是否滿足card(P)>2,若滿足,則按照最大值pmax1和最小值pmin1、次大值pmax2和次小值pmin2歸并的方式,將集合P中每個(gè)值pr對(duì)應(yīng)的行段歸并為子塊,得到最終的歸并結(jié)果。否則,直接將集合P當(dāng)作最終歸并結(jié)果。

經(jīng)過以上步驟后,每個(gè)子塊的操作數(shù)個(gè)數(shù)將趨近于平衡,使子塊間操作數(shù)個(gè)數(shù)的方差變小,從而實(shí)現(xiàn)了GPU線程塊間的負(fù)載平衡。基于迭代歸并的矩陣重構(gòu)算法如算法1所示。

算法1 基于迭代歸并的矩陣重構(gòu)算法輸入:稀疏矩陣B,分塊閾值nrows。輸出:重構(gòu)后的操作數(shù)集合P。1:P ← Getoperands(B,nrows);2:sort(P.begin(),P.end(),descending);3:while pmax1> 2*pmin1 and pmax1> 2*pmin2do4: pmin2+=pmin1;5: P.pop_back();6: sort(P.begin(),P.end(),descending);7:end while8:int len = P.size();9:if len mod 2=0 then10: for r=0 → lendo11: if r < (len -1- r)then12: pr+=plen-1-r;13: P.pop_back();14: end if15: end for16:else17: plen-2+=plen-1;18: P.pop_back();19: sort(P.begin(),P.end(),descending);20: len=P.size();21: for r=0 → lendo22: if r < len -1- rthen23: pr+=plen-1-r;24: P.pop_back();25: end if26: end for27:end if

算法1的第1～2行計(jì)算矩陣劃分后每個(gè)行段中的操作數(shù)個(gè)數(shù),將其存入集合P,并對(duì)集合P中元素進(jìn)行遞減排序。第3～7行循環(huán)歸并滿足條件的行段。第8行統(tǒng)計(jì)集合P的大小,若為偶數(shù)(第9～15行),直接將集合P中對(duì)應(yīng)的行段歸并為子塊;若為奇數(shù)(第16～27行),則先歸并最后2個(gè)行段,并進(jìn)行一次遞減排序,再將剩余的行段歸并為子塊。

以矩陣C為例,行段迭代歸并過程如圖7所示。

Figure 7 Iterative merging process of line segments圖7 行段迭代歸并過程

首先利用式(6)計(jì)算每個(gè)行段中的操作數(shù)個(gè)數(shù),并將其存入集合P。然后,對(duì)集合P進(jìn)行遞減排序,取出P中的最大值8和最小值2,由于8>2×2=4,滿足循環(huán)執(zhí)行步驟②,取出集合P的次小值4,由于8=2×4,不滿足循環(huán)執(zhí)行步驟③,則執(zhí)行下一步。最后,判斷card(P)是奇數(shù)還是偶數(shù),由于card(P)=4,且大于2,則根據(jù)pr的大小,將集合P中的行段兩兩歸并為子塊,得到最終集合P={10,10}。

基于迭代歸并的矩陣重構(gòu)結(jié)果如圖8所示,其中,B1和B2為最終歸并的子塊,每個(gè)子塊含有2個(gè)行段。通過式(13)計(jì)算可知,子塊間操作數(shù)個(gè)數(shù)的方差為0,表示稀疏矩陣已達(dá)到完全的負(fù)載平衡狀態(tài)。此時(shí),在GPU上進(jìn)行SpMV運(yùn)算時(shí)會(huì)有效平衡線程塊間計(jì)算負(fù)載,提升運(yùn)算的效率和性能。

Figure 8 Matrix reconstruction results based on iterative merging圖8 基于迭代歸并的矩陣重構(gòu)結(jié)果

4.3 DRM存儲(chǔ)格式

DRM存儲(chǔ)格式首先基于固定閾值的矩陣劃分策略,將矩陣劃分為多個(gè)行段,再結(jié)合4.2節(jié)中的迭代歸并策略,將多個(gè)行段歸并為子塊,最終實(shí)現(xiàn)子塊間opei的最大均衡,有效解決了稀疏對(duì)角矩陣在并行SpMV運(yùn)算時(shí)的零填充和負(fù)載不平衡問題。

在GPU上進(jìn)行并行SpMV運(yùn)算之前,需要將矩陣壓縮的數(shù)據(jù)存儲(chǔ)到數(shù)組中。因此,首先建立與計(jì)算相關(guān)的數(shù)組。數(shù)組drm_dnum以累加的方式存儲(chǔ)子塊中每個(gè)行段中的對(duì)角線個(gè)數(shù),用于行段中操作數(shù)個(gè)數(shù)的計(jì)算。數(shù)組drm_index存儲(chǔ)每個(gè)行段歸并后的索引用于SpMV運(yùn)算時(shí)矩陣行索引。數(shù)組drm_blo以累加的方式存儲(chǔ)子塊中包含的行段個(gè)數(shù),用于對(duì)每個(gè)行段進(jìn)行索引。數(shù)組drm_off和drm_val分別存儲(chǔ)子塊中每一行段內(nèi)對(duì)角線的偏移量和元素值。

對(duì)圖5中的矩陣B,利用上述5個(gè)數(shù)組存儲(chǔ)的結(jié)果如圖9所示。

Figure 9 DRM storage format representation of matrix B圖9 矩陣B的DRM存儲(chǔ)格式表示

4.4 基于DRM存儲(chǔ)格式的并行SpMV算法

本文在4.2節(jié)中提供了一個(gè)矩陣重構(gòu)策略,實(shí)現(xiàn)了GPU線程塊間計(jì)算量的最大相似性,在本節(jié)進(jìn)一步提出一種新的基于DRM存儲(chǔ)格式的并行SpMV算法。由于DRM將稀疏對(duì)角矩陣劃分為行段,并對(duì)行段進(jìn)行了迭代歸并,導(dǎo)致子塊間的行段是亂序排列的。為了與劃分結(jié)果相對(duì)應(yīng),需要對(duì)矩陣進(jìn)行有效的重排。因此,在SpMV運(yùn)算時(shí)需要準(zhǔn)確地獲取當(dāng)前行的索引,避免最終計(jì)算結(jié)果出現(xiàn)錯(cuò)誤。此外,為了體現(xiàn)負(fù)載平衡對(duì)SpMV計(jì)算的加速效果,GPU處理時(shí)采用一個(gè)線程塊處理一個(gè)子塊的方式,在每個(gè)線程塊中一個(gè)線程處理矩陣的一行。若子塊中的行數(shù)超出硬件支持的線程數(shù)目(1 024),則在預(yù)處理階段對(duì)該子塊進(jìn)行分解處理。DRM存儲(chǔ)格式的并行SpMV算法如算法2所示。在算法2的第1～2行,使用bid和tid分別記錄線程塊的塊號(hào)和線程號(hào),便于索引子塊。第4～5行使用nrows和drm_blo數(shù)組創(chuàng)建變量rnum和rinx,計(jì)算子塊中的行段索引。第6行使用drm_index數(shù)組來計(jì)算行索引,以確保乘積結(jié)果相加的正確性。第7～14行是本文算法的核心,利用上述索引以及drm_dnum、drm_off和drm_val數(shù)組,采用一個(gè)線程計(jì)算一行的方式,對(duì)矩陣每一行的數(shù)據(jù)進(jìn)行乘積運(yùn)算和累加操作,得到最終的SpMV結(jié)果,并將結(jié)果傳入result數(shù)組。

算法2 DRM存儲(chǔ)格式的并行SpMV算法輸入:drm_dnum,drm_index,drm_blo,drm_off,drm_val,x,nrows。輸出:結(jié)果向量result。1:int bid=blockIdx.x;2:int tid=threadIdx.x;3:double res=0;4:int rnum=tid /nrows;5:int rinx=rnum + drm_blo [bid];6:int row=tid +(drm_index [rinx]- rnum)*nrows;7:if rnum < (drm_blo [bid +1] - drm_blo [bid]) then8: for i =drm_dnum [rinx]→drm_dnum [rinx +1] do9: int val_index = (tid -rnum*nrows) +i*nrows;10: int x_index=row+drm_off [i];11: res+= drm_val [val_index] *x[x_index];12: end for13: result[row]=res;14:end if

算法2利用5個(gè)矩陣壓縮數(shù)組對(duì)每個(gè)非零元進(jìn)行精準(zhǔn)索引,實(shí)現(xiàn)了DRM存儲(chǔ)格式的并行SpMV內(nèi)核。由于在4.1節(jié)和4.2節(jié)對(duì)矩陣進(jìn)行了劃分,并對(duì)劃分后的行段進(jìn)行了歸并,使并行SpMV過程中每個(gè)線程塊間的計(jì)算量達(dá)到了平衡,有效減少了線程的空閑等待,提升了整個(gè)運(yùn)算過程的效率和性能。

5 實(shí)驗(yàn)評(píng)估

本文在第4節(jié)提出了一種新的存儲(chǔ)格式DRM,并實(shí)現(xiàn)了基于此格式的并行SpMV算法。本節(jié)利用零填充數(shù)量、操作數(shù)個(gè)數(shù)方差、時(shí)間和GFLOPS等性能指標(biāo),通過與當(dāng)前主流的DIA、HDC、HDIA和DIA-Adaptive存儲(chǔ)格式進(jìn)行對(duì)比實(shí)驗(yàn),驗(yàn)證所提并行SpMV算法的有效性。

5.1 實(shí)驗(yàn)環(huán)境

本文實(shí)驗(yàn)采用的GPU型號(hào)為NVIDIA?Tesla?V100-PCIE 32 GB,包含5 120個(gè)CUDA核心SP(Stream Processor);內(nèi)存容量為32 GB,含有80個(gè)流多處理器SM(Stream Multiprocessor);每個(gè)線程塊的最大線程數(shù)是1 024;CPU為Intel?Xeon?Gold 6240R,主頻為2.4 GHz,24核;服務(wù)器系統(tǒng)為Ubuntu 16.04.7,CUDA版本是11.4。

本文實(shí)驗(yàn)采用的稀疏矩陣為方陣,全部來自佛羅里達(dá)大學(xué)公開矩陣數(shù)據(jù)集(SuiteSparse Collection)[25]。表1描述了這些稀疏矩陣的相關(guān)信息。主要包括矩陣名、矩陣縮略圖、維數(shù)和非零元個(gè)數(shù)。所選擇的稀疏矩陣代表了不同類型的實(shí)際應(yīng)用,包括電路仿真問題、計(jì)算流體力學(xué)問題、結(jié)構(gòu)問題和電磁學(xué)問題等。

在實(shí)驗(yàn)中,將HDIA和DIA-Adaptive中BRCSD-Ⅱ格式的分塊閾值設(shè)為32,將BRCSD-Ⅰ格式分塊后的再分塊閾值設(shè)為64。

5.2 零填充和負(fù)載平衡

本節(jié)從零填充和負(fù)載平衡2方面對(duì)DRM格式進(jìn)行分析評(píng)估。首先驗(yàn)證該格式在減少稀疏對(duì)角矩陣零填充量方面的有效性。由于實(shí)驗(yàn)中HDIA和DRM格式的分塊閾值相同,所以具有相同的零填充數(shù)量。因此,本節(jié)只對(duì)DIA、HDC、DIA-Adaptive和DRM格式進(jìn)行了對(duì)比。如圖10所示,在所有測(cè)試矩陣中,HDC格式取得了最少的零元素填充量,DRM格式的零填充數(shù)量小于DIA和DIA-Adaptive格式的。

對(duì)于稀疏對(duì)角矩陣來說,將非零元素按照對(duì)角線進(jìn)行存儲(chǔ)有利于維持并行SpMV過程中的訪存連續(xù)性。而HDC格式對(duì)矩陣進(jìn)行存儲(chǔ)時(shí),僅將少量對(duì)角線按照DIA的形式進(jìn)行存儲(chǔ),剩余非零元素全部按照CSR的形式進(jìn)行存儲(chǔ),這樣雖然在零填充方面表現(xiàn)較好,但會(huì)引起SpMV過程中不連續(xù)的內(nèi)存訪問,導(dǎo)致降低整個(gè)運(yùn)算過程的性能。

除HDC格式外,其余格式均按照對(duì)角線的形式對(duì)稀疏矩陣進(jìn)行存儲(chǔ)。其中,相比于DIA格式,DRM零元素填充數(shù)量減少地最為顯著,在矩陣add32、lns_3937、ex11、rw5151和shyy161中,DIA格式填充零元素的個(gè)數(shù)分別為DRM格式的112.80,87.70,107.30,112.60和292.80倍。其中,add32和ex11填充零元素的個(gè)數(shù)是矩陣中非零元總數(shù)的940.30和82.70倍,而在DRM格式下只有8.30和0.07倍。通過觀察矩陣的結(jié)構(gòu)發(fā)現(xiàn),其原因在于這些矩陣中存在大量的散點(diǎn)、長(zhǎng)斷行和長(zhǎng)偏移,從而導(dǎo)致DIA格式存儲(chǔ)的對(duì)角線條數(shù)急劇增加,從而提高了零元素的填充比例。以矩陣add32為例,采用DIA格式存儲(chǔ)的對(duì)角線有3 767條。

Table 1 Information of sparse matrix dataset表1 稀疏矩陣數(shù)據(jù)集相關(guān)信息

Figure 10 Number of zero padding for different storage formats圖10 不同存儲(chǔ)格式的零填充數(shù)量

相比于DIA-Adaptive格式,DRM雖然不如相比DIA格式那么顯著,但每個(gè)測(cè)試矩陣的零填充數(shù)量都有所減少。在矩陣lns_3937、rw5151、shipsec5和add32中,DIA-Adaptive格式的零元素填充數(shù)量分別是DRM格式的3.42,3.05,2.62和2.63倍。由此可知,DRM存儲(chǔ)格式經(jīng)過一次固定閾值的矩陣劃分后,能夠有效減少稀疏對(duì)角矩陣存儲(chǔ)時(shí)的零填充數(shù)量。

然后,驗(yàn)證該格式在解決稀疏對(duì)角矩陣負(fù)載不平衡方面的有效性。為了更直觀地展現(xiàn)并行SpMV計(jì)算時(shí)矩陣的負(fù)載是否平衡,將子塊間操作數(shù)個(gè)數(shù)的方差作為評(píng)估負(fù)載平衡的標(biāo)準(zhǔn)。如式(13)所示,方差越小,說明子塊間操作數(shù)個(gè)數(shù)的差異性越小,GPU線程塊間的計(jì)算量越平衡。

本文將DRM格式與HDIA、DIA-Adaptive格式進(jìn)行對(duì)比,統(tǒng)計(jì)矩陣分塊后每塊操作數(shù)個(gè)數(shù)的方差,作為最終的實(shí)驗(yàn)結(jié)果。如圖11所示,DRM存儲(chǔ)格式在所有測(cè)試矩陣中的子塊間操作數(shù)個(gè)數(shù)的方差均小于DIA-Adaptive格式的。在除nemeth22外的其余測(cè)試矩陣中,DRM的子塊間操作數(shù)個(gè)數(shù)方差均小于HDIA格式的。最小方差出現(xiàn)在矩陣dw2048中,因?yàn)樵摼仃嚲S數(shù)較小且對(duì)角結(jié)構(gòu)明顯,非零元分布比較規(guī)律,有利于DRM格式的行段歸并。最大方差出現(xiàn)在矩陣ex11中,在該矩陣中,非零元分布不規(guī)則使矩陣對(duì)角線數(shù)量驟增,DRM行段歸并后的結(jié)果欠佳,導(dǎo)致子塊間操作數(shù)個(gè)數(shù)有一定的差異,但仍優(yōu)于HDIA和DIA-Adaptive格式的。

Figure 11 Comparison of operand variance between sub-blocks圖11 子塊間操作數(shù)個(gè)數(shù)的方差對(duì)比

在矩陣nemeth22中,非零元素沿著主對(duì)角線周圍密集排列,HDIA格式對(duì)其分塊時(shí),每塊操作數(shù)個(gè)數(shù)已經(jīng)達(dá)到了一定的平衡,因此DRM格式對(duì)行段進(jìn)行迭代歸并后,子塊間操作數(shù)個(gè)數(shù)的方差略大于HDIA格式的。此外,矩陣Kuu、shyy161、Goodwin_095和Goodwin_127中包含很多散點(diǎn)和長(zhǎng)斷行,在采用DIA-Adaptive格式存儲(chǔ)時(shí),被自適應(yīng)地改為BRCSD-Ⅱ格式進(jìn)行存儲(chǔ)。因?yàn)锽RCSD-Ⅱ和HDIA具有相同的分塊閾值,所以在這4個(gè)矩陣中,DIA-Adaptive和HDIA格式的每塊操作數(shù)個(gè)數(shù)方差相同,但都大于DRM格式的。

圖12顯示了HDIA和DIA-Adaptive的方差與DRM方差的比值。從圖12中可以看出,DRM格式相比與另外2種格式,在負(fù)載平衡方面有很大的改進(jìn)。DIA-Adaptive格式的方差是DRM格式的平均66.70倍(最高197.20倍),HDIA的方差是DRM格式的平均8.20倍(最高24.80倍)。這表明,DRM格式利用矩陣重構(gòu)策略,對(duì)劃分的行段進(jìn)行迭代歸并,極大地減少了子塊間操作數(shù)個(gè)數(shù)的方差,實(shí)現(xiàn)了GPU線程塊間計(jì)算量的最大平衡。

Figure 12 Variance ratios of HDIA and DIA-Adaptive to DRM圖12 HDIA和DIA-Adaptive與DRM方差的比值

5.3 SpMV執(zhí)行時(shí)間和GFLOPS

基于本文提出的DRM格式,本節(jié)采用時(shí)間和GFLOPS作為性能指標(biāo),評(píng)價(jià)DRM格式的性能。在此過程中,分別采用單精度和雙精度進(jìn)行比較。

首先對(duì)時(shí)間性能進(jìn)行對(duì)比,為了方便比較不同存儲(chǔ)格式的性能差異,避免初始化、數(shù)據(jù)拷貝和邏輯判斷等過程中由于矩陣大小或數(shù)據(jù)傳輸速率不同而造成的誤差,本節(jié)對(duì)比的時(shí)間為SpMV內(nèi)核執(zhí)行時(shí)間,不包括數(shù)據(jù)傳輸時(shí)間。同時(shí),為了消除偶然性誤差導(dǎo)致的時(shí)間不準(zhǔn)確等問題,統(tǒng)計(jì)結(jié)果為矩陣執(zhí)行10 000次SpMV的時(shí)間,且每次測(cè)試都迭代50次,取平均值作為最終實(shí)驗(yàn)結(jié)果。

在雙精度下對(duì)不同存儲(chǔ)格式的SpMV執(zhí)行時(shí)間進(jìn)行對(duì)比,結(jié)果如圖13所示。其中,DIA和HDC格式不進(jìn)行分塊,HDIA和DRM格式利用固定閾值進(jìn)行分塊,DIA-Adaptive根據(jù)矩陣非零元分布隨機(jī)分塊。從圖13中可以看出,在所有測(cè)試矩陣中,DRM格式均優(yōu)于DIA和DIA-Adaptive格式。在除shipsec5外的其余測(cè)試矩陣中,DRM格式優(yōu)于HDC格式。其原因在于,矩陣shipsec5主對(duì)角線中全為非零元素,主對(duì)角線周圍存在大量散點(diǎn),在采用HDC存儲(chǔ)時(shí),DIA只存儲(chǔ)主對(duì)角線,CSR存儲(chǔ)剩余非零元,整體零元素填充個(gè)數(shù)為0,減少了計(jì)算資源的浪費(fèi),縮短了SpMV執(zhí)行時(shí)間。因此,DRM格式略慢于HDC格式。在除nemeth22外的其余測(cè)試矩陣中,DRM格式優(yōu)于HDIA格式。其原因在于,矩陣nemeth22中非零元密集分布在主對(duì)角線周圍,使用HDIA格式分塊存儲(chǔ)時(shí),子塊間操作數(shù)負(fù)載更加平衡(如圖11所示),因此其SpMV執(zhí)行時(shí)間優(yōu)于DRM格式的。

Figure 13 Execution time comparison of different storage formats under double precision圖13 雙精度下不同存儲(chǔ)格式的執(zhí)行時(shí)間對(duì)比

表2統(tǒng)計(jì)了DRM格式在雙精度下的SpMV執(zhí)行時(shí)間相比于其他格式的最小、最大和平均加速比?？梢钥闯?DRM格式與DIA格式相比有很大的性能提升,獲得了最高155.40倍、平均20.76倍的加速。最高加速比發(fā)生在矩陣add32中,這是因?yàn)樵谠摼仃囍蟹橇阍貨]有按照對(duì)角線形式排列,矩陣存在大量散點(diǎn),DIA格式填充零元素的數(shù)量急劇增加(如圖10所示),導(dǎo)致性能降低。此外,和HDC、HDIA和DIA-Adaptive格式相比,DRM格式也取得了不錯(cuò)的性能提升,分別獲得了平均1.94,1.13和2.26倍,最高3.01,1.47和4.07倍的加速。

Table 2 Speedups of SpMV execution time under double precision表2 雙精度下SpMV執(zhí)行時(shí)間加速比

除了和以上格式對(duì)比外,本文還在雙精度下將DRM與CSR格式進(jìn)行了對(duì)比。實(shí)驗(yàn)結(jié)果表明,DRM在20個(gè)稀疏矩陣中的執(zhí)行時(shí)間都小于CSR格式的,取得了平均2.46倍,最高4.82倍的加速。因此,針對(duì)更加一般的稀疏對(duì)角矩陣,DRM相比于非基于DIA構(gòu)造的格式(如CSR等),仍能取得較好的性能加速。

圖14展示了單精度下不同存儲(chǔ)格式的SpMV執(zhí)行時(shí)間?？梢钥闯?相比于雙精度,單精度下所有格式在20個(gè)稀疏矩陣上的SpMV執(zhí)行時(shí)間均有所下降。在除shipsec5外的所有測(cè)試矩陣中,DRM格式相比于DIA、HDC、HDIA和DIA- Adaptive格式,同樣獲得了顯著的性能增益,執(zhí)行時(shí)間最短。在矩陣shipsec5中,DRM的執(zhí)行時(shí)間雖然慢于HDIA和HDC格式的,但也取得了與之相近的性能。

Figure 14 Execution time comparison of different storage formats under single precision圖14 單精度下不同存儲(chǔ)格式的執(zhí)行時(shí)間對(duì)比

表3統(tǒng)計(jì)了DRM格式在單精度下的SpMV執(zhí)行時(shí)間相比于其他格式的最小、最大和平均加速比?？梢钥闯?相比于雙精度,DRM格式相對(duì)DIA、HDC和DIA-Adaptive格式的加速都有所提升,分別取得了平均23.94,2.30和2.36倍,最高188.76,6.54和5.86倍的加速。對(duì)HDIA格式的加速比減少了0.01,但仍取得了平均1.12倍(最高1.53倍)的加速。

Table 3 Speedups of SpMV execution time under single precision表3 單精度下SpMV執(zhí)行時(shí)間加速比

然后對(duì)浮點(diǎn)運(yùn)算性能進(jìn)行對(duì)比,采用GFLOPS(即每秒10億次的浮點(diǎn)運(yùn)算數(shù))作為評(píng)價(jià)指標(biāo),其計(jì)算如式(18)所示:

(18)

其中,OPE為矩陣中的操作數(shù)總數(shù),N為矩陣的維數(shù),T為SpMV的內(nèi)核執(zhí)行時(shí)間。

對(duì)雙精度下不同存儲(chǔ)格式的GFLOPS進(jìn)行對(duì)比,結(jié)果如圖15所示。從圖15中可以看出,浮點(diǎn)計(jì)算性能隨著矩陣規(guī)模的增大而增加(DRM格式的最小GFLOPS為8.2,出現(xiàn)在規(guī)模最小的矩陣dw2048中;最大GFLOPS為228.07,出現(xiàn)在規(guī)模最大的矩陣Goodwin_127中),這是因?yàn)樵诰仃囈?guī)模較小時(shí),GPU中有一部分SM沒有滿載運(yùn)行,而隨著矩陣規(guī)模的增大,SM的利用率越來越高,在數(shù)據(jù)上體現(xiàn)為計(jì)算性能越來越高。在所有測(cè)試矩陣中,DRM格式的GFLOPS優(yōu)于HDC和DIA-Adaptive格式的。在除nemeth22外的其余測(cè)試矩陣中,DRM格式的GFLOPS優(yōu)于HDIA格式的;在除lns_3937、ex11、dw2048、rw5151和shy161外的其余測(cè)試矩陣中,DRM格式的GFLOPS優(yōu)于DIA格式的。其原因在于,矩陣nemeth22中,HDIA格式的負(fù)載更加均衡,SpMV執(zhí)行時(shí)間更短(如圖13所示),因此其GFLOPS略高于DRM格式的。矩陣lns_3937、ex11、dw2048、rw5151和shy161中,盡管DRM的SpMV執(zhí)行時(shí)間優(yōu)于DIA格式的,但DIA格式的零填充量分別是DRM格式的87.70,107.30,5.50,112.60和292.80倍(如圖11所示),在SpMV運(yùn)算時(shí),零填充量也屬于操作數(shù),需要進(jìn)行乘法和加法操作,因此DIA取得了較高的GFLOPS。

Figure 15 GFLOPS comparison of different storage formats under double precision圖15 雙精度下不同存儲(chǔ)格式的GFLOPS對(duì)比

Figure 16 GFLOPS speedups of different formats under double precision圖16 雙精度下不同格式的GFLOPS加速比

為了更直觀地展示DRM格式的性能提升,本文統(tǒng)計(jì)了雙精度下DRM的GFLOPS相比于其他格式的加速比,如圖16所示?？梢钥闯?在4種格式中,相比于HDC格式,DRM獲得了最好的性能加速,提供了平均5.28倍,最大11.68倍的加速;相比于DIA、HDIA、DIA-Adaptive格式,DRM也取得了較好的性能加速,分別提供了平均1.54,1.13和1.94倍,最高3.57,1.47和2.93倍的加速。

對(duì)單精度下不同存儲(chǔ)格式的GFLOPS進(jìn)行對(duì)比,結(jié)果如圖17所示?？梢钥闯?相比于雙精度,單精度的整體浮點(diǎn)計(jì)算性能都有所提高。同時(shí),和雙精度類似,DRM格式在大部分測(cè)試矩陣中的GFLOPS優(yōu)于其余4種格式的,最高GFLOPS達(dá)到了386.604(在矩陣Goodwin_095中),相比之下,DIA格式的最高性能GFLOPS不超過300,DIA-Adaptive格式的最高性能GFLOPS不超過240,而HDC格式的最高性能GFLOPS只有120左右?？傮w來說,與其他方法相比,本文的方法具有更好的浮點(diǎn)計(jì)算能力。

Figure 17 GFLOPS comparison of different storage formats under single precision圖17 單精度下不同存儲(chǔ)格式的GFLOPS對(duì)比

圖18展示了單精度下DRM格式的GFLOPS與其他格式的加速比。從圖18中可以看出,與雙精度相比,DRM對(duì)DIA、HDC和DIA-Adaptive的加速比都有所提升,分別取得了平均1.83,5.80和2.00倍,最高4.15,13.27和3.80倍的加速。HDIA格式的加速比略有下降,其主要原因在于,HDIA格式在矩陣shipsec5和Goodwin_127上的GFLOPS在單雙精度下差異較大,單精度相比雙精度,分別提升了75.9%和77.6%,而DRM格式只提升了65.7%和62.5%,因此,DRM對(duì)HDIA的GFLOPS加速比略有降低,但仍然取得了平均1.12倍,最高1.53倍的加速。

Figure 18 GFLOPS speedups of different formats under single precision圖18 單精度下不同格式的GFLOPS加速比

5.4 預(yù)處理時(shí)間分析

Figure 19 Preprocessing time for the 14 test matrices圖19 14個(gè)測(cè)試矩陣的預(yù)處理時(shí)間

本節(jié)對(duì)所提算法的預(yù)處理時(shí)間進(jìn)行分析。預(yù)處理階段包括矩陣劃分和矩陣重構(gòu)2個(gè)部分,其中,時(shí)間成本主要來自矩陣重構(gòu)部分。如算法1所示,重構(gòu)過程中會(huì)出現(xiàn)多次數(shù)據(jù)重排和歸并操作,相比于矩陣劃分部分,會(huì)帶來較大的時(shí)間開銷。本文分別計(jì)算14個(gè)測(cè)試矩陣的預(yù)處理時(shí)間,如圖19所示。從圖19中可以看出,所有測(cè)試矩陣的預(yù)處理時(shí)間開銷平均為8.41 ms,最大為38.32 ms(Kuu),最小為0.045 ms(dw2048)。其中,對(duì)于維數(shù)9 509,包含684 169個(gè)非零元的矩陣nemeth22,只需花費(fèi)0.37 ms就能完成矩陣的預(yù)處理過程。由此可以看出,所有測(cè)試矩陣的預(yù)處理時(shí)間開銷都能控制在毫秒范圍內(nèi)。

在有限元算法中,SpMV操作是求解線性系統(tǒng)的主要時(shí)間開銷[26]。而在求解過程中,通常需要多次迭代才能將結(jié)果收斂。如在解決大型稀疏線性系統(tǒng)的迭代方法(如共軛梯度和廣義最小殘差方法)中,求解器反復(fù)調(diào)用SpMV,通常需要數(shù)百到數(shù)千次迭代才能收斂[27]。在使用SPAI預(yù)條件的迭代求解器中,SpMV操作是最耗時(shí)的部分[28,29]。因此,在實(shí)際科學(xué)計(jì)算應(yīng)用中,這種預(yù)處理成本是可以接受的。

6 結(jié)束語

針對(duì)現(xiàn)有存儲(chǔ)格式對(duì)稀疏對(duì)角矩陣進(jìn)行存儲(chǔ)時(shí)出現(xiàn)的零填充和負(fù)載不平衡問題,本文提出了一種新的存儲(chǔ)格式DRM,并實(shí)現(xiàn)了基于此格式的并行SpMV算法。DRM利用基于固定閾值的矩陣劃分策略和基于迭代歸并的矩陣重構(gòu)策略,將輸入矩陣劃分為大小相等的行段,并按照行段中的操作數(shù)個(gè)數(shù)進(jìn)行重新排序,將滿足歸并要求的行段歸并為子塊,實(shí)現(xiàn)了子塊間操作數(shù)個(gè)數(shù)的最大相近。在NVIDIA?Tesla?V100平臺(tái)上的實(shí)驗(yàn)結(jié)果表明,在幾乎所有的測(cè)試矩陣中,DRM有效減少了零元素的填充量,實(shí)現(xiàn)了子塊間操作數(shù)個(gè)數(shù)方差的最小化。與當(dāng)前主流的DIA、HDC、HDIA和DIA-Adaptive格式相比,DRM格式在SpMV執(zhí)行時(shí)間上分別取得了20.76,1.94,1.10和2.26倍的加速,在浮點(diǎn)計(jì)算性能方面,分別取得了1.54,5.28,1.13和1.90倍的加速。但是,本文研究目前只適用于稀疏對(duì)角矩陣,針對(duì)部分非零元分布不規(guī)則的矩陣,仍存在零元素填充過多的現(xiàn)象。未來,將對(duì)DRM的劃分策略和非零元存儲(chǔ)方式進(jìn)行研究,以進(jìn)一步減少零元素的填充量,提升處理不同類型稀疏矩陣的能力,并將其應(yīng)用于相關(guān)領(lǐng)域的實(shí)際問題中。