use FFTW, as described in @ref{Calling FFTW from Fortran}.@ref{Installation and Customization} explains how to install FFTW inyour computer system and how to adapt FFTW to your needs.  License andcopyright information is given in @ref{License and Copyright}.  Finally,we thank all the people who helped us in @ref{Acknowledgments}.@c ************************************************************@node  Tutorial, FFTW Reference, Introduction, Top@chapter Tutorial@cindex TutorialThis chapter describes the basic usage of FFTW, i.e., how to compute theFourier transform of a single array.  This chapter tells the truth, butnot the @emph{whole} truth. Specifically, FFTW implements additionalroutines and flags, providing extra functionality, that are notdocumented here.  @xref{FFTW Reference}, for more complete information.(Note that you need to compile and install FFTW before you can use it ina program.  @xref{Installation and Customization}, for the details ofthe installation.)Here, we assume a default installation of FFTW.  In some installations(particulary from binary packages), the FFTW header files and librariesare prefixed with @samp{@codeqkmoiiy} or @samp{@code{s}} to indicateversions in double or single precision, respectively.  The usage of FFTWin that case is the same, except that @code{#include} directives andlink commands must use the appropriate prefix.  @xref{Installing FFTW inboth single and double precision}, for more information.This tutorial chapter is structured as follows.  @ref{ComplexOne-dimensional Transforms Tutorial} describes the basic usage of theone-dimensional transform of complex data.  @ref{ComplexMulti-dimensional Transforms Tutorial} describes the basic usage of themulti-dimensional transform of complex data.  @ref{Real One-dimensionalTransforms Tutorial} describes the one-dimensional transform of realdata and its inverse.  Finally, @ref{Real Multi-dimensional TransformsTutorial} describes the multi-dimensional transform of real data and itsinverse.  We recommend that you read these sections in the order thatthey are presented.  We then discuss two topics in detail.  In@ref{Multi-dimensional Array Format}, we discuss the variousalternatives for storing multi-dimensional arrays in memory.  @ref{Wordsof Wisdom} shows how you can save FFTW's plans for future use.@menu* Complex One-dimensional Transforms Tutorial::  * Complex Multi-dimensional Transforms Tutorial::  * Real One-dimensional Transforms Tutorial::  * Real Multi-dimensional Transforms Tutorial::  * Multi-dimensional Array Format::  * Words of Wisdom::             @end menu@node  Complex One-dimensional Transforms Tutorial, Complex Multi-dimensional Transforms Tutorial, Tutorial, Tutorial@section Complex One-dimensional Transforms Tutorial@cindex complex one-dimensional transform@cindex complex transformThe basic usage of FFTW is simple.  A typical call to FFTW looks like:@example#include <fftw.h>...@{     fftw_complex in[N], out[N];     fftw_plan p;     ...     p = fftw_create_plan(N, FFTW_FORWARD, FFTW_ESTIMATE);     ...     fftw_one(p, in, out);     ...     fftw_destroy_plan(p);  @}@end exampleThe first thing we do is to create a @dfn{plan}, which is an object@cindex planthat contains all the data that FFTW needs to compute the FFT, using thefollowing function:@examplefftw_plan fftw_create_plan(int n, fftw_direction dir, int flags);@end example@findex fftw_create_plan@findex fftw_direction@tindex fftw_planThe first argument, @code{n}, is the size of the transform you aretrying to compute.  The size @code{n} can be any positive integer, butsizes that are products of small factors are transformed mostefficiently.  The second argument, @code{dir}, can be either@code{FFTW_FORWARD} or @code{FFTW_BACKWARD}, and indicates the directionof the transform you@ctindex FFTW_FORWARD@ctindex FFTW_BACKWARDare interested in.  Alternatively, you can use the sign of the exponentin the transform, @math{-1} or @math{+1}, which corresponds to@code{FFTW_FORWARD} or @code{FFTW_BACKWARD} respectively.  The@code{flags} argument is either @code{FFTW_MEASURE} or@cindex flags@code{FFTW_ESTIMATE}.  @code{FFTW_MEASURE} means that FFTW actually runs@ctindex FFTW_MEASUREand measures the execution time of several FFTs in order to find thebest way to compute the transform of size @code{n}.  This may take sometime, depending on your installation and on the precision of the timerin your machine.  @code{FFTW_ESTIMATE}, on the contrary, does not runany computation, and just builds a@ctindex FFTW_ESTIMATEreasonable plan, which may be sub-optimal.  In other words, if yourprogram performs many transforms of the same size and initializationtime is not important, use @code{FFTW_MEASURE}; otherwise use theestimate.  (A compromise between these two extremes exists.  @xref{Wordsof Wisdom}.)Once the plan has been created, you can use it as many times as you likefor transforms on arrays of the same size.  When you are done with theplan, you deallocate it by calling @code{fftw_destroy_plan(plan)}.@findex fftw_destroy_planThe transform itself is computed by passing the plan along with theinput and output arrays to @code{fftw_one}:@examplevoid fftw_one(fftw_plan plan, fftw_complex *in, fftw_complex *out);@end example@findex fftw_oneNote that the transform is out of place: @code{in} and @code{out} mustpoint to distinct arrays. It operates on data of type@code{fftw_complex}, a data structure with real (@code{in[i].re}) andimaginary (@code{in[i].im}) floating-point components.  The @code{in}and @code{out} arrays should have the length specified when the plan wascreated.  An alternative function, @code{fftw}, allows you toefficiently perform multiple and/or strided transforms (@pxref{FFTWReference}).@tindex fftw_complexThe DFT results are stored in-order in the array @code{out}, with thezero-frequency (DC) component in @code{out[0]}.@cindex frequencyThe array @code{in} is not modified.  Users should note that FFTWcomputes an unnormalized DFT, the sign of whose exponent is given by the@code{dir} parameter of @code{fftw_create_plan}.  Thus, computing aforward followed by a backward transform (or vice versa) results in theoriginal array scaled by @code{n}.  @xref{What FFTW Really Computes},for the definition of DFT.@cindex normalizationA program using FFTW should be linked with @code{-lfftw -lm} on Unixsystems, or with the FFTW and standard math libraries in general.@cindex linking on Unix@node Complex Multi-dimensional Transforms Tutorial, Real One-dimensional Transforms Tutorial, Complex One-dimensional Transforms Tutorial, Tutorial@section Complex Multi-dimensional Transforms Tutorial@cindex complex multi-dimensional transform@cindex multi-dimensional transformFFTW can also compute transforms of any number of dimensions(@dfn{rank}).  The syntax is similar to that for the one-dimensional@cindex ranktransforms, with @samp{fftw_} replaced by @samp{fftwnd_} (which standsfor ``@code{fftw} in @code{N} dimensions'').As before, we @code{#include <fftw.h>} and create a plan for thetransforms, this time of type @code{fftwnd_plan}:@examplefftwnd_plan fftwnd_create_plan(int rank, const int *n,                               fftw_direction dir, int flags);@end example@tindex fftwnd_plan@tindex fftw_direction@findex fftwnd_create_plan@code{rank} is the dimensionality of the array, and can be anynon-negative integer.  The next argument, @code{n}, is a pointer to aninteger array of length @code{rank} containing the (positive) sizes ofeach dimension of the array.  (Note that the array will be stored inrow-major order. @xref{Multi-dimensional Array Format}, for informationon row-major order.)  The last two parameters are the same as in@code{fftw_create_plan}.  We now, however, have an additional possibleflag, @code{FFTW_IN_PLACE}, since @code{fftwnd} supports true in-place@cindex flags@ctindex FFTW_IN_PLACE@findex fftwndtransforms.  Multiple flags are combined using a bitwise @dfn{or}(@samp{|}).  (An @dfn{in-place} transform is one in which the outputdata overwrite the input data.  It thus requires half as much memoryas---and is often faster than---its opposite, an @dfn{out-of-place}transform.)@cindex in-place transform@cindex out-of-place transformFor two- and three-dimensional transforms, FFTWND provides alternativeroutines that accept the sizes of each dimension directly, rather thanindirectly through a rank and an array of sizes.  These are otherwiseidentical to @code{fftwnd_create_plan}, and are sometimes moreconvenient:@examplefftwnd_plan fftw2d_create_plan(int nx, int ny,                               fftw_direction dir, int flags);fftwnd_plan fftw3d_create_plan(int nx, int ny, int nz,                               fftw_direction dir, int flags);@end example@findex fftw2d_create_plan@findex fftw3d_create_planOnce the plan has been created, you can use it any number of times fortransforms of the same size.  When you do not need a plan anymore, youcan deallocate the plan by calling @code{fftwnd_destroy_plan(plan)}.@findex fftwnd_destroy_planGiven a plan, you can compute the transform of an array of data bycalling:@examplevoid fftwnd_one(fftwnd_plan plan, fftw_complex *in, fftw_complex *out);@end example@findex fftwnd_oneHere, @code{in} and @code{out} point to multi-dimensional arrays inrow-major order, of the size specified when the plan was created.  Inthe case of an in-place transform, the @code{out} parameter is ignoredand the output data are stored in the @code{in} array.  The results arestored in-order, unnormalized, with the zero-frequency component in@code{out[0]}.@cindex frequencyA forward followed by a backward transform (or vice-versa) yields theoriginal data multiplied by the size of the array (i.e. the product ofthe dimensions).  @xref{What FFTWND Really Computes}, for a discussionof what FFTWND computes.@cindex normalizationFor example, code to perform an in-place FFT of a three-dimensionalarray might look like:@example#include <fftw.h>...@{     fftw_complex in[L][M][N];     fftwnd_plan p;     ...     p = fftw3d_create_plan(L, M, N, FFTW_FORWARD,                            FFTW_MEASURE | FFTW_IN_PLACE);     ...     fftwnd_one(p, &in[0][0][0], NULL);     ...     fftwnd_destroy_plan(p);  @}@end exampleNote that @code{in} is a statically-declared array, which isautomatically in row-major order, but we must take the address of thefirst element in order to fit the type expected by @code{fftwnd_one}.(@xref{Multi-dimensional Array Format}.)@node Real One-dimensional Transforms Tutorial, Real Multi-dimensional Transforms Tutorial, Complex Multi-dimensional Transforms Tutorial, Tutorial@section Real One-dimensional Transforms Tutorial@cindex real transform@cindex complex to real transform@cindex RFFTWIf the input data are purely real, you can save roughly a factor of twoin both time and storage by using the @dfn{rfftw} transforms, which areFFTs specialized for real data.  The output of a such a transform is a@dfn{halfcomplex} array, which consists of only half of the complex DFTamplitudes (since the negative-frequency amplitudes for real data arethe complex conjugate of the positive-frequency amplitudes).@cindex halfcomplex arrayIn exchange for these speed and space advantages, the user sacrificessome of the simplicity of FFTW's complex transforms.  First of all, toallow maximum performance, the output format of the one-dimensional realtransforms is different from that used by the multi-dimensionaltransforms.  Second, the inverse transform (halfcomplex to real) has theside-effect of destroying its input array.  Neither of theseinconveniences should pose a serious problem for users, but it isimportant to be aware of them.  (Both the inconvenient output formatand the side-effect of the inverse transform can be ameliorated forone-dimensional transforms, at the expense of some performance, by usinginstead the multi-dimensional transform routines with a rank of one.)The computation of the plan is similar to that for the complextransforms.  First, you @code{#include <rfftw.h>}.  Then, you create aplan (of type @code{rfftw_plan}) by calling:@examplerfftw_plan rfftw_create_plan(int n, fftw_direction dir, int flags);@end example@tindex rfftw_plan@tindex fftw_direction@findex rfftw_create_plan@code{n} is the length of the @emph{real} array in the transform (evenfor halfcomplex-to-real transforms), and can be any positive integer(although sizes with small factors are transformed more efficiently).@code{dir} is either @code{FFTW_REAL_TO_COMPLEX} or@code{FFTW_COMPLEX_TO_REAL}.@ctindex FFTW_REAL_TO_COMPLEX@ctindex FFTW_COMPLEX_TO_REALThe @code{flags} parameter is the same as in @code{fftw_create_plan}.Once created, a plan can be used for any number of transforms, and isdeallocated when you are done with it by calling@code{rfftw_destroy_plan(plan)}.@findex rfftw_destroy_planGiven a plan, a real-to-complex or complex-to-real transform is computedby calling:@examplevoid rfftw_one(rfftw_plan plan, fftw_real *in, fftw_real *out);@end example@findex rfftw_one(Note that @code{fftw_real} is an alias for the floating-point type forwhich FFTW was compiled.)  Depending upon the direction of the plan,either the input or the output array is halfcomplex, and is stored inthe following format:@cindex halfcomplex array@tex$$r_0, r_1, r_2, \ldots, r_{n/2}, i_{(n+1)/2-1}, \ldots, i_2, i_1$$@end tex@ifinfor0, r1, r2, r(n/2), i((n+1)/2-1), ..., i2, i1@end ifinfo