uv-faceting成像并行算法研究?

勞保強(qiáng) 安 濤 于 昂 郭紹光

(1 中國(guó)科學(xué)院上海天文臺(tái) 上海 200030)

(2 中國(guó)科學(xué)院射電天文重點(diǎn)實(shí)驗(yàn)室 南京 210008)

(3 桂林電子科技大學(xué)信息與通信學(xué)院 桂林 541004)

1 引言

為了獲得更高分辨率和靈敏度,綜合孔徑干涉陣列設(shè)備得到了飛速的發(fā)展.國(guó)際上,從上世紀(jì)70年代至90年代末,主要建成了厘米波、毫米波和亞毫米波的綜合孔徑干涉陣列.如美國(guó)的甚大陣(Very Large Array,VLA)[1]、澳大利亞望遠(yuǎn)鏡致密陣(Australia Telescope Compact Array,ATCA)、英國(guó)多天線微波相連干涉儀網(wǎng)(Multi-Element Radio-Linked Interferometer Network,MERLIN)[2]、印度巨米波射電望遠(yuǎn)鏡(Giant Metrewave Radio Telescope,GMRT)[3]等.進(jìn)入21世紀(jì),分米波、米波甚至更長(zhǎng)波段的低頻綜合孔徑干涉陣列重新成為了熱點(diǎn).如荷蘭的低頻陣列(Low Frequency Array,LOFAR)[4]、澳大利亞的默奇森寬視場(chǎng)陣列(Murchison Widefield Array,MWA)[5]、美國(guó)的長(zhǎng)波長(zhǎng)陣列(Long Wavelength Array,LWA)以及即將建造的平方公里陣列(Square Kilometre Array,SKA)[6]等.在國(guó)內(nèi),射電天文研究起步相對(duì)較晚,綜合孔徑干涉陣列屈指可數(shù).目前的設(shè)備有新疆的21 cm陣(21CMA)、內(nèi)蒙古的日像儀和新疆的天籟陣.從國(guó)際形式來(lái)看,射電天文的研究將進(jìn)入以SKA為核心的新時(shí)代.

隨著綜合孔徑干涉陣列設(shè)備的發(fā)展,原始數(shù)據(jù)的規(guī)模隨之增大,同時(shí)數(shù)據(jù)處理也面臨著巨大的挑戰(zhàn).SKA第1階段(SKA1)的科學(xué)數(shù)據(jù)預(yù)處理所需要的計(jì)算能力為260 PFLOPS (即每秒260千萬(wàn)億次浮點(diǎn)運(yùn)算),相當(dāng)于我國(guó)超級(jí)計(jì)算機(jī)天河二號(hào)(33 PFLOPS)的8倍,比目前世界排名第2的中國(guó)神威太湖之光的計(jì)算能力(93 PFLPOS)還高2倍,相當(dāng)于約1億臺(tái)個(gè)人計(jì)算機(jī)的處理能力[7].近年,由于高性能計(jì)算和人工智能等技術(shù)的不斷提高以及天文大數(shù)據(jù)處理的需求,使得利用這些技術(shù)來(lái)優(yōu)化綜合孔徑干涉陣列數(shù)據(jù)處理方法成為了熱門.

大視場(chǎng)成像是SKA等低頻綜合孔徑干涉陣列數(shù)據(jù)處理中較為關(guān)鍵的技術(shù)之一.目前,比較有影響力并且將用于SKA的算法有:uv-faceting[8]、w-projection[9]、wstacking[10]、w-snapshots[11].其中,Romein[12]對(duì)w-projection算法進(jìn)行了加速優(yōu)化,研發(fā)了基于Compute Unified Device Architecture (CUDA)和Open Computing Language(OpenCL)的w-projection成像算法.該方法支持單個(gè)或多個(gè)GPU節(jié)點(diǎn)運(yùn)行,并完成了相關(guān)性能測(cè)試,結(jié)果表明提出的加速方法比以往的方法在性能上提高了一個(gè)數(shù)量級(jí).w-snapshots、w-projection和w-stacking成像算法均是一次性對(duì)整個(gè)視場(chǎng)的圖像進(jìn)行校正處理,這導(dǎo)致了優(yōu)化后的算法只能有效地重建較低分辨率天空?qǐng)D像[13].相反,uv-faceting將視場(chǎng)分割成多個(gè)分面,每個(gè)分面采用傳統(tǒng)的2維傅里葉變換方法重建天空?qǐng)D像,最后將所有分面的圖像進(jìn)行合并.在重建大像素天空?qǐng)D像上,有較好的拓展性,較適合用于SKA這種需要產(chǎn)生萬(wàn)兆像素圖像項(xiàng)目中.此外,與其他方法相比,uv-faceting方法的另一個(gè)優(yōu)點(diǎn)是可以在柵格化和成像之前直接對(duì)數(shù)據(jù)進(jìn)行方向依賴效應(yīng)的校正[14].因此,本文主要對(duì)uv-faceting算法進(jìn)行并行優(yōu)化改造.

首先,介紹了uv-faceting成像算法具體的工程實(shí)現(xiàn)方法,從而分析了該算法各個(gè)處理步驟的時(shí)間消耗情況,并確立對(duì)柵格化(gridding)和數(shù)據(jù)讀取這兩個(gè)部分進(jìn)行并行優(yōu)化.然后,提出了基于Message Passing Interface (MPI)+Open Multi-Processing(OpenMP)的uv-faceting成像算法和基于MPI+CUDA的uv-faceting成像算法.最后,對(duì)這兩種算法進(jìn)行正確性的驗(yàn)證與分析和性能測(cè)試與分析.

