1@node Calling FFTW from Modern Fortran, Calling FFTW from Legacy Fortran, Distributed-memory FFTW with MPI, Top
2@chapter Calling FFTW from Modern Fortran
3@cindex Fortran interface
4
5Fortran 2003 standardized ways for Fortran code to call C libraries,
6and this allows us to support a direct translation of the FFTW C API
7into Fortran.  Compared to the legacy Fortran 77 interface
8(@pxref{Calling FFTW from Legacy Fortran}), this direct interface
9offers many advantages, especially compile-time type-checking and
10aligned memory allocation.  As of this writing, support for these C
11interoperability features seems widespread, having been implemented in
12nearly all major Fortran compilers (e.g. GNU, Intel, IBM,
13Oracle/Solaris, Portland Group, NAG).
14@cindex portability
15
16This chapter documents that interface.  For the most part, since this
17interface allows Fortran to call the C interface directly, the usage
18is identical to C translated to Fortran syntax.  However, there are a
19few subtle points such as memory allocation, wisdom, and data types
20that deserve closer attention.
21
22@menu
23* Overview of Fortran interface::
24* Reversing array dimensions::
25* FFTW Fortran type reference::
26* Plan execution in Fortran::
27* Allocating aligned memory in Fortran::
28* Accessing the wisdom API from Fortran::
29* Defining an FFTW module::
30@end menu
31
32@c -------------------------------------------------------
33@node Overview of Fortran interface, Reversing array dimensions, Calling FFTW from Modern Fortran, Calling FFTW from Modern Fortran
34@section Overview of Fortran interface
35
36FFTW provides a file @code{fftw3.f03} that defines Fortran 2003
37interfaces for all of its C routines, except for the MPI routines
38described elsewhere, which can be found in the same directory as
39@code{fftw3.h} (the C header file).  In any Fortran subroutine where
40you want to use FFTW functions, you should begin with:
41
42@cindex iso_c_binding
43@example
44  use, intrinsic :: iso_c_binding
45  include 'fftw3.f03'
46@end example
47
48This includes the interface definitions and the standard
49@code{iso_c_binding} module (which defines the equivalents of C
50types).  You can also put the FFTW functions into a module if you
51prefer (@pxref{Defining an FFTW module}).
52
53At this point, you can now call anything in the FFTW C interface
54directly, almost exactly as in C other than minor changes in syntax.
55For example:
56
57@findex fftw_plan_dft_2d
58@findex fftw_execute_dft
59@findex fftw_destroy_plan
60@example
61  type(C_PTR) :: plan
62  complex(C_DOUBLE_COMPLEX), dimension(1024,1000) :: in, out
63  plan = fftw_plan_dft_2d(1000,1024, in,out, FFTW_FORWARD,FFTW_ESTIMATE)
64  ...
65  call fftw_execute_dft(plan, in, out)
66  ...
67  call fftw_destroy_plan(plan)
68@end example
69
70A few important things to keep in mind are:
71
72@itemize @bullet
73
74@item
75@tindex fftw_complex
76@ctindex C_PTR
77@ctindex C_INT
78@ctindex C_DOUBLE
79@ctindex C_DOUBLE_COMPLEX
80FFTW plans are @code{type(C_PTR)}.  Other C types are mapped in the
81obvious way via the @code{iso_c_binding} standard: @code{int} turns
82into @code{integer(C_INT)}, @code{fftw_complex} turns into
83@code{complex(C_DOUBLE_COMPLEX)}, @code{double} turns into
84@code{real(C_DOUBLE)}, and so on. @xref{FFTW Fortran type reference}.
85
86@item
87Functions in C become functions in Fortran if they have a return value,
88and subroutines in Fortran otherwise.
89
90@item
91The ordering of the Fortran array dimensions must be @emph{reversed}
92when they are passed to the FFTW plan creation, thanks to differences
93in array indexing conventions (@pxref{Multi-dimensional Array
94Format}).  This is @emph{unlike} the legacy Fortran interface
95(@pxref{Fortran-interface routines}), which reversed the dimensions
96for you.  @xref{Reversing array dimensions}.
97
98@item
99@cindex alignment
100@cindex SIMD
101Using ordinary Fortran array declarations like this works, but may
102yield suboptimal performance because the data may not be not aligned
103to exploit SIMD instructions on modern proessors (@pxref{SIMD
104alignment and fftw_malloc}). Better performance will often be obtained
105by allocating with @samp{fftw_alloc}. @xref{Allocating aligned memory
106in Fortran}.
107
108@item
109@findex fftw_execute
110Similar to the legacy Fortran interface (@pxref{FFTW Execution in
111Fortran}), we currently recommend @emph{not} using @code{fftw_execute}
112but rather using the more specialized functions like
113@code{fftw_execute_dft} (@pxref{New-array Execute Functions}).
114However, you should execute the plan on the @code{same arrays} as the
115ones for which you created the plan, unless you are especially
116careful.  @xref{Plan execution in Fortran}.  To prevent
117you from using @code{fftw_execute} by mistake, the @code{fftw3.f03}
118file does not provide an @code{fftw_execute} interface declaration.
119
120@item
121@cindex flags
122Multiple 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.)
123
124@end itemize
125
126@menu
127* Extended and quadruple precision in Fortran::
128@end menu
129
130@node Extended and quadruple precision in Fortran,  , Overview of Fortran interface, Overview of Fortran interface
131@subsection Extended and quadruple precision in Fortran
132@cindex precision
133
134If FFTW is compiled in @code{long double} (extended) precision
135(@pxref{Installation and Customization}), you may be able to call the
136resulting @code{fftwl_} routines (@pxref{Precision}) from Fortran if
137your compiler supports the @code{C_LONG_DOUBLE_COMPLEX} type code.
138
139Because some Fortran compilers do not support
140@code{C_LONG_DOUBLE_COMPLEX}, the @code{fftwl_} declarations are
141segregated into a separate interface file @code{fftw3l.f03}, which you
142should include @emph{in addition} to @code{fftw3.f03} (which declares
143precision-independent @samp{FFTW_} constants):
144
145@cindex iso_c_binding
146@example
147  use, intrinsic :: iso_c_binding
148  include 'fftw3.f03'
149  include 'fftw3l.f03'
150@end example
151
152We also support using the nonstandard @code{__float128}
153quadruple-precision type provided by recent versions of @code{gcc} on
15432- and 64-bit x86 hardware (@pxref{Installation and Customization}),
155using the corresponding @code{real(16)} and @code{complex(16)} types
156supported by @code{gfortran}.  The quadruple-precision @samp{fftwq_}
157functions (@pxref{Precision}) are declared in a @code{fftw3q.f03}
158interface file, which should be included in addition to
159@code{fftw3.f03}, as above.  You should also link with
160@code{-lfftw3q -lquadmath -lm} as in C.
161
162@c -------------------------------------------------------
163@node Reversing array dimensions, FFTW Fortran type reference, Overview of Fortran interface, Calling FFTW from Modern Fortran
164@section Reversing array dimensions
165
166@cindex row-major
167@cindex column-major
168A minor annoyance in calling FFTW from Fortran is that FFTW's array
169dimensions are defined in the C convention (row-major order), while
170Fortran's array dimensions are the opposite convention (column-major
171order). @xref{Multi-dimensional Array Format}.  This is just a
172bookkeeping difference, with no effect on performance.  The only
173consequence of this is that, whenever you create an FFTW plan for a
174multi-dimensional transform, you must always @emph{reverse the
175ordering of the dimensions}.
176
177For example, consider the three-dimensional (@threedims{L,M,N}) arrays:
178
179@example
180  complex(C_DOUBLE_COMPLEX), dimension(L,M,N) :: in, out
181@end example
182
183To plan a DFT for these arrays using @code{fftw_plan_dft_3d}, you could do:
184
185@findex fftw_plan_dft_3d
186@example
187  plan = fftw_plan_dft_3d(N,M,L, in,out, FFTW_FORWARD,FFTW_ESTIMATE)
188@end example
189
190That is, from FFTW's perspective this is a @threedims{N,M,L} array.
191@emph{No data transposition need occur}, as this is @emph{only
192notation}.  Similarly, to use the more generic routine
193@code{fftw_plan_dft} with the same arrays, you could do:
194
195@example
196  integer(C_INT), dimension(3) :: n = [N,M,L]
197  plan = fftw_plan_dft_3d(3, n, in,out, FFTW_FORWARD,FFTW_ESTIMATE)
198@end example
199
200Note, by the way, that this is different from the legacy Fortran
201interface (@pxref{Fortran-interface routines}), which automatically
202reverses the order of the array dimension for you.  Here, you are
203calling the C interface directly, so there is no ``translation'' layer.
204
205@cindex r2c/c2r multi-dimensional array format
206An important thing to keep in mind is the implication of this for
207multidimensional real-to-complex transforms (@pxref{Multi-Dimensional
208DFTs of Real Data}).  In C, a multidimensional real-to-complex DFT
209chops the last dimension roughly in half (@threedims{N,M,L} real input
210goes to @threedims{N,M,L/2+1} complex output).  In Fortran, because
211the array dimension notation is reversed, the @emph{first} dimension of
212the complex data is chopped roughly in half.  For example consider the
213@samp{r2c} transform of @threedims{L,M,N} real input in Fortran:
214
215@findex fftw_plan_dft_r2c_3d
216@findex fftw_execute_dft_r2c
217@example
218  type(C_PTR) :: plan
219  real(C_DOUBLE), dimension(L,M,N) :: in
220  complex(C_DOUBLE_COMPLEX), dimension(L/2+1,M,N) :: out
221  plan = fftw_plan_dft_r2c_3d(N,M,L, in,out, FFTW_ESTIMATE)
222  ...
223  call fftw_execute_dft_r2c(plan, in, out)
224@end example
225
226@cindex in-place
227@cindex padding
228Alternatively, for an in-place r2c transform, as described in the C
229documentation we must @emph{pad} the @emph{first} dimension of the
230real input with an extra two entries (which are ignored by FFTW) so as
231to leave enough space for the complex output. The input is
232@emph{allocated} as a @threedims{2[L/2+1],M,N} array, even though only
233@threedims{L,M,N} of it is actually used.  In this example, we will
234allocate the array as a pointer type, using @samp{fftw_alloc} to
235ensure aligned memory for maximum performance (@pxref{Allocating
236aligned memory in Fortran}); this also makes it easy to reference the
237same memory as both a real array and a complex array.
238
239@findex fftw_alloc_complex
240@findex c_f_pointer
241@example
242  real(C_DOUBLE), pointer :: in(:,:,:)
243  complex(C_DOUBLE_COMPLEX), pointer :: out(:,:,:)
244  type(C_PTR) :: plan, data
245  data = fftw_alloc_complex(int((L/2+1) * M * N, C_SIZE_T))
246  call c_f_pointer(data, in, [2*(L/2+1),M,N])
247  call c_f_pointer(data, out, [L/2+1,M,N])
248  plan = fftw_plan_dft_r2c_3d(N,M,L, in,out, FFTW_ESTIMATE)
249  ...
250  call fftw_execute_dft_r2c(plan, in, out)
251  ...
252  call fftw_destroy_plan(plan)
253  call fftw_free(data)
254@end example
255
256@c -------------------------------------------------------
257@node FFTW Fortran type reference, Plan execution in Fortran, Reversing array dimensions, Calling FFTW from Modern Fortran
258@section FFTW Fortran type reference
259
260The following are the most important type correspondences between the
261C interface and Fortran:
262
263@itemize @bullet
264
265@item
266@tindex fftw_plan
267Plans (@code{fftw_plan} and variants) are @code{type(C_PTR)} (i.e. an
268opaque pointer).
269
270@item
271@tindex fftw_complex
272@cindex precision
273@ctindex C_DOUBLE
274@ctindex C_FLOAT
275@ctindex C_LONG_DOUBLE
276@ctindex C_DOUBLE_COMPLEX
277@ctindex C_FLOAT_COMPLEX
278@ctindex C_LONG_DOUBLE_COMPLEX
279The C floating-point types @code{double}, @code{float}, and @code{long
280double} correspond to @code{real(C_DOUBLE)}, @code{real(C_FLOAT)}, and
281@code{real(C_LONG_DOUBLE)}, respectively.  The C complex types
282@code{fftw_complex}, @code{fftwf_complex}, and @code{fftwl_complex}
283correspond in Fortran to @code{complex(C_DOUBLE_COMPLEX)},
284@code{complex(C_FLOAT_COMPLEX)}, and
285@code{complex(C_LONG_DOUBLE_COMPLEX)}, respectively.
286Just as in C
287(@pxref{Precision}), the FFTW subroutines and types are prefixed with
288@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}).
289
290@item
291@tindex ptrdiff_t
292@ctindex C_INT
293@ctindex C_INTPTR_T
294@ctindex C_SIZE_T
295@findex fftw_malloc
296The C integer types @code{int} and @code{unsigned} (used for planner
297flags) 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)}.
298
299@item
300@tindex fftw_r2r_kind
301@ctindex C_FFTW_R2R_KIND
302The @code{fftw_r2r_kind} type (@pxref{Real-to-Real Transform Kinds})
303becomes @code{integer(C_FFTW_R2R_KIND)}.  The various constant values
304of the C enumerated type (@code{FFTW_R2HC} etc.) become simply integer
305constants of the same names in Fortran.
306
307@item
308@ctindex FFTW_DESTROY_INPUT
309@cindex in-place
310@findex fftw_flops
311Numeric array pointer arguments (e.g. @code{double *})
312become @code{dimension(*), intent(out)} arrays of the same type, or
313@code{dimension(*), intent(in)} if they are pointers to constant data
314(e.g. @code{const int *}).  There are a few exceptions where numeric
315pointers refer to scalar outputs (e.g. for @code{fftw_flops}), in which
316case they are @code{intent(out)} scalar arguments in Fortran too.
317For the new-array execute functions (@pxref{New-array Execute Functions}),
318the input arrays are declared @code{dimension(*), intent(inout)}, since
319they can be modified in the case of in-place or @code{FFTW_DESTROY_INPUT}
320transforms.
321
322@item
323@findex fftw_alloc_real
324@findex c_f_pointer
325Pointer @emph{return} values (e.g @code{double *}) become
326@code{type(C_PTR)}.  (If they are pointers to arrays, as for
327@code{fftw_alloc_real}, you can convert them back to Fortran array
328pointers with the standard intrinsic function @code{c_f_pointer}.)
329
330@item
331@cindex guru interface
332@tindex fftw_iodim
333@tindex fftw_iodim64
334@cindex 64-bit architecture
335The @code{fftw_iodim} type in the guru interface (@pxref{Guru vector
336and transform sizes}) becomes @code{type(fftw_iodim)} in Fortran, a
337derived data type (the Fortran analogue of C's @code{struct}) with
338three @code{integer(C_INT)} components: @code{n}, @code{is}, and
339@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)}.
340
341@item
342@ctindex C_FUNPTR
343Using the wisdom import/export functions from Fortran is a bit tricky,
344and is discussed in @ref{Accessing the wisdom API from Fortran}.  In
345brief, 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)}.
346
347@end itemize
348
349@cindex portability
350You may be wondering if you need to search-and-replace
351@code{real(kind(0.0d0))} (or whatever your favorite Fortran spelling
352of ``double precision'' is) with @code{real(C_DOUBLE)} everywhere in
353your program, and similarly for @code{complex} and @code{integer}
354types.  The answer is no; you can still use your existing types.  As
355long as these types match their C counterparts, things should work
356without a hitch.  The worst that can happen, e.g. in the (unlikely)
357event of a system where @code{real(kind(0.0d0))} is different from
358@code{real(C_DOUBLE)}, is that the compiler will give you a
359type-mismatch error.  That is, if you don't use the
360@code{iso_c_binding} kinds you need to accept at least the theoretical
361possibility of having to change your code in response to compiler
362errors on some future machine, but you don't need to worry about
363silently compiling incorrect code that yields runtime errors.
364
365@c -------------------------------------------------------
366@node Plan execution in Fortran, Allocating aligned memory in Fortran, FFTW Fortran type reference, Calling FFTW from Modern Fortran
367@section Plan execution in Fortran
368
369In C, in order to use a plan, one normally calls @code{fftw_execute},
370which executes the plan to perform the transform on the input/output
371arrays passed when the plan was created (@pxref{Using Plans}).  The
372corresponding subroutine call in modern Fortran is:
373@example
374 call fftw_execute(plan)
375@end example
376@findex fftw_execute
377
378However, we have had reports that this causes problems with some
379recent optimizing Fortran compilers.  The problem is, because the
380input/output arrays are not passed as explicit arguments to
381@code{fftw_execute}, the semantics of Fortran (unlike C) allow the
382compiler to assume that the input/output arrays are not changed by
383@code{fftw_execute}.  As a consequence, certain compilers end up
384repositioning the call to @code{fftw_execute}, assuming incorrectly
385that it does nothing to the arrays.
386
387There are various workarounds to this, but the safest and simplest
388thing is to not use @code{fftw_execute} in Fortran.  Instead, use the
389functions described in @ref{New-array Execute Functions}, which take
390the input/output arrays as explicit arguments.  For example, if the
391plan is for a complex-data DFT and was created for the arrays
392@code{in} and @code{out}, you would do:
393@example
394 call fftw_execute_dft(plan, in, out)
395@end example
396@findex fftw_execute_dft
397
398There are a few things to be careful of, however:
399
400@itemize @bullet
401
402@item
403@findex fftw_execute_dft_r2c
404@findex fftw_execute_dft_c2r
405@findex fftw_execute_r2r
406You must use the correct type of execute function, matching the way
407the plan was created.  Complex DFT plans should use
408@code{fftw_execute_dft}, Real-input (r2c) DFT plans should use use
409@code{fftw_execute_dft_r2c}, and real-output (c2r) DFT plans should
410use @code{fftw_execute_dft_c2r}.  The various r2r plans should use
411@code{fftw_execute_r2r}.  Fortunately, if you use the wrong one you
412will get a compile-time type-mismatch error (unlike legacy Fortran).
413
414@item
415You should normally pass the same input/output arrays that were used when
416creating the plan.  This is always safe.
417
418@item
419@emph{If} you pass @emph{different} input/output arrays compared to
420those used when creating the plan, you must abide by all the
421restrictions of the new-array execute functions (@pxref{New-array
422Execute Functions}).  The most tricky of these is the
423requirement that the new arrays have the same alignment as the
424original arrays; the best (and possibly only) way to guarantee this
425is to use the @samp{fftw_alloc} functions to allocate your arrays (@pxref{Allocating aligned memory in Fortran}). Alternatively, you can
426use the @code{FFTW_UNALIGNED} flag when creating the
427plan, in which case the plan does not depend on the alignment, but
428this may sacrifice substantial performance on architectures (like x86)
429with SIMD instructions (@pxref{SIMD alignment and fftw_malloc}).
430@ctindex FFTW_UNALIGNED
431
432@end itemize
433
434@c -------------------------------------------------------
435@node Allocating aligned memory in Fortran, Accessing the wisdom API from Fortran, Plan execution in Fortran, Calling FFTW from Modern Fortran
436@section Allocating aligned memory in Fortran
437
438@cindex alignment
439@findex fftw_alloc_real
440@findex fftw_alloc_complex
441In order to obtain maximum performance in FFTW, you should store your
442data in arrays that have been specially aligned in memory (@pxref{SIMD
443alignment and fftw_malloc}).  Enforcing alignment also permits you to
444safely use the new-array execute functions (@pxref{New-array Execute
445Functions}) to apply a given plan to more than one pair of in/out
446arrays.  Unfortunately, standard Fortran arrays do @emph{not} provide
447any alignment guarantees.  The @emph{only} way to allocate aligned
448memory in standard Fortran is to allocate it with an external C
449function, like the @code{fftw_alloc_real} and
450@code{fftw_alloc_complex} functions.  Fortunately, Fortran 2003 provides
451a simple way to associate such allocated memory with a standard Fortran
452array pointer that you can then use normally.
453
454We therefore recommend allocating all your input/output arrays using
455the following technique:
456
457@enumerate
458
459@item
460Declare a @code{pointer}, @code{arr}, to your array of the desired type
461and dimensions.  For example, @code{real(C_DOUBLE), pointer :: a(:,:)}
462for a 2d real array, or @code{complex(C_DOUBLE_COMPLEX), pointer ::
463a(:,:,:)} for a 3d complex array.
464
465@item
466The number of elements to allocate must be an
467@code{integer(C_SIZE_T)}.  You can either declare a variable of this
468type, e.g. @code{integer(C_SIZE_T) :: sz}, to store the number of
469elements to allocate, or you can use the @code{int(..., C_SIZE_T)}
470intrinsic function. e.g. set @code{sz = L * M * N} or use
471@code{int(L * M * N, C_SIZE_T)} for an @threedims{L,M,N} array.
472
473@item
474Declare a @code{type(C_PTR) :: p} to hold the return value from
475FFTW'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.
476
477@item
478@findex c_f_pointer
479Associate your pointer @code{arr} with the allocated memory @code{p}
480using the standard @code{c_f_pointer} subroutine: @code{call
481c_f_pointer(p, arr, [...dimensions...])}, where
482@code{[...dimensions...])} are an array of the dimensions of the array
483(in the usual Fortran order). e.g. @code{call c_f_pointer(p, arr,
484[L,M,N])} for an @threedims{L,M,N} array.  (Alternatively, you can
485omit the dimensions argument if you specified the shape explicitly
486when declaring @code{arr}.)  You can now use @code{arr} as a usual
487multidimensional array.
488
489@item
490When you are done using the array, deallocate the memory by @code{call
491fftw_free(p)} on @code{p}.
492
493@end enumerate
494
495For example, here is how we would allocate an @twodims{L,M} 2d real array:
496
497@example
498  real(C_DOUBLE), pointer :: arr(:,:)
499  type(C_PTR) :: p
500  p = fftw_alloc_real(int(L * M, C_SIZE_T))
501  call c_f_pointer(p, arr, [L,M])
502  @emph{...use arr and arr(i,j) as usual...}
503  call fftw_free(p)
504@end example
505
506and here is an @threedims{L,M,N} 3d complex array:
507
508@example
509  complex(C_DOUBLE_COMPLEX), pointer :: arr(:,:,:)
510  type(C_PTR) :: p
511  p = fftw_alloc_complex(int(L * M * N, C_SIZE_T))
512  call c_f_pointer(p, arr, [L,M,N])
513  @emph{...use arr and arr(i,j,k) as usual...}
514  call fftw_free(p)
515@end example
516
517See @ref{Reversing array dimensions} for an example allocating a
518single array and associating both real and complex array pointers with
519it, for in-place real-to-complex transforms.
520
521@c -------------------------------------------------------
522@node Accessing the wisdom API from Fortran, Defining an FFTW module, Allocating aligned memory in Fortran, Calling FFTW from Modern Fortran
523@section Accessing the wisdom API from Fortran
524@cindex wisdom
525@cindex saving plans to disk
526
527As explained in @ref{Words of Wisdom-Saving Plans}, FFTW provides a
528``wisdom'' API for saving plans to disk so that they can be recreated
529quickly.  The C API for exporting (@pxref{Wisdom Export}) and
530importing (@pxref{Wisdom Import}) wisdom is somewhat tricky to use
531from Fortran, however, because of differences in file I/O and string
532types between C and Fortran.
533
534@menu
535* Wisdom File Export/Import from Fortran::
536* Wisdom String Export/Import from Fortran::
537* Wisdom Generic Export/Import from Fortran::
538@end menu
539
540@c =========>
541@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
542@subsection Wisdom File Export/Import from Fortran
543
544@findex fftw_import wisdom_from_filename
545@findex fftw_export_wisdom_to_filename
546The easiest way to export and import wisdom is to do so using
547@code{fftw_export_wisdom_to_filename} and
548@code{fftw_wisdom_from_filename}.  The only trick is that these
549require you to pass a C string, which is an array of type
550@code{CHARACTER(C_CHAR)} that is terminated by @code{C_NULL_CHAR}.
551You can call them like this:
552
553@example
554  integer(C_INT) :: ret
555  ret = fftw_export_wisdom_to_filename(C_CHAR_'my_wisdom.dat' // C_NULL_CHAR)
556  if (ret .eq. 0) stop 'error exporting wisdom to file'
557  ret = fftw_import_wisdom_from_filename(C_CHAR_'my_wisdom.dat' // C_NULL_CHAR)
558  if (ret .eq. 0) stop 'error importing wisdom from file'
559@end example
560
561Note that prepending @samp{C_CHAR_} is needed to specify that the
562literal string is of kind @code{C_CHAR}, and we null-terminate the
563string by appending @samp{// C_NULL_CHAR}.  These functions return an
564@code{integer(C_INT)} (@code{ret}) which is @code{0} if an error
565occurred during export/import and nonzero otherwise.
566
567It is also possible to use the lower-level routines
568@code{fftw_export_wisdom_to_file} and
569@code{fftw_import_wisdom_from_file}, which accept parameters of the C
570type @code{FILE*}, expressed in Fortran as @code{type(C_PTR)}.
571However, you are then responsible for creating the @code{FILE*}
572yourself.  You can do this by using @code{iso_c_binding} to define
573Fortran intefaces for the C library functions @code{fopen} and
574@code{fclose}, which is a bit strange in Fortran but workable.
575
576@c =========>
577@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
578@subsection Wisdom String Export/Import from Fortran
579
580@findex fftw_export_wisdom_to_string
581Dealing with FFTW's C string export/import is a bit more painful.  In
582particular, the @code{fftw_export_wisdom_to_string} function requires
583you to deal with a dynamically allocated C string.  To get its length,
584you must define an interface to the C @code{strlen} function, and to
585deallocate it you must define an interface to C @code{free}:
586
587@example
588  use, intrinsic :: iso_c_binding
589  interface
590    integer(C_INT) function strlen(s) bind(C, name='strlen')
591      import
592      type(C_PTR), value :: s
593    end function strlen
594    subroutine free(p) bind(C, name='free')
595      import
596      type(C_PTR), value :: p
597    end subroutine free
598  end interface
599@end example
600
601Given these definitions, you can then export wisdom to a Fortran
602character array:
603
604@example
605  character(C_CHAR), pointer :: s(:)
606  integer(C_SIZE_T) :: slen
607  type(C_PTR) :: p
608  p = fftw_export_wisdom_to_string()
609  if (.not. c_associated(p)) stop 'error exporting wisdom'
610  slen = strlen(p)
611  call c_f_pointer(p, s, [slen+1])
612  ...
613  call free(p)
614@end example
615@findex c_associated
616@findex c_f_pointer
617
618Note that @code{slen} is the length of the C string, but the length of
619the array is @code{slen+1} because it includes the terminating null
620character.  (You can omit the @samp{+1} if you don't want Fortran to
621know about the null character.) The standard @code{c_associated} function
622checks whether @code{p} is a null pointer, which is returned by
623@code{fftw_export_wisdom_to_string} if there was an error.
624
625@findex fftw_import_wisdom_from_string
626To import wisdom from a string, use
627@code{fftw_import_wisdom_from_string} as usual; note that the argument
628of this function must be a @code{character(C_CHAR)} that is terminated
629by the @code{C_NULL_CHAR} character, like the @code{s} array above.
630
631@c =========>
632@node Wisdom Generic Export/Import from Fortran,  , Wisdom String Export/Import from Fortran, Accessing the wisdom API from Fortran
633@subsection Wisdom Generic Export/Import from Fortran
634
635The most generic wisdom export/import functions allow you to provide
636an arbitrary callback function to read/write one character at a time
637in any way you want.  However, your callback function must be written
638in a special way, using the @code{bind(C)} attribute to be passed to a
639C interface.
640
641@findex fftw_export_wisdom
642In particular, to call the generic wisdom export function
643@code{fftw_export_wisdom}, you would write a callback subroutine of the form:
644
645@example
646  subroutine my_write_char(c, p) bind(C)
647    use, intrinsic :: iso_c_binding
648    character(C_CHAR), value :: c
649    type(C_PTR), value :: p
650    @emph{...write c...}
651  end subroutine my_write_char
652@end example
653
654Given such a subroutine (along with the corresponding interface definition), you could then export wisdom using:
655
656@findex c_funloc
657@example
658  call fftw_export_wisdom(c_funloc(my_write_char), p)
659@end example
660
661@findex c_loc
662@findex c_f_pointer
663The standard @code{c_funloc} intrinsic converts a Fortran
664@code{bind(C)} subroutine into a C function pointer.  The parameter
665@code{p} is a @code{type(C_PTR)} to any arbitrary data that you want
666to pass to @code{my_write_char} (or @code{C_NULL_PTR} if none).  (Note
667that you can get a C pointer to Fortran data using the intrinsic
668@code{c_loc}, and convert it back to a Fortran pointer in
669@code{my_write_char} using @code{c_f_pointer}.)
670
671Similarly, to use the generic @code{fftw_import_wisdom}, you would
672define a callback function of the form:
673
674@findex fftw_import_wisdom
675@example
676  integer(C_INT) function my_read_char(p) bind(C)
677    use, intrinsic :: iso_c_binding
678    type(C_PTR), value :: p
679    character :: c
680    @emph{...read a character c...}
681    my_read_char = ichar(c, C_INT)
682  end function my_read_char
683
684  ....
685
686  integer(C_INT) :: ret
687  ret = fftw_import_wisdom(c_funloc(my_read_char), p)
688  if (ret .eq. 0) stop 'error importing wisdom'
689@end example
690
691Your function can return @code{-1} if the end of the input is reached.
692Again, @code{p} is an arbitrary @code{type(C_PTR} that is passed
693through to your function.  @code{fftw_import_wisdom} returns @code{0}
694if an error occurred and nonzero otherwise.
695
696@c -------------------------------------------------------
697@node Defining an FFTW module,  , Accessing the wisdom API from Fortran, Calling FFTW from Modern Fortran
698@section Defining an FFTW module
699
700Rather than using the @code{include} statement to include the
701@code{fftw3.f03} interface file in any subroutine where you want to
702use FFTW, you might prefer to define an FFTW Fortran module.  FFTW
703does not install itself as a module, primarily because
704@code{fftw3.f03} can be shared between different Fortran compilers while
705modules (in general) cannot.  However, it is trivial to define your
706own FFTW module if you want.  Just create a file containing:
707
708@example
709  module FFTW3
710    use, intrinsic :: iso_c_binding
711    include 'fftw3.f03'
712  end module
713@end example
714
715Compile this file into a module as usual for your compiler (e.g. with
716@code{gfortran -c} you will get a file @code{fftw3.mod}).  Now,
717instead of @code{include 'fftw3.f03'}, whenever you want to use FFTW
718routines you can just do:
719
720@example
721  use FFTW3
722@end example
723
724as usual for Fortran modules.  (You still need to link to the FFTW
725library, of course.)
726