趙玉文,敖玉龍,楊 超1,,劉芳芳,尹萬旺，林蓉芬

1(中國科學院 軟件研究所 并行軟件與計算科學實驗室,北京 100190)

2(北京大學 數(shù)學科學學院,北京 100871)

3(計算機科學國家重點實驗室(中國科學院 軟件研究所),北京 100190)

4(中國科學院大學,北京 100049)

5(國家并行計算機工程技術(shù)研究中心,北京 100190)

快速傅里葉變換(fast Fourier transform,簡稱FFT)算法被列為“20 世紀十大算法[1]”之一,在數(shù)字信號處理、圖像處理、微分方程求解和分子動力學等很多學科和科學計算領(lǐng)域具有重要地位.另外,FFT 作為 HPC Challenge 基準測試程序的組成部分[2],可以用于評估超級計算機的體系結(jié)構(gòu)和整體性能.

目前,超級計算機的體系結(jié)構(gòu)已經(jīng)從多核向眾核乃至異構(gòu)眾核發(fā)展,“神威·太湖之光”[3]是世界上首臺峰值運算速度超過10 億億次量級、國內(nèi)首臺采用國產(chǎn)處理器構(gòu)建的世界第一的超級計算機.截至2017 年11 月,已連續(xù)4 次榮登全球超級計算機500 強榜首.目前部署在國家超級計算無錫中心,使用的處理器是國產(chǎn)“申威26010”異構(gòu)眾核處理器.該處理器不同于現(xiàn)有的純CPU、CPU-MIC、CPU-GPU 等架構(gòu),其采用主-從核架構(gòu),單處理器峰值計算能力為3TFlops/s,訪存帶寬為130GB/s,相比計算能力,其訪存能力偏弱.然而,FFT 算法具有計算密集型和訪存密集型的特點,需要根據(jù)處理器體系結(jié)構(gòu)的特點設(shè)計出適用于此處理器的并行算法和優(yōu)化策略,這就給FFT 的高效實現(xiàn)帶來了巨大挑戰(zhàn).

本文根據(jù)申威26010 處理器體系結(jié)構(gòu)的特點,提出了基于兩層分解的一維FFT 眾核并行算法.該算法基于迭代的Stockham FFT 計算框架,根據(jù)Cooley-Tukey FFT 算法進行分解,將大規(guī)模FFT 分解成一系列的小規(guī)模FFT 來計算.其通過構(gòu)造合適的任務(wù)劃分方式,來滿足并行度要求和解決負載均衡問題,并采用寄存器通信、訪存計算重疊以及向量化等優(yōu)化方法進行優(yōu)化,以有效提高FFT 的計算性能.

1 申威26010 處理器

國產(chǎn)申威26010 處理器[3]采用異構(gòu)眾核架構(gòu),每個芯片由4 個核組(core group,簡稱CG)組成,每個核組由管理核心(management processing element,簡稱MPE)、運算核心陣列(computing processing elements cluster,簡稱CPE cluster)、協(xié)議處理部件(PPU)和存儲控制器(memory controller,簡稱MC)組成,核組間由片上網(wǎng)絡(luò)(Noc)連接,如圖1 所示.

Fig.1 The general architecture of the Sunway 26010 processor圖1 申威26010 處理器的總體結(jié)構(gòu)