2 uv-faceting成像

本節(jié)主要從工程實(shí)現(xiàn)的角度描述uv-faceting成像算法,我們之前的研究給出了詳細(xì)的uv-faceting算法理論推導(dǎo)過(guò)程,最終給出了第j個(gè)分面的亮度分布函數(shù)I與其可見(jiàn)度函數(shù)Vfacet的關(guān)系式[15?16]:

上式中,第j個(gè)分面的可見(jiàn)度函數(shù)可以通過(guò)改變?cè)伎梢?jiàn)度函數(shù)V的相位得到:

其中,lj、mj和nj是第j個(gè)分面中心的方向余弦,e表示指數(shù)函數(shù),?lj=l?lj和?mj=m?mj,l和m是原始觀測(cè)方向的方向余弦,u′和v′由原始的u、v和w坐標(biāo)計(jì)算得到:

這里,圖像坐標(biāo)l和m在工程實(shí)現(xiàn)中可以利用其與相位中心的關(guān)系得到[17],因此lj、mj和nj可以通過(guò)下面的公式計(jì)算得到:

其中,nra和ndec是分面相位中心的赤經(jīng)和赤緯,odec和ora是原始相位中心的赤經(jīng)和赤緯.

從上述的分析,我們很容易得出理論上只要給出需要求解的分面的相位中心,就可以利用傳統(tǒng)的2維傅里葉變換方法求出該分面的亮度分布函數(shù).因?yàn)樵嫉目梢?jiàn)度數(shù)據(jù)、相位中心和u、v和w坐標(biāo)可以在觀測(cè)數(shù)據(jù)中獲取.

圖1是實(shí)際觀測(cè)數(shù)據(jù)進(jìn)行uv-faceting成像的流程圖.這里我們使用的輸入數(shù)據(jù)必須為校準(zhǔn)后的數(shù)據(jù)(calibrated data),格式是Measurement Set (即ms文件)1http://www.astron.nl/casacore/trunk/casacore/doc/notes/229.html,它是目前比較通用的數(shù)據(jù)存儲(chǔ)格式.該文件格式目前版本號(hào)是V2.0,SKA科學(xué)數(shù)據(jù)處理器團(tuán)隊(duì)正在將其升級(jí)為V3.0,并將作為SKA1數(shù)據(jù)的文件存儲(chǔ)格式.圖1中uv-faceting成像主要由3層嵌套循環(huán)來(lái)完成,最外層循環(huán)是按照分面的個(gè)數(shù)執(zhí)行,使用facet-id變量作分面的循環(huán)計(jì)數(shù),num-facets是分面的總數(shù).次外層循環(huán)按照可見(jiàn)度數(shù)據(jù)的行數(shù)來(lái)執(zhí)行,這里使用變量row作行數(shù)的循環(huán)計(jì)數(shù),可見(jiàn)度數(shù)據(jù)的總行數(shù)為row-count.內(nèi)層循環(huán)按照頻率通道數(shù)來(lái)執(zhí)行,變量c作為頻率通道數(shù)的循環(huán)計(jì)數(shù),總的頻率通道數(shù)為channel-count.1次循環(huán)即完成1個(gè)分面1行可見(jiàn)數(shù)據(jù)中1個(gè)頻率通道的柵格化(gridding)處理任務(wù).gridding的實(shí)質(zhì)是可見(jiàn)度函數(shù)與卷積核函數(shù)(卷積核濾波器)進(jìn)行卷積,這里的gridding算法我們將采用casacore grider的設(shè)計(jì)思路實(shí)現(xiàn).循環(huán)結(jié)束后,對(duì)每個(gè)分面得到的gridding結(jié)果進(jìn)行2維快速傅里葉逆變換(IFFT),得到各個(gè)分面的圖像結(jié)果.最后,馬賽克(mosaic)處理所有分面圖像得到最終的目標(biāo)圖像.其中,步驟phase rotation利用(2)式計(jì)算完成;步驟baseline rotation由(3)式和(4)式聯(lián)立計(jì)算得到;步驟scaleuvcoverage根據(jù)傅里葉變換的展縮(尺度變換)特性獲得.例如當(dāng)想獲得的圖像像素大小為nx×ny,并且像元大小(cellsize)設(shè)定為xcellsize和ycellsize時(shí),u、v坐標(biāo)需要分別乘以一個(gè)尺度常量,分別為:

圖1 uv-faceting成像流程圖.圖中N表示否,Y表示是Fig.1 The flow chart of uv-faceting imaging.N indicates no,and Y indicates yes

3 uv-faceting成像并行算法

3.1 基于MPI+OpenMP的uv-faceting成像并行算法

