Wie mache ich eine fftw3 MPI "transponierte" 2D Transformation wenn möglich?

Betrachten Sie eine 2D-Transformation der Form L x M (Spaltenhauptaufbau) aus einem komplexen Arraysrc zu einem realen Arraytgt. Oder in Fortranese,

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

Entsprechende Hinweise sind

type(C_PTR) :: csrc,ctgt   .

Ich würde sie folgendermaßen zuordnen:

  ! 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])

Der Plan würde nun wie folgt erstellt:

! 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))

Schließlich würde die Transformation durchgeführt werden als:

call fftw_mpi_execute_dft_c2r(plan, src, tgt)

Dieses Rezept funktioniert jedoch nicht. Der letzte Aufruf verursacht einen Segmentierungsfehler. Zuerst dachte ich, das könnte etwas damit zu tun haben, wie ich es zuordnetesrc undtgt Arrays, aber spielen mit unterschiedlicher Menge an Speicher zugewiesentgt gab kein Ergebnis. Also mache ich entweder etwas wirklich Dummes, oder das ist überhaupt nicht möglich.

EDIT: MINIMALISTISCHES KOMPILIERBARES BEISPIEL

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

Antworten auf die Frage(1)

Ihre Antwort auf die Frage