MT Occam反演的CPU/GPU異構(gòu)混合并行算法研究

劉 羽,熊壬浩,肖 熠

(桂林理工大學(xué)信息科學(xué)與工程學(xué)院,廣西桂林541004)

大量實(shí)際計(jì)算經(jīng)驗(yàn)表明,電磁法正反演計(jì)算量大,非常耗時(shí),計(jì)算速度難以滿足當(dāng)前地質(zhì)勘查的要求。文獻(xiàn)[1]利用單臺(tái)主流PC機(jī)對(duì)101×50網(wǎng)格、20觀測(cè)頻率、21測(cè)點(diǎn)的二維航空電磁法數(shù)據(jù)進(jìn)行正演模擬,耗時(shí)4.0～4.5h。文獻(xiàn)[2]采用曙光TC5000A單節(jié)點(diǎn)對(duì)34×34×33網(wǎng)格、16觀測(cè)頻率、25測(cè)點(diǎn)的三維MT數(shù)據(jù)進(jìn)行反演,運(yùn)行時(shí)間長(zhǎng)達(dá)47.6h。影響計(jì)算速度的主要因素有:①顯式的雅可比矩陣計(jì)算耗時(shí)太多;②當(dāng)采用有限單元法求解正演響應(yīng)時(shí),密集的剖分網(wǎng)格會(huì)形成巨大的系數(shù)矩陣導(dǎo)致計(jì)算量陡增;③有些反演方法采用自適應(yīng)算法求取正則因子,導(dǎo)致迭代次數(shù)和總體計(jì)算量增加。圍繞電磁法計(jì)算速度的提升方法,國(guó)內(nèi)外地球物理工作者做了很多算法上的改進(jìn)[3-5],特別是在并行計(jì)算方面,進(jìn)行了很多探索。2006年,譚捍東等[6]基于頻率劃分,采用MPI平臺(tái)實(shí)現(xiàn)了大地電磁三維正演的并行計(jì)算;同年,劉羽等[7-8]基于頻率和拉格朗日乘子劃分,在并行計(jì)算機(jī)(PVM)平臺(tái)上,實(shí)現(xiàn)了二維MT Occam并行反演;陳露軍[9]實(shí)現(xiàn)了基于MPI的井地三維電磁數(shù)值模擬;2009年,SIRIPUNVARAPORN等[10]利用頻率劃分方式實(shí)現(xiàn)了包含垂直磁場(chǎng)轉(zhuǎn)換函數(shù)的三維數(shù)據(jù)空間并行反演;2010—2012年,胡祥云等將并行計(jì)算應(yīng)用于二維大地電磁正演,然后擴(kuò)展到三維數(shù)據(jù)空間反演,取得了較好的并行效果[2,11]。

前述并行計(jì)算成果均采用了消息傳遞模型和基于頻率劃分的大粒度并行方法,當(dāng)觀測(cè)頻率數(shù)較少時(shí),此種方式不具備可擴(kuò)展并行性,計(jì)算性能難以通過(guò)增加節(jié)點(diǎn)數(shù)得到相應(yīng)提高。例如,文獻(xiàn)[6]在4節(jié)點(diǎn)集群上獲得2.42的加速比(2D prism模型);文獻(xiàn)[2]和文獻(xiàn)[8]在8節(jié)點(diǎn)集群上分別獲得5.69和5.65的加速比;文獻(xiàn)[10]在8節(jié)點(diǎn)和16節(jié)點(diǎn)集群上分別獲得4.5和7.3的加速比(采用PCG求解器)?？梢?jiàn),要想獲得更高計(jì)算性能,僅靠擴(kuò)大集群規(guī)模很難奏效。