管理核心又稱為主核,用于運行操作系統(tǒng)和用戶程序,其采用自主指令集的64 位全功能RISC 核心,實現(xiàn)了SIMD 擴展指令集.運算核心陣列采用緊耦合的方式進行管理,包含8×8 方式排列的64 個運算核心和DMA 控制器(DMA controlller).每個運算核心之間采用拓撲為8×8 Mesh 的簇通信網(wǎng)絡(luò)連接.運算核心又稱為從核,為精簡的64 位RISC 核心,可提供強大的計算能力,支持256 位整數(shù)和浮點向量化操作.DMA 控制器用于解析和管理來自運算核心的數(shù)據(jù)流傳輸命令,其可實現(xiàn)主存與LDM 間的快速數(shù)據(jù)傳輸.存儲層次上,主核采用一級數(shù)據(jù)和指令cache 分離、二級指令數(shù)據(jù)共享的兩級片上存儲層次,其中,L1 指令cache 為32KB,數(shù)據(jù)cache 為32KB,L2的cache 容量為256KB.每個從核采用了私有的16KB 的L1 指令Cache 和64KB 的SPM(scratch pad memory).SPM 可看作用戶可控的局部數(shù)據(jù)存儲LDM(local data memory),并利用DMA(direct memory access)控制器實現(xiàn)解析和管理主存與LDM 間的數(shù)據(jù)傳輸.此外,從核還提供高效的寄存器通信機制,能夠使核組內(nèi)同行/列的從核進行數(shù)據(jù)交換,能夠讓每個從核向它的同行/列的其他從核發(fā)送(源方)或者接收(目標方)數(shù)據(jù),源方和目標方利用生產(chǎn)者-消費者模式完成數(shù)據(jù)交換.寄存器通信由PUT 和GET 指令對組成,源方使用PUT 指令將通用寄存器中的一個向量長度的數(shù)據(jù)經(jīng)過從核通信網(wǎng)絡(luò),送入一個或者多個目標方的接收緩沖中.目標方使用GET 指令從接收緩沖中讀取數(shù)據(jù)并送入本地通用寄存器.并行編程模型上,對于一個核組,應(yīng)用程序由主核啟動,借助高性能線程庫Athread 將計算任務(wù)異步加載到從核執(zhí)行,雙方通過同步接口協(xié)同.由于申威26010 處理器復雜的架構(gòu)特性,針對FFT 算法的特點,我們主要從LDM、DMA 和寄存器通信機制這3 個體系結(jié)構(gòu)特性來提高帶寬的利用率,從而獲得較好的性能.

2 相關(guān)工作

近年來,有很多研究工作致力于CPU、GPU 和MIC 等平臺上對FFT 算法進行并行優(yōu)化,并開發(fā)出多個優(yōu)秀的軟件包.CPU 上比較著名的軟件包主要包括FFTW(fastest Fourier transform in the west)[4,5]、UHFFT[6]、SPIRAL[7]、ESSL 和MKL 等.FFTW 由MIT 的Frigo 和Johnson 開發(fā),是目前使用最為廣泛、運算速度較快的離散傅立葉變換計算庫.其不是采用固定算法計算變換,而是在運行時根據(jù)具體硬件和變換參數(shù)自動地選擇最佳的分解方案和不同算法,以期達到最佳效果.本文在測試算法的整體性能時會與FFTW 庫的性能進行對比.GPU 上主要包括CUFFT、GPUFFTW、AccFFT、FFTE 和P3DFFT[8]等.其中,CUFFT 是由NVIDIA 提供的基于GPU 的快速傅里葉變換函數(shù)庫,它基于Cooley-Tukey FFT 和Bluestein 算法并采用分治的算法.GPUFFTW 項目由美國北卡羅萊納大學提出,使用Stockham 的Autosort 方法實現(xiàn)了基于GPU 的一維傅立葉變換,運算速度最高可達29GFlops.AccFFT 采用兩種數(shù)據(jù)分解方式并基于CUFFT 和FFTW 庫實現(xiàn)分布式FFT 計算,支持CPU和GPU.FFTE 是HPC challenge benchmark 的組成部分,主要利用了cache 分塊技術(shù)和Stockham 的自動排序技術(shù).MIC 上主要是Intel 公司隨MIC 協(xié)處理器發(fā)布的MKL 數(shù)學庫.

