馬 現(xiàn)，王勇獻(xiàn)，朱小謙，屠厚旺，李 朋，顏愷壯

(國(guó)防科技大學(xué)氣象海洋學(xué)院，湖南 長(zhǎng)沙 410073)

1 引言

近年來，水聲技術(shù)被廣泛應(yīng)用于水下通訊[1]、海洋環(huán)境測(cè)量[2]和海底測(cè)繪[3]等各個(gè)方面。海洋中聲波的傳播滿足基本的波動(dòng)方程，但由于海洋環(huán)境的時(shí)空復(fù)雜性，聲場(chǎng)分布易受環(huán)境的影響，導(dǎo)致聲波在海水中的傳播異常復(fù)雜[4]。波動(dòng)方程是所有的聲傳播數(shù)學(xué)模型的理論基礎(chǔ)，至今已發(fā)展出了多種傳播模型[5]，目前利用數(shù)值模型進(jìn)行水聲傳播計(jì)算已經(jīng)成為最常用的研究手段之一。

經(jīng)典的水聲傳播計(jì)算模型包括簡(jiǎn)正波近似法、拋物方程近似法等。每一種模型均是對(duì)原始波動(dòng)方程在特定條件下的近似，需求解一組微分方程。在傳統(tǒng)水聲傳播數(shù)值模擬中，有限差分方法是最常用的離散方法之一。例如，Porter等人[6]在簡(jiǎn)正波近似法中使用有限差分離散的方法，開發(fā)了Kraken程序；Collins等人[7]在拋物方程近似法中使用有限差分離散的方法，擴(kuò)展了針對(duì)二維寬水平傳播角情況的RAM(Range dependent Acoustic Model)程序；Lee等人[8]利用有限差分法開發(fā)了三維拋物模型的FOR3D程序;石鈴林等人[9]進(jìn)一步對(duì)拋物方程模型和FOR3D程序的聲傳播規(guī)律進(jìn)行了研究。盡管有限差分離散在傳統(tǒng)水聲傳播計(jì)算中發(fā)揮了重要作用，但它仍然存在很多不足，如處理復(fù)雜的邊值問題時(shí)不夠靈活，構(gòu)造高精度的差分格式困難等。

在計(jì)算流體力學(xué)、地震波傳播等其它方面，譜方法由于具有精度高、收斂速度快等優(yōu)點(diǎn)[10-12]，利用其進(jìn)行數(shù)值離散也得到了學(xué)者們的青睞。20世紀(jì)80年代，韋達(dá)[13-18]對(duì)譜方法的理論進(jìn)行了系統(tǒng)研究，發(fā)現(xiàn)譜方法無(wú)窮階的收斂特性；譜方法在物理方面也得到了廣泛的應(yīng)用，如大氣環(huán)流[19]、數(shù)值渦流[20]等；近年來，一些學(xué)者嘗試將譜方法引入到水聲傳播的數(shù)值計(jì)算中，取得了較好的效果，Tu等人[21]最近提出了一種可以解決不連續(xù)分層問題的簡(jiǎn)正波-譜方法數(shù)值模擬新手段，用于求解經(jīng)典的簡(jiǎn)正波模型，可以處理聲速、密度和衰減剖面不連續(xù)的問題，并開發(fā)了相應(yīng)的NM-CT(Normal Mode model program based on the Chebyshev-Tau spectral method)程序，計(jì)算結(jié)果具有較高的精度。

