大點(diǎn)數(shù)FFT 在“申威26010”上的并行優(yōu)化

郭俊，劉鵬，楊昕遙，張魯飛，吳東

(1.湖州職業(yè)技術(shù)學(xué)院 信息工程與物聯(lián)網(wǎng)學(xué)院，浙江 湖州 313000；2.湖州職業(yè)技術(shù)學(xué)院 湖州市物聯(lián)網(wǎng)智能系統(tǒng)集成技術(shù)重點(diǎn)實(shí)驗(yàn)室，浙江 湖州 313000；3.浙江大學(xué) 信息與電子工程學(xué)院，浙江 杭州 310027；4.螞蟻科技集團(tuán)股份有限公司，浙江 杭州 310013；5.數(shù)學(xué)工程與先進(jìn)計(jì)算國家重點(diǎn)實(shí)驗(yàn)室，江蘇 無錫 214125)

隨著數(shù)字計(jì)算機(jī)和超大規(guī)模集成電路的飛速發(fā)展，數(shù)字信號處理技術(shù)得到了越來越廣泛的應(yīng)用，快速傅里葉變換（fast Fourier transform，F(xiàn)FT）是其中最基本、使用最頻繁的核心算法之一.其他諸多信號處理算法，例如卷積、濾波、頻譜分析等，都可以轉(zhuǎn)化為FFT 來實(shí)現(xiàn).

超級計(jì)算機(jī)為解決當(dāng)今海量級數(shù)據(jù)的科學(xué)與工程計(jì)算問題提供了良好的平臺.“神威·太湖之光”是世界上首臺峰值每秒浮點(diǎn)運(yùn)算次數(shù)（FLOPs）超過100×1015的超級計(jì)算機(jī)，也是中國第一臺全部采用自主技術(shù)構(gòu)建的排名世界第一的超級計(jì)算機(jī)[1].“神威·太湖之光”由國產(chǎn)“申威26010”處理器組成，單處理器峰值FLOPs 可達(dá)3.1×1012.相比之下，訪存帶寬僅有136.5 GB/s 則略顯遜色，計(jì)算訪存比接近22.7.在面對諸如FFT 這樣兼具計(jì)算密集型和訪存密集型特點(diǎn)的程序時(shí)，需要設(shè)計(jì)特定的并行實(shí)現(xiàn)策略，以優(yōu)化其性能.

本文的主要工作是根據(jù)“申威26010”處理器的架構(gòu)特點(diǎn)和編程規(guī)范，提出針對大點(diǎn)數(shù)FFT 的眾核優(yōu)化方案.該方案源自經(jīng)典的Cooley-Tukey FFT 算法，將一維數(shù)據(jù)迭代分解為二維FFT 進(jìn)行加速.為了解決“列FFT”的讀寫和計(jì)算問題，降低矩陣轉(zhuǎn)置的影響，提出“列均分-行連續(xù)”結(jié)合DMA跨步傳輸?shù)淖x寫策略.該方法還可以推廣至其他需要進(jìn)行大量列數(shù)據(jù)操作的應(yīng)用.計(jì)算部分采用旋轉(zhuǎn)因子優(yōu)化和向量化操作，傳輸部分采用雙緩沖策略和寄存器通信，能夠充分利用眾核系統(tǒng)的計(jì)算資源和傳輸帶寬，達(dá)到了良好的加速效果.

1 相關(guān)工作及研究現(xiàn)狀

1.1 快速傅里葉變換

FFT 是離散傅里葉變換（discrete Fourier transform，DFT）的快速算法，由Cooley 等[2]提出，經(jīng)過近60 年的發(fā)展，形成了龐大的算法體系，陸續(xù)衍生出包括高基[3]、混合基[4]、分裂基[5]、多維度[6]等多種FFT 算法，以應(yīng)對各種數(shù)據(jù)類型和數(shù)據(jù)量.

在具體實(shí)現(xiàn)方面，F(xiàn)FT 的硬件優(yōu)化主要針對FPGA、ASIC、DSP 等專用處理器，大致可以分為基于存儲器和基于流水線2 類.基于存儲器的結(jié)構(gòu)利用數(shù)據(jù)在計(jì)算單元和存儲單元之間往復(fù)運(yùn)算和循環(huán)路由的特征來降低硬件復(fù)雜度[7-8].基于流水線的結(jié)構(gòu)通過設(shè)計(jì)單路/多路的前饋/反饋路徑，提升蝶形單元的吞吐率[9-10].

