program main implicit none integer::m double precision::dx,Xrange,xa,xb double precision,allocatable::x(:) integer::n double precision::dk,Krange double precision,allocatable::k(:) integer::Np complex(kind(0d0)),allocatable::f(:) complex(kind(0d0)),external::func Np = 200 ! Number of DFT points xa = -exp(1d0) xb = exp(1d0) Xrange = xb-xa ! x=[xa, xa + Xrange] ! Initialize and allocation for DFT allocate(x(0:Np-1),k(0:Np-1),f(0:Np-1)) x = 0d0 k = 0d0 f = dcmplx(0d0,0d0) ! Define position x dx = Xrange / Np call dft_abscissa_shift(Np, x, dx, xa) ! x=[xa, xa + Xrange] ! Calculate k abscissa Krange = 1d0/dx dk = Krange / Np call dft_abscissa(Np, k, dk) ! k=[0, Krange] ! Calculate function f at x=x(m) do m = 0,Np-1 f(m) = func(x(m)) enddo ! Write down f(x) call dft_write("before_fx.d",Np,x,f) ! Forward dft call dft(Np,f,"forward") ! Write down f(k) call dft_writefk_shift("fk.d",Np,k,f,xa) ! output k=[0,Krange] call dft_writefk_center0_shift("fk_center0.d",Np,k,f,xa) ! output k=[-Krange/2,Krange/2] ! Backward dft call dft(Np,f,"backward") ! Write down f(x) call dft_write("after_fx.d",Np,x,f) stop end program main complex(kind(0d0)) function func(x) implicit none double precision,intent(in)::x double precision::pi=dacos(-1d0) func=dcmplx(dcos(2d0*pi*3d0*x),0d0) return end function func !---------------------------------- subroutine dft_abscissa(N, w, dw) implicit none integer,intent(in)::N double precision,intent(in)::dw double precision,intent(out)::w(0:N-1) ! !w=[0, W] ! w(0) = 0 ! w(N-1) = (N-1)*dw = W-dw ! integer::j do j = 0,N-1 w(j) = j*dw enddo return end subroutine dft_abscissa subroutine dft_abscissa_shift(N, w, dw, wa) implicit none integer,intent(in)::N double precision,intent(in)::dw,wa double precision,intent(out)::w(0:N-1) ! !w=[wa, wa + W] ! w(0) = wa ! w(N-1) = (N-1)*dw = W - dw + wa ! integer::j do j = 0,N-1 w(j) = j*dw + wa enddo return end subroutine dft_abscissa_shift include "/opt/mkl/include/mkl_dfti.f90" subroutine dft(N,f,FB) use MKL_DFTI implicit none integer,intent(in)::N complex(kind(0d0)),intent(inout)::f(0:N-1) character(*),intent(in)::FB ! !sikinote ! author : sikino ! date : 2022/06/11 ! !n: Number of DFT points. !f(i): value of data at i !FB: "forward" -> Forward DFT ! "backward" -> Backward DFT integer::Status TYPE(DFTI_DESCRIPTOR),POINTER::hand Status = DftiCreateDescriptor(hand,DFTI_DOUBLE,DFTI_COMPLEX,1,N) Status = DftiSetValue(hand,DFTI_FORWARD_SCALE,1d0) Status = DftiSetValue(hand,DFTI_BACKWARD_SCALE,1d0/dble(N)) Status = DftiCommitDescriptor(hand) if(trim(FB)=="forward")then Status = DftiComputeForward(hand,f) elseif(trim(FB)=="backward")then Status = DftiComputeBackward(hand,f) else write(6,*)"DFT string different" stop endif Status = DftiFreeDescriptor(hand) return end subroutine dft subroutine dft_write(fname,N,w,f) implicit none character(*),intent(in)::fname integer,intent(in)::N double precision,intent(in)::w(0:N-1) complex(kind(0d0)),intent(in)::f(0:N-1) ! w = [wa,W] ! w(0) = wa ! --> output range w=[wa,wa+W] integer::j open(21,file=trim(fname)) do j = 0,N-1 write(21,'(3e26.16e3)')w(j),dble(f(j)),dimag(f(j)) enddo close(21) return end subroutine dft_write subroutine dft_writefk_shift(fname,N,k,f,xa) implicit none character(*),intent(in)::fname integer,intent(in)::N double precision,intent(in)::k(0:N-1),xa complex(kind(0d0)),intent(in)::f(0:N-1) ! k = [0,K] ! k(0) = 0 ! --> output range k=[0,K] integer::nn complex(kind(0d0))::tmp double precision,parameter::pi=dacos(-1d0) open(21,file=trim(fname)) do nn = 0,N-1 tmp = exp(dcmplx(0d0,2d0*pi*k(nn)*xa))*f(nn) write(21,'(3e26.16e3)')k(nn),dble(tmp),dimag(tmp) enddo close(21) return end subroutine dft_writefk_shift subroutine dft_writefk_center0_shift(fname,N,k,f,xa) implicit none character(*),intent(in)::fname integer,intent(in)::N double precision,intent(in)::k(0:N-1),xa complex(kind(0d0)),intent(in)::f(0:N-1) ! k = [0,K] ! k(0) = 0 ! --> output range k=[-K/2,K/2] integer::nn,Nhalf double precision::dk,Krange complex(kind(0d0))::tmp double precision,parameter::pi=dacos(-1d0) if(mod(N,2)==0)then Nhalf = N/2 else Nhalf = (N-1)/2 endif dk = abs(k(1)-k(0)) Krange = N*dk open(21,file=trim(fname)) do nn = Nhalf,N-1 tmp = exp(dcmplx(0d0,2d0*pi*k(nn)*xa))*f(nn) write(21,'(3e26.16e3)')k(nn)-Krange,dble(tmp),dimag(tmp) enddo do nn = 0,Nhalf-1 tmp = exp(dcmplx(0d0,2d0*pi*k(nn)*xa))*f(nn) write(21,'(3e26.16e3)')k(nn),dble(tmp),dimag(tmp) enddo close(21) return end subroutine dft_writefk_center0_shift