annotate src/fftw-3.3.8/doc/modern-fortran.texi @ 169:223a55898ab9 tip default

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