ARL中Gridding算法的并行化實現(xiàn)

吳懷廣 劉琳琳 石永生 李代祎 謝鵬杰

中圖分類號：F222.3文獻(xiàn)標(biāo)識碼：ADOI：10.3969/j.issn.2096-1553.2019.02.011

文章編號：2096-1553（2019）02-0082-06

關(guān)鍵詞：ARL;并行化算法;Gridding算法;CUDA

Key words：ARL;parallelization algorithm;Gridding algorithm;CUDA

摘要：針對海量天文數(shù)據(jù)實時性處理效率低的問題，通過對SKA圖像采集及成像ARL算法庫中耗時較長的Gridding算法進(jìn)行耗時分析，找出了該算法中調(diào)用頻率高且運行時間長的兩個函數(shù)convolutional-grid和convolutional-degrid，利用GPU的多線程并行化處理降低兩個函數(shù)的循環(huán)迭代，實現(xiàn)了Gridding算法在GPU和CPU上的協(xié)同運行.驗證實驗結(jié)果表明，在相同的數(shù)據(jù)量下，改進(jìn)后的Gridding算法運行時間大大縮短，特別是在處理海量數(shù)據(jù)時，有效提高了ARL的整體運行效率.

Abstract：Aiming at the low real-time processing efficiency of massive astronomical data， through time-consuming analysis of gridding algorithm in SKA image acquisition and imaging ARL library， two functions of convolutional-grid and convolutional-degrid with high frequency and long running time were found out in this algorithm. Then， two functions were parallelized on GPU by multi-threading to realize the cooperative operation of gridding algorithm on GPU and CPU. The experimental results showed that under the same amount of data， the running time of the improved gridding algorithm was greatly shortened， especially when dealing with massive data， the overall running efficiency of ARL was effectively improved.

0 引言

在過去的幾十年里，射電望遠(yuǎn)鏡的靈敏度和圖像分辨率均有很大的提升.SKA[1-2]作為世界上最大綜合孔徑的射電望遠(yuǎn)鏡，采集數(shù)據(jù)的速率非?？?，數(shù)據(jù)采集量非常大，其設(shè)計目標(biāo)是要大于12 TB/s.海量數(shù)據(jù)的產(chǎn)生和天文成像本身對實時處理的嚴(yán)格要求，給計算機(jī)的計算能力帶來了巨大的挑戰(zhàn).海量數(shù)據(jù)的科學(xué)處理一般需要百億億次量級的超級計算機(jī)來完成.目前，通過提高主頻來提高CPU處理能力的傳統(tǒng)方式已受到集成電路集成度的制約，因此利用多核CPU和多核加速器（如GPU，Cell/BE等）處理大量時效性數(shù)據(jù)成為發(fā)展趨勢.在天文成像過程中，網(wǎng)格化和去網(wǎng)格化是最耗時的兩個操作.如果處理的可見度數(shù)據(jù)是EB量級或以上，則消耗的時間不能通過調(diào)整計算機(jī)性能來緩和.傳統(tǒng)的網(wǎng)格化一般都是在CPU上運行，加速器在每一個浮點上的帶寬很小，應(yīng)用的時效性也相應(yīng)較低.所以，并行化算法作為提高計算速度的一個重要途徑，在射電天文學(xué)數(shù)據(jù)處理[3]方面越來越受到關(guān)注.

網(wǎng)格化（Gridding）算法的作用是將落在笛卡爾坐標(biāo)外的數(shù)據(jù)點插值到網(wǎng)格上.W.N.Brouw[4]最早提出了Gridding方法，用于實現(xiàn)極坐標(biāo)網(wǎng)格采樣的離散傅里葉變換.針對極坐標(biāo)網(wǎng)格采樣，J.D.OSullivan[5]提出了一種快速的sinc函數(shù)網(wǎng)格算法，并將其應(yīng)用在CT圖像重建上，該算法利用有限范圍的卷積函數(shù)得到sinc函數(shù)插值以減小計算時間，但是計算量太大.C.H.Meyer等[6]將Gridding算法應(yīng)用于磁共振成像（MRI）中的螺旋采樣，采用的卷積核是Kaiser-Bessel窗口函數(shù)，但是Kaiser-Bessel窗口函數(shù)不是最優(yōu)的卷積函數(shù).勞保強等[7]使用卷積函數(shù)網(wǎng)格化不規(guī)則的微波全息數(shù)據(jù)，分析了不同卷積函數(shù)的抗混疊性能，驗證了球體函數(shù)是最具抗混淆性能的卷積函數(shù).A.L.Varbanescu等[8]以增加附加計算為代價，通過對數(shù)據(jù)的排序和搜索來改善空間局部性，進(jìn)而改善內(nèi)存的性能，但其結(jié)果尚未達(dá)到GPU運算峰值的14%.

