分解在FORTRAN77

Question

[编辑1]加数字来显示原始数据和所获得的数据

[编辑2]我发现我的错误，我用来代替fftw_estimate fftw_measure在的呼叫dfftw_plan_many_dft

中的代码（用dfftw_execute_dft_r2c U2D替换U）

我试图执行使用多个1D FFT阵列的2D FFT而不是使用

[编辑3]修正笔误2D fft函数已经存在于fftw库中。随后，我需要执行逆2D fft。 这样做的原因是（将来）我的数组将会太大而无法一次加载来执行2D fft。

我的代码的第一个草案看上去或多或少像这样的时刻：

double precision u2d(nx,ny),u2d2(nx,ny) 
    double complex qhat2d(nx/2+1,ny),qhat2d2(nx/2+1,ny) 
    integer N(1) 
    integer howmany, idist, odist, istride, ostride 
    integer inembed(2), onembed(2) 
    integer rank 

    ! some function to read the data into u2d 

    ! perform x-fft 
    N(1) = NX 
    howmany = NY 
    inembed(1) = NX 
    inembed(2) = NY 
    istride = 1 
    idist = NX 
    ostride = 1 
    odist = (NX/2+1) 
    onembed(1) = (NX/2+1) 
    onembed(2) = NY 
    rank = 1 
    write(*,*) 'u', u2d(1,1) 
    CALL dfftw_plan_many_dft_r2c(PLAN,rank,N(1),howmany, 
&        u2d,inembed, 
&        istride,idist, 
&        qhat2d,onembed, 
&        ostride,odist,FFTW_ESTIMATE) ! 
    CALL dfftw_execute_dft_r2c(PLAN,u2d,qhat2d) ! x-fft 
    CALL dfftw_destroy_plan(PLAN) 


    ! perform y-fft 
    N(1) = NY 
    howmany = (NX/2+1) 
    inembed(1) = (NX/2+1) 
    inembed(2) = NY 
    istride = (NX/2+1) 
    idist = 1 
    ostride = (NX/2+1) 
    odist = 1 
    onembed(1) = (NX/2+1) 
    onembed(2) = NY 
    rank = 1 
    CALL dfftw_plan_many_dft(PLAN,rank,N(1),howmany, 
&        qhat2d,inembed, 
&        istride,idist, 
&        qhat2d2,onembed, 
&        ostride,odist,FFTW_FORWARD, 
&        FFTW_MEASURE) ! 
    CALL dfftw_execute_dft(PLAN,qhat2d,qhat2d2) ! y-fft 
    CALL dfftw_destroy_plan(PLAN) 

    ! normally here, perform some filtering operation 
    ! but at the moment, I do nothing 

    ! perform inv-y-fft 
    N(1) = NY 
    howmany = (NX/2+1) 
    inembed(1) = (NX/2+1) 
    inembed(2) = NY 
    istride = (NX/2+1) 
    idist = 1 
    ostride = (NX/2+1) 
    odist = 1 
    onembed(1) = (NX/2+1) 
    onembed(2) = NY 
    rank = 1 
    CALL dfftw_plan_many_dft(PLAN,rank,N(1),howmany, 
&        qhat2d2,inembed, 
&        istride,idist, 
&        qhat2d,onembed, 
&        ostride,odist,FFTW_BACKWARD, 
&        FFTW_MEASURE) ! 
    CALL dfftw_execute_dft(PLAN,qhat2d2,qhat2d) ! inv-y-fft 
    CALL dfftw_destroy_plan(PLAN) 

    ! perform inv-x-fft 
    N(1) = NX ! I'm not too sure about this value here 
    howmany = NY 
    inembed(1) = (NX/2+1) 
    inembed(2) = NY 
    istride = 1 
    idist = (NX/2+1) 
    ostride = 1 
    odist = NX 
    onembed(1) = NX 
    onembed(2) = NY 
    rank = 1 
    CALL dfftw_plan_many_dft_c2r(PLAN,rank,N(1),howmany, 
&        qhat2d,inembed, 
&        istride,idist, 
&        u2d2,onembed, 
&        ostride,odist,FFTW_ESTIMATE) ! 
    CALL dfftw_execute_dft_c2r(PLAN,qhat2d,u2d2) ! x-fft 
    CALL dfftw_destroy_plan(PLAN) 
    write(*,*) 'u2d2', u2d2(1,1) 

    do i=1,nx 
    do j=1,ny 
    u2d2(i,j) = u2d2(i,j)/(nx*ny) 
    enddo 
    enddo 
    write(*,*) 'u2d2', u2d2(1,1) ! here the values u2d2(1,1) is different from u2d(1,1) 

    ! some action to write u2d2 to file 
    end 

我期待U2D和u2d2是相同的，但我获得相对不同的值。我在某个地方犯了错吗？

原文和结果如下所示。形状看起来相似，但数值相对不同（例如最小值和最大值）。

Original field

Field obtained after fft and i-fft

Answer 1

清除混淆。会发生什么情况是FFTW c2r和r2c例程不能保证保留输入数组。他们可以用垃圾覆盖结果。

现在，你可以很幸运，只需使用FFTW_ESTIMATE而不是FFTW_MEASURE，但这不是一个好主意。 FFTW_MEASURE尝试了很多算法，因此可能会尝试一种不保留输入的算法。 FFTW_ESTIMATE将不会尝试计算任何内容并且不会覆盖输入。

问题是，执行转换（执行计划）时可能会覆盖您的输入。只有当FFT_ESTIMATE选择一种为您保留输入的算法时，您才会幸运。这是一个彩票。

为确保输入得到保留，除了FFT_ESTIMATE或FFTW_MEASURE之外，还应使用FFTW_INPUT_PRESERVE。

您也可以不使用它，而是将输入保存在某处。我认为这通常更好，因为FFTW_INPUT_PRESERVE可以（很可能）选择较慢的算法。

Answer 2

，我发现我的错误，我用fftw_measure代替fftw_estimate在dfftw_plan_many_dft的召唤。

纠正，给了我适当的原始字段。