大海域、高頻聲源等復(fù)雜場(chǎng)景中的高精度聲場(chǎng)模擬仍然面臨著計(jì)算量較大、模擬速度慢和實(shí)時(shí)性不足等問題，難以滿足水下實(shí)際應(yīng)用場(chǎng)景中聲場(chǎng)快速分析的需求。隨著高性能計(jì)算技術(shù)的迅速發(fā)展，利用高性能平臺(tái)研究水聲傳播優(yōu)化與并行算法，為解決這個(gè)問題提供了新的途徑。對(duì)于簡(jiǎn)正波模型，吉虹宇[22]基于并行應(yīng)用框架OpenFOAM實(shí)現(xiàn)了水聲傳播并行數(shù)值模擬；對(duì)于拋物方程模型，范培勤等人[23]實(shí)現(xiàn)了弱三維情況下FOR3D模型的并行計(jì)算，徐閩等人[24]在高性能平臺(tái)天河二號(hào)上實(shí)現(xiàn)了FOR3D模型的并行計(jì)算，王魯軍等人[25]使用共享存儲(chǔ)并行編程(OpenMP)方法在多核計(jì)算機(jī)上實(shí)現(xiàn)了RAM程序的并行，均取得了較好的加速效果；對(duì)于射線模型，Xiao等人[26]在高性能計(jì)算平臺(tái)上，利用OpenMP對(duì)三維楔形海底的水聲傳播模型進(jìn)行并行加速與優(yōu)化；Zhu等人[27]綜合利用串行優(yōu)化和并行加速的方法，在天河二號(hào)高性能計(jì)算平臺(tái)上對(duì)三維水聲傳播模型進(jìn)行了優(yōu)化和并行加速，同樣也取得了較好的效果。

本文擬針對(duì)譜方法求解大規(guī)模水聲傳播問題過程中，計(jì)算速度慢、實(shí)時(shí)性差等問題，面向主流眾核平臺(tái)開展程序的優(yōu)化與并行，加速數(shù)值模擬。

各個(gè)章節(jié)主要內(nèi)容如下：第1節(jié)主要介紹譜方法求解水聲傳播問題的現(xiàn)狀與存在的問題；第2節(jié)主要介紹譜方法求解水聲傳播問題的原理及其計(jì)算流程；第3節(jié)介紹譜方法求解水聲傳播問題的優(yōu)化與并行方法；第4節(jié)給出了在天河二號(hào)高性能計(jì)算平臺(tái)上的數(shù)值測(cè)試結(jié)果，并對(duì)優(yōu)化效果進(jìn)行詳細(xì)分析和評(píng)估；第5節(jié)對(duì)全文的工作進(jìn)行總結(jié)，展望后續(xù)研究。

2 譜方法求解水聲傳播問題的原理及計(jì)算流程

考慮一個(gè)關(guān)于深度z和水平距離r的二維水聲傳播問題，設(shè)聲源角頻率為ω、深度為zs，聲壓為p(r,z)，密度為ρ，聲速為c，聲音傳播介質(zhì)(海水及海底沉積底)是水平分層的，則在聲源點(diǎn)以外的空間，聲壓滿足Helmholtz方程，如式(1)所示：

(1)

采用分離變量法求解該方程，則有p(r,z)=ψ(z)R(r)，其中R(r)為漢克爾函數(shù)，ψ(z)滿足如式(2)所示的常微分方程：

(2)

(3)

其中i為虛數(shù)單位。利用譜方法求解式(2)的基本思路是把解函數(shù)在一簇光滑的基函數(shù)上作近似展開(即譜展開)，將原始物理空間中的問題轉(zhuǎn)化為譜空間中求解展開系數(shù)的問題。下面以Tu等人[21]提出的Chebyshev-Tau譜方法計(jì)算簡(jiǎn)正波模型和其開發(fā)的NM-CT程序(用Fortran語(yǔ)言開發(fā)的開源代碼，可從https://oalib-acoustics.org/Modes/index.html下載)為例，簡(jiǎn)要介紹譜方法求解水聲傳播模型的流程。

(4)

其中k表示波數(shù)。利用算子￡代替式(4)左側(cè)作用于ψ(x)的算子，則式(4)簡(jiǎn)化為式(5)：

(5)

其中,對(duì)定義在[-1，1]上的任意光滑函數(shù)ψ(x)使用Chebyshev變換，即利用Chebyshev多項(xiàng)式Tmp(x)展開(mp表示譜方法的截?cái)嚯A數(shù))并對(duì)無(wú)窮項(xiàng)作有限的N階截?cái)嘟铺幚?，如?6)所示：

(6)

(7)

