在 Fortran 中使用 fftw3 的大矩阵问题示例代码

问题描述

这个问题来自 my last question，但现在是代码。

我在 Fortran 中实现的“西方最快的傅立叶变换”（link）有问题，特别是计算 fft 的倒数。当我用小矩阵测试时，结果是完美的，但从 8x8 开始，结果是错误的。

这是my code here。我在里面写了评论。示例矩阵在文件 ex1.dat,... ex5.dat 中，因此很容易测试（我使用 intel 编译器，我不确定是否使用 gfortran 运行）。示例 ex2 和 ex3 工作完美（5x5 和 7x7），但其他示例给出了错误的结果，所以我无法理解错误或在哪里寻找。

代码内部：为了验证一切是否正确，我计算了

PRINT*,MINVAL( AA - AA_new ),MAXVAL( AA-AA_new )

其中AA是原始矩阵，AA_new是计算fft后的矩阵AA，然后是fft的逆矩阵。我们期望 AA==AA_new，因此我们期望 MINVAL( AA - AA_new ),MAXVAL( AA-AA_new ) 为零。但是对于更大的矩阵，这些数字很大，因此 AA 和 AA_new 非常不同。

我还将结果与 matlab 的命令 fft2 进行比较，因为根据其文档使用相同的库。

关于代码：

主文件是 fft_test.f90，
所有对应于使用fftw3的fft的演算都在fft_modules.f90中。 - _dp（精度）的定义为decimal.f90。

你可以从最后一个链接下载代码，但我也写在下面：

PROGRAM fftw_test
  !
  USE decimal
  USE fftw_module
  !
  IMPLICIT NONE
  !
  REAL(KIND=dp),ALLOCATABLE :: AA(:,:),AA_new(:,:)
  COMPLEX(KIND=dp),ALLOCATABLE :: ATF(:,:)
  INTEGER                       :: NN,MM,ii,jj
  !
  OPEN(unit=10,file='ex2.dat',status='old',action='read')
  READ(10,*) NN,MM ! <- NNxMM is the size of the matrix inside the file vel1.dat
  !
  ALLOCATE( AA(NN,MM),AA_new(NN,ATF(NN,MM) )
  AA = 0.0_dp;  AA_new = 0.0_dp;  ATF = 0.0_dp
  !
  DO ii=1,NN
     READ(10,*) ( AA(ii,jj),jj=1,NN ) ! AA is the original matrix (real)
  END DO
  CLOSE(10)
  !
  PRINT*,'MIN and MAX value of AA' !<- just to verify that is not null
  PRINT*,MINVAL( AA ),MAXVAL( AA )
  !
  CALL fft(NN,AA,ATF) ! ATF (complex) is the fft of AA (real)
  !
  CALL ifft(NN,ATF,AA_new) ! AA_new is the inverse fft AA,so must be == AA
  !
  PRINT*,'MIN and MAX value of AA_new' !<- just to verify that is not null
  PRINT*,MINVAL( AA_new ),MAXVAL( AA_new )
  !
  PRINT*,'Max and min val of (AA-AA_new)'!<- just to compare some way the result
  PRINT*,MAXVAL( AA-AA_new )
  PRINT*,' ------------------------------------------------------------------- '
  !
  DEALLOCATE( AA,ATF )
  !
END PROGRAM fftw_test

MODULE fftw_module
  !
  ! Contains the forward and inverse discrete fourier transform of a matrix
  ! using the library fftw3
  !
  USE decimal
  !
CONTAINS
  !
  SUbroUTINE fft(NX,NY,ATF)
    !
    ! Calculate the forward discrete Fourier transform ATF of the matriz AA
    ! Both matrices have the same size,but AA (the input) is real
    ! and ATF (the output) is complex
    !
    USE,INTRINSIC :: iso_c_binding
    IMPLICIT none
    INCLUDE 'fftw3.f03'
    !   include 'aslfftw3.f03'
    !
    INTEGER(C_INT),INTENT(in)  :: Nx,Ny
    COMPLEX(C_DOUBLE_COMPLEX),ALLOCATABLE :: zin(:),zout(:)
    TYPE(C_PTR)                            :: planf
    !
    INTEGER                                :: ii,yy,ix,iy
    REAL(KIND=dp)                          :: AA(:,:)
    COMPLEX(KIND=dp)                       :: ATF(:,:)
    !
    ALLOCATE( zin(NX * NY),zout(NX * NY) )
    !
    ! Plan Creation (out-of-place forward and backward FFT)
    planf = fftw_plan_dft_2d(NY,NX,zin,zout,FFTW_FORWARD,FFTW_ESTIMATE)
    IF ( .NOT. C_ASSOCIATED(planf) ) THEN
       PRINT*,"plan creation error!!"
       STOP
    END IF
    !
    zin = CMPLX( RESHAPE( AA,(/ NX*NY /) ),KIND=dp )
    !
    ! FFT Execution (forward)
    CALL FFTW_EXECUTE_DFT(planf,zout)
    !
    DO iy = 1,Ny
       DO ix = 1,Nx
          ii = ix + nx*(iy-1)
          ATF(ix,iy) = zout(ii)
       END DO
    END DO
    !
    ! Plan Destruction
    CALL FFTW_DESTROY_PLAN(planf)
    CALL FFTW_CLEANUP
    !
    DEALLOCATE( zin,zout )
    !
  END SUbroUTINE fft
  !
  !
  SUbroUTINE ifft(NX,AA)
    !
    ! Calculate the inverse discrete Fourier transform ATF of the matriz AA
    ! Both matrices have the same size,but AA ATF (the input) is complex
    ! and AA (the output) is real
    !
    USE,zout(:)
    TYPE(C_PTR)                            :: planb
    !
    INTEGER                                :: ii,zout(NX * NY) )
    !
    ! Plan Creation (out-of-place forward and backward FFT)
   planb = fftw_plan_dft_2d(NY,FFTW_BACKWARD,FFTW_ESTIMATE)
    IF ( .NOT. C_ASSOCIATED(planb) ) THEN
       PRINT*,"plan creation error!!"
       STOP
    END IF
    !
    zin = RESHAPE( ATF,(/ NX*NY /) )
    !
    ! FFT Execution (backward)