FFT 的軟件優(yōu)化主要是通用平臺上高性能軟件包（庫）的開發(fā).常見的有CPU 平臺上Frigo 等[11]開發(fā)的FFTW、Intel 公司[12]的MKL FFTW 等、GPU平臺上NVIDIA 公司[13]的CUFFT、Ayala 等[14]提出的heFFTe 等.其中，F(xiàn)FTW 是目前使用最廣泛、綜合性能較好的FFT 庫.張明[15]在龍芯3B 處理器上、郭金鑫[16]在ARM V8 平臺上、操廬寧[17]在X86-64 平臺上均對FFTW 庫進(jìn)行了移植和優(yōu)化.

1.2 申威26010 與神威·太湖之光

國產(chǎn)“申威26010”處理器由上海高性能集成電路設(shè)計(jì)中心通過自主技術(shù)研制，采用片上計(jì)算陣列集群和分布式共享存儲相結(jié)合的異構(gòu)眾核體系架構(gòu)，使用64 位自主申威指令系統(tǒng).單芯片集成了4 個(gè)核組共260 個(gè)核心，每個(gè)核組包含1 個(gè)控制核心（主核）和由64 個(gè)運(yùn)算核心（從核）組成的核心陣列.如圖1（a）所示為處理器的結(jié)構(gòu)示意圖，如圖1（b）所示為從核與主存、局存、寄存器之間的簡化通信模型.在編程規(guī)范上，核組內(nèi)部采用Athread 多線程庫或OpenACC 并行編程，核組間使用MPI 或OpenMP 編程接口.新一代“申威26010-Pro”處理器于2021 年發(fā)布，單核組運(yùn)算性能和訪存帶寬分別提升了3.1 倍和1.8 倍.

圖1 “申威26010”處理器的示意圖Fig.1 Schematic of Sunway 26010 processor

“神威·太湖之光”由中國國家并行計(jì)算機(jī)工程技術(shù)研究中心于2015 年12 月完成研制，投入運(yùn)行，已成功解決眾多的科學(xué)工程計(jì)算問題.目前已升級為配備有“申威26010-Pro”處理器的新一代超算系統(tǒng).代表性應(yīng)用包括Yang 等[18]的“千萬核可擴(kuò)展大氣動力學(xué)全隱式模擬”、Fu 等[19]的“非線性地震模擬”、Liu 等[20]的“超大規(guī)模量子隨機(jī)電路實(shí)時(shí)模擬”等，分別于2016 年、2017 年、2021 年三度榮獲“戈登貝爾獎”.

在“神威·太湖之光”上開發(fā)了高性能擴(kuò)展數(shù)學(xué)庫xMath 并移植了FFTW 庫，其中面向科學(xué)計(jì)算BLAS 庫的優(yōu)化已十分成熟[21-22]，但是缺少針對FFT 的研究.趙玉文等[23]在“申威26010”上實(shí)現(xiàn)了迭代Stockham FFT 計(jì)算框架，與單主核FFTW相比，獲得了平均44.5x、最高56.3x 的加速比.Stockham FFT 是Cooley-Tukey FFT 的一種常見變化算法，主要特點(diǎn)是消除了輸入輸出的倒位排序，可以避免傳統(tǒng)FFT 算法中的顯式轉(zhuǎn)置，在大數(shù)據(jù)量下能夠?qū)崿F(xiàn)一定的性能提升，代價(jià)是失去了蝶形結(jié)構(gòu)的對稱性和原位計(jì)算的特征，迭代框架的實(shí)現(xiàn)相對比較復(fù)雜.

針對FFT 優(yōu)化的開發(fā)需求，根據(jù)“申威26010”處理器的架構(gòu)特點(diǎn)與FFT 算法的內(nèi)在聯(lián)系，提出“列均分-行連續(xù)”讀寫策略.結(jié)合其他優(yōu)化手段，對數(shù)據(jù)進(jìn)行巧妙的劃分、重排和交換，可以減少行列轉(zhuǎn)置帶來的影響，完成了經(jīng)典Cooley-Tukey FFT 的眾核實(shí)現(xiàn)，達(dá)到了大點(diǎn)數(shù)下平均60x 以上的加速比.

2 大點(diǎn)數(shù)FFT 算法

離散傅里葉變換DFT 的計(jì)算公式為

