郭金鑫，張廣婷，張?jiān)迫悵扇A，賈海鵬

1.太原理工大學(xué) 大數(shù)據(jù)學(xué)院，太原 030024

2.中國科學(xué)院 計(jì)算技術(shù)研究所 計(jì)算機(jī)體系結(jié)構(gòu)國家重點(diǎn)實(shí)驗(yàn)室，北京 100190

快速傅里葉變換（fast Fourier transform，F(xiàn)FT）是處理器基礎(chǔ)軟件生態(tài)最關(guān)鍵的算法之一，是計(jì)算離散傅里葉變換（discrete Fourier transform，DFT）或其逆運(yùn)算的快速算法，并將算法復(fù)雜度由()降為了(lb)。FFT 算法被用于物理、天文學(xué)、工程、應(yīng)用數(shù)學(xué)、密碼學(xué)和計(jì)算金融等許多不同的領(lǐng)域。如在國際大科學(xué)工程——平方公里陣列射電望遠(yuǎn)鏡（square kilometer array，SKA）項(xiàng)目中，F(xiàn)FT 是數(shù)據(jù)處理的五大算法之一，其計(jì)算量占總計(jì)算量的40%。由于各應(yīng)用領(lǐng)域的急速發(fā)展及其實(shí)時(shí)性能的要求持續(xù)提高，F(xiàn)FT 在ARM（特別是ARMv8）和X86-64 架構(gòu)平臺(tái)上高性能的實(shí)現(xiàn)和優(yōu)化有著重要的研究意義和應(yīng)用價(jià)值。

雖然FFT算法在ARM和X86-64平臺(tái)上已經(jīng)有比較成熟的實(shí)現(xiàn)，如ARMPL（ARM performance library）、Intel MKL（math kernel library）和FFTW（fastest Fourier transform in the West）。但是由于FFT 算 法的復(fù)雜性和多樣性，依然有許多工作值得深入研究。例如在庫利-圖基FFT 算法這一目前應(yīng)用最為廣泛和流行的快速傅里葉算法中，依然存在蝶形網(wǎng)絡(luò)復(fù)雜、蝶形計(jì)算復(fù)雜多樣等問題。特別是對(duì)于大基的實(shí)現(xiàn)，雖然大基通過減少訪存提升性能，但是大基蝶形實(shí)現(xiàn)依然存在匯編實(shí)現(xiàn)復(fù)雜、寄存器不夠用等問題。本文針對(duì)這些問題，研究FFT 算法在不同架構(gòu)CPU 上的高性能實(shí)現(xiàn)方法，突破以上問題導(dǎo)致的性能瓶頸，從而實(shí)現(xiàn)了一個(gè)高性能FFT 算法庫。

在本文的研究中，F(xiàn)FT 算法的實(shí)現(xiàn)和優(yōu)化主要從如下三方面進(jìn)行：（1）蝶形網(wǎng)絡(luò)重構(gòu)，優(yōu)化不同基特別是一些大的基，降低蝶形網(wǎng)絡(luò)級(jí)數(shù)，減少訪存提升蝶形網(wǎng)絡(luò)性能；（2）利用DFT 矩陣性質(zhì)，提取蝶形計(jì)算公共項(xiàng)，將大基蝶形計(jì)算化到最簡(jiǎn)；（3）蝶形計(jì)算匯編實(shí)現(xiàn)，匯編SIMD（single instruction multiple data）優(yōu)化，寄存器復(fù)用策略制定和堆棧內(nèi)存使用解決寄存器不夠用等。通過以上優(yōu)化方法的使用，本文在ARMv8和X86-64 計(jì)算平臺(tái)上突破了大基寄存器不夠用的性能瓶頸，實(shí)現(xiàn)和優(yōu)化了一個(gè)高性能快速傅里葉變換算法庫。實(shí)驗(yàn)結(jié)果表明本文實(shí)現(xiàn)的FFT 算法庫相較FFTW、ARMPL 以及Intel MKL 性能有較大提升，相較算法中小基性能也有較大提升。

本文的主要貢獻(xiàn)如下：

（1）總結(jié)和重構(gòu)蝶形網(wǎng)絡(luò)，同時(shí)利用DFT 矩陣的對(duì)稱性和周期性，大幅降低了大基蝶形計(jì)算的復(fù)雜度；

（2）總結(jié)設(shè)計(jì)了基14、基20 等大基FFT 蝶形計(jì)算方法特別是寄存器使用策略，解決了由于寄存器不夠用導(dǎo)致的性能瓶頸；

（3）提出了一套FFT 算法在ARMv8 及X86-64 架構(gòu)上的實(shí)現(xiàn)策略和優(yōu)化方案，并構(gòu)建了一個(gè)可跨平臺(tái)移植的高性能FFT 算法庫。

