中国科学技术大学计算机科学与技术系 国家高性能计算中心(合肥) 2004年12月 并 行 计 算 中国科学技术大学计算机科学与技术系 国家高性能计算中心(合肥) 2004年12月

第三篇 并行数值算法 第八章 基本通讯操作 第九章 稠密矩阵运算 第十章 线性方程组的求解 第十一章 快速傅里叶变换 第三篇 并行数值算法 第八章 基本通讯操作 第九章 稠密矩阵运算 第十章 线性方程组的求解 第十一章 快速傅里叶变换

第十一章 快速傅里叶变换 11.1 快速傅里叶变换 11.2 并行FFT算法

11. 1 快速傅里叶变换(FFT) 11. 1. 1 离散傅里叶变换(DFT) 11. 1. 2 DFT的顺序代码 11. 1 11.1 快速傅里叶变换(FFT) 11.1.1 离散傅里叶变换(DFT) 11.1.2 DFT的顺序代码 11.1.3 串行FFT递归算法 11.1.4 串行FFT非递归算法

离散傅里叶变换(DFT) 定义 给定向量A=(a0,a1,…,an-1)T,DFT将A变换为B=(b0,b1,…,bn-1)T 国家高性能计算中心（合肥） 2019/1/1

11. 1 快速傅里叶变换(FFT) 11. 1. 1 离散傅里叶变换(DFT) 11. 1. 2 DFT的顺序代码 11. 1 11.1 快速傅里叶变换(FFT) 11.1.1 离散傅里叶变换(DFT) 11.1.2 DFT的顺序代码 11.1.3 串行FFT递归算法 11.1.4 串行FFT非递归算法

DFT的顺序代码 代码1 代码2 注：代码1需要计算ωk*j b[j]=0 for k=0 to n-1 do for j=0 to n-1 do b[j]=0 for k=0 to n-1 do b[j]=b[j]+ωk*ja[k] end for 注：代码1需要计算ωk*j 代码2的复杂度为O(n2) 代码2 w=ω0 for j=0 to n-1 do b[j]=0, s=ω0 for k=0 to n-1 do b[j]=b[j]+s*a[k] s=s*w end for w=w*ω 国家高性能计算中心（合肥） 2019/1/1

11. 1 快速傅里叶变换(FFT) 11. 1. 1 离散傅里叶变换(DFT) 11. 1. 2 DFT的顺序代码 11. 1 11.1 快速傅里叶变换(FFT) 11.1.1 离散傅里叶变换(DFT) 11.1.2 DFT的顺序代码 11.1.3 串行FFT递归算法 11.1.4 串行FFT非递归算法

串行FFT递归算法 蝶式递归计算原理 令 为n/2次单位元根，则有 . 将b向量的偶数项 和奇数项 分别记为 和 注意推导中反复使用 国家高性能计算中心（合肥） 2019/1/1

串行FFT递归算法 国家高性能计算中心（合肥） 2019/1/1

串行FFT递归算法 国家高性能计算中心（合肥） 2019/1/1

串行FFT递归算法 FFT的蝶式递归计算图(由计算原理推出) 国家高性能计算中心（合肥） 2019/1/1

串行FFT递归算法 特别地，n=8的FFT蝶式计算图(展开的) 国家高性能计算中心（合肥） 2019/1/1

串行FFT递归算法 SISD上的FFT分治递归算法 Procedure RFFT(a,b) begin 输入: a=(a0,a1,…,an-1); 输出: B=(b0,b1,…,bn-1) Procedure RFFT(a,b) begin if n=1 then b0=a0 else (1)RFFT(a0,a2,…,an-2, u0,u1,…,un/2-1) (2)RFFT(a1,a3,…,an-1, v0,v1,…,vn/2-1) (3)z=1 (4)for j=0 to n-1 do (4.1)bj=uj mod n/2+zvj mod n/2 (4.2)z=zω endfor 注: (1)算法时间复杂度t(n)=2t(n/2)+O(n) endif 解得 t(n)=O(nlogn) end (2)算法原理？ 国家高性能计算中心（合肥） 2019/1/1

11. 1 快速傅里叶变换(FFT) 11. 1. 1 离散傅里叶变换(DFT) 11. 1. 2 DFT的顺序代码 11. 1 11.1 快速傅里叶变换(FFT) 11.1.1 离散傅里叶变换(DFT) 11.1.2 DFT的顺序代码 11.1.3 串行FFT递归算法 11.1.4 串行FFT非递归算法

串行FFT非递归算法 蝶式计算示例 国家高性能计算中心（合肥） 2019/1/1

串行FFT非递归算法 蝶式计算流图 国家高性能计算中心（合肥） 2019/1/1

串行FFT非递归算法 行0 行1 行2 行3 行4 行5 行6 行7 如: b6=[(a0+a4)-(a2+a6)]-[(a1+a5)-(a3+a7)]ω2 注: ①下行线结点处的权因子的确定问题； ②bi的下标确定:取行号的位序反。如，行3: 3=(011)2 ==>(110)2=6, ==> 行3的输出为b6 国家高性能计算中心（合肥） 2019/1/1

