李 琨，賈海鵬，曹 婷，張云泉

1.中國科學院 計算技術(shù)研究所 計算機體系結(jié)構(gòu)國家重點實驗室，北京 100190

2.中國科學院大學 計算機與控制學院，北京 100190

大規(guī)模集群上多維FFT算法的實現(xiàn)與優(yōu)化研究*

李 琨1,2+，賈海鵬1，曹 婷1，張云泉1

1.中國科學院 計算技術(shù)研究所 計算機體系結(jié)構(gòu)國家重點實驗室，北京 100190

2.中國科學院大學 計算機與控制學院，北京 100190

LI Kun,JIA Haipeng,CAO Ting,et al.Implementation and optimization of multidimensional FFT algorithm on large-scale clusters.Journal of Frontiers of Computer Science and Technology,2017,11(6)：863-874.

快速傅里葉變換（fast Fourier transform，F(xiàn)FT）是用于計算離散傅里葉變換（discrete Fourier transform，DFT）或其逆運算的快速算法，在工程、科學和數(shù)學領(lǐng)域的應(yīng)用非常廣泛，例如信號分解、數(shù)字濾波、圖像處理等。因此，在實際應(yīng)用中對FFT算法進行細粒度優(yōu)化是非常重要的。研究了FFT算法常用的分解策略以及FFT算法在大規(guī)模集群系統(tǒng)上的并行實現(xiàn)，并提出了相關(guān)的優(yōu)化策略。在此基礎(chǔ)上，對多種FFT算法在不同平臺上進行了性能評估，并分析了各算法的實現(xiàn)、優(yōu)缺點及其在大規(guī)模計算時的可擴展性。實驗結(jié)果表明，相關(guān)研究有助于對現(xiàn)有的FFT算法進行進一步的優(yōu)化，以及指導如何在大規(guī)模CPU+GPU的異構(gòu)系統(tǒng)上根據(jù)不同需求選擇實現(xiàn)性能更優(yōu)的FFT算法。

集群；快速傅里葉變換（FFT）；消息傳遞接口（MPI）；性能優(yōu)化

1 引言

快速傅里葉變換（fast Fourier transform，F(xiàn)FT），是用于計算序列的離散傅里葉變換（discrete Fourier transform，DFT）或其逆變換的一種高效算法。FFT被廣泛應(yīng)用于工程、科學和數(shù)學領(lǐng)域，如在國際大科學工程平方公里陣列射電望遠鏡（square kilometer array，SKA）[1]的數(shù)據(jù)處理中，F(xiàn)FT占總計算量的40%，且由于其每秒TB級的數(shù)據(jù)量及底層高性能計算機系統(tǒng)的規(guī)模和復雜度，對FFT算法的計算效能提出了更高的要求，需針對性地研究相應(yīng)的方法。鑒于此，本文簡要介紹了FFT的經(jīng)典算法——庫利圖基算法，針對集群系統(tǒng)研究了并行FFT算法現(xiàn)存的瓶頸以及相應(yīng)的優(yōu)化策略，最后完整地將FFTW庫（faster Fourier transform in the West）、PFFT庫（parallel FFT library）與MPFFT庫（massive parallel FFT library）3種算法在不同平臺上進行了性能評估與分析。

本文的主要創(chuàng)新點如下：

（1）分析了PFFT和MPFFT算法的優(yōu)化現(xiàn)狀，指出目前大規(guī)模集群上FFT算法存在的兩點可優(yōu)化方向，一是數(shù)據(jù)在主存和設(shè)備內(nèi)存間的不斷通信，使得PCI-E總線帶寬成為了性能瓶頸；二是在數(shù)據(jù)輸出的同時可以進行細節(jié)上的處理，以更細致的劃分方式避免最后數(shù)據(jù)的整體調(diào)整。

（2）在深入分析目前主流的幾種FFT算法基礎(chǔ)上，研究了一維分解、二維分解的優(yōu)化策略，并提出了通信優(yōu)化和輸出優(yōu)化兩種改進方案。通信優(yōu)化將數(shù)據(jù)在CPU和GPU間的傳輸進行了更細致的分解，將有效地提高整體的通信效率，減少全局轉(zhuǎn)置所需的時間；輸出優(yōu)化在每一塊數(shù)據(jù)進行FFT計算的同時確定好最優(yōu)的內(nèi)存分布位置，便于對計算后的數(shù)據(jù)進行排列，如此可使得FFT計算的同時進行數(shù)據(jù)的整理，減少額外進行轉(zhuǎn)置的開銷。

（3）在單節(jié)點平臺和天河一號、天河二號超級計算機上通過對各種規(guī)模輸入數(shù)據(jù)進行大量的測試和統(tǒng)計，詳細地評估與分析了3種主流FFT算法的性能、可擴展性以及各算法的最佳適應(yīng)環(huán)境，為FFT在大規(guī)模計算時的算法選擇提供了參考。實驗表明，在大多數(shù)情況下PFFT算法與MPFFT算法要優(yōu)于FFTW算法，且當輸入數(shù)據(jù)量較小時，MPFFT算法性能優(yōu)勢相對明顯；輸入數(shù)據(jù)量較大時，PFFT算法性能優(yōu)勢則相對明顯。另外，這種性能優(yōu)勢會隨著進程數(shù)量的增加而更為突出，當進程數(shù)量較小時，3種算法的性能優(yōu)劣對比并不明顯。