1 相關(guān)背景

1.1 Cooley-Tukey FFT 算法

離散傅里葉變換是一種用于進(jìn)行傅里葉分析的基本離散變換，定義如下：

Cooley-Tukey 算法是在許多實(shí)際應(yīng)用中應(yīng)用最廣泛的快速傅里葉變換（FFT）算法。采取分而治之的方法，通過遞歸將大的DFT 分解為小的DFT。

為簡(jiǎn)化DFT 運(yùn)算，利用DFT 矩陣的對(duì)稱性周期性將計(jì)算時(shí)間復(fù)雜度由()降到(lb)。

1.2 ARMv8 架構(gòu)

ARM 是一種負(fù)載存儲(chǔ)體系結(jié)構(gòu)，是RISC 處理器的典型。ARMv8 是ARMv7 之后的下一個(gè)旗艦架構(gòu)，向后兼容ARMv7，是首次支持64 位指令集的ARM處理器架構(gòu)，引入64 位體系結(jié)構(gòu)的同時(shí)保持了與現(xiàn)有32 位系統(tǒng)結(jié)構(gòu)的兼容性。

ARMv8 提供了31×64 bit 通用寄存器及32 個(gè)128 bit 浮點(diǎn)寄存器V0～V31（如圖1），浮點(diǎn)寄存器在執(zhí)行指令時(shí)一條指令可以操作多個(gè)操作數(shù)，可以提高指令的執(zhí)行效率，提高性能。在SIMD 優(yōu)化中，浮點(diǎn)寄存器起著重要作用。

圖1 ARMv8 架構(gòu)浮點(diǎn)寄存器圖Fig.1 ARMv8 architecture float register

1.3 X86-64 架構(gòu)

X84-64 又 稱Intel 64 或AMD64,是X86 指令 集64位版本。Haswell 架構(gòu)是因特爾公司首個(gè)支持AVX2的X86 架構(gòu)。X86-64 處理器架構(gòu)提供了17×64 bit通用寄存器（RDI、RSI、RDX、RCX、R8-R15、RAX、RBX、RBP、RSP、RIP），其中RDI、RSI、RDX、RCX、R8、R9作為函數(shù)輸入?yún)?shù)；提供了16個(gè)256 bit浮點(diǎn)寄存器YMM0～YMM15。Haswell 向后兼容，寄存器低128 bit 可作為128 bit 浮點(diǎn)寄存器XMM0～XMM15。其中XMM0～XMM7 可作為函數(shù)輸入?yún)?shù)。XMM 寄存器每條指令可同時(shí)處理4 個(gè)float 浮點(diǎn)數(shù)，YMM 寄存器每條指令可處理8 個(gè)float 浮點(diǎn)數(shù)。在SIMD 優(yōu)化，浮點(diǎn)寄存器起著重要作用。

2 相關(guān)工作

FFT 算法的研究以主流FFT 庫為主，在特定的硬件架構(gòu)上實(shí)現(xiàn)高性能，包括FFTW、ARMPL、Intel MKL、AOCL（AMD optimizing CPU libraries）、CUFFT（CUDA fast Fourier transform library）等。

2.1 FFTW

FFTW 是MIT 的Frigo 和Johnson 開發(fā)的自適應(yīng)優(yōu)化FFT 軟件包，用于計(jì)算一維或多個(gè)維度、任意輸入大小、實(shí)數(shù)和復(fù)數(shù)數(shù)據(jù)以及偶數(shù)和奇數(shù)數(shù)據(jù)的離散傅里葉變換DFT。同時(shí)支持共享存儲(chǔ)多線程并行和MP（Imessage passing interface）并行，其運(yùn)算性能遠(yuǎn)遠(yuǎn)領(lǐng)先目前已有的其他FFT 軟件。FFTW 性能可移植，在大多數(shù)架構(gòu)上性能良好，并且自FFTW3.3.1 開始針對(duì)ARM 平臺(tái)實(shí)現(xiàn)了較高的性能。

2.2 ARMPL

ARMPL 性能庫是ARM 公司針對(duì)ARMv8 平臺(tái)推出的高性能商業(yè)庫。為ARM 處理器上的高性能計(jì)算應(yīng)用程序提供了標(biāo)準(zhǔn)核心數(shù)學(xué)庫，包含優(yōu)化的BLAS（basic linear algebra subprograms）、LAPACK（linear algebra package）和FFT，為FFT 計(jì)算提供了與FFTW3 相同的接口。

2.3 MKL

