cannam@95: @node Calling FFTW from Modern Fortran, Calling FFTW from Legacy Fortran, Distributed-memory FFTW with MPI, Top cannam@95: @chapter Calling FFTW from Modern Fortran cannam@95: @cindex Fortran interface cannam@95: cannam@95: Fortran 2003 standardized ways for Fortran code to call C libraries, cannam@95: and this allows us to support a direct translation of the FFTW C API cannam@95: into Fortran. Compared to the legacy Fortran 77 interface cannam@95: (@pxref{Calling FFTW from Legacy Fortran}), this direct interface cannam@95: offers many advantages, especially compile-time type-checking and cannam@95: aligned memory allocation. As of this writing, support for these C cannam@95: interoperability features seems widespread, having been implemented in cannam@95: nearly all major Fortran compilers (e.g. GNU, Intel, IBM, cannam@95: Oracle/Solaris, Portland Group, NAG). cannam@95: @cindex portability cannam@95: cannam@95: This chapter documents that interface. For the most part, since this cannam@95: interface allows Fortran to call the C interface directly, the usage cannam@95: is identical to C translated to Fortran syntax. However, there are a cannam@95: few subtle points such as memory allocation, wisdom, and data types cannam@95: that deserve closer attention. cannam@95: cannam@95: @menu cannam@95: * Overview of Fortran interface:: cannam@95: * Reversing array dimensions:: cannam@95: * FFTW Fortran type reference:: cannam@95: * Plan execution in Fortran:: cannam@95: * Allocating aligned memory in Fortran:: cannam@95: * Accessing the wisdom API from Fortran:: cannam@95: * Defining an FFTW module:: cannam@95: @end menu cannam@95: cannam@95: @c ------------------------------------------------------- cannam@95: @node Overview of Fortran interface, Reversing array dimensions, Calling FFTW from Modern Fortran, Calling FFTW from Modern Fortran cannam@95: @section Overview of Fortran interface cannam@95: cannam@95: FFTW provides a file @code{fftw3.f03} that defines Fortran 2003 cannam@95: interfaces for all of its C routines, except for the MPI routines cannam@95: described elsewhere, which can be found in the same directory as cannam@95: @code{fftw3.h} (the C header file). In any Fortran subroutine where cannam@95: you want to use FFTW functions, you should begin with: cannam@95: cannam@95: @cindex iso_c_binding cannam@95: @example cannam@95: use, intrinsic :: iso_c_binding cannam@95: include 'fftw3.f03' cannam@95: @end example cannam@95: cannam@95: This includes the interface definitions and the standard cannam@95: @code{iso_c_binding} module (which defines the equivalents of C cannam@95: types). You can also put the FFTW functions into a module if you cannam@95: prefer (@pxref{Defining an FFTW module}). cannam@95: cannam@95: At this point, you can now call anything in the FFTW C interface cannam@95: directly, almost exactly as in C other than minor changes in syntax. cannam@95: For example: cannam@95: cannam@95: @findex fftw_plan_dft_2d cannam@95: @findex fftw_execute_dft cannam@95: @findex fftw_destroy_plan cannam@95: @example cannam@95: type(C_PTR) :: plan cannam@95: complex(C_DOUBLE_COMPLEX), dimension(1024,1000) :: in, out cannam@95: plan = fftw_plan_dft_2d(1000,1024, in,out, FFTW_FORWARD,FFTW_ESTIMATE) cannam@95: ... cannam@95: call fftw_execute_dft(plan, in, out) cannam@95: ... cannam@95: call fftw_destroy_plan(plan) cannam@95: @end example cannam@95: cannam@95: A few important things to keep in mind are: cannam@95: cannam@95: @itemize @bullet cannam@95: cannam@95: @item cannam@95: @tindex fftw_complex cannam@95: @ctindex C_PTR cannam@95: @ctindex C_INT cannam@95: @ctindex C_DOUBLE cannam@95: @ctindex C_DOUBLE_COMPLEX cannam@95: FFTW plans are @code{type(C_PTR)}. Other C types are mapped in the cannam@95: obvious way via the @code{iso_c_binding} standard: @code{int} turns cannam@95: into @code{integer(C_INT)}, @code{fftw_complex} turns into cannam@95: @code{complex(C_DOUBLE_COMPLEX)}, @code{double} turns into cannam@95: @code{real(C_DOUBLE)}, and so on. @xref{FFTW Fortran type reference}. cannam@95: cannam@95: @item cannam@95: Functions in C become functions in Fortran if they have a return value, cannam@95: and subroutines in Fortran otherwise. cannam@95: cannam@95: @item cannam@95: The ordering of the Fortran array dimensions must be @emph{reversed} cannam@95: when they are passed to the FFTW plan creation, thanks to differences cannam@95: in array indexing conventions (@pxref{Multi-dimensional Array cannam@95: Format}). This is @emph{unlike} the legacy Fortran interface cannam@95: (@pxref{Fortran-interface routines}), which reversed the dimensions cannam@95: for you. @xref{Reversing array dimensions}. cannam@95: cannam@95: @item cannam@95: @cindex alignment cannam@95: @cindex SIMD cannam@95: Using ordinary Fortran array declarations like this works, but may cannam@95: yield suboptimal performance because the data may not be not aligned cannam@95: to exploit SIMD instructions on modern proessors (@pxref{SIMD cannam@95: alignment and fftw_malloc}). Better performance will often be obtained cannam@95: by allocating with @samp{fftw_alloc}. @xref{Allocating aligned memory cannam@95: in Fortran}. cannam@95: cannam@95: @item cannam@95: @findex fftw_execute cannam@95: Similar to the legacy Fortran interface (@pxref{FFTW Execution in cannam@95: Fortran}), we currently recommend @emph{not} using @code{fftw_execute} cannam@95: but rather using the more specialized functions like cannam@95: @code{fftw_execute_dft} (@pxref{New-array Execute Functions}). cannam@95: However, you should execute the plan on the @code{same arrays} as the cannam@95: ones for which you created the plan, unless you are especially cannam@95: careful. @xref{Plan execution in Fortran}. To prevent cannam@95: you from using @code{fftw_execute} by mistake, the @code{fftw3.f03} cannam@95: file does not provide an @code{fftw_execute} interface declaration. cannam@95: cannam@95: @item cannam@95: @cindex flags cannam@95: Multiple planner flags are combined with @code{ior} (equivalent to @samp{|} in C). e.g. @code{FFTW_MEASURE | FFTW_DESTROY_INPUT} becomes @code{ior(FFTW_MEASURE, FFTW_DESTROY_INPUT)}. (You can also use @samp{+} as long as you don't try to include a given flag more than once.) cannam@95: cannam@95: @end itemize cannam@95: cannam@95: @menu cannam@95: * Extended and quadruple precision in Fortran:: cannam@95: @end menu cannam@95: cannam@95: @node Extended and quadruple precision in Fortran, , Overview of Fortran interface, Overview of Fortran interface cannam@95: @subsection Extended and quadruple precision in Fortran cannam@95: @cindex precision cannam@95: cannam@95: If FFTW is compiled in @code{long double} (extended) precision cannam@95: (@pxref{Installation and Customization}), you may be able to call the cannam@95: resulting @code{fftwl_} routines (@pxref{Precision}) from Fortran if cannam@95: your compiler supports the @code{C_LONG_DOUBLE_COMPLEX} type code. cannam@95: cannam@95: Because some Fortran compilers do not support cannam@95: @code{C_LONG_DOUBLE_COMPLEX}, the @code{fftwl_} declarations are cannam@95: segregated into a separate interface file @code{fftw3l.f03}, which you cannam@95: should include @emph{in addition} to @code{fftw3.f03} (which declares cannam@95: precision-independent @samp{FFTW_} constants): cannam@95: cannam@95: @cindex iso_c_binding cannam@95: @example cannam@95: use, intrinsic :: iso_c_binding cannam@95: include 'fftw3.f03' cannam@95: include 'fftw3l.f03' cannam@95: @end example cannam@95: cannam@95: We also support using the nonstandard @code{__float128} cannam@95: quadruple-precision type provided by recent versions of @code{gcc} on cannam@95: 32- and 64-bit x86 hardware (@pxref{Installation and Customization}), cannam@95: using the corresponding @code{real(16)} and @code{complex(16)} types cannam@95: supported by @code{gfortran}. The quadruple-precision @samp{fftwq_} cannam@95: functions (@pxref{Precision}) are declared in a @code{fftw3q.f03} cannam@95: interface file, which should be included in addition to cannam@95: @code{fftw3l.f03}, as above. You should also link with cannam@95: @code{-lfftw3q -lquadmath -lm} as in C. cannam@95: cannam@95: @c ------------------------------------------------------- cannam@95: @node Reversing array dimensions, FFTW Fortran type reference, Overview of Fortran interface, Calling FFTW from Modern Fortran cannam@95: @section Reversing array dimensions cannam@95: cannam@95: @cindex row-major cannam@95: @cindex column-major cannam@95: A minor annoyance in calling FFTW from Fortran is that FFTW's array cannam@95: dimensions are defined in the C convention (row-major order), while cannam@95: Fortran's array dimensions are the opposite convention (column-major cannam@95: order). @xref{Multi-dimensional Array Format}. This is just a cannam@95: bookkeeping difference, with no effect on performance. The only cannam@95: consequence of this is that, whenever you create an FFTW plan for a cannam@95: multi-dimensional transform, you must always @emph{reverse the cannam@95: ordering of the dimensions}. cannam@95: cannam@95: For example, consider the three-dimensional (@threedims{L,M,N}) arrays: cannam@95: cannam@95: @example cannam@95: complex(C_DOUBLE_COMPLEX), dimension(L,M,N) :: in, out cannam@95: @end example cannam@95: cannam@95: To plan a DFT for these arrays using @code{fftw_plan_dft_3d}, you could do: cannam@95: cannam@95: @findex fftw_plan_dft_3d cannam@95: @example cannam@95: plan = fftw_plan_dft_3d(N,M,L, in,out, FFTW_FORWARD,FFTW_ESTIMATE) cannam@95: @end example cannam@95: cannam@95: That is, from FFTW's perspective this is a @threedims{N,M,L} array. cannam@95: @emph{No data transposition need occur}, as this is @emph{only cannam@95: notation}. Similarly, to use the more generic routine cannam@95: @code{fftw_plan_dft} with the same arrays, you could do: cannam@95: cannam@95: @example cannam@95: integer(C_INT), dimension(3) :: n = [N,M,L] cannam@95: plan = fftw_plan_dft_3d(3, n, in,out, FFTW_FORWARD,FFTW_ESTIMATE) cannam@95: @end example cannam@95: cannam@95: Note, by the way, that this is different from the legacy Fortran cannam@95: interface (@pxref{Fortran-interface routines}), which automatically cannam@95: reverses the order of the array dimension for you. Here, you are cannam@95: calling the C interface directly, so there is no ``translation'' layer. cannam@95: cannam@95: @cindex r2c/c2r multi-dimensional array format cannam@95: An important thing to keep in mind is the implication of this for cannam@95: multidimensional real-to-complex transforms (@pxref{Multi-Dimensional cannam@95: DFTs of Real Data}). In C, a multidimensional real-to-complex DFT cannam@95: chops the last dimension roughly in half (@threedims{N,M,L} real input cannam@95: goes to @threedims{N,M,L/2+1} complex output). In Fortran, because cannam@95: the array dimension notation is reversed, the @emph{first} dimension of cannam@95: the complex data is chopped roughly in half. For example consider the cannam@95: @samp{r2c} transform of @threedims{L,M,N} real input in Fortran: cannam@95: cannam@95: @findex fftw_plan_dft_r2c_3d cannam@95: @findex fftw_execute_dft_r2c cannam@95: @example cannam@95: type(C_PTR) :: plan cannam@95: real(C_DOUBLE), dimension(L,M,N) :: in cannam@95: complex(C_DOUBLE_COMPLEX), dimension(L/2+1,M,N) :: out cannam@95: plan = fftw_plan_dft_r2c_3d(N,M,L, in,out, FFTW_ESTIMATE) cannam@95: ... cannam@95: call fftw_execute_dft_r2c(plan, in, out) cannam@95: @end example cannam@95: cannam@95: @cindex in-place cannam@95: @cindex padding cannam@95: Alternatively, for an in-place r2c transform, as described in the C cannam@95: documentation we must @emph{pad} the @emph{first} dimension of the cannam@95: real input with an extra two entries (which are ignored by FFTW) so as cannam@95: to leave enough space for the complex output. The input is cannam@95: @emph{allocated} as a @threedims{2[L/2+1],M,N} array, even though only cannam@95: @threedims{L,M,N} of it is actually used. In this example, we will cannam@95: allocate the array as a pointer type, using @samp{fftw_alloc} to cannam@95: ensure aligned memory for maximum performance (@pxref{Allocating cannam@95: aligned memory in Fortran}); this also makes it easy to reference the cannam@95: same memory as both a real array and a complex array. cannam@95: cannam@95: @findex fftw_alloc_complex cannam@95: @findex c_f_pointer cannam@95: @example cannam@95: real(C_DOUBLE), pointer :: in(:,:,:) cannam@95: complex(C_DOUBLE_COMPLEX), pointer :: out(:,:,:) cannam@95: type(C_PTR) :: plan, data cannam@95: data = fftw_alloc_complex(int((L/2+1) * M * N, C_SIZE_T)) cannam@95: call c_f_pointer(data, in, [2*(L/2+1),M,N]) cannam@95: call c_f_pointer(data, out, [L/2+1,M,N]) cannam@95: plan = fftw_plan_dft_r2c_3d(N,M,L, in,out, FFTW_ESTIMATE) cannam@95: ... cannam@95: call fftw_execute_dft_r2c(plan, in, out) cannam@95: ... cannam@95: call fftw_destroy_plan(plan) cannam@95: call fftw_free(data) cannam@95: @end example cannam@95: cannam@95: @c ------------------------------------------------------- cannam@95: @node FFTW Fortran type reference, Plan execution in Fortran, Reversing array dimensions, Calling FFTW from Modern Fortran cannam@95: @section FFTW Fortran type reference cannam@95: cannam@95: The following are the most important type correspondences between the cannam@95: C interface and Fortran: cannam@95: cannam@95: @itemize @bullet cannam@95: cannam@95: @item cannam@95: @tindex fftw_plan cannam@95: Plans (@code{fftw_plan} and variants) are @code{type(C_PTR)} (i.e. an cannam@95: opaque pointer). cannam@95: cannam@95: @item cannam@95: @tindex fftw_complex cannam@95: @cindex precision cannam@95: @ctindex C_DOUBLE cannam@95: @ctindex C_FLOAT cannam@95: @ctindex C_LONG_DOUBLE cannam@95: @ctindex C_DOUBLE_COMPLEX cannam@95: @ctindex C_FLOAT_COMPLEX cannam@95: @ctindex C_LONG_DOUBLE_COMPLEX cannam@95: The C floating-point types @code{double}, @code{float}, and @code{long cannam@95: double} correspond to @code{real(C_DOUBLE)}, @code{real(C_FLOAT)}, and cannam@95: @code{real(C_LONG_DOUBLE)}, respectively. The C complex types cannam@95: @code{fftw_complex}, @code{fftwf_complex}, and @code{fftwl_complex} cannam@95: correspond in Fortran to @code{complex(C_DOUBLE_COMPLEX)}, cannam@95: @code{complex(C_FLOAT_COMPLEX)}, and cannam@95: @code{complex(C_LONG_DOUBLE_COMPLEX)}, respectively. cannam@95: Just as in C cannam@95: (@pxref{Precision}), the FFTW subroutines and types are prefixed with cannam@95: @samp{fftw_}, @code{fftwf_}, and @code{fftwl_} for the different precisions, and link to different libraries (@code{-lfftw3}, @code{-lfftw3f}, and @code{-lfftw3l} on Unix), but use the @emph{same} include file @code{fftw3.f03} and the @emph{same} constants (all of which begin with @samp{FFTW_}). The exception is @code{long double} precision, for which you should @emph{also} include @code{fftw3l.f03} (@pxref{Extended and quadruple precision in Fortran}). cannam@95: cannam@95: @item cannam@95: @tindex ptrdiff_t cannam@95: @ctindex C_INT cannam@95: @ctindex C_INTPTR_T cannam@95: @ctindex C_SIZE_T cannam@95: @findex fftw_malloc cannam@95: The C integer types @code{int} and @code{unsigned} (used for planner cannam@95: flags) become @code{integer(C_INT)}. The C integer type @code{ptrdiff_t} (e.g. in the @ref{64-bit Guru Interface}) becomes @code{integer(C_INTPTR_T)}, and @code{size_t} (in @code{fftw_malloc} etc.) becomes @code{integer(C_SIZE_T)}. cannam@95: cannam@95: @item cannam@95: @tindex fftw_r2r_kind cannam@95: @ctindex C_FFTW_R2R_KIND cannam@95: The @code{fftw_r2r_kind} type (@pxref{Real-to-Real Transform Kinds}) cannam@95: becomes @code{integer(C_FFTW_R2R_KIND)}. The various constant values cannam@95: of the C enumerated type (@code{FFTW_R2HC} etc.) become simply integer cannam@95: constants of the same names in Fortran. cannam@95: cannam@95: @item cannam@95: @ctindex FFTW_DESTROY_INPUT cannam@95: @cindex in-place cannam@95: @findex fftw_flops cannam@95: Numeric array pointer arguments (e.g. @code{double *}) cannam@95: become @code{dimension(*), intent(out)} arrays of the same type, or cannam@95: @code{dimension(*), intent(in)} if they are pointers to constant data cannam@95: (e.g. @code{const int *}). There are a few exceptions where numeric cannam@95: pointers refer to scalar outputs (e.g. for @code{fftw_flops}), in which cannam@95: case they are @code{intent(out)} scalar arguments in Fortran too. cannam@95: For the new-array execute functions (@pxref{New-array Execute Functions}), cannam@95: the input arrays are declared @code{dimension(*), intent(inout)}, since cannam@95: they can be modified in the case of in-place or @code{FFTW_DESTROY_INPUT} cannam@95: transforms. cannam@95: cannam@95: @item cannam@95: @findex fftw_alloc_real cannam@95: @findex c_f_pointer cannam@95: Pointer @emph{return} values (e.g @code{double *}) become cannam@95: @code{type(C_PTR)}. (If they are pointers to arrays, as for cannam@95: @code{fftw_alloc_real}, you can convert them back to Fortran array cannam@95: pointers with the standard intrinsic function @code{c_f_pointer}.) cannam@95: cannam@95: @item cannam@95: @cindex guru interface cannam@95: @tindex fftw_iodim cannam@95: @tindex fftw_iodim64 cannam@95: @cindex 64-bit architecture cannam@95: The @code{fftw_iodim} type in the guru interface (@pxref{Guru vector cannam@95: and transform sizes}) becomes @code{type(fftw_iodim)} in Fortran, a cannam@95: derived data type (the Fortran analogue of C's @code{struct}) with cannam@95: three @code{integer(C_INT)} components: @code{n}, @code{is}, and cannam@95: @code{os}, with the same meanings as in C. The @code{fftw_iodim64} type in the 64-bit guru interface (@pxref{64-bit Guru Interface}) is the same, except that its components are of type @code{integer(C_INTPTR_T)}. cannam@95: cannam@95: @item cannam@95: @ctindex C_FUNPTR cannam@95: Using the wisdom import/export functions from Fortran is a bit tricky, cannam@95: and is discussed in @ref{Accessing the wisdom API from Fortran}. In cannam@95: brief, the @code{FILE *} arguments map to @code{type(C_PTR)}, @code{const char *} to @code{character(C_CHAR), dimension(*), intent(in)} (null-terminated!), and the generic read-char/write-char functions map to @code{type(C_FUNPTR)}. cannam@95: cannam@95: @end itemize cannam@95: cannam@95: @cindex portability cannam@95: You may be wondering if you need to search-and-replace cannam@95: @code{real(kind(0.0d0))} (or whatever your favorite Fortran spelling cannam@95: of ``double precision'' is) with @code{real(C_DOUBLE)} everywhere in cannam@95: your program, and similarly for @code{complex} and @code{integer} cannam@95: types. The answer is no; you can still use your existing types. As cannam@95: long as these types match their C counterparts, things should work cannam@95: without a hitch. The worst that can happen, e.g. in the (unlikely) cannam@95: event of a system where @code{real(kind(0.0d0))} is different from cannam@95: @code{real(C_DOUBLE)}, is that the compiler will give you a cannam@95: type-mismatch error. That is, if you don't use the cannam@95: @code{iso_c_binding} kinds you need to accept at least the theoretical cannam@95: possibility of having to change your code in response to compiler cannam@95: errors on some future machine, but you don't need to worry about cannam@95: silently compiling incorrect code that yields runtime errors. cannam@95: cannam@95: @c ------------------------------------------------------- cannam@95: @node Plan execution in Fortran, Allocating aligned memory in Fortran, FFTW Fortran type reference, Calling FFTW from Modern Fortran cannam@95: @section Plan execution in Fortran cannam@95: cannam@95: In C, in order to use a plan, one normally calls @code{fftw_execute}, cannam@95: which executes the plan to perform the transform on the input/output cannam@95: arrays passed when the plan was created (@pxref{Using Plans}). The cannam@95: corresponding subroutine call in modern Fortran is: cannam@95: @example cannam@95: call fftw_execute(plan) cannam@95: @end example cannam@95: @findex fftw_execute cannam@95: cannam@95: However, we have had reports that this causes problems with some cannam@95: recent optimizing Fortran compilers. The problem is, because the cannam@95: input/output arrays are not passed as explicit arguments to cannam@95: @code{fftw_execute}, the semantics of Fortran (unlike C) allow the cannam@95: compiler to assume that the input/output arrays are not changed by cannam@95: @code{fftw_execute}. As a consequence, certain compilers end up cannam@95: repositioning the call to @code{fftw_execute}, assuming incorrectly cannam@95: that it does nothing to the arrays. cannam@95: cannam@95: There are various workarounds to this, but the safest and simplest cannam@95: thing is to not use @code{fftw_execute} in Fortran. Instead, use the cannam@95: functions described in @ref{New-array Execute Functions}, which take cannam@95: the input/output arrays as explicit arguments. For example, if the cannam@95: plan is for a complex-data DFT and was created for the arrays cannam@95: @code{in} and @code{out}, you would do: cannam@95: @example cannam@95: call fftw_execute_dft(plan, in, out) cannam@95: @end example cannam@95: @findex fftw_execute_dft cannam@95: cannam@95: There are a few things to be careful of, however: cannam@95: cannam@95: @itemize @bullet cannam@95: cannam@95: @item cannam@95: @findex fftw_execute_dft_r2c cannam@95: @findex fftw_execute_dft_c2r cannam@95: @findex fftw_execute_r2r cannam@95: You must use the correct type of execute function, matching the way cannam@95: the plan was created. Complex DFT plans should use cannam@95: @code{fftw_execute_dft}, Real-input (r2c) DFT plans should use use cannam@95: @code{fftw_execute_dft_r2c}, and real-output (c2r) DFT plans should cannam@95: use @code{fftw_execute_dft_c2r}. The various r2r plans should use cannam@95: @code{fftw_execute_r2r}. Fortunately, if you use the wrong one you cannam@95: will get a compile-time type-mismatch error (unlike legacy Fortran). cannam@95: cannam@95: @item cannam@95: You should normally pass the same input/output arrays that were used when cannam@95: creating the plan. This is always safe. cannam@95: cannam@95: @item cannam@95: @emph{If} you pass @emph{different} input/output arrays compared to cannam@95: those used when creating the plan, you must abide by all the cannam@95: restrictions of the new-array execute functions (@pxref{New-array cannam@95: Execute Functions}). The most tricky of these is the cannam@95: requirement that the new arrays have the same alignment as the cannam@95: original arrays; the best (and possibly only) way to guarantee this cannam@95: is to use the @samp{fftw_alloc} functions to allocate your arrays (@pxref{Allocating aligned memory in Fortran}). Alternatively, you can cannam@95: use the @code{FFTW_UNALIGNED} flag when creating the cannam@95: plan, in which case the plan does not depend on the alignment, but cannam@95: this may sacrifice substantial performance on architectures (like x86) cannam@95: with SIMD instructions (@pxref{SIMD alignment and fftw_malloc}). cannam@95: @ctindex FFTW_UNALIGNED cannam@95: cannam@95: @end itemize cannam@95: cannam@95: @c ------------------------------------------------------- cannam@95: @node Allocating aligned memory in Fortran, Accessing the wisdom API from Fortran, Plan execution in Fortran, Calling FFTW from Modern Fortran cannam@95: @section Allocating aligned memory in Fortran cannam@95: cannam@95: @cindex alignment cannam@95: @findex fftw_alloc_real cannam@95: @findex fftw_alloc_complex cannam@95: In order to obtain maximum performance in FFTW, you should store your cannam@95: data in arrays that have been specially aligned in memory (@pxref{SIMD cannam@95: alignment and fftw_malloc}). Enforcing alignment also permits you to cannam@95: safely use the new-array execute functions (@pxref{New-array Execute cannam@95: Functions}) to apply a given plan to more than one pair of in/out cannam@95: arrays. Unfortunately, standard Fortran arrays do @emph{not} provide cannam@95: any alignment guarantees. The @emph{only} way to allocate aligned cannam@95: memory in standard Fortran is to allocate it with an external C cannam@95: function, like the @code{fftw_alloc_real} and cannam@95: @code{fftw_alloc_complex} functions. Fortunately, Fortran 2003 provides cannam@95: a simple way to associate such allocated memory with a standard Fortran cannam@95: array pointer that you can then use normally. cannam@95: cannam@95: We therefore recommend allocating all your input/output arrays using cannam@95: the following technique: cannam@95: cannam@95: @enumerate cannam@95: cannam@95: @item cannam@95: Declare a @code{pointer}, @code{arr}, to your array of the desired type cannam@95: and dimensions. For example, @code{real(C_DOUBLE), pointer :: a(:,:)} cannam@95: for a 2d real array, or @code{complex(C_DOUBLE_COMPLEX), pointer :: cannam@95: a(:,:,:)} for a 3d complex array. cannam@95: cannam@95: @item cannam@95: The number of elements to allocate must be an cannam@95: @code{integer(C_SIZE_T)}. You can either declare a variable of this cannam@95: type, e.g. @code{integer(C_SIZE_T) :: sz}, to store the number of cannam@95: elements to allocate, or you can use the @code{int(..., C_SIZE_T)} cannam@95: intrinsic function. e.g. set @code{sz = L * M * N} or use cannam@95: @code{int(L * M * N, C_SIZE_T)} for an @threedims{L,M,N} array. cannam@95: cannam@95: @item cannam@95: Declare a @code{type(C_PTR) :: p} to hold the return value from cannam@95: FFTW's allocation routine. Set @code{p = fftw_alloc_real(sz)} for a real array, or @code{p = fftw_alloc_complex(sz)} for a complex array. cannam@95: cannam@95: @item cannam@95: @findex c_f_pointer cannam@95: Associate your pointer @code{arr} with the allocated memory @code{p} cannam@95: using the standard @code{c_f_pointer} subroutine: @code{call cannam@95: c_f_pointer(p, arr, [...dimensions...])}, where cannam@95: @code{[...dimensions...])} are an array of the dimensions of the array cannam@95: (in the usual Fortran order). e.g. @code{call c_f_pointer(p, arr, cannam@95: [L,M,N])} for an @threedims{L,M,N} array. (Alternatively, you can cannam@95: omit the dimensions argument if you specified the shape explicitly cannam@95: when declaring @code{arr}.) You can now use @code{arr} as a usual cannam@95: multidimensional array. cannam@95: cannam@95: @item cannam@95: When you are done using the array, deallocate the memory by @code{call cannam@95: fftw_free(p)} on @code{p}. cannam@95: cannam@95: @end enumerate cannam@95: cannam@95: For example, here is how we would allocate an @twodims{L,M} 2d real array: cannam@95: cannam@95: @example cannam@95: real(C_DOUBLE), pointer :: arr(:,:) cannam@95: type(C_PTR) :: p cannam@95: p = fftw_alloc_real(int(L * M, C_SIZE_T)) cannam@95: call c_f_pointer(p, arr, [L,M]) cannam@95: @emph{...use arr and arr(i,j) as usual...} cannam@95: call fftw_free(p) cannam@95: @end example cannam@95: cannam@95: and here is an @threedims{L,M,N} 3d complex array: cannam@95: cannam@95: @example cannam@95: complex(C_DOUBLE_COMPLEX), pointer :: arr(:,:,:) cannam@95: type(C_PTR) :: p cannam@95: p = fftw_alloc_complex(int(L * M * N, C_SIZE_T)) cannam@95: call c_f_pointer(p, arr, [L,M,N]) cannam@95: @emph{...use arr and arr(i,j,k) as usual...} cannam@95: call fftw_free(p) cannam@95: @end example cannam@95: cannam@95: See @ref{Reversing array dimensions} for an example allocating a cannam@95: single array and associating both real and complex array pointers with cannam@95: it, for in-place real-to-complex transforms. cannam@95: cannam@95: @c ------------------------------------------------------- cannam@95: @node Accessing the wisdom API from Fortran, Defining an FFTW module, Allocating aligned memory in Fortran, Calling FFTW from Modern Fortran cannam@95: @section Accessing the wisdom API from Fortran cannam@95: @cindex wisdom cannam@95: @cindex saving plans to disk cannam@95: cannam@95: As explained in @ref{Words of Wisdom-Saving Plans}, FFTW provides a cannam@95: ``wisdom'' API for saving plans to disk so that they can be recreated cannam@95: quickly. The C API for exporting (@pxref{Wisdom Export}) and cannam@95: importing (@pxref{Wisdom Import}) wisdom is somewhat tricky to use cannam@95: from Fortran, however, because of differences in file I/O and string cannam@95: types between C and Fortran. cannam@95: cannam@95: @menu cannam@95: * Wisdom File Export/Import from Fortran:: cannam@95: * Wisdom String Export/Import from Fortran:: cannam@95: * Wisdom Generic Export/Import from Fortran:: cannam@95: @end menu cannam@95: cannam@95: @c =========> cannam@95: @node Wisdom File Export/Import from Fortran, Wisdom String Export/Import from Fortran, Accessing the wisdom API from Fortran, Accessing the wisdom API from Fortran cannam@95: @subsection Wisdom File Export/Import from Fortran cannam@95: cannam@95: @findex fftw_import wisdom_from_filename cannam@95: @findex fftw_export_wisdom_to_filename cannam@95: The easiest way to export and import wisdom is to do so using cannam@95: @code{fftw_export_wisdom_to_filename} and cannam@95: @code{fftw_wisdom_from_filename}. The only trick is that these cannam@95: require you to pass a C string, which is an array of type cannam@95: @code{CHARACTER(C_CHAR)} that is terminated by @code{C_NULL_CHAR}. cannam@95: You can call them like this: cannam@95: cannam@95: @example cannam@95: integer(C_INT) :: ret cannam@95: ret = fftw_export_wisdom_to_filename(C_CHAR_'my_wisdom.dat' // C_NULL_CHAR) cannam@95: if (ret .eq. 0) stop 'error exporting wisdom to file' cannam@95: ret = fftw_import_wisdom_from_filename(C_CHAR_'my_wisdom.dat' // C_NULL_CHAR) cannam@95: if (ret .eq. 0) stop 'error importing wisdom from file' cannam@95: @end example cannam@95: cannam@95: Note that prepending @samp{C_CHAR_} is needed to specify that the cannam@95: literal string is of kind @code{C_CHAR}, and we null-terminate the cannam@95: string by appending @samp{// C_NULL_CHAR}. These functions return an cannam@95: @code{integer(C_INT)} (@code{ret}) which is @code{0} if an error cannam@95: occurred during export/import and nonzero otherwise. cannam@95: cannam@95: It is also possible to use the lower-level routines cannam@95: @code{fftw_export_wisdom_to_file} and cannam@95: @code{fftw_import_wisdom_from_file}, which accept parameters of the C cannam@95: type @code{FILE*}, expressed in Fortran as @code{type(C_PTR)}. cannam@95: However, you are then responsible for creating the @code{FILE*} cannam@95: yourself. You can do this by using @code{iso_c_binding} to define cannam@95: Fortran intefaces for the C library functions @code{fopen} and cannam@95: @code{fclose}, which is a bit strange in Fortran but workable. cannam@95: cannam@95: @c =========> cannam@95: @node Wisdom String Export/Import from Fortran, Wisdom Generic Export/Import from Fortran, Wisdom File Export/Import from Fortran, Accessing the wisdom API from Fortran cannam@95: @subsection Wisdom String Export/Import from Fortran cannam@95: cannam@95: @findex fftw_export_wisdom_to_string cannam@95: Dealing with FFTW's C string export/import is a bit more painful. In cannam@95: particular, the @code{fftw_export_wisdom_to_string} function requires cannam@95: you to deal with a dynamically allocated C string. To get its length, cannam@95: you must define an interface to the C @code{strlen} function, and to cannam@95: deallocate it you must define an interface to C @code{free}: cannam@95: cannam@95: @example cannam@95: use, intrinsic :: iso_c_binding cannam@95: interface cannam@95: integer(C_INT) function strlen(s) bind(C, name='strlen') cannam@95: import cannam@95: type(C_PTR), value :: s cannam@95: end function strlen cannam@95: subroutine free(p) bind(C, name='free') cannam@95: import cannam@95: type(C_PTR), value :: p cannam@95: end subroutine free cannam@95: end interface cannam@95: @end example cannam@95: cannam@95: Given these definitions, you can then export wisdom to a Fortran cannam@95: character array: cannam@95: cannam@95: @example cannam@95: character(C_CHAR), pointer :: s(:) cannam@95: integer(C_SIZE_T) :: slen cannam@95: type(C_PTR) :: p cannam@95: p = fftw_export_wisdom_to_string() cannam@95: if (.not. c_associated(p)) stop 'error exporting wisdom' cannam@95: slen = strlen(p) cannam@95: call c_f_pointer(p, s, [slen+1]) cannam@95: ... cannam@95: call free(p) cannam@95: @end example cannam@95: @findex c_associated cannam@95: @findex c_f_pointer cannam@95: cannam@95: Note that @code{slen} is the length of the C string, but the length of cannam@95: the array is @code{slen+1} because it includes the terminating null cannam@95: character. (You can omit the @samp{+1} if you don't want Fortran to cannam@95: know about the null character.) The standard @code{c_associated} function cannam@95: checks whether @code{p} is a null pointer, which is returned by cannam@95: @code{fftw_export_wisdom_to_string} if there was an error. cannam@95: cannam@95: @findex fftw_import_wisdom_from_string cannam@95: To import wisdom from a string, use cannam@95: @code{fftw_import_wisdom_from_string} as usual; note that the argument cannam@95: of this function must be a @code{character(C_CHAR)} that is terminated cannam@95: by the @code{C_NULL_CHAR} character, like the @code{s} array above. cannam@95: cannam@95: @c =========> cannam@95: @node Wisdom Generic Export/Import from Fortran, , Wisdom String Export/Import from Fortran, Accessing the wisdom API from Fortran cannam@95: @subsection Wisdom Generic Export/Import from Fortran cannam@95: cannam@95: The most generic wisdom export/import functions allow you to provide cannam@95: an arbitrary callback function to read/write one character at a time cannam@95: in any way you want. However, your callback function must be written cannam@95: in a special way, using the @code{bind(C)} attribute to be passed to a cannam@95: C interface. cannam@95: cannam@95: @findex fftw_export_wisdom cannam@95: In particular, to call the generic wisdom export function cannam@95: @code{fftw_export_wisdom}, you would write a callback subroutine of the form: cannam@95: cannam@95: @example cannam@95: subroutine my_write_char(c, p) bind(C) cannam@95: use, intrinsic :: iso_c_binding cannam@95: character(C_CHAR), value :: c cannam@95: type(C_PTR), value :: p cannam@95: @emph{...write c...} cannam@95: end subroutine my_write_char cannam@95: @end example cannam@95: cannam@95: Given such a subroutine (along with the corresponding interface definition), you could then export wisdom using: cannam@95: cannam@95: @findex c_funloc cannam@95: @example cannam@95: call fftw_export_wisdom(c_funloc(my_write_char), p) cannam@95: @end example cannam@95: cannam@95: @findex c_loc cannam@95: @findex c_f_pointer cannam@95: The standard @code{c_funloc} intrinsic converts a Fortran cannam@95: @code{bind(C)} subroutine into a C function pointer. The parameter cannam@95: @code{p} is a @code{type(C_PTR)} to any arbitrary data that you want cannam@95: to pass to @code{my_write_char} (or @code{C_NULL_PTR} if none). (Note cannam@95: that you can get a C pointer to Fortran data using the intrinsic cannam@95: @code{c_loc}, and convert it back to a Fortran pointer in cannam@95: @code{my_write_char} using @code{c_f_pointer}.) cannam@95: cannam@95: Similarly, to use the generic @code{fftw_import_wisdom}, you would cannam@95: define a callback function of the form: cannam@95: cannam@95: @findex fftw_import_wisdom cannam@95: @example cannam@95: integer(C_INT) function my_read_char(p) bind(C) cannam@95: use, intrinsic :: iso_c_binding cannam@95: type(C_PTR), value :: p cannam@95: character :: c cannam@95: @emph{...read a character c...} cannam@95: my_read_char = ichar(c, C_INT) cannam@95: end function my_read_char cannam@95: cannam@95: .... cannam@95: cannam@95: integer(C_INT) :: ret cannam@95: ret = fftw_import_wisdom(c_funloc(my_read_char), p) cannam@95: if (ret .eq. 0) stop 'error importing wisdom' cannam@95: @end example cannam@95: cannam@95: Your function can return @code{-1} if the end of the input is reached. cannam@95: Again, @code{p} is an arbitrary @code{type(C_PTR} that is passed cannam@95: through to your function. @code{fftw_import_wisdom} returns @code{0} cannam@95: if an error occurred and nonzero otherwise. cannam@95: cannam@95: @c ------------------------------------------------------- cannam@95: @node Defining an FFTW module, , Accessing the wisdom API from Fortran, Calling FFTW from Modern Fortran cannam@95: @section Defining an FFTW module cannam@95: cannam@95: Rather than using the @code{include} statement to include the cannam@95: @code{fftw3.f03} interface file in any subroutine where you want to cannam@95: use FFTW, you might prefer to define an FFTW Fortran module. FFTW cannam@95: does not install itself as a module, primarily because cannam@95: @code{fftw3.f03} can be shared between different Fortran compilers while cannam@95: modules (in general) cannot. However, it is trivial to define your cannam@95: own FFTW module if you want. Just create a file containing: cannam@95: cannam@95: @example cannam@95: module FFTW3 cannam@95: use, intrinsic :: iso_c_binding cannam@95: include 'fftw3.f03' cannam@95: end module cannam@95: @end example cannam@95: cannam@95: Compile this file into a module as usual for your compiler (e.g. with cannam@95: @code{gfortran -c} you will get a file @code{fftw3.mod}). Now, cannam@95: instead of @code{include 'fftw3.f03'}, whenever you want to use FFTW cannam@95: routines you can just do: cannam@95: cannam@95: @example cannam@95: use FFTW3 cannam@95: @end example cannam@95: cannam@95: as usual for Fortran modules. (You still need to link to the FFTW cannam@95: library, of course.)