首先,我們探討一下uv-faceting成像算法處理過(guò)程的耗時(shí)情況.從傳統(tǒng)的射電成像中得知,在成像的所有處理步驟中g(shù)ridding是耗時(shí)最多且是較為關(guān)鍵的一步.另外,隨著輸入數(shù)據(jù)的增大,數(shù)據(jù)的讀取所需的時(shí)間也會(huì)急劇增大.因?yàn)閙s文件存儲(chǔ)底層的數(shù)據(jù)輸入輸出(IO)使用的是casacore庫(kù)中的Casacore Table Data System (CTDS)[18],目前CTDS只支持?jǐn)?shù)據(jù)的串行讀寫.所以當(dāng)數(shù)據(jù)量非常大的時(shí)候,數(shù)據(jù)讀取消耗的時(shí)間將是巨大的,這將嚴(yán)重影響后續(xù)的科學(xué)數(shù)據(jù)處理與分析任務(wù)的開展.如圖2是uv-faceting成像關(guān)鍵步驟隨著數(shù)據(jù)量增大所耗時(shí)間的變化曲線,關(guān)鍵步驟依次為卷積核濾波器的建立(Conv-filter Create)、數(shù)據(jù)加載(Data Loading)、柵格化(gridding)和傅里葉逆變換(Fourier Inversion).數(shù)據(jù)大小分別是0.36、3.7、36和349,單位GB.圖中也充分證明了數(shù)據(jù)讀取和gridding是uv-faceting成像算法最耗時(shí)間的部分,其中g(shù)ridding耗時(shí)最長(zhǎng),并且這個(gè)增大呈線性性.基于此,本文將重點(diǎn)對(duì)gridding和數(shù)據(jù)讀取這兩個(gè)部分進(jìn)行并行化改造.

圖2 uv-faceting成像關(guān)鍵步驟耗時(shí)曲線.Conv-filter表示卷積核濾波器Fig.2 The consuming time curve of the key step of uv-faceting imaging.Convfilter represents the convolution kernel filter

在gridding的并行改造方面,從圖1流程可以看出,uv-faceting算法實(shí)際上提供了一種粗粒度并行的方法.這種方法只需對(duì)最外層循環(huán)進(jìn)行并行化,無(wú)需考慮任何內(nèi)層循環(huán)的處理過(guò)程和同步鎖問(wèn)題,最為簡(jiǎn)單方便.這個(gè)粗粒度并行可以使用OpenMP編程實(shí)現(xiàn),本文則是在最外層循環(huán)外加“#pragma omp parallel for schedule (static)”語(yǔ)句實(shí)現(xiàn)多線程并行.這相當(dāng)于將分面的處理過(guò)程分為多個(gè)線程并行實(shí)現(xiàn).例如有n個(gè)分面需要處理,并使用t個(gè)CPU核數(shù)或線程處理時(shí),那么將給每個(gè)線程靜態(tài)分配n/t個(gè)循環(huán).當(dāng)使用的CPU核數(shù)或線程數(shù)接近于分面數(shù),該方法所消耗的時(shí)間最短.

在數(shù)據(jù)讀取方面,大致有兩種方法作并行優(yōu)化:底層存儲(chǔ)數(shù)據(jù)IO并行化和應(yīng)用層數(shù)據(jù)讀取并行化.第1種并行化方法是利用并行IO軟件包,如Adaptable IO System(ADIOS)等,開發(fā)基于該并行IO的數(shù)據(jù)管理系統(tǒng).它可以作為新的管理系統(tǒng)與傳統(tǒng)的串行數(shù)據(jù)管理系統(tǒng)并存,同時(shí)將復(fù)雜的并行化代碼封裝起來(lái)并提供一種簡(jiǎn)潔的接口給用戶使用,因此使用時(shí)避開了繁瑣的編程實(shí)現(xiàn)過(guò)程.這種方法能夠從根本上解決海量數(shù)據(jù)讀寫上的瓶頸,但是目前仍處于驗(yàn)證性實(shí)驗(yàn)階段,還沒(méi)得到廣泛應(yīng)用[19].第2種并行優(yōu)化方法是保留存儲(chǔ)底層數(shù)據(jù)IO的讀寫方式,在用戶應(yīng)用層通過(guò)編寫并行程序?qū)?shù)據(jù)進(jìn)行分割同時(shí)讀取.這種方式需要對(duì)數(shù)據(jù)格式和內(nèi)容有充分的了解,且具備較強(qiáng)的編程功底.本節(jié)的并行算法將使用第2種方法.

圖3是本文提出的基于MPI+OpenMP的uv-faceting成像并行算法.首先,根據(jù)數(shù)據(jù)的存儲(chǔ)格式特點(diǎn)將輸入的校準(zhǔn)后的數(shù)據(jù)利用MPI進(jìn)程按行分割成多個(gè)分塊,分割數(shù)等于MPI進(jìn)程數(shù),圖3使用4個(gè)進(jìn)程或4個(gè)節(jié)點(diǎn)為例;其次,使用MPI將分割之后的數(shù)據(jù)分配到不同的CPU節(jié)點(diǎn)上進(jìn)行并行處理,每個(gè)進(jìn)程或CPU節(jié)點(diǎn)將完成gridding等數(shù)據(jù)處理.同時(shí)在每個(gè)進(jìn)程使用OpenMP多線程并行完成多個(gè)分面的gridding處理;然后,將所有進(jìn)程的gridding結(jié)果進(jìn)行求和并發(fā)送至主節(jié)點(diǎn);最后,在主節(jié)點(diǎn)對(duì)最終的gridding結(jié)果進(jìn)行2維傅里葉逆變換得到目標(biāo)多個(gè)分面的圖像并輸出相應(yīng)的fits格式文件,再將多個(gè)分面的圖像進(jìn)行馬賽克處理得到最終的目標(biāo)圖像.本文中馬賽克的處理軟件采用Swarp[20].

圖3 基于MPI+OpenMP的uv-faceting成像并行策略.k表示觀測(cè)時(shí)間戳的開始索引,m表示觀測(cè)時(shí)間編號(hào),pol表示偏振Fig.3 The uv-faceting imaging parallel strategy based on the MPI and OpenMP.k is the start index of observation time stamp,m is the numbering of observation time,and pol is the polarization

