基于Intel至強眾核加速積分法疊前時間偏移技術的實現

陳 維

(東方地球物理公司 油藏地球物理研究中心,涿州 072751)

基于Intel至強眾核加速積分法疊前時間偏移技術的實現

陳 維

(東方地球物理公司 油藏地球物理研究中心,涿州 072751)

為進一步提升積分法疊前時間偏移(PSTM)程序的運行效率,利用Intel至強眾核(MIC)技術,將該程序移植到至強眾核平臺上,以實現對該程序的加速。 GeoEast系統下的PSTM模塊CPU版本采用MPI+PThread多線程方式實現程序運行的并行計算。在CPU+MIC協同計算方式下，繼承上述的并行方式，同時，成像空間數據在節點間和節點內的CPU和MIC之間采用交替分割方式，并優化CPU與MIC之間的數據異步流水線式傳輸。 通過CPU+MIC協同計算優化，PSTM模塊運行效率相比CPU版本提升3.8倍，取得了預期的效果，也驗證了CPU+MIC協同計算加速PSTM的可行性。

疊前時間偏移; 眾核; 卸載; 移植; 多線程

0 引言

隨著地震勘探技術的進步，勘探采集得到的數據已經進入TB級時代，絕大多數三維項目都達到幾個TB，甚至幾十個TB的海量數據。同時由于勘探區域復雜的地下構造，要求每個項目都需要使用疊前時間偏移(PSTM)方法(常用的是Kirchhoff積分法PSTM)，對于地層傾角較大，地下橫向速度變化不劇烈的情況下，PSTM能夠獲得較好的成像效果[1]。 但由于疊前時間偏移是一個計算量極大的算法，是整個常規處理流程中最耗時的模塊，約占整個流程40%的處理時間。因此，國際上流行的幾大地震資料處理系統(GeoEast，Omega，CGG等)都在不斷優化加速PSTM模塊，以滿足海量地震資料數據處理需求。計算機硬件廠商面對高性能計算領域越來越大的計算需求，也在不斷推出新的技術。 面對計算量越來越大的高性能計算需求，陸續出現了英偉達公司的GPGPU加速器以及英特爾公司的MIC架構的至強融合Phi。GPGPU采用超輕量線程，使用CUDA并行開發環境對程序進行并行優化，在眾多高性能計算領域獲得應用，在地震資料處理領域，疊前時間偏移、疊前深度偏移等模塊都有GPGPU版本，并已經應用于實際生產。英特爾公司作為CPU市場的領導者，采用了與Intel CPU相同架構的x86技術于2011年推出開發代號為Knight corner的MIC(Many Integrated Core - 眾核)架構至強Phi眾核協處理器，每個處理器內部含62個核，每個核可以啟動4個硬件線程，用戶應用可使用60個核共240個線程。MIC的最大特點在于程序源代碼與CPU代碼兼容，減少了維護成本，并可以使用傳統編程模型，減少了學習成本[2]。MIC的出現，為高性能計算的并行程序優化提供了易于實現和易于維護的平臺，也為地震資料處理中類似疊前時間偏移等非常耗時模塊的加速優化提供了新的解決方案。 在GeoEast地震資料處理解釋一體化系統開發過程中，對于疊前時間偏移等非常耗時的模塊，其優化工作一直在持續進行。積分法疊前時間偏移模塊經歷了從進程版MPI(Message Passing Interface - 消息通訊接口，一種高性能計算領域常用的并行計算環境)并行，到MPI+多線程并行，以及CPU+GPU協同計算的持續不斷的并行優化。新的并行加速技術的出現與應用促進了地震處理高耗時模塊技術的進步及計算效率的不斷提升，不斷滿足海量數據處理生產的需求。這里所述的工作在Knight Corner版本的MIC上完成，通過在該平臺上對GeoEast系統原MPI+多線程CPU版本積分法疊前時間偏移并行模塊進行MIC平臺的移植加速，取得了比較好的加速效果，在相同硬件環境下通過增加兩個MIC卡，運行效率相比原來的模塊提升了3.8倍。

1 CPU版本疊前時間偏移并行算法

