annotate fft/fftw/fftw-3.3.4/doc/modern-fortran.texi @ 40:223f770b5341 kissfft-double tip

Try a double-precision kissfft

@node Calling FFTW from Modern Fortran, Calling FFTW from Legacy Fortran, Distributed-memory FFTW with MPI, Top
@chapter Calling FFTW from Modern Fortran
@cindex Fortran interface

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
(@pxref{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).
@cindex portability

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.

@menu
* Overview of Fortran interface::
* Reversing array dimensions::
* FFTW Fortran type reference::
* Plan execution in Fortran::
* Allocating aligned memory in Fortran::
* Accessing the wisdom API from Fortran::
* Defining an FFTW module::
@end menu

@c -------------------------------------------------------
@node Overview of Fortran interface, Reversing array dimensions, Calling FFTW from Modern Fortran, Calling FFTW from Modern Fortran
@section Overview of Fortran interface

FFTW provides a file @code{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
@code{fftw3.h} (the C header file).  In any Fortran subroutine where
you want to use FFTW functions, you should begin with:

@cindex iso_c_binding
@example
  use, intrinsic :: iso_c_binding 
  include 'fftw3.f03'
@end example

This includes the interface definitions and the standard
@code{iso_c_binding} module (which defines the equivalents of C
types).  You can also put the FFTW functions into a module if you
prefer (@pxref{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:

@findex fftw_plan_dft_2d
@findex fftw_execute_dft
@findex fftw_destroy_plan
@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)
@end example

A few important things to keep in mind are:

@itemize @bullet

@item
@tindex fftw_complex
@ctindex C_PTR
@ctindex C_INT
@ctindex C_DOUBLE
@ctindex C_DOUBLE_COMPLEX
FFTW plans are @code{type(C_PTR)}.  Other C types are mapped in the
obvious way via the @code{iso_c_binding} standard: @code{int} turns
into @code{integer(C_INT)}, @code{fftw_complex} turns into
@code{complex(C_DOUBLE_COMPLEX)}, @code{double} turns into
@code{real(C_DOUBLE)}, and so on. @xref{FFTW Fortran type reference}.

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

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

@item
@cindex alignment
@cindex SIMD
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 (@pxref{SIMD
alignment and fftw_malloc}). Better performance will often be obtained
by allocating with @samp{fftw_alloc}. @xref{Allocating aligned memory
in Fortran}.

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

@item
@cindex flags
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.)

@end itemize

@menu
* Extended and quadruple precision in Fortran::
@end menu

@node Extended and quadruple precision in Fortran,  , Overview of Fortran interface, Overview of Fortran interface
@subsection Extended and quadruple precision in Fortran
@cindex precision

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

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

@cindex iso_c_binding
@example
  use, intrinsic :: iso_c_binding 
  include 'fftw3.f03'
  include 'fftw3l.f03'
@end example

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

@c -------------------------------------------------------
@node Reversing array dimensions, FFTW Fortran type reference, Overview of Fortran interface, Calling FFTW from Modern Fortran
@section Reversing array dimensions

@cindex row-major
@cindex column-major
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). @xref{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 @emph{reverse the
ordering of the dimensions}.

For example, consider the three-dimensional (@threedims{L,M,N}) arrays:

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

To plan a DFT for these arrays using @code{fftw_plan_dft_3d}, you could do:

@findex fftw_plan_dft_3d
@example
  plan = fftw_plan_dft_3d(N,M,L, in,out, FFTW_FORWARD,FFTW_ESTIMATE)
@end example

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

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

Note, by the way, that this is different from the legacy Fortran
interface (@pxref{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.

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

@findex fftw_plan_dft_r2c_3d
@findex fftw_execute_dft_r2c
@example
  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)
@end example

@cindex in-place
@cindex padding
Alternatively, for an in-place r2c transform, as described in the C
documentation we must @emph{pad} the @emph{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
@emph{allocated} as a @threedims{2[L/2+1],M,N} array, even though only
@threedims{L,M,N} of it is actually used.  In this example, we will
allocate the array as a pointer type, using @samp{fftw_alloc} to
ensure aligned memory for maximum performance (@pxref{Allocating
aligned memory in Fortran}); this also makes it easy to reference the
same memory as both a real array and a complex array.

@findex fftw_alloc_complex
@findex c_f_pointer
@example
  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)
@end example

@c -------------------------------------------------------
@node FFTW Fortran type reference, Plan execution in Fortran, Reversing array dimensions, Calling FFTW from Modern Fortran
@section FFTW Fortran type reference

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

@itemize @bullet

@item
@tindex fftw_plan
Plans (@code{fftw_plan} and variants) are @code{type(C_PTR)} (i.e. an
opaque pointer).

@item
@tindex fftw_complex
@cindex precision
@ctindex C_DOUBLE
@ctindex C_FLOAT
@ctindex C_LONG_DOUBLE
@ctindex C_DOUBLE_COMPLEX
@ctindex C_FLOAT_COMPLEX
@ctindex C_LONG_DOUBLE_COMPLEX
The C floating-point types @code{double}, @code{float}, and @code{long
double} correspond to @code{real(C_DOUBLE)}, @code{real(C_FLOAT)}, and
@code{real(C_LONG_DOUBLE)}, respectively.  The C complex types
@code{fftw_complex}, @code{fftwf_complex}, and @code{fftwl_complex}
correspond in Fortran to @code{complex(C_DOUBLE_COMPLEX)},
@code{complex(C_FLOAT_COMPLEX)}, and
@code{complex(C_LONG_DOUBLE_COMPLEX)}, respectively.  
Just as in C
(@pxref{Precision}), the FFTW subroutines and types are prefixed with
@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}).

@item
@tindex ptrdiff_t
@ctindex C_INT
@ctindex C_INTPTR_T
@ctindex C_SIZE_T
@findex fftw_malloc
The C integer types @code{int} and @code{unsigned} (used for planner
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)}.