利用Chebyshev變換將式(5)由原始的物理空間變換到譜空間，最終形成一個(gè)線性代數(shù)特征值系統(tǒng)，寫為矩陣形式如式(8)所示：

(8)

輸入:初始聲速、密度、聲源位置和頻率等物理量。

步驟1利用Chebyshev多項(xiàng)式展開ψ(x)，離散式(2);

步驟2求解式(8)特征值和特征函數(shù);

步驟3利用Chebyshev變換求得矩陣L的特征解;

步驟4計(jì)算漢克爾函數(shù)R(r);

步驟5計(jì)算聲壓p;

步驟6計(jì)算TL。

算法包含如下3個(gè)主要的計(jì)算函數(shù)，第1個(gè)計(jì)算函數(shù)即算法第1～2行，命名為EIGsolve。EIGsolve函數(shù)主要計(jì)算式(8)矩陣L的特征值。第2個(gè)計(jì)算函數(shù)即算法第3行，命名為Genemodes。Genemodes函數(shù)主要對(duì)矩陣L的特征值和特征函數(shù)進(jìn)行Chebyshev變換，得到每個(gè)簡(jiǎn)正波模態(tài)函數(shù)ψm,m=1,2,…。第3個(gè)計(jì)算函數(shù)即算法第4～6行，也是算法計(jì)算量最大的部分，命名為Syn。Syn函數(shù)主要計(jì)算聲壓p和傳播損失TL。

譜方法與有限差分方法相比，是一種高精度的計(jì)算方法，兩者的水聲計(jì)算精度的對(duì)比在文獻(xiàn)[21]中有詳細(xì)分析。以文獻(xiàn)[21]分析的存在解析解的算例1為例，與有限差分方法相比，譜方法在垂直方向離散點(diǎn)更少的情況下取得更高的精度，適合于對(duì)聲場(chǎng)計(jì)算精度需求較高時(shí)的場(chǎng)景。下面以NM-CT程序?yàn)槔?，分析譜方法求解水聲傳播問題的優(yōu)化方法和并行方案，詳細(xì)描述每種優(yōu)化方法的原理并對(duì)測(cè)試結(jié)果進(jìn)行分析。

3 程序的優(yōu)化方法與并行方案

為了使程序性能優(yōu)化與并行工作更有針對(duì)性，本文利用Vtune工具測(cè)試和分析原始NM-CT串行程序，各個(gè)函數(shù)的耗時(shí)比例為：EIGsolve約2.91%，Genemodes約0.23%，Syn約54.95%，顯然，程序的主要時(shí)間開銷集中在Syn函數(shù)。因此，后續(xù)著重針對(duì)Syn函數(shù)進(jìn)行優(yōu)化與并行。首先對(duì)原始串行程序進(jìn)行不同層次的優(yōu)化，主要包括：編譯器調(diào)整和優(yōu)化、調(diào)用高性能數(shù)學(xué)庫(kù)MKL、優(yōu)化訪存和精簡(jiǎn)計(jì)算等；其次針對(duì)高性能眾核計(jì)算平臺(tái)，對(duì)調(diào)優(yōu)后的代碼進(jìn)行多線程并行加速，以充分利用多核心的計(jì)算資源。

3.1 串行優(yōu)化

串行代碼性能優(yōu)化與并行優(yōu)化同等重要，且極有可能獲得大幅度的加速效果，因此首先對(duì)代碼進(jìn)行串行優(yōu)化是非常必要的。

3.1.1 編譯器優(yōu)化

由于編譯器首先會(huì)對(duì)代碼整體進(jìn)行優(yōu)化，首先測(cè)試GNU 的gfortran和Intel的ifort 2種編譯器對(duì)代碼性能的影響。在保證正確性的基礎(chǔ)上選用優(yōu)化級(jí)別更高的編譯器-O3、添加選項(xiàng)-ipo過程間優(yōu)化和-funroll-all-loops循環(huán)展開選項(xiàng)對(duì)代碼進(jìn)行優(yōu)化。

