FFTW3 real-to-complex and inverse transforms

[Navdeep_Rana]

I would like to perform inplace real-to-complex and inverse transforms! The function has a signature

fftw_plan fftw_mpi_plan_dft_r2c_2d(ptrdiff_t n0, ptrdiff_t n1,
double *in, fftw_complex *out,
MPI_Comm comm, unsigned flags);

For an inplace complex-to-complex transform *in,*out are pointers of same type. However, here the pointer type differs. In C you can possibly typecast

a double pointer to a fftw_complex. I am confused how to do it in fortran? Any suggestions?

[Steve_Lionel]

From what I read, fftw_complex is the same as COMPLEX(C_DOUBLE) (or COMPLEX(8)) in Intel Fortran. It is not a double pointer.  http://www.fftw.org/doc/Complex-numbers.html

Note that COMPLEX(8) is not the same as COMPLEX*8!

[Navdeep_Rana]

Steve Lionel (Ret.) wrote:

From what I read, fftw_complex is the same as COMPLEX(C_DOUBLE) (or COMPLEX(8)) in Intel Fortran. It is not a double pointer.  http://www.fftw.org/doc/Complex-numbers.html

Note that COMPLEX(8) is not the same as COMPLEX*8!

This I understand but an inplace transform means the input (double *in) and output (fftw_complex *out) are the same. In C you can do (fftw_complex *) double_pointer to achieve it. How do I do it in fortran?

[Steve_Lionel]

Why not just pass the same variable as both arguments? Unless my understanding of C is much worse than I think it is (always a possibility), you're not being asked to pass a pointer to a pointer, unless I am seriously misreading http://www.fftw.org/doc/Real_002ddata-DFTs.html

The problem you'll run into in Fortran is type checking if you have declared an explicit interface for the FFTW routine that has the arguments of different types. The easy solution is to give both arguments the same type in your interface. A more complicated solution would be to declare a pointer to an array of double and use C_LOC and C_F_POINTER to "cast" the complex to double (or vice-versa.) The easiest solution is to not use an explicit interface at all, though I find that distasteful.

[Navdeep_Rana]

fftw_complex is defined as

typedef double fftw_complex[2]

so It is not a pointer to a pointer as you mentioned.

Yes I have an explicit interface and I plan to keep it. As I understand it, I have to define two pointers

  real(C_DOUBLE), pointer, dimension(:,:,:) :: u_r
  complex(C_DOUBLE_COMPLEX), pointer, dimension(:,:,:) :: u_c

one for real and other for complex data. I allocate my data using fftw_alloc_real i.e.

type(C_PTR) :: ptr_u
ptr_u = fftw_alloc_real(2*alloc_size)

and then typecast as

call c_f_pointer(ptr_u,u_r,[2*(N(1)/2+1),N(2),local_n3])
call c_f_pointer(ptr_u,u_c,[(N(1)/2+1),N(2),local_n3])

the shapes come out correctly, which I checked but when I go on creating a plan i.e

pfor = fftw_mpi_plan_dft_r2c_3d(N(3),N(2),N(1),u_r,u_c,MPI_COMM_WORLD,FFTW_FORWARD,FFTW_ESTIMATE)

It spits out an error i.e.

Error: More actual than formal arguments in procedure call at (1)

Any ideas on this?

[Steve_Lionel]

That's not an ifort error (gfortran?), and I don't know what your interface looks like. It does seem a simple error to resolve if you compare the call with the interface.

[FortranFan]

Navdeep R. wrote:

.. Any ideas on this?

Assuming the function prototype is

fftw_plan fftw_mpi_plan_dft_r2c_3d(ptrdiff_t n0, ptrdiff_t n1, ptrdiff_t n2,
                                   double *in, fftw_complex *out,
                                   MPI_Comm comm, unsigned flags);

and with typedef  fftw_complex as you indicate, your code is inconsistent when it comes to 5th parameter in this function i.e., with respect to fftw_complex *out:

• complex(C_DOUBLE_COMPLEX), pointer is not the same as fftw_complex * because
• fftw_complex * essentially reduces to double * given the above typedef

You can just have both your u_r and u_c objects in Fortran declared as real(C_DOUBLE), pointer, dimension(:,:,:) and i think your code will work ok.

[Navdeep_Rana]

FortranFan wrote:

Quote:

You can just have both your u_r and u_c objects in Fortran declared as real(C_DOUBLE), pointer, dimension(:,:,:) and i think your code will work ok.

That's the first thing I tried, but it gives a type mismatch error.

Error: Type mismatch in argument 'out' at (1); passed REAL(8) to COMPLEX(8)

Hence the question!

[Steve_Lionel]

I'm just curious why you are asking in this forum when you're not using ifort, not that we're unwilling to help.

Please show your explicit interface.

[Navdeep_Rana]

I have a student license for ifort and related tools, which is installed on my main machine.

For the time being I am on my laptop which only has gfortran.

