并 行 计 算中国科学技术大学计算机科学与技术系国家高性能计算中心 (合肥 )
2003年 9月第三篇 并行数值算法第八章 基本通讯操作第九章 稠密矩阵运算第十章 线性方程组的求解第十一章 快速傅里叶变换第十一章 快速傅里叶变换
11.1快速傅里叶变换
11.2并行 FFT算法
11.1 快速傅里叶变换 (FFT)
11.1.1 离散傅里叶变换 (DFT)
11.1.2 DFT的顺序代码
11.1.3 串行 FFT递归算法
11.1.4串行 FFT非递归算法国家高性能计算中心(合肥) 52009-7-24
离散傅里叶变换 (DFT)
定义给定向量 A=(a0,a1,…,an-1)T,DFT将 A变换为 B=(b0,b1,…,bn-1)T
1
1
0
)1)(1()1(210
1210
0000
1
1
0
/2
1
0;1
10
n
nnnn
n
n
ni
n
k
kj
kj
a
a
a
b
b
b
ine
njab
写成矩阵形式为次单位元根,为=这里即
11.1 快速傅里叶变换 (FFT)
11.1.1 离散傅里叶变换 (DFT)
11.1.2 DFT的顺序代码
11.1.3 串行 FFT递归算法
11.1.4串行 FFT非递归算法国家高性能计算中心(合肥) 72009-7-24
DFT的顺序代码
代码 1
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
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*ω
end for
11.1 快速傅里叶变换 (FFT)
11.1.1 离散傅里叶变换 (DFT)
11.1.2 DFT的顺序代码
11.1.3 串行 FFT递归算法
11.1.4串行 FFT非递归算法国家高性能计算中心(合肥) 92009-7-24
串行 FFT递归算法
蝶式递归计算原理令 为 n/2次单位元根,则有,
将 b向量的偶数项 和奇数项 分别记为和注意推导中反复使用
)2//(2~ nie = 2~=
Tnbbb ),...,,( 131?
Tnbbb ),...,,( 110
2
Tnbbb ),...,,( 110
2
,,1,1,1 2/ ppsnlnnn
ω
8
= ( 1,0 )
ω
6
= ( 0,- i )
ω
2
= ( 0,i )
ω
4
= ( - 1,0 )
ω
7
ω
3
ω
5
ω
1
图 1 1,1
Tnbbb ),...,,( 220?
国家高性能计算中心(合肥) 102009-7-24
串行 FFT递归算法
D F Taaaaaabbb
laa
aa
aaaaaa
aa
aaaaaa
aaaa
aaaa
abb
T
n
T
n
n
k
kk
lk
n
l
ll
n
l
ll
n
lll
lll
n
k
k
lk
ll
nnn
n
n
n
n
nnn
n
n
nnn
n
nnn
n
n
的是向量因此,
偶数时:
),.,,,,(),.,,,,(
1,,1,0)(
~
)(
~
)(
~
)(
~
)(
)(
)()()(
11110220
2
1
0
11
1
22
2
110
11
12
22
4
11
2
0
1
12
2
4
1
2
1
12
2
4
1
2
0
1
0
2
2
222
2
2
2
2
222
2
2
222
2
222
2
2
国家高性能计算中心(合肥) 112009-7-24
串行 FFT递归算法
D F Taaaaaabbb
laa
aaaaaaaa
aaaaaaaa
aaa
aaaa
aaa
aaaa
abb
T
n
T
n
n
k
kk
klk
n
lll
n
lll
n
ll
lll
n
lnll
lll
n
k
k
kl
ll
n
n
nn
n
n
n
nn
nnn
n
nn
nnn
nn
nn
n
nn
n
n
n
n
n
n
的是因此,向量奇数时:
))(),.,,,(),((),.,,,,(
1,,1,0)(
~
)(
~
)(
~
)(
~
)(
)()()()(
11
1
110131
2
1
0
11
11
22
22
110
11
112
22
24
11
2
0
1
112
1
2
1
112
2
24
1
2
0
1
)12(1
1
)12(1)12(
1
)12(1
2
)12(2
1
12
0
1
0
)12(
12
2
2
22
2
2
2
22
222
2
22
222
22
22
2
22
2
2
2
2
2
2
国家高性能计算中心(合肥) 122009-7-24
串行 FFT递归算法
FFT的蝶式递归计算图 (由计算原理推出 )
.
.
.
a
0
a
1
a
n - 1
D F T ( n / 2 )
D F T ( n / 2 )
.
.
.
.
.
.
.
.
.
+
+
+
-
-
-
1
- a ( a
)
ω
n
2
- 1
n
2
- 1
n - 1
- a ( a
)
ω
n
2
+ 1
1
- aa
n
2
0
+ aa
n
2
- 1
n - 1
+ aa
n
2
0
+ aa
n
2
+ 1
1
a
n
2
- 1
a
n
2
a
n
2
+ 1
.
.
.
.
.
.
ω
n
2
- 1
ω
国家高性能计算中心(合肥) 132009-7-24
串行 FFT递归算法特别地,n=8的 FFT蝶式计算图 (展开的 )
a
0
a
1
a
2
a
3
a
4
a
5
a
6
a
7
b
0
b
4
b
2
b
6
b
1
b
5
b
3
7
1
ω
2
ω
ω
3
ω
2
ω
2
1
1
1
1
1
1
+
+
+
+
+
+
+
+
+
+
+
+
_
_
_
_
_
_
_
_
_
_
_
_
b
图 1 1,4
国家高性能计算中心(合肥) 142009-7-24
串行 FFT递归算法
SISD上的 FFT分治递归算法输入,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)算法原理?
11.1 快速傅里叶变换 (FFT)
11.1.1 离散傅里叶变换 (DFT)
11.1.2 DFT的顺序代码
11.1.3 串行 FFT递归算法
11.1.4串行 FFT非递归算法国家高性能计算中心(合肥) 162009-7-24
串行 FFT非递归算法
蝶式计算示例
)()(
)()(
)()(
)()(
00
0101
1010
0101
1100
1100
0011
0011
11
1111
11
1111
1
1
1
1111
,4
11
11
1
11
,2
31203
31201
31202
31200
3
2
1
0
3
1
2
0
21
3
2
1
0
3
2
1
0
963
642
32
3
2
1
0
101
100
1
0
1
0
1
0
aaaab
aaaab
aaaab
aaaab
a
a
a
a
b
b
b
b
bb
a
a
a
a
a
a
a
a
b
b
b
b
n
aab
aab
a
a
a
a
b
b
n
-
-
-
-
=
,和对调
=时当时当国家高性能计算中心(合肥) 172009-7-24
串行 FFT非递归算法
蝶式计算流图
a
0
a
1
a
2
a
3
a
0
a
1
+ a
1
a
0
- a
1
( a
0
) ω
0
( a )
b
0
b
2
b
1
b
3
0
1
0
0
( b )
0
国家高性能计算中心(合肥) 182009-7-24
串行 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
a
0
a
1
a
2
a
3
a
4
a
5
a
6
a
7
b
0
b
4
b
2
b
6
b
1
b
5
b
3
b
7
0
0
0
0
0
0
2
2
1
2
3
0
第十一章 快速傅里叶变换
11.1快速傅里叶变换
11.2并行 FFT算法
11.2 并行 FFT算法
11.2.1 SIMD-MC2上的 FFT算法
11.2.2 SIMD-BF上的 FFT算法国家高性能计算中心(合肥) 212009-7-24
SIMD-MC2上的 FFT算法
算法描述
n个处理器组成 n1/2× n1/2的方阵,处理器以行主序编号
P
0
P
1
P
2
P
4
P
5
P
6
P
8
P
9
P
1 0
P
1 2
P
1 3
P
1 4
P
3
P
7
P
1 1
P
1 5
n
1 / 2
n
1 / 2
图 1 1,5
国家高性能计算中心(合肥) 222009-7-24
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
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
国家高性能计算中心(合肥) 232009-7-24
SIMD-MC2上的 FFT算法
示例,P271例 11.5,n=4
第 (1)步:
c ← a
0
0
c ← a
2
2
c ← a
1
1
c ← a
3
3
P
0
P
2
P
1
P
3
国家高性能计算中心(合肥) 242009-7-24
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
c ← a + a
0
0 2
c ← a + a
1
1 3
c ← a - a
2
0 2
c ← a - a
3
1 3
②
①
①
②
c
2
= a
2
c
3
= a
3
c
2
= a
0
- a
2
c
3
= a
1
- a
3
国家高性能计算中心(合肥) 252009-7-24
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)ω
①
②
c
1
= a
1
+ a
3
c
1
= ( a
0
+ a )
2
- ( a
1
+ a )
3
c ← ( a + a ) + ( a + a )
0
1 32
0
c ← ( a - a ) + ( a - a ) ω
0
1 32
2
c ← ( a + a ) - ( a + a )
0
1 32
1
c ← ( a - a ) - ( a - a ) ω
0
1 32
3
①
②
c
3
= a
1
- a
3
c
3
= ( a
0
- a )
2
- ( a
1
- a ) ω
3
国家高性能计算中心(合肥) 262009-7-24
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)
11.2 并行 FFT算法
11.2.1 SIMD-MC2上的 FFT算法
11.2.2 SIMD-BF上的 FFT算法国家高性能计算中心(合肥) 282009-7-24
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
国家高性能计算中心(合肥) 292009-7-24
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位不同行 0
行 2
行 3
行 1
行 2
行 3
ω
0
ω
0
ω
0
ω
0
ω
4
ω
4
ω
4
ω
4
ω
0
ω
0
ω
4
ω
4
ω
2
ω
2
ω
6
ω
6
ω
0
ω
4
ω
2
ω
6
ω
1
ω
5
ω
3
ω
7
行 1
ω
0
ω
0
ω
0
ω
0
ω
0
ω
0
ω
0
ω
0
列 0 列 1 列 2 列 3 列 4 列 5 列 6 列 7
0 0 0 0 0 1 0 1 0 0 1 1 1 0 0 1 0 1 1 1 0 1 1 1
国家高性能计算中心(合肥) 302009-7-24
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)
2003年 9月第三篇 并行数值算法第八章 基本通讯操作第九章 稠密矩阵运算第十章 线性方程组的求解第十一章 快速傅里叶变换第十一章 快速傅里叶变换
11.1快速傅里叶变换
11.2并行 FFT算法
11.1 快速傅里叶变换 (FFT)
11.1.1 离散傅里叶变换 (DFT)
11.1.2 DFT的顺序代码
11.1.3 串行 FFT递归算法
11.1.4串行 FFT非递归算法国家高性能计算中心(合肥) 52009-7-24
离散傅里叶变换 (DFT)
定义给定向量 A=(a0,a1,…,an-1)T,DFT将 A变换为 B=(b0,b1,…,bn-1)T
1
1
0
)1)(1()1(210
1210
0000
1
1
0
/2
1
0;1
10
n
nnnn
n
n
ni
n
k
kj
kj
a
a
a
b
b
b
ine
njab
写成矩阵形式为次单位元根,为=这里即
11.1 快速傅里叶变换 (FFT)
11.1.1 离散傅里叶变换 (DFT)
11.1.2 DFT的顺序代码
11.1.3 串行 FFT递归算法
11.1.4串行 FFT非递归算法国家高性能计算中心(合肥) 72009-7-24
DFT的顺序代码
代码 1
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
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*ω
end for
11.1 快速傅里叶变换 (FFT)
11.1.1 离散傅里叶变换 (DFT)
11.1.2 DFT的顺序代码
11.1.3 串行 FFT递归算法
11.1.4串行 FFT非递归算法国家高性能计算中心(合肥) 92009-7-24
串行 FFT递归算法
蝶式递归计算原理令 为 n/2次单位元根,则有,
将 b向量的偶数项 和奇数项 分别记为和注意推导中反复使用
)2//(2~ nie = 2~=
Tnbbb ),...,,( 131?
Tnbbb ),...,,( 110
2
Tnbbb ),...,,( 110
2
,,1,1,1 2/ ppsnlnnn
ω
8
= ( 1,0 )
ω
6
= ( 0,- i )
ω
2
= ( 0,i )
ω
4
= ( - 1,0 )
ω
7
ω
3
ω
5
ω
1
图 1 1,1
Tnbbb ),...,,( 220?
国家高性能计算中心(合肥) 102009-7-24
串行 FFT递归算法
D F Taaaaaabbb
laa
aa
aaaaaa
aa
aaaaaa
aaaa
aaaa
abb
T
n
T
n
n
k
kk
lk
n
l
ll
n
l
ll
n
lll
lll
n
k
k
lk
ll
nnn
n
n
n
n
nnn
n
n
nnn
n
nnn
n
n
的是向量因此,
偶数时:
),.,,,,(),.,,,,(
1,,1,0)(
~
)(
~
)(
~
)(
~
)(
)(
)()()(
11110220
2
1
0
11
1
22
2
110
11
12
22
4
11
2
0
1
12
2
4
1
2
1
12
2
4
1
2
0
1
0
2
2
222
2
2
2
2
222
2
2
222
2
222
2
2
国家高性能计算中心(合肥) 112009-7-24
串行 FFT递归算法
D F Taaaaaabbb
laa
aaaaaaaa
aaaaaaaa
aaa
aaaa
aaa
aaaa
abb
T
n
T
n
n
k
kk
klk
n
lll
n
lll
n
ll
lll
n
lnll
lll
n
k
k
kl
ll
n
n
nn
n
n
n
nn
nnn
n
nn
nnn
nn
nn
n
nn
n
n
n
n
n
n
的是因此,向量奇数时:
))(),.,,,(),((),.,,,,(
1,,1,0)(
~
)(
~
)(
~
)(
~
)(
)()()()(
11
1
110131
2
1
0
11
11
22
22
110
11
112
22
24
11
2
0
1
112
1
2
1
112
2
24
1
2
0
1
)12(1
1
)12(1)12(
1
)12(1
2
)12(2
1
12
0
1
0
)12(
12
2
2
22
2
2
2
22
222
2
22
222
22
22
2
22
2
2
2
2
2
2
国家高性能计算中心(合肥) 122009-7-24
串行 FFT递归算法
FFT的蝶式递归计算图 (由计算原理推出 )
.
.
.
a
0
a
1
a
n - 1
D F T ( n / 2 )
D F T ( n / 2 )
.
.
.
.
.
.
.
.
.
+
+
+
-
-
-
1
- a ( a
)
ω
n
2
- 1
n
2
- 1
n - 1
- a ( a
)
ω
n
2
+ 1
1
- aa
n
2
0
+ aa
n
2
- 1
n - 1
+ aa
n
2
0
+ aa
n
2
+ 1
1
a
n
2
- 1
a
n
2
a
n
2
+ 1
.
.
.
.
.
.
ω
n
2
- 1
ω
国家高性能计算中心(合肥) 132009-7-24
串行 FFT递归算法特别地,n=8的 FFT蝶式计算图 (展开的 )
a
0
a
1
a
2
a
3
a
4
a
5
a
6
a
7
b
0
b
4
b
2
b
6
b
1
b
5
b
3
7
1
ω
2
ω
ω
3
ω
2
ω
2
1
1
1
1
1
1
+
+
+
+
+
+
+
+
+
+
+
+
_
_
_
_
_
_
_
_
_
_
_
_
b
图 1 1,4
国家高性能计算中心(合肥) 142009-7-24
串行 FFT递归算法
SISD上的 FFT分治递归算法输入,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)算法原理?
11.1 快速傅里叶变换 (FFT)
11.1.1 离散傅里叶变换 (DFT)
11.1.2 DFT的顺序代码
11.1.3 串行 FFT递归算法
11.1.4串行 FFT非递归算法国家高性能计算中心(合肥) 162009-7-24
串行 FFT非递归算法
蝶式计算示例
)()(
)()(
)()(
)()(
00
0101
1010
0101
1100
1100
0011
0011
11
1111
11
1111
1
1
1
1111
,4
11
11
1
11
,2
31203
31201
31202
31200
3
2
1
0
3
1
2
0
21
3
2
1
0
3
2
1
0
963
642
32
3
2
1
0
101
100
1
0
1
0
1
0
aaaab
aaaab
aaaab
aaaab
a
a
a
a
b
b
b
b
bb
a
a
a
a
a
a
a
a
b
b
b
b
n
aab
aab
a
a
a
a
b
b
n
-
-
-
-
=
,和对调
=时当时当国家高性能计算中心(合肥) 172009-7-24
串行 FFT非递归算法
蝶式计算流图
a
0
a
1
a
2
a
3
a
0
a
1
+ a
1
a
0
- a
1
( a
0
) ω
0
( a )
b
0
b
2
b
1
b
3
0
1
0
0
( b )
0
国家高性能计算中心(合肥) 182009-7-24
串行 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
a
0
a
1
a
2
a
3
a
4
a
5
a
6
a
7
b
0
b
4
b
2
b
6
b
1
b
5
b
3
b
7
0
0
0
0
0
0
2
2
1
2
3
0
第十一章 快速傅里叶变换
11.1快速傅里叶变换
11.2并行 FFT算法
11.2 并行 FFT算法
11.2.1 SIMD-MC2上的 FFT算法
11.2.2 SIMD-BF上的 FFT算法国家高性能计算中心(合肥) 212009-7-24
SIMD-MC2上的 FFT算法
算法描述
n个处理器组成 n1/2× n1/2的方阵,处理器以行主序编号
P
0
P
1
P
2
P
4
P
5
P
6
P
8
P
9
P
1 0
P
1 2
P
1 3
P
1 4
P
3
P
7
P
1 1
P
1 5
n
1 / 2
n
1 / 2
图 1 1,5
国家高性能计算中心(合肥) 222009-7-24
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
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
国家高性能计算中心(合肥) 232009-7-24
SIMD-MC2上的 FFT算法
示例,P271例 11.5,n=4
第 (1)步:
c ← a
0
0
c ← a
2
2
c ← a
1
1
c ← a
3
3
P
0
P
2
P
1
P
3
国家高性能计算中心(合肥) 242009-7-24
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
c ← a + a
0
0 2
c ← a + a
1
1 3
c ← a - a
2
0 2
c ← a - a
3
1 3
②
①
①
②
c
2
= a
2
c
3
= a
3
c
2
= a
0
- a
2
c
3
= a
1
- a
3
国家高性能计算中心(合肥) 252009-7-24
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)ω
①
②
c
1
= a
1
+ a
3
c
1
= ( a
0
+ a )
2
- ( a
1
+ a )
3
c ← ( a + a ) + ( a + a )
0
1 32
0
c ← ( a - a ) + ( a - a ) ω
0
1 32
2
c ← ( a + a ) - ( a + a )
0
1 32
1
c ← ( a - a ) - ( a - a ) ω
0
1 32
3
①
②
c
3
= a
1
- a
3
c
3
= ( a
0
- a )
2
- ( a
1
- a ) ω
3
国家高性能计算中心(合肥) 262009-7-24
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)
11.2 并行 FFT算法
11.2.1 SIMD-MC2上的 FFT算法
11.2.2 SIMD-BF上的 FFT算法国家高性能计算中心(合肥) 282009-7-24
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
国家高性能计算中心(合肥) 292009-7-24
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位不同行 0
行 2
行 3
行 1
行 2
行 3
ω
0
ω
0
ω
0
ω
0
ω
4
ω
4
ω
4
ω
4
ω
0
ω
0
ω
4
ω
4
ω
2
ω
2
ω
6
ω
6
ω
0
ω
4
ω
2
ω
6
ω
1
ω
5
ω
3
ω
7
行 1
ω
0
ω
0
ω
0
ω
0
ω
0
ω
0
ω
0
ω
0
列 0 列 1 列 2 列 3 列 4 列 5 列 6 列 7
0 0 0 0 0 1 0 1 0 0 1 1 1 0 0 1 0 1 1 1 0 1 1 1
国家高性能计算中心(合肥) 302009-7-24
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)