目前,對FFT 算法的優(yōu)化策略主要包括兩種[9]:(1)數(shù)據(jù)分塊.利用分塊技術(shù)對數(shù)據(jù)進行處理,將子塊的數(shù)據(jù)傳輸?shù)矫總€處理器(線程),一般針對三維FFT 計算;文獻[9]、AccFFT、PFFT[10]、P3DFFT[8]、Takahashi[11]、2DECOMP & FFT 和Song[12]等策略均利用數(shù)據(jù)分塊策略進行并行優(yōu)化,其中,AccFFT、PFFT、2DECOMP & FFT和Song 考慮了通信優(yōu)化問題.另外,文獻[13,14]研究了天河1A 大型GPU 異構(gòu)集群系統(tǒng)上基于Parray 語言的數(shù)組維度類型程序設(shè)計方法.(2)算法分解.利用算法(如Cooley-Tukey FFT)將一個較大規(guī)模的一維FFT 分解成一系列較小規(guī)模的一維FFT,每個處理器(線程)并行地完成較小規(guī)模的FFT.文獻[15,16]嘗試著在IBM Cylops 64上使用Cooley-Tukey FFT 算法優(yōu)化FFT 計算,并且提出了一個針對該框架的性能模型.文獻[17]是在GPU 上基于Stockham FFT 框架和Cooley-Tukey FFT 算法實現(xiàn),并在文獻[18]中提出了一個用于生成GPU 上優(yōu)化的FFT內(nèi)核的自動優(yōu)化框架.文獻[19]提出了一種基于Cooley-Tukey FFT 算法和一套高效小規(guī)模codelet 代碼的計算多維FFT 框架.文獻[20]在單線程MKL 的基礎(chǔ)上,利用Cooley-Tukey FFT 算法和Intel Cilk Plus 的計算框架實現(xiàn)多線程并行.Nukada 等人也在文獻[21-24]中對3D FFT 算法在GPU 架構(gòu)上進行了一系列的研究.文獻[25]在單MIC 節(jié)點上基于Cooley-Tukey FFT 算法利用兩種分解算法和兩層并行方法實現(xiàn)了訪存高效的兩遍3D FFT算法.文獻[26]在Intel Xeon Phi 協(xié)處理器集群上實現(xiàn)了分布式1D FFT 的計算并進行了訪存與計算的重疊優(yōu)化.另外,文獻[27]主要對FFT 的性能進行了分析,也有一些文獻[28-32]對稀疏型的FFT 進行了研究,文獻[33]采用快速多極子方法來減少通信需求,文獻[34]在FFT 算法的容錯方面進行了研究.

在“神威·太湖之光”上已完成的各類運算任務(wù)中,其優(yōu)化主要從LDM、DMA、向量化、計算訪存重疊和寄存器通信等方面進行設(shè)計,主要包括:楊超團隊[35]的大氣動力模擬應(yīng)用項目,其獲得了國際高性能計算應(yīng)用領(lǐng)域的最高獎“戈登貝爾獎”;付昊桓團隊完成的“基于神威太湖之光的非線性地震模擬”[36],其獲得2017 年的“戈登貝爾獎”;文獻[37]主要是對stencil 計算進行了優(yōu)化,作者提出了一種求解歐拉大氣動力學方程的有效數(shù)據(jù)流方法[38];文獻[39]主要是對CAM(community atmosphere model)進行了重構(gòu)和優(yōu)化.

3 算法背景

離散傅里葉變換(discrete Fourier transform,簡稱DFT)是傅里葉變換在時域和頻域上都呈離散的形式,將信號的時域采樣變換為其離散時間傅里葉變換的頻域采樣.使用離散傅里葉變換可以將自然科學和工程技術(shù)中連續(xù)而復雜的問題轉(zhuǎn)化為離散且簡單的運算.對于一維長度為N的輸入序列x=[x0,x1,...,xN-1],其DFT 的計算公式為

DFT 可以表示成矩陣向量乘的形式XT=FN xT,其中,

DFT 的算法復雜度為O(N2),通常,x(n)和X(k)是復數(shù),當已知時,計算X(k)(k=0,1,...,N-1)需要N2次復數(shù)乘法和N(N-1)次復數(shù)加法,共需要8N2次浮點運算.

快速傅里葉變換FFT 是離散傅里葉變換的快速算法,其中,最經(jīng)典的是1965 年由Cooley 和Tukey 發(fā)表的Cooley-Tukey FFT 算法[40],首次將離散傅里葉變換的計算復雜度從O(N2)減少到O(NlogN),該算法主要依據(jù)旋轉(zhuǎn)因子的對稱性和周期性,采用分而治之的策略遞歸地進行計算.

Cooley-Tukey FFT算法.對于輸入長度為N的序列,若N可以分解為N1和N2的乘積,即N=N1×N2,則可以把輸入數(shù)據(jù)按行優(yōu)先自然地映射成N1×N2二維矩陣的形式.根據(jù)索引映射,令n=(n1,n2)=n1?N2+n2,k=(k2,k1)=k1+k2?N1,則式(1)可表示為

即輸出數(shù)據(jù)按行優(yōu)先映射到一個N2×N1的二維矩陣形式.因此,N點一維FFT 問題的計算可以分解成N1點和N2點的小規(guī)模一維FFT 問題來完成,具體計算步驟如圖2 所示.(1)N2個N1點一維FFT 計算;(2)每個元素乘旋轉(zhuǎn)因子;(3)N1個N2點一維FFT 計算;(4)二維數(shù)組轉(zhuǎn)置,得到正確且有序的計算結(jié)果.

