Calling FFTW from Modern Fortran
********************************

Fortran 2003 standardized ways for Fortran code to call C libraries, and this allows us to support a direct translation of the FFTW C API into Fortran.  Compared to the legacy Fortran 77 interface (see Calling FFTW from Legacy Fortran), this direct interface offers many advantages, especially compile-time type-checking and aligned memory allocation.  As of this writing, support for these C interoperability features seems widespread, having been implemented in nearly all major Fortran compilers (e.g. GNU, Intel, IBM, Oracle/Solaris, Portland Group, NAG).

This chapter documents that interface.  For the most part, since this interface allows Fortran to call the C interface directly, the usage is identical to C translated to Fortran syntax.  However, there are a few subtle points such as memory allocation, wisdom, and data types that deserve closer attention.

Overview of Fortran interface
=============================

FFTW provides a file fftw3.f03 that defines Fortran 2003 interfaces for all of its C routines, except for the MPI routines described elsewhere, which can be found in the same directory as fftw3.h (the C header file).  In any Fortran subroutine where you want to use FFTW functions, you should begin with:

       use, intrinsic :: iso_c_binding
       include 'fftw3.f03'

This includes the interface definitions and the standard iso_c_binding module (which defines the equivalents of C types).  You can also put the FFTW functions into a module if you prefer (see Defining an FFTW module).

At this point, you can now call anything in the FFTW C interface directly, almost exactly as in C other than minor changes in syntax. For example:

       type(C_PTR) :: plan
       complex(C_DOUBLE_COMPLEX), dimension(1024,1000) :: in, out
       plan = fftw_plan_dft_2d(1000,1024, in,out, FFTW_FORWARD,FFTW_ESTIMATE)
       ...
       call fftw_execute_dft(plan, in, out)
       ...
       call fftw_destroy_plan(plan)