@item
@tindex fftw_r2r_kind
@ctindex C_FFTW_R2R_KIND
The @code{fftw_r2r_kind} type (@pxref{Real-to-Real Transform Kinds})
becomes @code{integer(C_FFTW_R2R_KIND)}.  The various constant values
of the C enumerated type (@code{FFTW_R2HC} etc.) become simply integer
constants of the same names in Fortran.

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

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

@item
@cindex guru interface
@tindex fftw_iodim
@tindex fftw_iodim64
@cindex 64-bit architecture
The @code{fftw_iodim} type in the guru interface (@pxref{Guru vector
and transform sizes}) becomes @code{type(fftw_iodim)} in Fortran, a
derived data type (the Fortran analogue of C's @code{struct}) with
three @code{integer(C_INT)} components: @code{n}, @code{is}, and
@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)}.

@item
@ctindex C_FUNPTR
Using the wisdom import/export functions from Fortran is a bit tricky,
and is discussed in @ref{Accessing the wisdom API from Fortran}.  In
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)}.

@end itemize

@cindex portability
You may be wondering if you need to search-and-replace
@code{real(kind(0.0d0))} (or whatever your favorite Fortran spelling
of ``double precision'' is) with @code{real(C_DOUBLE)} everywhere in
your program, and similarly for @code{complex} and @code{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 @code{real(kind(0.0d0))} is different from
@code{real(C_DOUBLE)}, is that the compiler will give you a
type-mismatch error.  That is, if you don't use the
@code{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.

@c -------------------------------------------------------
@node Plan execution in Fortran, Allocating aligned memory in Fortran, FFTW Fortran type reference, Calling FFTW from Modern Fortran
@section Plan execution in Fortran

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

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
@code{fftw_execute}, the semantics of Fortran (unlike C) allow the
compiler to assume that the input/output arrays are not changed by
@code{fftw_execute}.  As a consequence, certain compilers end up
repositioning the call to @code{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 @code{fftw_execute} in Fortran.  Instead, use the
functions described in @ref{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
@code{in} and @code{out}, you would do:
@example
 call fftw_execute_dft(plan, in, out)
@end example
@findex fftw_execute_dft

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

@itemize @bullet

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

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

@item
@emph{If} you pass @emph{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 (@pxref{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 @samp{fftw_alloc} functions to allocate your arrays (@pxref{Allocating aligned memory in Fortran}). Alternatively, you can
use the @code{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 (@pxref{SIMD alignment and fftw_malloc}).
@ctindex FFTW_UNALIGNED

@end itemize

@c -------------------------------------------------------
@node Allocating aligned memory in Fortran, Accessing the wisdom API from Fortran, Plan execution in Fortran, Calling FFTW from Modern Fortran
@section Allocating aligned memory in Fortran

@cindex alignment
@findex fftw_alloc_real
@findex fftw_alloc_complex
In order to obtain maximum performance in FFTW, you should store your
data in arrays that have been specially aligned in memory (@pxref{SIMD
alignment and fftw_malloc}).  Enforcing alignment also permits you to
safely use the new-array execute functions (@pxref{New-array Execute
Functions}) to apply a given plan to more than one pair of in/out
arrays.  Unfortunately, standard Fortran arrays do @emph{not} provide
any alignment guarantees.  The @emph{only} way to allocate aligned
memory in standard Fortran is to allocate it with an external C
function, like the @code{fftw_alloc_real} and
@code{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:

@enumerate

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

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

@item
Declare a @code{type(C_PTR) :: p} to hold the return value from
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.

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

@item
When you are done using the array, deallocate the memory by @code{call
fftw_free(p)} on @code{p}.

@end enumerate

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

@example
  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])
  @emph{...use arr and arr(i,j) as usual...}
  call fftw_free(p)
@end example

and here is an @threedims{L,M,N} 3d complex array:

@example
  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])
  @emph{...use arr and arr(i,j,k) as usual...}
  call fftw_free(p)
@end example

See @ref{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.

@c -------------------------------------------------------
@node Accessing the wisdom API from Fortran, Defining an FFTW module, Allocating aligned memory in Fortran, Calling FFTW from Modern Fortran
@section Accessing the wisdom API from Fortran
@cindex wisdom
@cindex saving plans to disk

As explained in @ref{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 (@pxref{Wisdom Export}) and
importing (@pxref{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.

@menu
* Wisdom File Export/Import from Fortran::
* Wisdom String Export/Import from Fortran::
* Wisdom Generic Export/Import from Fortran::
@end menu

@c =========>
@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
@subsection Wisdom File Export/Import from Fortran

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

@example
  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'
@end example

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

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

@c =========>
@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
@subsection Wisdom String Export/Import from Fortran

@findex fftw_export_wisdom_to_string
Dealing with FFTW's C string export/import is a bit more painful.  In
particular, the @code{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 @code{strlen} function, and to
deallocate it you must define an interface to C @code{free}:

@example
  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
@end example

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

@example
  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)
@end example
@findex c_associated
@findex c_f_pointer

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

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

@c =========>
@node Wisdom Generic Export/Import from Fortran,  , Wisdom String Export/Import from Fortran, Accessing the wisdom API from Fortran
@subsection 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 @code{bind(C)} attribute to be passed to a
C interface.

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

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

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

@findex c_funloc
@example
  call fftw_export_wisdom(c_funloc(my_write_char), p)
@end example

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

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

@findex fftw_import_wisdom
@example
  integer(C_INT) function my_read_char(p) bind(C)
    use, intrinsic :: iso_c_binding
    type(C_PTR), value :: p
    character :: c
    @emph{...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'
@end example

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

@c -------------------------------------------------------
@node Defining an FFTW module,  , Accessing the wisdom API from Fortran, Calling FFTW from Modern Fortran
@section Defining an FFTW module

Rather than using the @code{include} statement to include the
@code{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
@code{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:

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

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

@example
  use FFTW3
@end example

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