如果嚴格按照上述4 個步驟來執(zhí)行,總共需要對主存中數(shù)組讀寫4 次,即總訪存量為8N(不統(tǒng)計訪問旋轉(zhuǎn)因子所增加的訪存量).因而在實現(xiàn)中為了減少訪存量,通常把乘旋轉(zhuǎn)因子和全局數(shù)組轉(zhuǎn)置分別合并到N1點FFT計算或者N2點FFT 計算中,使訪存量減少到4N.而N1點和N2點一維FFT 的計算通常遞歸地應(yīng)用Cooley-Tukey FFT 繼續(xù)分解,直到規(guī)模足夠小可以直接計算為止.

Fig.2 Flowchart of the Cooley-Tukey FFT algorithm圖2 Cooley-Tukey FFT 算法流程圖

對于Cooley-tukey FFT 算法中最后一步的轉(zhuǎn)置存回操作(對應(yīng)于基-2 Cooley-Tukey FFT 算法中的位倒序bit-reversal),其涉及到對輸入數(shù)據(jù)的全局轉(zhuǎn)置過程,雖然這個操作只需要O(n)時間,但仍然不可忽略,特別是對于變換規(guī)模較大的情況.我們利用Stockham FFT 框架進行優(yōu)化,避免了轉(zhuǎn)置操作.

Stockham FFT 計算框架.基于分解N=N1×...×Nr×...×Nm(1≤r≤m)使用迭代的方法將規(guī)模為N的一維FFT 轉(zhuǎn)化成一系列規(guī)模為Nr的一維FFT 計算.其算法流程如圖3 所示.

Fig.3 The flow chart of Stockham FFT algorithm圖3 Stockham FFT 算法流程圖

其通過在每步分解計算中穿插多維轉(zhuǎn)置,避免了一個顯式的全局轉(zhuǎn)置的過程.具體地,Nr點一維FFT 計算前后,數(shù)據(jù)在主存中的組織和維度為

其中,Nr表示當前需要計算的FFT 規(guī)模,P表示已經(jīng)完成的FFT 計算部分,T表示未計算的部分.舉例來說,若N可以分解為3 個因子的乘積,即N=N1×N2×N3,在計算每個Ni(1≤i≤3)點一維FFT 時,需要對整個數(shù)據(jù)讀寫遍歷1 遍,3 次遍歷數(shù)組中數(shù)據(jù)維度的變化過程為

4 基于兩層分解的一維FFT 眾核并行算法的設(shè)計與實現(xiàn)

FFT 算法具有計算密集型和訪存密集型的特點,其并行度有限、訪存局部性低,在26010 眾核平臺上很難充分利用眾多的計算資源.依據(jù)計算平臺的結(jié)構(gòu)特性和存儲層次特點,本文提出了基于兩層分解的一維FFT 眾核并行算法.該算法基于迭代的Stockham FFT 計算框架,將大規(guī)模FFT 分解成一系列的小規(guī)模FFT 來計算,分解規(guī)則基于Cooley-Tukey FFT 算法.通過設(shè)計合理的任務(wù)劃分方案、數(shù)據(jù)重用機制和計算與訪存重疊的雙緩沖優(yōu)化來有效地解決并行度和訪存量等問題.本文主要致力于雙精度復數(shù)到復數(shù)的2 的冪次一維FFT 計算,其他情況不作為本文研究的重點.

4.1 并行方案設(shè)計和任務(wù)劃分

