Calling FFTW from Legacy Fortran
********************************

This chapter describes the interface to FFTW callable by Fortran code in older compilers not supporting the Fortran 2003 C interoperability features (see Calling FFTW from Modern Fortran).  This interface has the major disadvantage that it is not type-checked, so if you mistake the argument types or ordering then your program will not have any compiler errors, and will likely crash at runtime.  So, greater care is needed.  Also, technically interfacing older Fortran versions to C is nonstandard, but in practice we have found that the techniques used in this chapter have worked with all known Fortran compilers for many years.

The legacy Fortran interface differs from the C interface only in the prefix (dfftw_ instead of fftw_ in double precision) and a few other minor details.  This Fortran interface is included in the FFTW libraries by default, unless a Fortran compiler isn't found on your system or --disable-fortran is included in the configure flags.  We assume here that the reader is already familiar with the usage of FFTW in C, as described elsewhere in this manual.

The MPI parallel interface to FFTW is not currently available to legacy Fortran.

Fortran-interface routines
==========================

Nearly all of the FFTW functions have Fortran-callable equivalents. The name of the legacy Fortran routine is the same as that of the corresponding C routine, but with the fftw_ prefix replaced by dfftw_. (Technically, Fortran 77 identifiers are not allowed to have more than 6 characters, nor may they contain underscores.  Any compiler that enforces this limitation doesn't deserve to link to FFTW.)  The single and long-double precision versions use sfftw_ and lfftw_, respectively, instead of fftwf_ and fftwl_; quadruple precision (real*16) is available on some systems as fftwq_ (see Precision). (Note that long double on x86 hardware is usually at most 80-bit extended precision, not quadruple precision.)

For the most part, all of the arguments to the functions are the same, with the following exceptions:

  * plan variables (what would be of type fftw_plan in C), must be declared as a type that is at least as big as a pointer (address) on your machine.  We recommend using integer*8 everywhere, since this should always be big enough.

  * Any function that returns a value (e.g. fftw_plan_dft) is converted into a subroutine.  The return value is converted into an additional first parameter of this subroutine. (The reason for this is that some Fortran implementations seem to have trouble with C function return values, and vice versa.)

  * The Fortran routines expect multi-dimensional arrays to be in column-major order, which is the ordinary format of Fortran arrays (see Multi-dimensional Array Format).  They do this transparently and costlessly simply by reversing the order of the dimensions passed to FFTW, but this has one important consequence for multi-dimensional real-complex transforms, discussed below.

  * Wisdom import and export is somewhat more tricky because one cannot easily pass files or strings between C and Fortran; see Wisdom of Fortran?.

  * Legacy Fortran cannot use the fftw_malloc dynamic-allocation routine. If you want to exploit the SIMD FFTW (see SIMD alignment and fftw_malloc), you'll need to figure out some other way to ensure that your arrays are at least 16-byte aligned.

  * Since Fortran 77 does not have data structures, the fftw_iodim structure from the guru interface (see Guru vector and transform sizes) must be split into separate arguments.  In particular, any fftw_iodim array arguments in the C guru interface become three integer array arguments (n, is, and os) in the Fortran guru interface, all of whose lengths should be equal to the corresponding rank argument.

  * The guru planner interface in Fortran does not do any automatic translation between column-major and row-major; you are responsible for setting the strides etcetera to correspond to your Fortran arrays. However, as a slight bug that we are preserving for backwards compatibility, the plan_guru_r2r in Fortran does reverse the order of its kind array parameter, so the kind array of that routine should be in the reverse of the order of the iodim arrays (see above).

In general, you should take care to use Fortran data types that correspond to (i.e. are the same size as) the C types used by FFTW. In practice, this correspondence is usually straightforward (i.e. integer corresponds to int, real corresponds to float, etcetera).  The native Fortran double/single-precision complex type should be compatible with fftw_complex/fftwf_complex.  Such simple correspondences are assumed in the examples below.

FFTW Constants in Fortran
=========================

When creating plans in FFTW, a number of constants are used to specify options, such as FFTW_MEASURE or FFTW_ESTIMATE.  The same constants must be used with the wrapper routines, but of course the C header files where the constants are defined can't be incorporated directly into Fortran code.

Instead, we have placed Fortran equivalents of the FFTW constant definitions in the file fftw3.f, which can be found in the same directory as fftw3.h.  If your Fortran compiler supports a preprocessor of some sort, you should be able to include or #include this file; otherwise, you can paste it directly into your code.

In C, you combine different flags (like FFTW_PRESERVE_INPUT and FFTW_MEASURE) using the | operator; in Fortran you should just use +.  (Take care not to add in the same flag more than once, though.  Alternatively, you can use the ior intrinsic function standardized in Fortran 95.)

FFTW 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 legacy Fortran is:

             call dfftw_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 dfftw_execute, the semantics of Fortran (unlike C) allow the compiler to assume that the input/output arrays are not changed by dfftw_execute.  As a consequence, certain compilers end up optimizing out or repositioning the call to dfftw_execute, assuming incorrectly that it does nothing.

There are various workarounds to this, but the safest and simplest thing is to not use dfftw_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 dfftw_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 dfftw_execute_dft, Real-input (r2c) DFT plans should use use dfftw_execute_dft_r2c, and real-output (c2r) DFT plans should use dfftw_execute_dft_c2r.  The various r2r plans should use dfftw_execute_r2r.

  * 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 difficult of these, in Fortran, is the requirement that the new arrays have the same alignment as the original arrays, because there seems to be no way in legacy Fortran to obtain guaranteed-aligned arrays (analogous to fftw_malloc in C).  You can, of course, 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).

