劉 仲，陳海燕，向宏衛(wèi)

(國防科技大學(xué)計(jì)算機(jī)學(xué)院，湖南長沙 410073)

現(xiàn)有的許多處理器體系結(jié)構(gòu)提供融合乘加(Fused Multiply-Add，F(xiàn)MA)指令，如 Intel的Itanium，IBM的Power處理器等。給定3個(gè)輸入操作數(shù)a，b，c，使用一條FMA指令即可實(shí)現(xiàn)y=±a±(b×c)中的任何一個(gè)操作。FMA指令對(duì)于數(shù)值計(jì)算亦具有重要的意義。因?yàn)椋谲浖?，一條FMA指令在執(zhí)行時(shí)間上和一條乘法或加法指令幾乎一樣快;在硬件上，F(xiàn)MA單元通常比分開的乘法器和加法器耗費(fèi)少;在計(jì)算上，F(xiàn)MA指令通過減少舍入操作可以提高數(shù)值計(jì)算精度。在提供FMA指令的處理器平臺(tái)上，對(duì)于有些數(shù)值計(jì)算應(yīng)用，轉(zhuǎn)換到FMA指令是比較直接的，如普通的矩陣和向量乘法、矩陣和矩陣乘法以及成對(duì)出現(xiàn)的乘法和加法操作計(jì)算。但是，對(duì)更多的其他應(yīng)用就沒有那么簡(jiǎn)單了，需要平衡加法和乘法操作，改造原有的計(jì)算方法或流程以適合FMA指令，才能充分發(fā)揮FMA指令的作用，加速計(jì)算的性能。

向量處理器在單個(gè)芯片上集成多個(gè)向量處理單元(Vector Processing Element，VPE)，每個(gè) VPE包含相同的MAC，ALU，BP等多個(gè)功能部件，能夠在降低面積和功耗的情況下大幅提高整芯片的計(jì)算能力，適合處理高密集運(yùn)算的任務(wù)，如矩陣計(jì)算、FFT運(yùn)算、濾波運(yùn)算等。然而，現(xiàn)有的大量程序和算法是基于單核處理器設(shè)計(jì)的，很多高密集運(yùn)算的任務(wù)由于算法本身的特性，向量化處理困難，如何針對(duì)向量處理器的多運(yùn)算部件的向量計(jì)算的體系結(jié)構(gòu)特點(diǎn)，充分開發(fā)各個(gè)層次的并行性，高效地向量化這些算法是當(dāng)前面臨的主要困難［1-2］。

快速傅里葉變換(Fast Fourier Transform，F(xiàn)FT)作為時(shí)域和頻域的基本變換工具，在現(xiàn)代信號(hào)處理系統(tǒng)領(lǐng)域應(yīng)用廣泛，如雷達(dá)、聲吶、地震信號(hào)分析、頻譜分析、語音識(shí)別、3G通信和圖像處理等。FFT能顯著減少離散傅里葉變換(Discrete Fourier Transform，DFT)計(jì)算的復(fù)雜度，對(duì)于實(shí)時(shí)性要求很高的信號(hào)處理應(yīng)用來說，DFT計(jì)算效率越高，信號(hào)處理的實(shí)時(shí)性就越好，因此FFT算法的優(yōu)化方法一直是國內(nèi)外的研究熱點(diǎn)。Pease［3］采用數(shù)學(xué)工具Kronecker積描述FFT算法，方便將FFT算法高效地映射到并行計(jì)算機(jī)上;Linzer［4］，Goedecker［5］和 Karner［6］等針對(duì) FMA 結(jié)構(gòu)的處理器，分別研究了如何利用FMA指令優(yōu)化FFT計(jì)算;Voronenko［7］依據(jù)有向無環(huán)圖和具有一定結(jié)構(gòu)的矩陣因子化算法，提出一種將離散傅里葉變換、離散余弦變換等線性變換轉(zhuǎn)化為FMA算法的一般方法;FFTW(Fastest Fourier Transform in theWest)［8］是通用CPU平臺(tái)上廣泛使用的FFT數(shù)學(xué)庫，移植性好，可以為任意長度和維數(shù)的FFT自動(dòng)生成一種最有效的實(shí)現(xiàn)方式;Loberiras［9］等針對(duì)GPU提出小點(diǎn)數(shù)FFT的優(yōu)化技術(shù)，受限于紋理存儲(chǔ)器的大小，大點(diǎn)數(shù) FFT的性能不理想;Li［10］等針對(duì) GPU提出一種 FFT自適應(yīng)優(yōu)化框架，當(dāng)點(diǎn)數(shù)大于1K時(shí)，獲得接近CUDA快速傅里葉變換(CUDA Fast Fourier Transform，CUFFT)的90%性能;文獻(xiàn)［11］針對(duì)GPU的大點(diǎn)數(shù)FFT優(yōu)化方法達(dá)到CUFFT性能的2．1倍。