本文將傳統(tǒng)的CPU集群擴(kuò)展為CPU/GPU異構(gòu)集群,以MT Occam反演為例,基于不同粒度并行分量的挖掘,實(shí)現(xiàn)計(jì)算任務(wù)的不同層次映射,從而在大粒度并行的基礎(chǔ)上進(jìn)一步提升計(jì)算性能。首先給出MT Occam反演的理論方法,然后分析反演方法的并行特征和現(xiàn)有并行方法存在的性能瓶頸,并針對(duì)性地提出混合并行方案,利用2種規(guī)模的5個(gè)理論模型進(jìn)行了反演測(cè)試,給出了方案的實(shí)際效果,最后討論了方案的適用性并給出結(jié)論。

1 反演方法

MT Occam反演目標(biāo)函數(shù)可表示為[12-13]:

(1)

式中:μ為拉格朗日乘子;d為觀測(cè)數(shù)據(jù)向量;m為模型參數(shù)向量;F(m)為模型正演響應(yīng);W為歸一化矩陣;‖?m‖2為模型粗糙度,對(duì)二維構(gòu)造,同時(shí)考慮了模型的橫向和垂向光滑問(wèn)題,?=[?y,?z]T,其中?y和?z分別表示橫向和縱向相鄰的模型參量的粗糙度矩陣;X2為d與F(m)的擬合差,表示為:

(2)

在迭代尋優(yōu)中,需要在每一步使用局部線性化[2]。設(shè)當(dāng)前模型為mk,則mk+1的響應(yīng)為:

(3)

式中:Δm=mk+1-mk為模型的修改量;Jk=(?F/?m)k為模型mk處的雅可比矩陣。將(3)式代入(1)式對(duì)m求導(dǎo)置零,可得到以μ為變量的模型迭代表達(dá)式:

(4)

式中:dk=d-F(mk)+Jkmk。μ值通常使用掃描法或者一維搜索法求得,其取值應(yīng)使得殘差平方和‖Wd-WF[mk+1(μ)]‖2最小。

MT Occam反演是迭代算法,通過(guò)不斷修改給定的初始模型,使其逐步趨近真實(shí)模型。每次迭代通過(guò)綜合判斷擬合差和粗糙度來(lái)找出最優(yōu)解,滿足一定條件時(shí),輸出模型并結(jié)束計(jì)算。

2 并行性分析

如前所述,反演迭代計(jì)算之間嚴(yán)格順序依賴關(guān)系,無(wú)法直接并行實(shí)現(xiàn),因此并行計(jì)算只能在一個(gè)迭代步內(nèi)進(jìn)行。一個(gè)迭代步包含模型正演響應(yīng)F(m)、雅可比矩陣Jk、叉積矩陣、楚列斯基(Cholesky)分解、拉格朗日乘子μ等主要計(jì)算項(xiàng),不同計(jì)算項(xiàng)通常具有不同的并行特征,因而需要采用不同的并行策略。

2.1 粗粒度及細(xì)粒度并行性

MT Occam反演的一個(gè)迭代步內(nèi),含有粗粒度和細(xì)粒度并行性。非常耗時(shí)的雅可比矩陣計(jì)算和運(yùn)行次數(shù)極多的模型正演響應(yīng)計(jì)算都是針對(duì)頻率分別進(jìn)行的,各頻率計(jì)算之間無(wú)數(shù)據(jù)耦合,容易基于頻率劃分實(shí)現(xiàn)大粒度并行。當(dāng)采用掃描方式求取拉格朗日乘子時(shí),μ值點(diǎn)預(yù)先按固定間隔選取,其計(jì)算相互獨(dú)立,也可以實(shí)現(xiàn)大粒度分解。μ值一維搜索方式下,采用黃金分割法尋找凹區(qū)間,新μ值點(diǎn)的確定依賴前點(diǎn),屬于迭代方式,無(wú)法直接并行,但可以在正演響應(yīng)計(jì)算部分以分頻方式實(shí)現(xiàn)并行計(jì)算。