MKL 是一個(gè)用于科學(xué)、工程和金融應(yīng)用程序的包含快速傅里葉變換的優(yōu)化數(shù)學(xué)例程庫。Intel MKL FFTW 是因特爾公司在FFTW 基礎(chǔ)上二次開發(fā)的商業(yè)FFT 性能庫，是目前X86 平臺(tái)上性能最好的FFT 商業(yè)庫，但其只用于X86 架構(gòu)，可移植性差。

3 FFT 算法的實(shí)現(xiàn)和優(yōu)化

3.1 蝶形網(wǎng)絡(luò)優(yōu)化

在Cooley-Tukey FFT 算法中，蝶形網(wǎng)絡(luò)決定了數(shù)據(jù)訪問模式和蝶形計(jì)算執(zhí)行順序。DFT 是逐級(jí)求解的，每級(jí)重復(fù)處理蝶形計(jì)算，因此蝶形網(wǎng)絡(luò)的組織方式影響整個(gè)優(yōu)化。相同的蝶形網(wǎng)絡(luò)，不同的實(shí)現(xiàn)和優(yōu)化可能導(dǎo)致不同的性能。依照蝶形因子在計(jì)算中出現(xiàn)的不同位置，實(shí)現(xiàn)該算法有兩種方式：時(shí)域抽取（decimation-in-time，DIT）和頻域抽取（decimation-infrequency，DIF）。時(shí)域抽取時(shí)，蝶形因子在計(jì)算輸入端，輸入向量需位反轉(zhuǎn)，輸出自然順序；頻域抽取則相反，蝶形因子在計(jì)算輸出端，輸入向量自然順序，輸出需位反轉(zhuǎn)。

傳統(tǒng)蝶形網(wǎng)絡(luò)存在位反轉(zhuǎn)操作，如圖2 所示，增加了額外的內(nèi)存成本，還增加了混合基建立的困難度，對(duì)SIMD 不友好。本文采用了如圖3所示的Stockham蝶形網(wǎng)絡(luò)結(jié)構(gòu)。

圖2 時(shí)域抽取Fig.2 DIT network

圖3 Stockham 蝶形網(wǎng)絡(luò)Fig.3 Stockham butterfly network

Stockham 蝶形網(wǎng)絡(luò)結(jié)構(gòu)相比傳統(tǒng)蝶形網(wǎng)絡(luò)：（1）去除了位反轉(zhuǎn)排列，DIT 時(shí)域抽取需要將輸入序列重新排序?yàn)槲环崔D(zhuǎn)順序，位反轉(zhuǎn)排列引入了額外的不連續(xù)的內(nèi)存訪問，不連續(xù)的內(nèi)存訪問會(huì)導(dǎo)致統(tǒng)一輸入輸出的內(nèi)存訪問困難。Stockham 蝶形網(wǎng)絡(luò)各級(jí)計(jì)算的輸入輸出都是自然順序，消除了位反轉(zhuǎn)排列，統(tǒng)一了蝶形網(wǎng)絡(luò)的訪存行為。（2）SIMD 友好，SIMD 為一條指令作用在多個(gè)數(shù)據(jù)操作上，為了有效使用SIMD，從內(nèi)存中加載和存儲(chǔ)到內(nèi)存中的數(shù)據(jù)應(yīng)是連續(xù)的，Stockham 蝶形網(wǎng)絡(luò)結(jié)構(gòu)中，蝶形網(wǎng)絡(luò)的輸入輸出是連續(xù)定位的，在同一級(jí)內(nèi)SIMD 并行化友好。（3）混合基友好，由于每級(jí)的輸入輸出都是自然順序，不同RADIX 算法采取統(tǒng)一的方式，可以完美融合在一起。

為了得到更好的性能，將第一級(jí)蝶形網(wǎng)絡(luò)單獨(dú)優(yōu)化。第一級(jí)蝶形網(wǎng)絡(luò)的旋轉(zhuǎn)因子為1，沒有必要從內(nèi)存中讀取和計(jì)算，降低了不必要的內(nèi)存訪問和計(jì)算成本。第一級(jí)輸出結(jié)果的寫入并不連續(xù)，在匯編優(yōu)化時(shí)需要再進(jìn)行數(shù)據(jù)重組和轉(zhuǎn)置。

蝶形網(wǎng)絡(luò)帶寬依賴較高，由于每一級(jí)內(nèi)存訪存寫入，蝶形網(wǎng)絡(luò)級(jí)數(shù)過多會(huì)增加數(shù)據(jù)訪問量。使用大基參與FFT 計(jì)算可降低蝶形網(wǎng)絡(luò)的級(jí)數(shù)，減少數(shù)據(jù)訪問量。雖然大基參與蝶形網(wǎng)絡(luò)計(jì)算，寄存器不夠用，一定程度上降低了蝶形網(wǎng)絡(luò)性能，但使用大基參與計(jì)算帶來的性能增益優(yōu)于性能損耗。