本文組織結(jié)構(gòu)如下：第2章介紹了相關(guān)工作；第3章介紹了FFT的庫利圖基算法；第4章對FFTW、PFFT、MPFFT算法的分解策略和實現(xiàn)方法進行了分析，針對3種算法的性能瓶頸提出了兩種優(yōu)化策略；第5章在單節(jié)點和大規(guī)模集群上使用不同的輸入規(guī)模對3種算法進行了性能評估和性能分析；第6章進行了總結(jié)討論。

2 相關(guān)工作

目前常用的FFT庫有FFTW[2]、PFFT[3]、MPFFT[4]、CUFFT[5]（CUDA fast Fourier transform library）、PKUFFT[6-7]（Peking University fast Fourier transform library）等。Frigo和Johnson開發(fā)的FFTW庫可以計算一維或者多維實數(shù)據(jù)、復數(shù)據(jù)的離散傅里葉變換，其具備自適應(yīng)性，并為各種領(lǐng)域內(nèi)所需的FFT計算提供了基本的運行框架。PFFT算法是在FFTW庫上的拓展，其加強了FFTW軟件庫的MPI（message passing interface）部分以適于多維的數(shù)據(jù)分解，即規(guī)模為Nd的d維FFT可以至多被Nd-1個MPI進程并行計算[8]。中國科學院李焱等人提出了FFT針對GPU平臺的自適應(yīng)性能優(yōu)化框架和針對CPU和GPU異構(gòu)集群的MPFFT軟件包，融合了異構(gòu)平臺的硬件異同和FFT算法的特性，使得FFT在這些平臺都能取得良好的加速效果[4]。CUFFT是NVIDIA SDK中提供的FFT庫，它的API接口與FFTW比較相像，包括FFT plan的搜索階段和執(zhí)行階段。另外陳一峯等人通過將數(shù)據(jù)各維度進行細粒度劃分，利用GPU通信緩沖區(qū)與GPU Pinned Memory之間拷貝的機會對數(shù)據(jù)進行粒度上的調(diào)整來優(yōu)化AlltoAll通信，使得PKUFFT取得了較高的性能[6]。

與以上工作不同，本文分析了各種庫的并行分解策略以及各自的優(yōu)勢和不足，從通信、輸出等多方面對FFTW、PFFT、MPFFT算法提出了可能的改進方案，并對3種算法在大規(guī)模集群上的性能進行了全面的評測和分析，所得結(jié)果可為3種算法今后的優(yōu)化方向及在大規(guī)模計算時FFT算法的選擇提供參考。

3 背景介紹

3.1 快速傅里葉變換

離散傅里葉變換在計算機科學、物理學等學科中扮演著非常關(guān)鍵的角色。對于長度為N的復數(shù)序列x，其中x=x0,x1,…,xn-1，DFT的計算公式[9]為：

根據(jù)式（1），可以發(fā)現(xiàn)Xk共有N次輸出結(jié)果，每一次的輸出都有N次的累加求和運算，若按照DFT的定義進行計算求值需O(n2)次運算。在實際應(yīng)用中，F(xiàn)FT算法通常被用于DFT變換。原因有兩點：一是FFT算法在求解過程中會有舍入誤差的存在，因此許多FFT算法運行后的結(jié)果甚至比用定義進行求解的結(jié)果要精確很多；二是FFT算法的運算速度更快，所有的FFT算法的復雜度均可在對數(shù)時間復雜度內(nèi)完成[10]。FFT已經(jīng)有很多成熟的算法，在所有的FFT算法中，庫利圖基算法應(yīng)用最廣、最流行，是眾多算法庫實現(xiàn)的基礎(chǔ)，故下文將簡要地介紹庫利圖基算法。

3.2 庫利圖基快速傅里葉變換算法

庫利圖基快速傅里葉變換算法（Cooley-Tukey算法）[11]是目前應(yīng)用最廣泛、最流行的FFT算法。利用式（1），將公式中的項目在時域上進行重新分組，并將進行替換，其中替換后的被稱為“旋轉(zhuǎn)因子”（twiddle factor），亦稱為“蝶形因子”。根據(jù)旋轉(zhuǎn)因子在計算過程中出現(xiàn)的不同位置，可以將FFT算法分為時域抽?。╠ecimation-in-time，DIT）和頻域抽?。╠ecimation-in-frequency，DIF）兩大類。DIT的旋轉(zhuǎn)因子出現(xiàn)在計算的輸入端，而DIF的旋轉(zhuǎn)因子則出現(xiàn)在計算的輸出端。

以DIT為例，將旋轉(zhuǎn)因子替換后，DIT是將長度為N的復數(shù)序列x根據(jù)數(shù)據(jù)的偶數(shù)項和奇數(shù)項進行分解，分解后的式（1）將成為偶序列與奇序列兩部分。其中偶序列為x0,x2,x4,…,xn-2，奇序列為x1,x3, x5,…,xn-1，運算如式（2）：

其中k=0,1,…,N-1。

因為Gk與Hk具有周期性，WN具有對稱性和周期性，所以可得：

按照同樣的方式，對式（3）中的Gk與Hk也繼續(xù)以相同的方式分解，則可將一個N點的離散傅里葉變換最終用一組兩點的DFT來計算。在基數(shù)為2的快速傅里葉變換中，共有l(wèi)bN級運算，而且在每一級的計算中有(N/2)個兩點FFT蝶形運算[12]。因此采用該FFT算法的復雜度將直接降為O(nlbn)。

4 并行FFT優(yōu)化分析及策略

4.1 優(yōu)化分析

對于許多應(yīng)用，F(xiàn)FT算法的輸入并不能完全只在一個單節(jié)點上進行，因此對輸入數(shù)據(jù)勢必要劃分到各個節(jié)點上運算以提高效率。每個子任務(wù)計算自己的FFT部分，之后再將數(shù)據(jù)與其他子任務(wù)進行交換。實現(xiàn)該操作的策略有一維分解和二維分解。