在傳統(tǒng)的FFT計(jì)算中，蝶形單元的最終計(jì)算往往轉(zhuǎn)化為實(shí)數(shù)的乘法和加法操作，如時(shí)域抽取基2 FFT算法的一個(gè)蝶形單元計(jì)算需要4次實(shí)數(shù)乘法和6次實(shí)數(shù)加法，即需要10次實(shí)數(shù)乘(加)操作;時(shí)域抽取基4 FFT算法的一個(gè)蝶形單元計(jì)算需要12次實(shí)數(shù)乘法和22次實(shí)數(shù)加法，即需要34次實(shí)數(shù)乘(加)操作。顯然，F(xiàn)FT計(jì)算中的乘、加操作是不平衡的，不能直接利用FMA指令來實(shí)現(xiàn)蝶形單元的計(jì)算，所以傳統(tǒng)的FFT計(jì)算方法不能有效發(fā)揮具有FMA指令的處理器的計(jì)算效率。劉仲等針對(duì)FMA結(jié)構(gòu)的向量處理器，提出一種基于FMA加速FFT計(jì)算的向量化方法，通過變換FFT的蝶形單元運(yùn)算流程，將傳統(tǒng)計(jì)算方式中的獨(dú)立的乘法和加法操作組合成次數(shù)更少的FMA操作，使得計(jì)算一個(gè)DIT基2 FFT算法的蝶形單元所需要的浮點(diǎn)操作減少到6次FMA操作，計(jì)算一個(gè)DIT基4 FFT算法的蝶形單元所需要的浮點(diǎn)操作減少到24次FMA操作。實(shí)驗(yàn)結(jié)果表明，提出的方法能夠顯著加速FFT的計(jì)算性能。

1 向量處理器Matrix的體系結(jié)構(gòu)

如圖1所示，Matrix是一款面向高密度計(jì)算應(yīng)用的高性能浮點(diǎn)向量處理器，內(nèi)核是超長指令字(Very Long Instruction Word，VLIW)體系結(jié)構(gòu)，包括標(biāo)量處理部件(Scalar Processing Unit，SPU)和向量處理部件(Vector Processing Unit，VPU)。SPU負(fù)責(zé)標(biāo)量任務(wù)計(jì)算和流控，SPU和VPU可通過共享寄存器交換數(shù)據(jù)。Matrix每時(shí)鐘周期發(fā)射11條指令，包括5條標(biāo)量指令和6條向量指令。指令派發(fā)單元對(duì)執(zhí)行包進(jìn)行識(shí)別，并將其中的指令派發(fā)到相應(yīng)的功能單元中執(zhí)行。VPU負(fù)責(zé)向量計(jì)算，包括16個(gè)向量處理單元(Vector Processing Element，VPE)，每個(gè)VPE含一個(gè)局部寄存器文件，以及3個(gè)浮點(diǎn)乘加單元(Float Multiply and Accumulate，F(xiàn)MAC)、1個(gè)BP和2個(gè)L/S共6個(gè)并行功能部件，3個(gè)FMAC均支持FMA指令。局部寄存器文件包含64個(gè)64位寄存器，所有VPE的同一編號(hào)的局部寄存器在邏輯上又組成一個(gè)1024位的向量寄存器。功能部件支持定點(diǎn)和浮點(diǎn)操作，向量指令在各個(gè)VPE上同時(shí)獨(dú)立運(yùn)行。向量數(shù)據(jù)訪問單元支持向量數(shù)據(jù)的Load/Store，提供大容量陣列向量存儲(chǔ)器(Array Memory，AM)，每周期同時(shí)支持2個(gè)Load/Store指令。