以上研究主要是在CPU環(huán)境下對Gridding算法的改進(jìn)，鮮見在ARL（SKA的算法參考庫，用以實現(xiàn)整體天文數(shù)據(jù)成像流程）中針對Gridding算法在GPU實現(xiàn)時的運行效率的研究.鑒于此，本文擬分析CPU環(huán)境下ARL中耗時最長的Gridding算法，利用性能分析工具找出Gridding算法中耗時最長的函數(shù)模塊，然后對函數(shù)模塊使用CUDA進(jìn)行GPU并行執(zhí)行，以期利用GPU在并行運算方面的優(yōu)勢，實現(xiàn)Gridding算法的CPU和GPU協(xié)同處理，從而提高ARL的整體運行效率.

1 Gridding算法在ARL網(wǎng)格化中的應(yīng)用

在射電天文學(xué)中，所采集到的信號通常受到各類因素影響，數(shù)據(jù)呈現(xiàn)出不規(guī)律性.在利用射電望遠(yuǎn)鏡采集到的數(shù)據(jù)重建天空圖像時，需要將不規(guī)則采樣的數(shù)據(jù)映射到標(biāo)準(zhǔn)的二維網(wǎng)格中，這稱為網(wǎng)格化.只有經(jīng)歷這一步，網(wǎng)格才可以通過快速傅里葉變換來重建天空圖像.在射電望遠(yuǎn)鏡數(shù)據(jù)處理中，網(wǎng)格化是計算量最大、耗費時間最長的步驟.

在天文成像中，對于天空亮度I，初始光束模式A和能見度值V，存在以下傅里葉變換關(guān)系：

其中，S表示（u，v）采樣函數(shù)，V′表示觀測值.

對式②作數(shù)值估計的方法有直接傅里葉變換DFT（direct Fourier transform）和快速傅里葉變換FFT（fast Fourier transform）兩種方法.其中，DFT對N×N網(wǎng)格的每個點都做計算，通過蠻力計算求和來估計ID（l，m），其表達(dá)式為

ID（l，m）=1M∑Mk=1V′（uk，vk）e2πi（ukl+vkm）

FFT則通過利用DFT運算中的對稱性和周期性，減少DFT的運算量，從而更快地估計ID（l，m）.但是利用FFT需要先將數(shù)據(jù)插值到矩形網(wǎng)格點上，然后才可以應(yīng)用FFT.插值的過程就稱為網(wǎng)格化.ARL中Gridding算法的作用是將落在笛卡爾坐標(biāo)外的數(shù)據(jù)點插值到網(wǎng)格上，以便對網(wǎng)格化的數(shù)據(jù)進(jìn)行FFT計算.

天文學(xué)中的Gridding算法步驟如下.

1）密度補償函數(shù)與采樣數(shù)據(jù)相乘.Gridding算法首先要用密度補償函數(shù)與采樣數(shù)據(jù)相乘來彌補采樣數(shù)據(jù)的不均勻.網(wǎng)格算法中的采樣密度補償函數(shù)和插值函數(shù)（卷積函數(shù)）的選擇對圖像重建呈現(xiàn)的最后效果影響較大，也是網(wǎng)格化算法研究的核心問題.

2）在網(wǎng)格點（uc，vc）處計算出卷積，即

3）對數(shù)據(jù)進(jìn)行重采樣，使數(shù)據(jù)落到網(wǎng)格點上.在所有網(wǎng)格點上對C×VW作采樣處理，即

4）應(yīng)用FFT對VR作傅里葉變換，得出經(jīng)過簡單修正的初步圖像.VR是規(guī)則間隔的δ函數(shù)的線性組合，可以由離散傅里葉變換計算出其傅里葉變換VR的采樣矩陣，

5）為了消除卷積函數(shù)的影響，需要對重建后的圖像作進(jìn)一步的修正，即除以卷積函數(shù)的傅里葉變換.

2 Gridding算法的耗時分析