3.1.2 利用高性能數(shù)學(xué)庫(kù)MKL

較大規(guī)模的矩陣乘法計(jì)算耗時(shí)比較突出，利用Intel公司開發(fā)的計(jì)算速度更快的MKL庫(kù)函數(shù)GEMM來替換原始的矩陣乘計(jì)算，以提升運(yùn)算速度。

3.1.3 訪存優(yōu)化

在代碼的優(yōu)化過程中，訪存也是影響代碼運(yùn)行速度的一個(gè)重要因素，保證訪存的連續(xù)性提升Cache命中率并盡可能減少內(nèi)存的使用，對(duì)代碼的性能有較為明顯的提升。

在提升Cache命中率方面，在讀取數(shù)據(jù)的過程中，按照存放順序讀取，可以保證較高的緩存命中率。以NM-CT代碼計(jì)算矩陣L為例，EIGsolve函數(shù)的部分代碼如下所示：

簡(jiǎn)化前的函數(shù):

1dok=1,size(v)

2doi=1,size(v)

3doj=1,size(v)

4if((i-1+j-1)==(k-1))

5Co(k，i)=Co(k，i)+v(j)*0.5;

6if(abs(i-j)==(k-1))

7Co(k，i)=Co(k，i)+v(j)*0.5;

8enddo

9enddo

10enddo

簡(jiǎn)化后的函數(shù):

1n=size(v);

2doi=1,n

3dok=1,n

4j=k-i+1;

5if(1 ≤j.and.j≤n)then

6Co(k,i)=Co(k,i)+v(j);

7endif

8j=i-k+1

9if(j≤n.and.j≥1)then

10Co(k,i)=Co(k,i)+v(j);

11endif

12Co(k,i)=Co(k,i)* 0.5;

13enddo

14enddo

簡(jiǎn)化前的函數(shù)為原始循環(huán)方式，代碼第1、2行改進(jìn)后，按照簡(jiǎn)化后的函數(shù)中第2、3行所示的方式，大大提高了Cache命中率，對(duì)程序性能的提升有較為明顯的作用。

在減少內(nèi)存使用方面，訪存優(yōu)化更直接有效的方式是減少內(nèi)存的使用。通過分析可以刪減部分不必要的數(shù)組，減少內(nèi)存空間的使用。以NM-CT代碼計(jì)算psi數(shù)組為例，計(jì)算psi數(shù)組的具體過程如下所示：

改進(jìn)前的計(jì)算:

1doi=1,nmodes

2psizs(i,i)=psi(…,i)*(…)+psi(…,i)*(…);

3enddo

4psi=matmul(psi,psizs);

改進(jìn)后的計(jì)算:

1doi=1,nmodes

2psi(:,i)=psi(:,i)*[psi(…,i)*(…)+

3psi(…,i)*(…)]

4enddo

改進(jìn)前的計(jì)算中，對(duì)psi數(shù)組的更新僅僅是將其每一列元素乘一個(gè)相同的數(shù)，即psizs對(duì)應(yīng)的元素值，因此直接將數(shù)組psi的每一列元素乘原數(shù)組,如改進(jìn)后的計(jì)算。通過分析發(fā)現(xiàn)，改進(jìn)后的計(jì)算可以少開辟一個(gè)維度nmodes×nmodes的數(shù)組，減少了內(nèi)存的使用，加快了其運(yùn)行速度。

3.1.4 精簡(jiǎn)計(jì)算

程序中包含有大量的計(jì)算過程，若能夠有效去除冗余計(jì)算，整個(gè)程序會(huì)有較大幅度的性能提升。對(duì)NM-CT的簡(jiǎn)化計(jì)算主要是減少分支判斷語(yǔ)句。以EIGsolve函數(shù)的計(jì)算為例，利用譜方法計(jì)算矩陣L時(shí)，部分代碼如簡(jiǎn)化前的EIGsolve函數(shù)所示。三重嵌套循環(huán)最內(nèi)層第4、5行賦值語(yǔ)句所需j索引值可以用i和k表示，從而極大地減少生成指令中的循環(huán)分支語(yǔ)句數(shù)目，有利于串行代碼的指令級(jí)調(diào)度與優(yōu)化。

