Как выполнить транспонированное 2D-преобразование MPI fftw3, если это вообще возможно?

Рассмотрим двумерное преобразование формы L x M (основная настройка столбца) из сложного массива src в реальный массив 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)

Однако этот рецепт не работает. Последний вызов вызывает ошибку сегментации. Сначала я подумал, что это может иметь какое-то отношение к тому, как я выделяю массивы src и 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

c fortran mpi fftw
person Mali Remorker    schedule 18.06.2014    source источник
comment
У вас есть автономный фрагмент кода, который мы можем попытаться скомпилировать и запустить?   -  person Vladimir F    schedule 18.06.2014
comment
У меня его пока нет, но я постараюсь его предоставить. А пока я надеюсь, что кто-то, кто действительно использовал транспонированные реальные преобразования FFTW аналогичным образом, может дать пару советов.   -  person Mali Remorker    schedule 18.06.2014
comment
Подпрограмма планирования возвращает NULL вместо действительного плана. Теперь вопрос почему?   -  person Vladimir F    schedule 18.06.2014
comment
О да, проверка возвращаемых значений :)   -  person Mali Remorker    schedule 18.06.2014

Ответы (1)


arrow_upward
4
arrow_downward

И ответ таков:

Для реальных преобразований mpi допустимы только две комбинации транспозиций и направлений:

  • реальное в сложное преобразование и FFTW_MPI_TRANSPOSED_OUT
  • сложное в реальное преобразование и FFTW_MPI_TRANSPOSED_IN

Я нашел это, копаясь внутри версии fftw3. 3.3.4 код, файл "rdft2-problem.c", комментарий к строке 120.

РЕДАКТИРОВАТЬ:

МИНИМАЛЬНЫЙ КОМПИЛЯЦИОННЫЙ И ПРИМЕР РАБОТЫ:

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(L/2+1,M, MPI_COMM_WORLD, &
       local_l, local_offset1)

  print *, id, "alloc complex=",alloc_local, local_l

  csrc = fftw_alloc_complex(alloc_local)
  call c_f_pointer(csrc, src, [M,local_l])

  !Caveat: Must partition the real storage according to complex layout, this is why
  ! I am using M and L/2+1 instead of M, 2*(L/2+1) as it was done in the original post
  alloc_local = fftw_mpi_local_size_2d(M,L/2+1, MPI_COMM_WORLD, &
      &                               local_M, local_offset2)

  print *, id, "alloc real=",alloc_local, local_m
  ! Two reals per complex
  ctgt = fftw_alloc_real(2*alloc_local)
  ! Only the first L are relevant, the rest is just dangling space (see fftw3 docs) 
  !caveat: since the padding is in the first index, the 2d data is laid out non-contiguously 
  !(L sensible reals, padding, padding, L sensible reals, padding, padding, ....)
  call c_f_pointer(ctgt, tgt, [2*(L/2+1),local_m])


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

  ! Should be non-null
  print *, 'plan:', plan

  src(3,2)=(1.,0)
  call fftw_mpi_execute_dft_c2r(plan, src, tgt) 

  call mpi_finalize(ierr)
end program thrashingfftw
person Mali Remorker    schedule 23.06.2014
comment
Пожалуйста, также примите ваш ответ, если вы не ожидаете лучшего в ближайшее время (я не жду), я думаю, что это уже возможно. - person Vladimir F; 24.06.2014
comment
Я был готов немного подождать, пока кто-нибудь не подтвердит ответ. Теперь, когда это произошло, я принимаю это. - person Mali Remorker; 24.06.2014
comment
Следующий вопрос неправильный вывод через fftw mpi r2c 2d и fftw mpi c2r 2d">stackoverflow.com/questions/41730864/ , этот код требует исправления для учета для fftw требуется дополнительное заполнение. - person francis; 20.01.2017
comment
Да, с тех пор я придумал для этого совершенно другую схему и всегда подозревал, что должен обновить ответ :). Кроме того, следует быть осторожным с тем, как глобальные измерения распределены по вычислительным узлам. Может возникнуть проблема с балансировкой нагрузки (которая иногда легко решается) - person Mali Remorker; 21.01.2017