Gridding中一共含有14個函數(shù)，對這些函數(shù)進(jìn)行分析，選出調(diào)用頻率較高、運行時間最長的4個函數(shù)：1）grdsf函數(shù)，其功能是獲取網(wǎng)格函數(shù)及對網(wǎng)格進(jìn)行修正;2）anti_aliasing_calculate 函數(shù)，其功能是計算prolate球體反混疊函數(shù);3）convolutional_degrid 函數(shù)，其功能是利用被采樣的gcf部分uv坐標(biāo)值，進(jìn)行卷積解網(wǎng)格;4）convolutional_grid函數(shù)，其功能是由頻率和偏振獨立的gcf進(jìn)行卷積網(wǎng)格.Gridding算法中其他函數(shù)因其調(diào)用次數(shù)和單個函數(shù)耗時都很小，在整個過程中可以忽略不計.

選出的4個函數(shù)被調(diào)用次數(shù)如圖1所示，執(zhí)行時間如圖2所示.分析中可見性數(shù)據(jù)量為104條.

由圖1可以看出，grdsf函數(shù)的被調(diào)用次數(shù)最高，達(dá)到370次，anti_aliasing_calculate函數(shù)被調(diào)用185次，convolutional_degrid函數(shù)被調(diào)用47次，convolutional_grid函數(shù)被調(diào)用44次.

由圖2可以看出，convolutional_grid和convolutional_degrid兩個函數(shù)運行時間最長，而grdsf和anti_aliasing_calculate兩個函數(shù)運行時間較少.

綜合來看，需要對convolutional_grid和convolutional_degrid兩個函數(shù)進(jìn)行并行化處理.

3 Gridding算法的并行化實現(xiàn)

在整個SKA圖像采集和最終成像的過程中，Gridding算法是計算量大、耗時長的算法之一.ARL算法庫是基于CPU運行的，利用PyCUDA 進(jìn)行并行化加速可以使整個ARL的運行時間得到有效縮短.

通常，GPU并行加速方式兩種：一種是在串行代碼的基礎(chǔ)上通過降低循環(huán)進(jìn)行加速;另一種是重新設(shè)計算法，將總模塊劃分為多個獨立的子模塊，通過分配線程來并行處理子任務(wù).本文采用第一種優(yōu)化措施.

這是因為Gridding算法在處理數(shù)據(jù)的時候迭代次數(shù)很多，通過對Gridding算法的調(diào)用次數(shù)和耗時的分析可以看出，計算量最多的是convolutional_grid和convolutional_degrid函數(shù)，而這兩個函數(shù)運算時for循環(huán)較多，且數(shù)據(jù)之間不存在依賴關(guān)系，所以可以直接利用多線程實現(xiàn)其在GPU下的并行計算，從而達(dá)到降低for循環(huán)次數(shù)，加快運算速度的目的.

由于GPU與CPU之間不能直接通信，所以進(jìn)行并行化處理之前需將數(shù)據(jù)從主機(jī)端（CPU）拷貝到設(shè)備端（GPU），然后才能執(zhí)行相關(guān)的核函數(shù).同樣的，對于GPU并行化處理后的結(jié)果也需要從GPU拷貝到CPU.

本文利用全局存儲器進(jìn)行并行化處理.實驗所用顯卡的靜態(tài)存儲器的大小是11 440 MB，因此所有的線程模塊中的線程不會產(chǎn)生訪問沖突.利用CUDA[9]中的cudaMemcpy（）將數(shù)據(jù)載入到GPU的靜態(tài)存儲器中，以便所有的線程共享.根據(jù)線程索引獲取可見性數(shù)據(jù)進(jìn)行求和.等到線程都完成計算，最后將結(jié)果載入到內(nèi)存.進(jìn)行GPU并行處理的流程圖如圖3所示.

4 驗證實驗結(jié)果與分析

實驗環(huán)境：CPU 為Intel Xeon E5-2620 V3，內(nèi)存為816GB DDR4，內(nèi)存最大支持為768 GB，操作系統(tǒng)是CentOS 7.0，python版本為3.5;GPU采用Tesla K80，為雙核GK210架構(gòu)，主板型號為MG50-G20.