對于異構(gòu)平臺上的應(yīng)用，一個比較普遍的問題是哪一種任務(wù)適于使用GPU設(shè)備進行加速。GPU適合處理計算密集型或者訪存密集型的任務(wù)，其他任務(wù)，例如Internet帶寬，存儲設(shè)備的讀、寫帶寬以及通信延遲等都是不能被GPU加速的。

對于小規(guī)模的FFT，如果數(shù)據(jù)能夠被全部地置于GPU設(shè)備上，那么此時的計算就會從設(shè)備內(nèi)存的高帶寬中受益。這種方案是有特定情景的，即大部分的數(shù)據(jù)都置于設(shè)備的內(nèi)存上，并且FFT在這些數(shù)據(jù)上運算了多次。這時數(shù)據(jù)在主存和設(shè)備內(nèi)存間的通信開銷與GPU的計算時間相比就可以忽略了。

然而，對于大規(guī)模的FFT，數(shù)據(jù)需要在主存和設(shè)備內(nèi)存間不斷通信，PCI-E總線帶寬就會成為瓶頸。因為此時通過PCI-E總線的數(shù)據(jù)傳輸時間要比在GPU的計算時間長得多。這種方案也是有特定情景的，即針對的是在節(jié)點上計算大規(guī)模FFT，所有數(shù)據(jù)必須在不同的節(jié)點之間、在主存和設(shè)備內(nèi)存之間來回傳輸?shù)臅r候，如果能對各節(jié)點間的通信、GPU與CPU間的通信進行優(yōu)化，將會大幅提高整體的計算效率。

4.2 優(yōu)化策略

4.2.1 一維分解

一維分解是基于轉(zhuǎn)置的并行FFT分解策略，該分解策略已在FFTW算法庫[8]中實現(xiàn)。以三維輸入數(shù)據(jù)為例，假設(shè)三維的輸入數(shù)據(jù)大小n0×n1×n2，其中n0≥n1≥n2。假定MPI進程的數(shù)量為P，則首先沿著n0方向進行分解，將數(shù)據(jù)切分到P個MPI進程上，這P個進程并行地計算n0/P組大小為n1×n2的二維FFT。此時，n0維的數(shù)據(jù)已被切分分布于所有的P個進程上，故n0維數(shù)據(jù)的計算需使用AlltoAll的MPI函數(shù)對數(shù)據(jù)交換轉(zhuǎn)置以完成n0維數(shù)據(jù)的計算[2]。

MPI_ALLTOALL函數(shù)是對MPI_ALLGATHER的擴展，兩者的區(qū)別是調(diào)用MPI_ALLTOALL時每個接收者可接收來自別的發(fā)送者的數(shù)目不同的數(shù)據(jù)。例如，第m個進程發(fā)送的第n塊數(shù)據(jù)將被第n個進程接收并存放在其接收消息緩沖區(qū)recv_buffer的第m塊上。

一維分解較為簡單，然而采用該策略進行實現(xiàn)的庫在可擴展性方面存在很大的不足，它要求P=n0，n0=max{n0,n1,n2}。當P>n0時，一維分解將不能被使用。另外當P>n1或P>n2時，數(shù)據(jù)分解所使用的進程數(shù)P將受限于n1或n2。數(shù)據(jù)轉(zhuǎn)置之后，多余的進程數(shù)將會空閑得不到利用。因此，對數(shù)據(jù)的劃分需要采用更多維來進行計算[4]。圖1為一維分解示意圖。

Fig.1 One-dimension decomposition strategy圖1 一維分解策略

4.2.2 二維分解

二維分解[3]相較于一維分解很好地解決了可擴展性不足這一嚴重的弊端，故在大規(guī)模集群上也應(yīng)用得比較普遍。圖2即為二維分解的效果圖。該分解策略首先由Ding等人提出，后來Eleftheriou等人引入了體區(qū)域分解（volumetric domain decomposition）的概念，其在本質(zhì)上使用的還是二維分解，即將輸入數(shù)組首先沿n0和n1進行分解，因此這種分解策略允許增加的MPI進程數(shù)目最多可達n0×n1。因為輸入的數(shù)據(jù)僅有一維數(shù)據(jù)位于本地進程中，所以每次變換需進行多次的數(shù)據(jù)交換轉(zhuǎn)置，包括本地的數(shù)據(jù)轉(zhuǎn)置和全局的數(shù)據(jù)轉(zhuǎn)置。

Fig.2 Tow-dimension decomposition strategy圖2 二維分解策略

對于規(guī)模為n0×n1×n2的輸入，二維分解會將其并行地分布于P=P0×P1個進程上，每個進程首先會得到大小為n2的二維條，圖3表示的就是將輸入數(shù)據(jù)分布到4×4的進程網(wǎng)格上。為了計算三維FFT，算法將沿著n2維進行一批FFT計算，之后進行全局轉(zhuǎn)置操作，以使得每個進程在本地得到第二維的數(shù)據(jù)。類似的操作會再次進行以得到第一維的數(shù)據(jù)，最終的輸出分布形式為此時輸出的內(nèi)存排列與輸入時是不一致的，差別為輸入數(shù)組沿著第一和第二個維度分布，而輸出數(shù)組則沿著第二和第三維度分布，如圖4所示。PFFT、MPFFT算法在此之后會在輸出數(shù)組上進行內(nèi)存布局的調(diào)整，以使最后的輸出維度的分布也與輸入時一致。

Fig.3 Input data distributed on 16 processes initially(4×4 process grid)圖3 初始時輸入數(shù)據(jù)分布于16個進程上（4×4進程網(wǎng)格）