FFT 作為DFT 的快速算法，核心思想是將長輸入逐次分解成短序列DFT進(jìn)行遞歸計(jì)算，通過分而治之的策略，將計(jì)算復(fù)雜度降至O(Nlog2N).大點(diǎn)數(shù)FFT 的加速算法沿用了上述思想，通過將一維數(shù)據(jù)轉(zhuǎn)化為二維或更高維，減少了每一維上的數(shù)據(jù)量，更有利于進(jìn)行多核/眾核并行.

將輸入規(guī)模為N=N1×N2的一維數(shù)組映射為N1行、N2列的二維矩陣.令輸入、輸出序號為

其中，0 ≤n1、k1≤N1-1，0 ≤n2、k2≤N2-1.有如下關(guān)系：

一維轉(zhuǎn)二維FFT 的計(jì)算流程（見圖2）如下.

圖2 一維轉(zhuǎn)二維FFT 的計(jì)算流程Fig.2 Calculation process of 1D to 2D FFT

1）將一維數(shù)據(jù)按行優(yōu)先順序排列為N1行、N2列的二維矩陣.

2）進(jìn)行每一列N1點(diǎn) FFT，共進(jìn)行N2次，按原位置存回.

3）對于矩陣中的每一個(gè)數(shù)據(jù)，根據(jù)其坐標(biāo)，乘以旋轉(zhuǎn)因子 e xp(-j2πn2k1/N).

4）進(jìn)行每一行N2點(diǎn)FFT，共進(jìn)行N1次，按原位置存回.

5）按列優(yōu)先順序?qū)?shù)據(jù)讀出，得到N=N2N1的最終結(jié)果.

3 大點(diǎn)數(shù)FFT 在“申威26010”上的并行優(yōu)化

將大點(diǎn)數(shù)FFT 在“申威26010”上并行實(shí)現(xiàn)時(shí)，主要面臨如下難題.1）主核與從核之間列數(shù)據(jù)的讀寫問題.2）從核間行列數(shù)據(jù)的分布與轉(zhuǎn)置問題.3）從核內(nèi)小點(diǎn)數(shù)FFT 的高效實(shí)現(xiàn)問題.

下面給出具體的并行優(yōu)化解決方案.

3.1 基本存儲和分解策略

使用雙精度浮點(diǎn)復(fù)數(shù)，包括8 B 的實(shí)部和8 B 的虛部.“申威26010”從核對局存LDM 的訪問延遲為4 個(gè)CPU 周期，僅為對主存訪問延遲的1/40，應(yīng)盡可能將主存中的數(shù)據(jù)放入LDM 內(nèi)再計(jì)算.考慮到單從核LDM 大小僅為64 kB，采用雙緩沖策略需要2 份輸入/輸出緩沖區(qū)（3.6 節(jié)），將可單次寫入LDM 的最大數(shù)據(jù)量規(guī)定為512 個(gè)復(fù)數(shù).輸入和輸出共占512×16×2×2=32 kB，余下一半空間存儲中間結(jié)果和臨時(shí)變量.

“申威26010”單個(gè)核組陣列包括8 行8 列共64 個(gè)運(yùn)算核心.為了方便維度分解和FFT 計(jì)算，總數(shù)據(jù)量采用2 的正整數(shù)冪次，每一維上的最小尺寸設(shè)置為8.下面給出基本的分解策略.

1）將總點(diǎn)數(shù)分解為N=N1N2，通常N1和N2較大，需要繼續(xù)分解.以N1為例，可以按以下a）、b）、c）共3 種情況進(jìn)行分解.

a）N1＜64，即8、16、32 點(diǎn)，停止分解，直接在單核上進(jìn)行小點(diǎn)數(shù)FFT 計(jì)算.

b）64 ≤N1≤256，將N1進(jìn)行二維分解，包括8×8、16×8、32×8 共3 種方案.前一位數(shù)（8、16、32）代表單核上的FFT 點(diǎn)數(shù)，后一位數(shù)（8）代表啟動一行（或一列）8 個(gè)從核進(jìn)行計(jì)算.

c）512 ≤N1≤2 048，將N1進(jìn)行三維分解，包括8×8×8、16×8×8、32×8×8 共3 種方案.前一位數(shù)（8、16、32）代表單核上的FFT 點(diǎn)數(shù)，后兩位數(shù)（8×8）代表啟動64 個(gè)從核進(jìn)行計(jì)算.

2）同理，N2按以上3 條進(jìn)行再分解.當(dāng)N1、N2≥4 096 時(shí)，對總點(diǎn)數(shù)N進(jìn)行三維及以上分解，并啟動多個(gè)核組.

