数值计算讨论组

欢迎来到数值计算讨论组。 本论坛旨在增强师生交流的时效性,主要讨论数值计算相关问题,包括初等数值分析、矩阵计算、非线性方程(组)数值解法、微分方程数值解法等方面的内容。


    分数阶扩散方程的fast差分格式-O(NlogN)

    分享

    Admin
    Admin

    帖子数 : 46
    注册日期 : 10-05-13

    分数阶扩散方程的fast差分格式-O(NlogN)

    帖子 由 Admin 于 周一 十二月 02, 2013 4:44 pm

    分数阶扩散方程的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官网建议










      目前的日期/时间是周五 十二月 14, 2018 2:39 am