基于兩層分解的一維FFT 眾核并行算法首先利用Cooley-Tukey FFT 算法對輸入數(shù)據(jù)規(guī)模為N(N=2n,n≥9)的一維FFT 進行第1 層分解,N=N1×...×Nr×...×Nm(1≤r≤m),即將原N點的一維FFT 轉(zhuǎn)化成m個相互獨立的計算步驟,輸入數(shù)據(jù)按行優(yōu)先自然地映射成m維矩陣的形式,其中,第Nm維度的數(shù)據(jù)連續(xù)存儲,其他維度非連續(xù)存儲.所有的分解因子Nr(1≤r≤m)滿足其規(guī)模的相應(yīng)數(shù)據(jù)可以完全加載到從核整個核組的LDM 中.如果Nr的規(guī)模足夠小,則不對Nr進行第2 次分解,直接計算Nr點FFT;如果Nr的規(guī)模較大,從高效的LDM 空間利用、DMA 傳輸方式和負載均衡3 個方面考慮,一般需要多個從核共同完成Nr點一維FFT 計算,此時利用Cooley-Tukey FFT 算法遞歸地對Nr進行第2 層分解Nr=f1×f2[×f3](f3可選),將數(shù)據(jù)平均分布到8 個或者64 個從核上.算法中,不進行第2 層分解的Nr點和fk(k=1,2,3)點一維FFT 直接調(diào)用底層高效的計算小因子FFT 的Codelets 庫完成.Codelets FFT 直接使用平臺特有的向量化指令編寫,并進行循環(huán)展開、仔細規(guī)劃寄存器資源和指令流水等優(yōu)化,使代碼更高效且具有更少的浮點運算量.另外,算法利用Stockham FFT 計算框架的自動排序功能來優(yōu)化Cooley-Tukey FFT 算法中將結(jié)果轉(zhuǎn)置存回的問題.

對于算法中m個相互獨立的計算步驟,每一步都需要利用DMA 將相應(yīng)數(shù)據(jù)傳輸?shù)胶私MLDM 中,并在核組內(nèi)完成Nr(1≤r≤m)點一維FFT 的計算.根據(jù)數(shù)據(jù)存儲是否連續(xù),可將m個計算步驟分成兩類.

(1)數(shù)據(jù)非連續(xù)存儲Nr(1≤r≤m),其分解可簡寫成P×Nr×T(1＜r＜m),P為已經(jīng)完成FFT 計算的部分(對于第N1維度,P=1),T為未計算的部分,此時將計算任務(wù)沿T的方向進行分塊,每個任務(wù)塊的大小為V*Nr,任務(wù)塊個數(shù)為(P*T)/V,其可以滿足DMA 傳輸數(shù)據(jù)的連續(xù)度要求,從而使DMA 傳輸更高效.但是由于單個從核的LDM 空間有限,能夠獨立處理的任務(wù)塊一般較小,需要利用多個從核協(xié)作以完成一個任務(wù)塊,此時需要對Nr進行第2 層分解Nr=f1×f2[×f3].如果Nr分解成Nr=f1×f2,則每個任務(wù)塊由一行/列從核完成;如果Nr分解成Nr=f1×f2×f3,則每個任務(wù)塊由整個核組共同完成.

(2)數(shù)據(jù)連續(xù)存儲Nm(r=m),分解可以簡寫成P×Nm,計算任務(wù)沿P的方向進行分塊.當Nm≤M(M為每個從核所能計算的最大數(shù)據(jù)量)時,按行優(yōu)先方式讀取數(shù)據(jù),將Nm點一維FFT 計算所需的數(shù)據(jù)可以直接加載到一個從核的LDM 上,每個從核的計算量為V=M/Nm個Nm點FFT.但當V較小時,由于結(jié)果需要以轉(zhuǎn)置的形式存回到不連續(xù)的主存中,單個從核直接存回只能保證V個數(shù)的連續(xù)度,為了滿足高效的DMA 傳輸對傳輸數(shù)據(jù)的連續(xù)度要求,可以通過對一行/列中的8 個從核的數(shù)據(jù)進行重排和轉(zhuǎn)置,從而保證存回的數(shù)據(jù)塊有8*V的連續(xù)度,此時需要將Nm進一步分解成Nm=f1×f2,由一行/列從核共同完成.但當Nm＞M時,單個從核不能容納Nm規(guī)模的數(shù)據(jù),需要多個從核協(xié)作完成計算任務(wù).一般由一行/列從核完成V=(8*M)/Nm個Nm點一維FFT 計算,此時需要將Nm進一步分解成Nm=f1×f2×f3,其計算方式與Nm≤M類似.對于以上兩種情況中的數(shù)據(jù)重排操作一般借助寄存器通信機制完成,以進一步減少訪存量.

