Performing split step fourier using MKL dfti

[Nadia_A_]

Hi,

I am trying to perform split step fourier method with the following code. There seems to be a problem with the Fourier portion because when I compare my data with the values obtained with another program, they match perfectly just until I begin the Fourier transform. Here is my code:

program ss

Use MKL_DFTI

implicit none

 real*8 , parameter :: distance = 100.0d0
 real*8 , parameter ::beta2 = 1.0d0
 real*8, parameter :: N = 1.0d0
 real*8, parameter :: nt = 1024.0d0
 real*8, parameter :: Tmax = 32.0d0
 real*8 :: step_num
 real*8 :: deltaz
 real*8 :: dtau
 real*8, dimension (1024) :: tau, omega, uu
 real*8, parameter :: pi = 3.141592653590d0
 type(DFTI_DESCRIPTOR), POINTER :: handle
 complex*16 :: hhz, C= (0,1.0d0)
 complex*16, dimension (1024):: dispersion, xy,temp,temp1,temp2
 integer :: i, status

 
  step_num = 100000
  deltaz = distance/step_num
  dtau = (2.0d0*Tmax)/nt

 do i = 1,1024
  tau (i) = dble(i-513)*dtau
 
  if (i<=512) then
   omega (i) = dble(i-1)*pi/Tmax
  else
   omega (i) = dble(i-1025)*pi/Tmax
  end if
 
  end do
 
  uu = 1/dcosh(tau)
 
  dispersion = exp(0.50d0*C*beta2*(omega**2)*deltaz)
 
  hhz = 1*(0,1.0d0)*N**2*deltaz
 
  temp = uu*exp(abs(uu)**2*hhz/2)
 
 
  status = DftiCreateDescriptor(handle, DFTI_DOUBLE,DFTI_COMPLEX,1,1024)
  status = DftiSetValue(handle,DFTI_FORWARD_SCALE,1.0d0)
  status = DftiSetValue(handle,DFTI_BACKWARD_SCALE,1/dble(1024))
  status = DftiCommitDescriptor(handle)

do j=1,step_num
 
 status = DftiComputeBackward(handle,temp)
 
 do i= 1,1024
 
  temp1(i) = temp(i)*dispersion(i)
 
 end do

 Status = DftiComputeForward(handle,temp1)
 
 do i= 1,1024
 
  temp2(i) = temp1(i)*exp(abs(temp1(i))**2*hhz)
 
 end do
 
end do
 

 status = DftiFreeDescriptor(handle)
 
 
  xy = temp2*exp(-abs(temp1)**2*hhz/2)
 
   do i=1,1024
 
  open (10, file = 'spec.dat')
  write (10,*) tau(i), abs(xy(i)**2)
 
  end do
 
end program ss

[Zhen_Z_Intel]

Hi Nadia,

You compute temp backward for 100000 times. I am afraid the result for 100000 later would be small enough which be seen as 0 value. Thus, temp1 and temp2 also be 0. I tried to print the original temp value, and result for calculating once, the result is fine. And I also found the result after 30 times computing already extend the boundary of smallest double floating number (around 10E-60). Please reduce step_num to see the result, Thanks.

Best regards,
Fiona