基于FFTW库分步傅里叶变换算法并行方案研究

刘 帅， 李 智， 王 晶

（1.装备学院 研究生管理大队，北京101416； 2.装备学院 重点实验室，北京101416； 3.91336部队）

在求解抛物型波动方程时，Hardin和Tappert于1973年提出了一种分步傅里叶变换算法［1］，该算法是一种采用FFT技术的频域算法，其计算不受波长限制，在求解大区域电波传播问题中得到了广泛的应用。该算法属于一种步进迭代算法，当计算区域采样点数或迭代次数增加时，计算速度便会下降。

FFTW是由麻省理工学院（MIT）的Frigo和Johnson基于C语言开发的程序库［2］，能够高效地进行任意维数的离散傅里叶变换（DFT）、离散正弦变换（DST）以及离散余弦变换（DCT），是目前公认的最有效的变换工具。它能够根据用户的底层硬件属性如缓存、内存、寄存器大小等，自动调整算法参数，以获得最优方案。FFTW2.1.5和最新的FFTW3.3Alpha测试版还能够支持共享存储多线程并行和分布式存储MPI（消息传递接口）并行［3］。

为解决适应波动方程计算量大、计算速度慢的问题，作者借助于FFTW库研究了SSFT算法的并行方案，测试结果表明，所设计方案是有效的。

1 SSFT算法介绍

为阐述方便，先对SSFT算法的计算过程做简要的介绍。无论是三维还是二维情况，采用Feit－Fleck近 似 的 抛 物 型 波 动 方 程［4－6］都 可 以 化为如下形式

这样抛物方程的解可以近似为

称B为折射算子，A为绕射算子。在SSFT算法中，折射部分的计算在空间域内进行，绕射部分，也就是exp（ΔxA）u部分的计算在变换域即频域内进行，最终结果为

在一个计算步长Δx内，计算被分解成4个步骤：

1）计算初始场分布的傅氏变换U0＝F（u）；

2）绕射步计算T1＝exp（ΔxA）F（u）；

3）对T1计 算 其 傅 氏 反 变 换T2＝F－1（exp（ΔxA）F（u））；

4）折射步计算u（x＋Δx）＝exp（ΔxB）·F－1（exp（ΔxA）F（u））。

本文提出的并行方案正是立足于上述计算步骤，通过将数组拆分成较小规模的子数组，分布到多个进程分别计算，实现循环级的并行粒度。对于第2）和第4）步骤，由于只涉及对应点之间简单的复数乘法运算，在各个进程内部的数组上便可完成，无需进程间的通信，故有着直观的独立性或并发性。而第1）、第3）步的傅氏变换，依据定义，计算一个点需要变换前数组的所有点参与，将数组分发到各个进程后，这种计算必然涉及通信，事实上，通信正是影响傅氏变换并行算法的关键所在，因而第1）、第3）步构成了整个SSFT算法的计算核心。

2 基于FFTW库的并行方案

在介绍并行方案之前先给出评价并行算法的一个关键指标，即并行加速比，该值越大表明并行算法越高效。

并行加速比［7－8］：设串行执行时间为Ts，使用P个进程（或处理器）并行执行的时间为Tp（P），则加速比为

本文测试环境是实验室配备的银河集群机的一个机组，包含10个节点。每个节点集成2个Intel新一代Xeon高性能64位微处理器，配置有24G内 存、146G硬 盘。

文中所提出的SSFT算法并行方案核心，就是利用支持分布式存储的FFTW并行库实现计算过程中傅氏变换的计算，具体有2种思路。

1）串－并－串模式：即只对第1）、第3）步傅氏变换部分采用并行计算，其余步骤仍然在单个节点计算，算法流程图如图1，其中黑色双向箭头表示进程之间的通信。

图1 串－并－串模式流程图

2）分布式模式：此种模式下，除了傅氏变换的步骤外，绕射项和折射项的计算过程也分布到各个进程内进行，即预先将绕射项、折射项的指数数组计算好并分发到各个进程。此模式由于省去了各进程与主进程间的通信，整个算法的效率要高于第一种模式，流程图如图2。

图2 分布式模式流程图

本文所提并行方案具备通用性，在任何集群环境下都可以适用。

3 测试结果

为了测试方案的效果，需要不同的测试条件，包括采样点数、启动的进程数以及循环的步数。在不影响评价算法性能的前提下，本文取定循环的步数为Ns＝1 000。所有计时工作都是基于墙上时钟进行的。表1、表2为二维情况测试结果，表3为三维情况测试结果，其中P栏表示启动的进程数，P＝1即表示串行执行。

表1 串－并－串模式并行方案测试结果

表2 分布式模式并行方案测试结果

表3 3DSSFT基于FFTW库并行方案测试结果

由表1中数据可知，当高度方向采样点数N一定时，随着启动的进程数的增加，加速比呈上升趋势；当进程数一定时，随着采样点数的增加，加速比也呈上升趋势，这符合预期的设想。

表2亦表明，随着问题规模（N）的扩大和进程数（P）的增加，加速比也随之增大。对比表2中的数据，同等测试条件下分布式模式的加速比要优于串－并－串模式，这得益于第二种模式省去了一些全局收集操作，节约了通信开支。

由表3可看出，并行加速比是进程数和问题规模的增函数，随着问题规模的增大和进程数的增加，并行加速比呈上升趋势，且由于三维情况涉及的计算量更大，采用并行计算后速度提升更明显。

4 结 论

抛物型波动方程在电波传播预测、电磁环境数据生成等领域有着重要的应用，而用于求解此类方程的SSFT算法逐渐成为研究的热点。为了使求解过程更加快速，本文将并行计算技术引入到SSFT算法中，基于FFTW库研究了2种并行方案：串－并－串模式和分布式模式。测试结果表明，本文方案能够有效实现SSFT求解过程的加速，为相关问题研究提供思路。下一步工作将继续完成多种条件下方案的测试，以期得到方案的性能曲线，找到最佳配置状态。

（

）

［1］HARDIN R H，TAPPERT F D.Applications of the splitstep Fourier method to the numerical solution of nonlinear and variable coefficient wave equations［J］.SIAM Review Chronicle，1973，15：423.

［2］FRIGO M，JOHNSON S.FFTW：An adaptive software architecture for the FFT［C］／／IEEE.Proceedings of the IEEE International Conference on Acoustics，Speech，and Signal Processing.Seattle，WA，USA：IEEE，1998：1381－1384.

［3］FRIGO M.FFTW home page［EB／OL］.［2010－10－20］.http：／／www.fftw.org.

［4］胡绘斌.预测复杂环境下电波传播特性的算法研究［D］.长沙：国防科技大学，2006：18－19.

［5］LEVY M F.Parabolic equation methods for electromagnetic wave propagation［M］.London：London IEE Press，2000：21－22.

［6］THIEM K B.A 3Dparabolic equation（PE）based technique for predicting propagation path loss in an urban area［D］.California：Naval Postgraduate School，2001：6－9.

［7］陈国良.并行算法的设计与分析［M］.北京：高等教育出版社，1994：16.

［8］李晓梅，吴建平.数值并行算法与软件［M］.北京：科学出版社，2007：22.