!    CALL FFTW_EXECUTE_DFT(planb,zin)
    CALL FFTW_EXECUTE_DFT(planb,Nx
          ii = ix + nx*(iy-1)
          AA(ix,iy) = dble( zout(ii) )
       END DO
    END DO
    !
    ! Plan Destruction
    CALL FFTW_DESTROY_PLAN(planb)
    CALL FFTW_CLEANUP
    !
    DEALLOCATE( zin,zout )
    !
  END SUbroUTINE ifft
  !
END MODULE fftw_module

MODULE decimal
  !
  ! Precision in the whole program
  !
  ! dp : double precision
  ! sp : simple precision
  !
  IMPLICIT NONE
  !
  INTEGER,ParaMETER ::  dp = KIND(1.D0)
  INTEGER,ParaMETER ::  sp = KIND(1.0)
  !
END MODULE decimal

一个矩阵示例（因为小而有效），对应于 ex2.dat 输入文件：

5 5
1    2     3    4    5.3
6    7     8    9    10
11   12    3    5    3
0.1 -0.32  0.4  70   12
0 1  0    -1    0    -70

解决方法

当您对某些数据执行前向和后向离散傅立叶变换时，结果的归一化是传统的，通常您要么将数据恢复原样（达到浮点精度），要么如果您使用的是未归一化的转换数据将按数据集中的点数进行缩放，或者您提供归一化作为参数。找出您将阅读用于进行转换的任何软件的文档； fftw uses unnormalised transforms。因此，在您的代码中，您需要进行适当的缩放。如果您在数据集上运行您的代码，您会发现缩放比例如所述 - 在 10x10 数据集上，最终数据是原始数据的 100 倍。

我无法重现您的说法，即给定的代码适用于较小的数据集。我得到了预期的缩放比例。

fftw fortran