2

Multithreaded FFTW can be implemented as in this from FFTW homepage. Instead, we want to call the serial FFTW routines within OpenMP parallel environment using multiple threads. We want variable and fourier_variable to be thread-safe. This could be done by using PURE subroutines and declaring variable and fourier_variable inside it. The question here is related to calling FFTW routines like fftw_execute_dft_r2c from within a PURE subroutine.

A stripped-down version of the code is presented here just for reference (the full code is an optimisation solver involving many FFTW calls).

PROGRAM main
USE fft_call
REAL(8), DIMENSION(1:N, 1:N) :: variable
COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(1:N/2+1, 1:N) :: fourier_variable
INTEGER :: JJ

!$OMP PARALLEL
!$OMP DO 
 DO JJ = 1, 5
   call fourier_to_physical(fourier_variable, variable)
 END DO
!$OMP END DO
!$OMP END PARALLEL
END PROGRAM main


MODULE fft_call

contains
PURE SUBROUTINE fourier_to_physical( fourier_variable, variable)
IMPLICIT NONE
REAL(8), DIMENSION(1:N, 1:N) :: variable
COMPLEX(C_DOUBLE_COMPLEX), DIMENSION(1:N/2+1, 1:N), INTENT(OUT) :: fourier_variable
    CALL fftw_execute_dft_r2c (plan_fwd, variable, fourier_variable)
END SUBROUTINE fourier_to_physical

END MODULE fft_call

The error while calling fftw_plan_dft_r2c_2d from the PURE subroutine fourier_to_physical:

Error: Function reference to 'fftw_plan_dft_r2c' at (1) is to a non-PURE procedure within a PURE procedure

The question: is there a way to call FFTW routines like fftw_execute_dft_r2c from within a PURE subroutine in Fortran90?

Or, in other words, are their PURE versions of fftw_execute_dft_r2c such that we can call them from PURE procedures? We are beginners to OpenMP.

RTh
  • 107
  • 1
  • 1
  • 4
  • 1
    No. FFTW routines are not pure (in the Fortran sense) so can not be called from a pure Fortran subprogram. But why do you need pure? I can't see anything in what you are doing that requires it. Oh, and please don't use the non-portable Real(8), see any number of questions here which discuss that. – Ian Bush Aug 26 '20 at 14:46
  • As already mentioned above, we need `variable` and `fourier_variable` to be thread-safe. – RTh Aug 26 '20 at 14:52
  • 3
    Thread safety does not require called routines to be pure. I agree it is very nice to do so, but the use of pure is severely limited by use of external libraries such as FFTW, BLAS or MPI, few of which provide the appropriate pure interfaces. Thus in practice more often than not you just have to forget pure. If properly ceded you don't need to declare the routine pure to make it work. – Ian Bush Aug 26 '20 at 15:03
  • Thanks. We will work on improving our knowledge about thread-safety. – RTh Aug 26 '20 at 15:08

0 Answers0