Блог

У меня возникли проблемы с использованием cufftPlanMany. После создания планов и принятия прямого и обратного БПФ я не смог вернуть исходные данные. Пожалуйста, найдите прилагаемую минимальную версию кода.

 program test_cufft
  use cudafor
  use cufft

  integer :: plan_r2c
  integer :: plan_c2r
  real,allocatable,dimension(:,:,:,:), device :: eta_d
  complex,allocatable,dimension(:,:,:,:), device :: etak_d

  nv = 4
  nx = 256
  ny = 512
  nz = 512
  nx21 = nx/2 1

  allocate( eta_d(nv,nx,ny,nz) )
  allocate( etak_d(nv,nx21,ny,nz) )

  batch = nv;
  rank = 3;
  n = (/ nx, ny, nz /);
  idist = nx*ny*nz;
  odist = nx21*ny*nz;
  inembed = (/ nx, ny, nz /);
  onembed = (/ nx21, ny, nz /);
  istride = 1;
  ostride = 1;

  istat = cufftPlanMany( plan_r2c, rank, n, inembed, istride, idist, amp;
                         onembed, ostride, odist, CUFFT_R2C, batch )
  istat = cufftPlanMany( plan_c2r, rank, n, onembed, ostride, odist, amp;
                         inembed, istride, idist, CUFFT_C2R, batch )

  ! Initialize eta_d

  istat = cufftExecR2C( plan_r2c, eta_d, etak_d )
  istat = cufftExecC2R( plan_c2r, etak_d, eta_d )
  eta_d = eta_d/idist

end program test_cufft
 

Проблема в том, что после того, как я выполнил прямое и обратное БПФ, я не смог вернуть исходные данные.
Пожалуйста, что я делаю не так?
Должен ли порядок данных быть eta_d(batch,nx,ny,nz) or eta_d(nx,ny,nz,batch) ?

Я бы сказал, что правильный порядок (nz, ny, nx, batch)

Но важно также соотнести их с индексацией вашего массива и порядком хранения.

В терминологии CUFFT, для 3D-преобразования (*) nz направление является наиболее быстро изменяющимся индексом, при этом типичное использование (шаг = 1) является смежными данными в памяти, соответствующими смежным элементам в преобразовании.

Это направление (я думаю о нем как об элементах вдоль строки, т. Е. Индекс «z» является индексом столбца) также является направлением многомерного преобразования, которое «уменьшается» в сложной области для типов преобразований R2C / C2R.

Имея это в виду, я бы переписал ваш код таким образом:

 $ cat t4.cuf
program test_cufft
  use cudafor
  use cufft

  integer :: plan_r2c
  integer :: plan_c2r
  real,allocatable,dimension(:,:,:,:), managed :: eta_d
  complex,allocatable,dimension(:,:,:,:), managed :: etak_d
  integer :: n(3), inembed(3), onembed(3),rank,istride,idist,ostride,odist,batch

  nv = 4
  nx = 8
  ny = 8
  nz = 4
  nz21 = nz/2 1

  allocate( eta_d(nz,ny,nx,nv) )
  allocate( etak_d(nz21,ny,nx,nv) )

  batch = nv;
  rank = 3;
  n = (/ nx, ny, nz /);
  idist = nx*ny*nz;
  odist = nx*ny*nz21;
  inembed = (/ nx, ny, nz /);
  onembed = (/ nx, ny, nz21 /);
  istride = 1;
  ostride = 1;

  istat = cufftPlanMany( plan_r2c, rank, n, inembed, istride, idist, onembed, ostride, odist, CUFFT_R2C, batch )
  istat = cufftPlanMany( plan_c2r, rank, n, onembed, ostride, odist, inembed, istride, idist, CUFFT_C2R, batch )

  ! Initialize eta_d
  eta_d(:,:,:,:) = 1.0
  eta_d(1,1,1,2) = 2.0
  istat = cufftExecR2C( plan_r2c, eta_d, etak_d )
  istat = cudaDeviceSynchronize()
  eta_d(:,:,:,:) = 0.0
  istat = cufftExecC2R( plan_c2r, etak_d, eta_d )
  istat = cudaDeviceSynchronize()
  eta_d = eta_d/idist
  print *,eta_d(1,1,1,1)
  print *,eta_d(1,1,1,2)
end program test_cufft
$ nvfortran t4.cuf -lcufft
$ ./a.out
    1.000000
    2.000000
$
 

(NVIDIA HPC SDK 20.9, графический процессор Tesla V100)

и, похоже, это дает ожидаемые результаты для моего простого тестового примера.

(*) Для 2D-преобразования ny размерность будет наиболее быстро изменяться, а для одномерного преобразования nx размерность (конечно) будет наиболее быстро изменяться.

Разделы многомерных преобразований и расширенной компоновки данных руководства CUFFT также могут быть полезны для чтения.