随着数字计算机和超大规模集成电路的飞速发展，数字信号处理技术得到了越来越广泛的应用，快速傅里叶变换（fast Fourier transform，FFT）是其中最基本、使用最频繁的核心算法之一. 其他诸多信号处理算法，例如卷积、滤波、频谱分析等，都可以转化为FFT来实现.

超级计算机为解决当今海量级数据的科学与工程计算问题提供了良好的平台. “神威·太湖之光”是世界上首台峰值每秒浮点运算次数（FLOPs）超过100×10¹⁵的超级计算机，也是中国第一台全部采用自主技术构建的排名世界第一的超级计算机^[1]. “神威·太湖之光”由国产“申威26010”处理器组成，单处理器峰值FLOPs可达3.1×10¹². 相比之下，访存带宽仅有136.5 GB/s则略显逊色，计算访存比接近22.7. 在面对诸如FFT这样兼具计算密集型和访存密集型特点的程序时，需要设计特定的并行实现策略，以优化其性能.

本文的主要工作是根据“申威26010”处理器的架构特点和编程规范，提出针对大点数FFT的众核优化方案. 该方案源自经典的Cooley-Tukey FFT算法，将一维数据迭代分解为二维FFT进行加速. 为了解决“列FFT”的读写和计算问题，降低矩阵转置的影响，提出“列均分-行连续”结合DMA跨步传输的读写策略. 该方法还可以推广至其他需要进行大量列数据操作的应用. 计算部分采用旋转因子优化和向量化操作，传输部分采用双缓冲策略和寄存器通信，能够充分利用众核系统的计算资源和传输带宽，达到了良好的加速效果.

FFT是离散傅里叶变换（discrete Fourier transform，DFT）的快速算法，由Cooley等^[2]提出，经过近60年的发展，形成了庞大的算法体系，陆续衍生出包括高基^[3]、混合基^[4]、分裂基^[5]、多维度^[6]等多种FFT算法，以应对各种数据类型和数据量.

在具体实现方面，FFT的硬件优化主要针对FPGA、ASIC、DSP等专用处理器，大致可以分为基于存储器和基于流水线2类. 基于存储器的结构利用数据在计算单元和存储单元之间往复运算和循环路由的特征来降低硬件复杂度^[7-8]. 基于流水线的结构通过设计单路/多路的前馈/反馈路径，提升蝶形单元的吞吐率^[9-10].

FFT的软件优化主要是通用平台上高性能软件包（库）的开发. 常见的有CPU平台上Frigo等^[11]开发的FFTW、Intel公司^[12]的MKL FFTW等、GPU平台上NVIDIA公司^[13]的CUFFT、Ayala等^[14]提出的heFFTe等. 其中，FFTW是目前使用最广泛、综合性能较好的FFT库. 张明^[15]在龙芯3B处理器上、郭金鑫^[16]在ARM V8平台上、操庐宁^[17]在X86-64平台上均对FFTW库进行了移植和优化.

国产“申威26010”处理器由上海高性能集成电路设计中心通过自主技术研制，采用片上计算阵列集群和分布式共享存储相结合的异构众核体系架构，使用64位自主申威指令系统. 单芯片集成了4个核组共260个核心，每个核组包含1个控制核心（主核）和由64个运算核心（从核）组成的核心阵列. 如图1（a）所示为处理器的结构示意图，如图1（b）所示为从核与主存、局存、寄存器之间的简化通信模型. 在编程规范上，核组内部采用Athread多线程库或OpenACC并行编程，核组间使用MPI或OpenMP编程接口. 新一代“申威26010-Pro”处理器于2021年发布，单核组运算性能和访存带宽分别提升了3.1倍和1.8倍.

“神威·太湖之光”由中国国家并行计算机工程技术研究中心于2015年12月完成研制，投入运行，已成功解决众多的科学工程计算问题. 目前已升级为配备有“申威26010-Pro”处理器的新一代超算系统. 代表性应用包括Yang等^[18]的“千万核可扩展大气动力学全隐式模拟”、Fu等^[19]的“非线性地震模拟”、Liu等^[20]的“超大规模量子随机电路实时模拟”等，分别于2016年、2017年、2021年三度荣获“戈登贝尔奖”.