實驗中選用的GPU一共含有4992個CUDA Cores，每個塊最多分配1024個線程，對于不同的數(shù)據(jù)量可以選擇不同的線程和線程塊.本實驗使用的數(shù)據(jù)量分別為103條，104條，105條，106條，107條，108條，針對不同的數(shù)據(jù)量，線程塊統(tǒng)一采用（32，32，1）的形式，采用的網(wǎng)格塊的分配分別是（1，1，1），（10，1，1），（10，10，1），（100，10，1），（100，100，1），（1000，100，1）.

ARL中convolutional_grid函數(shù)和convolutional_degrid函數(shù)在CPU上運行的時間與在GPU加速后的運行時間見表1，其中加速比指在CPU上運行時間和GPU上運行時間的比值.

從表1可以看出，隨著數(shù)據(jù)量的增加，兩個函數(shù)的加速比不斷增加，但是加速比增加的速度逐漸減小.為了更直觀地表示加速比和數(shù)據(jù)量之間的關(guān)系，采用直方圖進(jìn)行表示，結(jié)果如圖4和圖5所示.

由圖4和圖5可以看出，當(dāng)數(shù)據(jù)量較少（103條）時，兩個函數(shù)的加速比分別是0.22和0.38，利用GPU并行處理的運行時間反而增加，這是因為數(shù)據(jù)量相對較少時，GPU器件的啟動時間、內(nèi)存和顯存之間數(shù)據(jù)交互等都需要時間.因此在數(shù)據(jù)量比較小時，GPU因其內(nèi)部調(diào)度相對來說耗時較長，其并行加速并不占優(yōu)勢.隨著數(shù)據(jù)量的增加，GPU加速的優(yōu)勢逐漸

顯現(xiàn).當(dāng)數(shù)據(jù)量大于104條時，加速比不斷增加.因此，GPU并行加速適用于大量數(shù)據(jù)的運算情況，可在一定程度上提高ARL的整體運行效率.

4 結(jié)語

本文針對SKA海量天文數(shù)據(jù)實時性處理效率較低的問題，對ARL算法庫中Gridding算法進(jìn)行了耗時分析與優(yōu)化：通過分析Gridding算法中每個函數(shù)的調(diào)用次數(shù)和運行時間，找出了需要優(yōu)化的兩個函數(shù)，即convolutional_grid和convolutional_degrid函數(shù);然后利用GPU的多線程并行化處理降低兩個函數(shù)的循環(huán)迭代，使得ARL中Gridding算法在GPU和CPU上實現(xiàn)了協(xié)同運行.驗證實驗結(jié)果顯示，在數(shù)據(jù)量不影響輸入輸出數(shù)據(jù)格式和大小的情況下，對Gridding算法的并行化處理縮短了ARL的整體運算時間，提高了運行的效率，且數(shù)據(jù)量越大，加速優(yōu)勢越明顯.

本文所做的并行化研究對于天文圖像的成像過程具有一定的參考意義，對Gridding算法在多GPU下實現(xiàn)算法加速將是下一步的研究方向.

參考文獻(xiàn)：

[1] RAZAVI-GHODS N，ACEDO E D L，EL-MAKA-DEMA A，et al.Analysis of sky contributions to system temperature for low frequency SKA aperture array geometries[J].Experimental Astronomy，2012，33（1）：141.

[2] ZHANG Y，BROWN A K.Bunny ear combline antennas for compact wide-band dual-polarized aperture array[J].IEEE Transactions on Antennas and Propagation，2011，59（8）：3071.

[3] 彭曉明，郭浩然，龐建民.多核處理器——技術(shù)、趨勢和挑戰(zhàn)[J].計算機(jī)科學(xué)，2012，39（S3）：320.

[4] BROUW W N.Aperture synthesis[J].Methods in Computational Physics，1975，14：131.

[5] OSULLIVAN J D.A fast sinc function gridding algorithm for fourier inversion in computer tomography[J].IEEE Transactions on Medical Imaging，1985，4（4）：200.

[6] MEYER C H，HU B S，NISHIMURA D G，et al.Fast spiral coronary artery imaging[J].Magnetic Resonance in Medicine，1992，28（2）：202.

[7] 勞保強，王俊義，王錦清，等.基于卷積核網(wǎng)格化二維近程微波全息[J].微波學(xué)報，2014，30（5）：82.

[8] VARBANESCU A L.On the effective parallel programming of multi-core processors[D].Romania：Universitatea Politehnica Bucuresti，2010.

[9] CHENG J.CUDA by example：an introduction to general-purpose GPU programming[M].Boston：Addison-Wesley Professional，2010.