3）每多分解一次，需要將當(dāng)前FFT 替換為圖2所示的“列FFT-乘旋轉(zhuǎn)因子-行FFT”流程，直至迭代計(jì)算完畢.

下面以輸入數(shù)據(jù)量為262 144=512×512 的實(shí)例，闡述各創(chuàng)新實(shí)現(xiàn)技術(shù).

3.2 “列均分-行連續(xù)”策略

在行優(yōu)先的存儲系統(tǒng)中，對列數(shù)據(jù)的讀寫是離散非連續(xù)的，傳輸效率很低.一維轉(zhuǎn)二維FFT并行方案的主要難點(diǎn)是對矩陣“列FFT”的處理.大數(shù)據(jù)量加劇了FFT 訪存密集型特點(diǎn)帶來的上述影響，若在主核上先進(jìn)行顯式轉(zhuǎn)置再傳輸，則會增加不可忽略的額外開銷.

為了提升“列FFT”的并行效率，結(jié)合FFT 算法、處理器架構(gòu)以及數(shù)據(jù)在從核間的分布特點(diǎn)，提出“列均分-行連續(xù)”讀寫策略，即并非用單個(gè)從核去讀取主存中的列數(shù)據(jù)，而是把每一列的數(shù)據(jù)循環(huán)平分給64 個(gè)從核.如圖3 所示，對于“大矩陣”的每一列，0 號核讀第0 行，1 號核讀第1 行，······，63 號核讀第63 行.然后0 號核讀第64 行，······，直至63 號核讀第511 行.從列的方向來看，每一列循環(huán)平分到64 個(gè)從核內(nèi).每個(gè)核一次可以讀取多列，從行的方向來看，數(shù)據(jù)是連續(xù)讀取的，即為“列均分-行連續(xù)”策略.

圖3 “列均分-行連續(xù)”策略Fig.3 Column-sharing,row-continuity strategy

每個(gè)從核分到512/64=8 行的數(shù)據(jù).由于LDM最大單次輸入量為512，可以連續(xù)讀入“大矩陣”每一行中512/8=64 個(gè)數(shù).由于“申威26010”主存與局存之間的DMA 帶寬性能在傳輸粒度≥256 B（32 個(gè)double 數(shù)據(jù)）時(shí)達(dá)到峰值，采用“列均分-行連續(xù)”策略的第1 個(gè)優(yōu)點(diǎn)，即在保證一列數(shù)據(jù)能夠完整讀入從核陣列的同時(shí)，符合行優(yōu)先的存儲結(jié)構(gòu)，且充分利用了傳輸帶寬，從而有效解決了列讀寫的不連續(xù)性.

3.3 從核間數(shù)據(jù)分布與寄存器通信轉(zhuǎn)置

“列均分-行連續(xù)”策略的第2 個(gè)優(yōu)點(diǎn)，即從核間的數(shù)據(jù)分布與大點(diǎn)數(shù)分解方案吻合，有利于“列FFT”的后續(xù)計(jì)算.如圖4 所示為“大矩陣”列數(shù)據(jù)在從核間的分布，方格中的數(shù)字代表該數(shù)據(jù)所在的從核號.將“大矩陣”任意一列512 個(gè)數(shù)（見圖4（a））排列成64×8 的“中矩陣”（見圖4（b）），再將“中矩陣”的任意一列繼續(xù)轉(zhuǎn)換成8×8 的“小矩陣”（圖4（c）以圖4（b）中的第0 列為例）.此時(shí)，“小矩陣”中的每一列數(shù)據(jù)均保存在同一個(gè)從核中，因此可以直接計(jì)算“小矩陣”“列FFT”.

圖4 從核數(shù)據(jù)分布Fig.4 Data distribution among slave cores