Fig.4 Input data distributed on 16 processes at the end(4×4 process grid)圖4 結(jié)束時輸入數(shù)據(jù)分布于16個進程上（4×4進程網(wǎng)格）

4.2.3 通信優(yōu)化

基于上面的分解策略，F(xiàn)FTW庫采用了一維分解，PFFT和MPFFT庫采用了二維分解。FFTW和PFFT算法未進行通信方面的優(yōu)化。MPFFT算法在每維FFT計算時直接調(diào)用NVIDIA CUFFT庫，之后進行本地轉(zhuǎn)置，由于此時轉(zhuǎn)置所需的數(shù)據(jù)全部位于GPU中，故其是借助GPU上的矩陣轉(zhuǎn)置函數(shù)來提高性能。之后與CPU通信時的全局轉(zhuǎn)置是通過封裝FFTW中相關(guān)接口直接實現(xiàn)，并未進行更進一步的優(yōu)化。而在大規(guī)模計算FFT時，顯然這將會成為一個限制效率的瓶頸。

基于此可以考慮將需要傳送拷貝的數(shù)據(jù)劃分為相同大小的數(shù)據(jù)塊[13]。通過異步通信的方式，一次將一個數(shù)據(jù)塊拷貝到CPU內(nèi)存中，此時CPU間的通信便可以開始，同時其他的數(shù)據(jù)塊也以異步的方式傳輸?shù)紺PU內(nèi)存。也可以調(diào)用異步接收指令，這樣數(shù)據(jù)接收操作與數(shù)據(jù)塊從設(shè)備向主存的拷貝操作同時進行。同時，每個接收的數(shù)據(jù)塊之后可以異步地拷回到GPU設(shè)備中，這樣便可有效地提高通信效率，減少全局轉(zhuǎn)置所需的時間。圖5展示了對于P=4時all-to-all的過程。4×4矩陣中的每個元素代表了存儲的數(shù)據(jù)塊。表示從進程i到進程j的一個發(fā)送指令，表示進程j接收來自進程i的消息。

Fig.5 Communication optimization圖5 通信優(yōu)化示意圖

4.2.4 輸出優(yōu)化

假設(shè)用T0代表本地轉(zhuǎn)置，用T1代表全局轉(zhuǎn)置，以0代表經(jīng)FFT計算后的n0。在運算過程中，初始的輸入分布可表示為采用二維分解策略，經(jīng)過一系列計算后得到的輸出分布為此后，一部分FFT算法，例如PFFT、MPFFT算法，將在此輸出形式上進行多次轉(zhuǎn)置計算以改變其內(nèi)存布局[14-15]，轉(zhuǎn)換過程可表示如下：

經(jīng)過多次的本地轉(zhuǎn)置和全局轉(zhuǎn)置，可以將輸出的分布形式完全轉(zhuǎn)化為和初始的輸入分布一致，即為然而這種轉(zhuǎn)化是在已經(jīng)得到正確的輸出結(jié)果之上所做的進一步的整理工作，這種額外的數(shù)據(jù)整理又會增加一些計算和通信上的開銷，而這種開銷理論上也可以進行進一步的分析優(yōu)化。在不可避免使用全局轉(zhuǎn)置的地方盡量以塊劃分的異步通信策略進行優(yōu)化，而在每一維度上進行轉(zhuǎn)置來計算FFT的時候可以進行更為細致的粒度劃分，即在每一塊數(shù)據(jù)進行FFT計算的同時確定好最優(yōu)的內(nèi)存分布位置，便于對計算后的數(shù)據(jù)進行排列。這樣，便可以在FFT計算的同時進行數(shù)據(jù)整理，減少額外進行轉(zhuǎn)置的開銷。另外，在某些不需要得到和輸入相同數(shù)據(jù)分布的輸出的情況下，對輸出數(shù)據(jù)的整理就是一步冗余操作。此時可以直接取消輸出轉(zhuǎn)化的過程，輸出計算后的結(jié)果或按實際需求相應(yīng)地進行小規(guī)模的數(shù)據(jù)整理，這也可以減少算法運行的時間。

5 性能評估

本文對FFTW、PFFT和MPFFT庫在異構(gòu)單節(jié)點、天河一號和天河二號3種平臺上的性能和可擴展性進行評估，其中FFTW庫的版本為3.3.4，具體測試評估將分別進行介紹。

5.1 單節(jié)點上算法性能分析

本節(jié)將在固定的異構(gòu)平臺A上分別對FFTW、PFFT以及MPFFT算法進行常規(guī)測試。

本次實驗平臺A為CPU和GPU的異構(gòu)平臺，其中CPU配置的詳細信息如表1所示。平臺上共有3塊GPU，分別是一塊GeForce GTX Titan Black，內(nèi)存為6 GB；兩塊Tesla K80，內(nèi)存均為11 GB。使用的CUDA版本為7.5，且運行時調(diào)用了3塊GPU進行計算。

Table 1 CPU configuration表1CPU配置信息

實驗分為兩種輸入：一是1 024×1 024×1 024的規(guī)模下，采用FFTW和PFFT算法得到的性能；二是輸入規(guī)模為256×256×256，采用FFTW、PFFT和MPFFT算法得到的性能。具體實驗所得數(shù)據(jù)如圖6~圖10所示。

圖6和圖7為1 024×1 024×1 024的輸入時，采用FFTW和PFFT算法得到的性能。從中可以看出PFFT比FFTW算法可擴展性更強，當輸入為1 024×1 024× 1 024時，其性能一直優(yōu)于FFTW算法。隨著進程數(shù)量的增加，算法性能也在增長，但由于該平臺是單節(jié)點服務(wù)器，該測試并未能真正體現(xiàn)出在大規(guī)模集群上的運行情況。另外平臺中GPU內(nèi)存大小的限制使得MPFFT算法不能夠接收如此大規(guī)模的輸入。

