YAKL FFTs

YAKL provides FFT routines to perform batched 1-D real-to-complex forward FFTs and complex-to-real inverse FFTs for multi-dimensional styleC, memDevice Array objects with float or double underlying types. Complex numbers are stored in interleaved format using float, meaning they are real,imag,real,imag format. All FFTs are in-place and overwrite the given data.

Your Array object needs to have room for the transformed data in each dimension you plan to transform. Assuming N is the number of elements you plan to transform in a given dimension, since this is a real-to-complex forward transform, you'll need N+2 elements available when N is even and N+1 elements available when N is odd. For this reason, YAKL cannot tell from the dimensions size alone how many elements you plan to transform for that dimension. Therefore, you must explicitly give the number of real-valued data elements you plan to transform in the forward and inverse transforms.

Because the numbers are interleaved, the result you get from the forward transform still has either float or double type, depending upon what you specify.

Example workflow

// nx, ny, and nz are the number of elements in each dimension
int nz = 100;
int ny = 50;
int nx = 41;

// We plan to transform the "y" and "x" dimensions of this array
// For even numbers of elements, you need two additional elements for the forward transform
// For odd numbers of elements, you need one additional element for the forward transform
int ny_trans = ny % 2 == 0 ? ny+2 : ny+1;
int nx_trans = nx % 2 == 0 ? nx+2 : nx+1;
Array<float,3,memDevice,styleC> data("data",nz,ny_trans,nx_trans);

RealFFT1D<float> fft_y;  // You need a separate FFT object for *each dimension you plan to transform*
RealFFT1D<float> fft_x;  // Even if the dimensions are the same size, the batching strategy is different
// init(Array const &arr, int tr_dim, int tr_size);
// This allocates and initializes temporary storage and data you need like twiddle and Chirp factors
// This routine is optional. If you don't call init(), then you must pass the
//   optional tr_dim and tr_size parameters to your forward transform, and allocation
//   and initialization will be performed on the first transform.
// The first parameter is the Array object you plan to transform. It need not be initialized with any
//   meaningful data. It's just there to tell you the sizes of the *other dimensions* to create
//   a batching strategy.
// The second parameter is the index (zero-based!!!) of the dimension you plan to transform. Again,
//   this is just to create the batching strategy.
// The third parameter is the number of real-valued elements you want to transform. Again, it's not
//   possible to know a priori what this number is because both even and odd sized transforms end up
//   with an even number of complex-valued elements in the forward transform.
fft_x.init(data , 2 , nx );
fft_y.init(data , 1 , ny );


// Initialize data in the array with some meaningful data
//   like momentum divergence for an incompressible solver using cyclic reduction
// Keep in mind you do not need to initialize any data in the array before calling init().
// init() only uses the size of the other dimensions to create a batching strategy.


// Perform the in-place forward transform, overwriting the array's original data
// forward_real( Array &arr , int tr_dim = -1 , int tr_size = -1 )
//  Allocation and initialization of temporary data is only performed if you did not call init()
//    first or if tr_dim or tr_size have changed since calling init() or your last forward_real() or
//    inverse_real() calls. Therefore, as long as tr_dim and tr_size do not change, there is no harm
//    in specifying them anyway if you want.
// If you called init(), and tr_dim and tr_size haven't changed (they typically don't need to change),
//   then you do not need to specify these optional parameters.
// This will transform N real values into either (N+1)/2 or (N+2)/2 complex numbers in interleaved
//   format if N is odd or even, respectively.
// The one or two extra elements of space needed for the forward transform are not used by this routine.
fft_x.forward_real( data );  // Batched over y and z dimensions
fft_y.forward_real( data );  // Batched over x and z dimensions

// Do fancy stuff in Fourier space

// Perform the in-place inverse transform, overwriting the array's original data
// inverse_real( Arra &arr , int tr_dim = -1 , int tr_size = -1 )
// Again, you only need to specify tr_dim or tr_size if they've changed since calling forward_real(),
//   and I cannot imagine any scenario where this would happen. Allocation and initialization of
//   temporary data is only performed if tr_dim or tr_size differ from their previous values from
//   forward_real(), so there is no harm in specifying them anyway if you want.
// This will transform the appropriate number of interleaved complex numbers into N real values.
fft_y.inverse_real( data );  // Batched over x and z dimensions
fft_x.inverse_real( data );  // Batched over y and z dimensions

// IMPORTANT: After the inverse transform, the extra one or two elements used for holding
//   complex data for the forward transform is undefined and must not be used.

// Do fancy stuff in physical space