3.2 蝶形計(jì)算優(yōu)化

在FFT 計(jì)算過程中，蝶形計(jì)算反復(fù)調(diào)用，蝶形計(jì)算的性能將直接影響FFT 算法的最終性能。因此，本節(jié)將介紹如何將FFT 蝶形計(jì)算的復(fù)雜度降到最低。根據(jù)離散傅里葉變換式（1），基(Radix-)的蝶形，本質(zhì)上就是長(zhǎng)度為的DFT 計(jì)算，而DFT 計(jì)算的實(shí)質(zhì)即DFT 矩陣與輸入矩陣向量乘法。

由于當(dāng)基較小時(shí)已經(jīng)有了成熟的計(jì)算方案，本文將研究大基的高性能實(shí)現(xiàn)方法，如基14 和基20。

下面將詳細(xì)地分析Radix-14 的蝶形計(jì)算方法。Radix-14 蝶形計(jì)算本質(zhì)上是數(shù)據(jù)規(guī)模為14 的DFT計(jì)算。

圖4 Radix-14 旋轉(zhuǎn)因子復(fù)平面分布圖Fig.4 Radix-14 twiddles complex plane distribution

根據(jù)Radix-14 在如圖4 所示復(fù)平面上的分布旋轉(zhuǎn)因子關(guān)于軸和軸對(duì)稱：實(shí)部相同虛部互為相反數(shù)；虛部相同實(shí)部互為相反數(shù)。

式（1）旋轉(zhuǎn)因子具有如下性質(zhì)：

通過提取和預(yù)計(jì)算公因子，可以減少浮點(diǎn)計(jì)算和代碼冗余，將蝶形計(jì)算的時(shí)間復(fù)雜度降到最低。

3.3 SIMD 優(yōu)化

ARMv8 提供了32 個(gè)128 位的浮點(diǎn)寄存器，每個(gè)浮點(diǎn)寄存器可以存儲(chǔ)4 個(gè)32 位的單精度float 浮點(diǎn)數(shù)或2 個(gè)64 位雙精度double浮點(diǎn)數(shù)，一條指令最多可以同時(shí)并行處理4 個(gè)數(shù)據(jù)。在ARMv8 體系結(jié)構(gòu)中，使用2個(gè)128位寄存器分別容納4個(gè)復(fù)數(shù)的實(shí)部和虛部。

X86-64 架構(gòu)提供了16 個(gè)256 位的浮點(diǎn)寄存器，X86 架構(gòu)SSE 指令可以用128 位通路XMM 浮點(diǎn)寄存器處理4 個(gè)32 位的運(yùn)算或處理2 個(gè)64 位的運(yùn)算。X86-64 AVX 指令是SSE 的兩倍，可操作16 個(gè)YMM 256 位浮點(diǎn)寄存器。并行操作1～8 個(gè)單精度float 浮點(diǎn)數(shù)，1～4 個(gè)雙精度double 浮點(diǎn)數(shù)。在Haswell X86-64 體系結(jié)構(gòu)中，使用一個(gè)256 位寄存器交錯(cuò)容納4 個(gè)復(fù)數(shù)的實(shí)部和虛部。

FFT 計(jì)算時(shí)，每個(gè)蝶形計(jì)算相互獨(dú)立。如圖5 所示SIMD 優(yōu)化同時(shí)處理4 個(gè)蝶形計(jì)算，提高了程序的并行效率。

圖5 SIMD 優(yōu)化Fig.5 SIMD optimization

高級(jí)語言程序，一般由編譯器負(fù)責(zé)寄存器使用。為提高FFT 算法性能，本文蝶形計(jì)算過程采用匯編語言。通過寄存器使用優(yōu)化，可提高算法的性能。寄存器使用的主要思想是寄存器分組。浮點(diǎn)寄存器的使用分為輸入寄存器in，旋轉(zhuǎn)因子寄存器twiddles，中間結(jié)果寄存器scratch 及輸出結(jié)果寄存器out。ARM 架構(gòu)復(fù)數(shù)實(shí)部虛部需要×2 個(gè)in 寄存器，(-1)×2 個(gè)tw 寄存器，×2 個(gè)輸出out寄存器；X86-64架構(gòu)需要個(gè)in寄存器，-1 個(gè)tw寄存器，個(gè)out寄存器。隨著基數(shù)的增長(zhǎng)，4 組寄存器需要更多的向量寄存器，寄存器資源無法獨(dú)立完成蝶形計(jì)算。