下面以4×4 的從核陣列來展示N=N1×N2的任務(wù)劃分方式,其中圖4 展示第N1維度,圖5 展示第N2維度,圖中×表示分解操作,*表示任務(wù)劃分中的乘法操作,Z方向表示數(shù)據(jù)連續(xù)存儲,Y和X分別次之.圖4 中,計算任務(wù)沿N2的方向進行分塊,每個任務(wù)塊大小為V*N1,任務(wù)塊個數(shù)為TN=N2/V=4,分別為T0、T1、T2 和T3,每個任務(wù)塊由一行/列從核完成.將N1分解為N1=f1×f2,首先計算f1點一維FFT,對每個任務(wù)塊從f2的方向進行劃分得到4 個小任務(wù)塊,并將每個小任務(wù)塊分配給一個從核;然后計算f2點一維FFT,其任務(wù)劃分方式與計算f1點FFT 過程類似,但其任務(wù)劃分沿f1方向.從核中f1和f2的計算直接調(diào)用Codelets FFT 來完成.

Fig.4 The diagram of two-layer decomposition 1-D FFT multi-core parallel algorithm (N1=f1×f2)圖4 基于兩層分解的一維FFT 眾核并行算法圖(N1=f1×f2)

圖5 中,N2≤M,N2=f1×f2,其中,f1的取值一般與一行/列從核的個數(shù)相同,紅框部分表示一行/列從核的計算任務(wù).首先計算f1點一維FFT,將計算任務(wù)沿N1的方向進行分塊,其任務(wù)的個數(shù)為TN=16,任務(wù)大小為V*N2.每個任務(wù)按行優(yōu)先方式讀取數(shù)據(jù)到單個從核LDM 中;然后計算f2點一維FFT,因計算結(jié)果需轉(zhuǎn)置后存回,計算前需要對數(shù)據(jù)進行重排和轉(zhuǎn)置,以達到DMA 數(shù)據(jù)傳輸中f1*V的連續(xù)度.

Fig.5 The diagram of two-layer decomposition 1-D FFT multi-core parallel algorithm (N2=f1×f2)圖5 基于兩層分解的一維FFT 眾核并行算法圖(N2=f1×f2)

4.2 規(guī)模決策

對于規(guī)模為N的一維FFT,按照基于兩層分解的一維FFT 眾核并行算法原理,其存在多種分解策略,但是并非所有的分解策略都能得到較好的計算性能.本文算法雖然具有很好的并行度,但其帶來的是訪存量的大幅增加.對于分解N=N1×N2×...×N m(m≥2),當總數(shù)據(jù)規(guī)模大于LDM 空間時,至少需要對主存數(shù)組讀寫m次,其訪存量為2mN.對于第2 層分解Nr=f1×f2[×f3](1≤r≤m),因遞歸地應(yīng)用Cooley-Tukey FFT,每一個分解又至少增加2N訪存量,因此在數(shù)據(jù)規(guī)模允許的情況下,m的取值應(yīng)盡量小,第2 層分解類似.

對于分解規(guī)模的選取,本算法采用down-up 型兩層策略樹的方式進行選取,其特點是首先考慮第2 層分解中Nr(1≤r≤m)的分解(down 方向),然后再考慮第1 層中N的分解(up 方向).算法中Nr的取值一般為64、128、256、512 和1024 等,目前可通過分析和實驗的方式確定性能最優(yōu)的Nr分解結(jié)果.已知第2 層的分解結(jié)果后,根據(jù)具體規(guī)模得到性能最優(yōu)的第1 層的分解.具體內(nèi)容如圖6 所示,該圖僅用部分實例來描述down-up 型兩層策略樹,并不代表最優(yōu)的分解策略.圖中連接線的數(shù)字代表分解因子的個數(shù),例如節(jié)點64 與節(jié)點8 之間的數(shù)字2 表示64 分解為64=8×8.

Fig.6 Down-up policy tree圖6 Down-up 型策略樹

4.3 數(shù)據(jù)重用機制