由于當前的機群系統同時具有共享存儲和分布式存儲兩級結構，節點內采用OpenMP多線程，而節點間采用MPI多進程的混合模式可以更好地利用機群的特性[3]；為了解決I/O對偏移效率的影響，王霖等[4]提出了流水線方式的PSTM并行計算方法；在此基礎上，馬召貴等[5]進一步提出了PSTM多級并行優化技術，采用MPI多進程并行，數據I/O與計算異步并行及單進程內多線程成像計算方法。 GeoEast系統MPI+多線程CPU版本疊前時間偏移并行模塊同樣采用了節點內多線程以及節點間使用MPI通訊的并行方式。程序使用了基于輸出并行的方法，其基本思想是把輸出成像空間進行劃分并分配給各個計算節點，各個計算節點在各自所分配的成像空間內對所有輸入數據進行偏移，處理完畢后將所有結果按序進行輸出。其過程如下：

循環所有地震道：

for (all input traces)

主節點廣播地震道到各個計算節點

各個計算節點預處理地震道

對節點內成像空間所有輸出道

for (all output traces )

image(out) += 輸入道偏移結果

為了提高數據分發效率，本模塊采用分組數據交換方式，每組包含若干個節點，并為每組分配一個頭節點，該節點負責從磁盤讀取地震數據，并將數據分發給該組的各個節點。在每個節點內，偏移作業包含一個MPI主進程，多個FFT/反FFT線程(AntiDraper)以及多個偏移線程(PSTMKernel)，一個數據接收線程，每組頭節點還包含一個數據讀取線程。程序采用多緩沖區異步計算方式，實現數據準備與偏移計算同時進行，當一個數據已經準備好的緩沖區開始用于偏移計算時，另外一個緩沖區立即開始準備新的用于偏移的數據。

MPI主進程主要負責以下工作：①偏移成像空間分割和最終成像空間的拼接輸出；②線程啟動與結束；③數據獲取與任務控制。

1.1 輸入道的數據分發過程

1) 頭節點分批從磁盤讀入數據，每批的數據大小是10 000道數據，并分發給組內其他節點。

2) AntiDraper線程將10 000道數據分割成更小的數據塊，一塊一塊地依次處理。本模塊設定每個數據塊的道數是CPU數cpu_num。以cpu_num個數據道為最小單位，經AntiDraper線程處理完后，接著由PSTMKernel進行偏移處理。

1.2 輸出成像空間的分割策略

將成像空間交替分配給nnode個節點，在每個節點上再將成像空間交替分配給cpu_num個偏移線程，分割規律如圖1所示(以3個節點，每個節點5個偏移線程為例)。在圖1中，VELS表示成像空間，Rk(i)表示第i節點上的MPI主進程，Thd(i)表示第i個線程。此處在MPI進程之間的成像空間分割的特點是：

圖1 MPI+多線程CPU版本成像空間分割示意圖Fig.1 Image map for CPU MPI+multi thread PSMT

對于total-trace個輸出道大小的成像空間，第0號 ～cpu_num-1號分配給0號進程，第cpu_num號～2*cpu_num-1號則分給1號進程，第2*cpu_num號～3*cpu_num-1則分給2號進程；第3*cpu_num號～4*cpu_num-1又分給0號進程，如此直到將所有total_trace個輸出道交替分割完畢。程序這樣交替分割的原因是為了每個節點的計算負載平衡。

1.3 AntiDraper線程與PSTMKernel線程的協同工作

作業啟動時，創建cpu_num個AntiDraper線程和cpu_num個PSTMKernel線程。每個AntiDraper線程負責一個輸入道數據的處理，每執行完一次輸入道數據塊的FFT/iFFT的計算工作之后，就向每個PSTMKernel線程的輸入緩沖區隊列中寫入自己的索引號，通知PSTMKernel線程這個緩沖區的一個輸入道數據轉換完畢，可以開始偏移計算。同時對該緩沖區禁止寫，直到其被cpu_num個PSTMKernel線程各訪問一次后，才允許寫。每個PSTMKernel線程查詢自己緩沖區隊列中可以計算的輸入道數據。每計算完一個緩沖區，相應的PSTMKernel線程就釋放一次該緩沖區。待到該緩沖區被釋放cpu_num次后，AntiDraper線程又獲得向該緩沖區寫的權利。如此彼此協同直至所有輸入道數據都處理完畢。

cpu_num個AntiDraper線程與cpu_num個PSTMKernel線程之間協作的過程如下：當某個AntiDraper線程處理向其對應的緩沖區中寫完一個輸入道數據時，該AntiDraper線程即在所有的PSTMKernel線程的輸入緩沖區索引隊列中寫入該緩沖區索引，表示該緩沖區準備完畢。每個PSTMKernel線程則通過讀自己的緩沖區索引隊列來確定哪個緩沖區可以使用了。每個PSTMKernel使用完某個緩沖區后即釋放一次資源，待該緩沖區被釋放cpu_num次后，該緩沖區就重新被對應的AntiDraper線程獲取，再次向其寫入新的數據，供PSTMKernel使用。如此往復，直到所有輸入道數據都處理完畢。