圖4（c）中，“小矩陣”的每一行數(shù)據(jù)分散在同列8 個(gè)從核中，可以利用“申威26010”的寄存器通信機(jī)制來實(shí)現(xiàn)從核間的數(shù)據(jù)交換.使用putc 和getc 指令開啟寄存器列通信模式，同列8 個(gè)核之間彼此將各行數(shù)據(jù)傳到對應(yīng)的從核，并接收其他從核傳給自己的數(shù)據(jù)，即可完成“小矩陣”的列與行轉(zhuǎn)置（見圖4（d）），繼續(xù)執(zhí)行“小矩陣”“行FFT”計(jì)算.在“小矩陣”FFT 計(jì)算完畢后，核間數(shù)據(jù)按圖4（d）中的列優(yōu)先順序排列，實(shí)際上與圖4（b）保持一致.“中矩陣”的轉(zhuǎn)置過程與“小矩陣”類似，但須使用putr 和getr 指令開啟寄存器行通信模式（見圖4（e））.在完成“中矩陣”FFT 后，將數(shù)據(jù)按列優(yōu)先順序?qū)懟刂鞔妫ㄒ妶D4（f））.LDM 內(nèi)的64 列數(shù)據(jù)均按圖4（f）的順序排列，因此寫回過程符合“列均分-行連續(xù)”原則，DMA 傳輸可達(dá)最大速度.

3.4 SIMD FFT

FFT 計(jì)算采用單指令多數(shù)據(jù)流操作（single instruction multiple data，SIMD）進(jìn)行加速.“申威26010”支持寬度為4 個(gè)double 數(shù)據(jù)的向量運(yùn)算（doublev4），但須解決以下2 個(gè)問題.1）若將同組FFT 數(shù)據(jù)放在一個(gè)向量內(nèi)，則不同間距的蝶形會引起數(shù)據(jù)依賴問題（見圖5（a））.2）若分散在不同向量內(nèi)，則離散數(shù)據(jù)裝載/卸載向量的效率較低.

圖5 SIMD FFT 數(shù)據(jù)排列Fig.5 SIMD FFT data arrangement

當(dāng)前每個(gè)從核LDM 內(nèi)的數(shù)據(jù)排列為8×64 的矩陣，行方向上為從“大矩陣”中讀入的一行64 個(gè)數(shù)據(jù)（見圖3），列方向上是“小矩陣”中一列（行）（見圖4（c）、（d））或“中矩陣”中一行（見圖4（e））的8 個(gè)數(shù)據(jù).采用doublev4 數(shù)組進(jìn)行裝載，每個(gè)doublev4 內(nèi)的4 個(gè)子元素互不相關(guān)，不同doublev4內(nèi)相同位置上的子元素來自同一組FFT 數(shù)據(jù)，如圖5（b）所示，既消除了向量內(nèi)數(shù)據(jù)依賴，可以一次完成4 組FFT 計(jì)算，又使得數(shù)據(jù)在LDM 內(nèi)的存儲連續(xù)且邊界對齊，向量的裝載和卸載效率最高.這是“列均分-行連續(xù)”策略的第3 個(gè)優(yōu)點(diǎn).

3.5 旋轉(zhuǎn)因子計(jì)算優(yōu)化

1）和差化積.所有旋轉(zhuǎn)因子系數(shù)均為2 π/N的nk倍（N為矩陣規(guī)模，n和k為坐標(biāo)）.對基本單元cos(2π/N) 和 sin(2π/N) 使用和差化積，用復(fù)數(shù)乘法代替三角函數(shù)，大幅減少三角函數(shù)的調(diào)用次數(shù).

2）數(shù)據(jù)復(fù)用.每個(gè)從核讀入的“大矩陣”多列數(shù)據(jù)，實(shí)際上每一列在圖4 所示的“小矩陣”和“中矩陣”中所處的位置都是一樣的，即坐標(biāo)相同，因此只須計(jì)算第1 列旋轉(zhuǎn)因子并復(fù)用.

3）提前計(jì)算.所有的旋轉(zhuǎn)因子均由矩陣尺寸，即輸入數(shù)據(jù)量決定，可以提前完成計(jì)算.和差化積縮短了計(jì)算時(shí)間，數(shù)據(jù)復(fù)用減少了計(jì)算次數(shù)，將旋轉(zhuǎn)因子計(jì)算與第一次DMA 傳輸同步進(jìn)行，不會給存儲空間和計(jì)算時(shí)間帶來額外負(fù)擔(dān).

此外，所有旋轉(zhuǎn)因子的計(jì)算和相乘采用SIMD操作，可以與FFT 結(jié)合.

3.6 雙緩沖策略和跨步讀寫

“申威26010”為異步DMA 傳輸，當(dāng)數(shù)據(jù)需要在主存與局存之間進(jìn)行多次傳輸時(shí)，可以采用雙緩沖策略.如圖6 所示，數(shù)字表示循環(huán)輪次，奇數(shù)輪和偶數(shù)輪分別占用2 份不同的輸入/輸出緩沖區(qū).可以看出，除了第1 輪讀入和最后1 輪寫回外，在進(jìn)行中間輪次計(jì)算的同時(shí)，可以進(jìn)行下一輪讀入及上一輪寫回，以計(jì)算時(shí)間隱藏通信開銷.