圖1 Matrix的體系結(jié)構(gòu)Fig．1 Architecture of Matrix

2 基于FMA的FFT向量化優(yōu)化方法

序列x(n)(n=0，…，N-1)的離散傅里葉變換 X(k)(k=0，…，N-1)定義為:

依據(jù)離散傅里葉變換公式計(jì)算的DFT復(fù)數(shù)乘加運(yùn)算的計(jì)算復(fù)雜度是N2。當(dāng)N比較大時(shí)，運(yùn)算量增長得非常快，使得直接采用DFT計(jì)算方法難以滿足很多實(shí)際應(yīng)用的需求。1965年Cooley和Turkey提出一種快速傅里葉變換算法［12］，計(jì)算復(fù)雜度由原來的O(N2)降到O(N log2N)，可顯著地減少運(yùn)算量。FFT算法分為時(shí)間抽取法(Decimation In Time，DIT)和頻率抽取法(Decimation In Frequency，DIF)兩種，現(xiàn)以 DIT 方法為例闡述FFT的向量化方法。

2．1 傳統(tǒng)的FFT蝶形單元計(jì)算方法

2．1．1 DIT基2FFT的蝶形單元計(jì)算方法

當(dāng)N是2的整數(shù)次方時(shí)，DIT基2FFT將輸入數(shù)據(jù)序列x(n)進(jìn)行奇、偶分組:

由旋轉(zhuǎn)因子的周期性特性易知:

令 a(l)=x(2l)，b(l)=x(2l+1)，則序列X(k)劃分為2個(gè)長度為N/2的子序列:

圖2是DIT基2FFT的蝶形單元運(yùn)算流程圖，DIT基2FFT的每次蝶形單元運(yùn)算需要1次復(fù)數(shù)乘法，2次復(fù)數(shù)加法，轉(zhuǎn)變實(shí)數(shù)計(jì)算即為4次實(shí)數(shù)乘法和6次實(shí)數(shù)加法，即需要10次實(shí)數(shù)乘(加)操作。

圖2 DIT基2 FFT的蝶形單元運(yùn)算流程圖Fig．2 Radix-2 DIT FFT butterfly diagram

2．1．2 DIT基4FFT的蝶形單元計(jì)算方法

當(dāng)N是4的整數(shù)次方時(shí)，DIT基4FFT將輸入數(shù)據(jù)序列x(n)按模4后的余數(shù)分組:

由旋轉(zhuǎn)因子的周期性特性易知:

令 a(l)=x(4l)，b(l)=x(4l+1)，c(l)=x(4l+2)，d(l)=x(4l+3)，則序列 X(k)劃分為4個(gè)長度為N/4的子序列:

圖3是DIT基4FFT的蝶形單元運(yùn)算流程圖，DIT基4FFT的每次蝶形單元運(yùn)算需要3次復(fù)數(shù)乘法，8次復(fù)數(shù)加法，轉(zhuǎn)變實(shí)數(shù)計(jì)算即為12次實(shí)數(shù)乘法和22次實(shí)數(shù)加法，即需要34次實(shí)數(shù)乘(加)操作。

圖3 DIT基4 FFT的蝶形單元運(yùn)算流程圖Fig．3 Radix-4 DIT FFT butterfly diagram

2．2 FMA優(yōu)化的FFT蝶形單元運(yùn)算

2．2．1 FMA優(yōu)化的DIT基2FFT蝶形單元計(jì)算

假定DIT基2FFT的一個(gè)蝶形單元的兩個(gè)輸入分別為A和B，輸出分別為Y1和Y2，蝶形因子為W，下標(biāo)r和i分別表示復(fù)數(shù)的實(shí)部和虛部，根據(jù)式(1)，則有:

根據(jù)式(3)，DIT基2 FFT的蝶形單元運(yùn)算可以分解成兩步完成:

1)計(jì)算中間結(jié)果ˉA，ˉB:

2)計(jì)算輸出結(jié)果Y1和Y2:

其中，步驟1的計(jì)算由2條FMA指令完成，步驟2的計(jì)算由4條FMA指令完成。從而計(jì)算一個(gè)DIT基2 FFT蝶形單元僅需要6條FMA指令操作，相比傳統(tǒng)的計(jì)算方式減少了4次。在傳統(tǒng)的計(jì)算方法下，計(jì)算N點(diǎn)的DIT基2 FFT需要的浮點(diǎn)操作次數(shù)為5N log2N;FMA優(yōu)化后需要的浮點(diǎn)操作次數(shù)減少到3N log2N，減少了40%。

2．2．2 FMA優(yōu)化的DIT基4FFT蝶形單元計(jì)算

假定DIT基4 FFT的一個(gè)蝶形單元的4個(gè)輸入分別為 A，B，C，D，輸出分別為 Y1，Y2，Y3，Y4，3個(gè)蝶形因子分別為W1，W2，W3，下標(biāo)r和i分別表示復(fù)數(shù)的實(shí)部和虛部，根據(jù)式(2)，則有:

根據(jù)式(4)，DIT基4 FFT的蝶形單元運(yùn)算可以分解成兩步完成:

1)計(jì)算中間結(jié)果ˉA，ˉB，ˉC，ˉD:

其中，步驟1的計(jì)算由12條FMA指令完成，步驟2的計(jì)算由12條FMA指令完成。從而計(jì)算一個(gè)DIT基4FFT蝶形單元僅需要24條FMA指令操作，相比傳統(tǒng)的計(jì)算方式減少了10次。在傳統(tǒng)的計(jì)算方法下，計(jì)算N點(diǎn)的DIT基4 FFT需要的浮點(diǎn)操作次數(shù)為8．5N log4N;FMA優(yōu)化后需要的浮點(diǎn)操作次數(shù)減少到6N log4N，減少了29．4%。

2．3 混合基4和基2的FFT計(jì)算方法

如圖4所示，對(duì)于任意的N=2n，可采用混合基4和基2的FFT方法加速FFT的計(jì)算。若n是偶數(shù)，令n=2m，則N=22m=4m;若n是奇數(shù)，令n=2m+1，則 N=22m+1=4m×2。因此，N 點(diǎn)的FFT計(jì)算可轉(zhuǎn)化為m級(jí)基4FFT或者m級(jí)基4FFT和最后一級(jí)的基2 FFT計(jì)算。

圖4 混合基4和基2的FFT計(jì)算流程Fig．4 FFT computation diagram ofmixed radix-2/4

2．4 旋轉(zhuǎn)因子向量訪問的優(yōu)化

傳統(tǒng)的基4 FFT計(jì)算方法中，每個(gè)蝶形單元運(yùn)算需要乘以3個(gè)不同的旋轉(zhuǎn)因子:和FMA優(yōu)化后的蝶形單元運(yùn)算只需計(jì)算2個(gè)不同旋轉(zhuǎn)因子:WiN和W2Ni，并且預(yù)先計(jì)算出旋轉(zhuǎn)因子虛部/實(shí)部(Wi/Wr)的值，然后把旋轉(zhuǎn)因子按實(shí)部(Wr)、虛部/實(shí)部(Wi/Wr)交叉存儲(chǔ)放置。這樣，能減少旋轉(zhuǎn)因子存儲(chǔ)的個(gè)數(shù)。由于向量處理器只能連續(xù)存取向量數(shù)據(jù)，所以每級(jí)FFT運(yùn)算的不同旋轉(zhuǎn)因子都要單獨(dú)放置。

同樣，DIT基2 FFT算法需要預(yù)先計(jì)算旋轉(zhuǎn)因子WiN和虛部/實(shí)部(Wi/Wr)的值，并且按實(shí)部(Wr)、虛部/實(shí)部(Wi/Wr)交叉的方式存放。

以1024點(diǎn)的DIT基4 FFT為例，算法由5級(jí)蝶形運(yùn)算實(shí)現(xiàn)，且每一級(jí)不同旋轉(zhuǎn)因子所需的數(shù)目分別為:1個(gè)、4個(gè)、16個(gè)、64個(gè)和256個(gè)。由此可知，第1級(jí)和第2級(jí)的不同旋轉(zhuǎn)因子的數(shù)目少于VPE的個(gè)數(shù)，并且第1級(jí)需要乘的旋轉(zhuǎn)因子的值為1，所以可以省略。為了提高運(yùn)算的并行度，可以把第2級(jí)所需的旋轉(zhuǎn)因子進(jìn)行冗余存放4次。