每個AntiDraper線程在自己的輸入緩沖區中讀取存放10 000道輸入道數據的地址，獲取緩沖區中的道數traceNum。之后每個AntiDraper線程就只訪問自己負責的輸入道的數據，計算后放入對應的緩沖區之中。第_cpuIndex號AntiDraper線程將緩沖區索引號_cpuIndex寫入所有的PstmKernel線程的輸入緩沖區索引隊列中，告訴所有的PstmKernel線程第_cpuIndex緩沖區數據準備完畢。之后第_cpuIndex號AntiDraper線程就進入等待狀態，等待所有PstmKernel線程都使用了第_cpuIndex緩沖區中的數據，開始計算新的一道輸入道數據。待traceNum道數據處理完畢時，再批量讀入下一個10 000道數據。每個PstmKernel線程則是不斷的訪問自己的緩沖區索引隊列，當有有效值時，返回值表示那個緩沖區的數據已經準備好，可以進行偏移計算。偏移計算完成之后，將該緩沖區釋放一次，信號量增加1，接著進入下一個緩沖區數據的處理。

2 CPU與MIC協同計算PSTM加速

2.1 MIC環境下的編程

與GPU類似，MIC以PCI卡的方式插在主機的PCI插槽上，通過Intel提供的驅動程序實現主機與MIC卡之間的數據通訊與控制。基本的編程方式采用類似OpenMP的指導語句offload實現程序在MIC上運行。 例如：

float reduction(float *data,size_t size)

{

float ret=0.0;

#pragma offload target(mic)

in(data:length(size))

for (int i=0; i

ret += data[i];

}

return ret;

}

在上面的代碼中，#pragma offload 指定了緊隨其后的程序段要在MIC上運行，target(mic)指定在哪個MIC上運行。這里所述在MIC上程序的并行計算使用OpenMP的多線程方式，線程啟動方式如下：

#pragma omp parallel default(shared) num_threads(線程個數)

{

并行計算疊前時間偏移線程

}

2.2 CPU/MIC協同工作下程序的優化

由于MIC架構的至強融合處理器以卡的形式插在PCI插槽上，與CPU之間的數據交換受到PCI帶寬的限制，同時，MIC上的內存相對主機要小得多，而地震處理中大部分的高耗時處理模塊都需要比較大的內存，因此不僅要考慮主機與MIC卡之間的數據傳輸與MIC計算的協同問題，同時也要對原來的CPU程序進行一定的修改，以減少對內存的需求。另外，MIC架構支持512位數據帶寬，可以同時進行16個浮點數據的計算，因此，如何更好地發揮這一優勢，也是程序優化重點關注的地方。 與此同時，由于CPU的主頻與核數與MIC上的CPU主頻與核數的差異，會導致計算效率的不同。 綜上所述，主機CPU與MIC之間要實現高效的協同工作，需要在以下幾個方面進行程序的優化：

1)合理的任務分割，以達到較好的任務均衡。

2)采用“軟流水”技術，使數據傳輸與計算重疊，盡可能避免數據等待。

3)程序循環的向量化，即使用#pragma simd指導語句對程序循環進行向量化，充分利用MIC的512位數據帶寬的優勢。

4)優化數據傳輸，盡量較少數據傳輸次數。

2.3 積分法疊前時間偏移的CPU和MIC的協同計算優化

根據該程序的特點，對其優化主要進行以下各項工作：①成像空間的分割，也即進行任務的劃分；②數據傳輸的設計，包含每次傳送多少地震道用于計算，計算與數據傳輸如何協同；③偏移程序核心計算部分的向量化。

這里對積分法疊前時間偏移進行MIC加速優化的基本思想就是將以前CPU的計算任務分出一部分給MIC處理，將CPU+MIC作為一個節點來看待。如果使用雙MIC卡計算，則在一個物理節點上啟動兩個MPI進程，進程號即是其可以利用的MIC設備號，在機群的雙路節點情況下，將CPU也分成兩部分均勻分給兩個MPI進程。

在CPU+MIC的情況下，將CPU+MIC認為是一個節點，該節點上有cpu_num個AntiDraper線程用于輸入道數據的FFT/iFFT變換(iFFT即反變換)，cpu_num-1個在CPU上運行的PstmKernel線程，mic_num個在MIC上運行的PstmKernel線程，mic_num個PstmKernel的MIC線程由一個CPU線程管理，所以程序實際在HOST端創建了2*cpu_num個線程。