細(xì)粒度并行性分布于反演過(guò)程中不同的計(jì)算項(xiàng)。例如,一個(gè)頻率的模型正演響應(yīng)計(jì)算,需要采用高斯消元法求解形如ax=b的大型線性方程組,其中a為有限單元?jiǎng)偠汝?b為常數(shù)項(xiàng),x為待求向量,消元過(guò)程中的行操作可以分解,形成極細(xì)粒度計(jì)算任務(wù),這些計(jì)算可以通過(guò)解耦實(shí)現(xiàn)并行處理。雅可比矩陣計(jì)算也具有同樣的特征。此外,每次迭代中,需要完成多個(gè)矩陣叉乘運(yùn)算和多次Chelosky分解,這些部分也都具有極細(xì)粒度的并行性。

2.2 傳統(tǒng)并行方式的性能瓶頸及解決方案

基于頻率和μ值掃描點(diǎn)的分解具有天然的并行性,粒度也可以通過(guò)任務(wù)組合加以調(diào)整,適合消息傳遞模式下的大粒度并行,是電磁法正反演傳統(tǒng)并行計(jì)算的主要分解方式,文獻(xiàn)[2]、文獻(xiàn)[6]和文獻(xiàn)[10]都采用了頻率劃分策略,文獻(xiàn)[8]則結(jié)合采用了頻率劃分和μ值掃描點(diǎn)劃分。但一般來(lái)說(shuō),電磁法勘探中的觀測(cè)頻率數(shù)有限,通常不會(huì)超過(guò)40個(gè),反演中μ乘子掃描點(diǎn)數(shù)一般也不會(huì)超過(guò)30個(gè),因此,基于頻率和μ值掃描點(diǎn)的大粒度分解方式并行度是受限的[14]。以f表示程序中的串行部分,T1為串行執(zhí)行時(shí)間,假定并行工作由N份大小固定的工作tu組成,每份工作可以與其它工作并行執(zhí)行,但不再繼續(xù)細(xì)分(一個(gè)包含N個(gè)觀測(cè)頻率的模型正演及偏導(dǎo)數(shù)計(jì)算就具有此種特征),則在一個(gè)含有p個(gè)處理機(jī)的并行機(jī)上執(zhí)行程序的時(shí)間可以表示為:

(5)

公式(5)中N實(shí)際為并行度,當(dāng)P超過(guò)這一并行度時(shí),并行運(yùn)行時(shí)間便不再減少。這暴露出MT Occam反演傳統(tǒng)并行方式的性能瓶頸:

1) 當(dāng)并行度較低時(shí),增加處理機(jī)數(shù)沒(méi)有意義。

2) 由于并行度并非遠(yuǎn)大于處理機(jī)數(shù),很難保證所有處理機(jī)能幾乎平均地分配計(jì)算任務(wù),這容易導(dǎo)致短板效應(yīng),降低效率。

由此可見(jiàn),在消息傳遞模式下,MT Occam反演的并行計(jì)算性能由于并行度受限而難以提高。

要繼續(xù)減小并行計(jì)算時(shí)間Tp,必須細(xì)分tu,即充分利用算法中的細(xì)粒度并行性。但不同并行特征的計(jì)算只有適應(yīng)計(jì)算機(jī)集群的體系結(jié)構(gòu)才能獲得最佳的計(jì)算效率[15]。在消息傳遞模式下,tu的細(xì)分沒(méi)有意義,因?yàn)橛?jì)算每個(gè)細(xì)粒度任務(wù)需要進(jìn)行繁雜的通信,每一次的通信都包含啟動(dòng)時(shí)間和傳輸時(shí)間,因此計(jì)算/通信比會(huì)很小,程序的運(yùn)行效率非常低下,極端情況下甚至?xí)?lái)災(zāi)難性后果。因此,細(xì)粒度并行不適合消息傳遞模式,需要運(yùn)行在共享內(nèi)存架構(gòu)上,這需要對(duì)傳統(tǒng)結(jié)構(gòu)進(jìn)行擴(kuò)展。