在“神威·太湖之光”上开发了高性能扩展数学库xMath并移植了FFTW库，其中面向科学计算BLAS库的优化已十分成熟^[21-22]，但是缺少针对FFT的研究. 赵玉文等^[23]在“申威26010”上实现了迭代Stockham FFT计算框架，与单主核FFTW相比，获得了平均44.5x、最高56.3x的加速比. Stockham FFT是Cooley-Tukey FFT的一种常见变化算法，主要特点是消除了输入输出的倒位排序，可以避免传统FFT算法中的显式转置，在大数据量下能够实现一定的性能提升，代价是失去了蝶形结构的对称性和原位计算的特征，迭代框架的实现相对比较复杂.

针对FFT优化的开发需求，根据“申威26010”处理器的架构特点与FFT算法的内在联系，提出“列均分-行连续”读写策略. 结合其他优化手段，对数据进行巧妙的划分、重排和交换，可以减少行列转置带来的影响，完成了经典Cooley-Tukey FFT的众核实现，达到了大点数下平均60x以上的加速比.

离散傅里叶变换DFT的计算公式为

(1) $ \begin{array}{*{20}{c}} {X\left( k \right) = \displaystyle\sum\limits_{n = 0}^{N - 1} {x\left( n \right)W_N^{nk}} };&{k = 0,1, \cdots ,N - 1} \end{array}. $

式中： $W_N^{nk} $为旋转因子， $W_N^{nk} = {\rm{exp}} ({ - {\rm{j}}{{2{\text{π}} nk} / N}})$，其中j为虚数单位；N为序列长度.

FFT作为DFT的快速算法，核心思想是将长输入逐次分解成短序列DFT进行递归计算，通过分而治之的策略，将计算复杂度降至 $\mathcal{O}\left( {N{{\log }_2}N} \right)$. 大点数FFT的加速算法沿用了上述思想，通过将一维数据转化为二维或更高维，减少了每一维上的数据量，更有利于进行多核/众核并行.

将输入规模为 $N = {N_1} \times {N_2}$的一维数组映射为 ${N_1}$行、 ${N_2}$列的二维矩阵. 令输入、输出序号为

(2) $ n = \left( {{n_1},{n_2}} \right) = {n_1} {N_2}+{n_2}, $

(3) $ k = \left( {{k_2},{k_1}} \right) = {k_2} {N_1}+{k_1}. $

其中， $ 0 \leqslant {n_1}、{k_1} \leqslant {N_1} - 1 $， $ 0 \leqslant {n_2}、{k_2} \leqslant {N_2} - 1 $. 有如下关系：

(4) $ \begin{split} & X\left( {{k_2}{N_1}+{k_1}} \right) = \displaystyle\sum\limits_{n = 0}^{N - 1} {\left[ {x\left( {{n_1}{N_2}+{n_2}} \right)} \right.} \cdot {{{\rm{exp}}\;{( - {\rm{j}}({{2{\text{π}} }}/{N})\left( {{n_1}{N_2}+{n_2}} \right)\left( {{k_2}{N_1}+{k_1}} \right))}}} ] = \\ & \displaystyle\sum\limits_{{n_2} = 0}^{{N_2} - 1} {\displaystyle\sum\limits_{{n_1} = 0}^{{N_1} - 1} {\left[ {x\left( {{n_1}{N_2}+{n_2}} \right)} \right.} } {{{\rm{exp}}\;{( - {\rm{j}}({{2{\text{π}} }}/{N})\left( {{n_1}{k_2}{N_1}{N_2}+{n_1}{k_1}{N_2}+{n_2}{k_2}{N_1}+{n_2}{k_1}} \right))}}} ] = \\ & \displaystyle\sum\limits_{{n_2} = 0}^{{N_2} - 1} {\left\{ {\left[ {\displaystyle\sum\limits_{{n_1} = 0}^{{N_1} - 1} {x\left( {{n_1}{N_2}+{n_2}} \right)} } \right.} \right.} \Biggl.{20} {\Biggl.{20} {{{\rm{exp}}\;{( - {\rm{j}}2{\text{π}} {{{n_1}{k_1}}}/{{{N_1}}})}}} \Biggl]{20} {{\rm{exp}}\;{( - {\rm{j}}2{\text{π}} {{{n_2}{k_1}}}/{N})}}} \Biggl\}{20} {{\rm{exp}}\;{( - {\rm{j}}2{\text{π}} {{{n_2}{k_2}}}/{{{N_2}}})}} \end{split} . $

