Spectral derivatives D2Z transform

Accelerated Computing CUDA CUDA Programming and Performance
1

Good morning everyone,

I have started using cuda fortran for two weeks now, and I must recognize I am stuck with a strange problem.

I want to compute the derivative of a 2D field using spectral derivatives according to the following program :

program example
  
  use cudafor
  use precision_m
  use cufft_m
  
  implicit none 
  
  integer,          parameter :: N = 8
  integer,          parameter :: Nh = N/2 + 1
  double precision, parameter :: Ns = real(N*N)
  double precision, parameter :: pi = 3.141592653589793238d0
  
  double precision, dimension(N,N)  :: signal, dx_signal, dy_signal
  double precision, dimension(N)    :: kx
  double precision, dimension(Nh)   :: ky
  
  double precision, dimension(N,N),  device :: d_signal, d_dx_signal, d_dy_signal
  double complex,   dimension(N,Nh), device :: d_fft,    d_dx_fft,    d_dy_fft
  double precision, dimension(N),    device :: d_kx
  double precision, dimension(Nh),   device :: d_ky

  integer :: i, j, plan
  double complex :: temp
  
  !--- Initialisation
  
  do i = 1, N
    do j = 1, N
      signal(i,j) = cos( 2.*pi * ( 0.*(i-1) + 1.*(j-1) ) / N)
    end do
  end do
  
  kx(1:N/2)   = (/(dfloat(i-1)   * 2. * pi, i=1,N/2,1)/)
  kx(N/2+1:N) = (/(dfloat(i-1-N) * 2. * pi, i=N/2+1,N,1)/)
  ky(1:Nh)    = (/(dfloat(i-1)   * 2. * pi, i=1,Nh,1)/)

  ! Move to device
  d_signal = signal
  d_kx = kx
  d_ky = ky
  
  !--- Computations
  
  ! We compute the Fourier transform
  call cufftPlan2D(plan, N, N, CUFFT_D2Z)
  call cufftExecD2Z(plan, d_signal, d_fft)
  call cufftDestroy(plan)
  
  ! We compute the x derivative
  !$ cuf kernel do(2) <<<*,*>>>
  do i = 1, N
    do j = 1, Nh
      temp = cmplx(0., 1.) * d_kx(i)
      temp = temp * d_fft(i,j)
      d_dx_fft(i,j) = temp / Ns
    end do
  end do
  
  ! We compute the y derivative
  !$ cuf kernel do(2) <<<*,*>>>
  do i = 1, N
    do j = 1, Nh
      temp = cmplx(0., 1.) * d_ky(j)
      temp = temp * d_fft(i,j)
      d_dy_fft(i,j) = temp / Ns
    end do
  end do
  
  ! We compute the Fourier inverse transform
  call cufftPlan2D(plan, N, N, CUFFT_Z2D)
  call cufftExecZ2D(plan, d_dx_fft, d_dx_signal)
  call cufftExecZ2D(plan, d_dy_fft, d_dy_signal)
  call cufftDestroy(plan)
  
  ! Move to host
  dx_signal = d_dx_signal
  dy_signal = d_dy_signal
  
  !--- Print dx_signal and dy_signal
  
  write(*,*) "dx_signal"
  do i = 1, N
    write(*,*) ( dx_signal(i,j), j=1,N )
  end do
  
  write(*,*) "dy_signal"
  do i = 1, N
    write(*,*) ( dy_signal(i,j), j=1,N )
  end do
  
end program example

What is strange is that the derivative according to x should be zero, but this is not what we observe. On the contrary, we can almost find the derivative according to y.

I compile using nvfortran and the following command:
nvfortran -fast -Mpreprocess -Mcuda -c precision_m.f90 -o precision_m.o -lcufft
nvfortran -fast -Mpreprocess -Mcuda -c cufft_m.f90 -o cufft_m.o -lcufft
nvfortran -fast -Mpreprocess -Mcuda -c example.f90 -o example.o -lcufft
nvfortran -fast -Mpreprocess -Mcuda -O3 -o example example.o cufft_m.o precision_m.o -lcufft

And I execute using:
./example

Thank you very much for your help

B.

2

fortran is column major for array storage. CUFFT expects row-major inputs (that is, the first dimension for CUFFT transforms will have elements that are adjacent to each other in memory). Have you accounted for that?

The complex domain (i.e. output) from a D2Z transform will exhibit hermitian symmetry. This has specific implications for the data that will actually be stored - it is not identical to the data storage for an equivalent Z2Z forward transform. Have you accounted for that?

The complex domain must have the same kind of hermitian symmetry for the inverse (Z2D) transform (to be valid). Are you certain that your complex domain arithmetic has that implication?