申威26010眾核處理器上一維FFT實(shí)現(xiàn)與優(yōu)化?

趙玉文,敖玉龍,楊 超1,,劉芳芳,尹萬(wàn)旺，林蓉芬

1(中國(guó)科學(xué)院 軟件研究所 并行軟件與計(jì)算科學(xué)實(shí)驗(yàn)室,北京 100190)

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

3(計(jì)算機(jī)科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室(中國(guó)科學(xué)院 軟件研究所),北京 100190)

4(中國(guó)科學(xué)院大學(xué),北京 100049)

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

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

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

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

1 申威26010 處理器

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

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

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

2 相關(guān)工作

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

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

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

3 算法背景

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

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

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

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

Cooley-Tukey FFT算法.對(duì)于輸入長(zhǎng)度為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)先映射到一個(gè)N2×N1的二維矩陣形式.因此,N點(diǎn)一維FFT 問題的計(jì)算可以分解成N1點(diǎn)和N2點(diǎn)的小規(guī)模一維FFT 問題來(lái)完成,具體計(jì)算步驟如圖2 所示.(1)N2個(gè)N1點(diǎn)一維FFT 計(jì)算;(2)每個(gè)元素乘旋轉(zhuǎn)因子;(3)N1個(gè)N2點(diǎn)一維FFT 計(jì)算;(4)二維數(shù)組轉(zhuǎn)置,得到正確且有序的計(jì)算結(jié)果.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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的取值一般與一行/列從核的個(gè)數(shù)相同,紅框部分表示一行/列從核的計(jì)算任務(wù).首先計(jì)算f1點(diǎn)一維FFT,將計(jì)算任務(wù)沿N1的方向進(jìn)行分塊,其任務(wù)的個(gè)數(shù)為TN=16,任務(wù)大小為V*N2.每個(gè)任務(wù)按行優(yōu)先方式讀取數(shù)據(jù)到單個(gè)從核LDM 中;然后計(jì)算f2點(diǎn)一維FFT,因計(jì)算結(jié)果需轉(zhuǎn)置后存回,計(jì)算前需要對(duì)數(shù)據(jù)進(jìn)行重排和轉(zhuǎn)置,以達(dá)到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ī)模決策

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

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

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

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

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

實(shí)現(xiàn)時(shí)每個(gè)從核中的數(shù)據(jù)可以看作是一個(gè)m0×m1的矩陣,首先每個(gè)從核將本從核中已有的數(shù)據(jù)放到指定的位置中,然后使用PUT 指令向同行/列中其他7 個(gè)從核發(fā)送數(shù)據(jù),同時(shí)使用GET 指令接收其他從核發(fā)來(lái)的數(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ù)傳輸與從核上的計(jì)算可異步執(zhí)行,即DMA 訪存與計(jì)算重疊.本文算法在實(shí)現(xiàn)時(shí)利用了此特性,以提升性能.在從核函數(shù)中,其訪存與計(jì)算過(guò)程主要包括4 個(gè)操作,分別為從主存讀取數(shù)據(jù)(DMA read)、計(jì)算操作(computation)、寄存器通信(register communication)和將數(shù)據(jù)寫回主存(DMA write).對(duì)于分解N1=f1×f2,首先使用DMA Read 將數(shù)據(jù)從主存讀取到從核LDM 中,計(jì)算f1點(diǎn)FFT,寄存器通信完成同行/列從核的數(shù)據(jù)交換,計(jì)算f2點(diǎn)FFT,最后將計(jì)算結(jié)果寫回主存.當(dāng)寄存器通信完成之后從核就可以從主存中預(yù)取數(shù)據(jù),從而可以隱藏一部分計(jì)算操作.另外,對(duì)于DMA Write 操作,當(dāng)讀取操作完成時(shí),就可以進(jìn)行下一個(gè)循環(huán)的計(jì)算操作,而不需要等到DMA Write 操作完成后再執(zhí)行,也可隱藏一部分計(jì)算操作.具體如圖8 所示.

5 實(shí)驗(yàn)結(jié)果和性能分析

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

5.1 總體性能

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

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

5.2 訪存帶寬

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

對(duì)于第2 層分解,Nr=f1×f2[×f3](1≤r≤m),tw=max{N1,N2,...,Nm}.圖10 給出了各數(shù)據(jù)規(guī)模的訪存帶寬以及與帶寬利用率(總帶寬按照實(shí)測(cè)帶寬27.5Gb/s 計(jì)算),最高可達(dá)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 計(jì)算的時(shí)間復(fù)雜度由O(N2)降到了O(Nlog2N),FFT 性能通常由如下公式計(jì)算獲得:

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

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

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

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

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

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

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

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

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