在 Fortran 中使用 fftw3 的大矩阵问题（示例代码）
Problem with big matrices using fftw3 in Fortran (example code)
这个问题来自 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.
你可以下载上次的代码link，也可以写在下面：
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,MM), 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,MM,AA,ATF) ! ATF (complex) is the fft of AA (real)
!
CALL ifft(NN,MM,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*, MINVAL( AA - AA_new ), 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,AA,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, zin, 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,NY,ATF,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, 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) :: planb
!
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)
planb = fftw_plan_dft_2d(NY, NX, zin, zout, 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, zout, zin)
CALL FFTW_EXECUTE_DFT(planb, zin, zout)
!
DO iy = 1, Ny
DO ix = 1, 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 倍。
我无法重现您关于给定代码适用于较小数据集的说法。我得到了预期的缩放比例。
这个问题来自 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.
你可以下载上次的代码link，也可以写在下面：
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,MM), 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,MM,AA,ATF) ! ATF (complex) is the fft of AA (real)
!
CALL ifft(NN,MM,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*, MINVAL( AA - AA_new ), 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,AA,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, zin, 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,NY,ATF,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, 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) :: planb
!
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)
planb = fftw_plan_dft_2d(NY, NX, zin, zout, 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, zout, zin)
CALL FFTW_EXECUTE_DFT(planb, zin, zout)
!
DO iy = 1, Ny
DO ix = 1, 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 倍。
我无法重现您关于给定代码适用于较小数据集的说法。我得到了预期的缩放比例。