由于對(duì)細(xì)粒度并行的良好支持,GPU的通用計(jì)算技術(shù)(GPGPU)發(fā)展極為迅速并得到廣泛應(yīng)用[16-17]。CPU和GPU能夠容易地構(gòu)成混合異構(gòu)系統(tǒng),利用傳統(tǒng)消息傳遞模型、OpenMP編譯制導(dǎo)和最新GPU編程接口CUDA可以實(shí)現(xiàn)混合并行算法,從而靈活處理具有不同并行特征的計(jì)算問(wèn)題。由于MT Occam反演中存在的層次并行特征,本文采用異構(gòu)混合并行方案來(lái)解決反演中不同粒度的并行計(jì)算,利用GPU與CPU的協(xié)同工作,實(shí)現(xiàn)系統(tǒng)整體計(jì)算能力的最大化。

3 混合并行方案設(shè)計(jì)

圖1給出了混合并行模型框架,硬件自上向下分別是集群、節(jié)點(diǎn)和GPU計(jì)算設(shè)備,對(duì)應(yīng)的軟件環(huán)境為MPI,OpenMP和CUDA。

上層分布式并行模型MPI用于實(shí)現(xiàn)計(jì)算的節(jié)點(diǎn)間并行,本文采用全局靜態(tài)通信,在映射過(guò)程中由主節(jié)點(diǎn)進(jìn)程(主進(jìn)程)負(fù)責(zé)對(duì)子節(jié)點(diǎn)進(jìn)程(子進(jìn)程)進(jìn)行主動(dòng)發(fā)送和接收,子進(jìn)程被動(dòng)接收和發(fā)送,在確認(rèn)主進(jìn)程接收前持續(xù)阻塞狀態(tài)。圖2給出了MPI層主進(jìn)程的通信過(guò)程的實(shí)現(xiàn)代碼,其中nproc為子進(jìn)程數(shù)。

OpenMP層對(duì)上層發(fā)送的任務(wù)繼續(xù)細(xì)分,并對(duì)GPU進(jìn)行管理。反演的主要計(jì)算由GPU完成,節(jié)點(diǎn)CPU只負(fù)責(zé)少量的串行計(jì)算,因此沒(méi)有設(shè)計(jì)利用CPU多核的算法來(lái)完成GPU之外的工作,加入OpenMP層的主要目的是為將來(lái)預(yù)留GPU擴(kuò)展機(jī)制。

CUDA層負(fù)責(zé)完成有限單元法計(jì)算、矩陣叉乘、Chelosky分解等具有極細(xì)粒度重復(fù)計(jì)算特征的部分。主機(jī)端需要將計(jì)算數(shù)據(jù)拷貝至設(shè)備端顯存,調(diào)用Kernel函數(shù)通過(guò)多線程方式完成并行處理,CPU最后將結(jié)果拷貝回內(nèi)存。

圖3為MT Occam CPU/GPU并行反演簡(jiǎn)化流程:主進(jìn)程讀入數(shù)據(jù)、對(duì)頻率分組并分發(fā)數(shù)據(jù)到各子進(jìn)程;子進(jìn)程接收任務(wù),調(diào)用GPU完成模型正演和雅可比矩陣計(jì)算;主進(jìn)程接收子進(jìn)程計(jì)算結(jié)果并對(duì)模型正演響應(yīng)和雅可比矩陣進(jìn)行匯總,調(diào)用GPU矩陣算法計(jì)算叉積矩陣,針對(duì)不同的μ值,調(diào)用GPU Cholesky分解算法計(jì)算相應(yīng)模型并發(fā)送給子進(jìn)程;子進(jìn)程接收模型數(shù)據(jù),調(diào)用GPU計(jì)算正演響應(yīng);主進(jìn)程接收模型正演響應(yīng)數(shù)據(jù),計(jì)算其擬合差和粗糙度,評(píng)判模型,確定是進(jìn)入下次迭代還是輸出最終模型。

圖1 CPU/GPU混合并行模型框架

圖2 MPI層主進(jìn)程通信代碼