3.2 基于MPI+CUDA的uv-faceting成像并行算法

就單節(jié)點(diǎn)而言,CPU的核數(shù)是非常有限的,而GPU核處理器則通常有成千上萬(wàn)個(gè)處理核.因此,上述的gridding處理步驟,使用GPU加速會(huì)更加有效.本文主要采用CUDA編程實(shí)現(xiàn)gridding加速,利用了由相同基線進(jìn)行連續(xù)測(cè)量之間的空間局部性.由于連續(xù)測(cè)量之間的積分時(shí)間足夠短,相同基線的連續(xù)測(cè)量數(shù)據(jù)將會(huì)被gridding到相同的網(wǎng)格點(diǎn)上.我們把相同基線下連續(xù)測(cè)量的多個(gè)時(shí)間戳數(shù)據(jù)的gridding結(jié)果進(jìn)行累加并緩存在寄存器中,當(dāng)處理下一個(gè)基線數(shù)據(jù)時(shí),再將寄存器的數(shù)據(jù)寫入全局存儲(chǔ)器(global memory).當(dāng)global memory的寫訪問(wèn)是原子性(atomic)的,多線程之間的不存在競(jìng)爭(zhēng),多個(gè)基線就可以實(shí)現(xiàn)并行g(shù)ridding.在gridding處理中,我們給每個(gè)卷積核濾波器的抽頭分配一個(gè)線程,一共需要CFull support-size×CFull-support-size個(gè)線程,其中CFull-support-size是卷積核濾波器的全支撐大小.因此,基于該方法的uv-faceting成像算法需要的最小線程數(shù)量為Nfacets×NBaselines×CFull-support size×CFull-support-size,其中Nfacets是分面總數(shù),NBaselines是基線總數(shù).

圖4是基于CUDA的單個(gè)GPU節(jié)點(diǎn)uv-faceting成像算法流程.圖中主機(jī)先進(jìn)行輸入?yún)?shù)的初始化,同時(shí)對(duì)GPU device進(jìn)行初始化,包括device的重置和device的選擇,因?yàn)橛行〨PU節(jié)點(diǎn)會(huì)掛載多個(gè)GPU加速卡.隨后,主機(jī)開始讀數(shù)據(jù),再將讀入的數(shù)據(jù)進(jìn)行重新排序和打包,最后將重排后的數(shù)據(jù)傳送至GPU的global memory.GPU device經(jīng)過(guò)初始化后,將會(huì)事先給從主機(jī)傳送來(lái)的數(shù)據(jù)開辟相應(yīng)的數(shù)組memory.數(shù)據(jù)傳送到GPU global memory后,在GPU device中完成數(shù)據(jù)的gridding處理任務(wù).首先,獲取CUDA線程索引(thread-id),這里線程總數(shù)設(shè)置為上述提到的最小線程數(shù)量,線程塊大小為256;其次,根據(jù)這個(gè)線程索引、設(shè)置的線程總數(shù)和線程塊數(shù)量,計(jì)算出當(dāng)前的基線起始索引、卷積核濾波器支撐索引(u-support和v-support)、每個(gè)基線的時(shí)間戳數(shù)量(num-t)和分面ID號(hào)(facet-id);然后,開始完成當(dāng)前基線和分面下,多個(gè)時(shí)間戳的gridding并累加,累加的結(jié)果緩存到CUDA寄存器中,其中,這個(gè)gridding累加只計(jì)算當(dāng)前卷積核支撐索引下的卷積結(jié)果.當(dāng)需要進(jìn)行下一個(gè)基線數(shù)據(jù)處理前,當(dāng)前基線gridding結(jié)果將從寄存器拷貝到global memory;最后,把所有基線和分面計(jì)算的gridding結(jié)果傳送給主機(jī).主機(jī)完成IFFT和mosaic,得到最終的圖像.

在數(shù)據(jù)讀取方面,我們同樣采用3.1節(jié)中的第2種方法,不同的是分割數(shù)據(jù)時(shí)按照觀測(cè)時(shí)間的順序來(lái)進(jìn)行切割.如圖5所示,是基于MPI與CUDA的uv-faceting成像并行策略.首先,將輸入數(shù)據(jù)利用MPI進(jìn)程按觀測(cè)時(shí)間分割成多個(gè)ms文件,文件數(shù)等于MPI進(jìn)程數(shù),圖5中同樣是以4個(gè)GPU節(jié)點(diǎn)為例;其次,使用MPI將分割之后的文件數(shù)據(jù)分配到不同的GPU節(jié)點(diǎn)上進(jìn)行并行處理,每個(gè)GPU節(jié)點(diǎn)將完成數(shù)據(jù)的重排序和gridding等數(shù)據(jù)處理,其中,gridding采用上述的CUDA編程進(jìn)行加速處理;然后,將所有GPU節(jié)點(diǎn)的gridding結(jié)果進(jìn)行求和并發(fā)送至主節(jié)點(diǎn)(GPU node0);最后,在主節(jié)點(diǎn)完成分面成像與馬賽克處理獲得最終的圖像.