表1 1024點(diǎn)DIT基4 FFT的旋轉(zhuǎn)因子Tab．1 Twiddle factors of 1024-point radix-4 DIT FFT

3 性能測(cè)試與分析

在向量處理器Matrix上對(duì)不同點(diǎn)數(shù)的FFT計(jì)算性能進(jìn)行了測(cè)試(稱 MatrixFFT)，并與 Intel，AMD，IBM 平臺(tái)上的 MKL，ACML，ESSL 和 FFTW 3．0算法庫的性能進(jìn)行比較。實(shí)驗(yàn)平臺(tái)中，Matrix的主頻為1GHz(雙精度峰值性能為96GFLOPS)，MatrixFFT的測(cè)試結(jié)果為RTL級(jí)仿真環(huán)境下的測(cè)試性能數(shù)據(jù);比對(duì)實(shí)驗(yàn)的 CPU中，Intel選擇主頻3．0GHz Intel Xeon Core Duo(Woodcrest)(峰值性能為12GFLOPS);AMD選擇主頻2．2 GHz Dual Core AMD Opteron(峰值性能為4．4GFLOPS);IBM選擇主頻2GHz PowerPC 970(峰值性能為10GFLOPS)，其性能數(shù)據(jù)選自FFTW網(wǎng)站報(bào)告的評(píng)測(cè)數(shù)據(jù)［8］。

圖5和圖6分別給出了在Matrix上測(cè)試的不同點(diǎn)數(shù)的基2和基4FFT的單精度和雙精度計(jì)算性能。從圖5、圖6中可以看出，在點(diǎn)數(shù)較小時(shí)，F(xiàn)FT的性能較低，隨著點(diǎn)數(shù)的增大，F(xiàn)FT的計(jì)算性能顯著提高，16K點(diǎn)的單精度基2和基4FFT性能分別達(dá)到74GFLOPS和98GFLOPS，計(jì)算效率分別為38．5%和51．04%。16K點(diǎn)的雙精度基2和基4FFT 性能分別達(dá)到48．83GFLOPS 和49．05GFLOPS，計(jì)算效率分別為50．08%和51．09%。

圖5 基2FFT的計(jì)算性能Fig．5 Performance of radix-2 FFT

圖6 基4FFT的計(jì)算性能Fig．6 Performance of radix-4 FFT

圖7 對(duì)比了不同點(diǎn)數(shù)的雙精度FFT在 FMA優(yōu)化前后的計(jì)算性能，從圖7中可以看出，無論是DIT基2FFT，還是 DIT基4FFT，經(jīng)過 FMA優(yōu)化后，F(xiàn)FT的計(jì)算性能均顯著提高，其中 DIT基2FFT的計(jì)算性能平均提高30．4%，DIT基4FFT的計(jì)算性能平均提高30．3%。

圖7 FMA優(yōu)化的FFT性能對(duì)比Fig．7 Performance comparison of FMA optimized FFT

表2 不同處理器平臺(tái)下FFT性能(MFLOPS)Tab．2 FFT performance of different processors(MFLOPS)

表3 不同處理器平臺(tái)下FFT計(jì)算效率Tab．3 FFT efficiency of different processors

表2和表3分別從計(jì)算性能和效率兩方面給出了不同點(diǎn)數(shù)的雙精度FFT分別在Matrix，Intel，AMD，IBM 平 臺(tái) 上 的 MKL，ACML，ESSL 和FFTW3．0算法庫測(cè)試結(jié)果。