一维转二维FFT的计算流程（见图2）如下.

1）将一维数据按行优先顺序排列为 ${N_1}$行、 ${N_2}$列的二维矩阵.

2）进行每一列 ${N_1}$点FFT，共进行 ${N_2}$次，按原位置存回.

3）对于矩阵中的每一个数据，根据其坐标，乘以旋转因子 ${{\rm{exp}}\;{( - {\rm{j}}{{2{\text{π}} {n_2}{k_1}} / N})}}$.

4）进行每一行 ${N_2}$点FFT，共进行 ${N_1}$次，按原位置存回.

5）按列优先顺序将数据读出，得到 $N = {N_2} {N_1}$的最终结果.

将大点数FFT在“申威26010”上并行实现时，主要面临如下难题. 1）主核与从核之间列数据的读写问题. 2）从核间行列数据的分布与转置问题. 3）从核内小点数FFT的高效实现问题.

下面给出具体的并行优化解决方案.

使用双精度浮点复数，包括8 B的实部和8 B的虚部. “申威26010”从核对局存LDM的访问延迟为4个CPU周期，仅为对主存访问延迟的1/40，应尽可能将主存中的数据放入LDM内再计算. 考虑到单从核LDM大小仅为64 kB，采用双缓冲策略需要2份输入/输出缓冲区（3.6节），将可单次写入LDM的最大数据量规定为512个复数. 输入和输出共占512×16×2×2=32 kB，余下一半空间存储中间结果和临时变量.

“申威26010”单个核组阵列包括8行8列共64个运算核心. 为了方便维度分解和FFT计算，总数据量采用2的正整数幂次，每一维上的最小尺寸设置为8. 下面给出基本的分解策略.

1）将总点数分解为 $N = {N_1} {N_2}$，通常 ${N_1}$和 ${N_2}$较大，需要继续分解. 以 ${N_1}$为例，可以按以下a）、b）、c）共3种情况进行分解.

a） ${N_1} < 64$，即8、16、32点，停止分解，直接在单核上进行小点数FFT计算.

b） $64 \leqslant {N_1} \leqslant 256$，将 ${N_1}$进行二维分解，包括8×8、16×8、32×8共3种方案. 前一位数（8、16、32）代表单核上的FFT点数，后一位数（8）代表启动一行（或一列）8个从核进行计算.

c） $512 \leqslant {N_1} \leqslant 2\;048$，将 ${N_1}$进行三维分解，包括8×8×8、16×8×8、32×8×8共3种方案. 前一位数（8、16、32）代表单核上的FFT点数，后两位数（8×8）代表启动64个从核进行计算.

2）同理， ${N_2}$按以上3条进行再分解. 当 ${N_1}、 {N_2} \geqslant 4\;096$时，对总点数 $N$进行三维及以上分解，并启动多个核组.

3）每多分解一次，需要将当前FFT替换为图2所示的“列FFT-乘旋转因子-行FFT”流程，直至迭代计算完毕.

下面以输入数据量为262144=512×512的实例，阐述各创新实现技术.

在行优先的存储系统中，对列数据的读写是离散非连续的，传输效率很低. 一维转二维FFT并行方案的主要难点是对矩阵“列FFT”的处理. 大数据量加剧了FFT访存密集型特点带来的上述影响，若在主核上先进行显式转置再传输，则会增加不可忽略的额外开销.

为了提升“列FFT”的并行效率，结合FFT算法、处理器架构以及数据在从核间的分布特点，提出“列均分-行连续”读写策略，即并非用单个从核去读取主存中的列数据，而是把每一列的数据循环平分给64个从核. 如图3所示，对于“大矩阵”的每一列，0号核读第0行，1号核读第1行，······，63号核读第63行. 然后0号核读第64行，······，直至63号核读第511行. 从列的方向来看，每一列循环平分到64个从核内. 每个核一次可以读取多列，从行的方向来看，数据是连续读取的，即为“列均分-行连续”策略.