3.2 并行優(yōu)化

當(dāng)程序中含有大量的循環(huán)和數(shù)據(jù)計(jì)算時(shí)，采用并行計(jì)算是最有效的提速方法。以NM-CT中的Syn函數(shù)計(jì)算為例，其包含有大量的循環(huán)和數(shù)據(jù)計(jì)算，因此這個(gè)函數(shù)中的代碼是可以并行計(jì)算的熱點(diǎn)代碼段。

(1)OpenMP多線程并行方案設(shè)計(jì)。

在循環(huán)計(jì)算中，若各個(gè)計(jì)算之間相互獨(dú)立，則采用任務(wù)并行的策略，直接利用OpenMP多線程并行，還可以直接利用collapse對(duì)兩重循環(huán)進(jìn)行展開,從而增大了并行度，并且可以保證線程間的負(fù)載均衡；在調(diào)度的過程中，嘗試靜態(tài)和動(dòng)態(tài)等多種不同的調(diào)度方式，尋找最優(yōu)的調(diào)度方式。以Syn函數(shù)中的計(jì)算為例，嵌套循環(huán)的并行過程如下所示：

1 !$omp parallel default(none)&

2 !$omp shared(…)private(…)

3 !$omp do collapse(2)schedule(static,…)

4doi=1，nr

5dok=1，nmodes

6bessel0=r(i)*kr(k);

7 callZBESH(…);

8bessel(k，i)=cmplx(CYR，CYI)

9enddo

10enddo

11 !$omp end do

12 !$omp end parallel

上述第3行代碼，循環(huán)內(nèi)的計(jì)算沒有依賴關(guān)系，直接利用OpenMP多線程并行，用collapse對(duì)兩重循環(huán)進(jìn)行展開，并且嘗試多種不同的調(diào)度方式，尋找效率最高的計(jì)算方式。

(2)利用多線程高性能計(jì)算庫(kù)MKL。

與OpenMP多線程相適應(yīng)，通過調(diào)用多線程版本的高性能計(jì)算庫(kù)MKL，可實(shí)現(xiàn)矩陣相乘GEMM函數(shù)的多線程加速，進(jìn)一步提升計(jì)算速度。

4 實(shí)驗(yàn)與結(jié)果分析

4.1 高性能平臺(tái)簡(jiǎn)介

為了檢驗(yàn)本文所用的優(yōu)化方法的效果，在高性能平臺(tái)天河二號(hào)上進(jìn)行了測(cè)試。天河二號(hào)高性能計(jì)算平臺(tái)的CPU 為Intel(R)Xeon(R)CPU E5-2692，每個(gè)結(jié)點(diǎn)包含2個(gè)CPU，每個(gè)CPU包含12個(gè)核。在測(cè)試過程中，采用gcc/6.3.0版本和Intel-compilers/mkl-15 2種編譯器進(jìn)行測(cè)試。在實(shí)際測(cè)試過程中，對(duì)于每一次優(yōu)化都進(jìn)行5次測(cè)試，將時(shí)間最短的測(cè)試時(shí)間作為最優(yōu)時(shí)間。

4.2 算例

譜方法求解水聲傳播模型能夠?qū)Χ喾N復(fù)雜的算例和模型進(jìn)行數(shù)值模擬。選取具有代表性的Munk波導(dǎo)算例進(jìn)行測(cè)試(本文所有圖表的彩色版本可從https://gitee.com/maxian-nudt/nm-ct_-parallel/issues下載)。

Munk剖面是一個(gè)理想化的聲速剖面，它可以刻畫深海聲場(chǎng)傳播的規(guī)律與特征，聲源位于1 000 m處，頻率為50 Hz，水層密度為1.0 g/cm3，沉積層密度為1.5 g/cm3，水層衰減系數(shù)為0，沉積層衰減系數(shù)為0.01,水層深度為2 500 m，沉積層深度為5 000 m。Munk波導(dǎo)的聲速值為c(z)=1500.0[1.0+其中,圖1a展示了利用譜方法畫出的傳播損失場(chǎng)。圖1b展示了海深1 000 m處利用不同方法計(jì)算的TL值，此算例不存在解析解，因此用傳統(tǒng)有限差分方法的程序代碼Kraken的計(jì)算結(jié)果作為參考解。從圖1中可以看出，截?cái)嚯A數(shù)N值越大，計(jì)算結(jié)果越精確?？紤]實(shí)際復(fù)雜海洋環(huán)境對(duì)計(jì)算精度的要求以及算法整體計(jì)算量的大小，選擇N=1 000進(jìn)行計(jì)算和后續(xù)的測(cè)試分析。

