Home

FFTW User's Manual

1. FFTW IN PLACE produce a plan assuming that you want to perform the trans form in place Unlike the one dimensional transform this really performs the transform in place Note that if you want to perform in place transforms you must use a plan created with this option The use of this option has important implications for the size of the input output array see Section 3 5 2 Computing the Real Multi dimensional Transform page 30 The default mode of operation is FFTW OUT OF PLACE FFTW USE WISDOM use any wisdom that is available to help in the creation of the plan See Section 2 6 Words of Wisdom page 13 This can greatly speed the creation of plans especially with the FFTW MEASURE option FFTW ESTIMATE plans can also take advantage of wisdom to produce a more optimal plan based on past measurements than the estimation heuristic would normally generate When the FFTW MEASURE option is used new wisdom will also be generated if the current transform size is not completely understood by existing wisdom Note that the same wisdom is shared between one dimensional and multi dimensional transforms 3 5 2 Computing the Real Multi dimensional Transform include lt rfftw h gt Chapter 3 FFTW Reference 31 void rfftwnd_real_to_complex rfftwnd_plan plan int howmany fftw_real in int istride int idist fftw complex out int ostride int odist void rfftwnd complex to real rfftwnd plan plan int howmany f
2. void fftw one fftw plan plan fftw complex in fftw complex out The function fftw computes the one dimensional Fourier transform using a plan created by fftw create plan See Section 3 2 1 Plan Creation for One dimensional Transforms page 18 The function fftw one provides a simplified interface for the common case of single input array of stride 1 Arguments e plan is the plan created by fftw create plan see Section 3 2 1 Plan Creation for One dimensional Transforms page 18 e howmany is the number of transforms fftw will compute It is faster to tell FFTW to compute many transforms instead of simply calling fftw many times e in istride and idist describe the input array s There are howmany input arrays the first one is pointed to by in the second one is pointed to by in idist and so on up to in howmany 1 idist Each input array consists of complex numbers see Section 3 1 Data Types page 17 which are not necessarily contiguous in memory Specifically in 0 is the first element of the first array in istride is the second element of the first array and so on In general the i th element of the j th input array will be in position in i istride j idist Chapter 3 FFTW Reference 21 e out ostride and odist describe the output array s The format is the same as for the input array In place transforms If the plan specifies an in place transform ostride and odist are always ignored
3. 8 License and Copyright 9 Concept Index 10 Library Index ii FFTW
4. The real to complex FFTW REAL TO COMPLEX transform of a real array X of size n computes an hermitian array Y where n 1l Y n 5 bir j 0 That Y is a hermitian array is not intended to be obvious although the proof is easy The hermitian array Y is stored in halfcomplex order see Section 3 1 Data Types page 17 Currently RFFTW provides no way to compute a real to complex transform with a positive sign in the exponent Chapter 3 FFTW Reference 29 The complex to real FFTW_COMPLEX_TO_REAL transform of a hermitian array X of size n computes a real array Y where n 1i Y K JU IV m j 0 That Y is a real array is not intended to be obvious although the proof is easy The hermitian input array X is stored in halfcomplex order see Section 3 1 Data Types page 17 Currently RFFTW provides no way to compute a complex to real transform with a negative sign in the exponent Like FFTW RFFTW computes an unnormalized transform In other words applying the real to complex forward and then the complex to real backward transform will multiply the input by n 3 5 Real Multi dimensional Transforms Reference The multi dimensional real routines are generally prefixed with rfftwnd_ Programs using RFFTWND should be linked with 1rfftw 1f ftw lm on Unix systems or with the FFTW RFFTW and standard math libraries in general 3 5 1 Plan Creation for Real Multi dimensional Transforms include rfftw h
5. index ij explicitly in the convolution product loop In order to normalize the convolution we must multiply by a scale factor we can do so either before or after the inverse transform and choose the former because it obviates the necessity of an additional loop Notice the limits of the loops and the dimensions of the various arrays As with the one dimensional RFFTW an out of place FFTW_COMPLEX_TO_REAL trans form has the side effect of overwriting its input array The real to complex transform on the other hand leaves its input array untouched If you use RFFTWND for a rank one transform however this side effect does not occur Because of this fact and the simpler output format users may find the RFFTWND interface more convenient than RFFTW for one dimensional transforms However RFFTWND in one dimension is slightly slower than RFFTW because RFFTWND uses an extra buffer array internally 2 5 Multi dimensional Array Format This section describes the format in which multi dimensional arrays are stored We felt that a detailed discussion of this topic was necessary since it is often a source of confusion among users and several different formats are common Although the comments below refer to fftwnd they are also applicable to the rfftwnd routines 2 5 1 Row major Format The multi dimensional arrays passed to fftwnd are expected to be stored as a single contiguous block in row major order sometimes called C order Basical
6. page 32 In particular the z dimension of the real input data is padded to a size 2 nz 2 1 and after the transform it contains nz 2 1 complex values Some other important things to know about the real MPI transforms e As for the uniprocessor rfftwnd create plan the dimensions passed for the FFTW_ COMPLEX TO REAL plan are those of the real data In particular even when FFTW_ TRANSPOSED ORDER is used as in this case the dimensions are those of the untrans posed real output not the transposed complex input For the complex MPI trans forms on the other hand the dimensions are always those of the input array The output ordering of the transform FFTW TRANSPOSED ORDER or FFTW_TRANSPOSED_ ORDER must be the same for both forward and backward transforms This is not required in the complex case total Local size is the required size in fftw real values not fftw complex values as it is for the complex transforms local nm after transpose and local y start after transpose describe the por tion of the array after the transform that is they are indices in the complex array for an FFTW REAL TO COMPLEX transform and in the real array for an FFTW COMPLEX TO REAL transform Chapter 4 Parallel FFTW 47 e rfftwnd_mpi always expects fftw_real array arguments but of course these pointers can refer to either real or complex arrays depending upon which side of the transform you are on Just as for in place uniprocesso
7. FETW IN PLACE clle e RE NA E ae 5 fftwqmalloG uiu seeni iaiia n des red dors 35 fftw malloc hoOk o o LR Reg 35 EPTW MEASURE dasc hee hes doe x REESE Fed 4 fftw mpi se o d rex per RR EFE E Rd AT fftw mpi create plan estrei 4T fftw mpi destroy Glan 4T fftw mpi local siz68 eriiic ret 48 fftw mpi plan es doce dag ace eg d 4T fftw mpi EE 41 48 FFTW NORMAL ORDER EELER I RB 42 fftw OnQu leLererer49ed esee deeg 4 20 FETW OUT OF PLACE oed e EG E Y EET Exc x 19 fftw pleno e ee eae kde ee eee RO ae 3 18 fftw eech 2 dere ebe rege pees be Ge Sech kg FETW READ AE BE Se ERS EI 18 FFTW REAL TO COMPLEX ees 6 26 FETW SCRAMBLED INPUT 2 cor t RR 4T FETW SCRAMBLED QUTPUT ees EE e n 4T fftw test sud de aa E ehh ES 57 58 fftw threads INTC ios fence Bw deans 38 fftw threads test eee cea haere ss 39 FEIW THREADSAPE EE ie AN ANE 0 36 FETW TIMB MIW 2c I eda tee eek eas wy ocn 58 FFTW TRANSPOSED ORDER 42 44 46 48 67 FFIW USE WISDOM si bci diac desde dakar Rer dees 13 fftw2d create plan o e t s 5 22 fftw2d create plan specific 22 fftw2d mpi create plan 42 fftw3d create plan cece ee 5 22 fftw3d create plan specific 22 fftw3d f77 create plan 53 fftw3d_mpi_create_plan 42 PE CWI EE 5 24 fftw d threadS 2 EE Poke gege de qn 38 ftwnd threads One s IESEL eee gaa dees
8. FF T W comes with a shared memory implementation on top of POSIX and similar threads as well as a distributed memory implementation based on MPI We also provide an experimental parallel implemen 2 FFTW tation written in Cilk the superior programming tool of choice for discriminating hackers Olin Shivers See the Cilk home page For more information regarding FFTW see the paper The Fastest Fourier Transform in the West by M Frigo and S G Johnson which is the technical report MIT LCS TR 728 Sep 97 See also FFTW An Adaptive Software Architecture for the FFT by M Frigo and S G Johnson which appeared in the 23rd International Conference on Acoustics Speech and Signal Processing Proc ICASSP 1998 3 p 1381 The code generator is described in the paper A Fast Fourier Transform Compiler by M Frigo to appear in the Proceedings of the 1999 ACM SIGPLAN Conference on Programming Language Design and Implementation PLDI Atlanta Georgia May 1999 These papers along with the latest version of FFTW the FAQ benchmarks and other links are available at the FF TW home page The current version of FFTW incorporates many good ideas from the past thirty years of FFT literature In one way or another FFTW uses the Cooley Tukey algorithm the Prime Factor algorithm Rader s algorithm for prime sizes and the split radix algorithm with a variation due to Dan Bernstein Our code generator also produces new algor
9. The nthreads parameter specifies the number of threads of execution to use when performing the transform actually the maximum number of threads For example to parallelize a single one dimensional transform of complex data instead of calling the uniprocessor fftw_one plan in out you would call fftw_threads_ one nthreads plan in out Passing an nthreads of 1 means to use only one thread the main thread and is equivalent to calling the uniprocessor routine Passing an nthreads of 2 means that the transform is potentially parallelized over two threads and two processors if you have them and so on These are the only changes you need to make to your source code Calls to all other FFTW routines plan creation destruction wisdom etcetera are not parallelized and remain the same The same plans and wisdom are used by both uniprocessor and multi threaded transforms Your arrays are allocated and formatted in the same way and so on Programs using the parallel complex transforms should be linked with 1fftw threads lfftw 1m on Unix Programs using the parallel real transforms should be linked with lrfftw threads lfftw threads lrfftw lfftw 1m You will also need to link with whatever library is responsible for threads on your system e g 1pthread on Linux 4 1 3 How Many Threads to Use There is a fair amount of overhead involved in spawning and synchronizing threads so the optimal number of threads to use depends
10. Then for any equal or smaller power of two FFTW can create a plan with the same direction and flags quickly using the 14 FFTW precomputed wisdom Even for larger powers of two or sizes that are a power of two times some other prime factors plans will be computed more quickly than they would otherwise although some measurements still have to be made The wisdom is cumulative and is stored in a global private data structure managed internally by FFTW The storage space required is minimal proportional to the logarithm of the sizes the wisdom was generated from The wisdom can be forgotten and its associated memory freed by a call to fftw_forget_wisdom otherwise it is remembered until the program terminates It can also be exported to a file a string or any other medium using fftw_export_wisdom and restored during a subsequent execution of the program or a different program using fftw_import_wisdom these functions are described below Because wisdom is incorporated into FFTW at a very low level the same wisdom can be used for one dimensional transforms multi dimensional transforms and even the parallel extensions to FFTW Just include FFTW_USE_WISDOM in the flags for whatever plans you create i e always plan wisely Plans created with the FFTW_ESTIMATE plan can use wisdom but cannot generate it only FFTW_MEASURE plans actually produce wisdom Also plans can only use wisdom generated from plans created with the same
11. and the local y extent will be given by Local ny after transpose and local y start after transpose Otherwise the output has the same dimensions and layout as the input 44 FFTW For instance suppose you want to transform three dimensional data of size nx x ny x nz Then the current process will store a subset of this data of size local_nx x ny x nz where the x indices correspond to the range local_x_start to local_x_start local_nx 1 in the real i e logical array If fftwnd mpi is called with FFTW_TRANSPOSED_ORDER output then the result will be a ny x nx x nz array of which a 1ocal ny after transpose x nx x nz subset is stored on the current process corresponding to y values starting at local y start after transpose The following is an example of allocating such a three dimensional array array local_ data before the transform and initializing it to some function f x y z fftwnd mpi local sizes plan amp local nx amp local x start amp local ny after transpose amp local y start after transpose amp total local size local data fftw complex malloc sizeof fftw complex total local size for x 0 x lt local nx x for y 0 y lt ny y for z 0 z lt nz z local data x ny y nz z f x local x start y z Some important things to remember e Although the local data is of dimensions Local nx x ny x nz in the above example do not allocate the array to be of size
12. because the transform data is distributed across the processes involved in the transform It is discussed in detail by the next section see Section 4 2 3 MPI Data Layout page 43 The actual computation of the transform is performed by the function fftwnd mpi which differs somewhat from its uniprocessor equivalent and is described by void fftwnd mpi fftwnd mpi plan p int n fields fftw complex local data fftw complex work fftwnd mpi output order output order There are several things to notice here e First of all all fftw mpi transforms are in place the output is in the Local data parameter and there is no need to specify FFTW IN PLACE in the plan flags e The MPI transforms also only support a limited subset of the hownany stride dist functionality of the uniprocessor routines the n fields parameter is equivalent to howmany n fields stride n fields and dist 1 Conceptually the n fields pa rameter allows you to transform an array of contiguous vectors each with length n_ fields n fields is 1 if you are only transforming a single ordinary array e The work parameter is an optional workspace If it is not NULL it should be exactly the same size as the 1ocal data array If it is provided FFTW is able to use the built in MPI Alltoall primitive for often greater efficiency at the expense of extra storage space e Finally the last parameter specifies whether the output data has the same ordering as the input
13. over In particular the array is divided according to the rows first dimension of the data each process gets a subset of the rows of the data This is sometimes called a slab decom position One consequence of this is that you can t take advantage of more processors than you have rows e g 64x64x64 matrix can at most use 64 processors This isn t usually much of a limitation however as each processor needs a fair amount of data in order for the parallel computation benefits to outweight the communications costs Below the first dimension of the data will be referred to as x and the second dimension as y FFTW supplies a routine to tell you exactly how much data resides on the current process void fftwnd_mpi_local_sizes fftwnd_mpi_plan p int local_nx int local_x_start int local_ny_after_transpose int local_y_start_after_transpose int total local size Given a plan p the other parameters of this routine are set to values describing the required data layout described below total local size is the number of fftw complex elements that you must allocate for your local data and workspace if you choose This value should of course be multiplied by n fields if that parameter to fftwnd mpi is not 1 The data on the current process has 1ocal nx rows starting at row local x start If fftwnd mpi is called with FFTW TRANSPOSED ORDER output then y will be the first dimen sion of the output
14. plus one element of the last dimension of the data from the ordinary complex transform We could have instead taken half of any other dimension but implementation turns out to be simpler if the last contiguous dimension is used Since the complex data is slightly larger than the real data some complications arise for in place transforms In this case the final dimension of the real data must be padded with extra values to accommodate the size of the complex data two extra if the last dimension is even and one if it is odd That is the last dimension of the real data must physically contain 2 n4 2 1 fftw real values exactly enough to hold the complex data This physical array size does not however change the logical array size only n values are actually stored in the last dimension and n is the last dimension passed to rfftwnd create plan 3 5 4 Strides in In place RFFTWND The fact that the input and output datatypes are different for rfftwnd complicates the meaning of the stride and dist parameters of in place transforms are they in units of fftw_real or fftw_complex elements When reading the input they are interpreted in units of the datatype of the input data When writing the output the istride and idist are translated to the output datatype s units in one of two ways corresponding to the two most common situations in which stride and dist parameters are useful Below we refer to these translated parame
15. rfftwnd plan rfftwnd create plan int rank const int n fftw direction dir int flags rfftwnd plan rfftw2d create plan int nx int ny fftw direction dir int flags rfftwnd plan rfftw3d create plan int nx int ny int nz fftw direction dir int flags The function rfftwnd create plan creates a plan which is a data structure containing all the information that rfftwnd needs in order to compute a multi dimensional real Fourier transform You can create as many plans as you need but only one plan for a given array size is required a plan can be reused many times The functions rfftw2d create plan and rfftw3d create plan are optional alternative interfaces to rfftwnd create plan for two and three dimensions respectively rfftwnd create plan returns a valid plan or NULL if for some reason the plan can t be created This can happen if the arguments are invalid in some way e g if rank 0 Arguments e rank is the dimensionality of the arrays to be transformed It can be any non negative integer e nisa pointer to an array of rank integers giving the size of each dimension of the arrays to be transformed Note that these are always the dimensions of the real arrays the 30 FFTW complex arrays have different dimensions see Section 3 5 3 Array Dimensions for Real Multi dimensional Transforms page 32 These sizes which must be positive integers correspond to the dimensions of row major arrays i e n 0 is the size
16. 0 we found that the crossover point at which 2 threads became beneficial for complex trans forms was about 4k points while 4 threads became beneficial at 8k points 4 1 4 Using Multi threaded FFTW in a Multi threaded Program It is perfectly possible to use the multi threaded FFTW routines from a multi threaded program e g have multiple threads computing multi threaded transforms simultaneously If you have the processors more power to you However the same restrictions apply as for the uniprocessor FFTW routines see Section 3 8 Thread safety page 36 In particular you should recall that you may not create or destroy plans in parallel 4 1 5 Tips for Optimal Threading Not all transforms are equally well parallelized by the multi threaded FF TW routines This is merely a consequence of laziness on the part of the implementors and is not inherent to the algorithms employed Mainly the limitations are in the parallel one dimensional transforms The things to avoid if you want optimal parallelization are as follows 4 1 6 Parallelization deficiencies in one dimensional transforms e Large prime factors can sometimes parallelize poorly Of course you should avoid these anyway if you want high performance e Single in place transforms don t parallelize completely Multiple in place transforms i e howmany gt 1 are fine Again you should avoid these in any case if you want high performance as they require transforming to
17. Destroying a Multi dimensional Plan include lt fftw h gt void fftwnd_destroy_plan fftwnd_plan plan The function fftwnd destroy plan frees the plan plan and releases all the memory associated with it After destruction a plan is no longer valid 3 3 4 What FFTWND Really Computes The conventions that we follow for the multi dimensional transform are analogous to those for the one dimensional transform In particular the forward transform has a negative sign in the exponent and neither the forward nor the backward transforms will perform any normalization Computing the backward transform of the forward transform will multiply the array by the product of its dimensions The output is in order and the zeroth element of the output is the amplitude of the zero frequency component The exact mathematical definition of our multi dimensional transform follows Let X be a d dimensional complex array whose elements are X j1 j2 ja where 0 lt js lt ng for all s 1 2 d Let also w e V for all s 1 2 d The forward transform computes a complex array Y whose structure is the same as that of X defined by nj1 1 n3 1 ngq 1 Yen X Y Me jaor wg ees m ji 0 jo 0 ja 0 The backward transform computes nj1 1 n3 1 na l Vita cl 9 X Xl jz le a eu ji 0 j2 0 ja 0 Computing the forward transform followed by the backward transform will multiply the array by I2 4 na 26 FFTW 3 4
18. If out is NULL out is ignored too Otherwise out is interpreted as a pointer to an array of n complex numbers that FFTW will use as temporary space to perform the in place computation out is used as scratch space and its contents destroyed In this case out must be an ordinary array whose elements are contiguous in memory no striding The function fftw_one transforms a single contiguous input array to a contiguous output array By definition the call fftw_one plan in out is equivalent to fftw plan 1 in 1 0 out 1 0 3 2 4 Destroying a One dimensional Plan include lt fftw h gt void fftw_destroy_plan fftw_plan plan The function fftw_destroy_plan frees the plan plan and releases all the memory asso ciated with it After destruction a plan is no longer valid 3 2 5 What FFTW Really Computes In this section we define precisely what FFTW computes Please be warned that differ ent authors and software packages might employ different conventions than FFTW does The forward transform of a complex array X of size n computes an array Y where n 1i Y Se j 0 The backward transform computes n 1 Y 5 Ke j 0 FFTW computes an unnormalized transform that is the equation ZFFT FFT X nX holds In other words applying the forward and then the backward transform will multiply the input by n An FFTW FORWARD transform corresponds to a sign of 1 in the exponent of the DFT Note also that we
19. N call fftwnd_f77_create_plan plan 3 n FFTW_FORWARD FFTW_ESTIMATE FFTW_IN_PLACE call fftwnd_f77_one plan arr 0 call fftwnd f77 destroy plan plan Instead of calling fftwnd f77 create plan plan 3 n we could also have called fftw3d f77 create plan plan L M N Note that we pass the array dimensions in the natural order also note that the last argument to fftwnd f77 is ignored since the transform is FFTW IN PLACE 54 FFTW To transform a one dimensional real array in Fortran you might do double precision in out dimension in N out N integer plan call rfftw f77 create plan plan N FFTW REAL TO COMPLEX t FFTW ESTIMATE call rfftw f77 one plan in out call rfftw f77 destroy plan plan To transform a two dimensional real array out of place you might use the following double precision in double complex out dimension in M N out M 2 1 N integer plan call rfftw2d f77 create plan plan M N FFTW REAL TO COMPLEX FFTW_ESTIMATE call rfftwnd_f77_one_real_to_complex plan in out call rfftwnd f77 destroy plan plan Important Notice that it is the first dimension of the complex output array that is cut in half in Fortran rather than the last dimension as in C This is a consequence of the wrapper routines reversing the order of the array dimensions passed to FFTW so that the Fortran program can use its ordinary column major order Chapter 6 Installation and Customization 55 6 Installatio
20. Real One dimensional Transforms Reference The one dimensional real routines are generally prefixed with rfftw Programs using RFFTW should be linked with 1rfftw 1fftw 1m on Unix systems or with the RFFTW the FFTW and the standard math libraries in general 3 4 1 Plan Creation for Real One dimensional Transforms include lt rfftw h gt rfftw_plan rfftw_create_plan int n fftw_direction dir int flags rfftw_plan rfftw_create_plan_specific int n fftw_direction dir int flags fftw_real in int istride fftw real out int ostride The function rfftw_create_plan creates a plan which is a data structure containing all the information that rfftw needs in order to compute the 1D real Fourier transform You can create as many plans as you need but only one plan for a given array size is required a plan can be reused many times rfftw_create_plan returns a valid plan or NULL if for some reason the plan can t be created In the default installation this cannot happen but it is possible to configure RFFTW in such a way that some input sizes are forbidden and RFFTW cannot create a plan The rfftw_create_plan_specific variant takes as additional arguments specific in put output arrays and their strides For the last four arguments you should pass the arrays and strides that you will eventually be passing to rfftw The resulting plans will be optimized for those arrays and strides although they may be used on other arrays as
21. Y i5 i4 Y ni t1 ne do n4 ta holds for all 0 i lt n Because of this symmetry Y is stored in the peculiar way described in Section 3 5 3 Array Dimensions for Real Multi dimensional Transforms page 32 Complex to real backward transform Let X be a d dimensional complex array whose elements are X j1 j2 Ja where 0 js lt n for all s 1 2 d The array X must be hermitian that is the identity X j1 j2 ja X ni ji na ja na Aal must hold for all 0 j n Moreover X must be stored in memory in the peculiar way described in Section 3 5 3 Array Dimensions for Real Multi dimensional Transforms page 32 Let w e 7V V s for all s 1 2 d The complex to real transform computes a real array Y whose structure is the same as that of X defined by nj1 1 n3 1 ng 1 Y h i id SD So Y EE ji 0 j2 0 ja 0 That Y is real is not meant to be obvious although the proof is easy Computing the forward transform followed by the backward transform will multiply the array by T5 na 3 6 Wisdom Reference 3 6 1 Exporting Wisdom include lt fftw h gt void fftw_export_wisdom void emitter char c void void data void fftw export wisdom to file FILE output_file char fftw export wisdom to string void These functions allow you to export all currently accumulated wisdom in a form from which it can be later imported and restor
22. all gone fftw import wisdom from string wisdom string wisdom is back fftw free wisdom string deallocate it since we re done Chapter 3 FFTW Reference 17 3 FFTW Reference This chapter provides a complete reference for all sequential i e one processor FF TW functions We first define the data types upon which FFTW operates that is real com plex and halfcomplex numbers see Section 3 1 Data Types page 17 Then in four sec tions we explain the FFTW program interface for complex one dimensional transforms see Section 3 2 One dimensional Transforms Reference page 18 complex multi dimensional transforms see Section 3 3 Multi dimensional Transforms Reference page 22 and real one dimensional transforms see Section 3 4 Real One dimensional Transforms Reference page 26 real multi dimensional transforms see Section 3 5 Real Multi dimensional Trans forms Reference page 29 Section 3 6 Wisdom Reference page 34 describes the wisdom mechanism for exporting and importing plans Finally Section 3 7 Memory Allocator Ref erence page 35 describes how to change FFTW s default memory allocator For parallel transforms See Chapter 4 Parallel FFTW page 37 3 1 Data Types The routines in the FFTW package use three main kinds of data types Real and complex numbers should be already known to the reader We also use the term halfcomplex to describe complex arrays in a special packed format
23. and should be at least several hundred times the resolution of your clock The default constants are on the conservative side and may cause FF TW to take longer than necessary when you create a plan Set FFTW TIME MIN to whatever is appropriate on your system be sure to set the right FFTW TIME MIN there are several definitions in fftw int h corresponding to different platforms and timers As an aid in checking the resolution of your clock you can use the tests fftw test program with the t option c f tests README Remember the mere fact that your Chapter 6 Installation and Customization 59 clock reports times in say picoseconds does not mean that it is actually accurate to that resolution 6 6 Generating your own code If you know that you will only use transforms of a certain size say powers of 2 and want to reduce the size of the library you can reconfigure FF TW to support only those sizes you are interested in You may even generate code to enable efficient transforms of a size not supported by the default distribution The default distribution supports transforms of any size but not all sizes are equally fast The default installation of FFTW is best at handling sizes of the form 273 5 7411 13 where e f is either 0 or 1 and the other exponents are arbitrary Other sizes are computed by means of a slow general purpose routine However if you have an application that requires fast transforms of size say 17 ther
24. computes the actual transform as dictated by the plan The plan can be reused as many times as needed In typical high performance applications many transforms of the same size are computed and consequently a relatively expensive initialization of this sort is acceptable On the other hand if you need a single transform of a given size the one time cost of the planner becomes significant For this case FF TW provides fast planners based on heuristics or on previously computed plans The pattern of planning execution applies to all four operation modes of FFTW that is I one dimensional complex transforms FFTW II multi dimensional complex transforms FFTWND III one dimensional transforms of real data RFFTW IV multi dimensional transforms of real data RFFTWND Each mode comes with its own planner and executor Besides the automatic performance adaptation performed by the planner it is also pos sible for advanced users to customize FFTW for their special needs As distributed FF TW works most efficiently for arrays whose size can be factored into small primes 2 3 5 and 7 and uses a slower general purpose routine for other factors FF TW however comes with a code generator that can produce fast C programs for any particular array size you may care about For example if you need transforms of size 513 19 3 you can customize FFTW to support the factor 19 efficiently FFTW can exploit multiple processors if you have them
25. data FFTW NORMAL ORDER or if it is transposed FFTW TRANSPOSED ORDER Leaving the data transposed results in significant performance improvements due to a saved communication step needed to un transpose the data Specifically the first two dimensions of the array are transposed as is described in more detail by the next section Chapter 4 Parallel FFTW 43 The output of fftwnd mpi is identical to that of the corresponding uniprocessor trans form In particular you should recall our conventions for normalization and the sign of the transform exponent The same plan can be used to compute many transforms of the same size After you are done with it you should deallocate it by calling fftwnd mpi destroy plan Important The FFTW MPI routines must be called in the same order by all processes involved in the transform You should assume that they all are blocking as if each contained a call to MPI Barrier Programs using the FFTW MPI routines should be linked with 1fftw mpi lfftw 1m on Unix in addition to whatever libraries are required for MPI 4 2 3 MPI Data Layout The transform data used by the MPI FFTW routines is distributed a distinct portion of it resides with each process involved in the transform This allows the transform to be parallelized for example over a cluster of workstations each with its own separate memory so that you can take advantage of the total memory of all the processors you are parallelizing
26. direction and flags For example a size 42 FFTW_BACKWARD transform will not use wisdom produced by a size 42 FFTW_FORWARD transform The only exception to this rule is that FFTW_ESTIMATE plans can use wisdom from FFTW_MEASURE plans 2 6 1 Caveats in Using Wisdom For in much wisdom is much grief and he that increaseth knowledge increaseth sorrow Ecclesiastes 1 18 There are pitfalls to using wisdom in that it can negate FFTW s ability to adapt to changing hardware and other conditions For example it would be perfectly possible to export wisdom from a program running on one processor and import it into a program running on another processor Doing so however would mean that the second program would use plans optimized for the first processor instead of the one it is running on It should be safe to reuse wisdom as long as the hardware and program binaries remain unchanged Actually the optimal plan may change even between runs of the same binary on identical hardware due to differences in the virtual memory environment etcetera Users seriously interested in performance should worry about this problem too It is likely that if the same wisdom is used for two different program binaries even running on the same machine the plans may be sub optimal because of differing code alignments It is therefore wise to recreate wisdom every time an application is recompiled The more the underlying hardware and software changes between the cr
27. document ing the different parallel libraries that we have provided Users calling FFTW from a multi threaded program should also consult Section 3 8 Thread safety page 36 The FFTW package currently contains three parallel transform implementations that leverage the uniprocessor FFTW code e The first set of routines utilizes shared memory threads for parallel one and multi dimensional transforms of both real and complex data Any program using FFTW can be trivially modified to use the multi threaded routines This code can use any com mon threads implementation including POSIX threads POSIX threads are available on most Unix variants including Linux These routines are located in the threads directory and are documented in Section 4 1 Multi threaded FFTW page 37 e The mpi directory contains multi dimensional transforms of real and complex data for parallel machines supporting MPI It also includes parallel one dimensional transforms for complex data The main feature of this code is that it supports distributed memory transforms so it runs on everything from workstation clusters to massively parallel supercomputers More information on MPI can be found at the MPI home page The FFTW MPI routines are documented in Section 4 2 MPI FFTW page 40 e We also have an experimental parallel implementation written in Cilk a C like parallel language developed at MIT and currently available for several SMP platforms For more informati
28. function performs any one time initialization required to use threads on your system It returns zero if successful and a non zero value if there was an error in which case something is seriously wrong and you should probably exit the program Third when you want to actually compute the transform you should use one of the following transform routines instead of the ordinary FF TW functions fftw threads nthreads plan howmany in istride idist out ostride odist fftw threads one nthreads plan in out fftwnd threads nthreads plan howmany in istride idist out ostride odist fftwnd threads one nthreads plan in out rfftw threads nthreads plan howmany in istride idist out ostride odist rfftw threads one nthreads plan in out rfftwnd threads real to complex nthreads plan howmany in istride idist out ostride odist rfftwnd threads one real to complex nthreads plan in out rfftwnd threads complex to real nthreads plan howmany in Chapter 4 Parallel FFTW 39 istride idist out ostride odist rfftwnd_threads_one_real_to_complex nthreads plan in out rfftwnd_threads_one_complex_to_real nthreads plan in out All of these routines take exactly the same arguments and have exactly the same effects as their uniprocessor counterparts i e without the _threads except that they take one extra parameter nthreads of type int before the normal parameters
29. have for fftw 1 The basic problem is the resolution of the clock FF TW needs to run for a certain time for the clock to be reliable 20 FFTW 3 2 2 Discussion on Specific Plans We recommend the use of the specific planners even in cases where you will be transform ing arrays different from those passed to the specific planners as they confer the following advantages e The resulting plans will be optimized for your specific arrays and strides This may or may not make a significant difference but it certainly doesn t hurt The ordinary planner does its planning based upon a stride one temporary array that it allocates e Less intermediate storage is required during the planning process The ordinary plan ner uses O N temporary storage where N is the maximum dimension while it is cre ating the plan e For multi dimensional transforms new parameters become accessible for optimization by the planner Since multi dimensional arrays can be very large we don t dare to allocate one in the ordinary planner for experimentation This prevents us from doing certain optimizations that can yield dramatic improvements in some cases On the other hand note that the specific planner destroys the contents of the in and out arrays 3 2 3 Computing the One dimensional Transform include lt fftw h gt void fftw fftw_plan plan int howmany fftw complex sin int istride int idist fftw complex out int ostride int odist
30. input array s There are howmany multi dimensional input arrays the first one is pointed to by in the second one is pointed to by in idist and so on up to in howmany 1 idist Each multi dimensional input array consists of complex numbers see Section 3 1 Data Types page 17 stored in row major format see Section 2 5 Multi dimensional Array Format page 11 which are not necessarily contiguous in memory Specifically in 0 is the first element of the first array in istride is the second element of the first array and so on In general the i th element of the j th input array will be in position in i istride j idist Note that here i refers to an index into the row major format for the multi dimensional array rather than an index in any particular dimension In place transforms For plans created with the FFTW IN PLACE option the trans form is computed in place the output is returned in the in array using the same strides etcetera as were used in the input Chapter 3 FFTW Reference 25 e out ostride and odist describe the output array s The format is the same as for the input array In place transforms These parameters are ignored for plans created with the FFTW_IN_PLACE option The function fftwnd one transforms a single contiguous input array to a contiguous output array By definition the call fftwnd one plan in out is equivalent to fftwnd plan 1 in 1 0 out 1 0 3 3 3
31. local data contains the portion of the output local to the current process see below To find out what portion of the array is stored local to the current process you call the following routine void fftw mpi local sizes fftw mpi plan p int local n int local start int local n after transform int local start after transform int total local size total local size is the number of fftw complex elements you should actually allocate for local data and work local n and Local start indicate that the current process stores local n elements corresponding to the indices local start to local start local n 1 in the real array After the transform the process may store a different portion of the array The portion of the data stored on the process after the transform is given by local n after transform and local start after transform This data is exactly the same as a contiguous segment of the corresponding uniprocessor transform output i e an in order sequence of sequential frequency bins Note that if you compute both a forward and a backward transform of the same size the local sizes are guaranteed to be consistent That is the local size after the forward transform will be the same as the local size before the backward transform and vice versa Programs using the FFTW MPI routines should be linked with 1fftw mpi lfftw 1m on Unix in addition to whatever libraries are required for MPI 4 2 6 MPI Tips There ar
32. local_nx ny nz Use total local size instead e The amount of data on each process will not necessarily be the same in fact Local nx may even be zero for some processes For example suppose you are doing a 6x6 transform on four processors There is no way to effectively use the fourth processor in a slab decomposition so we leave it empty Proof left as an exercise for the reader e All arrays are of course in row major order see Section 2 5 Multi dimensional Array Format page 11 e If you want to compute the inverse transform of the output of fftwnd mpi the dimen sions of the inverse transform are given by the dimensions of the output of the forward transform For example if you are using FFTW TRANSPOSED ORDER output in the above example then the inverse plan should be created with dimensions ny x nx x nz e The data layout only depends upon the dimensions of the array not on the plan so you are guaranteed that different plans for the same size or inverse plans will use the same consistent data layouts 4 2 4 Usage of MPI FFTW for Real Multi dimensional Transforms MPI transforms specialized for real data are also available similiar to the uniprocessor rfftwnd transforms Just as in the uniprocessor case the real data MPI functions gain roughly a factor of two in speed and save a factor of two in space at the expense of more complicated data formats in the calling program Before reading this section you should
33. of the dimension whose indices vary most slowly and so on See Section 2 5 Multi dimensional Array Format page 11 for more information See Section 3 4 1 Plan Creation for Real One dimensional Transforms page 26 for more information regarding optimal array sizes nx and ny in rfftw2d create plan are positive integers specifying the dimensions of the rank 2 array to be transformed i e they specify that the transform will operate on nx x ny arrays in row major order where nx is the number of rows and ny is the number of columns nx ny and nz in rfftw3d create plan are positive integers specifying the dimensions of the rank 3 array to be transformed i e they specify that the transform will operate on nx x ny x nz arrays in row major order dir is the direction of the desired transform either FFTW REAL TO COMPLEX or FFTW COMPLEX TO REAL corresponding to FFTW FORWARD or FFTW_BACKWARD respectively flags is a boolean OR of zero or more of the following FFTW MEASURE this flag tells FF TW to find the optimal plan by actually computing several FFTs and measuring their execution time FFTW ESTIMATE do not run any FFT and provide a reasonable plan for a RISC processor with many registers If neither FFTW ESTIMATE nor FFTW MEASURE is provided the default is FFTW ESTIMATE FFTW OUT OF PLACE produce a plan assuming that the input and output arrays will be distinct this is the default
34. the wrapper subroutine The reason for this is that some Fortran implementations seem to have trouble with C function return values e When performing one dimensional FFTW IN PLACE transforms you don t have the op tion of passing NULL for the out argument since there is no way to pass NULL from Fortran Therefore when performing such transforms you must allocate and pass a contiguous scratch array of the same size as the transform Note that for in place multi dimensional r fftwnd transforms the out argument is ignored so you can pass anything for that parameter E Technically Fortran 77 identifiers are not allowed to have more than 6 characters nor may they contain underscores Any compiler that enforces this limitation doesn t deserve to link to FFTW 52 FFTW e The wrapper routines expect multi dimensional arrays to be in column major order which is the ordinary format of Fortran arrays They do this transparently and cost lessly simply by reversing the order of the dimensions passed to FFTW but this has one important consequence for multi dimensional real complex transforms discussed below In general you should take care to use Fortran data types that correspond to i e are the same size as the C types used by FFTW If your C and Fortran compilers are made by the same vendor the correspondence is usually straightforward i e integer corresponds to int real corresponds to float etcetera Such simple correspo
35. upon the size of the transform as well as on the number of processors you have As a general rule you don t want to use more threads than you have processors Using more threads will work but there will be extra overhead with no benefit In fact if the problem size is too small you may want to use fewer threads than you have processors You will have to experiment with your system to see what level of parallelization is best for your problem size Useful tools to help you do this are the test programs that are automatically compiled along with the threads libraries fftw threads test and rfftw_ threads test in the threads subdirectory These take the same arguments as the other FFTW test programs see tests README except that they also take the number of threads to use as a first argument and report the parallel speedup in speed tests For example fftw threads test 2 s 128x128 1 There is one exception when performing one dimensional in place transforms the out parameter is always ignored by the multi threaded routines instead of being used as a workspace if it is non NULL as in the uniprocessor routines The multi threaded routines always allocate their own workspace the size of which depends upon the number of threads 40 FFTW will benchmark complex 128x128 transforms using two threads and report the speedup relative to the uniprocessor transform For instance on a 4 processor 200MHz Pentium Pro system running Linux 2 2
36. well Note the contents of the in and out arrays are destroyed by the specific planner the initial contents are ignored so the arrays need not have been initialized See Section 3 2 2 Discussion on Specific Plans page 20 for a discussion on specific plans Arguments e nis the size of the transform It can be any positive integer RFFTW is best at handling sizes of the form 2 3 5 7411 13 where e f is either 0 or 1 and the other exponents are arbitrary Other sizes are computed by means of a slow general purpose routine reducing to O n performance for prime sizes It is possible to customize RFFTW for different array sizes See Chapter 6 In stallation and Customization page 55 for more information Transforms whose sizes are powers of 2 are especially fast If you have large prime factors it may be faster to switch over to the complex FFTW routines which have O nlog n performance even for prime sizes we don t know of a similar algorithm specialized for real data unfortunately e dir is the direction of the desired transform either FFTW REAL TO COMPLEX or FFTW_ COMPLEX TO REAL corresponding to FFTW FORWARD or FFTW_BACKWARD respectively The etymologically correct spelling would be frftw_ but it is hard to remember Chapter 3 FFTW Reference 27 e flags is a boolean OR of zero or more of the following FFTW_MEASURE this flag tells RFFTW to find the optimal plan by actually com puti
37. with POSIX threads avail able on most Unix variants including Linux Solaris threads BeOS threads tested on 38 FFTW BeOS DR8 2 Mach C threads reported to work by users and Win32 threads reported to work by users There is also untested code to use MacOS MP threads We also support using OpenMP or SGI MP compiler directives to launch threads enabled by using with openmp or with sgimp in addition to enable threads This is especially useful if you are employing that sort of directive in your own code in order to minimize conflicts If you have a shared memory machine that uses a different threads API it should be a simple matter of programming to include support for it see the file fftw threads int h for more detail SMP hardware is not required although of course you need multiple processors to get any benefit from the multithreaded transforms 4 1 2 Usage of Multi threaded FFTW Here it is assumed that the reader is already familiar with the usage of the uniprocessor FFTW routines described elsewhere in this manual We only describe what one has to change in order to use the multi threaded routines First instead of including fftw h or rfftw h you should include the files lt fftw_ threads h gt or rfftw threads h respectively Second before calling any FFTW routines you should call the function int fftw threads init void This function which should only be called once probably in your mainO
38. your compiler appropriately Be sure to enable the highest possible level of optimization in your compiler By default FFTW is compiled for double precision transforms To work in single preci sion rather than double precision define the symbol FFTW ENABLE FLOAT in fftw h in the fftw directory and re compile FFTW Chapter 6 Installation and Customization 57 These libraries should be linked with any program that uses the corresponding trans forms The required header files fftw h and rfftw h are located in the fftw and rfftw directories respectively you may want to put them with the libraries or wherever header files normally go on your system FFTW includes test programs fftw_test and rfftw_test in the tests directory These are compiled and linked like any program using FFTW except that they use addi tional header files located in the fftw and rfftw directories so you will need to set your compiler include paths appropriately fftw_test is compiled from fftw_test c and test_main c while rfftw_test is compiled from rfftw_test c and test_main c When you run these programs you will be prompted interactively for various possible tests to perform see also tests README for more information 6 3 Installing FF TW in both single and double precision It is often useful to install both single and double precision versions of the FF TW libraries on the same machine and we provide a convenient mechanism for achieving this on Unix
39. 11 20 24 27 315 32 42 T thread safety oi 3 eed oie eR eor okii 36 40 tbreads 2 vid tate br ts 1 36 37 56 timer customization of 58 Tutorial zen n erepeqet Euren sree 3 ELSE Sein ee 13 19 23 30 34 wisdom import and export sssssseuuuu 14 wisdom problems with 14 Chapter 10 Library Index 10 Library Index F bn P TTE 20 fftw threads ipee ere Ier ease wad 38 fftw threads one scere 38 FETW_BACKWARD NN dE use e hx tet eop ER 4 fftw comlnpl x c er e np heo ne ANT FEIW CHMPLEX orea mnp REG RACER e READS D ROSS 18 FFTW_COMPLEX_TO_REAL 2 02040 EEN 008 6 26 46 fftw create Dan ed EEN Ee 3 18 fftw create plan specific 18 fftw destroy plan sess 4 21 fftw direction 3 4 6 18 22 29 FEFTW ENABLE PLOAT aoaesacesa sr R9 IRR 18 FETW ESTIMATE 4elawdog m ra e RR ERI n Gene 4 fftw export wisdom sss 13 34 fftw export wisdom to file 14 34 fftw export wisdom to string 15 34 fftw f 7 create plan o ss ga es 53 fftw f77 destroy Glan 53 fftwfi Oon6 onsiosg caes eei goes up des 51 53 fftw f77 threads one 53 fftw forget wisdOom osse ee 35 FEIW PUBWABD A ut ku REX REA dimes She BOSS 4 fftw free hook fs godsend Ree Ene 35 fftw import wisdom SEELEN SEELEN 13 35 fftw import wisdom from file 14 35 fftw import wisdom from string 15 35
40. 2 40 42 51 installation coe che er Rer LE CERA MES 55 L linking on Unix 4 7 39 43 47 48 DISP deg ege geg eier 15 61 load balancing EE 48 M malloG s see ged codes bte e 12 35 UR RE Rr P bEICEEREERSTE PRESS 59 monadic programming 05 59 MPD ine ohana dence xe feuis 1 37 40 56 MPI Alte 41 42 48 MPT Barriere Eer Peces predam 43 MPI COMM WORLD 5 rank bm 42 47 MPI EE ER 23 0 2vLd4eeme Le RUE 42 MPD UNDG itis asa ie are ae wearer ned Icio dide 42 multi dimensional transform 4 22 29 N MSL AS asi sites psa duh Sob en VEL tnd dys e pe 42 47 nerd readable test 15 normalization 4 5 7 8 11 21 29 42 number of threads 00 0c e nipi eae 39 O out of place transform ee eee 5 P pad dilg ie den aries edge EP 8 32 46 parallel transform mee REED Berl ac 1 37 Pentium hack o2 2202 2505 Se EE ge R 57 Plaise steel need utei Aided che eed ees 1 3 42 II TE tud der ph Aaa Bea Re ed aieo eius i 66 PANG ae Ae pror EES 4 real multi dimensional transform 7 29 real number 42 cca b Gare e Sege dE Ee e real transfor ee EEN eee eee ee 6 26 SN ENEE E 6 26 REPTWND 19 4 Bb d EE 29 rfftwnd array format 8 32 46 54 TOW MAJOL EE 11 23 44 S saving plans to disk see os C rv EE PRA 13 FFTW slab decomposition 00 00 eee eee 43 specific planners 02 e eee eee eee 20 16016 er RR i
41. 38 fftwnd create plan Ea 4 22 fftwnd create plan specific 22 fftwnd destroy plan 5 25 fftwnd f77 create plan 58 ffitwid f77 destroy plan xy 59 fftwnd 71 ONE eaba ees RE RE a ERES 53 FEC WN DEEG 42 fftwnd_mpi_create_plan 42 fftwnd mpi destroy plan 43 fftwnd mpi local sizes 43 fftwnd mpi melen aisen ge de 42 fftwnrd Oone ss bs wees dar e RESO UpPERS heed 5 24 fftwnd plam ii sre r eg vem uS 4 22 G PONE EE 59 61 R EEL EE 2 rfftw threadS ood sop AR See wee bes 38 rfftw threads ne oae xd er 38 tfftw create plan ie M eI ER tri 6 26 rfftw create plan specific 26 rfftw destroy plam obo enep tas 6 28 rfftw f77 6reate plan d a eta 54 rfftw 77 destroy plam e rwn 54 ffftw f f One Joop esewd eR ET 54 rfftw Dpi te6Sto c ie d ace ER weeds 41 49 tfft NS o irena retiri que p diy 6 27 rffftw plani so aged er RLERRU Epe 6 26 Ifftw test 0 210 PUER a Rer rS 5T rfftw threads test sange cba e ees 39 rfftw2d create plan 7 29 47 rfftw2d f77 create plan 54 rfftw3d mpi create plan 45 rfftw3d create plan 7 29 rfftwnd mpi i o arip iar EA nade 46 rfftwnd mpi destroy plan 46 rfftwnd mpi local sizes 45 68 rfftwnd_threads_complex_to_real 39 rff
42. Chapter 4 Parallel FFTW 45 definitely understand how to call the uniprocessor rfftwnd functions and also the complex MPI FF TW functions The following is an example of a program using rfftwnd_mpi It computes the size nx x ny x nz transform of a real function f x y z multiplies the imaginary part by 2 for fun then computes the inverse transform We ll also use FFTW TRANSPOSED ORDER output for the transform and additionally supply the optional workspace parameter to rfftwnd_mpi just to add a little spice include lt rfftw_mpi h gt int main int argc char argv 1 const int nX ace Dy senn DZ As int local nx local x start local ny after transpose local y start after transpose total local size int x y Z rfftwnd mpi plan plan iplan fftw real data work fftw complex cdata MPI Init amp argc amp argv create the forward and backward plans plan rfftw3d mpi create plan MPI COMM WORLD nx ny nz FFTW REAL TO COMPLEX FFTW ESTIMATE iplan rfftw3d mpi create plan MPI COMM WORLD dim s of REAL data gt nx ny nz FFTW COMPLEX TO REAL FFTW ESTIMATE rfftwnd mpi local sizes plan amp local nx amp local x start amp local ny after transpose amp local y start after transpose amp total local size data fftw real malloc sizeof fftw real total local size workspace is the same size as the data work fftw real malloc sizeof f
43. D page 32 e out ostride and odist describe the output array s The format is the same as that for the input array See below for a discussion of the dimensions of the output array for real and complex data In place transforms These parameters are ignored for plans created with the FFTW_IN_PLACE option The function rfftwnd_one transforms a single contiguous input array to a contiguous output array By definition the call rfftwnd_one_ plan in out is equivalent to rfftwnd_ plan 1 in 1 0 out 1 0 3 5 3 Array Dimensions for Real Multi dimensional Transforms The output of a multi dimensional transform of real data contains symmetries that in principle make half of the outputs redundant see Section 3 5 6 What RFFTWND Really Computes page 33 In practice it is not possible to entirely realize these savings in an efficient and understandable format Instead the output of the rfftwnd transforms is slightly over half of the output of the corresponding complex transform We do not pack the data in any way but store it as an ordinary array of fftw complex values In fact this data is simply a subsection of what would be the array in the corresponding complex transform Specifically for a real transform of dimensions n x ny x x na the complex data is an n4 X n3 X X n4 2 4 1 array of ftw complex values in row major order with the division rounded down That is we only store the lower half
44. FFTW User s Manual For version 2 1 5 16 March 2003 Matteo Frigo Steven G Johnson Copyright 1997 1999 Massachusetts Institute of Technology Permission is granted to make and distribute verbatim copies of this manual provided the copyright notice and this permission notice are preserved on all copies Permission is granted to copy and distribute modified versions of this manual under the con ditions for verbatim copying provided that the entire resulting derived work is distributed under the terms of a permission notice identical to this one Permission is granted to copy and distribute translations of this manual into another lan guage under the above conditions for modified versions except that this permission notice may be stated in a translation approved by the Free Software Foundation Chapter 1 Introduction 1 1 Introduction This manual documents version 2 1 5 of FFTW the Fastest Fourier Transform in the West FFTW is a comprehensive collection of fast C routines for computing the discrete Fourier transform DFT in one or more dimensions of both real and complex data and of arbitrary input size FFTW also includes parallel transforms for both shared and distributed memory systems We assume herein that the reader is already familiar with the properties and uses of the DFT that are relevant to her application Otherwise see e g The Fast Fourier Transform by E O Brigham Prentice Hall Englewood Cliffs NJ 1974 Ou
45. FTW 026 ccancveveccedsenes EQ ES 37 4 1 Multi threaded HHIW ccc eens 37 4 1 1 Installation and Supported Hardware Software 37 4 1 2 Usage of Multi threaded FFTW 38 4 1 8 How Many Threads to Une 39 4 1 4 Using Multi threaded FFTW in a Multi threaded e EE 40 4 1 5 Tips for Optimal Threading 40 4 1 6 Parallelization deficiencies in one dimensional transforms A aia o1 ed rhe Ca aep io Y d OR es 40 42 MPE EET S cett esee ssh E S ini 40 4 2 1 MPI FFTW Installation 40 4 2 2 Usage of MPI FFTW for Complex Multi dimensional T ra sfOFll8 22 52 aut dec cedet a orte Rn e ve 41 4 2 3 MPI Data Dayot meh 43 4 2 4 Usage of MPI FFTW for Real Multi dimensional Transforms o 2223 ppt od fad eae LE eds 44 4 2 5 Usage of MPI FFTW for Complex One dimensional Transforms cei zies deae peterent depend is 4T 420 MPIPS 4 oO IER and et Postea is 48 5 Calling FFTW from Fortran 51 5 1 Wrapper Routnes 0 0 cece cee eee eee neces 51 5 2 FFTW Constants in Furtran eee eee 52 5 3 Fortran Examples 0 00000 cece cece eee es 52 6 Installation and Customization 55 6 4 Installation on Umts esses 55 6 2 Installation on non Unix Systems sese else 56 6 3 Installing FFTW in both single and double precision 57 6 4 gcc and Pentium hacks 2 seo bg esi 57 6 5 Customizing the mer 58 6 6 Generating your own code 59 7 Acknowledements eese 61 FFTW
46. Transform Compiler by M Frigo which is available from the FF T W home page and will appear in the Proceedings of the 1999 ACM SIGPLAN Conference on Programming Language Design and Implementation PLDI 60 FFTW Chapter 7 Acknowledgments 61 7 Acknowledgments Matteo Frigo was supported in part by the Defense Advanced Research Projects Agency DARPA under Grants N00014 94 1 0985 and F30602 97 1 0270 and by a Digital Equip ment Corporation Fellowship Steven G Johnson was supported in part by a DoD NDSEG Fellowship an MIT Karl Taylor Compton Fellowship and by the Materials Research Sci ence and Engineering Center program of the National Science Foundation under award DMR 9400334 Both authors were also supported in part by their respective girlfriends by the letters Q and R and by the number 12 We are grateful to SUN Microsystems Inc for its donation of a cluster of 9 8 processor Ultra HPC 5000 SMPs 24 Gflops peak These machines served as the primary platform for the development of earlier versions of FFTW We thank Intel Corporation for donating a four processor Pentium Pro machine We thank the Linux community for giving us a decent OS to run on that machine The genfft program was written using Objective Caml a dialect of ML Objective Caml is asmall and elegant language developed by Xavier Leroy The implementation is available from ftp inria fr in the directory lang caml light We used versions 1 07 an
47. _to_complex functions and FFTW COMPLEX TO REAL plans must be used with the complex to real functions It is an error to mismatch the plan direction and the transform function e howmany is the number of transforms to be computed e in istride and idist describe the input array s There are howmany input arrays the first one is pointed to by in the second one is pointed to by in idist and so on up to in howmany 1 idist Each input array is stored in row major format see Section 2 5 Multi dimensional Array Format page 11 and is not neces sarily contiguous in memory Specifically in 0 is the first element of the first array in istride is the second element of the first array and so on In general the i th element of the j th input array will be in position in i istride j idist Note that here i refers to an index into the row major format for the multi dimensional array rather than an index in any particular dimension The dimensions of the arrays are different for real and complex data and are discussed in more detail below see Section 3 5 3 Array Dimensions for Real Multi dimensional Transforms page 32 In place transforms For plans created with the FFTW IN PLACE option the trans form is computed in place the output is returned in the in array The meaning 32 FFTW of the stride and dist parameters in this case is subtle and is discussed below see Section 3 5 4 Strides in In place RFFTWN
48. a scratch array and copying back e Single real complex rfftw transforms don t parallelize completely This is unfortu nate but parallelizing this correctly would have involved a lot of extra code and a much larger library You still get some benefit from additional processors but if you have a very large number of processors you will probably be better off using the parallel complex fftw transforms Note that multi dimensional real transforms or multiple one dimensional real transforms are fine 4 2 MPI FFTW This section describes the MPI FFTW routines for distributed memory and shared memory machines supporting MPI Message Passing Interface The MPI routines are significantly different from the ordinary FFTW because the transform data here are dis tributed over multiple processes so that each process gets only a portion of the array Currently multi dimensional transforms of both real and complex data as well as one dimensional transforms of complex data are supported 4 2 1 MPI FFTW Installation The FFTW MPI library code is all located in the mpi subdirectoy of the FFTW package along with source code for test programs On Unix systems the FFTW MPI libraries and Chapter 4 Parallel FFTW 41 header files can be automatically configured compiled and installed along with the unipro cessor FF TW libraries simply by including enable mpi in the flags to the configure script see Section 6 1 Installation on Unix pag
49. ariety in place transforms are used for the forward FFTs and an out of place transform is used for the inverse transform include lt rfftw h gt fftw real a M 2 N 2 1 b M 2 N 2 1 c M N fftw complex A B C M N 2 1 rfftwnd plan p pinv fftw real scale 1 0 M N int i j rfftw2d create plan M N FFTW REAL TO COMPLEX FFTW ESTIMATE FFTW IN PLACE pinv rfftw2d create plan M N FFTW COMPLEX TO REAL FFTW ESTIMATE P aliases for accessing complex transform outputs A fftw complex amp a 0 0 B fftw complex amp b 0 0 for i 0 i lt M i for j 0 j lt N j f ali j b i j rfftwnd_one_real_to_complex p amp a 0 0 NULL rfftwnd_one_real_to_complex p amp b 0 0 NULL for i 0 i lt M i for j 0 j lt N 2 1 j int ij i N 2 1 j C i j re Alijl re B ijl re A ijl im B ijl im scale C ilLjl im A ijl re B ijl im A ij im B ijl re scale inverse transform to get c the convolution of a and b this has the side effect of overwriting C rfftwnd_one_complex_to_real pinv amp C 0 0 amp c 0 0 rfftwnd_destroy_plan p rfftwnd_destroy_plan pinv We access the complex outputs of the in place transforms by casting each real array to a fftw_complex pointer Because this is a flat pointer we have to compute the row major Chapter 2 Tutorial 11
50. ay is a rank three N x M x L matrix in column major order you should pass the dimensions of the array as if it were anLxMxN matrix which it is from the perspective of fftwnd This is done for you 12 FFTW automatically by the FFTW Fortran wrapper routines see Chapter 5 Calling FFTW from Fortran page 51 2 5 3 Static Arrays in C Multi dimensional arrays declared statically that is at compile time not necessarily with the static keyword in C are already in row major order You don t have to do any thing special to transform them See Section 2 2 Complex Multi dimensional Transforms Tutorial page 4 for an example of this sort of code 2 5 4 Dynamic Arrays in C Often especially for large arrays it is desirable to allocate the arrays dynamically at runtime This isn t too hard to do although it is not as straightforward for multi dimensional arrays as it is for one dimensional arrays Creating the array is simple using a dynamic allocation routine like malloc allocate an array big enough to store N fftw complex values where N is the product of the sizes of the array dimensions ie the total number of complex values in the array For example here is code to allocate a 5x12x27 rank 3 array fftw complex an array an array fftw complex malloc 5 12 27 sizeof fftw complex Accessing the array elements however is more tricky you can t simply use multiple applications of the operator like
51. be generated if the current transform size is not completely understood by existing wisdom e in out istride ostride only for rfftw_create_plan_specific see correspond ing arguments in the description of rfftw See Section 3 4 2 Computing the Real One dimensional Transform page 27 In particular the out and ostride parameters have the same special meaning for FFTW_IN_PLACE transforms as they have for rfftw 3 4 2 Computing the Real One dimensional Transform include lt rfftw h gt void rfftw rfftw_plan plan int howmany fftw_real in int istride int idist fftw_real out int ostride int odist void rfftw one rfftw plan plan fftw_real sin fftw_real out The function rfftw computes the Real One dimensional Fourier Transform using a plan created by rfftw create plan see Section 3 4 1 Plan Creation for Real One dimensional Transforms page 26 The function rfftw one provides a simplified interface for the common case of single input array of stride 1 Important When invoked for an out of place FFTW COMPLEX TO REAL transform the input array is overwritten with scratch values by these routines The input array is not modified for FFTW REAL TO COMPLEX transforms 28 FFTW Arguments e plan is the plan created by rfftw_create_plan see Section 3 4 1 Plan Creation for Real One dimensional Transforms page 26 e howmany is the number of transforms rfftw will compute It is faster to tell RFFTW to comp
52. by calling fftw mpi destroy plan plan If you don t care about the ordering of the input or output data of the transform you can include FFTW SCRAMBLED INPUT and or FFTW SCRAMBLED OUTPUT in the flags These save some communications at the expense of having the input and or output reordered in an undocumented way For example if you are performing an FF T based convolution you might use FFTW SCRAMBLED OUTPUT for the forward transform and FFTW SCRAMBLED INPUT for the inverse transform The transform itself is computed by void fftw mpi fftw mpi plan p int n fields fftw complex local data fftw complex work n fields as in fftwnd mpi is equivalent to howmany n fields stride n fields and dist 1 and should be 1 when you are computing the transform of a single array local data contains the portion of the array local to the current process described below work is either NULL or an array exactly the same size as local data in the latter case The 1D transforms require much more communication All the communication in our FFT routines takes the form of an all to all communication the multi dimensional transforms require two all to all communications or one if you use FFTW TRANSPOSED ORDER while the one dimensional transforms require three or two if you use scrambled input or output 48 FFTW FFTW can use the MPI_Alltoall communications primitive which is usually faster at the expense of extra storage Upon return
53. ct of the dimensions See Section 3 3 4 What FFTWND Really Computes page 25 for a discussion of what FFTWND computes For example code to perform an in place FFT of a three dimensional array might look like include lt fftw h gt 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 amp in 0 0 0 NULL fftwnd destroy plan p Note that in is a statically declared array which is automatically in row major order but we must take the address of the first element in order to fit the type expected by fftwnd_one See Section 2 5 Multi dimensional Array Format page 11 6 FFTW 2 3 Real One dimensional Transforms Tutorial If the input data are purely real you can save roughly a factor of two in both time and storage by using the rfftw transforms which are FFTs specialized for real data The output of a such a transform is a halfcomplex array which consists of only half of the complex DFT amplitudes since the negative frequency amplitudes for real data are the complex conjugate of the positive frequency amplitudes In exchange for these speed and space advantages the user sacrifices some of the sim plicity of FFTW s complex transforms First of all to allow maximum performance the output format of the one dimensional real transforms is different from that used by the multi dimensional transforms Second the inverse transform hal
54. d 2 00 of the software In previous releases of FFTW genfft was written in Caml Light by the same authors An even earlier implementation of genfft was written in Scheme but Caml is definitely better for this kind of application FFTW uses many tools from the GNU project including automake texinfo and libtool Prof Charles E Leiserson of MIT provided continuous support and encouragement This program would not exist without him Charles also proposed the name codelets for the basic FFT blocks Prof John D Joannopoulos of MIT demonstrated continuing tolerance of Steven s extra curricular computer science activities Steven s chances at a physics degree would not exist without him Andrew Sterian contributed the Windows timing code Didier Miras reported a bug in the test procedure used in FFTW 1 2 We now use a completely different test algorithm by Funda Ergun that does not require a separate FFT program to compare against Wolfgang Reimer contributed the Pentium cycle counter and a few fixes that help porta bility Ming Chang Liu uncovered a well hidden bug in the complex transforms of FFTW 2 0 and supplied a patch to correct it The FFTW FAQ was written in bfnn Bizarre Format With No Name and formatted using the tools developed by Ian Jackson for the Linux FAQ We are especially thankful to all of our users for their continuing support feedback and interest during our development of FFTW 62 FFTW Chap
55. d is stored in the following format To T1 T2 5 Tn 25 n 1 2 15 5 2 1 Here rx is the real part of the kth output and i is the imaginary part We follow here the C convention that integer division is rounded down e g 7 2 3 For a halfcomplex array hc the kth component has its real part in hc k and its imaginary part in hc n k with the exception of k 0 or n 2 the latter only if n is even in these two cases the imaginary part is zero due to symmetries of the real complex transform and is not stored Thus the transform of n real values is a halfcomplex array of length n and vice versa This is actually only half of the DFT spectrum of the data Although the other half can be obtained by complex conjugation it is not required by many applications such as convolution and filtering 1 The output for the multi dimensional rfftw is a more conventional array of fftw complex values but the format here permitted us greater efficiency in one dimension Chapter 2 Tutorial 7 Like the complex transforms the RFFTW transforms are unnormalized so a forward followed by a backward transform or vice versa yields the original data scaled by the length of the array n Let us reiterate here our warning that an FFTW COMPLEX TO REAL transform has the side effect of destroying its halfcomplex input The FFTW_REAL_TO_COMPLEX transform however leaves its real input untouched just as you would hope As an
56. ds in order to compute the 1D Fourier transform You can create as many plans as you need but only one plan for a given array size is required a plan can be reused many times fftw create plan returns a valid plan or NULL if for some reason the plan can t be created In the default installation this cannot happen but it is possible to configure FFTW in such a way that some input sizes are forbidden and FFTW cannot create a plan The fftw create plan specific variant takes as additional arguments specific in put output arrays and their strides For the last four arguments you should pass the arrays and strides that you will eventually be passing to fftw The resulting plans will be optimized for those arrays and strides although they may be used on other arrays as well Note the contents of the in and out arrays are destroyed by the specific planner the initial contents are ignored so the arrays need not have been initialized Arguments e nis the size of the transform It can be any positive integer FFTW is best at handling sizes of the form 273 5 7711 13 where e f is either 0 or 1 and the other exponents are arbitrary Other sizes are computed by means of a slow general purpose routine which nevertheless retains O nlog n performance Chapter 3 FFTW Reference 19 even for prime sizes It is possible to customize FFTW for different array sizes See Chapter 6 Installation and Customization page 55 for more infor
57. e 55 The only requirement of the FFTW MPI code is that you have the standard MPI 1 1 or later libraries and header files installed on your system A free implementation of MPI is available from the MPICH home page Previous versions of the FFTW MPI routines have had an unfortunate tendency to expose bugs in MPI implementations The current version has been largely rewritten and hopefully avoids some of the problems If you run into difficulties try passing the optional workspace to r fftwnd mpi see below as this allows us to use the standard and hopefully well tested MPI Alltoall primitive for communications Please let us know fftwOfftw org how things work out Several test programs are included in the mpi directory The ones most useful to you are probably the fftw mpi test and rfftw mpi test programs which are run just like an ordinary MPI program and accept the same parameters as the other FF TW test programs c f tests README For example mpirun params fftw mpi test r O will run non terminating complex transform correctness tests of random dimensions They can also do performance benchmarks 4 2 2 Usage of MPI FFTW for Complex Multi dimensional Transforms Usage of the MPI FFTW routines is similar to that of the uniprocessor FFTW We assume that the reader already understands the usage of the uniprocessor FF TW routines described elsewhere in this manual Some familiarity with MPI is also helpful A typical progra
58. e is a way to generate specialized code to handle that The directory gensrc contains all the programs and scripts that were used to generate FFTW In particular the program gensrc genfft ml was used to generate the code that FFTW uses to compute the transforms We do not expect casual users to use it genfft is a rather sophisticated program that generates directed acyclic graphs of FFT algorithms and performs algebraic simplifications on them genfft is written in Objective Caml a dialect of ML Objective Caml is described at http pauillac inria fr ocaml and can be downloaded from from ftp ftp inria fr lang caml light If you have Objective Caml installed you can type sh bootstrap sh in the top level directory to re generate the files If you change the gensrc config file you can optimize FFTW for sizes that are not currently supported efficiently say 17 or 19 We do not provide more details about the code generation process since we do not expect that users will need to generate their own code However feel free to contact us at fftwOfftw org if you are interested in the subject You might find it interesting to learn Caml and or some modern programming techniques that we used in the generator including monadic programming especially if you heard the rumor that Java and object oriented programming are the latest advancement in the field The internal operation of the codelet generator is described in the paper A Fast Fourier
59. e several things you should consider in order to get the best performance out of the MPI FFTW routines First if possible the first and second dimensions of your data should be divisible by the number of processes you are using If only one can be divisible then you should choose the first dimension This allows the computational load to be spread evenly among the processes and also reduces the communications complexity and overhead In the one dimensional transform case the size of the transform should ideally be divisible by the square of the number of processors Second you should consider using the FFTW_TRANSPOSED_ORDER output format if it is not too burdensome The speed gains from communications savings are usually substantial Third you should consider allocating a workspace for r fftw nd _mpi as this can often but not always improve performance at the cost of extra storage Fourth you should experiment with the best number of processors to use for your problem There comes a point of diminishing returns when the communications costs outweigh the computational benefits The fftw mpi test program can output helpful performance benchmarks It accepts the same parameters as the uniprocessor test programs c f tests README and is run like an ordinary MPI program For example mpirun np 3 An FFT is particularly hard on communications systems as it requires an all to all communication which is more or less the worst po
60. eation of wisdom and its use the greater grows the risk of sub optimal plans 2 6 2 Importing and Exporting Wisdom void fftw export wisdom to file FILE output_file fftw status fftw import wisdom from file FILE sinput file fftw export wisdom to file writes the wisdom to output file which must be a file open for writing fftw import wisdom from file reads the wisdom from input file Chapter 2 Tutorial 15 which must be a file open for reading and returns FFTW_SUCCESS if successful and FFTW_ FAILURE otherwise In both cases the file is left open and must be closed by the caller It is perfectly fine if other data lie before or after the wisdom in the file as long as the file is positioned at the beginning of the wisdom data before import char fftw export wisdom to string void fftw status fftw import wisdom from string const char input string fftw export wisdom to string allocates a string exports the wisdom to it in NULL terminated format and returns a pointer to the string If there is an error in allocating or writing the data it returns NULL The caller is responsible for deallocating the string with fftw free when she is done with it fftw import wisdom from string imports the wisdom from input string returning FEI SUCCESS if successful and FFTW FAILURE otherwise Exporting wisdom does not affect the store of wisdom Imported wisdom supplements the current store rather than replacing it except when there is c
61. ed even during a separate run of the program See Section 2 6 Words of Wisdom page 13 The current store of wisdom is not affected by calling any of these routines fftw export wisdom exports the wisdom to any output medium as specified by the callback function emitter emitter is a putc like function that writes the character c to some output its second parameter is the data pointer passed to fftw export wisdom For convenience the following two wrapper routines are provided fftw export wisdom to file writes the wisdom to the current position in output file which should be open with write permission Upon exit the file remains open and is positioned at the end of the wisdom data fftw export wisdom to string returns a pointer to a NULL terminated string holding the wisdom data This string is dynamically allocated and it is the responsibility of the caller to deallocate it with fftw free when it is no longer needed All of these routines export the wisdom in the same format which we will not document here except to say that it is LISP like ASCII text that is insensitive to white space Chapter 3 FFTW Reference 35 3 6 2 Importing Wisdom include lt fftw h gt fftw status fftw import wisdom int get_input void void data fftw status fftw import wisdom from file FILE input file fftw status fftw import wisdom from string const char input string These functions import wisdom into a program from data sto
62. example here is an outline of how you might use RFFTW to compute the power spectrum of a real array i e the squares of the absolute values of the DFT amplitudes include lt rfftw h gt fftw real in N out N power_spectrum N 2 1 rfftw plan p int k p rfftw create plan N FFTW REAL TO COMPLEX FFTW ESTIMATE rfftw one p in out power spectrum 0 out O out O DC component for k 1 k lt N 1 2 k k lt N 2 rounded up power spectrum k out k out k out N k out N k if N 2 0 N is even power spectrum N 2 out N 2 out N 2 Nyquist freq rfftw destroy plan p Programs using RFFTW should link with 1rfftw 1fftw lm on Unix or with the FFTW RFFTW and math libraries in general 2 4 Real Multi dimensional Transforms Tutorial FFTW includes multi dimensional transforms for real data of any rank As with the one dimensional real transforms they save roughly a factor of two in time and storage over complex transforms of the same size Also as in one dimension these gains come at the expense of some increase in complexity the output format is different from the one dimensional RFFTW and is more similar to that of the complex FFTW and the inverse complex to real transforms have the side effect of overwriting their input data To use the real multi dimensional transforms you first include rfftw h and then create a plan for the size and direction of trans
63. f the current transform size is not completely understood by existing wisdom Note that the same wisdom is shared between one dimensional and multi dimensional transforms e in out istride ostride only for the create plan specific variants see cor responding arguments in the description of fftwnd See Section 3 3 2 Computing the Multi dimensional Transform page 24 3 3 2 Computing the Multi dimensional Transform include lt fftw h gt void fftwnd fftwnd_plan plan int howmany fftw complex sin int istride int idist fftw complex out int ostride int odist void fftwnd one fftwnd plan p fftw complex in fftw complex out The function fftwnd computes one or more multi dimensional Fourier Transforms us ing a plan created by fftwnd create plan see Section 3 3 1 Plan Creation for Multi dimensional Transforms page 22 Note that the plan determines the rank and dimensions of the array to be transformed The function fftwnd one provides a simplified interface for the common case of single input array of stride 1 Arguments e plan is the plan created by fftwnd create plan see Section 3 3 1 Plan Creation for Multi dimensional Transforms page 22 In the case of two and three dimensional transforms it could also have been created by fftw2d create plan or fftw3d create plan respectively howmany is the number of multi dimensional transforms fftwnd will compute in istride and idist describe the
64. fcomplex to real has the side effect of destroying its input array Neither of these inconveniences should pose a serious problem for users but it is important to be aware of them Both the inconvenient output format and the side effect of the inverse transform can be ameliorated for one dimensional transforms at the expense of some performance by using instead the multi dimensional transform routines with a rank of one The computation of the plan is similar to that for the complex transforms First you include rfftw h Then you create a plan of type rfftw plan by calling rfftw plan rfftw create plan int n fftw direction dir int flags n is the length of the real array in the transform even for halfcomplex to real trans forms and can be any positive integer although sizes with small factors are transformed more efficiently dir is either FFTW REAL TO COMPLEX or FFTW COMPLEX TO REAL The flags parameter is the same as in fftw create plan Once created a plan can be used for any number of transforms and is deallocated when you are done with it by calling rfftw destroy plan plan Given a plan a real to complex or complex to real transform is computed by calling void rfftw one rfftw plan plan fftw real in fftw real out Note that fftw_real is an alias for the floating point type for which FFTW was com piled Depending upon the direction of the plan either the input or the output array is halfcomplex an
65. form that you are interested in rfftwnd plan rfftwnd create plan int rank const int n fftw direction dir int flags The first two parameters describe the size of the real data not the halfcomplex data which will have different dimensions The last two parameters are the same as those for rfftw create plan Just as for fftwnd there are two alternate versions of this routine rfftw2d create plan and rfftw3d create plan that are sometimes more convenient for two and three dimensional transforms Also as in fftwnd rfftwnd supports true in place transforms specified by including FFTW IN PLACE in the flags 8 FFTW Once created a plan can be used for any number of transforms and is deallocated by calling rfftwnd_destroy_plan plan Given a plan the transform is computed by calling one of the following two routines void rfftwnd_one_real_to_complex rfftwnd_plan plan fftw real sin fftw complex out void rfftwnd one complex to real rfftwnd plan plan fftw complex sin fftw real out As is clear from their names and parameter types the former function is for FFTW REAL TO COMPLEX transforms and the latter is for FFTW COMPLEX TO REAL transforms We could have used only a single routine since the direction of the transform is encoded in the plan but we wanted to correctly express the datatypes of the parameters The latter routine as we discuss elsewhere has the side effect of overwriting its input except when the ran
66. ftw complex in int istride int idist fftw real out int ostride int odist void rfftwnd one real to complex rfftwnd plan p fftw real in fftw complex out void rfftwnd one complex to real rfftwnd plan p fftw complex in fftw real out These functions compute the real multi dimensional Fourier Transform using a plan cre ated by rfftwnd create plan see Section 3 5 1 Plan Creation for Real Multi dimensional Transforms page 29 Note that the plan determines the rank and dimensions of the ar ray to be transformed The rfftwnd_one_ functions provide a simplified interface for the common case of single input array of stride 1 Unlike other transform routines in FFTW we here use separate functions for the two directions of the transform in order to correctly express the datatypes of the parameters Important When invoked for an out of place FFTW COMPLEX TO REAL transform with rank gt 1 the input array is overwritten with scratch values by these routines The input array is not modified for FFTW REAL TO COMPLEX transforms or for FFTW COMPLEX TO REAL with rank Arguments e planisthe plan created by rfftwnd create plan see Section 3 5 1 Plan Creation for Real Multi dimensional Transforms page 29 In the case of two and three dimensional transforms it could also have been created by rfftw2d create plan or rfftw3d_ create plan respectively FFTW REAL TO COMPLEX plans must be used with the real
67. ftw real total local size initialize data to f x y z for x 0 x local nx x for y 0 y lt ny ty for z 0 z nz z data x ny y 2 nz 2 1 z f x local x start y z 46 FFTW Now compute the forward transform rfftwnd_mpi plan 1 data work FFTW_TRANSPOSED_ORDER the data is now complex so typecast a pointer cdata fftw_complex data multiply imaginary part by 2 for fun note that the data is transposed for y 0 y lt local_ny_after_transpose y for x 0 x lt nx x for z 0 z lt nz 2 1 z cdata y nx x nz 241 z im 2 0 Finally compute the inverse transform the result is transposed back to the original data layout rfftwnd mpi iplan 1 data work FFTW TRANSPOSED ORDER free data free work rfftwnd mpi destroy plan plan rfftwnd mpi destroy plan iplan MPI_Finalize d There s a lot of stuff in this example but it s all just what you would have guessed right We replaced all the ftwnd_mpi functions by rfftwnd_mpi but otherwise the parameters were pretty much the same The data layout distributed among the processes just like for the complex transforms see Section 4 2 3 MPI Data Layout page 43 but in addition the final dimension is padded just like it is for the uniprocessor in place real transforms see Section 3 5 3 Array Dimensions for Real Multi dimensional Transforms
68. gn of whose exponent is given by the dir parameter of fftw create plan Thus computing a forward followed by a backward transform or vice versa results in the original array scaled by n See Section 3 2 5 What FFTW Really Computes page 21 for the definition of DFT A program using FFTW should be linked with 1fftw 1m on Unix systems or with the FFTW and standard math libraries in general 2 2 Complex Multi dimensional Transforms Tutorial FFTW can also compute transforms of any number of dimensions rank The syntax is similar to that for the one dimensional transforms with fftw_ replaced by fftwnd_ which stands for fftw in N dimensions As before we include lt fftw h gt and create a plan for the transforms this time of type fftwnd plan fftwnd plan fftwnd create plan int rank const int n fftw direction dir int flags rank is the dimensionality of the array and can be any non negative integer The next argument n is a pointer to an integer array of length rank containing the positive sizes of each dimension of the array Note that the array will be stored in row major order See Chapter 2 Tutorial 5 Section 2 5 Multi dimensional Array Format page 11 for information on row major order The last two parameters are the same as in fftw_create_plan We now however have an additional possible flag FFTW_IN_PLACE since fftwnd supports true in place transforms Multiple flags are combined using a b
69. have something like the following to transform a one dimensional complex array fftw complex in N out N fftw plan plan Chapter 5 Calling FFTW from Fortran 53 plan fftw create plan N FFTW FORWARD FFTW ESTIMATE fftw one plan in out fftw destroy plan plan In Fortran you use the following to accomplish the same thing double complex in out dimension in N out N integer plan call fftw f77 create plan plan N FFTW FORWARD FFTW ESTIMATE call fftw f77 one plan in out call fftw f77 destroy plan plan Notice how all routines are called as Fortran subroutines and the plan is returned via the first argument to fftw f77 create plan Important these examples assume that integer is the same size as a pointer and may need modification on a 64 bit machine See Section 5 1 Wrapper Routines page 51 above To do the same thing but using 8 threads in parallel see Section 4 1 Multi threaded FFTW page 37 you would simply replace the call to fftw f77 one with call fftw f77 threads one 8 plan in out To transform a three dimensional array in place with C you might do fftw complex arr L M N fftwnd_plan plan int n 3 L M N plan fftwnd_create_plan 3 n FFTW_FORWARD FFTW_ESTIMATE FFTW_IN_PLACE fftwnd_one plan arr 0 fftwnd destroy plan plan In Fortran you would use this instead double complex arr dimension arr L M N integer n dimension n 3 integer plan n 1 L n 2 M n 3
70. he scheme depends on the parity of n Let n 4 In this case we have Y Bel X Y Re X1 Yo Re X2 18 FFTW and Y Im X Let now n 5 In this case we have Y Re Xo Yi Re X1 Y Re X5 Y Im X2 and Y Im X By default the type fftw real equals the C type double To work in single precision rather than double precision define the symbol FFTW ENABLE FLOAT in fftw h and then recompile the library On Unix systems you can instead use configure enable float at installation time see Chapter 6 Installation and Customization page 55 In version 1 of FFTW the data types were called FFTW REAL and FFTW COMPLEX We changed the capitalization for consistency with the rest of FFTW s conventions The old names are still supported but their use is deprecated 3 2 One dimensional Transforms Reference The one dimensional complex routines are generally prefixed with fftw Programs using FFTW should be linked with 1fftw 1m on Unix systems or with the FFTW and standard math libraries in general 3 2 1 Plan Creation for One dimensional Transforms include lt fftw h gt fftw_plan fftw_create_plan int n fftw_direction dir int flags fftw_plan fftw_create_plan_specific int n fftw_direction dir int flags fftw_complex in int istride fftw complex out int ostride The function fftw create plan creates a plan which is a data structure containing all the information that fftw nee
71. here See Chapter 3 FFTW Reference page 17 for more complete informa tion Note that you need to compile and install FFTW before you can use it in a program See Chapter 6 Installation and Customization page 55 for the details of the installation Here we assume a default installation of FFTW In some installations particulary from binary packages the FFTW header files and libraries are prefixed with d or s to indicate versions in double or single precision respectively The usage of FF TW in that case is the same except that include directives and link commands must use the appropriate prefix See Section 6 3 Installing FFTW in both single and double precision page 57 for more information This tutorial chapter is structured as follows Section 2 1 Complex One dimensional Transforms Tutorial page 3 describes the basic usage of the one dimensional transform of complex data Section 2 2 Complex Multi dimensional Transforms Tutorial page 4 de scribes the basic usage of the multi dimensional transform of complex data Section 2 3 Real One dimensional Transforms Tutorial page 6 describes the one dimensional trans form of real data and its inverse Finally Section 2 4 Real Multi dimensional Transforms Tutorial page 7 describes the multi dimensional transform of real data and its inverse We recommend that you read these sections in the order that they are presented We then discuss two topics in detail In Secti
72. ike this technique and want to maximize convenience in accessing the array but still want to pass the array to FFTW you can use a hybrid method Allocate the array as one contiguous block but also declare an array of arrays of pointers that point to appropriate places in the block That sort of trick is beyond the scope of this documentation for more information on multi dimensional arrays in C see the comp lang c FAQ 2 6 Words of Wisdom FFTW implements a method for saving plans to disk and restoring them In fact what FFTW does is more general than just saving and loading plans The mechanism is called wisdom Here we describe this feature at a high level See Chapter 3 FFTW Reference page 17 for a less casual but more complete discussion of how to use wisdom in FFTW Plans created with the FFTW MEASURE option produce near optimal FFT performance but it can take a long time to compute a plan because FFTW must actually measure the runtime of many possible plans and select the best one This is designed for the situations where so many transforms of the same size must be computed that the start up time is irrelevant For short initialization times but slightly slower transforms we have provided FFTW ESTIMATE The wisdom mechanism is a way to get the best of both worlds There are however certain caveats that the user must be aware of in using wisdom For this reason wisdom is an optional feature which is not enabled by default At
73. int flags fftw complex in int istride fftw complex out int ostride The function fftwnd create plan creates a plan which is a data structure containing all the information that fftwnd needs in order to compute a multi dimensional Fourier transform You can create as many plans as you need but only one plan for a given array size is required a plan can be reused many times The functions fftw2d create plan and fftw3d create plan are optional alternative interfaces to fftwnd create plan for two and three dimensions respectively fftwnd create plan returns a valid plan or NULL if for some reason the plan can t be created This can happen if memory runs out or if the arguments are invalid in some way e g if rank lt 0 The create_plan_specific variants take as additional arguments specific input output arrays and their strides For the last four arguments you should pass the arrays and strides that you will eventually be passing to fftwnd The resulting plans will be optimized for those arrays and strides although they may be used on other arrays as well Note the Chapter 3 FFTW Reference 23 contents of the in and out arrays are destroyed by the specific planner the initial contents are ignored so the arrays need not have been initialized See Section 3 2 2 Discussion on Specific Plans page 20 for a discussion on specific plans Arguments rank is the dimensionality of the arrays to be transformed It can be any no
74. isdom causes all accumulated wisdom to be discarded and its associated memory to be freed New wisdom can still be gathered subsequently however 3 7 Memory Allocator Reference include lt fftw h gt void fftw_malloc_hook size_t n void fftw_free_hook void p Whenever it has to allocate and release memory FFTW ordinarily calls malloc and free If malloc fails FFTW prints an error message and exits This behavior may be undesirable in some applications Also special memory handling functions may be necessary in certain environments Consequently FFTW provides means by which you can install your own memory allocator and take whatever error correcting action you find appropriate The 36 FFTW variables fftw malloc hook and fftw_free_hook are pointers to functions and they are normally NULL If you set those variables to point to other functions then FFTW will use your routines instead of malloc and free fftw_malloc_hook must point to a malloc like function and fftw free hook must point to a free like function 3 8 Thread safety Users writing multi threaded programs must concern themselves with the thread safety of the libraries they use that is whether it is safe to call routines in parallel from multiple threads FFTW can be used in such an environment but some care must be taken because certain parts of FFTW use private global variables to share data between calls In particular the plan creation functions share t
75. ithms that we do not yet completely understand The reader is referred to the cited papers for the appropriate references The rest of this manual is organized as follows We first discuss the sequential one processor implementation We start by describing the basic features of FFTW in Chap ter 2 Tutorial page 3 This discussion includes the storage scheme of multi dimensional arrays Section 2 5 Multi dimensional Array Format page 11 and FFTW s mechanisms for storing plans on disk Section 2 6 Words of Wisdom page 13 Next Chapter 3 FFTW Reference page 17 provides comprehensive documentation of all FFTW s features Parallel transforms are discussed in their own chapter Chapter 4 Parallel FFTW page 37 Fortran programmers can also use FF TW as described in Chapter 5 Calling FFTW from Fortran page 51 Chapter 6 Installation and Customization page 55 explains how to install FFTW in your computer system and how to adapt FFTW to your needs License and copyright information is given in Chapter 8 License and Copyright page 63 Finally we thank all the people who helped us in Chapter 7 Acknowledgments page 61 Chapter 2 Tutorial 3 2 Tutorial This chapter describes the basic usage of FFTW i e how to compute the Fourier trans form of a single array This chapter tells the truth but not the whole truth Specifically FFTW implements additional routines and flags providing extra functionality that are not documented
76. its simplest wisdom provides a way of saving plans to disk so that they can be reused in other program runs You create a plan with the flags FFTW MEASURE and FFTW USE WISDOM and then save the wisdom using fftw export wisdom plan fftw create plan FFTW MEASURE FFTW USE WISDOM fftw export wisdom The next time you run the program you can restore the wisdom with fftw import wisdom and then recreate the plan using the same flags as before This time however the same optimal plan will be created very quickly without measurements FFTW still needs some time to compute trigonometric tables however The basic outline is fftw import wisdom plan fftw create plan FFTW USE WISDOM Wisdom is more than mere rote memorization however FF TW s wisdom encompasses all of the knowledge and measurements that were used to create the plan for a given size Therefore existing wisdom is also applied to the creation of other plans of different sizes Whenever a plan is created with the FEI MEASURE and FFTW USE WISDOM flags wisdom is generated Thereafter plans for any transform with a similar factorization will be com puted more quickly so long as they use the FFTW USE WISDOM flag In fact for transforms with the same factors and of equal or lesser size no measurements at all need to be made and an optimal plan can be created with negligible delay For example suppose that you create a plan for N 2 9
77. itwise or An in place transform is one in which the output data overwrite the input data It thus requires half as much memory as and is often faster than its opposite an out of place transform For two and three dimensional transforms FFTWND provides alternative routines that accept the sizes of each dimension directly rather than indirectly through a rank and an array of sizes These are otherwise identical to fftwnd_create_plan and are sometimes more convenient fftwnd_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 Once the plan has been created you can use it any number of times for transforms of the same size When you do not need a plan anymore you can deallocate the plan by calling fftwnd_destroy_plan plan Given a plan you can compute the transform of an array of data by calling void fftwnd one fftwnd plan plan fftw complex sin fftw complex out Here in and out point to multi dimensional arrays in row major order of the size specified when the plan was created In the case of an in place transform the out parameter is ignored and the output data are stored in the in array The results are stored in order unnormalized with the zero frequency component in out 0 A forward followed by a backward transform or vice versa yields the original data multiplied by the size of the array i e the produ
78. k of the array is one In both cases the out parameter is ignored for in place transforms The format of the complex arrays deserves careful attention Suppose that the real data has dimensions n X nz X x n4 in row major order Then after a real to complex transform the output is an n X na x x na 2 1 array of f ftw complex values in row major order corresponding to slightly over half of the output of the corresponding complex transform Note that the division is rounded down The ordering of the data is otherwise exactly the same as in the complex case In principle the output could be exactly half the size of the complex transform output but in more than one dimension this requires too complicated a format to be practical Note that unlike the one dimensional RFFTW the real and imaginary parts of the DFT amplitudes are here stored together in the natural way Since the complex data is slightly larger than the real data some complications arise for in place transforms In this case the final dimension of the real data must be padded with extra values to accommodate the size of the complex data two extra if the last dimension is even and one if it is odd That is the last dimension of the real data must physically contain 2 n4 2 1 fftw real values exactly enough to hold the complex data This physical array size does not however change the logical array size only nq values are actually stored in the last dimensio
79. le type prefix Adds a d or s prefix to all installed libraries and header files to indicate the floating point precision See Section 6 3 Installing FFTW in both single and double precision page 57 enable type prefix lt prefix gt lets you add an arbitrary prefix By default no prefix is used e enable threads Enables compilation and installation of the FFTW threads library see Section 4 1 Multi threaded FFTW page 37 which provides a simple interface to parallel transforms for SMP systems By default the threads routines are not compiled e with openmp with sgimp In conjunction with enable threads causes the multi threaded FFTW library to use OpenMP or SGI MP compiler directives in order to induce parallelism rather than spawning its own threads directly Useful especially for programs already employing such directives in order to minimize conflicts between different parallelization mechanisms e enable mpi Enables compilation and installation of the FFTW MPI library see Section 4 2 MPI FFTW page 40 which provides parallel transforms for distributed memory systems with MPI By default the MPI routines are not compiled e disable fortran Disables inclusion of Fortran callable wrapper routines see Chap ter 5 Calling FFTW from Fortran page 51 in the standard FFTW libraries These wrapper routines increase the library size by only a negligible amount so they are included by default as long a
80. ly this means that as you step through adjacent memory locations the first dimension s index varies most slowly and the last dimension s index varies most quickly To be more explicit let us consider an array of rank d whose dimensions are n X nz x na X X na Now we specify a location in the array by a sequence of zero based indices one for each dimension 4 i5 43 i4 If the array is stored in row major order then this element is located at the position tg nalia 1 nq 1 Netz Note that each element of the array must be of type fftw_complex i e a real imagi nary pair of double precision numbers Note also that in fftwnd the expression above is multiplied by the stride to get the actual array index this is useful in situations where each element of the multi dimensional array is actually a data structure or another array and you just want to transform a single field In most cases however you use a stride of 1 2 5 2 Column major Format Readers from the Fortran world are used to arrays stored in column major order some times called Fortran order This is essentially the exact opposite of row major order in that here the first dimension s index varies most quickly If you have an array stored in column major order and wish to transform it using fftwnd it is quite easy to do When creating the plan simply pass the dimensions of the array to fftwnd create plan in reverse order For example if your arr
81. m odist_t is idist 2 3 5 5 Destroying a Multi dimensional Plan include lt rfftw h gt void rfftwnd_destroy_plan rfftwnd_plan plan The function rfftwnd_destroy_plan frees the plan plan and releases all the memory associated with it After destruction a plan is no longer valid 3 5 6 What RFFTWND Really Computes The conventions that we follow for the real multi dimensional transform are analogous to those for the complex multi dimensional transform In particular the forward transform has a negative sign in the exponent and neither the forward nor the backward transforms will perform any normalization Computing the backward transform of the forward transform will multiply the array by the product of its dimensions that is the logical dimensions of the real data The forward transform is real to complex and the backward transform is complex to real The exact mathematical definition of our real multi dimensional transform follows Real to complex forward transform Let X be a d dimensional real array whose elements are X j1 Jo ja where 0 js lt n for all s 1 2 d Let also w e V s for all s 1 2 d The real to complex transform computes a complex array Y whose structure is the same as that of X defined by ni1 1 n2 1 ngq 1 AT E le e Xh jz jaor ug tros ms j1 0 j2 0 Ja 0 34 FFTW The output array Y enjoys a multidimensional hermitian symmetry that is the identity
82. m performing a complex two dimensional MPI transform might look something like include lt fftw_mpi h gt int main int argc char argv x const int NX NY fftwnd mpi plan plan fftw complex data MPI Init amp argc amp argv plan fftw2d mpi create plan MPI COMM WORLD NX NY FFTW FORWARD FFTW ESTIMATE allocate and initialize data fftwnd mpi p 1 data NULL FFTW NORMAL ORDER 42 FFTW fftwnd mpi destroy plan plan MPI Finalize The calls to MPI_Init and MPI_Finalize are required in all MPI programs see the MPI home page for more information Note that all of your processes run the program in parallel as a group there is no explicit launching of threads processes in an MPI program As in the ordinary FF TW the first thing we do is to create a plan of type fftwnd_ mpi_plan using fftwnd_mpi_plan fftw2d_mpi_create_plan MPI_Comm comm int nx int ny fftw_direction dir int flags Except for the first argument the parameters are identical to those of fftw2d create plan There are also analogous fftwnd mpi create plan and fftw3d mpi create plan functions Transforms of any rank greater than one are supported The first argument is an MPI communicator which specifies the group of processes that are to be involved in the transform the standard constant MPI COMM WORLD indicates all available processes Next one has to allocate and initialize the data This is somewhat tricky
83. mation Transforms whose sizes are powers of 2 are especially fast e dir is the sign of the exponent in the formula that defines the Fourier transform It can be 1 or 1 The aliases FFTW_FORWARD and FFTW_BACKWARD are provided where FFTW_FORWARD stands for 1 e flags is a boolean OR of zero or more of the following FFTW_MEASURE this flag tells FF TW to find the optimal plan by actually computing several FFTs and measuring their execution time Depending on the installation this can take some time FFTW_ESTIMATE do not run any FFT and provide a reasonable plan for a RISC processor with many registers If neither FFTW_ESTIMATE nor FFTW_MEASURE is provided the default is FFTW_ESTIMATE FFTW OUT OF PLACE produce a plan assuming that the input and output arrays will be distinct this is the default FFTW IN PLACE produce a plan assuming that you want the output in the input array The algorithm used is not necessarily in place FF TW is able to compute true in place transforms only for small values of n If FFTW is not able to compute the transform in place it will allocate a temporary array unless you provide one yourself compute the transform out of place and copy the result back Warning This option changes the meaning of some parameters of fftw see Section 3 2 3 Computing the One dimensional Transform page 20 The in place option is mainly provided for people who want to write their own in place m
84. n 64 FFTW Chapter 9 Concept Index 9 Concept Index A Ee EE 2 B benchfft i e pb Saad dace Sg 1 benchmark 22er exu 1 40 48 58 blocki gussoReekbe ter ERU r EE Ed RR 43 C C multi dimensional aas 12 EE 59 61 di EE 2 37 e 58 code generat ss ge R dE 1 59 column major EELER 11 52 54 COMPILER gasetenn ERR e tton 2 51 55 56 58 Compiler flags Ee ed nde mue 55 56 complex multi dimensional transform 4 22 complex nuniber 2222e 9 ltem RR Rede 17 complex one dimensional transform 3 complex to real transform 6 26 complex transform 00 eee eee ee eee 3 EE 18 37 41 55 57 conyolutions este chases eod aa ema eo tog eds 10 cyclic convolution SANS nr kerine e 10 D Discrete Fourier Transform 21 25 28 33 distributed array format 43 46 48 distributed memory suuus 37 40 43 E Ecclesiastes 00 0 ccc cece ccc nese ee eee 14 EXECUTOR aeos A Ar AE acd rete aoe gia J EREM etu tete ct et Ria 1 EID MUNIDY be elen ee 29 flags EEN NN 4 5 19 23 27 30 47 52 floating point precision 18 52 55 56 57 Fortran callable wrappers 12 51 56 frequenGey 2b aires E REOR RA Dd 4 5 21 65 gettimeofday 2 22 sdb a iqq d bREEd ads 58 ke 61 H halfcomplex array 0 ee eee eee 6 17 hermitian array ositos eser EN 17 33 34 I in place transform 5 19 21 24 30 3
85. n and n is the last dimension passed to rfftwnd create plan For example consider the transform of a two dimensional real array of size nx by ny The output of the rfftwnd transform is a two dimensional complex array of size nx by ny 2 1 where the y dimension has been cut nearly in half because of redundancies in the output Because fftw complex is twice the size of fftw real the output array is slightly bigger than the input array Thus if we want to compute the transform in place we must pad the input array so that it is of size nx by 2 ny 2 1 If ny is even then there are two padding elements at the end of each row which need not be initialized as they are only used for output Figure 1 depicts the input and output arrays just described for both the out of place and in place transforms with the arrows indicating consecutive memory locations The RFFTWND transforms are unnormalized so a forward followed by a backward transform will result in the original data scaled by the number of real data elements that is the product of the logical dimensions of the real data Below we illustrate the use of RFFTWND by showing how you might use it to compute the cyclic convolution of two dimensional real arrays a and b using the identity that a Chapter 2 Tutorial Figure 1 Illustration of the data layout for real to complex transforms 10 FFTW convolution corresponds to a pointwise product of the Fourier transforms For v
86. n and Customization This chapter describes the installation and customization of FFTW the latest version of which may be downloaded from the FFTW home page As distributed FFTW makes very few assumptions about your system All you need is an ANSI C compiler gcc is fine although vendor provided compilers often produce faster code However installation of FFTW is somewhat simpler if you have a Unix or a GNU system such as Linux In this chapter we first describe the installation of FFTW on Unix and non Unix systems We then describe how you can customize FFTW to achieve better performance Specifically you can I enable gcc x86 specific hacks that improve performance on Pentia and PentiumPro s II adapt FFTW to use the high resolution clock of your machine if any III produce code codelets to support fast transforms of sizes that are not supported efficiently by the standard FFTW distribution 6 1 Installation on Unix FFTW comes with a configure program in the GNU style Installation can be as simple as configure make make install This will build the uniprocessor complex and real transform libraries along with the test programs We strongly recommend that you use GNU make if it is available on some systems it is called gmake The make install command installs the fftw and rfftw libraries in standard places and typically requires root privileges unless you specify a different install directory with the prefix flag to co
87. n negative integer n is a pointer to an array of rank integers giving the size of each dimension of the arrays to be transformed These sizes which must be positive integers correspond to the dimensions of row major arrays i e n 0 is the size of the dimension whose indices vary most slowly and so on See Section 2 5 Multi dimensional Array Format page 11 for more information on row major storage See Section 3 2 1 Plan Creation for One dimensional Transforms page 18 for more information regarding optimal array sizes nx and ny in fftw2d create plan are positive integers specifying the dimensions of the rank 2 array to be transformed i e they specify that the transform will operate on nx x ny arrays in row major order where nx is the number of rows and ny is the number of columns nx ny and nz in fftw3d create plan are positive integers specifying the dimensions of the rank 3 array to be transformed i e they specify that the transform will operate on nx x ny x nz arrays in row major order dir is the sign of the exponent in the formula that defines the Fourier transform It can be 1 or 1 The aliases FFTW FORWARD and FFTW_BACKWARD are provided where FFTW FORWARD stands for 1 flags is a boolean OR of zero or more of the following FFTW MEASURE this flag tells FFTW to find the optimal plan by actually computing several FFTs and measuring their execution time FFTW ESTIMATE do not run any FFT a
88. nd provide a reasonable plan for a RISC processor with many registers If neither FFTW ESTIMATE nor FEI MEASURE is provided the default is FFTW ESTIMATE FFTW OUT OF PLACE produce a plan assuming that the input and output arrays will be distinct this is the default FFTW IN PLACE produce a plan assuming that you want to perform the trans form in place Unlike the one dimensional transform this really performs the transform in place Note that if you want to perform in place transforms you must use a plan created with this option The default mode of operation is FFTW OUT OF PLACE FFTW USE WISDOM use any wisdom that is available to help in the creation of the plan See Section 2 6 Words of Wisdom page 13 This can greatly speed the creation of plans especially with the FFTW MEASURE option FFTW ESTIMATE 2 fftwnd actually may use some temporary storage hidden in the plan but this storage space is only the size of the largest dimension of the array rather than being as big as the entire array Unless you use fftwnd to perform one dimensional transforms in which case the temporary storage required for in place transforms is as big as the entire array 24 FFTW plans can also take advantage of wisdom to produce a more optimal plan based on past measurements than the estimation heuristic would normally generate When the FFTW MEASURE option is used new wisdom will also be generated i
89. ndences are assumed in the examples below The examples also assume that FFTW was compiled in double precision the default 5 2 FFTW Constants in Fortran When creating plans in FFTW a number of constants are used to specify options such as FFTW FORWARD or FFTW USE WISDOM The same constants must be used with the wrapper routines but of course the C header files where the constants are defined can t be incorporated directly into Fortran code Instead we have placed Fortran equivalents of the FFTW constant definitions in the file fortran fftw f77 i of the FFTW package If your Fortran compiler supports a prepro cessor you can use that to incorporate this file into your code whenever you need to call FFTW Otherwise you will have to paste the constant definitions in directly They are integer FFTW FORWARD FFTW BACKWARD parameter FFTW FORWARD 1 FFTW BACKWARD 1 integer FFTW REAL TO COMPLEX FFTW COMPLEX TO REAL parameter FFTW REAL TO COMPLEX 1 FFTW COMPLEX TO REAL 1 integer FFTW_ESTIMATE FFTW_MEASURE parameter FFTW_ESTIMATE 0 FFTW_MEASURE 1 integer FFTW_OUT_OF_PLACE FFTW_IN_PLACE FFTW_USE_WISDOM parameter FFTW_OUT_OF_PLACE 0 parameter FFTW_IN_PLACE 8 FFTW_USE_WISDOM 16 integer FFTW_THREADSAFE parameter FFTW_THREADSAFE 128 In C you combine different flags like FFTW_USE_WISDOM and FFTW_MEASURE using the operator in Fortran you should just use 5 3 Fortran Examples In C you might
90. nding on your installation and on the precision of the timer in your machine FFTW_ESTIMATE on the contrary does not run any computation and just builds a reasonable plan which may be sub optimal In other words if your program performs many transforms of the same size and initialization time is not important use FFTW_MEASURE otherwise use the estimate A compromise between these two extremes exists See Section 2 6 Words of Wisdom page 13 Once the plan has been created you can use it as many times as you like for transforms on arrays of the same size When you are done with the plan you deallocate it by calling fftw_destroy_plan plan The transform itself is computed by passing the plan along with the input and output arrays to fftw one void fftw one fftw plan plan fftw complex sin fftw complex out Note that the transform is out of place in and out must point to distinct arrays It operates on data of type fftw complex a data structure with real in i re and imaginary in i im floating point components The in and out arrays should have the length specified when the plan was created An alternative function fftw allows you to efficiently perform multiple and or strided transforms see Chapter 3 FFTW Reference page 17 The DFT results are stored in order in the array out with the zero frequency DC component in out 0 The array in is not modified Users should note that FFTW computes an unnormalized DFT the si
91. nfigure You can also type make check to put the FFTW test programs through their paces If you have problems during configuration or compilation you may want to run make distclean before trying again this ensures that you don t have any stale files left over from previous compilation attempts The configure script knows good CFLAGS C compiler flags for a few systems If your system is not known the configure script will print out a warning In this case you can compile FFTW with the command make CFLAGS lt write your CFLAGS here gt If you do find an optimal set of CFLAGS for your system please let us know what they are along with the output of config guess so that we can include them in future releases The configure program supports all the standard flags defined by the GNU Coding Standards see the INSTALL file in FF TW or the GNU web page Note especially help to list all flags and enable shared to create shared rather than static libraries configure also accepts a few FF TW specific flags particularly e enable float Produces a single precision version of FFTW float instead of the default double precision double See Section 6 3 Installing FFTW in both single and double precision page 57 1 Each version of cc seems to have its own magic incantation to get the fastest code most of the time you d think that people would have agreed upon some convention e g Omax by now 56 FFTW e enab
92. ng several FFTs and measuring their execution time Depending on the in stallation this can take some time FFTW_ESTIMATE do not run any FFT and provide a reasonable plan for a RISC processor with many registers If neither FFTW_ESTIMATE nor FFTW_MEASURE is provided the default is FFTW_ESTIMATE FFTW_OUT_OF_PLACE produce a plan assuming that the input and output arrays will be distinct this is the default FFTW_IN_PLACE produce a plan assuming that you want the output in the input array The algorithm used is not necessarily in place RFF TW is able to compute true in place transforms only for small values of n If RFFTW is not able to compute the transform in place it will allocate a temporary array unless you provide one yourself compute the transform out of place and copy the result back Warning This option changes the meaning of some parameters of rfftw see Section 3 4 2 Computing the Real One dimensional Transform page 27 The default mode of operation is FFTW_OUT_OF_PLACE FFTW_USE_WISDOM use any wisdom that is available to help in the creation of the plan See Section 2 6 Words of Wisdom page 13 This can greatly speed the creation of plans especially with the FFTW_MEASURE option FFTW_ESTIMATE plans can also take advantage of wisdom to produce a more optimal plan based on past measurements than the estimation heuristic would normally generate When the FFTW_MEASURE option is used new wisdom will also
93. ns of egcs since 1 0 3 These optimizations affect only the performance and not the correctness of FF TW i e it is always safe to try them out 58 FFTW These hacks provide a workaround to the incorrect alignment of local double variables in gcc The compiler aligns these variables to multiples of 4 bytes but execution is much faster on Pentium and PentiumPro if doubles are aligned to a multiple of 8 bytes By carefully counting the number of variables allocated by the compiler in performance critical regions of the code we have been able to introduce dummy allocations using alloca that align the stack properly The hack depends crucially on the compiler flags that are used For example it won t work without fomit frame pointer In principle these hacks are no longer required under gcc versions 2 95 and later which automatically align the stack correctly see mpreferred stack boundary in the gcc man ual However we have encountered a bug in the stack alignment of versions 2 95 012 that causes FFTW s stack to be misaligned under some circumstances The configure script automatically detects this bug and disables gcc s stack alignment in favor of our own hacks when enable i386 hacks is used The fftw test program outputs speed measurements that you can use to see if these hacks are beneficial The configure option enable pentium timer enables the use of the Pentium and PentiumPro cycle counter for timing purposes In orde
94. on 2 5 Multi dimensional Array Format page 11 we discuss the various alternatives for storing multi dimensional arrays in memory Section 2 6 Words of Wisdom page 13 shows how you can save FFTW s plans for future use 2 1 Complex One dimensional Transforms Tutorial The basic usage of FFTW is simple A typical call to FFTW looks like include lt fftw h gt fftw complex in N out N fftw plan p ps fftw_create_plan N FFTW_FORWARD FFTW_ESTIMATE SEN in out eve daderoycniantes d The first thing we do is to create a plan which is an object that contains all the data that FFTW needs to compute the FFT using the following function fftw_plan fftw_create_plan int n fftw_direction dir int flags The first argument n is the size of the transform you are trying to compute The size n can be any positive integer but sizes that are products of small factors are transformed most 4 FFTW efficiently The second argument dir can be either FFTW_FORWARD or FFTW_BACKWARD and indicates the direction of the transform you are interested in Alternatively you can use the sign of the exponent in the transform 1 or 1 which corresponds to FFTW_FORWARD or FFTW_BACKWARD respectively The flags argument is either FFTW_MEASURE or FFTW_ ESTIMATE FFTW_MEASURE means that FFTW actually runs and measures the execution time of several FFTs in order to find the best way to compute the transform of size n This may take some time depe
95. on on Cilk see the Cilk home page The FFTW Cilk code can be found in the cilk directory with parallelized one and multi dimensional transforms of complex data The Cilk FFTW routines are documented in cilk README 4 1 Multi threaded FFTW In this section we document the parallel FFTW routines for shared memory threads on SMP hardware These routines which support parallel one and multi dimensional trans forms of both real and complex data are the easiest way to take advantage of multiple processors with FFTW They work just like the corresponding uniprocessor transform rou tines except that they take the number of parallel threads to use as an extra parameter Any program that uses the uniprocessor FF T W can be trivially modified to use the multi threaded FFTW 4 1 1 Installation and Supported Hardware Software A of the FFTW threads code is located in the threads subdirectory of the FFTW package On Unix systems the FFTW threads libraries and header files can be auto matically configured compiled and installed along with the uniprocessor FFTW libraries simply by including enable threads in the flags to the configure script see Section 6 1 Installation on Unix page 55 Note also that the threads routines when enabled are automatically tested by the make check self tests The threads routines require your operating system to have some sort of shared memory threads support Specifically the FFTW threads package works
96. on on Specific Dans 20 3 2 3 Computing the One dimensional Transform 20 3 2 4 Destroying a One dimensional Plan 21 3 25 What FFTW Really Computes 21 3 3 Multi dimensional Transforms Reference 22 3 83 1 Plan Creation for Multi dimensional Transforms 22 3 3 2 Computing the Multi dimensional Transform 24 3 3 3 Destroying a Multi dimensional Plan 25 3 3 4 What FFTWND Really Computes 25 3 4 Real One dimensional Transforms Reference 26 3 4 1 Plan Creation for Real One dimensional Transforms EE 26 3 4 2 Computing the Real One dimensional Transform SE E Ee 2T 3 4 3 Destroying a Real One dimensional Plan 28 3 4 4 What RFFTW Really Computes 28 3 5 Real Multi dimensional Transforms Reference 29 3 5 1 Plan Creation for Real Multi dimensional Transforms m 29 3 5 3 Array Dimensions for Real Multi dimensional Transforms dest Sg Maoh Ree Sender a 32 3 5 4 Strides in In place RFFTWND 32 3 5 5 Destroying a Multi dimensional Plan 33 3 5 6 What RFFTWND Really Computes 33 3 6 Wisdom Reference 00 cece eee eect ees 34 3 6 1 Exporting Wisdom 34 3 6 2 Importing Wisdom 35 3 6 3 Forgetting Wisdom 35 3 7 Memory Allocator Reference 0000 eee eee eens 35 3 8 Thread safety 0 cece cence nee eee eee eee 36 4 Parallel F
97. onflicting wisdom in which case the older wisdom is discarded The format of the exported wisdom is nerd readable LISP like ASCII text we will not document it here except to note that it is insensitive to white space interested users can contact us for more details See Chapter 3 FFTW Reference page 17 for more information and for a description of how you can implement wisdom import export for other media besides files and strings The following is a brief example in which the wisdom is read from a file a plan is created possibly generating more wisdom and then the wisdom is exported to a string and printed to stdout 1 fftw plan plan char wisdom string FILE input file open file to read wisdom from input file fopen sample wisdom r if FFTW FAILURE fftw import wisdom from file input file printf Error reading wisdom Nn fclose input_file be sure to close the file create a plan for N 64 possibly creating and or using wisdom plan fftw_create_plan 64 FFTW_FORWARD FFTW MEASURE FFTW_USE_WISDOM do some computations with the plan always destroy plans when you are done fftw_destroy_plan plan write the wisdom to a string wisdom string fftw_export_wisdom_to_string if wisdom string NULL printf Accumulated wisdom s n wisdom_string 16 FFTW Just for fun destroy and restore the wisdom fftw_forget_wisdom
98. r real transforms and also in the example above this is most easily handled by typecasting to a complex pointer when handling the complex data e As with the complex transforms there are also rfftwnd_create_plan and rfftw2d_ create_plan functions and any rank greater than one is supported Programs using the MPI FF TW real transforms should link with 1rfftw mpi lfftw mpi lrfftw lfftw 1m on Unix 4 2 5 Usage of MPI FFTW for Complex One dimensional Transforms The MPI FFTW also includes routines for parallel one dimensional transforms of complex data only Although the speedup is generally worse than it is for the multi dimensional routines these distributed memory one dimensional transforms are especially useful for performing one dimensional transforms that don t fit into the memory of a single machine The usage of these routines is straightforward and is similar to that of the multi dimensional MPI transform functions You first include the header fftw mpi h and then create a plan by calling fftw mpi plan fftw mpi create plan MPI Comm comm int n fftw direction dir int flags The last three arguments are the same as for fftw create plan except that all MPI transforms are automatically FFTW IN PLACE The first argument specifies the group of processes you are using and is usually MPI COMM WORLD all processes A plan can be used for many transforms of the same size and is destroyed when you are done with it
99. r to get correct results you must define FFTW CYCLES DER SEC in fftw config h to be the clock speed of your processor the resulting FFTW library will be nonportable The use of this option is deprecated On serious operating systems such as Linux FFTW uses gettimeofday O which has enough resolution and is portable Note that Win32 has its own high resolution timing routines as well FFTW contains unsupported code to use these routines 6 5 Customizing the timer FFTW needs a reasonably precise clock in order to find the optimal way to compute a transform On Unix systems configure looks for gettimeofday and other system specific timers If it does not find any high resolution clock it defaults to using the clock function which is very portable but forces FF T W to run for a long time in order to get reliable measurements If your machine supports a high resolution clock not recognized by FF TW it is therefore advisable to use it You must edit fftw fftw int h There are a few macros you must redefine The code is documented and should be self explanatory By the way fftw int stands for fftw internal but for some inexplicable reason people are still using primitive systems with 8 3 filenames Even if you don t install high resolution timing code we still recommend that you look at the FFTW TIME MIN constant in fftw fftw int h This constant holds the minimum time interval in seconds required to get accurate timing measurements
100. r web page also has links to FFT related information online FFTW is usually faster and sometimes much faster than all other freely available Fourier transform programs found on the Net For transforms whose size is a power of two it compares favorably with the FFT codes in Sun s Performance Library and IBM s ESSL library which are targeted at specific machines Moreover FF TW s performance is portable Indeed FFTW is unique in that it automatically adapts itself to your machine your cache the size of your memory the number of registers and all the other factors that normally make it impossible to optimize a program for more than one machine An extensive comparison of FFTW s performance with that of other Fourier transform codes has been made The results are available on the Web at the benchFF T home page In order to use FF TW effectively you need to understand one basic concept of FFTW s internal structure FFTW does not used a fixed algorithm for computing the transform but it can adapt the DFT algorithm to details of the underlying hardware in order to achieve best performance Hence the computation of the transform is split into two phases First FFTW s planner is called which learns the fastest way to compute the transform on your machine The planner produces a data structure called a plan that contains this information Subsequently the plan is passed to FFTW s executor along with an array of input data The executor
101. red by the fftw export wisdom functions above See Section 2 6 Words of Wisdom page 13 The imported wisdom supplements rather than replaces any wisdom already accumulated by the running program except when there is conflicting wisdom in which case the existing wisdom is replaced fftw import wisdom imports wisdom from any input medium as specified by the call back function get input get input is a getc like function that returns the next character in the input its parameter is the data pointer passed to fftw import wisdom If the end of the input data is reached which should never happen for valid data it may return either NULL ASCII 0 or EOF as defined in stdio h For convenience the following two wrapper routines are provided fftw import wisdom from file reads wisdom from the current position in input file which should be open with read permission Upon exit the file remains open and is positioned at the end of the wisdom data fftw import wisdom from string reads wisdom from the NULL terminated string input string The return value of these routines is FFTW SUCCESS if the wisdom was read successfully and FFTW FAILURE otherwise Note that in all of these functions any data in the input stream past the end of the wisdom data is simply ignored it is not even read if the wisdom data is well formed 3 6 3 Forgetting Wisdom include lt fftw h gt void fftw_forget_wisdom void Calling fftw forget w
102. rigonometric tables and accumulated wisdom Users should note that these comments only apply to programs using shared memory threads Parallelism using MPI or forked processes involves a separate address space and global variables for each process and is not susceptible to problems of this sort The central restriction of FFTW is that it is not safe to create multiple plans in parallel You must either create all of your plans from a single thread or instead use a semaphore mutex or other mechanism to ensure that different threads don t attempt to create plans at the same time The same restriction also holds for destruction of plans and import ing forgetting wisdom Once created a plan may safely be used in any thread The actual transform routines in FFTW fftw one etcetera are re entrant and thread safe so it is fine to call them simultaneously from multiple threads Another question arises however is it safe to use the same plan for multiple transforms in parallel It would be unsafe if for example the plan were modified in some way by the transform We address this question by defining an additional planner flag FFTW THREADSAFE When included in the flags for any of the plan creation routines FFTW THREADSAFE guarantees that the resulting plan will be read only and safe to use in parallel by multiple threads Chapter 4 Parallel FFTW 37 4 Parallel FFTW In this chapter we discuss the use of FFTW in a parallel environment
103. s the configure script finds a Fortran compiler on your system e with gcc Enables the use of gcc By default FFTW uses the vendor supplied cc compiler if present Unfortunately gcc produces slower code than cc on many systems e enable i386 hacks See Section 6 4 gcc and Pentium hacks page 57 below e enable pentium timer See Section 6 4 gcc and Pentium hacks page 57 below To force configure to use a particular C compiler instead of the default usually cc set the environment variable CC to the name of the desired compiler before running configure you may also need to set the flags via the variable CFLAGS 6 2 Installation on non Unix Systems It is quite straightforward to install FF TW even on non Unix systems lacking the niceties of the configure script The FFTW Home Page may include some FFTW packages pre configured for particular systems compilers and also contains installation notes sent in by users All you really need to do though is to compile all of the c files in the appropriate directories of the FFTW package You needn t worry about the many extraneous files lying around For the complex transforms compile all of the c files in the fftw directory and link them into a library Similarly for the real transforms compile all of the c files in the rfftw directory into a library Note that these sources include various files in the fftw and rfftw directories so you may need to set up the include paths for
104. ssible case Chapter 4 Parallel FFTW 49 4 fftw mpi test s 128x128x128 will benchmark a 128x128x128 transform on four pro cessors reporting timings and parallel speedups for all variants of f ftwnd mpi transposed with workspace etcetera Note also that there is the rfftw mpi test program for the real transforms 50 FFTW Chapter 5 Calling FFTW from Fortran 51 5 Calling FFTW from Fortran The standard FFTW libraries include special wrapper functions that allow Fortran pro grams to call FFTW subroutines This chapter describes how those functions may be employed to use FFTW from Fortran We assume here that the reader is already familiar with the usage of FFTW in C as described elsewhere in this manual In general it is not possible to call C functions directly from Fortran due to Fortran s inability to pass arguments by value and also because Fortran compilers typically expect identifiers to be mangled somehow for linking However if C functions are written in a special way they are callable from Fortran and we have employed this technique to create Fortran callable wrapper functions around the main FFTW routines These wrapper functions are included in the FFTW libraries by default unless a Fortran compiler isn t found on your system or disable fortran is included in the configure flags As a result calling FFTW from Fortran requires little more than appending Er to the function names and then linking normall
105. systems When the enable type prefix option of configure is used the FFTW libraries and header files are installed with a prefix of d or s depending upon whether you compiled in double or single precision Then instead of linking your program with 1rfftw lfftw for example you would link with 1drfftw ldfftw to use the double precision version or with lsrfftw lsfftw to use the single precision version Also you would include lt drfftw h gt or lt srfftw h gt instead of lt rfftw h gt and so on The names of FFTW functions data types and constants remain unchanged You still call for instance fftw one and not dfftw one Only the names of header files and libraries are modified One consequence of this is that you cannot use both the single and double precision FFTW libraries in the same program simultaneously as the function names would conflict So to install both the single and double precision libraries on the same machine you would do configure enable type prefix other options make make install make clean configure enable float enable type prefix other options make make install 6 4 gcc and Pentium hacks The configure option enable i386 hacks enables specific optimizations for the Pen tium and later x86 CPUs under gcc which can significantly improve performance of double precision transforms Specifically we have tested these hacks on Linux with gcc 2 789 and versio
106. ter 8 License and Copyright 63 8 License and Copyright FFTW is copyright 1997 1999 Massachusetts Institute of Technology FFTW is free software you can redistribute it and or modify it under the terms of the GNU General Public License as published by the Free Software Foundation either version 2 of the License or at your option any later version This program is distributed in the hope that it will be useful but WITHOUT ANY WARRANTY without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE See the GNU General Public License for more details You should have received a copy of the GNU General Public License along with this program if not write to the Free Software Foundation Inc 59 Temple Place Suite 330 Boston MA 02111 1307 USA You can also find the GPL on the GNU web site In addition we kindly ask you to acknowledge FFTW and its authors in any program or publication in which you use FFTW You are not required to do so it is up to your common sense to decide whether you want to comply with this request or not Non free versions of FFTW are available under terms different than the General Public License e g they do not require you to accompany any object code using FFTW with the corresponding source code For these alternate terms you must purchase a license from MIT s Technology Licensing Office Users interested in such a license should contact us fftwOfftw org for more informatio
107. ters as ostride_t and odist_t Note that these are Chapter 3 FFTW Reference 33 computed internally by rfftwnd the actual ostride and odist parameters are ignored for in place transforms First there is the case where you are transforming a number of contiguous arrays located one after another in memory In this situation istride is 1 and idist is the product of the physical dimensions of the array ostride_t and odist_t are then chosen so that the output arrays are contiguous and lie on top of the input arrays ostride_t is therefore 1 For a real to complex transform odist_t is idist 2 for a complex to real transform odist_t is idist 2 The second case is when you have an array in which each element has nc components e g a structure with nc numeric fields and you want to transform all of the components at once Here istride is nc and idist is 1 For this case it is natural to want the output to also have nc consecutive components now of the output data type this is exactly what rfftwnd does Specifically it uses an ostride_t equal to istride and an odist_t of 1 Astute readers will realize that some extra buffer space is required in order to perform such a transform this is handled automatically by rfftwnd The general rule is as follows ostride_t equals istride If idist is 1 and idist is less than istride then odist_t is 1 Otherwise for a real to complex transform odist_t is idist 2 and for a complex to real transfor
108. twnd_threads_one_complex_to_real 39 rfftwnd_threads_one_real_to_complex 38 39 rfftwnd_threads_real_to_complex 38 rfftwnd complex to real 31 rfftwnd create plan 7 29 47 rfftwnd destroy plan sss 33 FFTW rfftwnd f77 destroy plan 54 rfftwnd f77 one real to complex 54 rfftwnd one complex to real 8 31 rfftwnd_one_real_to_complex 8 31 rfftwnd page oo sa Er RR ER 7 29 rfftwnd_real_to_complex 31 Table of Contents l Introductions ser r4 R4 PERCY RERUM Ey CER Y 1 2 SOM BEE 3 2 1 Complex One dimensional Transforms Tutorial 3 2 2 Complex Multi dimensional Transforms Tutorial 4 2 3 Real One dimensional Transforms Tutorial 6 2 4 Real Multi dimensional Transforms Tutorial 7 2 5 Multi dimensional Array boat 11 2 5 4 Row major Format 0 000002 eee 11 2 5 2 Column major boamat elles else 11 2 5 9 Static Arrays RE 12 2 5 4 Dynamic Arrays mm 12 2 5 5 Dynamic Arrays in C The Wrong Way 12 2 6 Words of WisdOm bae er epp beo Oops d 13 2 6 1 Caveats in Using Wisdom 14 2 6 2 Importing and Exporting Wisdom 14 3 FR TW Reference oc e NEE ces he 17 Ol BR TEE 17 3 2 One dimensional Transforms Reference aaaaaaaaaa 18 3 2 1 Plan Creation for One dimensional Transforms 18 3 2 2 Discussi
109. ulti dimensional Fourier transform using FFTW as a base For example consider a three dimensional n n n transform An out of place algorithm will need another array which may be huge However FFTW can compute the in place transform along each dimension using only a temporary array of size n Moreover if FFTW happens to be able to compute the transform truly in place no temporary array and no copying are needed As distributed FFTW knows how to compute in place transforms of size 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 32 and 64 The default mode of operation is FFTW OUT OF PLACE FFTW USE WISDOM use any wisdom that is available to help in the creation of the plan See Section 2 6 Words of Wisdom page 13 This can greatly speed the creation of plans especially with the FFTW MEASURE option FFTW ESTIMATE plans can also take advantage of wisdom to produce a more optimal plan based on past measurements than the estimation heuristic would normally generate When the FFTW MEASURE option is used new wisdom will also be generated if the current transform size is not completely understood by existing wisdom e in out istride ostride only for fftw create plan specific see correspond ing arguments in the description of fftw See Section 3 2 3 Computing the One dimensional Transform page 20 In particular the out and ostride parameters have the same special meaning for FFTW IN PLACE transforms as they
110. use the standard in order output ordering the k th output corresponds to the frequency k n or k T where T is your total sampling period For those who like to think in terms of positive and negative frequencies this means that the positive frequencies are stored in the first half of the output and the negative frequencies are stored in backwards order in the second half of the output The frequency k n is the same as the frequency n k n 22 FFTW 3 3 Multi dimensional Transforms Reference The multi dimensional complex routines are generally prefixed with fftwnd_ Programs using FFTWND should be linked with 1fftw lm on Unix systems or with the FF TW and standard math libraries in general 3 3 1 Plan Creation for Multi dimensional Transforms include lt fftw h gt fftwnd plan fftwnd_create_plan int rank const int n fftw direction dir int flags fftwnd 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 fftwnd plan fftwnd create plan specific int rank const int n fftw direction dir int flags fftw complex in int istride fftw complex out int ostride fftwnd plan fftw2d create plan specific int nx int ny fftw direction dir int flags fftw complex sin int istride fftw complex out int ostride fftwnd plan fftw3d create plan specific int nx int ny int nz fftw direction dir
111. used by the one dimensional real transforms taking advantage of the hermitian symmetry that arises in those cases By including lt fftw h gt or lt rfftw h gt you will have access to the following definitions typedef double fftw real typedef struct fftw_real re im fftw complex define c_re c c re define c_im c c im All FFTW operations are performed on the fftw real and fftw complex data types For fftw complex numbers the two macros c re and c im retrieve respectively the real and imaginary parts of the number A real array is an array of real numbers A complex array is an array of complex numbers A one dimensional array X of n complex numbers is hermitian if the following property holds for all 0 i lt n we have X X where x denotes the complex conjugate of x Hermitian arrays are relevant to FF TW because the Fourier transform of a real array is hermitian Because of its symmetry a hermitian array can be stored in half the space of a complex array of the same size FF TW s one dimensional real transforms store hermitian arrays as halfcomplex arrays A halfcomplex array of size n is a one dimensional array of n fftw_ real numbers A hermitian array X in stored into a halfcomplex array Y as follows For all integers i such that 0 i n 2 we have Y Re X For all integers i such that 0 i n 2 we have Y Im X We now illustrate halfcomplex storage for n 4 and n 5 since t
112. ute many transforms instead of simply calling rfftw many times e in istride and idist describe the input array s There are two cases If the plan defines a FFTW_REAL_TO_COMPLEX transform in is a real array Otherwise for FFTW_COMPLEX_TO_REAL transforms in is a halfcomplex array whose contents will be destroyed e out ostride and odist describe the output array s and have the same meaning as the corresponding parameters for the input array In place transforms If the plan specifies an in place transform ostride and odist are always ignored If out is NULL out is ignored too Otherwise out is interpreted as a pointer to an array of n complex numbers that FFTW will use as temporary space to perform the in place computation out is used as scratch space and its contents destroyed In this case out must be an ordinary array whose elements are contiguous in memory no striding The function rfftw one transforms a single contiguous input array to a contiguous output array By definition the call rfftw one plan in out is equivalent to rfftw plan 1 in 1 O out 1 O 3 4 3 Destroying a Real One dimensional Plan include lt rfftw h gt void rfftw_destroy_plan rfftw_plan plan The function rfftw_destroy_plan frees the plan plan and releases all the memory associated with it After destruction a plan is no longer valid 3 4 4 What RFFTW Really Computes In this section we define precisely what RFFTW computes
113. y with the FFTW libraries There are a few wrinkles however as we shall discuss below 5 1 Wrapper Routines All of the uniprocessor and multi threaded transform routines have Fortran callable wrappers except for the wisdom import export functions since it is not possible to ex change string and file arguments portably with Fortran and the specific planner routines see Section 3 2 2 Discussion on Specific Plans page 20 The name of the wrapper routine is the same as that of the corresponding C routine but with fftw fftwnd rfftw rfftwnd replaced by fftw f77 fftwnd f77T rfftw f 7 rfftwnd f77 For example in Fortran in stead of calling fftw one you would call fftw f77 one For the most part all of the arguments to the functions are the same with the following exceptions e plan variables what would be of type fftw plan rfftwnd plan etcetera in C must be declared as a type that is the same size as a pointer address on your machine Fortran has no generic pointer type The Fortran integer type is usually the same size as a pointer but you need to be wary especially on 64 bit machines You could also use integer 4 on a 32 bit machine and integer 8 on a 64 bit machine Ugh g77 has a special type integer kind 7 that is defined to be the same size as a pointer e Any function that returns a value e g fftw create plan is converted into a subrou tine The return value is converted into an additional first parameter of
114. you could for static arrays Instead you have to explicitly compute the offset into the array using the formula given earlier for row major arrays For example to reference the i j k th element of the array allocated above you would use the expression an array k 27 j 12 i This pain can be alleviated somewhat by defining appropriate macros or in C cre ating a class and overloading the operator 2 5 5 Dynamic Arrays in C The Wrong Way A different method for allocating multi dimensional arrays in C is often suggested that is incompatible with fftwnd using it will cause FFTW to die a painful death We discuss the technique here however because it is so commonly known and used This method is to create arrays of pointers of arrays of pointers of etcetera For example the analogue in this method to the example above is int i j fftw complex a bad array another way to make a 5x12x27 array a bad array fftw complex malloc 5 sizeof fftw complex for i 0 i lt 5 i a bad array i fftw complex malloc 12 sizeof fftw complex for j 0 j lt 12 j a bad array i jl fftw complex malloc 27 sizeof fftw complex Chapter 2 Tutorial 13 As you can see this sort of array is inconvenient to allocate and deallocate On the other hand it has the advantage that the i j k th element can be referenced simply by a_bad_array i j k If you l

FFTW User's Manual

Contents

Download Pdf Manuals

Related Search

Related Contents