譜方法求解水聲傳播問題的優(yōu)化與并行*

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

(國防科技大學氣象海洋學院，湖南 長沙 410073)

1 引言

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

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

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

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

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

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

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

考慮一個關(guān)于深度z和水平距離r的二維水聲傳播問題，設(shè)聲源角頻率為ω、深度為zs，聲壓為p(r,z)，密度為ρ，聲速為c，聲音傳播介質(zhì)(海水及海底沉積底)是水平分層的，則在聲源點以外的空間，聲壓滿足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譜方法計算簡正波模型和其開發(fā)的NM-CT程序(用Fortran語言開發(fā)的開源代碼，可從https://oalib-acoustics.org/Modes/index.html下載)為例，簡要介紹譜方法求解水聲傳播模型的流程。

(4)

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

(5)

其中,對定義在[-1，1]上的任意光滑函數(shù)ψ(x)使用Chebyshev變換，即利用Chebyshev多項式Tmp(x)展開(mp表示譜方法的截斷階數(shù))并對無窮項作有限的N階截斷近似處理，如式(6)所示：

(6)

(7)

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

(8)

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

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

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

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

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

步驟5計算聲壓p;

步驟6計算TL。

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

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

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

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

3.1 串行優(yōu)化

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

3.1.1 編譯器優(yōu)化

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

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

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

3.1.3 訪存優(yōu)化

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

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

簡化前的函數(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

簡化后的函數(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

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

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

改進前的計算:

1doi=1,nmodes

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

3enddo

4psi=matmul(psi,psizs);

改進后的計算:

1doi=1,nmodes

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

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

4enddo

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

3.1.4 精簡計算

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

3.2 并行優(yōu)化

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

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

在循環(huán)計算中，若各個計算之間相互獨立，則采用任務(wù)并行的策略，直接利用OpenMP多線程并行，還可以直接利用collapse對兩重循環(huán)進行展開,從而增大了并行度，并且可以保證線程間的負載均衡；在調(diào)度的過程中，嘗試靜態(tài)和動態(tài)等多種不同的調(diào)度方式，尋找最優(yōu)的調(diào)度方式。以Syn函數(shù)中的計算為例，嵌套循環(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)的計算沒有依賴關(guān)系，直接利用OpenMP多線程并行，用collapse對兩重循環(huán)進行展開，并且嘗試多種不同的調(diào)度方式，尋找效率最高的計算方式。

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

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

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

4.1 高性能平臺簡介

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

4.2 算例

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

Munk剖面是一個理想化的聲速剖面，它可以刻畫深海聲場傳播的規(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展示了利用譜方法畫出的傳播損失場。圖1b展示了海深1 000 m處利用不同方法計算的TL值，此算例不存在解析解，因此用傳統(tǒng)有限差分方法的程序代碼Kraken的計算結(jié)果作為參考解。從圖1中可以看出，截斷階數(shù)N值越大，計算結(jié)果越精確?？紤]實際復(fù)雜海洋環(huán)境對計算精度的要求以及算法整體計算量的大小，選擇N=1 000進行計算和后續(xù)的測試分析。

Figure 1 TL schematic diagram of Munk

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

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

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

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

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

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

4.3.3 訪存優(yōu)化

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

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

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

4.3.4 精簡計算

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

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

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)的并行效果。在眾核平臺上，若可以使用的最大線程數(shù)為nt，程序代碼中可并行部分的執(zhí)行時間所占百分比為q，則理想加速比S=1/(1-q+q/nt)。

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

Table 2 Program running time under different thread numbers

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

圖3顯示了Syn函數(shù)中3個子計算步驟的加速比情況。由于各網(wǎng)格點之間的計算沒有依賴，因此多線程并行計算具有良好的加速效果。隨著線程數(shù)的增加，加速比基本呈線性增長。當使用24個線程時，計算R(r)的代碼段加速比可達17.59。圖3中虛線表示對應(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

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

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

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

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

Figure 5 Comparison of calculation results before and after optimization

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

4.6 實驗結(jié)果分析

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

Table 3 Changes in program performance under different optimization techniques

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

5 結(jié)束語

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