I inherited a pseudo-spectral CFD code which uses FFTW2 interface, and I would like to switch to FFTW3

program fftw3_demo
  use mpi
  use, intrinsic :: iso_c_binding
  implicit none
  include "fftw3-mpi.f03"

  integer(C_INTPTR_T), dimension(3) :: N = [32,32,32]
  real(C_DOUBLE), pointer, dimension(:,:,:) :: u_r
  complex(C_DOUBLE_COMPLEX), pointer, dimension(:,:,:) :: u_c
  integer(C_INTPTR_T) :: alloc_size,local_n3,local_n3_start
  type(C_PTR) :: ptr_u,pfor,pinv

  integer :: mpi_error_flag,mpi_rank,mpi_procs
  real(kind=8) :: u_sum_old
  
  call mpi_init(mpi_error_flag)
  call mpi_comm_rank(MPI_COMM_WORLD, mpi_rank,mpi_error_flag)
  call mpi_comm_size(MPI_COMM_WORLD, mpi_procs, mpi_error_flag)

  call fftw_mpi_init()

  alloc_size = fftw_mpi_local_size_3d(N(3),N(2),N(1)/2+1,MPI_COMM_WORLD,local_n3,local_n3_start)
  ptr_u = fftw_alloc_real(2*alloc_size)
  
  call c_f_pointer(ptr_u,u_r,[2*(N(1)/2+1),N(2),local_n3])
  call c_f_pointer(ptr_u,u_c,[(N(1)/2+1),N(2),local_n3])
  
  write(*,*) shape(u_r), shape(u_c)
  !! compiles and runs correctly till here.
  pfor = fftw_mpi_plan_dft_r2c_3d(N(3),N(2),N(1),u_r,u_c,MPI_COMM_WORLD,FFTW_FORWARD,FFTW_ESTIMATE)
  pinv = fftw_mpi_plan_dft_c2r_3d(N(3),N(2),N(1),u_c,u_r,MPI_COMM_WORLD,FFTW_BACKWARD,FFTW_ESTIMATE)

  !!u_r(:,:,:) = 0.d0 
  !!u_r(:N(1),:,:) = 1.d0 
  !!u_sum_old = sum(u_r(:N(1),:,:))
  !!call fftw_execute(pfor)
  !!call fftw_execute(pinv)

  !!write(*,*) u_sum_old,sum(u_r(:N(1),:,:))/(N(1)*N(2)*N(3))
  !!call fftw_free(ptr_u)
  call fftw_mpi_cleanup()
  call mpi_finalize(mpi_error_flag)
  
end program fftw3_demo

.

About the interface, right now it's just single program.

[FortranFan · Accepted solution]

Navdeep R. wrote:

..

program fftw3_demo
  ..
  !! compiles and runs correctly till here.
  pfor = fftw_mpi_plan_dft_r2c_3d(N(3),N(2),N(1),u_r,u_c,MPI_COMM_WORLD,FFTW_FORWARD,FFTW_ESTIMATE)
..

According to this site, the interface for fftw_mpi_plan_dft_r2c_3d is

    type(C_PTR) function fftw_mpi_plan_dft_r2c_3d(n0,n1,n2,in,out,comm,flags) bind(C, name='fftw_mpi_plan_dft_r2c_3d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INTPTR_T), value :: n2
      real(C_DOUBLE), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_MPI_FINT), value :: comm
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_dft_r2c_3d

which has 7 dummy arguments, however your code has 8 which would explain, "Error: More actual than formal arguments in procedure call at (1)". from gfortran.

Is one of FFTW_FORWARD or FFTW_ESTIMATE superfluous?

[Navdeep_Rana]

FortranFan wrote:

According to this site, the interface for fftw_mpi_plan_dft_r2c_3d is

    type(C_PTR) function fftw_mpi_plan_dft_r2c_3d(n0,n1,n2,in,out,comm,flags) bind(C, name='fftw_mpi_plan_dft_r2c_3d_f03')
      import
      integer(C_INTPTR_T), value :: n0
      integer(C_INTPTR_T), value :: n1
      integer(C_INTPTR_T), value :: n2
      real(C_DOUBLE), dimension(*), intent(out) :: in
      complex(C_DOUBLE_COMPLEX), dimension(*), intent(out) :: out
      integer(C_MPI_FINT), value :: comm
      integer(C_INT), value :: flags
    end function fftw_mpi_plan_dft_r2c_3d

which has 7 dummy arguments, however your code has 8 which would explain, "Error: More actual than formal arguments in procedure call at (1)". from gfortran.

Is one of FFTW_FORWARD or FFTW_ESTIMATE superfluous?

Yes indeed. When I wrote the line, I pulled in FFTW_FORWARD/ BACKWARD from memory, which is not required here.

I will tweak the code and see if it works. Thanks a ton.

Update : It works. So the solution is to use two different pointers, one real and other complex.