寄存器優(yōu)化分為兩部分：一部分為寄存器復(fù)用；一部分為極大基臨時(shí)存入堆棧或內(nèi)存。寄存器復(fù)用策略：復(fù)用旋轉(zhuǎn)因子臨時(shí)寄存器tmptw，輸入和旋轉(zhuǎn)因子數(shù)據(jù)相互獨(dú)立，寄存器使用緊張時(shí)分批加載旋轉(zhuǎn)因子，完成復(fù)數(shù)乘法后，繼續(xù)復(fù)用tmptw，在全部完成旋轉(zhuǎn)因子輸入乘法后，釋放tmptw；復(fù)用臨時(shí)輸入寄存器tmpin，在加載部分輸入后，計(jì)算中間變量，釋放tmpin 以復(fù)用；復(fù)用臨時(shí)寄存器tmp，在乘法等運(yùn)算時(shí)需要臨時(shí)寄存器，這時(shí)的臨時(shí)寄存器需及時(shí)復(fù)用處理運(yùn)算；復(fù)用tmpout 輸出寄存器，在獲得輸出后，立即存儲(chǔ)在內(nèi)存，釋放tmpout以復(fù)用。

（1）在大基的情況下如Radix-14，浮點(diǎn)寄存器的合理充分利用，直接影響FFT計(jì)算程序性能。ARM架構(gòu)寄存器合理復(fù)用的同時(shí)存在不夠用的情況。需將公共因子臨時(shí)存入堆棧，計(jì)算輸出結(jié)果時(shí)載入寄存器。

ARM 架構(gòu)中，Radix-14 的Kernel 計(jì)算：中間結(jié)果scratch 實(shí)部虛部分別需要26 個(gè)浮點(diǎn)寄存器。輸入數(shù)據(jù)和旋轉(zhuǎn)因子成對(duì)載入，輸入數(shù)據(jù)乘以旋轉(zhuǎn)因子后載入下一對(duì)數(shù)據(jù)時(shí)計(jì)算中間結(jié)果，復(fù)用旋轉(zhuǎn)因子和輸入數(shù)據(jù)寄存器；計(jì)算中間結(jié)果后，+類中間結(jié)果寄存器可復(fù)用計(jì)算；計(jì)算過程中將通過寄存器分配無法分配溢出的中間數(shù)據(jù)+，-存入堆棧，計(jì)算輸出out 時(shí)取出；蝶形計(jì)算結(jié)果及時(shí)輸出釋放復(fù)用輸出寄存器。寄存器復(fù)用中間處理足夠多的指令可消除相鄰指令寄存器依賴，達(dá)到蝶形計(jì)算寄存器最大化合理利用。

X86-64 架構(gòu)中，Radix-14 的Kernel 計(jì)算:中間結(jié)果scratch 需要26 個(gè)YMM，旋轉(zhuǎn)因子寄存器在乘以輸入后重復(fù)復(fù)用；輸入寄存器計(jì)算中間結(jié)果后釋放復(fù)用為中間結(jié)果寄存器。通過寄存器合理復(fù)用即可實(shí)現(xiàn)Radix-14 蝶形計(jì)算。

（2）在極大基的情況下，如Radix-20，X86-64 架構(gòu)提供的16 個(gè)YMM 寄存器，ARM 架構(gòu)提供的32 個(gè)V0～V31 浮點(diǎn)寄存器遠(yuǎn)不夠用。通過復(fù)用寄存器，寄存器仍不夠用，此時(shí)需要使用堆?；騼?nèi)存指令暫時(shí)保存相關(guān)數(shù)據(jù)，蝶形計(jì)算需要時(shí)取出載入寄存器。

ARM 架構(gòu)中，Radix-20 的Kernel 計(jì)算：中間結(jié)果scratch 實(shí)部虛部分別需要50 個(gè)浮點(diǎn)寄存器。數(shù)據(jù)載入時(shí)分奇偶序列分開載入；輸入寄存器和旋轉(zhuǎn)因子乘運(yùn)算后，復(fù)用旋轉(zhuǎn)因子寄存器；輸入數(shù)據(jù)載入后及時(shí)計(jì)算，復(fù)用輸入數(shù)據(jù)占用的寄存器；計(jì)算過程中寄存器占滿32 個(gè)，通過寄存器分配策略將計(jì)算輸出頻繁用到的S 類中間數(shù)據(jù)存入堆棧釋放所需寄存器，需要使用時(shí)再取出；蝶形計(jì)算結(jié)果及時(shí)輸出釋放復(fù)用寄存器。