從上述的描述知,使用本文提出的并行算法的前提是觀測(cè)數(shù)據(jù)以基線為一組并按照時(shí)間排序方式存儲(chǔ),即數(shù)據(jù)按照基線索引來(lái)排序的.然而,ms格式的數(shù)據(jù)一般不是按照這種方式存儲(chǔ)的,因此在數(shù)據(jù)處理之前需要進(jìn)行數(shù)據(jù)重新打包和排序.此外,需要注意的是:所有基線的觀測(cè)時(shí)步的數(shù)量(時(shí)間戳的個(gè)數(shù))并不是都相等,有些基線的某些時(shí)步數(shù)據(jù)有可能在標(biāo)記和校準(zhǔn)的時(shí)候就從Measurement Set中分離出來(lái),所以需要一個(gè)可變長(zhǎng)度的數(shù)組來(lái)存儲(chǔ)每個(gè)基線的數(shù)據(jù).如此,這個(gè)可變長(zhǎng)度的數(shù)組想要傳輸?shù)紾PU上進(jìn)行計(jì)算就非常困難.本文則使用一個(gè)一維數(shù)組存儲(chǔ)每個(gè)基線第1個(gè)時(shí)步位置(即基線的起始索引)來(lái)替代.然而,Measurement Set沒(méi)有存儲(chǔ)基線索引的信息,只存儲(chǔ)了天線的ID號(hào).因此本文給出了一種計(jì)算基線索引的方法.首先,按照基線對(duì)數(shù)據(jù)進(jìn)行分組,每組對(duì)應(yīng)一個(gè)編號(hào).這個(gè)編號(hào)可以使用天線的ID號(hào)來(lái)計(jì)算得到:

其中,B表示天線1 (ant1)和天線2 (ant2)組成的基線;S表示ant1和ant2中最小ID號(hào);Nant表示天線的個(gè)數(shù);Iant1和Iant2分別表示ant1和ant2的ID號(hào).根據(jù)(6)式計(jì)算出基線的編號(hào),然后統(tǒng)計(jì)相同基線編號(hào)下時(shí)步的個(gè)數(shù)(相同基線的采樣數(shù)據(jù)個(gè)數(shù)),并存儲(chǔ)到數(shù)組C中.最后,利用前綴求和方法計(jì)算出每個(gè)基線的起始索引Bstart-index,其數(shù)學(xué)表達(dá)式如下:

其中,ks表示相同基線下時(shí)步的編號(hào),范圍為0 ～I(xiàn)B?1.

圖4 基于CUDA的單GPU節(jié)點(diǎn)uv-faceting成像算法流程圖Fig.4 The flow chart of the uv-faceting imaging algorithm based on the CUDA on a single GPU node

圖5 基于MPI+CUDA的uv-faceting成像并行策略.k表示觀測(cè)時(shí)間戳的開始索引,m表示觀測(cè)時(shí)間編號(hào),pol表示偏振Fig.5 The uv-faceting imaging parallel strategy based on the MPI and CUDA.k is the start index of observation time stamp,m is the numbering of observation time,and pol is the polarization

結(jié)合上述計(jì)算得到的基線索引起始值和相同基線索引下的時(shí)步個(gè)數(shù),我們就可以將讀入的數(shù)據(jù)按照基線索引進(jìn)行重新排序.

4 驗(yàn)證性實(shí)驗(yàn)與分析

本次實(shí)驗(yàn)的數(shù)據(jù)采用甚大陣VLA B和C構(gòu)型組成的望遠(yuǎn)鏡陣列在74 MHz波段觀測(cè)后發(fā)座星系團(tuán)(Coma cluster)得到的數(shù)據(jù),觀測(cè)時(shí)間為:1998-11-22UTC11:19:57.5—1998-12-02UTC22:59:55.0.該觀測(cè)的時(shí)間總長(zhǎng)為905998 s,分為3組觀測(cè).

Coma cluster是一個(gè)超大的星系團(tuán),它包含超過(guò)1000個(gè)已經(jīng)被識(shí)別的星系,視場(chǎng)范圍較大,因此非常適用于本文算法的驗(yàn)證測(cè)試.

實(shí)驗(yàn)的測(cè)試平臺(tái)采用廣州超算中心的天河二號(hào)超級(jí)計(jì)算機(jī),由于測(cè)試需要用到GPU,因此實(shí)驗(yàn)最終在GPU節(jié)點(diǎn)上進(jìn)行.天河二號(hào)目前共有25個(gè)GPU節(jié)點(diǎn)可用,每個(gè)GPU節(jié)點(diǎn)各指標(biāo)參數(shù)如表1所示.基于MPI+CUDA的算法測(cè)試采用10個(gè)節(jié)點(diǎn)進(jìn)行測(cè)試,每個(gè)節(jié)點(diǎn)使用1個(gè)GPU加速卡.基于MPI+OpenMP的算法測(cè)試同樣采用10個(gè)節(jié)點(diǎn)進(jìn)行測(cè)試,但每個(gè)節(jié)點(diǎn)只使用CPU核,每個(gè)節(jié)點(diǎn)共使用20個(gè)CPU核(即20個(gè)線程).