針對GPU內(nèi)存大小的限制問題，本文使用256× 256×256的輸入在該平臺上進行新的測試，圖8和圖9即為測試結(jié)果。整體上看，3種算法的運算性能大致為MPFFT優(yōu)于PFFT優(yōu)于FFTW算法，且3種算法的運算性能隨進程數(shù)量的增加都在提升。更確切地說，在并行的進程數(shù)小于8時，使用GPU的MPFFT異構(gòu)算法和使用并行計算的PFFT算法要優(yōu)于FFTW算法，而隨著進程數(shù)量的增加，兩種算法相對于FFTW算法的優(yōu)勢則不再明顯。主要原因可歸為以下3點：

（1）隨著進程數(shù)量的增加，并行算法中進程間的通信和數(shù)據(jù)交換將明顯地影響到算法的性能。

Fig.6 Running time comparison of different algorithms when input size is 1 024×1 024×1 024(PlatformA)圖6 數(shù)據(jù)規(guī)模為1 024×1 024×1 024時不同算法運行時間對比（平臺A）

Fig.8 Running time comparison of different algorithms when input size is 256×256×256(PlatformA)圖8 數(shù)據(jù)規(guī)模為256×256×256時不同算法運行時間對比（平臺A）

Fig.7 Performance comparison of different algorithms when input size is 1 024×1 024×1 024(PlatformA)圖7 數(shù)據(jù)規(guī)模為1 024×1 024×1 024時不同算法運行性能對比（平臺A）

Fig.9 Performance comparison of different algorithms when input size is 256×256×256(PlatformA)圖9 數(shù)據(jù)規(guī)模為256×256×256時不同算法運行性能對比（平臺A）

（2）PFFT與MPFFT算法針對的都是集群優(yōu)化的算法，平臺A只是一個單節(jié)點服務(wù)器，故當進程數(shù)量、輸入量較小時，在并行和GPU方面進行優(yōu)化的MPFFT算法與PFFT算法會大幅地改善性能。

（3）當進程數(shù)量逐漸增加，此時測試的數(shù)據(jù)與之前的大規(guī)模輸入相比較小，因此計算任務(wù)比之前要小，在這種情況下GPU的加速性能并沒有充分利用，使得優(yōu)化效果減弱，且較小規(guī)模的計算任務(wù)時PFFT算法的可擴展性將不如FFTW算法，從而當進程增加時MPFFT算法與PFFT算法的總體性能將會出現(xiàn)暫時不如FFTW算法的現(xiàn)象。

圖10為1 024×1 024×1 024與256×256×256的輸入時各FFT算法在平臺A上的加速情況對比。據(jù)圖可知，各算法在接收大輸入數(shù)據(jù)量時加速比率要高于接收小輸入數(shù)據(jù)量時的加速比率，并且FFTW的加速效果較優(yōu)，MPFFT的加速效果較差。原因可歸為兩點：

Fig.10 Speedup comparison of different algorithms when input size is 1 024×1 024×1 024 and 256×256×256(PlatformA)圖10 數(shù)據(jù)規(guī)模為1 024×1 024×1 024與256×256×256時不同算法加速情況對比（平臺A）

（1）各算法在接收大規(guī)模數(shù)據(jù)輸入時的性能表現(xiàn)能較好地排除接收小規(guī)模數(shù)據(jù)時所產(chǎn)生的偶然性，更能反映出其真正的性能情況，并且PFFT與MPFFT算法針對接收大規(guī)模數(shù)據(jù)時的操作進行了相關(guān)的優(yōu)化。

（2）PFFT與MPFFT算法由于進行了優(yōu)化，使得其在較小進程數(shù)時已產(chǎn)生相對FFTW較高的性能，而隨進程數(shù)升高，上述分析已說明其優(yōu)化效果將受平臺環(huán)境的制約而減弱，故其加速情況也將受影響而小于FFTW的加速比。

5.2 大規(guī)模集群上算法性能分析

本節(jié)探討的是在大規(guī)模集群上的FFT算法測試。針對平臺A是單節(jié)點服務(wù)器不能進行大規(guī)模集群上的運算從而無法獲得算法可擴展性的問題，本文使用了平臺B——天河一號超級計算機（TH-1A）。

TH-1A以持續(xù)性能每秒2 507萬億次列世界超級計算機排名第五位，其配備了14 336顆至強X5670處理器、7 168塊基于英偉達的Tesla M2050計算卡、2 048顆國防科大的飛騰處理器以及5 PB存儲設(shè)備。

本文使用了256×256×256、512×512×512、1 024× 1 024×1 024的輸入數(shù)據(jù)分別在TH-1A上進行測試，且每個節(jié)點進行多個進程的并行計算獲得如下結(jié)果。

由于GPU的內(nèi)存大小限制，MPFFT算法不能在TH-1A上接收輸入數(shù)據(jù)規(guī)模為512×512×512及以上的輸入。圖11、圖12分別是輸入為1 024×1 024× 1 024和512×512×512時FFTW與PFFT算法運行時間對比。由圖可以分析出，在大規(guī)模集群上接收512×512×512及以上規(guī)模的輸入數(shù)據(jù)時，PFFT算法耗時始終要少于FFTW算法，且進程數(shù)量達到256及以上時FFTW算法的性能將達到瓶頸，而PFFT算法至多可擴展到2 048個進程上繼續(xù)進行計算，且性能也在不斷提升。