Figure 1 TL schematic diagram of Munk

4.3 串行優(yōu)化結(jié)果與分析

4.3.1 編譯器和編譯選項(xiàng)優(yōu)化

表1為原始串行程序在天河二號(hào)高性能平臺(tái)上不同編譯器和編譯選項(xiàng)的測(cè)試結(jié)果：選用Intel編譯器并增加合適的編譯選項(xiàng)可使性能進(jìn)一步提升。以gfortran編譯器下的-O2選項(xiàng)為基準(zhǔn)，通過更換編譯器及編譯選項(xiàng)優(yōu)化，最終獲得了1.67倍的加速效果。

Table 1 Effects of different compilers and compilation options on running time of serial programs

4.3.2 利用高性能數(shù)學(xué)庫(kù)MKL

測(cè)試計(jì)算聲壓大規(guī)模矩陣乘法代碼段在使用MKL庫(kù)前后的運(yùn)行時(shí)間，時(shí)間由201.2 s減少至56.83 s，使用MKL庫(kù)后運(yùn)行時(shí)間縮短144.37 s，提速3.54倍。

4.3.3 訪存優(yōu)化

在提升Cache命中率方面，對(duì)程序進(jìn)行訪存優(yōu)化，首先對(duì)計(jì)算矩陣L的兩重嵌套循環(huán)交換循環(huán)嵌套的順序，使訪存滿足空間的連續(xù)性，該代碼段優(yōu)化后運(yùn)行時(shí)間由3.03 s減少至2.79 s。

在減少內(nèi)存使用方面，通過分析具體計(jì)算過程，省略psizs數(shù)組，將運(yùn)算簡(jiǎn)化，該代碼段運(yùn)行時(shí)間由0.58 s減少至0.003 s。

根據(jù)訪存連續(xù)性的原理，對(duì)整個(gè)程序的嵌套循環(huán)和數(shù)組進(jìn)行檢查，最大限度地保證訪存的連續(xù)性，優(yōu)化后進(jìn)行測(cè)試，運(yùn)行時(shí)間由304.77 s減少至288.13 s，運(yùn)行時(shí)間減少了16 s，加速了1.05倍。訪存優(yōu)化不僅提升了程序性能，且精簡(jiǎn)了代碼，有利于后續(xù)對(duì)代碼的維護(hù)。

4.3.4 精簡(jiǎn)計(jì)算

對(duì)計(jì)算矩陣L的部分代碼減少分支判斷語(yǔ)句后，計(jì)算量大大下降，測(cè)試優(yōu)化前后每次調(diào)用此函數(shù)所需要的時(shí)間由3.03 s減少至0.02 s。由于調(diào)用矩陣L次數(shù)較多，測(cè)試簡(jiǎn)化運(yùn)算后對(duì)整個(gè)程序性能的影響。測(cè)試結(jié)果表明，程序運(yùn)行時(shí)間由304.77 s減少至242.19 s，對(duì)代碼進(jìn)行簡(jiǎn)化后，程序運(yùn)行時(shí)間減少了62 s，加速了1.26倍。