從絕對(duì)計(jì)算性能上看，MatrixFFT的性能遠(yuǎn)超其他幾種算法庫，主要原因是Matrix的峰值性能遠(yuǎn)超對(duì)比CPU的峰值性能;另一方面也體現(xiàn)出提出的基于FMA優(yōu)化的FFT計(jì)算效率較高，能夠充分發(fā)揮Matrix的計(jì)算性能。這一點(diǎn)從表3可以看出，在點(diǎn)數(shù)較小時(shí)，MatrixFFT的計(jì)算效率較低，這是因?yàn)镸atrix是向量處理器，點(diǎn)數(shù)較小時(shí)，軟件流水中的循環(huán)填充和排空的開銷在整個(gè)計(jì)算中的占比較高，影響計(jì)算效率。在點(diǎn)數(shù)超過1024點(diǎn)以后，MatrixFFT的計(jì)算效率顯著提高，在4K點(diǎn)以后，計(jì)算效率都是最高的，其中16K點(diǎn)效率達(dá)50．86%，而對(duì)比的其他算法庫效率分別是47．88%，36．69%，26．83%，29．27%，12．28%。

4 結(jié)論

本文針對(duì) FMA結(jié)構(gòu)的向量處理器，提出FMA加速FFT計(jì)算的向量化方法。通過優(yōu)化和重組FFT算法的蝶形單元運(yùn)算流程，將原本不平衡的乘法和加法操作組合成融合乘加操作，利用FMA指令減少FFT計(jì)算的浮點(diǎn)操作指令次數(shù)，進(jìn)而提高FFT算法的計(jì)算性能和向量處理器的計(jì)算效率，提高硬件資源的利用率，加速FFT算法的計(jì)算性能。

References)

［1］ 劉仲，陳躍躍，陳海燕．支持任意系數(shù)長度和數(shù)據(jù)類型的FIR濾波器向量化方法［J］．電子學(xué)報(bào)，2013，41(2):346-351．LIU Zhong，CHEN Yueyue，CHEN Haiyan．A vectorization of FIR filter supporting any length and data types of coefficients［J］．Acta Electronica Sinica，2013，41(2):346-351．(in Chinese)

［2］ 劉仲，邢彬朝，陳躍躍．一種面向多核處理器的高效并行PCA-SIFT算法［J］．國防科技大學(xué)學(xué)報(bào)，2012，34(4):103-107．LIU Zhong，XING Binchao，CHEN Yueyue．An efficient parallel PCA-SIFT algorithm for multi-core processor［J］．Journal of National University of Defense Technology，2012，34(4):103-107．(in Chinese)

［3］ Pease M C．An adaptation of the fast Fourier transform for parallel processing［J］．Journal of the ACM，1968，15(2):252-264．

［4］ Linzer EN，F(xiàn)eig E．Implementation ofefficient FFT algorithms on fused multiply-add architectures［J］．IEEE Transactions on Signal Processing，1993，41(1):93-107．

［5］ Goedecker S．Fast radix 2，3，4，and 5 kernels for fast Fourier transformations on computers with overlapping multiply-add instructions［J］．SIAM Journal on Scientific Computing，1997，18(6):1605-1611．

［6］ Karner H，Auer M，Ueberhuber CW．Multiply-add optimized FFT kernels［J］．MathematicalModels and Methods in Applied Sciences，2001，11(1):105-117．

［7］ Voronenko Y，Puschel M．Mechanical derivation of fused multiply-add algorithms for linear transforms［J］． IEEE Transactions on Signal Processing，2007，55(9):4458-4473．

［8］ Frigo M，Johnson S G．BenchFFT［EB/OL］．［2014-03-15］．http://www．fftw．org/benchfft/．

［9］ Lobeiras J，Amor M，Doallo R．Influence of memory access patterns to small-scale FFT performance［J］． Journal of Supercomputing，2013，64(1):120-131．

［10］ Li Y，Zhang Y Q，Liu Y Q，et al．MPFFT:An auto-tuning FFT library for OpenCL GPUs［J］．Journal of Computer Science and Technology，2013，28(1):90-105．

［11］ 何濤，朱岱寅．大點(diǎn)數(shù)一維FFT的GPU設(shè)計(jì)實(shí)現(xiàn)［J］．計(jì)算機(jī)工程與科學(xué)，2013，35(11):34-41．HE Tao，ZHU Daiyin．Design and implementation of largepoint1D FFT on GPU［J］．Computer Engineering ＆ Science，2013，35(11):34-41．

［12］ Cooley JW，Turkey JW．An algorithm for the machine calculation of complex Fourier series［J］．Mathematics of Computation，1965，19:297-301．