Fig.11 Running time comparison of different algorithms when input size is 1 024×1 024×1 024(TH-1A)圖11 數(shù)據(jù)規(guī)模為1 024×1 024×1 024時不同算法運行時間對比（TH-1A）

Fig.12 Running time comparison of different algorithms when input size is 512×512×512(TH-1A)圖12 數(shù)據(jù)規(guī)模為512×512×512時不同算法運行時間對比（TH-1A）

圖13展示了FFTW與PFFT算法在接收512× 512×512和1 024×1 024×1 024的輸入數(shù)據(jù)時的性能對比。對比表示，在進程數(shù)量較少時，同一種算法在接收不同規(guī)模的輸入時性能差異不大，并且PFFT相對于FFTW算法的優(yōu)勢也不明顯；然而，PFFT算法的性能在接收1 024×1 024×1 024的輸入反而比接收512×512×512的輸入要高很多，且隨著進程數(shù)量的增加性能優(yōu)勢越明顯。主要原因為PFFT算法的設(shè)計考慮到了多進程上的并行優(yōu)化，所以輸入數(shù)據(jù)量越大，使用的進程數(shù)量越多，利用PFFT算法進行優(yōu)化計算得到的效果就越明顯。這充分說明了PFFT算法在大規(guī)模輸入多進程上的可擴展性非常強。另外，F(xiàn)FTW算法在使用多進程計算時優(yōu)化性能至多支持到256個進程，且在接收億級輸入數(shù)據(jù)時在256個進程上已經(jīng)開始出現(xiàn)明顯的性能下降現(xiàn)象，這也說明了FFTW算法在多進程上的可擴展性不高。

Fig.13 Performance comparison of different algorithms when input size is 1 024×1 024×1 024 and 512×512×512(TH-1A)圖13 數(shù)據(jù)規(guī)模為1 024×1 024×1 024和512×512×512時不同算法運行性能對比（TH-1A）

圖14與圖15分別是在接收256×256×256輸入數(shù)據(jù)量時3種算法的運算時間對比和性能對比。由圖分析可知，在進程數(shù)量較小的情況下，MPFFT算法性能優(yōu)勢明顯；隨著進程數(shù)量的增多，在進程數(shù)小于64時，3種算法的性能相差不大；當FFTW算法運行的進程數(shù)大于32，MPFFT算法運行的進程數(shù)大于128時，F(xiàn)FTW與MPFFT算法的性能都開始呈下降趨勢，且這兩種算法隨著進程數(shù)的翻倍增長性能卻沒有明顯的增長，而此時PFFT算法仍有較好的性能。

Fig.14 Running time comparison of different algorithms when input size is 256×256×256(TH-1A)圖14 數(shù)據(jù)規(guī)模為256×256×256時不同算法運行時間對比（TH-1A）

Fig.15 Performance comparison of different algorithms when input size is 256×256×256(TH-1A)圖15 數(shù)據(jù)規(guī)模為256×256×256時不同算法運行性能對比（TH-1A）

綜合分析，MPFFT算法主要是針對異構(gòu)平臺上的FFT進行優(yōu)化，當數(shù)據(jù)量較小，不斷增大進程數(shù)時，計算時間上的減小不再明顯，相對于較小的數(shù)據(jù)計算時間，增加的節(jié)點間的通信以及多個GPU的啟用時間會使得整體的性能下降；PFFT算法著重多進程上的并行優(yōu)化和改進，當進程數(shù)量增加時其優(yōu)化效果將會逐漸顯露；而FFTW算法并未有上述兩種算法的優(yōu)化，故進程數(shù)漸增，MPFFT算法優(yōu)化性能開始下降而PFFT算法優(yōu)化性能效果還不明顯的時候，其將會有稍優(yōu)于MPFFT與PFFT算法的性能表現(xiàn)。

圖16為在接收256×256×256、512×512×512和1 024×1 024×1 024輸入數(shù)據(jù)量時各FFT算法的加速情況對比。為了使加速情況對比更為清晰，圖中橫軸使用的是以2為底進程數(shù)的對數(shù)，縱軸使用的是以2為底加速比的對數(shù)。在線性加速比的對比下，各加速算法在接收更大規(guī)模輸入數(shù)據(jù)時的加速比是相對更高的，這也印證了之前分析的正確性。未經(jīng)優(yōu)化的FFTW算法隨著進程數(shù)增加其加速比呈現(xiàn)了波動的趨勢，而PFFT和MPFFT算法的波動相對要小得多。其中，以接收輸入數(shù)據(jù)量為1 024×1 024×1 024時PFFT算法加速比最穩(wěn)定且最貼近線性加速比，說明PFFT算法在大規(guī)模計算時的擴展性還是較高的。

Fig.16 Speedup comparison of different algorithms under different input size(TH-1A)圖16 多種輸入數(shù)據(jù)規(guī)模下不同算法加速情況對比（TH-1A）

之所以出現(xiàn)加速比波動的情況，主要原因為隨著進程數(shù)不斷地呈兩倍提升，額外帶來的通信開銷卻將并行效果相應(yīng)地弱化，兩倍的計算資源不但不能提供兩倍的加速效果，反而還會帶來性能的降低。

圖17為3種FFT算法在TH-1A上的性能評估。橫軸使用以2為底進程數(shù)量的對數(shù)值來表示進程的規(guī)模，縱軸使用的是以2為底輸入數(shù)據(jù)第一維度的對數(shù)值來表示輸入的規(guī)模。氣泡的大小代表在指定的進程數(shù)量和指定的輸入規(guī)模下的性能值，氣泡越大，在該環(huán)境下的性能越好，反之則越差。

