FFTW Fortran Real2Complex R2C output array structure? fftw_plan_dft_r2c_3d?

Hi,

I'm using the FFTW package and am using the fortan wrappers
to compute a 3D transform. This I seem to do successfully, however
I can not seem to disentangle the raw output.
>From studying the figure on page 8 of the fftw 3.1.2 manual
this is what I think is going on:

For the output complex array I do the following shuffle:

-------------------------------------------------

  subroutine FFTWshuffle_R2C_3D(array)

    integer :: ngrid1,ngrid2,ngrid3
    integer :: i,j,k,ii,jj,kk

    complex(DPC), dimension(:,:,:) :: array
    complex(DPC), allocatable, dimension(:,:,:) :: image

    ngrid1 = size(array,dim=1)
    ngrid2 = size(array,dim=2)
    ngrid3 = size(array,dim=3)

    allocate (image(ngrid1,ngrid2,ngrid3))

    do j=1,ngrid2
    do i=1,ngrid1

       ii = i + ngrid1/2
       jj = j + ngrid2/2

       if (ii.gt.ngrid1) ii=ii-ngrid1
       if (jj.gt.ngrid2) jj=jj-ngrid2

       image(i,j,:)=array(ii,jj,:)

    enddo
    enddo

    array=image

    deallocate (image)

  end subroutine FFTWshuffle_R2C_3D

------------------------------------------

(I don't shuffle the 3rd dimension since this runs from 1---ngrid3/2+1
and already
contains just the +ve frequncies.)

---------This does not give me the correct answer,
so any help as to what is going wrong would be welcome.

Thanks
Robert

On Jun 12, 1:16 pm, "Robert S." <r...@astro.upenn.edu> wrote:
> I'm using the FFTW package and am using the fortan wrappers
> to compute a 3D transform. This I seem to do successfully, however
> I can not seem to disentangle the raw output.>From studying the figure on page 8 of the fftw 3.1.2 manual
>
> [...]
>
> (I don't shuffle the 3rd dimension since this runs from 1---ngrid3/2+1
> and already
> contains just the +ve frequncies.)

Remember that a Fortran array corresponds to a C array with *reversed*
dimensions (column-major vs. row-major format).  So, in the C FFTW
interface, for an r2c (real-input) transform the last dimension is cut
in ~half, but in the Fortran interface the *first* dimension is cut in
~half.

This is discussed in the FFTW manual.  See the third bullet point of:
   http://www.fftw.org/doc/Fortran_002dinterface-routines.html
and the last example of:
    http://www.fftw.org/doc/Fortran-Examples.html

Regards,
Steven G. Johnson

On Jun 13, 5:40 am, "Steven G. Johnson" <stev...@alum.mit.edu> wrote:
> On Jun 12, 1:16 pm, "Robert S." <r...@astro.upenn.edu> wrote:
>
> > I'm using the FFTW package and am using the fortan wrappers
> > to compute a 3D transform. This I seem to do successfully, however
> > I can not seem to disentangle the raw output.>From studying the figure on page 8 of the fftw 3.1.2 manual
>
> > [...]
>
> > (I don't shuffle the 3rd dimension since this runs from 1---ngrid3/2+1
> > and already
> > contains just the +ve frequncies.)
>
> Remember that a Fortran array corresponds to a C array with *reversed*
> dimensions (column-major vs. row-major format).  So, in the C FFTW
> interface, for an r2c (real-input) transform the last dimension is cut
> in ~half, but in the Fortran interface the *first* dimension is cut in
> ~half.
>
> This is discussed in the FFTW manual.  See the third bullet point of:
>    http://www.fftw.org/doc/Fortran_002dinterface-routines.html
> and the last example of:
>    http://www.fftw.org/doc/Fortran-Examples.html
>
> Regards,
> Steven G. Johnson


Dear Steven,
Many thanks!  The routine works pefectly.
Best,
Robert

AND just for completeness:
The working frequency shuffle routine, based on Stephen's comments, is
presented below.
This places the zero frequency value at the position  array(1, ng2/2 ,
ng3/2 )

-----------------------------------------------------------------------------------------------

  subroutine FFTWshuffle_R2C_3D(array)

    integer :: ng1,ng2,ng3
    integer :: j,k,jj,kk

    complex(DPC), dimension(:,:,:) :: array
    complex(DPC), allocatable, dimension(:,:,:) :: image

    ng1 = size(array,dim=1)
    ng2 = size(array,dim=2)
    ng3 = size(array,dim=3)

    allocate (image(ng1,ng2,ng3))

    do k=1,ng3
    do j=1,ng2

       jj = j + ng2/2 + 1
       kk = k + ng3/2 + 1

       if (jj.gt.ng2) jj=jj-ng2
       if (kk.gt.ng3) kk=kk-ng3

       image(:,j,k)=array(:,jj,kk)

    enddo
    enddo

    array=image

    deallocate (image)

  end subroutine FFTWshuffle_R2C_3D

-------------------------------------------------------------------