圖6是Coma cluster數(shù)據(jù)進(jìn)行傳統(tǒng)的2維傅里葉變換成像得到的結(jié)果.圖7是使用通用天文應(yīng)用軟件(Common Astronomy Software Applications,CASA)進(jìn)行faceting成像得到的結(jié)果,其中圖像的像素大小設(shè)為1800 × 1800,赤經(jīng)和赤緯方向的圖像單元大小(cellsize)均設(shè)為30′′(arcsecond),分面總數(shù)為25個(gè)分面(赤經(jīng)和赤緯方向各為5個(gè)分面).圖8是本文基于MPI+CUDA的uv-faceting成像得到的結(jié)果,圖9是本文基于MPI+OpenMP的uv-faceting成像得到的結(jié)果,這兩個(gè)測(cè)試的cellsize大小和分面數(shù)量的設(shè)置均與CASA軟件的設(shè)置相同.從視覺(jué)效果可以明顯看出,遠(yuǎn)離相位中心的星系或射電源在圖7、圖8和圖9均能夠得到較好重構(gòu).這說(shuō)明本文提出的這兩種并行算法能夠完成射電干涉陣列大視場(chǎng)成像,結(jié)果與廣泛使用的CASA軟件的結(jié)果基本一致.此外,圖8和圖9在相同的流量范圍內(nèi),重構(gòu)的圖像更加明亮.說(shuō)明我們提出的方法,在提高了運(yùn)行速度的同時(shí)也提高了圖像的信噪比.

表1 天河二號(hào)超級(jí)計(jì)算機(jī)GPU節(jié)點(diǎn)各指標(biāo)參數(shù)Table 1 The indicator parameters of the GPU node of the Tianhe-2 supercomputer

圖6 2維傅里葉變換成像結(jié)果Fig.6 The results of two-dimensional fourier transform imaging

圖7 CASA軟件的成像結(jié)果Fig.7 The imaging results of CASA software

圖8 基于MPI+CUDA的uv-faceting成像結(jié)果Fig.8 The results of uv-faceting imaging based on the MPI and CUDA

圖9 基于MPI+OpenMP的uv-faceting成像結(jié)果Fig.9 The results of uv-faceting imaging based on the MPI and OpenMP

為了更進(jìn)一步驗(yàn)證我們提出的算法的正確性和有效性,我們使用源查找軟件對(duì)圖7、圖8和圖9進(jìn)行源查找,對(duì)比分析它們的圖像結(jié)果的一致性.這里我們使用的源查找軟件是由Paul Hancolk開發(fā)的Aegean軟件[21],該軟件已經(jīng)廣泛應(yīng)用于默奇森廣角陣列MWA的數(shù)據(jù)處理當(dāng)中.源查找的對(duì)比結(jié)果如圖10、圖11和圖12所示,這里我們采用的波束大小(beam size)均為193.20′′× 142.96′′,位置角為?64.46?,這些值由原始觀測(cè)數(shù)據(jù)計(jì)算得到.表2為詳細(xì)的分析對(duì)比結(jié)果.從表中可知,以CASA得到的結(jié)果為基準(zhǔn),MPI+CUDA方法在1個(gè)波束大小的誤差范圍內(nèi)與CASA得到的結(jié)果位置相同的射電源有140顆,相同率約為94.6%,比CASA的結(jié)果多恢復(fù)了11顆源.MPI+OpenMP的方法得到的結(jié)果有141顆射電源與CASA的結(jié)果位置相同,相同率為95.27%,比CASA多恢復(fù)了42顆射電源.這表明本文提出的這兩種方法基本正確.存在差異的原因:我們的方法是將數(shù)據(jù)分割成多份,每份數(shù)據(jù)單獨(dú)進(jìn)行柵格化,再進(jìn)行加權(quán)求和得到最終的柵格化結(jié)果.在加權(quán)求和時(shí),邊界位置部分網(wǎng)格點(diǎn)出現(xiàn)了重復(fù)累加,導(dǎo)致差異的產(chǎn)生.這些差異可以通過(guò)以下方法進(jìn)行修正:以CASA圖像結(jié)果的流量大小為基準(zhǔn),對(duì)我們提出方法的圖像結(jié)果進(jìn)行加權(quán)重,修正它們的流量大小,從而獲得與CASA完全一致的圖像結(jié)果.這項(xiàng)工作是我們下一步的工作內(nèi)容.另一方面,MPI+OpenMP與MPI+CUDA的結(jié)果有159顆射電源的位置相同,說(shuō)明MPI+OpenMP的方法多恢復(fù)了31顆源.表明在正確性上,MPI+CUDA的方法比MPI+OpenMP的方法更加準(zhǔn)確.

在運(yùn)行速度方面,MPI+CUDA的方法消耗的總運(yùn)行時(shí)間最短,CASA消耗的總運(yùn)行時(shí)間最長(zhǎng),相比CASA、MPI+CUDA的方法提高了將近13倍.雖然MPI+OpenMP與MPI+CUDA消耗的總運(yùn)行時(shí)間差不多,它們只相差不到3 s.但在gridding方面,MPI+CUDA的方法消耗的時(shí)間為1.16 s,MPI+OpenMP消耗的時(shí)間為3.92 s,說(shuō)明MPI+OpenMP是MPI+CUDA的3.38倍左右.綜上分析,表明本文提出的算法在運(yùn)行速度上有很大的提高,不管是在運(yùn)行速度還是正確性上MPI+CUDA方法更優(yōu),但也有上升的空間.

圖10 圖7與圖8源搜尋結(jié)果對(duì)比圖Fig.10 Comparison of source finding results of Fig.7 and Fig.8

圖11 圖7與圖9源搜尋結(jié)果對(duì)比圖Fig.11 Comparison of source finding results of Fig.7 and Fig.9

圖12 圖8與圖9源搜尋結(jié)果對(duì)比圖Fig.12 Comparison of source finding results of Fig.8 and Fig.9

表2 源查找結(jié)果對(duì)比Table 2 Comparison of source finding results

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

5.1 單節(jié)點(diǎn)