計(jì)算流程中正演響應(yīng)和雅可比矩陣?yán)糜邢迒卧ㄓ?jì)算,其線性方程組通過(guò)高斯消元法求解,此部分計(jì)算占比最大,是提高計(jì)算效率的核心部分。本文針對(duì)其系數(shù)矩陣的二維等帶寬壓縮存儲(chǔ)結(jié)構(gòu)和MT Occam反演的特點(diǎn),設(shè)計(jì)了高斯消元GPU算法。算法以一個(gè)線性方程組的消元過(guò)程(一個(gè)工作三角形)為計(jì)算任務(wù)映射到CUDA模型的BLOCK線程塊,然后由CUDA二次映射到流處理器上完成,總的映射任務(wù)數(shù)與節(jié)點(diǎn)計(jì)算頻率數(shù)和GPU個(gè)數(shù)相關(guān)[18]。

叉積矩陣的求取和Cholesky分解計(jì)算量相對(duì)較小,不再分發(fā)到計(jì)算節(jié)點(diǎn),而由主進(jìn)程調(diào)用GPU單獨(dú)完成?？紤]到NVIDIA公司提供的函數(shù)庫(kù)已經(jīng)包含了部分高效算法,本文直接調(diào)用子程序cublasSsyrk_v2計(jì)算叉積矩陣,調(diào)用設(shè)備端函數(shù)cula_device_SPOSV實(shí)現(xiàn)Cholesky分解。

圖3 MT Occam CPU/GPU并行反演簡(jiǎn)化流程

4 理論模型試算

本文采用的異構(gòu)集群配置4個(gè)節(jié)點(diǎn),每個(gè)節(jié)點(diǎn)配置1個(gè)物理4核i5-3470 CPU,8GB DDR3雙通道內(nèi)存,1個(gè)具有2GB顯存的GTX 680 GPU(1536個(gè)流處理單元),節(jié)點(diǎn)間利用千兆局域網(wǎng)絡(luò)連接。操作系統(tǒng)為CentOS 6.4,開(kāi)發(fā)工具采用CDK,包含支持MPI,OpenMP和CUDA的C/C++和Fortran編譯器,同時(shí)提供了性能分析工具和調(diào)試工具。

4.1 反演實(shí)驗(yàn)?zāi)Ｐ?/h3>

為了評(píng)估算法的實(shí)際運(yùn)行效率,本文實(shí)驗(yàn)過(guò)程中設(shè)計(jì)了5個(gè)模型,在所構(gòu)建的異構(gòu)集群上進(jìn)行試算。每個(gè)模型均有2種采樣規(guī)模(即類(lèi)型1和類(lèi)型2),限于文章篇幅,這里只展示2個(gè)典型模型及其反演結(jié)果,即模型1(類(lèi)型1)、模型1(類(lèi)型2)、模型2(類(lèi)型1)和模型2(類(lèi)型2)。2種采樣規(guī)模的4項(xiàng)主要參數(shù)見(jiàn)表1,表中Nr為測(cè)點(diǎn)數(shù),Nm為模型參量,Ne為有限單元網(wǎng)格數(shù),Nf為觀測(cè)頻率數(shù)。最大反演次數(shù)定義為20。

表1 理論模型參數(shù)

模型1(圖4)設(shè)計(jì)為1000Ω·m基底,中部為相同電阻率的地壘,截面尺寸5.0km×12.5km,上部為100Ω·m均勻?qū)?。模?(圖5)設(shè)計(jì)為100Ω·m均勻半空間有2個(gè)不同電阻率棱鏡體,上部為高阻層。

圖4 反演模型1

圖5 反演模型2

類(lèi)型1模型測(cè)點(diǎn)數(shù)為21,模型參量為889,有限單元網(wǎng)格數(shù)為103×35,觀測(cè)頻率數(shù)為20。類(lèi)型2模型加密了各項(xiàng)參數(shù),4項(xiàng)參數(shù)分別增加至41,1549,193×35和36。反演初始模型定義為100Ω·m均勻半空間,理論模型數(shù)據(jù)均加入5%的隨機(jī)噪聲。