圖6 雙緩沖機(jī)制Fig.6 Double buffering scheme

實(shí)際上，DMA 還包括啟動athread_get（讀入）和athread_put（寫回）的函數(shù)調(diào)用開銷.該步驟為串行過程，無法被計(jì)算隱藏，當(dāng)一次完整讀寫需要多次調(diào)用這2 個(gè)函數(shù)時(shí)，耗時(shí)不可忽略.從核在每一輪“列均分”中共分到8 行數(shù)據(jù)，須調(diào)用get和put 函數(shù)各16 次（實(shí)部和虛部各1 次），因此無法被雙緩沖隱藏的DMA 時(shí)長T=16T1（設(shè)T1為函數(shù)調(diào)用時(shí)間），導(dǎo)致加速提升不到10%.

此時(shí)，每個(gè)從核在主存上的訪問空間呈現(xiàn)等數(shù)據(jù)量（行連續(xù)）、等間隔（列均分）的特點(diǎn)（見圖4（a）、（f））.通過配置athread_get(dma_mode mode,void *src,void *dest,int len,void *reply,char mask,int stride,int bsize)的最后2 個(gè)參數(shù)“間距”和“單次傳輸數(shù)據(jù)量”，可以實(shí)現(xiàn)跨步傳輸.跨步讀寫僅須調(diào)用get 和put 各2 次，無法隱藏的調(diào)用時(shí)間2T3＜＜16T1（設(shè)T3為跨步函數(shù)調(diào)用時(shí)間），達(dá)到了非常明顯的加速效果.

3.7 FFT 完整流程

如圖7 所示為“大矩陣”“列FFT”的完整計(jì)算流程，包括各環(huán)節(jié)采用的優(yōu)化方法.“大矩陣”“行FFT”的計(jì)算過程與“列FFT”類似，此處不再贅述.

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

4.1 實(shí)驗(yàn)運(yùn)行環(huán)境

測試了數(shù)據(jù)量為32 768～4 194 304 共8 組隨機(jī)生成的雙精度浮點(diǎn)復(fù)數(shù).眾核并行程序啟用“申威26010”1 個(gè)核組共64 個(gè)從核，運(yùn)行時(shí)間從athread_spawn 開始至athread_join 結(jié)束，包括傳輸和計(jì)算全過程.對比基準(zhǔn)為單主核上運(yùn)行FFTW 3.3.4庫以最優(yōu)plan 計(jì)算相同輸入的用時(shí)，僅統(tǒng)計(jì)fftw_execute 計(jì)算時(shí)間.所有實(shí)驗(yàn)均重復(fù)多次，取穩(wěn)定值求平均作為測試結(jié)果.

4.2 各優(yōu)化策略的加速效果

如圖8 所示為采用旋轉(zhuǎn)因子優(yōu)化、SIMD、雙緩沖以及改用跨步傳輸后再結(jié)合雙緩沖等各方案的加速比SP.

圖8 采用旋轉(zhuǎn)因子優(yōu)化、SIMD、跨步傳輸、雙緩沖方案的加速比Fig.8 Speedup using twiddle factor optimization,SIMD,stride transmission,double-buffering scheme

4.2.1 計(jì)算環(huán)節(jié) 如表1 所示為優(yōu)化前、后單個(gè)從核內(nèi)旋轉(zhuǎn)因子的計(jì)算量.其中，小點(diǎn)數(shù)FFT 所需的旋轉(zhuǎn)因子數(shù)量分別為4、8、16，可以一次完成計(jì)算，所以忽略不計(jì)，只統(tǒng)計(jì)一維轉(zhuǎn)二維方案中“列FFT”轉(zhuǎn)“行FFT”之間所須乘上的旋轉(zhuǎn)因子.3.5 節(jié)優(yōu)化策略的核心思想是以和差化積所需的復(fù)數(shù)乘法和加法（或4 次實(shí)數(shù)乘法和2 次實(shí)數(shù)加法）來代替耗時(shí)較長的三角函數(shù)調(diào)用.以262 144=512×512 為例，優(yōu)化前需要40 960 次三角函數(shù)，優(yōu)化后僅為24 次三角函數(shù)加上1 136 次浮點(diǎn)運(yùn)算.且該加速效果隨著計(jì)算規(guī)模的增大而提高.

