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