4.2 算法測(cè)試結(jié)果

對(duì)類(lèi)型1模型,分別用串行和并行算法計(jì)算5個(gè)模型的反演。定義偏差值為本次迭代串行算法所計(jì)算出的模型擬合差和粗糙度與相對(duì)應(yīng)并行算法計(jì)算值之間差值的絕對(duì)值。圖6給出了模型2反演過(guò)程中擬合差和粗造度的偏差值走勢(shì)曲線,偏差值均不大,說(shuō)明并行計(jì)算過(guò)程相對(duì)穩(wěn)定。

表2給出了類(lèi)型1和類(lèi)型2模型串行算法和并行算法的迭代次數(shù)、計(jì)算時(shí)間和加速比的統(tǒng)計(jì)結(jié)果。結(jié)果表明,規(guī)模較大的模型(類(lèi)型2)計(jì)算時(shí)間明顯增多,但加速比也更大。個(gè)別模型出現(xiàn)正演次數(shù)的差異,分析原因是CPU與GPU處理器在浮點(diǎn)運(yùn)算精度上的細(xì)微差異所致,最終反演結(jié)果沒(méi)有受到影響。圖7至圖10展示了模型1和模型2兩種規(guī)模下并行反演的結(jié)果(黑色實(shí)線為理論模型結(jié)構(gòu)輪廓),對(duì)模型1來(lái)說(shuō),兩種規(guī)模反演結(jié)果無(wú)太大差別,但從模型2的反演結(jié)果看,類(lèi)型2與理論模型明顯更為接近。

本混合并行方案反演結(jié)果正確,與純消息傳遞方案相比,加速比有較大提升,說(shuō)明方案合理可行。不同模型加速比有一定差異,這取決于模型規(guī)模、復(fù)雜度和映射到不同層次計(jì)算量的變化?？傮w而言,計(jì)算密度大的模型(如類(lèi)型2)有著更大的加速比。在本文4節(jié)點(diǎn)環(huán)境下,類(lèi)型1模型平均加速比為6.11,最高加速比為7.11,類(lèi)型2模型平均加速比為16.46,最高加速比達(dá)23.49。

圖6 模型2(類(lèi)型1)誤差曲線

表2 理論模型反演統(tǒng)計(jì)表

圖7 模型1(類(lèi)型1)反演結(jié)果

圖8 模型1(類(lèi)型2)反演結(jié)果

圖9 模型2(類(lèi)型1)反演結(jié)果

圖10 模型2(類(lèi)型2)反演結(jié)果

5 方案適用性分析

影響本文方案性能的因素主要源于兩方面:①通信延遲。本文測(cè)試結(jié)果表明,混合模型各層均會(huì)造成一定的通信延遲。由于OpenMP層和CUDA層是共享內(nèi)存方式,其延遲要遠(yuǎn)小于消息傳遞方式,因此,MPI層通信是主要的瓶頸。消息傳遞模型要求計(jì)算任務(wù)粒度盡可能大以提高通信效率,當(dāng)模型較小特別是頻點(diǎn)較少時(shí)(如類(lèi)型1模型),本方案上層的任務(wù)粒度就相對(duì)較小,盡管在下層采用了GPU處理,整體計(jì)算性能仍會(huì)受到影響。②GPU計(jì)算并行度。方案中,反演過(guò)程的大量計(jì)算最終將映射到GPU的流處理器上完成,細(xì)粒度并行度便成為影響性能的因素。GPU計(jì)算要求一次將盡可能多的數(shù)據(jù)投入到運(yùn)算中,當(dāng)并行度不夠時(shí),流處理器(SP)資源不能得到充分利用,計(jì)算效率將會(huì)降低。由于上述2方面的影響,本文混合并行方案更適合于規(guī)模較大的模型(多測(cè)點(diǎn)、多頻率、大密度網(wǎng)格)反演。