表1 旋轉(zhuǎn)因子優(yōu)化前、后的計(jì)算次數(shù)Tab.1 Calculation times of twiddle factor before and after optimization

實(shí)驗(yàn)結(jié)果表明，若未對旋轉(zhuǎn)因子計(jì)算進(jìn)行優(yōu)化，則多核并行總耗時(shí)比單主核FFTW 長，而通過3.5 節(jié)所述的和差化積、數(shù)據(jù)復(fù)用、提前計(jì)算等優(yōu)化策略，可以達(dá)到15 倍左右的加速比.因此，旋轉(zhuǎn)因子優(yōu)化可以認(rèn)為是計(jì)算環(huán)節(jié)中必須且最主要的加速部分.

在加入SIMD 向量化運(yùn)算后，可以在旋轉(zhuǎn)因子優(yōu)化的基礎(chǔ)上提升近2 倍的加速效果.SIMD 覆蓋了計(jì)算環(huán)節(jié)的幾乎所有部分，包括FFT 以及旋轉(zhuǎn)因子相乘.若單獨(dú)測量運(yùn)算部分，則可以達(dá)到約3.5 倍的加速，接近4 倍峰值，證明利用3.4 節(jié)的方案，降低了裝載/卸載向量引起的額外開銷.

4.2.2 通信環(huán)節(jié) 從核間寄存器單點(diǎn)通信的標(biāo)稱延遲為10 拍，但是在應(yīng)用中通常會疊加put/get 指令執(zhí)行開銷、顯式sync 同步開銷、網(wǎng)絡(luò)堵塞延遲、連續(xù)讀寫延遲等，效率明顯降低，需要采用匯編級流水線優(yōu)化進(jìn)行手動提升.在本實(shí)驗(yàn)中，每次寄存器通信共須傳輸512×7/8=448 個(gè)數(shù)，數(shù)量固定，因此耗時(shí)基本穩(wěn)定，優(yōu)于采用DMA 傳回主存，且改用DMA 傳輸后會導(dǎo)致計(jì)算時(shí)間無法覆蓋傳輸時(shí)間，使雙緩沖失效.寄存器通信是“小/中矩陣”“列/行”轉(zhuǎn)置的首選方案.

在“列均分-行連續(xù)”策略中，直接采用雙緩沖，加速比最大提升不到10%.與get 和put 調(diào)用次數(shù)有關(guān)，數(shù)據(jù)量越大，調(diào)用次數(shù)越多，雙緩沖效果越差.在改用跨步傳輸再加上雙緩沖后，達(dá)到了非常明顯的加速效果，大數(shù)據(jù)量時(shí)可以再提升1 倍以上，與3.6 節(jié)的傳輸過程分析基本吻合.

如圖9 所示為不同輸入點(diǎn)數(shù)時(shí)的DMA 平均傳輸帶寬BW及利用率UBW統(tǒng)計(jì).可以看出，“行連續(xù)”策略能夠保證連續(xù)傳輸字節(jié)數(shù)≥256 B，因此可以在大多數(shù)情況下保持對讀寫帶寬的充分利用.當(dāng)僅在最大數(shù)據(jù)量分解方案為2 048 列時(shí)，連續(xù)讀寫粒度減小為128 B，導(dǎo)致性能下降.以實(shí)測讀帶寬峰值27.9 GB/s 進(jìn)行折算，帶寬利用率最高可達(dá)89.6%，最低為61.6%，平均為79.8%.

圖9 DMA 傳輸帶寬及利用率Fig.9 DMA transmission bandwidth and utilization

如圖10 所示為“主-從核DMA 通信”、“從核間寄存器通信”、“從核FFT 計(jì)算”3 個(gè)環(huán)節(jié)各自在總時(shí)長中的占比TP，左邊為普通傳輸，右邊為跨步讀寫（圖例僅表示各部分功能的相對占比關(guān)系，不表示運(yùn)行時(shí)長的絕對數(shù)值）.可以看出，跨步傳輸+雙緩沖策略使得主/從核間DMA 傳輸時(shí)長占比明顯下降，不再是整體性能的瓶頸，各環(huán)節(jié)耗時(shí)更均衡.

圖10 傳輸和計(jì)算各部分功能的耗時(shí)占比Fig.10 Proportion of time consumed by transmission and calculation