Fortran Examples
================

In C, you might have something like the following to transform a one-dimensional complex array:

             fftw_complex in[N], out[N];
             fftw_plan plan;

             plan = fftw_plan_dft_1d(N,in,out,FFTW_FORWARD,FFTW_ESTIMATE);
             fftw_execute(plan);
             fftw_destroy_plan(plan);

In Fortran, you would use the following to accomplish the same thing:

             double complex in, out
             dimension in(N), out(N)
             integer*8 plan

             call dfftw_plan_dft_1d(plan,N,in,out,FFTW_FORWARD,FFTW_ESTIMATE)
             call dfftw_execute_dft(plan, in, out)
             call dfftw_destroy_plan(plan)

Notice how all routines are called as Fortran subroutines, and the plan is returned via the first argument to dfftw_plan_dft_1d. Notice also that we changed fftw_execute to dfftw_execute_dft (see FFTW Execution in Fortran).  To do the same thing, but using 8 threads in parallel (see Multi-threaded FFTW), you would simply prefix these calls with:

             integer iret
             call dfftw_init_threads(iret)
             call dfftw_plan_with_nthreads(8)

(You might want to check the value of iret: if it is zero, it indicates an unlikely error during thread initialization.)

To transform a three-dimensional array in-place with C, you might do:

             fftw_complex arr[L][M][N];
             fftw_plan plan;

             plan = fftw_plan_dft_3d(L,M,N, arr,arr,
                                     FFTW_FORWARD, FFTW_ESTIMATE);
             fftw_execute(plan);
             fftw_destroy_plan(plan);

In Fortran, you would use this instead:

             double complex arr
             dimension arr(L,M,N)
             integer*8 plan

             call dfftw_plan_dft_3d(plan, L,M,N, arr,arr,
            &                       FFTW_FORWARD, FFTW_ESTIMATE)
             call dfftw_execute_dft(plan, arr, arr)
             call dfftw_destroy_plan(plan)

Note that we pass the array dimensions in the ``natural'' order in both C and Fortran.

To transform a one-dimensional real array in Fortran, you might do:

             double precision in
             dimension in(N)
             double complex out
             dimension out(N/2 + 1)
             integer*8 plan

             call dfftw_plan_dft_r2c_1d(plan,N,in,out,FFTW_ESTIMATE)
             call dfftw_execute_dft_r2c(plan, in, out)
             call dfftw_destroy_plan(plan)

To transform a two-dimensional real array, out of place, you might use the following:

             double precision in
             dimension in(M,N)
             double complex out
             dimension out(M/2 + 1, N)
             integer*8 plan

             call dfftw_plan_dft_r2c_2d(plan,M,N,in,out,FFTW_ESTIMATE)
             call dfftw_execute_dft_r2c(plan, in, out)
             call dfftw_destroy_plan(plan)

Important: Notice that it is the first dimension of the complex output array that is cut in half in Fortran, rather than the last dimension as in C.  This is a consequence of the interface routines reversing the order of the array dimensions passed to FFTW so that the Fortran program can use its ordinary column-major order.

Wisdom of Fortran?
==================

In this section, we discuss how one can import/export FFTW wisdom (saved plans) to/from a Fortran program; we assume that the reader is already familiar with wisdom, as described in Words of Wisdom-Saving Plans.

The basic problem is that is difficult to (portably) pass files and strings between Fortran and C, so we cannot provide a direct Fortran equivalent to the fftw_export_wisdom_to_file, etcetera, functions.  Fortran interfaces are provided for the functions that do not take file/string arguments, however: dfftw_import_system_wisdom, dfftw_import_wisdom, dfftw_export_wisdom, and dfftw_forget_wisdom.

So, for example, to import the system-wide wisdom, you would do:

             integer isuccess
             call dfftw_import_system_wisdom(isuccess)

As usual, the C return value is turned into a first parameter; isuccess is non-zero on success and zero on failure (e.g. if there is no system wisdom installed).

If you want to import/export wisdom from/to an arbitrary file or elsewhere, you can employ the generic dfftw_import_wisdom and dfftw_export_wisdom functions, for which you must supply a subroutine to read/write one character at a time.  The FFTW package contains an example file doc/f77_wisdom.f demonstrating how to implement import_wisdom_from_file and export_wisdom_to_file subroutines in this way.  (These routines cannot be compiled into the FFTW library itself, lest all FFTW-using programs be required to link with the Fortran I/O library.)