Fig.17 Performance evaluation of different algorithms on TH-1A圖17 天河1號不同算法性能評估

由圖17可分析得知，在大多數(shù)的情況下PFFT算法與MPFFT算法要優(yōu)于FFTW算法，且當輸入數(shù)據(jù)量較小時，MPFFT算法性能優(yōu)勢相對明顯；輸入數(shù)據(jù)量較大時，PFFT算法性能優(yōu)勢則相對明顯。另外，這種性能優(yōu)勢會隨著進程數(shù)量的增加而更為突出，當進程數(shù)量較小時，3種算法的性能優(yōu)劣對比并不明顯。

5.3 大規(guī)模集群上算法性能分析

本節(jié)探討的是單一平臺（CPU）下的大規(guī)模集群上的FFT算法測試。在天河一號上的運行環(huán)境是異構(gòu)平臺多進程，而在本節(jié)中將使用天河二號超級計算機（TH-2）來重點評估多節(jié)點上的PFFT算法運行情況。

天河二號是由國防科技大學研制的異構(gòu)超級計算機，共有16 000個運算節(jié)點，每個節(jié)點配備兩顆Xeon E5 12核心的中央處理器，3個Xeon Phi 57核心的協(xié)處理器。本文采用了PFFT算法進行測試以分析其可擴展性。實驗中，為了更準確地對算法在大規(guī)模集群上的性能進行測量，在每個節(jié)點上均指定運行一個進程，輸入數(shù)據(jù)量為1 024×1 024×1 024。

數(shù)據(jù)結(jié)果如圖18和圖19所示，每個節(jié)點上運行一個進程，申請的最大節(jié)點數(shù)為256，因此該數(shù)據(jù)可以真實有效地反映在大規(guī)模集群上FFT算法的性能情況。從圖中可以看出，隨著節(jié)點數(shù)的增加，PFFT算法的性能也在逐步地提升，且加速比大致符合線性加速比，實際的加速情況與理想加速情況的最大差距不超過16%（16節(jié)點時）。因此，該數(shù)據(jù)可證明在大規(guī)模集群上（256節(jié)點）PFFT算法有著較為理想的加速實現(xiàn)，即該算法在多節(jié)點上的優(yōu)化已較為理想。

Fig.18 Performance evaluation of PFFT algorithm when input size is 1 024×1 024×1 024（TH-2）圖18 數(shù)據(jù)規(guī)模為1 024×1 024×1 024時PFFT算法運行性能（TH-2）

Fig.19 Speedup comparison of PFFT algorithm when input size is 1 024×1 024×1 024（TH-2）圖19 數(shù)據(jù)規(guī)模為1 024×1 024×1 024時PFFT算法加速情況對比（TH-2）

6 總結(jié)

多維FFT最終都是轉(zhuǎn)換成一維FFT來計算，目前分布式FFT算法主要有一維切片分解和二維束分解[14]。前者可擴展性比較差，目前只能應(yīng)用于小規(guī)模的集群系統(tǒng)上。而后者需要在計算中穿插更多的全局轉(zhuǎn)置和本地轉(zhuǎn)置，計算過程相對比較復雜。

本文針對CPU和GPU的異構(gòu)集群系統(tǒng)，研究了多維FFT并行算法進一步的優(yōu)化方法，從通信和輸出等方面提出的優(yōu)化策略可為現(xiàn)行FFT算法的改進提供參考。在此基礎(chǔ)上在不同的平臺上對3種FFT算法進行了全方位的性能評估。性能評估的結(jié)果也可以為FFT算法在大規(guī)模運算時提供可靠和全面的參考。

[1]SKA Telescope.The square kilometre array SKA home[EB/ OL].[2016-09-28].https://www.skatelescope.org/.

[2]Frigo M,Johnson S G.FFTW:an adaptive software architecture for the FFT[C]//Proceedings of the 1998 International Conference on Acoustics,Speech and Signal Processing, Seattle,USA,May 12-15,1998.Piscataway,USA:IEEE, 1998:1381-1384.

[3]Pippig M.PFFT:an extension of FFTW to massively parallel architectures[J].SIAM Journal on Scientific Computing, 2013,35(3):C213-C236.

[4]Li Yan,Zhang Yunquan,Liu Yiqun,et al.MPFFT:an autotuning FFT library for OpenCL GPUs[J].Journal of Computer Science and Technology,2013,28(1):90-105.

[5]Nvidia C.CUFFT library[DB].2010.

[6]Chen Yifeng,Cui Xiang,Mei Hong.Large-scale FFT on GPU clusters[C]//Proceedings of the 24th International Conference on Supercomputing,Tsukuba,Japan,Jun 2-4,2010.New York:ACM,2010:315-324.

[7]Cui Xiang,Chen Yifeng,Mei Hong.Improving performance of matrix multiplication and FFT on GPU[C]//Proceedings of the 15th International Conference on Parallel and Distributed Systems,Shenzhen,China,Dec 8-11,2009.Washington:IEEE Computer Society,2009:42-48.

[8]Pippig M.PFFT user manual.2014.

[9]Frigo M,Johnson S G.The design and implementation of FFTW3[J].Proceedings of the IEEE,2005,93(2):216-231.

[10]Bracewell R N.The Fourier transform and its application[M]. 3rd ed.New York:McGraw-Hill,1999.

[11]Cooley J W,Tukey J W.An algorithm for the machine calculation of complex Fourier series[J].Mathematics of Computation,1965,19(90):297-301.

[12]Cooley J W,Lewis P A W,Welch P D.Historical notes on the fast Fourier transform[J].IEEE Transactions on Audio &Electroacoustics,1967,15(2):76-79.