第十一章 快速傅里叶变换 11.1 快速傅里叶变换 11.2 并行FFT算法

11.2 并行FFT算法 11.2.1 SIMD-MC2上的FFT算法 11.2.2 SIMD-BF上的FFT算法

SIMD-MC2上的FFT算法 算法描述 n个处理器组成n1/2×n1/2的方阵, 处理器以行主序编号 国家高性能计算中心（合肥） 2019/1/1

SIMD-MC2上的FFT算法 算法11.3(P270): 输入: ak处于Pk中; 输出bk处于Pk中 Begin (1)for k=0 to n-1 par-do ck=ak endfor (2)for h=logn-1 to 0 do for k=0 to n-1 par-do (2.1)p=2h (2.2)q=n/p (2.3)z=ωp //先算出ωn/2,以后每次z=z1/2 (2.4)if ( k mod p = k mod 2p) then par-do //满足条件的处理器同时做 (i) ck=ck+ck+pzr(k)modq //(i)和(ii)同时执行 (ii)ck＋p=ck-ck+pzr(k)modq endif endfor (3)for k=0 to n-1 par-do bk=cr(k) endfor //r(k)为k的位序反 End 如：n=16 h=3, p=8, q=2, z=ω8 h=2, p=4, q=4, z=ω4 h=1, p=2, q=8, z=ω2 h=0, p=1, q=16,z=ω1 国家高性能计算中心（合肥） 2019/1/1

SIMD-MC2上的FFT算法 示例: P271例11.5, n=4 第(1)步： 国家高性能计算中心（合肥） 2019/1/1

SIMD-MC2上的FFT算法 第(2)步: 第1次迭代(h=1): p=2, q=2, z=ω2 满足k mod 2 = k mod 4的处理器为P0和P1, 同时计算 P0: c0=c0+(ω2)0c2=a0+a2 P1: c1=c1+(ω2)0c3=a1+a3 c2=c0-(ω2)0c2=a0-a2 c3=c1-(ω2)0c3=a1-a3 国家高性能计算中心（合肥） 2019/1/1

SIMD-MC2上的FFT算法 第(2)步: 第2次迭代(h=0): p=1, q=4, z=ω 满足k mod 1 = k mod 2的处理器为P0和P2, 同时计算 P0: c0=c0+ω0c1=(a0+a2)+(a1+a3) P2: c1=c1+ω1c3=(a0+a2)+(a1+a3)ω c1=c0-ω0c1=(a0+a2)-(a1+a3) c3=c1-ω1c3=(a0+a2)+(a1+a3)ω 国家高性能计算中心（合肥） 2019/1/1

SIMD-MC2上的FFT算法 算法分析 第(3)步: b0=c0, b1=c2, b2=c1, b3=c3 计算时间: tcomp=O(logn) 选路时间: trouting: 只涉及(2.4)和(3) (2.4): O(n1/2) (3): O(n1/2) 综上, 当n较大时t(n)=O(n1/2) 国家高性能计算中心（合肥） 2019/1/1

11.2 并行FFT算法 11.2.1 SIMD-MC2上的FFT算法 11.2.2 SIMD-BF上的FFT算法

SIMD-BF上的FFT算法 蝶形网络 处理器布局 有k+1层, 每层有n=2k个处理器, 共有n(1+logn)个处理器 第r行第i列的处理器记为Pr,i, i=(a1,a2,…,ak)2 互连方式 Pr,i与上层Pr-1,i, Pr-1,j相连, 这里i的第r位为0 Pr,j与上层Pr-1,i, Pr-1,j相连, 这里j与i仅在第r位不同 权因子ω在BF网络中的计算方法 Pr,i中ω的指数为j=exp(r,i) 这里exp(r,i)=(ar,…,a1,0,…,0) //即i的前r位取位序反,再后补0 k-r 国家高性能计算中心（合肥） 2019/1/1

SIMD-BF上的FFT算法 示例: n=8的BF网络表示 r,i与上层Pr-1,i, Pr-1,j相连, 这里i的第r位为0 Pr,j与上层Pr-1,i, Pr-1,j相连, 这里j与i仅在第r位不同 国家高性能计算中心（合肥） 2019/1/1

SIMD-BF上的FFT算法 算法描述: P272算法11.4 算法分析 算法的正确性可以用归纳法证明 时间分析 第(1)步时间: O(1) (2.1)和(2.2)的计算时间为O(1), (假定ωexp(r,i)已计算好) (2.1)和(2.2)的选路时间为O(1) ==>第(2)步时间: O(logn) 所以 t(n)=O(logn), p(n)=n(1+logn), c(n)=O(nlog2n) Sp(n)=O(n), Ep(n)=O(1/logn) //本算法的综合指标是较好的 ωexp(r,i)的计算 初始时: Pk,i读入ωexp(k,i)，k=logn 若Pr+1,i已有ωexp(r+1,i), 则Pr,i中的ωexp(r,i)=ω2exp(r+1,i) 所以, 经过logn步就可以计算出每个ωexp(r,i) 国家高性能计算中心（合肥） 2019/1/1