單節(jié)點(diǎn)的測(cè)試平臺(tái),我們采用上海天文臺(tái)區(qū)域中心原理樣機(jī)的GPU節(jié)點(diǎn)、澳大利亞Pawsey超算中心Galaxy超級(jí)計(jì)算機(jī)的GPU節(jié)點(diǎn)和上述天河二號(hào)的GPU節(jié)點(diǎn)進(jìn)行測(cè)試.上海天文臺(tái)GPU節(jié)點(diǎn)的各項(xiàng)指標(biāo)參數(shù)見(jiàn)表3,Galaxy超級(jí)計(jì)算機(jī)的GPU節(jié)點(diǎn)各項(xiàng)指標(biāo)參數(shù)見(jiàn)表4.測(cè)試數(shù)據(jù)采用第4節(jié)中Coma cluster的觀測(cè)數(shù)據(jù),分面總數(shù)設(shè)置為49個(gè)分面,即赤經(jīng)和赤緯方向各7個(gè)分面,其余的參數(shù)設(shè)置與第4節(jié)中的測(cè)試相同.我們主要通過(guò)改變環(huán)境變量OMP-NUM-THREADS來(lái)使用不同的線程數(shù)量測(cè)試基于MPI+OpenMP的uv-faceting成像算法的性能.然后,通過(guò)使用不同型號(hào)的GPU來(lái)測(cè)試基于MPI+CUDA的uv-faceting成像算法的性能.

如圖13所示,在56個(gè)線程范圍內(nèi),隨著線程數(shù)量的增加,gridding所消耗的時(shí)間隨之遞減.但沒(méi)有呈現(xiàn)良好的線性性,不管是單精度還是雙精度計(jì)算.極有可能是分面的數(shù)量沒(méi)有均勻分配給CPU核,導(dǎo)致了各個(gè)線程之間的負(fù)載不均衡.當(dāng)線程的數(shù)量接近分面總數(shù)時(shí),gridding所消耗的時(shí)間最少;不管是單精度還是雙精度情況下,相比沒(méi)有進(jìn)行并行改造的運(yùn)行速度均增加了20倍左右.表明本文提出的基于MPI+OpenMP的uvfaceting成像算法的并行優(yōu)化是有效的.另外,我們還對(duì)比測(cè)試了CPU和GPU之間以及不同GPU之間的性能.隨著GPU加速卡性能的提升,gridding消耗的時(shí)間變短.如圖13特斯拉(Tesla)V100所消耗的時(shí)間最短,無(wú)論是單精度還是雙精度相比沒(méi)有進(jìn)行并行改造的運(yùn)行速度增加了60倍左右.單精度情況下,Tesla K80 gridding所消耗的時(shí)間與CPU使用48個(gè)線程所消耗的時(shí)間相當(dāng),Tesla K20X gridding所消耗的時(shí)間與CPU使用40個(gè)線程的相當(dāng).雙精度情況下,Tesla K80 gridding所消耗的時(shí)間與CPU使用8個(gè)線程所消耗的時(shí)間相當(dāng),Tesla K20X gridding所消耗的時(shí)間與CPU使用4個(gè)線程的相當(dāng).表明本文提出的基于MPI+CUDA的uv-faceting成像算法是有效的.同時(shí)也表明,在一定的線程范圍內(nèi),基于MPI+CUDA的算法比基于MPI+OpenMP的算法性能要好,并且隨著GPU加速卡的提高,這個(gè)范圍的上限值會(huì)隨之增大.

表3 上海天文臺(tái)區(qū)域中心原理樣機(jī)GPU節(jié)點(diǎn)各指標(biāo)參數(shù)Table 3 The indicator parameters of the GPU node of the Shanghai Astronomical Observatory Regional Center Prototype

表4 Galaxy超級(jí)計(jì)算機(jī)GPU節(jié)點(diǎn)各指標(biāo)參數(shù)Table 4 The indicator parameters of the GPU node of the Galaxy supercomputer

5.2 多節(jié)點(diǎn)

多節(jié)點(diǎn)的測(cè)試數(shù)據(jù)我們采用擴(kuò)展甚大陣望遠(yuǎn)鏡(Expanded Very Large Array,EVLA)的觀測(cè)數(shù)據(jù),選擇了Sanjay Bhatnagar的超新星遺跡G55.7+3.4 L波段的觀測(cè)數(shù)據(jù).該觀測(cè)于2010-08-23UTC01:00:25開始,2010-08-23UTC08:25:16.5結(jié)束,歷時(shí)約7小時(shí)25分.首先我們使用CASA軟件對(duì)G55.7+3.4源的原始觀測(cè)數(shù)據(jù)進(jìn)行校準(zhǔn),并導(dǎo)出校準(zhǔn)后的數(shù)據(jù)作為本次測(cè)試的輸入數(shù)據(jù).

本次測(cè)試的平臺(tái)采用天河二號(hào)的GPU節(jié)點(diǎn),如表1所示.測(cè)試數(shù)據(jù)的參數(shù)設(shè)置:像素大小設(shè)為1024×1024,赤經(jīng)和赤緯方向的圖像單元大小(cellsize)均設(shè)為8′′,分面總數(shù)為16個(gè)分面(赤經(jīng)和赤緯方向各為4個(gè)分面).測(cè)試方法:計(jì)算節(jié)點(diǎn)從1到10自然數(shù)遞增,每個(gè)節(jié)點(diǎn)數(shù)測(cè)試10次,獲取數(shù)據(jù)加載的時(shí)間和gridding的時(shí)間,并取平均值和標(biāo)準(zhǔn)差.基于MPI+CUDA的方法每個(gè)節(jié)點(diǎn)使用1個(gè)GPU加速卡,基于MPI+OpenMP的方法每個(gè)節(jié)點(diǎn)使用20個(gè)CPU核.