2.3.1 節點進程間成像空間的分割

成像空間的分割是在節點進程之間交替分割。如圖2所示，一個連續的vels塊按照對3取模運算，均勻分割到所有的MPI進程上，在圖2中Rk(0)表示0號進程，Rk(1)表示1號進程，Rk(2)表示2號進程。

圖2 節點進程間VELS分割Fig.2 VELS for each MPI Process

2.3.2 AntiDraper與PstmKernel協同的過程

定義兩個可容納N*cpu_num道輸入數據的緩沖區，兩個緩沖區輪流向PstmKernel線程供應數據。cpu_num個AntiDraper OMP線程經過N次循環就可以填滿一個緩沖區。當該緩沖區填滿后，即使用異步Offload上傳數據到MIC上去。當異步Offload返回時，AntiDraper就向每個PstmKernel線程的緩沖區索引隊列中寫入相應的緩沖區號，表示該緩沖區數據更新完畢。AntiDraper線程開始等待所有PstmKernel線程釋放資源的信號，當有cpu_num次釋放后，AntiDraper中index++，即開始進行下一個緩沖區的填充工作(此處，index為緩沖區索引)。此時PstmKernel線程還在計算第index個緩沖區中的數據。實現了PstmKernel的計算與AntiDraper的計算及數據向MIC的傳遞重疊。當AntiDraper線程將index++組數據寫入第(index++)%2緩沖區時，AntiDraper即進入等待狀態，等待PstmKernel釋放資源的信號。

AntiDraper線程通過OpenMP派生線程，實現并行輸入道處理。如圖3所示，AntiDraper中OpenMP派生10個線程，每次處理10道輸入道數據，進行N次之后，即完成了N*10道輸入道數據的處理，填充滿一個緩沖區。AntiDraper之后將緩沖區中的數據復制到MIC上去，即進入等待狀態。當收到PstmKernel線程通知，該緩沖區正被處理時，AntiDraper開始下一個緩沖區的填充工作。

圖3 AntiDraper 和 PstmKernel線程協同Fig.3 Cooperation between AntiDraperand PstmKernel

PstmKernel線程中0號～8號線程直接在CPU上執行，第9號線程則負責MIC上運行程序段的管理工作，并將執行代碼Offload到MIC上去，創建mic_num個線程并行處理各自所負責的緩沖區。

2.3.3 節點內成像空間數據分割規律

節點內CPU和MIC之間的成像空間數據分割仍然保持交替分割的策略。以cpu_part+mic_part為最小單位來分割，cpu_part和mic_part分別表示每個節點負責的成像空間在CPU上和MIC上計算的比例，mic_part/(cpu_part+mic_part)表示MIC所占成像空間的比例。

1)每個節點內CPU和MIC之間的成像空間分割。對于mic_part=1，cpu_part=1的情況是以2為分割的最小單位，所有mod_index % (mic_part + cpu_part) < mic_part 的成像空間都給MIC(這里的mod_index為成像空間中的道索引)，由MIC線程負責計算；所有mod_index % (mic_part + cpu_part) >= mic_part && mod_index % (mic_part + cpu_part) < (mic_part + cpu_part)的成像空間都給CPU，由CPU線程負責計算。同時MIC相比CPU的計算比例為mic_part:cpu_part。

2)CPU和MIC上線程之間成像空間的分割。線程之間的成像空間的分割，則是直接采取取模的方式分割。MIC負責的第1個成像空間即給第0號MIC線程，依次第mic_num個成像空間則給第mic_num-1號線程；第mic_num+1個成像空間再給0號MIC線程，如此直到所有分割完畢；CPU 上成像空間的分割同樣按取模方式分割，如圖4所示。

圖4 MIC與CPU之間成像空間分割示意圖Fig.4 Image map for MIC and CPU(a)VELS;(b)VELS on MIC;(c)VELS on CPU

對于不能整除的情況，要向上去整數，保證不丟棄成像空間。圖4(a)的vels表示每個節點所負責的成像空間，圖4(b)表示MIC負責的部分，圖4(c)表示CPU負責的部分。mic_i表示第i個mic線程所負責的成像空間部分。cpu_i表示第i個cpu線程所負責的成像空間部分。

3 試驗結果

