i'm having issue computing fft of data in fortran. don't know if there's wrong algorithm, roundoff, lack of precision or what.
here code
module fft_mod public :: fft1d integer,parameter :: cp = selected_real_kind(14) real(cp),parameter :: pi = real(3.14159265358979,cp) contains subroutine fft1d(x) complex(cp), dimension(:), intent(inout) :: x complex(cp), dimension(:), allocatable :: temp complex(cp) :: s integer :: n complex(cp) :: j ! sqrt(-1) integer :: i,k n=size(x) allocate(temp(n)) j = cmplx(0.0_cp,1.0_cp,cp) s = cmplx(0.0_cp,0.0_cp,cp) temp = cmplx(0.0_cp,0.0_cp,cp) = 1,n k = 1,n s = s + x(k)*exp(real(-2.0,cp)*pi*j*real(k-1,cp)*real(i-1,cp)/real(n,cp)) enddo temp(i) = s s = cmplx(0.0_cp,0.0_cp,cp) enddo x = temp deallocate(temp) end subroutine end module program test use fft_mod implicit none complex(cp), dimension(10) :: data = (/1.0, 2.0, 3.0, 4.0, 5.0, 5.0, 4.0, 3.0, 2.0, 1.0/) integer :: call fft1d(data) i=1,10 write(*,*) data(i) end end program test running in fortran gives
c:\users\charlie\desktop\fft>gmake gfortran -j".\\mod" -fimplicit-none -wuninitialized -g -wall -wextra -fbacktrace -fcheck=all -o0 -fopenmp -d_quad_precision_ -cpp -c -o .\\obj\testfft.o testfft .f90 gfortran -o .\\test -j".\\mod" -fimplicit-none -wuninitialized -g -wall -wextra -fbacktrace -fcheck=all -o0 -fopenmp -d_quad_precision_ -cpp .\\obj\testfft.o c:\users\charlie\desktop\fft>test.exe ( 30.000000000000000 , 0.0000000000000000 ) ( -9.4721355260035178 , -3.0776825738331275 ) ( 1.2032715918097736e-007, 8.7422769579070803e-008) (-0.52786408204828272 ,-0.72654221835813126 ) ( 5.6810824045072650e-008, 1.7484555003832725e-007) ( 1.0325074129013956e-013, 2.6226834001115759e-007) ( -8.5216018574918451e-008, 2.6226836247200680e-007) (-0.52786395855490920 , 0.72654325051559143 ) ( -4.8130813040669906e-007, 3.4969128892559098e-007) ( -9.4721398159606647 , 3.0776922072585111 ) but running same dataset in matlab gives
format long ; g = [1:5 5:-1:1]; fft(g)' ans = 30.000000000000000 -9.472135954999580 + 3.077683537175253i 0 -0.527864045000420 + 0.726542528005361i 0 0 0 -0.527864045000420 - 0.726542528005361i 0 -9.472135954999580 - 3.077683537175253i i believe i'm using double precision using selected_real_kind(14), looks result single precision @ best. i'm sure of sprinkled real(,cp)'s not necessary, don't know where, how or why looks result single precision compared matlab.
any appreciated!
update:
using advice accepted answer, thing needed change here was:
real(cp),parameter :: pi = real(3.14159265358979,cp) to
real(cp),parameter :: pi = 3.14159265358979_cp
the problem how define real numbers, pi. when define
real(cp),parameter :: pi = real(3.14159265358979,cp) you're passing argument 3.14159265358979 function real. real numbers have default single precision real number cast single precision enters function. consider following example:
program main integer,parameter :: cp = selected_real_kind(14) real(cp),parameter :: pi = real(3.14159265358979,cp) real(cp),parameter :: pj = 3.14159265358979_cp write(*,*) pi write(*,*) pj end program main compiled pgfortran , no options, gives me:
3.141592741012573 3.141592653589790 when defining real number, should use []_cp assign kind instead of real([],cp).
edit: problem affects how define 0.0, 1.0, , 2.0, numbers may cast binary , not suffer same rounding error.
Comments
Post a Comment