此外,系統(tǒng)存儲(chǔ)性能制約了本文集群的適用范圍。由于MT Occam反演對(duì)存儲(chǔ)空間的需求較大,小規(guī)模集群通常適用于二維模型反演,測(cè)試結(jié)果表明本文設(shè)計(jì)的不同規(guī)模的多個(gè)二維模型反演完全可行。但就三維反演而言,由于有限單元節(jié)點(diǎn)數(shù)遠(yuǎn)遠(yuǎn)大于二維反演,其需要占用的存儲(chǔ)空間是小型集群存儲(chǔ)資源難以滿足的,因此,如果要用于小規(guī)模三維模型反演,節(jié)點(diǎn)內(nèi)存和GPU顯存都必須增加,并對(duì)存儲(chǔ)進(jìn)行優(yōu)化。

6 結(jié)論

MT Occam反演中傳統(tǒng)的并行方法不具擴(kuò)放性,計(jì)算性能的提升難以利用擴(kuò)大集群規(guī)模的方式來(lái)實(shí)現(xiàn)。采用CPU/GPU多層次混合并行架構(gòu),通過(guò)挖掘細(xì)粒度并行分量,可以很好地彌補(bǔ)這一不足。線性方程組Gauss消元求解是反演過(guò)程中最為核心的計(jì)算,在GPU上實(shí)現(xiàn)其細(xì)粒度并行,能使極為耗時(shí)的雅可比矩陣和模型正演響應(yīng)計(jì)算都得到加速,是提升反演速度的關(guān)鍵。矩陣叉乘、Cholesky分解也具有良好的細(xì)粒度并行性,并可直接利用CUDA并行函數(shù)庫(kù)高效計(jì)算,這進(jìn)一步提高了算法的綜合計(jì)算性能。

本方案在4節(jié)點(diǎn)的小型集群上對(duì)較大模型(類(lèi)型2)反演獲得了平均16.46的加速比,最高達(dá)23.49,相對(duì)于純消息傳遞方式有較大提高。較小模型(類(lèi)型1)的反演受MPI層計(jì)算粒度和GPU層并行度的影響,平均加速比為6.11,最高為7.11,效率相對(duì)要低,但仍高于純消息傳遞方式。混合并行模式下通信延遲是需要關(guān)注的問(wèn)題,在我們的實(shí)驗(yàn)中,三層模型通信耗時(shí)最高達(dá)總計(jì)算時(shí)間的30%左右,說(shuō)明數(shù)據(jù)通信已經(jīng)成為本方案效率的主要制約因素,需要進(jìn)行更多的優(yōu)化,特別需要考慮計(jì)算與通信的重疊(overlap)。

致謝:感謝MT Occam 反演作者Steven Constable先生提供的串行源代碼。

參 考 文 獻(xiàn)

[1] 李小康.基于MPI的頻率域航空電磁法有限元二維正演并行計(jì)算研究[D].北京:中國(guó)地質(zhì)大學(xué),2011

LI X K.A MPI based parallel calculation on two dimensional finite element modeling of AEM[D].Beijing:China University of Geosciences,2011

[2] 胡祥云,李焱,楊文采,等.大地電磁三維數(shù)據(jù)空間反演并行算法研究[J].地球物理學(xué)報(bào),2012,55(12):3969-3978

HU X Y,LI Y,YANG W C,et al.Three-dimensional magnetotelluric parallel inversion algorithm using data space method[J].Chinese Journal of Geophysics,2012,55(12):3969-3978

[3] 吳小平,徐果明.大地電磁數(shù)據(jù)的Occam反演改進(jìn)[J].地球物理學(xué)報(bào),1998,41(4):547-554

WU X P,XU G M.Improvement of Occam’s inversion for MT data[J].Chinese Journal of Geophysics,1998,41(4):547-554