X86-64 架構(gòu)中，Radix-20 的Kernel 計(jì)算:中間結(jié)果scratch 需要50 個(gè)YMM。數(shù)據(jù)成對(duì)載入，及時(shí)計(jì)算中間變量，復(fù)用輸入寄存器；寄存器分配策略將無法占用寄存器的相關(guān)數(shù)據(jù)存入內(nèi)存，計(jì)算輸出結(jié)果時(shí)載入寄存器；蝶形網(wǎng)絡(luò)計(jì)算第一級(jí)輸出結(jié)果存儲(chǔ)轉(zhuǎn)置時(shí)，還需將部分輸出結(jié)果存入內(nèi)存，數(shù)據(jù)操作時(shí)載入寄存器。

指令選擇時(shí)，選擇延遲低吞吐量高的指令。X86-64 選擇vfnmadd231ps 和vfmadd231ps 乘加指令，ARMv8 選擇fmla 和fmls 乘加乘減指令。AVX2 還提供了vaddsubps 指令來完成交錯(cuò)模式下的復(fù)數(shù)乘法。如圖6 列出了幾種運(yùn)算的指令對(duì)比。在ARM NEON中支持ld2、st2 高效的加載存取指令，因此ARMv8 使用ld2、faddp、st2 等指令進(jìn)行復(fù)數(shù)算數(shù)運(yùn)算。

圖6 計(jì)算指令對(duì)比Fig.6 Instruction comparison

圖7 順序執(zhí)行與指令重排對(duì)比Fig.7 Sequential execution and instruction rearrangement

相鄰指令間沒有依賴關(guān)系時(shí)可指令重排，避免了流水線stall。指令重排和順序指令相比并不影響計(jì)算結(jié)果，但性能會(huì)有一定提升。如圖7 是基14 指令重排對(duì)比圖，順序執(zhí)行過長(zhǎng)時(shí)間占用寄存器，寄存器還不夠使用，寄存器得不到有效利用，指令重排將載入數(shù)據(jù)的順序調(diào)整后及時(shí)計(jì)算中間結(jié)果。在中間結(jié)果計(jì)算后，釋放輸入寄存器的占用，供輸出結(jié)果復(fù)用。在Radix-14、20 這類大基計(jì)算時(shí)，指令重排配合寄存器的合理分配優(yōu)化，一定程度上提高了算法的計(jì)算性能。

4 性能評(píng)估

4.1 測(cè)試環(huán)境

本文采用華為鯤鵬920 CPU和IntelXeonCPU E5-2640 V4 作為性能測(cè)試平臺(tái)。華為鯤鵬920 CPU采用ARMv8 架構(gòu)，IntelXeonCPU E5-2640 V4 采用X86 架構(gòu)。本文實(shí)驗(yàn)條件如表1 所示。

由于FFTW、ARMPL、Intel MKL 是應(yīng)用最廣泛、最成熟的FFT 算法庫，將OpenFFT 的性能與這些庫進(jìn)行了比較。采用FFTW3.3.8 和ARM 公司的商業(yè)庫ARMPL20.0.0 在ARMv8 平臺(tái)上進(jìn)行性能對(duì)比；采用FFTW3.3.8 和Intel MKL 在X86-64 平臺(tái)進(jìn)行對(duì)比。本文實(shí)現(xiàn)的高性能FFT 算法庫為OpenFFT。

表1 實(shí)驗(yàn)環(huán)境Table 1 Experimental environment

4.2 性能分析

本文測(cè)試數(shù)據(jù)維度為一維，數(shù)據(jù)規(guī)模為14×20×（以下性能分析圖橫坐標(biāo)），輸入輸出均為復(fù)數(shù)序列。性能評(píng)估以每秒所執(zhí)行的浮點(diǎn)次數(shù)（giga floating-point operation per second，Gflops）為單位（以下性能分析圖縱坐標(biāo)）。

圖8給出了OpenFFT、ARMPL和FFTW在ARMv8體系結(jié)構(gòu)上的一維C2C 的性能。對(duì)于單精度和雙精度序列，OpenFFT 算法庫的性能整體高于FFTW 和ARMPL 兩個(gè)算法庫。

圖8 ARM 1D C2C FFT 性能對(duì)比Fig.8 Performance comparison of ARM 1D C2C FFT

（1）單精度Float

如圖8（a）所示，OpenFFT 在ARMv8 下一維C2C FFT 變換優(yōu)化結(jié)果中，相對(duì)于ARMPL 實(shí)現(xiàn)了平均31.90%的加速比，最大加速比為51.00%，最小加速比為3.58%；相對(duì)于FFTW 實(shí)現(xiàn)了平均95.50%的加速比，最大加速比為145.80%，最小加速比為31.00%。圖8（b）所示相對(duì)于ARMPL 實(shí)現(xiàn)了平均32.00%的加速比，最大加速比為51.00%，最小加速比為4.00%；相對(duì)于FFTW 實(shí)現(xiàn)了平均98.00%的加速比，最大加速比為155.00%，最小加速比為36.10%。