圖2主要展示了串行優(yōu)化前后3個(gè)主要函數(shù)的運(yùn)行時(shí)間變化。串行優(yōu)化效果最為明顯的是Syn函數(shù)，Syn函數(shù)利用了MKL替換原始矩陣乘計(jì)算、訪存優(yōu)化和精簡(jiǎn)計(jì)算3種方法。優(yōu)化后Syn函數(shù)運(yùn)行時(shí)間由320.91 s減少至80.43 s，加速比達(dá)到3.26，串行加速方法的有效性得到了證實(shí)。

Figure 2 Time changes of the three main functions and the total time before and after serial optimization

4.4 并行優(yōu)化

并行加速效果遵循Amdahl加速比定律[27]，據(jù)此可提前預(yù)判最優(yōu)的并行效果。在眾核平臺(tái)上，若可以使用的最大線程數(shù)為nt，程序代碼中可并行部分的執(zhí)行時(shí)間所占百分比為q，則理想加速比S=1/(1-q+q/nt)。

本文對(duì)3.2節(jié)的并行方法進(jìn)行了測(cè)試。首先對(duì)熱點(diǎn)Syn函數(shù)內(nèi)部的3個(gè)主要步驟(計(jì)算R(r)、計(jì)算聲壓、計(jì)算TL)進(jìn)行多線程并行，測(cè)試不同線程數(shù)目下的運(yùn)行時(shí)間，結(jié)果如表2所示。

Table 2 Program running time under different thread numbers

通過多線程并行加速后程序的總運(yùn)行時(shí)間由串行最優(yōu)版本的93.94 s減少為24.38 s，熱點(diǎn)函數(shù)Syn的運(yùn)行時(shí)間由原來的80.43 s減少至9.84 s。當(dāng)線程數(shù)為24時(shí)，整個(gè)程序的理想加速比S=5.68，實(shí)際并行后的絕對(duì)加速比為3.85，未達(dá)到理想值。究其原因，可能是由于串行程序中部分可并行的代碼之間夾雜部分不可并行的計(jì)算，需要多次進(jìn)行開關(guān)線程的操作，導(dǎo)致程序運(yùn)行時(shí)間增加，對(duì)程序的優(yōu)化效果產(chǎn)生不良影響。

圖3顯示了Syn函數(shù)中3個(gè)子計(jì)算步驟的加速比情況。由于各網(wǎng)格點(diǎn)之間的計(jì)算沒有依賴，因此多線程并行計(jì)算具有良好的加速效果。隨著線程數(shù)的增加，加速比基本呈線性增長(zhǎng)。當(dāng)使用24個(gè)線程時(shí)，計(jì)算R(r)的代碼段加速比可達(dá)17.59。圖3中虛線表示對(duì)應(yīng)線程數(shù)目下的最優(yōu)加速比。

Figure 3 Speedup of R(r),TL, and p under different thread numbers, and the optimal speedup under the corresponding number of threads

此外，為了進(jìn)一步測(cè)試多線程并行的任務(wù)調(diào)度對(duì)性能的影響,對(duì)Syn函數(shù)中的計(jì)算TL步驟進(jìn)行了測(cè)試，圖4展示了使用動(dòng)態(tài)調(diào)度以及不同粒度參數(shù)下的靜態(tài)調(diào)度時(shí)代碼的執(zhí)行時(shí)間。結(jié)果表明靜態(tài)調(diào)度隨著塊中的迭代次數(shù)chunksize值的增加，耗時(shí)逐漸減少，當(dāng)chunksize值達(dá)到1 024時(shí)，該代碼段的耗時(shí)達(dá)到最低。

Figure 4 Different scheduling methods to calculate the time change of TL under 24 threads

4.5 優(yōu)化前后結(jié)果對(duì)比

對(duì)程序進(jìn)行串行和并行優(yōu)化后，需要對(duì)程序改進(jìn)前后的計(jì)算結(jié)果進(jìn)行對(duì)比分析，比較計(jì)算結(jié)果精確度的變化，取聲源深度1 000 m的計(jì)算結(jié)果進(jìn)行比較，如圖5所示。

Figure 5 Comparison of calculation results before and after optimization