由第4.2 節(jié)可知,為了解決基于本文算法訪存量大幅增加的問題,還需考慮數(shù)據(jù)的重用問題.申威26010 處理器提供的寄存器通信機制可以實現(xiàn)核組內(nèi)同行/列從核之間的數(shù)據(jù)共享,能夠有效地減少利用主存進行數(shù)據(jù)交換和數(shù)據(jù)共享帶來的數(shù)據(jù)移動開銷.因此,本算法主要借助寄存器通信機制來避免對兩種訪存問題的開銷.(1)分解訪存開銷.對于第2 層分解Nr=f1×f2[×f3](1≤r≤m),在計算f1點一維FFT 時,數(shù)據(jù)已通過DMA 由主存?zhèn)鬏數(shù)綇暮薒DM 中,借助寄存器通信機制,使一行/列從核的數(shù)據(jù)進行交換,得到計算f2和f3點一維FFT 所需的數(shù)據(jù),以避免從主存讀取數(shù)據(jù),從而達到數(shù)據(jù)的重用.(2)數(shù)據(jù)重排.對于Nm點FFT 的計算,其結(jié)果需要以轉(zhuǎn)置的形式存回到不連續(xù)的主存中,如圖5 所示,故需借助寄存器通信機制完成對一行/列從核的數(shù)據(jù)重排,以保證DMA 寫回主存數(shù)據(jù)的連續(xù)度要求.

實現(xiàn)時每個從核中的數(shù)據(jù)可以看作是一個m0×m1的矩陣,首先每個從核將本從核中已有的數(shù)據(jù)放到指定的位置中,然后使用PUT 指令向同行/列中其他7 個從核發(fā)送數(shù)據(jù),同時使用GET 指令接收其他從核發(fā)來的數(shù)據(jù).如圖7 所示,圖中僅顯示了一行寄存器通信的內(nèi)容,且m0=f1=16,m1一般為算法中的V,f2=8,threadnum=8.

Fig.7 Row-wise register communication圖7 行寄存器通信

4.4 雙緩沖優(yōu)化

申威26010 處理器主從核之間的數(shù)據(jù)傳輸與從核上的計算可異步執(zhí)行,即DMA 訪存與計算重疊.本文算法在實現(xiàn)時利用了此特性,以提升性能.在從核函數(shù)中,其訪存與計算過程主要包括4 個操作,分別為從主存讀取數(shù)據(jù)(DMA read)、計算操作(computation)、寄存器通信(register communication)和將數(shù)據(jù)寫回主存(DMA write).對于分解N1=f1×f2,首先使用DMA Read 將數(shù)據(jù)從主存讀取到從核LDM 中,計算f1點FFT,寄存器通信完成同行/列從核的數(shù)據(jù)交換,計算f2點FFT,最后將計算結(jié)果寫回主存.當寄存器通信完成之后從核就可以從主存中預取數(shù)據(jù),從而可以隱藏一部分計算操作.另外,對于DMA Write 操作,當讀取操作完成時,就可以進行下一個循環(huán)的計算操作,而不需要等到DMA Write 操作完成后再執(zhí)行,也可隱藏一部分計算操作.具體如圖8 所示.

5 實驗結(jié)果和性能分析

本節(jié)利用申威26010 眾核的一個核組測試本文提出的算法的總體性能與訪存帶寬,并對其性能進行分析.具體信息見第1 節(jié).

5.1 總體性能

本節(jié)主要測試基于兩層分解的一維FFT 眾核并行算法的總體性能,并且與在申威26010 處理器單主核上運行的FFTW3.3.4 庫的性能進行對比.其中,Basic 是基于兩層分解的一維FFT 眾核并行算法的基礎(chǔ)版本,其沒有進行寄存器通信、SIMD 向量化和雙緩沖等優(yōu)化,Basic+Register 是在Basic 版本上進行寄存器通信優(yōu)化的版本,Basic+Register+Simd 是在 Basic+Register 版進行向量化優(yōu)化的版本,xMathFFT 為最終版本,其在 Basic+Register+Simd 基礎(chǔ)上進行了雙緩沖的優(yōu)化,其性能結(jié)果如圖9 所示.

總體來說,本文算法獲得了較好的性能,相比FFTW 獲得了平均44.53x 的加速比,最高加速比可達56.33x,最低加速比可達24.67x.另外,Basic 版獲得了平均13.76x 的加速比,Basic+Register 版獲得了平均24.94x 的加速比,Basic+Register+Simd 版獲得了平均33.52 的加速比.

5.2 訪存帶寬