通過比較OpenFFT、ARMPL、FFTW 性能曲線，三者在性能走勢(shì)上大致相同，從圖8（a）、圖8（b）可知，OpenFFT 性能在輸入規(guī)模為5 600 時(shí)，開始下降，主要原因在于，數(shù)據(jù)規(guī)模較小時(shí)，數(shù)據(jù)可存儲(chǔ)在Cache中，增加了Cache 命中率，減少了訪存開銷，導(dǎo)致小規(guī)模性能整體高于大規(guī)模。

（2）雙精度Double

如圖8（c）所示，OpenFFT 在ARMv8 下一維C2C Double FFT 變換優(yōu)化結(jié)果中，相對(duì)于ARMPL 實(shí)現(xiàn)了平均4.30%的加速比，最大加速比為19.10%，最小加速比為1.81%；相對(duì)于FFTW 實(shí)現(xiàn)了平均27.90%的加速比，最大加速比為42.80%，最小加速比為2.90%。圖8（d）所示相對(duì)于ARMPL 實(shí)現(xiàn)了平均5.36%的加速比，最大加速比為22.20%，最小加速比為0.70%；相對(duì)于FFTW 實(shí)現(xiàn)了平均35.00%的加速比，最大加速比為47.60%，最小加速比為8.00%。

從圖8 可知，雙精度加速性能相對(duì)于單精度要低。主要原因在于SIMD 優(yōu)化，雙精度Double 為64位，ARM 浮點(diǎn)寄存器為128 位，只能循環(huán)展開2 次，寄存器一次處理兩個(gè)數(shù)據(jù)；再有ARM 部分Double 浮點(diǎn)數(shù)指令性能低。圖8（c）、圖8（d）在數(shù)據(jù)規(guī)模為560、840 存在性能低于ARMPL 的情況，原因在于Cache命中率低，訪存開銷大且延遲高，再有數(shù)據(jù)預(yù)取不恰當(dāng)，這些都有可能造成性能不高。

圖9 給出了OpenFFT、MKL 和FFTW 在X86-64體系結(jié)構(gòu)上的一維C2C 的性能。對(duì)于單精度和雙精度序列，OpenFFT 算法庫的性能同樣整體高于FFTW和MKL 兩個(gè)算法庫。

圖9 X86-64 1D C2C FFT 性能對(duì)比Fig.9 Performance comparison of X86-64 1D C2C FFT

（1）單精度Float

如圖9（a）所示，OpenFFT 在X86-64 下一維C2C FFT變換優(yōu)化結(jié)果中，相對(duì)于MKL實(shí)現(xiàn)了平均26.00%的加速比，最大加速比為76.00%，最小加速比為0.92%；相對(duì)于FFTW 實(shí)現(xiàn)了平均70.00%的加速比，最大加速比為155.00%，最小加速比為3.60%。圖9（b）所示相對(duì)于MKL 實(shí)現(xiàn)了平均29.40%的加速比，最大加速比為55.20%，最小加速比為3.60%；相對(duì)于FFTW 實(shí)現(xiàn)了平均81.80%的加速比，最大加速比為175.00%，最小加速比為11.10%。

通過比較OpenFFT、MKL、FFTW 性能曲線，三者在性能走勢(shì)上大致相同，OpenFFT 性能走勢(shì)最高，其次是MKL，性能最低的是FFTW。從圖9（a）、圖9（b）可知，OpenFFT 性能在輸入規(guī)模為3 920 時(shí)達(dá)到最高點(diǎn)，隨后隨著規(guī)模的變大整體性能略下降后趨于穩(wěn)定。小規(guī)模性能高在于數(shù)據(jù)能存在Cache 中，增加了Cache命中率，因此小規(guī)模性能整體高于大規(guī)模。

（2）雙精度Double

如圖9（c）所示，OpenFFT 在X86-64 下一維C2C Double FFT 變換優(yōu)化結(jié)果中，相對(duì)于MKL 實(shí)現(xiàn)了平均45.00%的加速比，最大加速比為126.00%，最小加速比為18.80%；相對(duì)于FFTW 實(shí)現(xiàn)了平均50.50%的加速比，最大加速比為113.00%，最小加速比為6.50%。圖9（d）所示相對(duì)于MKL 實(shí)現(xiàn)了平均33.80%的加速比，最大加速比為96.00%，最小加速比為1.50%；相對(duì)于FFTW 實(shí)現(xiàn)了平均52.70%的加速比，最大加速比為110.00%，最小加速比為2.98%。