每个从核分到512/64=8行的数据. 由于LDM最大单次输入量为512，可以连续读入“大矩阵”每一行中512/8=64个数. 由于“申威26010”主存与局存之间的DMA带宽性能在传输粒度≥256 B（32个double数据）时达到峰值，采用“列均分-行连续”策略的第1个优点，即在保证一列数据能够完整读入从核阵列的同时，符合行优先的存储结构，且充分利用了传输带宽，从而有效解决了列读写的不连续性.

“列均分-行连续”策略的第2个优点，即从核间的数据分布与大点数分解方案吻合，有利于“列FFT”的后续计算. 如图4所示为“大矩阵”列数据在从核间的分布，方格中的数字代表该数据所在的从核号. 将“大矩阵”任意一列512个数（见图4（a））排列成64×8的“中矩阵”（见图4（b）），再将“中矩阵”的任意一列继续转换成8×8的“小矩阵”（图4（c）以图4（b）中的第0列为例）. 此时，“小矩阵”中的每一列数据均保存在同一个从核中，因此可以直接计算“小矩阵”“列FFT”.

图4（c）中，“小矩阵”的每一行数据分散在同列8个从核中，可以利用“申威26010”的寄存器通信机制来实现从核间的数据交换. 使用putc和getc指令开启寄存器列通信模式，同列8个核之间彼此将各行数据传到对应的从核，并接收其他从核传给自己的数据，即可完成“小矩阵”的列与行转置（见图4（d）），继续执行“小矩阵”“行FFT”计算. 在“小矩阵”FFT计算完毕后，核间数据按图4（d）中的列优先顺序排列，实际上与图4（b）保持一致. “中矩阵”的转置过程与“小矩阵”类似，但须使用putr和getr指令开启寄存器行通信模式（见图4（e））. 在完成“中矩阵”FFT后，将数据按列优先顺序写回主存（见图4（f））. LDM内的64列数据均按图4（f）的顺序排列，因此写回过程符合“列均分-行连续”原则，DMA传输可达最大速度.

FFT计算采用单指令多数据流操作（single instruction multiple data，SIMD）进行加速. “申威26010”支持宽度为4个double数据的向量运算（doublev4），但须解决以下2个问题. 1）若将同组FFT数据放在一个向量内，则不同间距的蝶形会引起数据依赖问题（见图5（a））. 2）若分散在不同向量内，则离散数据装载/卸载向量的效率较低.

当前每个从核LDM内的数据排列为8×64的矩阵，行方向上为从“大矩阵”中读入的一行64个数据（见图3），列方向上是“小矩阵”中一列（行）（见图4（c）、（d））或“中矩阵”中一行（见图4（e））的8个数据. 采用doublev4数组进行装载，每个doublev4内的4个子元素互不相关，不同doublev4内相同位置上的子元素来自同一组FFT数据，如图5（b）所示，既消除了向量内数据依赖，可以一次完成4组FFT计算，又使得数据在LDM内的存储连续且边界对齐，向量的装载和卸载效率最高. 这是“列均分-行连续”策略的第3个优点.

根据欧拉公式可知，FFT的旋转因子 $W_N^{nk}$可以表示为 ${{\rm{exp}}\;{({{ - {\rm{j}}2{\text{π}} nk} / N})}} = \cos \left( {{{2{\text{π}} nk} / N}} \right) - {\rm{j}}\sin \left( {{{2{\text{π}} nk} / N}} \right)$. “申威26010”上调用三角函数的用时为浮点乘法的近百倍，消耗了近80%的运算时间. 优化旋转因子的计算过程如下.

1）和差化积. 所有旋转因子系数均为 ${{2{\text{π}}} / N}$的 $n k$倍（ $N$为矩阵规模， $n$和 $k$为坐标）. 对基本单元 $\cos \left( {{{2{\text{π}} } / N}} \right)$和 $\sin \left( {{{2{\text{π}} } / N}} \right)$使用和差化积，用复数乘法代替三角函数，大幅减少三角函数的调用次数.

2）数据复用. 每个从核读入的“大矩阵”多列数据，实际上每一列在图4所示的“小矩阵”和“中矩阵”中所处的位置都是一样的，即坐标相同，因此只须计算第1列旋转因子并复用.