根據(jù)圖5優(yōu)化前后的計(jì)算結(jié)果，可以看出2條曲線基本重合，程序改進(jìn)前后計(jì)算結(jié)果非常接近，表明此次優(yōu)化沒有改變程序計(jì)算結(jié)果的精度。

4.6 實(shí)驗(yàn)結(jié)果分析

綜上所述，通過綜合使用串行優(yōu)化和多線程并行加速技術(shù)，分層海洋聲傳播Chebyshev-Tau譜方法程序NM-CT的性能得到了提升，計(jì)算效率大幅提升。以gfortran編譯器作為基準(zhǔn)版本，程序原始版本的運(yùn)行時(shí)間為584 s，優(yōu)化與并行的各階段的加速效果如表3所示。

Table 3 Changes in program performance under different optimization techniques

通過結(jié)果對(duì)比發(fā)現(xiàn)，首先選擇合適的編譯器對(duì)程序性能的提升有較大的作用，對(duì)于NM-CT程序，在天河二號(hào)的平臺(tái)上 Intel 的 ifort 編譯器對(duì)提升此程序的性能作用更加顯著。此外，加入合適的附加編譯選項(xiàng)、使用高性能MKL 函數(shù)庫(kù)、優(yōu)化訪存和合理簡(jiǎn)化計(jì)算等串行優(yōu)化方法，在單結(jié)點(diǎn)單線程下，充分利用了資源，取得了較好的加速效率。利用并行加速方法，并在此基礎(chǔ)上研究不同調(diào)度策略以及chunksize大小對(duì)程序的影響，發(fā)揮眾核潛能，進(jìn)一步提升了程序的運(yùn)行速度。針對(duì)天河二號(hào)單結(jié)點(diǎn)，本文提出的關(guān)于 NM-CT 程序的一系列優(yōu)化和并行加速方法非常有效，程序運(yùn)行時(shí)間從原始的584 s減少到最佳優(yōu)化版本的24 s，加速了23.98倍，優(yōu)化效果明顯，極大地提升了程序的性能，對(duì)解決計(jì)算大范圍海域聲傳播實(shí)時(shí)性問題做出了重要的貢獻(xiàn)。

5 結(jié)束語(yǔ)

為解決譜方法數(shù)值求解水聲傳播問題計(jì)算量大、實(shí)時(shí)性差等問題，本文充分利用高性能計(jì)算平臺(tái)對(duì)譜方法求解水聲程序進(jìn)行性能優(yōu)化與并行加速，首先從選擇合適的編譯器和優(yōu)化選項(xiàng)、調(diào)用高性能數(shù)學(xué)庫(kù)MKL函數(shù)、訪存優(yōu)化和精簡(jiǎn)計(jì)算4個(gè)方面對(duì)串行程序內(nèi)部進(jìn)行優(yōu)化;接下來在天河二號(hào)眾核平臺(tái)上對(duì)程序進(jìn)行更細(xì)粒度、更輕量級(jí)的多線程并行，在此基礎(chǔ)上研究不同調(diào)度策略和chunksize大小對(duì)程序性能的影響。以NM-CT程序?yàn)槔郎y(cè)試優(yōu)化與并行手段的加速效果，結(jié)果表明，本文所使用的優(yōu)化與并行方法，對(duì)于加速譜方法數(shù)值求解水聲傳播問題有非常明顯的加速效果，計(jì)算深海波導(dǎo)程序運(yùn)行時(shí)間由原始串行版本的584 s減少到24 s，加速了23.98倍，顯著提高了程序的運(yùn)行效率。進(jìn)一步分析程序，若截?cái)嚯A數(shù)更高，導(dǎo)致EIGsolve函數(shù)耗時(shí)較長(zhǎng)，因此在下一步的工作中，需要進(jìn)一步對(duì)程序在截?cái)嚯A數(shù)更高的情況下進(jìn)行優(yōu)化。盡管對(duì)于實(shí)際的大范圍海域聲場(chǎng)計(jì)算仍沒有達(dá)到實(shí)時(shí)處理的要求，但本文的工作為達(dá)成這一目標(biāo)邁進(jìn)了一大步。