本試驗在GeoEast2.6系統下運行MIC優化后的PSTM版本，硬件環境如下：CPU為SandyBridge 架構，雙路，16核,內存為64GB,2個MIC卡。試驗數據為：某工區三維數據，7s記錄，4 ms采樣,偏移參數如表1所示。

表1 CPU版本與MIC版本運行效率對比Tab.1 Run time compare between CPU and MIC version

試驗方式：在同一臺機器上，使用純CPU版本的PSTM模塊運行該作業，然后使用MIC優化的版本運行相同的作業，比較兩個不同版本運行同一作業的運行時間。 運行結果為，CPU+MIC協同計算,相比CPU版本效率提升3.8倍。

4 結束語

由于MIC架構是基于x86的架構，與主流Intel的CPU架構兼容，因此，針對MIC的代碼優化方法也基本上適用CPU代碼。(如#pragma simd的向量優化技術，OpenMP多線程技術等)，極大減少了代碼維護成本。然而，由于MIC卡仍然是PCI卡，數據需要在主存和MIC卡之間傳送，同時MIC卡上的內存相對較小，進行程序優化時，不僅要考慮數據傳輸對計算的影響，還要考慮如何修改程序，使其適合在小內存的MIC卡上運行，對于需要大內存的地震模塊的移植還是需要一定的工作量和技術，甚至可能需要對算法進行一定的修改。

從這里的優化移植試驗可以看出，MIC眾核技術程序優化對代碼相對較小的修改以及較好的代碼兼容性，對于地震數據處理軟件中那些高耗時模塊的優化具有較大的優勢，其在GeoEast系統下的移植也表明其很容易集成到系統中。本次研究也為對其它耗時模塊在眾核平臺上的移植提供了有益的經驗，該技術也將會隨著更多耗時模塊的移植在海量地震數據處理中得到更多的應用。

[1]王棣，王華忠，馬在田，等.疊前時間偏移方法綜述[J].勘探地球物理進展，2004 ，27(5) ：313-320.WANG D,WANG H Z,MA Z T,et al.Review of prestack time migration methods[J].Progress in Exploration Geophysics,2004,27(5):313-320.(In Chinese)

[2]王恩東,張清,沈鉑,等.MIC高性能計算編程指南[M].北京：中國水利水電出版社，2012.WANG E D,ZHANG Q,SHEN B,et al.MIC high performance computing programming guide[M].Beijing:China Water & Power Press,2012.(In Chinese)

[3]趙永華，遲學斌.基于SMP集群MPI+OpenMP混合編程模型及有效實現[J].微電子學與計算機，2005，22(10)：7-11.ZHAO Y H,CHI X B.MPI+OpenMP hybrid paradigms and efficient implementation base on SMP clusters[J].microelectronic & computer,2005,22(10):7-11.(In Chinese)

[4]王霖,晏海華.一種疊前時間偏移并行模式的流水線改進方法[J].微電子學與計算機,2008,25(10):54-57.WANG L,YAN H H.Method of improving performance of kirchhoff prestack time migration parallel pattern through pipelining[J].microelectronic & computer,2008,25(10):54-57.(In Chinese)

[5]馬召貴,趙改善,武港山等.起伏地表疊前時間偏移的多級并行優化技術[J].石油物探,2013,52(3):280-287.MA Z G,ZHAO G S,WU G S,et al.Multi-level parallel optimization technique for prestack time migration from rugged topography[J].geophysical prospecting for petroleum,2013,52(3):280-287.(In Chinese)

The implementation of MIC usage in accelerating Kirchhoff prestack time migration

CHEN Wei

(Research Center of Reservoir Geophysics,BGP Inc.CNPC,Zhuozhou 072751,China)

In order to improve Kirchhoff pre-stack time migration process efficiency,the program is ported onto Intel Xeon multiple integrated core(MIC) platform.The PSTM in GeoEast system uses MPI+Pthread method in parallel processing.In CPU+MIC cooperative computing,in addition to MPI+PThread,the migration imaging data are put on each node by turn,and also on CPU and MIC by turn.In addition,using asynchronous data transfer between CPU and MIC.Through the CPU+MIC cooperative computing optimization,compared with CPU version,the MIC version of PSTM got the 3.8 times shorter computing time.

PSTM; MIC; offload; porting; multithread

2015-10-09 改回日期：2015-10-30

國家科技重大專項(2011ZX05019-006)

陳維(1963-)，男，高級工程師，從事石油物探高性能計算并行軟件開發，E-mail：weic@cnpc.com.cn。

1001-1749(2016)06-0815-06

P 631.4