圖13 單節(jié)點(diǎn)性能測(cè)試結(jié)果Fig.13 The results of performance tests on a single node

數(shù)據(jù)讀取(或加載)的耗時(shí)情況如圖14所示.在GPU上運(yùn)行時(shí),在1–4個(gè)節(jié)點(diǎn)范圍內(nèi),無(wú)論是單精度情況下還是雙精度情況下,隨著節(jié)點(diǎn)數(shù)的增加所消耗的時(shí)間均接近線性遞減.當(dāng)節(jié)點(diǎn)數(shù)超過(guò)4個(gè)之后,隨著節(jié)點(diǎn)數(shù)的增加,所消耗的時(shí)間減少得比較緩慢,在8個(gè)節(jié)點(diǎn)數(shù)之后達(dá)到平衡.在CPU上運(yùn)行時(shí),只有在2個(gè)節(jié)點(diǎn)范圍內(nèi)接近線性遞減.當(dāng)節(jié)點(diǎn)數(shù)超過(guò)2個(gè)之后,變化緩慢,同樣在8個(gè)節(jié)點(diǎn)數(shù)之后達(dá)到平衡.表明本文提出的基于MPI+CUDA方法和基于MPI+OpenMP方法在一定計(jì)算節(jié)點(diǎn)范圍內(nèi)具有良好的可拓展性,基于MPI+CUDA方法的范圍更大,拓展性更優(yōu).另一方面,在CPU上運(yùn)行時(shí),所消耗的數(shù)據(jù)加載時(shí)間相對(duì)較少,單精度和雙精度所消耗的時(shí)間基本相同.此外,CPU上運(yùn)行的標(biāo)準(zhǔn)差較GPU上運(yùn)行的小,運(yùn)行相對(duì)穩(wěn)定.在GPU上運(yùn)行,所消耗的時(shí)間較大一些,且不夠穩(wěn)定.表明在數(shù)據(jù)加載方面,基于MPI+OpenMP的方法比基于MPI+CUDA的方法穩(wěn)定性強(qiáng)且耗時(shí)短.

gridding的耗時(shí)情況如圖15所示.從圖中可以看出,無(wú)論是CPU還是GPU上運(yùn)行,gridding所消耗的時(shí)間均隨著節(jié)點(diǎn)數(shù)的增加逐漸減少,但是均沒(méi)有呈現(xiàn)明顯的線性性,在8個(gè)節(jié)點(diǎn)數(shù)之后達(dá)到平衡.1個(gè)節(jié)點(diǎn)數(shù)所消耗的時(shí)間是10個(gè)節(jié)點(diǎn)數(shù)的9倍左右.表明本文提出的兩種方法在整體上具有良好的拓展性.在單精度情況下,在GPU上運(yùn)行所消耗的時(shí)間較短,標(biāo)準(zhǔn)差幾乎為零,相對(duì)較穩(wěn)定.在雙精度情況下,在CPU上的運(yùn)行時(shí)間較短,但GPU的運(yùn)行穩(wěn)定性更強(qiáng).表明在gridding方面,基于MPI+CUDA的方法比基于MPI+OpenMP的方法更加穩(wěn)定.

圖14 數(shù)據(jù)加載耗時(shí)曲線Fig.14 Data loading consuming time curve

圖15 gridding耗時(shí)曲線Fig.15 Gridding consuming time curve

6 總結(jié)與展望

本文主要介紹了uv-faceting的工程實(shí)現(xiàn)流程,進(jìn)而分析了uv-faceting主要步驟的耗時(shí)情況,明確了數(shù)據(jù)讀取和gridding這兩個(gè)步驟是uv-faceting中最耗時(shí)的部分.基于此,本文提出了兩種并行優(yōu)化算法:基于MPI+OpenMP的uv-faceting成像算法和基于MPI+CUDA的uv-faceting成像算法.驗(yàn)證性實(shí)驗(yàn)的結(jié)果表明本文提出的兩種算法均能夠正確實(shí)現(xiàn)uv-faceting成像,結(jié)果與目前廣泛使用的CASA軟件得到的結(jié)果基本一致.但還存在略微差異,這是我們將來(lái)工作需要完善的地方之一.性能測(cè)試方面,單節(jié)點(diǎn)的測(cè)試表明本文提出的兩種算法是有效的,在一定線程范圍內(nèi),基于MPI+CUDA的方法比基于MPI+OpenMP的方法性能更好.其中,反映了基于MPI+OpenMP的方法線程之間負(fù)載均衡不是很好,多線程的拓展性有待優(yōu)化和提高.多節(jié)點(diǎn)的測(cè)試表明,在數(shù)據(jù)讀取方面,本文的兩種方法在一定的節(jié)點(diǎn)數(shù)范圍內(nèi)具有較好的可擴(kuò)展性,基于MPI+OpenMP的方法比基于MPI+CUDA的方法更加穩(wěn)定.在gridding方面,表明本文提出的兩種方法在整體上有良好的可拓展性,但在局部上均沒(méi)有展現(xiàn)出明顯的線性性,說(shuō)明節(jié)點(diǎn)之間的負(fù)載均衡不是很好.因此,在后續(xù)的工作中,需要完善多節(jié)點(diǎn)之間的可拓展性性能.