3）提前计算. 所有的旋转因子均由矩阵尺寸，即输入数据量决定，可以提前完成计算. 和差化积缩短了计算时间，数据复用减少了计算次数，将旋转因子计算与第一次DMA传输同步进行，不会给存储空间和计算时间带来额外负担.

此外，所有旋转因子的计算和相乘采用SIMD操作，可以与FFT结合.

“申威26010”为异步DMA传输，当数据需要在主存与局存之间进行多次传输时，可以采用双缓冲策略. 如图6所示，数字表示循环轮次，奇数轮和偶数轮分别占用2份不同的输入/输出缓冲区. 可以看出，除了第1轮读入和最后1轮写回外，在进行中间轮次计算的同时，可以进行下一轮读入及上一轮写回，以计算时间隐藏通信开销.

实际上，DMA还包括启动athread_get（读入）和athread_put（写回）的函数调用开销. 该步骤为串行过程，无法被计算隐藏，当一次完整读写需要多次调用这2个函数时，耗时不可忽略. 从核在每一轮“列均分”中共分到8行数据，须调用get和put函数各16次（实部和虚部各1次），因此无法被双缓冲隐藏的DMA时长 $T = 16 {T_1}$（设 ${T_1}$为函数调用时间），导致加速提升不到10%.

此时，每个从核在主存上的访问空间呈现等数据量（行连续）、等间隔（列均分）的特点（见图4（a）、（f））. 通过配置athread_get(dma_mode mode, void *src, void *dest, int len, void *reply, char mask, int stride, int bsize)的最后2个参数“间距”和“单次传输数据量”，可以实现跨步传输. 跨步读写仅须调用get和put各2次，无法隐藏的调用时间 $2 {T_3} < < 16 {T_1}$（设 ${T_{\text{3}}}$为跨步函数调用时间），达到了非常明显的加速效果.

如图7所示为“大矩阵”“列FFT”的完整计算流程，包括各环节采用的优化方法. “大矩阵”“行FFT”的计算过程与“列FFT”类似，此处不再赘述.

测试了数据量为32 768~4 194 304共8组随机生成的双精度浮点复数. 众核并行程序启用“申威26010”1个核组共64个从核，运行时间从athread_spawn开始至athread_join结束，包括传输和计算全过程. 对比基准为单主核上运行FFTW 3.3.4库以最优plan计算相同输入的用时，仅统计fftw_execute计算时间. 所有实验均重复多次，取稳定值求平均作为测试结果.

如图8所示为采用旋转因子优化、SIMD、双缓冲以及改用跨步传输后再结合双缓冲等各方案的加速比S_P.

如表1所示为优化前、后单个从核内旋转因子的计算量. 其中，小点数FFT所需的旋转因子数量分别为4、8、16，可以一次完成计算，所以忽略不计，只统计一维转二维方案中“列FFT”转“行FFT”之间所须乘上的旋转因子. 3.5节优化策略的核心思想是以和差化积所需的复数乘法和加法（或4次实数乘法和2次实数加法）来代替耗时较长的三角函数调用. 以262 144=512×512为例，优化前需要40 960次三角函数，优化后仅为24次三角函数加上1 136次浮点运算. 且该加速效果随着计算规模的增大而提高.

实验结果表明，若未对旋转因子计算进行优化，则多核并行总耗时比单主核FFTW长，而通过3.5节所述的和差化积、数据复用、提前计算等优化策略，可以达到15倍左右的加速比. 因此，旋转因子优化可以认为是计算环节中必须且最主要的加速部分.

在加入SIMD向量化运算后，可以在旋转因子优化的基础上提升近2倍的加速效果. SIMD覆盖了计算环节的几乎所有部分，包括FFT以及旋转因子相乘. 若单独测量运算部分，则可以达到约3.5倍的加速，接近4倍峰值，证明利用3.4节的方案，降低了装载/卸载向量引起的额外开销.

从核间寄存器单点通信的标称延迟为10拍，但是在应用中通常会叠加put/get指令执行开销、显式sync同步开销、网络堵塞延迟、连续读写延迟等，效率明显降低，需要采用汇编级流水线优化进行手动提升. 在本实验中，每次寄存器通信共须传输512×7/8 = 448个数，数量固定，因此耗时基本稳定，优于采用DMA传回主存，且改用DMA传输后会导致计算时间无法覆盖传输时间，使双缓冲失效. 寄存器通信是“小/中矩阵”“列/行”转置的首选方案.

