Transcript FFTW

SCCAS
FFTW使用说明
王 武
中科院超级计算中心
Email: [email protected]
2010 年 5月
1
SCCAS
内容提要
一. 快速傅立叶变换
二. FFTW 的使用指南
三. FFTW 的技术特点
四. FFTW 的调用算例
五. FFTW在Unix上的安装
2
一. 快速傅立叶变换
SCCAS
1.1.1 Fourier Transform (FT)

F (k )   f ( x)e2 kx dx (Forward FT)


f ( x)   F (k )e
2 kx

dk (Inverse FT)
1.1.2 Discrete Fourier Transform (DFT)
One-dimension:
n 1
y  Fn ( x) , yr   xk wn rk , wn  e-2 i / n , r  0,1,
, n -1.
k 0
Multi-dimension:
n1 1 n2 1
y (i1 ,
, id )   
nd 1

k1  0 k2  0
kd  0
-2 i / ni
ni
i
with w  e
x(k1 ,
, r  0,1,
, kd ) wnr11k1 wnr22k2
, ni -1, i  1,
wnr11k1 ,
, d.
3
一. 快速傅立叶变换
SCCAS
1.1.3 DFT in Matrix-vector Form
Eg: DFT in 1D
y  Fx , x,y  R
F is a matrix as in right:
Complexity:
2
O(n )
It’s applied in
DSP ( Digital Signal Processing ):
analyse and manipulate real-world signals by computer
– Toys and consumer electronics (Spectrogram )
– Speech recognition – Audio/video compression
– Medical imaging
– Communications
– Radar
4
一. 快速傅立叶变换
1.2.1 Fast Fourier Transform (FFT)
Cooley & Tukey , 1965
FFT is a divide-and-conquer algorithm
using binary reversal sort
Eg: FFT in 1D
SCCAS
0
1
2
3
4
5
6
7
0
4
2
6
1
5
3
7
Order: even, odd
Complexity:
O( N log 2 N )
5
一. 快速傅立叶变换
SCCAS
1.2.3 Parallel FFT in 1D
Combine: using butterfly computation
1.2.4 FFT in Multi-dimension
2
Eg: DFT in 2D, Complexity : ( n1n2 )
n1 1 n2 1
y (i1 , i2 )    x( k1 , k2 ) w w
r1k1
n1
k1  0 k2  0
n1 1
 w
k1  0
r1k1
n1
n2 1
r2 k2
n2
 x(k , k ) w
k2  0
1
2
r2 k2
n2
[0, 4] [2, 6] [1, 5] [3, 7]
Log28 transform steps
reduced to 2 steps of DFT in 1D
Computed by FFT，Complexity :
O(n1n2 log n1n2 )
6
一. 快速傅立叶变换
SCCAS
1.3 Example: FFT for Convolution Form
N 1
y(k )  x(k )* h(k )   x( j )h(k  j ), k  0,..., N 1
j 0
Step 1: FFT for x(N) and h(N)
N 1
N 1
X ( n)   x ( k )e
2 ikn / N
Step 2: Vec-Vec Mult:
Y(n)=X(n)H(n)
k 0
k 0
Step 3: IFFT for Y(N):
Complexity:
, H ( n)   h( k ) e
2 ikn / N
O ( N log N )
1
y (k ) 
N
N 1
 Y ( n )e
2 ikn / N
n 0
7
二. FFTW的使用指南
SCCAS
2.1 FFTW简单概述
FFTW ( the Faster Fourier Transform in the West)
是一个快速计算离散傅里叶变换的标准C语言程序集
FFTW 由MIT的 M. Frigo 和 S. Johnson 开发,
可计算一维或多维, 实和复数据以及任意规模的DFT.
FFTW 还包含对共享和分布式存储系统的并行变换,
它可自动适应你的机器, 缓存, 存储器大小, 寄存器个数 .
FFTW 通常比目前其它开源Fourier变换程序都要快 !
FFTW 最新版本为fftw-3.2.2
FFTW下载及相关资料提供
http:// www.fftw.org
8
二. FFTW的使用指南
SCCAS
2.2. FFTW 计算:
1) 复(实)一维变换及其逆变换:
n 1
Yi   X j e
n 1
Yi   X j e
j 0
2 ij 1/ n
由于没有正规化: F
j 0
2) 复(实)多维变换及其逆变换:
n1 1 n2 1
Y [i1 , i2 ,
, id ]   
Y [i1 , i2 ,
, id ]   
j1  0 j2  0
n1 1 n2 1
j1  0 j2  0
2 ij 1/ n
nd 1