A few important things to keep in mind are:

  * FFTW plans are type(C_PTR).  Other C types are mapped in the obvious way via the iso_c_binding standard: int turns into integer(C_INT), fftw_complex turns into complex(C_DOUBLE_COMPLEX), double turns into real(C_DOUBLE), and so on. See FFTW Fortran type reference.

  * Functions in C become functions in Fortran if they have a return value, and subroutines in Fortran otherwise.

  * The ordering of the Fortran array dimensions must be reversed when they are passed to the FFTW plan creation, thanks to differences in array indexing conventions (see Multi-dimensional Array Format).  This is unlike the legacy Fortran interface (see Fortran-interface routines), which reversed the dimensions for you.  See Reversing array dimensions.

  * Using ordinary Fortran array declarations like this works, but may yield suboptimal performance because the data may not be not aligned to exploit SIMD instructions on modern proessors (see SIMD alignment and fftw_malloc). Better performance will often be obtained by allocating with fftw_alloc. See Allocating aligned memory in Fortran.

  * Similar to the legacy Fortran interface (see FFTW Execution in Fortran), we currently recommend not using fftw_execute but rather using the more specialized functions like fftw_execute_dft (see New-array Execute Functions). However, you should execute the plan on the same arrays as the ones for which you created the plan, unless you are especially careful.  See Plan execution in Fortran.  To prevent you from using fftw_execute by mistake, the fftw3.f03 file does not provide an fftw_execute interface declaration.

  * Multiple planner flags are combined with ior (equivalent to | in C).  e.g. FFTW_MEASURE | FFTW_DESTROY_INPUT becomes ior(FFTW_MEASURE, FFTW_DESTROY_INPUT).  (You can also use + as long as you don't try to include a given flag more than once.)

Extended and quadruple precision in Fortran
-------------------------------------------

If FFTW is compiled in long double (extended) precision (see Installation and Customization), you may be able to call the resulting fftwl_ routines (see Precision) from Fortran if your compiler supports the C_LONG_DOUBLE_COMPLEX type code.

Because some Fortran compilers do not support C_LONG_DOUBLE_COMPLEX, the fftwl_ declarations are segregated into a separate interface file fftw3l.f03, which you should include in addition to fftw3.f03 (which declares precision-independent FFTW_ constants):

       use, intrinsic :: iso_c_binding
       include 'fftw3.f03'
       include 'fftw3l.f03'

We also support using the nonstandard __float128 quadruple-precision type provided by recent versions of gcc on 32- and 64-bit x86 hardware (see Installation and Customization), using the corresponding real(16) and complex(16) types supported by gfortran.  The quadruple-precision fftwq_ functions (see Precision) are declared in a fftw3q.f03 interface file, which should be included in addition to fftw3l.f03, as above.  You should also link with -lfftw3q -lquadmath -lm as in C.

Reversing array dimensions
==========================

A minor annoyance in calling FFTW from Fortran is that FFTW's array dimensions are defined in the C convention (row-major order), while Fortran's array dimensions are the opposite convention (column-major order). See Multi-dimensional Array Format.  This is just a bookkeeping difference, with no effect on performance.  The only consequence of this is that, whenever you create an FFTW plan for a multi-dimensional transform, you must always reverse the ordering of the dimensions.

For example, consider the three-dimensional (L,M,N) arrays:

       complex(C_DOUBLE_COMPLEX), dimension(L,M,N) :: in, out

To plan a DFT for these arrays using fftw_plan_dft_3d, you could do:

       plan = fftw_plan_dft_3d(N,M,L, in,out, FFTW_FORWARD,FFTW_ESTIMATE)

That is, from FFTW's perspective this is a N,M,L array. No data transposition need occur, as this is only notation.  Similarly, to use the more generic routine fftw_plan_dft with the same arrays, you could do:

       integer(C_INT), dimension(3) :: n = [N,M,L]
       plan = fftw_plan_dft_3d(3, n, in,out, FFTW_FORWARD,FFTW_ESTIMATE)

Note, by the way, that this is different from the legacy Fortran interface (see Fortran-interface routines), which automatically reverses the order of the array dimension for you.  Here, you are calling the C interface directly, so there is no ``translation'' layer.

An important thing to keep in mind is the implication of this for multidimensional real-to-complex transforms (see Multi-Dimensional DFTs of Real Data).  In C, a multidimensional real-to-complex DFT chops the last dimension roughly in half (N,M,L real input goes to N,M,L/2+1 complex output).  In Fortran, because the array dimension notation is reversed, the first dimension of the complex data is chopped roughly in half.  For example consider the r2c transform of L,M,N real input in Fortran:

       type(C_PTR) :: plan
       real(C_DOUBLE), dimension(L,M,N) :: in
       complex(C_DOUBLE_COMPLEX), dimension(L/2+1,M,N) :: out
       plan = fftw_plan_dft_r2c_3d(N,M,L, in,out, FFTW_ESTIMATE)
       ...
       call fftw_execute_dft_r2c(plan, in, out)

Alternatively, for an in-place r2c transform, as described in the C documentation we must pad the first dimension of the real input with an extra two entries (which are ignored by FFTW) so as to leave enough space for the complex output. The input is allocated as a 2[L/2+1],M,N array, even though only L,M,N of it is actually used.  In this example, we will allocate the array as a pointer type, using fftw_alloc to ensure aligned memory for maximum performance (see Allocating aligned memory in Fortran); this also makes it easy to reference the same memory as both a real array and a complex array.

       real(C_DOUBLE), pointer :: in(:,:,:)
       complex(C_DOUBLE_COMPLEX), pointer :: out(:,:,:)
       type(C_PTR) :: plan, data
       data = fftw_alloc_complex(int((L/2+1) * M * N, C_SIZE_T))
       call c_f_pointer(data, in, [2*(L/2+1),M,N])
       call c_f_pointer(data, out, [L/2+1,M,N])
       plan = fftw_plan_dft_r2c_3d(N,M,L, in,out, FFTW_ESTIMATE)
       ...
       call fftw_execute_dft_r2c(plan, in, out)
       ...
       call fftw_destroy_plan(plan)
       call fftw_free(data)

FFTW Fortran type reference
===========================

The following are the most important type correspondences between the C interface and Fortran:

  * Plans (fftw_plan and variants) are type(C_PTR) (i.e. an opaque pointer).

  * The C floating-point types double, float, and long double correspond to real(C_DOUBLE), real(C_FLOAT), and real(C_LONG_DOUBLE), respectively.  The C complex types fftw_complex, fftwf_complex, and fftwl_complex correspond in Fortran to complex(C_DOUBLE_COMPLEX), complex(C_FLOAT_COMPLEX), and complex(C_LONG_DOUBLE_COMPLEX), respectively. Just as in C (see Precision), the FFTW subroutines and types are prefixed with fftw_, fftwf_, and fftwl_ for the different precisions, and link to different libraries (-lfftw3, -lfftw3f, and -lfftw3l on Unix), but use the same include file fftw3.f03 and the same constants (all of which begin with FFTW_).  The exception is long double precision, for which you should also include fftw3l.f03 (see Extended and quadruple precision in Fortran).

  * The C integer types int and unsigned (used for planner flags) become integer(C_INT).  The C integer type ptrdiff_t (e.g. in the 64-bit Guru Interface) becomes integer(C_INTPTR_T), and size_t (in fftw_malloc etc.) becomes integer(C_SIZE_T).

  * The fftw_r2r_kind type (see Real-to-Real Transform Kinds) becomes integer(C_FFTW_R2R_KIND).  The various constant values of the C enumerated type (FFTW_R2HC etc.) become simply integer constants of the same names in Fortran.

  * Numeric array pointer arguments (e.g. double *) become dimension(*), intent(out) arrays of the same type, or dimension(*), intent(in) if they are pointers to constant data (e.g. const int *).  There are a few exceptions where numeric pointers refer to scalar outputs (e.g. for fftw_flops), in which case they are intent(out) scalar arguments in Fortran too. For the new-array execute functions (see New-array Execute Functions), the input arrays are declared dimension(*), intent(inout), since they can be modified in the case of in-place or FFTW_DESTROY_INPUT transforms.

  * Pointer return values (e.g double *) become type(C_PTR).  (If they are pointers to arrays, as for fftw_alloc_real, you can convert them back to Fortran array pointers with the standard intrinsic function c_f_pointer.)

  * The fftw_iodim type in the guru interface (see Guru vector and transform sizes) becomes type(fftw_iodim) in Fortran, a derived data type (the Fortran analogue of C's struct) with three integer(C_INT) components: n, is, and os, with the same meanings as in C.  The fftw_iodim64 type in the 64-bit guru interface (see 64-bit Guru Interface) is the same, except that its components are of type integer(C_INTPTR_T).

  * Using the wisdom import/export functions from Fortran is a bit tricky, and is discussed in Accessing the wisdom API from Fortran.  In brief, the FILE * arguments map to type(C_PTR), const char * to character(C_CHAR), dimension(*), intent(in) (null-terminated!), and the generic read-char/write-char functions map to type(C_FUNPTR).

You may be wondering if you need to search-and-replace real(kind(0.0d0)) (or whatever your favorite Fortran spelling of ``double precision'' is) with real(C_DOUBLE) everywhere in your program, and similarly for complex and integer types.  The answer is no; you can still use your existing types.  As long as these types match their C counterparts, things should work without a hitch.  The worst that can happen, e.g. in the (unlikely) event of a system where real(kind(0.0d0)) is different from real(C_DOUBLE), is that the compiler will give you a type-mismatch error.  That is, if you don't use the iso_c_binding kinds you need to accept at least the theoretical possibility of having to change your code in response to compiler errors on some future machine, but you don't need to worry about silently compiling incorrect code that yields runtime errors.

Plan execution in Fortran
=========================

In C, in order to use a plan, one normally calls fftw_execute, which executes the plan to perform the transform on the input/output arrays passed when the plan was created (see Using Plans).  The corresponding subroutine call in modern Fortran is:

      call fftw_execute(plan)

However, we have had reports that this causes problems with some recent optimizing Fortran compilers.  The problem is, because the input/output arrays are not passed as explicit arguments to fftw_execute, the semantics of Fortran (unlike C) allow the compiler to assume that the input/output arrays are not changed by fftw_execute.  As a consequence, certain compilers end up repositioning the call to fftw_execute, assuming incorrectly that it does nothing to the arrays.

There are various workarounds to this, but the safest and simplest thing is to not use fftw_execute in Fortran.  Instead, use the functions described in New-array Execute Functions, which take the input/output arrays as explicit arguments.  For example, if the plan is for a complex-data DFT and was created for the arrays in and out, you would do:

      call fftw_execute_dft(plan, in, out)

There are a few things to be careful of, however:

  * You must use the correct type of execute function, matching the way the plan was created.  Complex DFT plans should use fftw_execute_dft, Real-input (r2c) DFT plans should use use fftw_execute_dft_r2c, and real-output (c2r) DFT plans should use fftw_execute_dft_c2r.  The various r2r plans should use fftw_execute_r2r.  Fortunately, if you use the wrong one you will get a compile-time type-mismatch error (unlike legacy Fortran).

  * You should normally pass the same input/output arrays that were used when creating the plan.  This is always safe.

  * If you pass different input/output arrays compared to those used when creating the plan, you must abide by all the restrictions of the new-array execute functions (see New-array Execute Functions).  The most tricky of these is the requirement that the new arrays have the same alignment as the original arrays; the best (and possibly only) way to guarantee this is to use the fftw_alloc functions to allocate your arrays (see Allocating aligned memory in Fortran). Alternatively, you can use the FFTW_UNALIGNED flag when creating the plan, in which case the plan does not depend on the alignment, but this may sacrifice substantial performance on architectures (like x86) with SIMD instructions (see SIMD alignment and fftw_malloc).

Allocating aligned memory in Fortran
====================================

In order to obtain maximum performance in FFTW, you should store your data in arrays that have been specially aligned in memory (see SIMD alignment and fftw_malloc).  Enforcing alignment also permits you to safely use the new-array execute functions (see New-array Execute Functions) to apply a given plan to more than one pair of in/out arrays.  Unfortunately, standard Fortran arrays do not provide any alignment guarantees.  The only way to allocate aligned memory in standard Fortran is to allocate it with an external C function, like the fftw_alloc_real and fftw_alloc_complex functions.  Fortunately, Fortran 2003 provides a simple way to associate such allocated memory with a standard Fortran array pointer that you can then use normally.

We therefore recommend allocating all your input/output arrays using the following technique:

 1. Declare a pointer, arr, to your array of the desired type and dimensions.  For example, real(C_DOUBLE), pointer :: a(:,:) for a 2d real array, or complex(C_DOUBLE_COMPLEX), pointer :: a(:,:,:) for a 3d complex array.

 2. The number of elements to allocate must be an integer(C_SIZE_T).  You can either declare a variable of this type, e.g. integer(C_SIZE_T) :: sz, to store the number of elements to allocate, or you can use the int(..., C_SIZE_T) intrinsic function. e.g. set sz = L * M * N or use int(L * M * N, C_SIZE_T) for an L,M,N array.

 3. Declare a type(C_PTR) :: p to hold the return value from FFTW's allocation routine.  Set p = fftw_alloc_real(sz) for a real array, or p = fftw_alloc_complex(sz) for a complex array.

 4. Associate your pointer arr with the allocated memory p using the standard c_f_pointer subroutine: call c_f_pointer(p, arr, [...dimensions...]), where [...dimensions...]) are an array of the dimensions of the array (in the usual Fortran order). e.g. call c_f_pointer(p, arr, [L,M,N]) for an L,M,N array.  (Alternatively, you can omit the dimensions argument if you specified the shape explicitly when declaring arr.)  You can now use arr as a usual multidimensional array.

 5. When you are done using the array, deallocate the memory by call fftw_free(p) on p.

For example, here is how we would allocate an L,M 2d real array:

       real(C_DOUBLE), pointer :: arr(:,:)
       type(C_PTR) :: p
       p = fftw_alloc_real(int(L * M, C_SIZE_T))
       call c_f_pointer(p, arr, [L,M])
       ...use arr and arr(i,j) as usual...
       call fftw_free(p)

and here is an L,M,N 3d complex array:

       complex(C_DOUBLE_COMPLEX), pointer :: arr(:,:,:)
       type(C_PTR) :: p
       p = fftw_alloc_complex(int(L * M * N, C_SIZE_T))
       call c_f_pointer(p, arr, [L,M,N])
       ...use arr and arr(i,j,k) as usual...
       call fftw_free(p)

See Reversing array dimensions for an example allocating a single array and associating both real and complex array pointers with it, for in-place real-to-complex transforms.

Accessing the wisdom API from Fortran
=====================================

As explained in Words of Wisdom-Saving Plans, FFTW provides a ``wisdom'' API for saving plans to disk so that they can be recreated quickly.  The C API for exporting (see Wisdom Export) and importing (see Wisdom Import) wisdom is somewhat tricky to use from Fortran, however, because of differences in file I/O and string types between C and Fortran.

Wisdom File Export/Import from Fortran
--------------------------------------

The easiest way to export and import wisdom is to do so using fftw_export_wisdom_to_filename and fftw_wisdom_from_filename.  The only trick is that these require you to pass a C string, which is an array of type CHARACTER(C_CHAR) that is terminated by C_NULL_CHAR. You can call them like this:

       integer(C_INT) :: ret
       ret = fftw_export_wisdom_to_filename(C_CHAR_'my_wisdom.dat' // C_NULL_CHAR)
       if (ret .eq. 0) stop 'error exporting wisdom to file'
       ret = fftw_import_wisdom_from_filename(C_CHAR_'my_wisdom.dat' // C_NULL_CHAR)
       if (ret .eq. 0) stop 'error importing wisdom from file'

Note that prepending C_CHAR_ is needed to specify that the literal string is of kind C_CHAR, and we null-terminate the string by appending // C_NULL_CHAR.  These functions return an integer(C_INT) (ret) which is 0 if an error occurred during export/import and nonzero otherwise.

It is also possible to use the lower-level routines fftw_export_wisdom_to_file and fftw_import_wisdom_from_file, which accept parameters of the C type FILE*, expressed in Fortran as type(C_PTR). However, you are then responsible for creating the FILE* yourself.  You can do this by using iso_c_binding to define Fortran intefaces for the C library functions fopen and fclose, which is a bit strange in Fortran but workable.

Wisdom String Export/Import from Fortran
----------------------------------------

Dealing with FFTW's C string export/import is a bit more painful.  In particular, the fftw_export_wisdom_to_string function requires you to deal with a dynamically allocated C string.  To get its length, you must define an interface to the C strlen function, and to deallocate it you must define an interface to C free:

       use, intrinsic :: iso_c_binding
       interface
         integer(C_INT) function strlen(s) bind(C, name='strlen')
           import
           type(C_PTR), value :: s
         end function strlen
         subroutine free(p) bind(C, name='free')
           import
           type(C_PTR), value :: p
         end subroutine free
       end interface

Given these definitions, you can then export wisdom to a Fortran character array:

       character(C_CHAR), pointer :: s(:)
       integer(C_SIZE_T) :: slen
       type(C_PTR) :: p
       p = fftw_export_wisdom_to_string()
       if (.not. c_associated(p)) stop 'error exporting wisdom'
       slen = strlen(p)
       call c_f_pointer(p, s, [slen+1])
       ...
       call free(p)

Note that slen is the length of the C string, but the length of the array is slen+1 because it includes the terminating null character.  (You can omit the +1 if you don't want Fortran to know about the null character.) The standard c_associated function checks whether p is a null pointer, which is returned by fftw_export_wisdom_to_string if there was an error.

To import wisdom from a string, use fftw_import_wisdom_from_string as usual; note that the argument of this function must be a character(C_CHAR) that is terminated by the C_NULL_CHAR character, like the s array above.

Wisdom Generic Export/Import from Fortran
-----------------------------------------

The most generic wisdom export/import functions allow you to provide an arbitrary callback function to read/write one character at a time in any way you want.  However, your callback function must be written in a special way, using the bind(C) attribute to be passed to a C interface.

In particular, to call the generic wisdom export function fftw_export_wisdom, you would write a callback subroutine of the form:

       subroutine my_write_char(c, p) bind(C)
         use, intrinsic :: iso_c_binding
         character(C_CHAR), value :: c
         type(C_PTR), value :: p
         ...write c...
       end subroutine my_write_char

Given such a subroutine (along with the corresponding interface definition), you could then export wisdom using:

       call fftw_export_wisdom(c_funloc(my_write_char), p)

The standard c_funloc intrinsic converts a Fortran bind(C) subroutine into a C function pointer.  The parameter p is a type(C_PTR) to any arbitrary data that you want to pass to my_write_char (or C_NULL_PTR if none).  (Note that you can get a C pointer to Fortran data using the intrinsic c_loc, and convert it back to a Fortran pointer in my_write_char using c_f_pointer.)

Similarly, to use the generic fftw_import_wisdom, you would define a callback function of the form:

       integer(C_INT) function my_read_char(p) bind(C)
         use, intrinsic :: iso_c_binding
         type(C_PTR), value :: p
         character :: c
         ...read a character c...
         my_read_char = ichar(c, C_INT)
       end function my_read_char

       ....

       integer(C_INT) :: ret
       ret = fftw_import_wisdom(c_funloc(my_read_char), p)
       if (ret .eq. 0) stop 'error importing wisdom'

Your function can return -1 if the end of the input is reached. Again, p is an arbitrary type(C_PTR that is passed through to your function.  fftw_import_wisdom returns 0 if an error occurred and nonzero otherwise.

Defining an FFTW module
=======================

Rather than using the include statement to include the fftw3.f03 interface file in any subroutine where you want to use FFTW, you might prefer to define an FFTW Fortran module.  FFTW does not install itself as a module, primarily because fftw3.f03 can be shared between different Fortran compilers while modules (in general) cannot.  However, it is trivial to define your own FFTW module if you want.  Just create a file containing:

       module FFTW3
         use, intrinsic :: iso_c_binding
         include 'fftw3.f03'
       end module

Compile this file into a module as usual for your compiler (e.g. with gfortran -c you will get a file fftw3.mod).  Now, instead of include 'fftw3.f03', whenever you want to use FFTW routines you can just do:

       use FFTW3

as usual for Fortran modules.  (You still need to link to the FFTW library, of course.)