在“列均分-行连续”策略中，直接采用双缓冲，加速比最大提升不到10%. 与get和put调用次数有关，数据量越大，调用次数越多，双缓冲效果越差. 在改用跨步传输再加上双缓冲后，达到了非常明显的加速效果，大数据量时可以再提升1倍以上，与3.6节的传输过程分析基本吻合.

如图9所示为不同输入点数时的DMA平均传输带宽B_W及利用率U_BW统计. 可以看出，“行连续”策略能够保证连续传输字节数≥256 B，因此可以在大多数情况下保持对读写带宽的充分利用. 当仅在最大数据量分解方案为2 048列时，连续读写粒度减小为128 B，导致性能下降. 以实测读带宽峰值27.9 GB/s进行折算，带宽利用率最高可达89.6%，最低为61.6%，平均为79.8%.

如图10所示为“主-从核DMA通信”、“从核间寄存器通信”、“从核FFT计算”3个环节各自在总时长中的占比T_P，左边为普通传输，右边为跨步读写（图例仅表示各部分功能的相对占比关系，不表示运行时长的绝对数值）. 可以看出，跨步传输+双缓冲策略使得主/从核间DMA传输时长占比明显下降，不再是整体性能的瓶颈，各环节耗时更均衡.

表2给出FFTW串行程序和最终多核并行方案各自的运行拍数C_S、C_P及两者的加速比测试结果，可以得到以下结论.

1）该并行优化方案的核心是“列均分-行连续”策略，优点如下. a）列数据读写连续非离散，能够充分利用DMA带宽. b）数据在从核间的分布能够直接进行后续FFT. c）数据在从核内的排列可以实现高效的SIMD. d）有助于旋转因子的快速运算. 以上4点保证了除必需的主存和局存之间的数据读写以及算法中的矩阵“列/行转置”操作以外，计算过程中基本没有额外的数据搬移，最大程度上减少了数据移动的次数，有效消除了FFT访存密集型特点带来的影响.

2）“列均分-行连续”策略的不足，即需要多次调用读写函数而带来的开销，可以采用跨步传输+双缓冲进行弥补. 合理的跨步传输能够有效地减少传输次数，从而提高异步传输时间在DMA总耗时中的占比，再利用双缓冲加以隐藏，有利于“列均分”策略的实施.

3）随着输入数据量的成倍增长，多核并行总时长的增加倍数与计算数据量的增长倍数基本相当，DMA传输粒度保持在256 B及以上. 测试结果显示，浮点运算性能及主、从核间的传输带宽基本可以达到并维持在峰值.

4）单主核FFTW在数据量较小时的运算性能较好，在数据量大于262 144之后基本保持成倍增长，应该是FFTW在不同数据量下会选择不同的最优算法所致.

5）本实验最终的加速比可以达到平均48x、最高65x以上的出色效果，尤其是在大点数情况下，能够始终保持在50x以上，平均接近60x.

6）利用“列均分-行连续+跨步传输”的策略，能够有效解决主、从核间“列”数据的读写问题，在矩阵应用场景下可以消除显式转置带来的影响，起到良好的加速效果. 在新一代“申威26010-Pro”处理器架构中，LDM扩大至256 kB，单次可传输数据量增加，更有利于DMA带宽的充分利用. 从核间的寄存器传输改为RMA通信机制，性能得到进一步的提升. 本方法重点关注处理器架构与算法的内在适配性，无须对软件算法框架进行大幅调整，因此可以作为一种通用优化策略进行推广使用，为“申威26010”高性能数学库中针对FFT算法的并行优化做出一定的贡献.

本文介绍了在国产“申威26010”众核处理器上对大点数FFT进行并行加速的方案. 为了有效地解决二维FFT对列数据的读写和计算问题，消除矩阵转置带来的影响，特别提出“列均分-行连续”的读写策略，结合SIMD向量操作、寄存器通信转置、旋转因子优化、双缓冲+跨步传输策略等优化方法，实现了经典Cooley-Tukey FFT算法在“申威26010”上的众核优化. 与单主核FFTW相比，可以达到平均48x以上、峰值65x以上的加速比.