[13]Gholami A,Hill J,Malhotra D,et al.AccFFT:a library for distributed-memory FFT on CPU and GPU architectures[J]. arXiv:1506.07933v2,2015.

[14]Li Yan,Zhang Yunquan,Jia Haipeng,et al.Automatic FFT performance tuning on OpenCL GPUs[C]//Proceedings of the 17th International Conference on Parallel and Distributed Systems,Tainan,China,Dec 7-9,2011.Washington:IEEE Computer Society,2011:228-235.

[15]Liu Yiqun,Li Yan,Zhang Yunquan,et al.Memory efficient two-pass 3D FFT algorithm for Intel?Xeon PhiTM,coprocessor[J].Journal of Computer Science and Technology, 2014,29(6):989-1002.

LI Kun was born in 1992.He is a Ph.D.candidate at University of Chinese Academy of Sciences.His research interests include high performance computing and parallel software,etc.

李琨（1992—），男，山東青島人，中國科學院計算技術(shù)研究所博士研究生，主要研究領(lǐng)域為高性能計算，并行軟件等。

賈海鵬（1983—），男，山東濰坊人，2013年于中國海洋大學獲得博士學位，現(xiàn)為中國科學院計算技術(shù)研究所助理研究員，主要研究領(lǐng)域為異構(gòu)計算，多核并行編程方法，多核上的計算機視覺算法等。

CAO Ting was born in 1984.She received the Ph.D.degree in computer system from The Australian National University in 2014.Now she is a post-doctor at Institute of Computing Technology,Chinese Academy of Sciences.Her research interests include high-level language implementation,novel computer architecture design,hybrid memory management and high performance computing,etc.

曹婷（1984—），女，遼寧撫順人，2014年于澳大利亞國立大學獲得博士學位，現(xiàn)為中國科學院計算技術(shù)研究所博士后，主要研究領(lǐng)域為高級語言實現(xiàn)，新型計算機體系結(jié)構(gòu)設(shè)計，異質(zhì)內(nèi)存管理，高性能計算等。

ZHANG Yunquan was born in 1973.He received the Ph.D.degree in parallel software from Institute of Software, Chinese Academy of Sciences in 2000.Now he is a professor and Ph.D.supervisor at Institute of Computing Technology,Chinese Academy of Sciences.His research interests include high performance parallel computing,with particular emphasis on large scale parallel computation and programming models,high-performance parallel numerical algorithms,performance modeling and evaluation for parallel programs,etc.

張云泉（1973—），男，山東聊城人，2000年于中國科學院軟件研究所獲得博士學位，現(xiàn)為中國科學院計算技術(shù)研究所研究員、博士生導師，主要研究領(lǐng)域為高性能并行計算，尤其是大規(guī)模并行計算及編程模型，高性能并行數(shù)值算法，并行程序的性能建模和評估等。

Implementation and Optimization of Multidimensional FFT Algorithm on Large-Scale Clusters*

LI Kun1,2+,JIAHaipeng1,CAO Ting1,ZHANG Yunquan1
1.State Key Laboratory of Computer Architecture,Institute of Computing Technology,Chinese Academy of Sciences, Beijing 100190,China
2.School of Computer and Control Engineering,University of ChineseAcademy of Sciences,Beijing 100190,China +Corresponding author:E-mail:likun@ict.ac.cn

Fast Fourier transform(FFT)is a fast algorithm which computes the discrete Fourier transform(DFT)of a sequence,or its inverse.Fast Fourier transform is widely used for many applications in engineering,science and mathematics,such as signal decomposition,digital filter and image processing.As a result,the fine-grained optimization of FFT is extremely important in application.This paper studies the typical decomposition strategies of FFT, the parallel implementation of FFT algorithms on large-scale clusters and proposes some optimization strategies.Onthat basis,this paper also evaluates a variety of FFT algorithms performance on different platforms,then analyzes the implementation,strength and weakness,and the computing scalability on large-scale clusters of these algorithms. The experimental results show that the research of this paper can contribute to the further optimization on existing FFT algorithms,and instruct to choose an FFT algorithm which performance is better on large-scale heterogeneous systems(CPU and GPU)in order to satisfy the specified requirements.

cluster;fast Fourier transforms(FFT);message passing interface(MPI);performance optimization

ng was born in 1983.He

the Ph.D.degree in parallel software from Ocean University of China in 2013.Now he is an assistant professor at Institute of Computing Technology,Chinese Academy of Sciences.His research interests include heterogeneous computing,many-core parallel programming method and computer vision algorithms on multi-/many-core processors,etc.

A

TP302.7

*The National Natural Science Foundation of China under Grant Nos.61432018,61133005,61272136,61521092,61502450(國家自然科學基金);the National Key Research and Development Program of China under Grant No.2016YFB0200803(國家重點研發(fā)計劃);the Postdoctoral Science Foundation of China under Grant No.2015T80139(中國博士后科學基金);the National High Technology Research and Development Program of China under Grant Nos.2015AA01A303,2015AA011505(國家高技術(shù)研究發(fā)展計劃(863計劃));the Key Technology Research and Development Programs of Guangdong Province under Grant No.2015B010108006 (廣東省科技計劃項目);the CAS Interdisciplinary Innovation Team of Efficient Space Weather Forecast Models(中科院高效空間天氣預報模式科技創(chuàng)新交叉與合作團隊項目).

Received 2016-11,Accepted 2017-01.

CNKI網(wǎng)絡(luò)優(yōu)先出版:2017-01-05,http://www.cnki.net/kcms/detail/11.5602.TP.20170105.0828.006.html