表2 總結(jié)了在ARMv8 和X86-64 架構(gòu)在數(shù)據(jù)規(guī)模為14×20×一維C2C 復(fù)數(shù)序列下，OpenFFT 分別與FFTW、ARMPL 和MKL 算法庫性能對(duì)比的平均加速和最大加速。實(shí)驗(yàn)表明OpenFFT 性能明顯優(yōu)于FFTW、ARMPL 和Intel MKL FFT 算法庫。

表2 OpenFFT 平均和最大加速Table 2 Average and maximum speedups of OpenFFT %

圖10、圖11 給出了OpenFFT 在ARMv8 體系結(jié)構(gòu)和X86-64 體系結(jié)構(gòu)上同一數(shù)據(jù)規(guī)模大基Radix-14、Radix-20 和中小基Radix-10、Radix-7 等一維C2C性能對(duì)比。對(duì)于單精度和雙精度序列，同一數(shù)據(jù)規(guī)模OpenFFT 大基性能總體優(yōu)于小基性能，大基對(duì)于總體性能的提升大于大基帶來的性能損耗。

（1）ARMv8 大小基性能對(duì)比

①Float 如圖10（a）所示，OpenFFT 在ARMv8 下單精度一維大基C2C FFT 變換，相對(duì)于小基實(shí)現(xiàn)了平均0.16%的加速比，最大加速比為8.70%；大基跟小基在性能走勢(shì)上大致相同，在數(shù)據(jù)規(guī)模達(dá)到54 880性能優(yōu)于中小基。從圖10（a）可知，大基在輸入規(guī)模為560、840、11 760、16 800，即為Radix-2 或Radix-3時(shí)，性能低于中小基，主要原因是ARM 體系結(jié)構(gòu)單精度下Radix-2、Radix-3 對(duì)性能的影響，數(shù)據(jù)規(guī)模分解方式不唯一，ARM 單精度下異常性能的數(shù)據(jù)規(guī)模Radix-20*Radix-14*Radix-2 或Radix-3 不是最優(yōu)分配造成的性能差異。

圖10 ARMv8 1D C2C FFT 性能對(duì)比Fig.10 Performance comparison of ARMv8 1D C2C FFT

圖11 X86-64 1D C2C FFT 性能對(duì)比Fig.11 Performance comparison of X86-64 1D C2C FFT

②Double 如圖10（b）所示，OpenFFT 在ARMv8下雙精度一維大基C2C FFT 變換，相對(duì)于中小基實(shí)現(xiàn)了平均10.00%的加速比，最大加速比為20.70%；大基跟中小基相比，性能優(yōu)于中小基。

（2）X86-64 大小基性能對(duì)比

①Float 如圖11（a）所示，OpenFFT 在X86-64 下單精度一維大基C2C FFT 變換，相對(duì)于小基實(shí)現(xiàn)了平均17.00%的加速比，最大加速比為35.90%；大基性能優(yōu)于中小基，在數(shù)據(jù)規(guī)模超過16 800，大基性能優(yōu)勢(shì)再次拉大。

②Double如圖11（b）所示，OpenFFT 在X86-64 下雙精度一維大基C2C FFT 變換，相對(duì)于中小基實(shí)現(xiàn)了平均33.30%的加速比，最大加速比為72.80%；大基跟中小基相比，性能優(yōu)于中小基，性能優(yōu)勢(shì)明顯。

實(shí)驗(yàn)結(jié)果表明，大基雖然存在寄存器不夠使用，計(jì)算復(fù)雜的問題，一定程度上降低了蝶形計(jì)算性能，但大基減少了計(jì)算過程中的訪存，一定程度提升了性能，大基性能增益明顯大于性能的損耗。

5 結(jié)束語

本文在原有FFT 基礎(chǔ)上，通過蝶形網(wǎng)絡(luò)優(yōu)化、大基網(wǎng)絡(luò)級(jí)數(shù)降低減少訪存、大基蝶形計(jì)算優(yōu)化、SIMD 優(yōu)化寄存器分配等優(yōu)化方式，突破了快速傅里葉變換在ARMv8 與X86-64 硬件平臺(tái)上的算法性能，形成了一套FFT 算法在ARMv8 及X86 架構(gòu)上的實(shí)現(xiàn)策略和優(yōu)化方案，實(shí)現(xiàn)了一個(gè)跨平臺(tái)的高性能FFT 算法庫。同時(shí)對(duì)ARMv8 及X86-64 平臺(tái)上程序優(yōu)化提供了思路。下一步的工作將實(shí)現(xiàn)和優(yōu)化快速傅里葉變換列主序，完善OpenFFT 高性能算法庫，形成一套實(shí)用完善的高性能FFT 算法庫。