因本文算法在申威26010 異構(gòu)眾核處理器上的訪存時間和計算時間(計算已采用向量化等技術(shù)充分優(yōu)化)相當,訪存時間略長,故采用帶寬利用率作為評價指標.本節(jié)測試了基于兩層分解的一維FFT 眾核并行算法的訪存帶寬情況.假設(shè)對于雙精度復數(shù)的一維FFT 計算,DMA 傳輸?shù)臄?shù)據(jù)量為td×N+tw(單位:雙精度復數(shù)),其中,N為一維FFT 的數(shù)據(jù)規(guī)模,td為整個計算過程中需要對主存中數(shù)據(jù)遍歷的次數(shù),tw表示旋轉(zhuǎn)因子的傳輸量.td和tw的結(jié)果均與本文算法中的分解結(jié)果有關(guān),若有分解N=N1×...×Nm,則有

對于第2 層分解,Nr=f1×f2[×f3](1≤r≤m),tw=max{N1,N2,...,Nm}.圖10 給出了各數(shù)據(jù)規(guī)模的訪存帶寬以及與帶寬利用率(總帶寬按照實測帶寬27.5Gb/s 計算),最高可達83.45%,最低為48.41%,平均帶寬利用率為68.78%.

Fig.10 The memory bandwidth of two-layer decomposition 1-D FFT multi-core parallel algorithm圖10 基于兩層分解的一維FFT 眾核并行算法的訪存帶寬

5.3 性能分析

FFT 算法將DFT 計算的時間復雜度由O(N2)降到了O(Nlog2N),FFT 性能通常由如下公式計算獲得:

其中,N為一維FFT 變換的規(guī)模;runtime為計算花費的時間,以秒(s)為單位.在申威26010 眾核上,從核計算可以與DMA 傳輸異步進行.一般來說,運行時間由計算時間和DMA 傳輸時間構(gòu)成,但是,由于FFT 是一種典型的帶寬密集型算法,好的算法設(shè)計方案可以將從核計算與DMA 傳輸重疊起來,從而隱藏從核計算的開銷.因而,我們可以根據(jù)申威26010 眾核的DMA 傳輸帶寬參數(shù)來估計實現(xiàn)能達到的性能上限.

由第5.2 節(jié)可知,DMA 傳輸?shù)臄?shù)據(jù)量為td×N+tw,由于所需加載的旋轉(zhuǎn)因子數(shù)組的規(guī)模與變換規(guī)模N以及其分解方案存在很大的關(guān)系,因此只考慮第1 層分解中的旋轉(zhuǎn)因子傳輸部分,即取td=N.

理想情況下,雙精度復數(shù)到復數(shù)FFT 的性能上限可由如下公式來計算:

表2 根據(jù)上述設(shè)計方案,對典型一維雙精度復數(shù)到復數(shù)2 的冪次FFT 的性能上限進行了估計.

Table 2 Comparison of ideal performance with actual performance at various scales表2 各規(guī)模理想性能與實際性能對比

從表2 可知,本文算法獲得了較好的性能,相比理想的Gflops 性能,其實際性能占比平均為78.5%,最高占比可達95.4%,最低占比為65.5%.

6 總結(jié)和下一步工作

FFT 算法具有計算密集型和訪存密集型的特點,其并行度有限、訪存局部性低,在申威26010 眾核處理器上很難充分利用眾多的計算和訪存資源.本文首先對DFT、Cooley-Tukey FFT 算法和Stockham FFT 計算框架以及申威26010 眾核處理器等進行了較為詳盡的介紹.然后根據(jù)申威26010 眾核處理器的結(jié)構(gòu)特點和存儲層次特點,提出了基于兩層分解的一維FFT 眾核并行算法.該算法采用迭代的Stockham FFT 計算框架,將大規(guī)模FFT分解成一系列的小規(guī)模FFT 來計算,而分解規(guī)則基于Cooley-Tukey FFT 算法.該算法利用寄存器通信、訪存計算重疊以及SIMD 向量化等優(yōu)化手段來有效地提高整體性能.最后對本文算法的性能和帶寬進行了測試并對其性能進行了剖析,其性能相比于單主核上運行的FFTW3.3.4 庫,獲得了平均44.53x 的加速比,最高加速比可達56.33x,且其帶寬利用率最高可達83.45%.本文雖然以申威26010 平臺來開展研究,但本文提出的算法對其他異構(gòu)眾核系統(tǒng)架構(gòu),如GPU 也具有很大的借鑒意義.

下一步工作可以從以下兩方面展開:一方面,可以對本文的工作進行擴展,功能上實現(xiàn)支持非2 的冪次一維FFT 計算以及任意規(guī)模的多維FFT;另一方面,將本文提出的算法移植到其他計算平臺.