X [ j1 , j2 ,
jd  0
nd 1

jd  0
X [ j1 , j2 ,
计算一个向前变换随后一个向后变换将用
1
||
F
F ( A)  A || 
变换的正确性检验:
1
F ( X )  nX
, jd ]1 i1 j1 2 i2 j2
 di
, jd ]1i1 j1 2i2 j2
 di j

d
d jd
d d
n 乘以输入(正规化).
s 1 s
3) 以上变换的并行实现
共享和分布式存储系统,一维和多维DFT, 实变换和复变换,多线程和MPI.
9
二. FFTW的使用指南
SCCAS
2.3.1 复一维变换
#include <fftw.h>
…
{ fftw_complex in[N], out[N]; //定义数据类型:复数组
fftw_plan p; //定义计划
…
//创建计划 ,向前变换
p = fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE);
… // ESTIMATE表示不运行任何计算
//而只是建立一个合理的计划
fftw_one(p, in, out); …
//一维变换的输入(出)数组
fftw_destroy_plan(p); //用完后,取消计划
}
// 在Unix上,一个使用FFTW的程序应该与-lfftw –lm相连接
10
二. FFTW的使用指南
SCCAS
2.3.2 复多维变换
#include <fftw.h>
…
//对于FFTW3.0以上版本则采用fftw3.h
{
fftw_complex in[L][M][N];
//数据类型
fftwnd_plan p; … // fftwnd代表N维fftw
p = fftwnd_create_plan(L, M, N, FFTW_FORWARD,
FFTW_MEASURE | FFTW_IN_PLACE);
… //在位(in-place)变换 用输出数组覆盖组
//MEASURE即FFTW实际执行和测量几个FFT的执行时间
以发现计算规模为N的变换的最好方式 ;
fftwnd_one(p, &in[0][0][0], NULL);
…
fftwnd_destroy_plan(p);
}
//程序编译时需与-lfftwnd –lfftw –lm相连接
11
二. FFTW的使用指南
SCCAS
2.3.3 实一维变换
#include <rfftw.h>
…
//专门针对实型数据的rfftw变换
{
// 数据类型
fftw_real in[N], out[N];
rfftw_plan p; //计划
int k …
//创建计划 ,实到复的变换,复到实则为逆变换
p = rfftw_create_plan(N, FFTW_REAL_TO_COMPLEX,
FFTW_ESTIMATE); …
rfftw_one(p, in, out);
rfftw_destroy_plan(p);
}
//程序编译时需与-lrfftw –lfftw –lm相连接
12
二. FFTW的使用指南
SCCAS
2.3.4 实多维变换
#include <rfftw.h>
…
//专门针对实型数据的rfftw变换
{
// 数据类型
fftw_real in[M, N], out[M, N];
rfftw_plan p; //计划
int k …
//创建计划 ,实到复的变换,复到实则为逆变换
p = rfftwnd_create_plan(M, N, FFTW_REAL_TO_COMPLEX,
FFTW_ESTIMATE); …
rfftwnd_one(p, in, out);
rfftwnd_destroy_plan(p);
} //程序编译时需与-lrfftwnd –lfftw –lm相连接
13
二. FFTW的使用指南
SCCAS
2.4.1. FFTW的多线程并行
1.头文件: <fftw_threads.h>或<rfftw_theads.h>
2.线程初始化: int fftw_threads_init(void);
3.用到的函数
fftw_threads_one(nthreads, plan, in, out); //一维复变换
fftwnd_threads_one(nthreads, plan, in, out); //n维复变换
rfftw_threads_one(nthreads, plan, in, out); //一维实变换
rfftwnd_threads_one(nthreads, plan, in, out); //n维实变换
以一维复变换为例
用 fftw_threads_one (nthreads, plan, in, out)
代替调用单机 fftw_one (plan, in, out)
在Unix上,使用并行复变换的程序应该与
-lfftw_threads - lfftw - lm相连接,
使用并行实变换的程序应与
-lrfftw_threads -lfftw_threads -lrfftw -lfftw –lm 相连接
14
二. FFTW的使用指南
SCCAS
2.4.2. FFTW的MPI并行
1. 调用头文件<fftw_mpi.h>
2. 创建计划:
fftw_mpi_plan fftw_mpi_create_plan(MPI_Comm comm, int n
FFTW_FORWARD, FFTW_ESTIMATE );
用完后通过fftw_mpi_destroy_plan(plan)来取消.
3. 变换函数为:
void fftw_mpi(fftw_mpi_plan p, int n_fields, fftw_complex *local_data,
fftw_complex *work);
返回时, local_data包含局限于当前进程的输出部分, 可调用:
void fftw_mpi_local_sizes(fftw_mpi_plan p,
int *local_n, int *local_start,int *local_n_after_transform,
int *local_start_after_transform,int *total_local_size);
在Unix上, FFTW MPI的程序应与MPI库和-lfftw_mpi -lfftw -lm连接.
15
三 FFTW的技术特点
SCCAS
3.1 The Planner
•
For a given N, there are many factorizations
–
–
•
not clear a priori which is best
Eg: 32768 =16x8x8x32=64x16x32
The planner tries them “all” and picks the best one
–
–
•
uses actual runtime timing measurements
result is encoded in a “plan”
Uses dynamic programming to reduce no. of possible plans
–
•
remembers optimal sub-plans for small sizes
The Runtime Planner
–
•
optimizes FFTW for your CPU, your cache size, etc.
Ideas
–
–
–
Modern architectures are invalidating conventional wisdom about
what is fast no new wisdom is emerging
In the name of performance, designers have sacrificed:
predictability, repeatability, composability
Hand-optimization of programs is becoming impractical
16
三 FFTW的技术特点
SCCAS
3.2 The Codelet Generator
•
Generates highly-optimized transforms of small sizes form the
base cases of the FFT recursion
•
Manipulates abstract syntax tree which is unparsed to C
knows about complex arithmetic, etc.
•
The codelets
–
composable blocks of optimized code computer generated
Advantages
•
Long, unrolled code takes advantage of: optimizing compilers
(instruction scheduling, etc.), large register sets
•
Applies tedious optimizations
•
Easy to experiment with different algorithms
–
–
prime factor, split-radix (transform sizes not a power of 2)
various optimization schemes
17
三 FFTW的技术特点
SCCAS
3.3 The Executor
•
•
Executes the plan by composing codelets
Explicit recursion divide-and-conquer
uses all levels of the memory hierarchy
fit in cache
•
Novel storage for the twiddle factors
store them in the order they are used
3.4 FFTW is Easy to Use
COMPLEX A[n], B[n];
fftw_plan plan;
plan = fftw_create_plan(n); /* create the plan */
fftw(plan, A);
/* use the plan */
fftw(plan, B);
/* re-use the plan */
fftw_destroy_plan(plan) /* destroy the plan*/
18
四 FFTW 的调用算例
SCCAS
4.1 调用说明
1.
2.
3.
•
•
•
FFTW 的C函数允许Fortran程序调用 .
fftw/fftwnd/rfftw/rfftwnd 由
fftw_f77/fftwnd_f77/rfftw_f77/rfftwnd_f77所替代
函数的大多数参数相同, 少数例外:
plan变量在C中为fftw_plan, rfftw_plan等类型，
Fortran上对64位机器用integer*8 类型
Fortran数组列存储,C为行存储, integer对应于int, real对应于float
fortran/fftw_f77.i中已将选项参数化
integer FFTW_FORWARD, FFTW_BACKWARD
parameter(FFTW_FORWARD=-1, FFTW_BACKWARD=1
integer FFTW_OUT_OF_PLACE, FFTW_IN_PLACE
parameter(FFTW_OUT_OF_PLACE=0)
parameter(FFTW_IN_PLACE=8)
integer FFTW_ESTIMATE, FFTW_MEASURE
parameter (FFTW_ESTIMATE=0, FFTW_MEASURE=1)
19
四 FFTW 的调用算例
SCCAS
4.2 串行一维复变换的例子
//C编译环境
fftw_complex in[N], *out[N];
fftw_plan plan;
plan = fftw_create_plan(N, FFTW_FORWARD,
FFTW_ESTIMATE);
fftw_one(plan, in, out);
fftw_destroy_plan(plan);
// Fortran编译环境
double complex in, out
dimension in(N), out(N)
integer plan
call fftw_f77_create_plan(plan, N, FFTW_FORWARD,
FFTW_ESTIMATE)
call fftw_f77_one(plan, in, out)
call fftw_f77_destroy_plan(plan)
20
四 FFTW 的调用算例
SCCAS
4.3. 并行一维复变换的例子
#include <fttw_mpi.h>
int main(int argc, char **argv)
{
const int N = 128; fftwnd_mpi_plan plan; fftw_complex *data;
MPI_Comm comm=MPI_COMM_WORLD; MPI_Init(&argc, &argv);
MPI_Comm_rank(comm,&rank);MPI_Comm_size(comm,&np);
plan = fftwnd_mpi_creat_plan(comm, N, -1, 8);
fftw_mpi_local_sizes(plan,&local_n,&nx_out,
&start_x,&start_x_out, &size);
…allocate and initialize data . . . //分配并初始化数据
data= (fftw_complex*) fftw_malloc(sizeof(fftw_complex)* size);
preal=(fftw_real *)in;;srand(rank*100); for(i=0;i<2*local_n;i++)
preal[i]=(fftw_real)rand()/RAND_MAX;
fftwnd_mpi(p, 1, data, NULL, FFTW_NORMAL_ORDER);
fftw_mpi_destroy_plan(plan);
MPI_Finalize( );
21
}
四 FFTW 的调用算例
SCCAS
4.4. FFTW的数值测试
[scwang@LB]$ locate fftw
/soft/fftw-215 fftw_test fftw_mpi_test
[scydf@LB270108 ~/fftw/fftw-f90/2ds_test]$
%Demo
计算时间：T(Forward_plan)+T(Forward_FFT)+T(Inverse_plan)+T(Inverse_FFT)
表1：串行性能（深腾7000, x86_64, 时间: 秒, 浮点性能: Mflops）
问题规模
1024x 1024
2048 x 2048
5120 x 5120
10240 x 10240
计算时间
0.16
5.40039
35.06543
146.90234
浮点性能
1 310 695
788 670
646 024
表2：并行效率（深腾7000, x86_64, 4096 x 4096 ）
647 917
处理机数
1
2
4
8
16
32
时间
6.7383
4.8354
3.2163
3.1541
1.5784
0.7881
效率
1
0.6967
0.5238
0..2670
0.2668
0.2672
22
五. FFTW在Unix上的安装
SCCAS
1) FFTW带有一个GNU形式的configure程序,安装很简单:
tar –xvfz fftw-3.0.1.tar.gz ./configure
make make install
这将建立单机形式的复和实变换库以及测试程序
2) 对于某些系统, configure脚本知道好的CFLAG.
如果你的系统不知道, 可用如下命令编译 FFTW:
./ configure CC=“<Compiler>”, CFLAGS=”<Compiler CFLAGS >”
configure 也可接受一些FFTW说明的标志 :
--enable-float 生成一个单精度版的FFTW
--enable-threads 使能够编译和安装FFTW线程库
--enable-mpi 使能够编译和安装FFTW的MPI库
--disable-fortran 禁止FFTW库中包含Fortran可调用程序
3) 比如调用 dfftw_one或 sfftw_one 之前的安装选项为:
./configure --prefix=/home_soft/soft/ia64/lib/Mathlib/fftw-2.1.5
CC=icc MPICC=icc CFLAGS=-O3 F77=ifort FFLAGS=-O3
--enable-double --enable-type-prefix --enable-mpi --enable-fortran
--enable-shared
(double/float对应于双精度或单精度)
make make install
23
SCCAS
24