4.2.3 實(shí)驗(yàn)小結(jié) 表2 給出FFTW 串行程序和最終多核并行方案各自的運(yùn)行拍數(shù)CS、CP及兩者的加速比測試結(jié)果，可以得到以下結(jié)論.

表2 并行FFT 的加速測試結(jié)果Tab.2 Accelerated test results of parallel FFT

1）該并行優(yōu)化方案的核心是“列均分-行連續(xù)”策略，優(yōu)點(diǎn)如下.a）列數(shù)據(jù)讀寫連續(xù)非離散，能夠充分利用DMA 帶寬.b）數(shù)據(jù)在從核間的分布能夠直接進(jìn)行后續(xù)FFT.c）數(shù)據(jù)在從核內(nèi)的排列可以實(shí)現(xiàn)高效的SIMD.d）有助于旋轉(zhuǎn)因子的快速運(yùn)算.以上4 點(diǎn)保證了除必需的主存和局存之間的數(shù)據(jù)讀寫以及算法中的矩陣“列/行轉(zhuǎn)置”操作以外，計(jì)算過程中基本沒有額外的數(shù)據(jù)搬移，最大程度上減少了數(shù)據(jù)移動的次數(shù)，有效消除了FFT 訪存密集型特點(diǎn)帶來的影響.

2）“列均分-行連續(xù)”策略的不足，即需要多次調(diào)用讀寫函數(shù)而帶來的開銷，可以采用跨步傳輸+雙緩沖進(jìn)行彌補(bǔ).合理的跨步傳輸能夠有效地減少傳輸次數(shù)，從而提高異步傳輸時(shí)間在DMA總耗時(shí)中的占比，再利用雙緩沖加以隱藏，有利于“列均分”策略的實(shí)施.

3）隨著輸入數(shù)據(jù)量的成倍增長，多核并行總時(shí)長的增加倍數(shù)與計(jì)算數(shù)據(jù)量的增長倍數(shù)基本相當(dāng)，DMA 傳輸粒度保持在256 B 及以上.測試結(jié)果顯示，浮點(diǎn)運(yùn)算性能及主、從核間的傳輸帶寬基本可以達(dá)到并維持在峰值.

4）單主核FFTW 在數(shù)據(jù)量較小時(shí)的運(yùn)算性能較好，在數(shù)據(jù)量大于262 144 之后基本保持成倍增長，應(yīng)該是FFTW 在不同數(shù)據(jù)量下會選擇不同的最優(yōu)算法所致.

5）本實(shí)驗(yàn)最終的加速比可以達(dá)到平均48x、最高65x 以上的出色效果，尤其是在大點(diǎn)數(shù)情況下，能夠始終保持在50x 以上，平均接近60x.

6）利用“列均分-行連續(xù)+跨步傳輸”的策略，能夠有效解決主、從核間“列”數(shù)據(jù)的讀寫問題，在矩陣應(yīng)用場景下可以消除顯式轉(zhuǎn)置帶來的影響，起到良好的加速效果.在新一代“申威26010-Pro”處理器架構(gòu)中，LDM 擴(kuò)大至256 kB，單次可傳輸數(shù)據(jù)量增加，更有利于DMA 帶寬的充分利用.從核間的寄存器傳輸改為RMA 通信機(jī)制，性能得到進(jìn)一步的提升.本方法重點(diǎn)關(guān)注處理器架構(gòu)與算法的內(nèi)在適配性，無須對軟件算法框架進(jìn)行大幅調(diào)整，因此可以作為一種通用優(yōu)化策略進(jìn)行推廣使用，為“申威26010”高性能數(shù)學(xué)庫中針對FFT算法的并行優(yōu)化做出一定的貢獻(xiàn).

5 結(jié)語

本文介紹了在國產(chǎn)“申威26010”眾核處理器上對大點(diǎn)數(shù)FFT 進(jìn)行并行加速的方案.為了有效地解決二維FFT 對列數(shù)據(jù)的讀寫和計(jì)算問題，消除矩陣轉(zhuǎn)置帶來的影響，特別提出“列均分-行連續(xù)”的讀寫策略，結(jié)合SIMD 向量操作、寄存器通信轉(zhuǎn)置、旋轉(zhuǎn)因子優(yōu)化、雙緩沖+跨步傳輸策略等優(yōu)化方法，實(shí)現(xiàn)了經(jīng)典Cooley-Tukey FFT 算法在“申威26010”上的眾核優(yōu)化.與單主核FFTW 相比，可以達(dá)到平均48x 以上、峰值65x 以上的加速比.