Как сделать «транспонированное» двумерное преобразование fftw3 MPI, если это вообще возможно?

Рассмотрим двумерное преобразование в форме L x M (настройка основного столбца) из сложного массиваЦСИ в реальный массивTGT, Или, на фортранском,

complex(C_DOUBLE_COMPLEX), pointer :: src(:,:)
real(8), pointer :: tgt(:,:)  .

Соответствующие указатели

type(C_PTR) :: csrc,ctgt   .

Я бы выделил их следующим образом:

  ! The complex array first
    alloc_local = fftw_mpi_local_size_2d(M,L/2+1,MPI_COMM_WORLD,local_M,local_offset1)
    csrc = fftw_alloc_complex(alloc_local)
    call c_f_pointer(csrc, src, [L/2,local_M])

    ! Now the real array
    alloc_local = fftw_mpi_local_size_2d(2*(L/2+1),M, &
                                   MPI_COMM_WORLD,local_L,local_offset2)
    ctgt = fftw_alloc_real(alloc_local)
    call c_f_pointer(ctgt, tgt, [M,local_L])

Теперь план будет создан как:

! Create c-->r transform with one transposition left out
plan =  fftw_mpi_plan_dft_c2r_2d(M,L,src,tgt, MPI_COMM_WORLD, & 
                                           ior(FFTW_MEASURE,FFTW_MPI_TRANSPOSED_OUT))

Наконец, преобразование будет выполнено как:

call fftw_mpi_execute_dft_c2r(plan, src, tgt)

Однако этот рецепт не работает. Последний вызов вызывает ошибку сегментации. Сначала я подумал, что это может быть связано с тем, как яЦСИ а такжеTGT массивы, но игра с другим объемом памяти, выделенной дляTGT не дал никакого результата. Итак, я либо делаю что-то действительно глупое, либо это невозможно вообще.

РЕДАКТИРОВАТЬ: МИНИМАЛИСТИЧЕСКИЙ СОВМЕСТИМЫЙ ПРИМЕР

program trashingfftw
  use, intrinsic :: iso_c_binding
  use MPI

  implicit none
  include 'fftw3-mpi.f03'

  integer(C_INTPTR_T), parameter :: L = 256
  integer(C_INTPTR_T), parameter :: M = 256

  type(C_PTR) :: plan, ctgt, csrc

  complex(C_DOUBLE_COMPLEX), pointer :: src(:,:)
  real(8), pointer :: tgt(:,:)

  integer(C_INTPTR_T) :: alloc_local, local_M, &
                         & local_L,local_offset1,local_offset2

  integer :: ierr,id


  call mpi_init(ierr)

  call mpi_comm_rank(MPI_COMM_WORLD,id,ierr)

  call fftw_mpi_init()


  alloc_local = fftw_mpi_local_size_2d(M,L/2+1, MPI_COMM_WORLD, &
       local_M, local_offset1)

  csrc = fftw_alloc_complex(alloc_local)
  call c_f_pointer(csrc, src, [L/2,local_M])


  alloc_local = fftw_mpi_local_size_2d(2*(L/2+1),M, MPI_COMM_WORLD, &
       &                               local_L, local_offset2)

  ctgt = fftw_alloc_real(alloc_local)
  call c_f_pointer(ctgt, tgt, [M,local_L])

  plan =  fftw_mpi_plan_dft_c2r_2d(M,L,src,tgt, MPI_COMM_WORLD, & 
       ior(FFTW_MEASURE, FFTW_MPI_TRANSPOSED_OUT))

  call fftw_mpi_execute_dft_c2r(plan, src, tgt)

  call mpi_finalize(ierr)


end program trashingfftw

Ответы на вопрос(1)

Ваш ответ на вопрос