[4] SIRIPUNVARAPORN W,EGBERT G.An efficient data-subspace inversion method for 2-D magnetotelluric data[J].Geophysics,2000,65(3):791-803

[5] 陳小斌.大地電磁正反演新算法研究及資料處理與解釋的可視化集成系統(tǒng)開(kāi)發(fā)[D].北京:中國(guó)地震局地質(zhì)研究所,2003

CHEN X B.New forward and inversion algorithms and a visual integrated system for MT data[D].Beijing:Institute of Geology,China Seismological Bureau,2003

[6] TAN H D,TONG T,LIN C H.The parallel 3D magnetotelluric forward modeling algorithm[J].Applied Geophysics,2006,3(4):197-202

[7] 劉羽,王家映,孟永良.基于PC機(jī)群的大地電磁Occam反演并行計(jì)算研究[J].石油物探,2006,45(3):311-315

LIU Y,WANG J Y,MENG Y L.PC cluster based magnetotelluric 2-D Occam’s inversion parallel calculation[J].Geophysical Prospecting for Petroleum,2006,45(3):311-315

[8] 劉羽.基于機(jī)群的二維大地電磁Occam反演的并行計(jì)算研究[D].武漢:中國(guó)地質(zhì)大學(xué),2006

LIU Y.PC-Cluster based parallel computation research on 2-D magnetotelluric occam inversion[D].Wuhan:China University of Geosciences,2006

[9] 陳露軍.基于MPI的三維井地電磁場(chǎng)并行計(jì)算研究與實(shí)現(xiàn)[D].成都:成都理工大學(xué),2006

CHEN L J.Research on parallel computation of 3-D borehole-surface EM based on MPI[D].Chengdu:Chengdu University of Technology,2006

[10] SIRIPUNVARAPORN W,EGBERT G.WSINV3D-MT:vertical magnetic field transfer function inversion and parallel implementation[J].Physics of the Earth and Planetary Interiors,2009,173(3/4):317-329

[11] 李焱,胡祥云,吳桂桔,等.基于MPI的二維大地電磁正演的并行計(jì)算[J].地震地質(zhì),2010,32(3):392-401

LI Y,HU X Y,WU G J,et al.Research of 1-D magnetotelluric parallel computation based on MPI[J].Seismology and Geology,2010,32(3):392-401

[12] CONSTABLE S C,PARKER R L,CONSTABLE C G.Occam’s inversion:a practical algorithm for generating smooth models from EM sounding data[J].Geophysics,1987,52(3):289-300

[13] GROOT-HEDLIN DE C,CONSTABLE S C.Occam’s inversion and the North American Central plains electrical anomaly[J].Journal of Geomagnetismc and Geoelectricity,1993,45(9):985-999

[14] 劉羽.MT Occam 并行反演方案及性能分析[J].武漢理工大學(xué)學(xué)報(bào),2007,29(12):136-140

LIU Y.Parallel MT occam inversion scheme and its performance analysis[J].Journal of Wuhan University of Technology,2007,29(12):136-140

[15] 馬召貴,趙改善,武港山,等.起伏地表疊前時(shí)間偏移的多級(jí)并行優(yōu)化技術(shù)[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

[16] JACOBSEN D A,SENOCAK I.Multi-level parallelism for incompressible flow computations on GPU clusters[J].Parallel Computing,2013,39(1):1-20

[17] 張猛,王華忠,任浩然,等.基于CPU/GPU異構(gòu)平臺(tái)的全波形反演及其實(shí)用化分析[J].石油物探,2014,53(4):461-467

ZHANG M,WANG H Z,REN H R,et al.Fall waveform inversion on the CPU/GPU heterogeneous platform and its application analysis[J].Geophysical Prospecting for Petroleum,2014,53(4):461-467

[18] XIAO Y,LIU Y.GPU acceleration for the gaussian elimination in magnetotelluric Occam inversion algorithm[J].Proceedings of the 4thInternational Conference on Computer Engineering and Networks,2015:123-131