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