分数阶扩散方程的fast差分格式-O(NlogN)
关键在于fftw的应用
1. 见D:\study\FFTW, 在lib目录下,有一个静态链接库,fftw3fortran.lib; 在include目录下,有一个头文件fftw3.f
(用于一些常值设定)
把这两个目录在tools->options-> 添加入lib and include搜索地址中去
在fftw3.f中可以看到
INTEGER FFTW_FORWARD
PARAMETER (FFTW_FORWARD=-1)
INTEGER FFTW_BACKWARD
PARAMETER (FFTW_BACKWARD=+1)
+1,-1分别表示正逆fft
2. 建立测试workspace,内含主程序x1.f90
3. 在此workspace下,添加 fftw3.dll 动态链接库
4. 在alt+F7 settings 里面,添加静态链接库 fftw3fortran.lib
5. compile, build and execute 得到结果
*******************************************************
! x1.f90
program main
include 'fftw3.f'
integer :: N
parameter (N=32)
double complex out(N)
double complex in(N)
double complex var1(N)
integer :: plan
integer :: i
!
do i=1,N
in(i)=i*(1.d0,0.d0)
end do
!
call dfftw_plan_dft_1d(plan,N,in,out,1,FFTW_ESTIMATE) ! in ---> out 1 means ifft
call dfftw_execute(plan)
call dfftw_destroy_plan(plan)
!
call dfftw_plan_dft_1d(plan,N,out,var1,-1,FFTW_ESTIMATE) ! out ---> var1 -1 means fft
call dfftw_execute(plan)
call dfftw_destroy_plan(plan)
!
do i = 1,N
write(*,*) i, in(i), out(i)
end do
end program main
**************************************************************************
其中,31行,dfftw_plan_dft_1d中的参数plan 应该是输出参数,用于fourier transform的计划
N 一维数组的长度
in 要做离散变换的原始数据
out 该离散变换的结果
1 也可以用FFTW_BACKWARD代替,表明是逆变换
FFTW_ESTIMATE: 这个属于flags参数,
“
flags参数一般情况下为FFTW_MEASURE 或 FFTW_ESTIMATE。FFTW_MEASURE表示FFTW会先计算一些FFT并测量所用的时间,以
便为大小为n的变换寻找最优的计算方法。依据 机器配置和变换的大小(n),这个过程耗费约数秒(时钟clock精度)。
FFTW_ESTIMATE则相反,它直接构造一个合理的但可能是次最优的方 案。
见 FFTW中文参考 http://blog.csdn.net/sjzcandy/article/details/5866409
”
*****************************************************************************
关于正逆变换与matlab里面的关系
FFTW中的逆变换 (标记为1) = matlab中的逆变换*N
FFTW中的正变换 (标记为-1) = matlab中的正变换
*****************************************************************************
2013.10.15
dfftw_plan_dft_r2c_1d(plan,n,x,y,fftw_measure)
dfftw_execute_dft_r2c(plan,x,y)
r2c版本:实输入数据,复Hermitian输出,正变换, 注意,这里的输出数据,只有1半的存储。
实型,双精度x -----> 复型,双精度y (若为双精度,子程序的名字第一个单词为dfftw,若为单精度,则fftw)
n = length(x) = length(y)
整型 fftw_measure = 0
整型 plan
由于只有一半的存储,所以y作为输出数据,仅仅存了第一个元素,与后面一半的元素。
比如 n = 5时, 只存y(1:3),因为y(4) = y(3)*, y(5) = y(2)*, *表示共轭
比如 n = 6时, 只存y(1:4),因为y(5) = y(3)*, y(6) = y(2)*, y(4)独立
如果想将y共轭对称出其他的元素,写一个subroutine :VectorHermit来处理
------------------------------------------------------
dfftw_plan_dft_c2r_1d(plan,n,y,x,fftw_measure)
dfftw_execute_dft_c2r(plan,y,x)
c2r版本:复Hermitian输入数据,实输出数据,逆变换
复型双精度 y ----> 实型 双精度 x
n = length(x) = length(y)
注意 y是输入变量,x为输出变量
------------------------------------------------------
测试程序见 fftw_test.f90
编译方式为 gfortran -o app fftw_test.f90 -lfftw3
------------------------------------------------------
发现了一个问题
c2r,r2c这两个子程序,对较长的向量处理的时候,结果都为0
通过查询fftw web examples for fftw.org
To transform a one-dimensional real array in Fortran, you might do:
double precision in
dimension in(N)
double complex out
dimension out(N/2 + 1)
integer*8 plan
call dfftw_plan_dft_r2c_1d(plan,N,in,out,FFTW_ESTIMATE)
call dfftw_execute_dft_r2c(plan, in, out)
call dfftw_destroy_plan(plan)
这说明, 进行r2c的fft时, 输出变量的长度为N/2+1, 这意味着 N 需要是偶数!
再重新测试!
2013.10.16
注意一个重要的问题,最后一个参数是FFTW_ESTIMATE,这个整型常量=64
而不是FFTW_MEASURE = 0
问题出现在这里,win下面因为include ‘fftw3.f’, 已经声明了这个常量,所以可以直接运行
但是linux下面需要自己定义这个常量,或者把fftw3.f直接copy到 目录下面去,然后在main.f90里面include一下
测试程序,见 fftw_test.f90
再强调一遍,在
dfftw_plan_dft_r2c_1d(plan,n,x,y,fftw_estimate)
或
dfftw_plan_dft_c2r_1d(plan,n,y,x,fftw_estimate)
中,y的长度为n/2+1, x 为double precision, y为 double complex
最后的测试,
比较一下,1024维的x再进行fft变换的时候需要多少时间,x为随机生成
两种方式,(I)为 每次变换都重新plan (II) 只做一次plan,用到每次的变换上面去
用时 (I) (II)
次数 100 3.1e-e 0.0e-2
次数 500 7.8e-e 1.5e-2
次数 1000 1.5e-1 4.7e-2
次数 5000 7.3e-1 1.7e-1
进行逆变换时,类似ifft。注意,应用时,plan是否应该作为全局变量
!----------------------------------------------------
2013.10.17
最后注意2件非常重要的事情
1. 在使用fftw进行FFT变换的时候,需要将输入参数(向量)的下标调整为 0:n-1,不然会出现问题。
fftw是区分这种下标移位的!!!! 包括r2c,c2r
!----------------------------------------------------
2013.10.18
2. 如果多次对相同维数的向量进行FFT或者iFFT(使用FFTW package),那么相应的plan可以只生成一次。
不要每次进行变换时都重新生成plan,这样会非常费时(经测试)。作为一个例子,见subroutine Fm in
FractionalPackage
!-----------------------------------------------------
2013.12.2补充
integer*8 :: FmFFTplan,FmiFFTplan
注意声明变量plan时,用integer*8,不然有时结果算不对,见fftw3官网建议
关键在于fftw的应用
1. 见D:\study\FFTW, 在lib目录下,有一个静态链接库,fftw3fortran.lib; 在include目录下,有一个头文件fftw3.f
(用于一些常值设定)
把这两个目录在tools->options-> 添加入lib and include搜索地址中去
在fftw3.f中可以看到
INTEGER FFTW_FORWARD
PARAMETER (FFTW_FORWARD=-1)
INTEGER FFTW_BACKWARD
PARAMETER (FFTW_BACKWARD=+1)
+1,-1分别表示正逆fft
2. 建立测试workspace,内含主程序x1.f90
3. 在此workspace下,添加 fftw3.dll 动态链接库
4. 在alt+F7 settings 里面,添加静态链接库 fftw3fortran.lib
5. compile, build and execute 得到结果
*******************************************************
! x1.f90
program main
include 'fftw3.f'
integer :: N
parameter (N=32)
double complex out(N)
double complex in(N)
double complex var1(N)
integer :: plan
integer :: i
!
do i=1,N
in(i)=i*(1.d0,0.d0)
end do
!
call dfftw_plan_dft_1d(plan,N,in,out,1,FFTW_ESTIMATE) ! in ---> out 1 means ifft
call dfftw_execute(plan)
call dfftw_destroy_plan(plan)
!
call dfftw_plan_dft_1d(plan,N,out,var1,-1,FFTW_ESTIMATE) ! out ---> var1 -1 means fft
call dfftw_execute(plan)
call dfftw_destroy_plan(plan)
!
do i = 1,N
write(*,*) i, in(i), out(i)
end do
end program main
**************************************************************************
其中,31行,dfftw_plan_dft_1d中的参数plan 应该是输出参数,用于fourier transform的计划
N 一维数组的长度
in 要做离散变换的原始数据
out 该离散变换的结果
1 也可以用FFTW_BACKWARD代替,表明是逆变换
FFTW_ESTIMATE: 这个属于flags参数,
“
flags参数一般情况下为FFTW_MEASURE 或 FFTW_ESTIMATE。FFTW_MEASURE表示FFTW会先计算一些FFT并测量所用的时间,以
便为大小为n的变换寻找最优的计算方法。依据 机器配置和变换的大小(n),这个过程耗费约数秒(时钟clock精度)。
FFTW_ESTIMATE则相反,它直接构造一个合理的但可能是次最优的方 案。
见 FFTW中文参考 http://blog.csdn.net/sjzcandy/article/details/5866409
”
*****************************************************************************
关于正逆变换与matlab里面的关系
FFTW中的逆变换 (标记为1) = matlab中的逆变换*N
FFTW中的正变换 (标记为-1) = matlab中的正变换
*****************************************************************************
2013.10.15
dfftw_plan_dft_r2c_1d(plan,n,x,y,fftw_measure)
dfftw_execute_dft_r2c(plan,x,y)
r2c版本:实输入数据,复Hermitian输出,正变换, 注意,这里的输出数据,只有1半的存储。
实型,双精度x -----> 复型,双精度y (若为双精度,子程序的名字第一个单词为dfftw,若为单精度,则fftw)
n = length(x) = length(y)
整型 fftw_measure = 0
整型 plan
由于只有一半的存储,所以y作为输出数据,仅仅存了第一个元素,与后面一半的元素。
比如 n = 5时, 只存y(1:3),因为y(4) = y(3)*, y(5) = y(2)*, *表示共轭
比如 n = 6时, 只存y(1:4),因为y(5) = y(3)*, y(6) = y(2)*, y(4)独立
如果想将y共轭对称出其他的元素,写一个subroutine :VectorHermit来处理
------------------------------------------------------
dfftw_plan_dft_c2r_1d(plan,n,y,x,fftw_measure)
dfftw_execute_dft_c2r(plan,y,x)
c2r版本:复Hermitian输入数据,实输出数据,逆变换
复型双精度 y ----> 实型 双精度 x
n = length(x) = length(y)
注意 y是输入变量,x为输出变量
------------------------------------------------------
测试程序见 fftw_test.f90
编译方式为 gfortran -o app fftw_test.f90 -lfftw3
------------------------------------------------------
发现了一个问题
c2r,r2c这两个子程序,对较长的向量处理的时候,结果都为0
通过查询fftw web examples for fftw.org
To transform a one-dimensional real array in Fortran, you might do:
double precision in
dimension in(N)
double complex out
dimension out(N/2 + 1)
integer*8 plan
call dfftw_plan_dft_r2c_1d(plan,N,in,out,FFTW_ESTIMATE)
call dfftw_execute_dft_r2c(plan, in, out)
call dfftw_destroy_plan(plan)
这说明, 进行r2c的fft时, 输出变量的长度为N/2+1, 这意味着 N 需要是偶数!
再重新测试!
2013.10.16
注意一个重要的问题,最后一个参数是FFTW_ESTIMATE,这个整型常量=64
而不是FFTW_MEASURE = 0
问题出现在这里,win下面因为include ‘fftw3.f’, 已经声明了这个常量,所以可以直接运行
但是linux下面需要自己定义这个常量,或者把fftw3.f直接copy到 目录下面去,然后在main.f90里面include一下
测试程序,见 fftw_test.f90
再强调一遍,在
dfftw_plan_dft_r2c_1d(plan,n,x,y,fftw_estimate)
或
dfftw_plan_dft_c2r_1d(plan,n,y,x,fftw_estimate)
中,y的长度为n/2+1, x 为double precision, y为 double complex
最后的测试,
比较一下,1024维的x再进行fft变换的时候需要多少时间,x为随机生成
两种方式,(I)为 每次变换都重新plan (II) 只做一次plan,用到每次的变换上面去
用时 (I) (II)
次数 100 3.1e-e 0.0e-2
次数 500 7.8e-e 1.5e-2
次数 1000 1.5e-1 4.7e-2
次数 5000 7.3e-1 1.7e-1
进行逆变换时,类似ifft。注意,应用时,plan是否应该作为全局变量
!----------------------------------------------------
2013.10.17
最后注意2件非常重要的事情
1. 在使用fftw进行FFT变换的时候,需要将输入参数(向量)的下标调整为 0:n-1,不然会出现问题。
fftw是区分这种下标移位的!!!! 包括r2c,c2r
!----------------------------------------------------
2013.10.18
2. 如果多次对相同维数的向量进行FFT或者iFFT(使用FFTW package),那么相应的plan可以只生成一次。
不要每次进行变换时都重新生成plan,这样会非常费时(经测试)。作为一个例子,见subroutine Fm in
FractionalPackage
!-----------------------------------------------------
2013.12.2补充
integer*8 :: FmFFTplan,FmiFFTplan
注意声明变量plan时,用integer*8,不然有时结果算不对,见fftw3官网建议