The cuFFT and cuFFTW libraries are available as shared libraries. They consist of compiled programs ready for users to incorporate into applications with the compiler and linker. cuFFT can be downloaded from http://developer.nvidia.com/cufft. By selecting Download CUDA Production Release users are all able to install the package containing the CUDA Toolkit, SDK code samples and development drivers. The CUDA Toolkit contains cuFFT and the samples include simplecuFFT.

The Linux release for simplecuFFT assumes that the root install directory is /usr/local/cuda and that the locations of the products are contained there as follows. Modify the Makefile as appropriate for your system.

The most common case is for developers to modify an existing CUDA routine (for example, filename.cu) to call cuFFT routines. In this case the include file cufft.h or cufftXt.h should be inserted into filename.cu file and the library included in the link line. A single compile and link line might appear as

• /usr/local/cuda/bin/nvcc [options] filename.cu … -I/usr/local/cuda/inc -L/usr/local/cuda/lib -lcufft

Of course there will typically be many compile lines and the compiler g++ may be used for linking so long as the library path is set correctly.

Users of the FFTW interface (see FFTW Interface to cuFFT) should include cufftw.h and link with both cuFFT and cuFFTW libraries.

For the best performance input data should reside in device memory. Therefore programs in the cuFFT library assume that the data is in GPU memory. For example, if one of the execution functions is called with data in host memory, the program will return CUFFT_EXEC_FAILED. Programs in the cuFFTW library assume that the input data is in host memory since this library is a porting tool for users of FFTW. If the data resides in GPU memory, the program will abort.

The first step in using the cuFFT Library is to create a plan using one of the following:

• cufftPlan1D() / cufftPlan2D() / cufftPlan3D() - Create a simple plan for a 1D/2D/3D transform respectively.
• cufftPlanMany() - Creates a plan supporting batched input and strided data layouts.

Among the plan creation functions, cufftPlanMany() allows use of more complicated data layouts and batched executions. Execution of a transform of a particular size and type may take several stages of processing. When a plan for the transform is generated, cuFFT derives the internal steps that need to be taken. These steps may include multiple kernel launches, memory copies, and so on. In addition, all the intermediate buffer allocations (on CPU/GPU memory) take place during planning. These buffers are released when the plan is destroyed. In the worst case, the cuFFT Library allocates space for 8*batch*n[0]*..*n[rank-1] cufftComplex or cufftDoubleComplex elements (where batch denotes the number of transforms that will be executed in parallel, rank is the number of dimensions of the input data (see Multidimensional Transforms) and n[] is the array of transform dimensions) for single and double-precision transforms respectively. Depending on the configuration of the plan, less memory may be used. In some specific cases, the temporary space allocations can be as low as 1*batch*n[0]*..*n[rank-1] cufftComplex or cufftDoubleComplex elements. This temporary space is allocated separately for each individual plan when it is created (i.e., temporary space is not shared between the plans).

The next step in using the library is to call an execution function such as cufftExecC2C() (see Parameter cufftType) which will perform the transform with the specifications defined at planning.

One can create a cuFFT plan and perform multiple transforms on different data sets by providing different input and output pointers. Once the plan is no longer needed, the cufftDestroy() function should be called to release the resources allocated for the plan.

Apart from the general complex-to-complex (C2C) transform, cuFFT implements efficiently two other types: real-to-complex (R2C) and complex-to-real (C2R). In many practical applications the input vector is real-valued. It can be easily shown that in this case the output satisfies Hermitian symmetry ( $X_k=X_{N-k}^{\ast }$, where the star denotes complex conjugation). The converse is also true: for complex-Hermitian input the inverse transform will be purely real-valued. cuFFT takes advantage of this redundancy and works only on the first half of the Hermitian vector.

Transform execution functions for single and double-precision are defined separately as:

• cufftExecC2C() / cufftExecZ2Z() - complex-to-complex transforms for single/double precision.
• cufftExecR2C() / cufftExecD2Z() - real-to-complex forward transform for single/double precision.
• cufftExecC2R() / cufftExecZ2D() - complex-to-real inverse transform for single/double precision.

Each of those functions demands different input data layout (see Data Layout for details).

In the cuFFT Library, data layout depends strictly on the configuration and the transform type. In the case of general complex-to-complex transform both the input and output data shall be a cufftComplex/cufftDoubleComplex array in single- and double-precision modes respectively. In C2R mode an input array $(x_1,x_2,\dots ,x_{\lfloor \frac{N}{2}\rfloor +1})$ of only non-redundant complex elements is required. The output array $(X_1,X_2,\dots ,X_N)$ consists of cufftReal/cufftDouble elements in this mode. Finally, R2C demands an input array $(X_1,X_2,\dots ,X_N)$ of real values and returns an array $(x_1,x_2,\dots ,x_{\lfloor \frac{N}{2}\rfloor +1})$ of non-redundant complex elements.

In real-to-complex and complex-to-real transforms the size of input data and the size of output data differ. For out-of-place transforms a separate array of appropriate size is created. For in-place transforms the user can specify one of two supported data layouts: padded ornative(deprecated). The default is padded for FFTW compatibility.

In the padded layout output signals begin at the same memory addresses as the input data. Therefore input data for real-to-complex and output data for complex-to-real must be padded. In the native(deprecated) layout no padding is required and both input and output data are formed as arrays of adequate types and sizes.

Expected sizes of input/output data for 1-d transforms are summarized in the table below:

The real-to-complex transform is implicitly a forward transform. For an in-place real-to-complex transform where FFTW compatible output is desired, the input size must be padded to $\left(\lfloor \frac{N}{2}\rfloor +1\right)$ complex elements. For out-of-place transforms, input and output sizes match the logical transform size $\mathrm{N}$ and the non-redundant size $\lfloor \frac{N}{2}\rfloor +1$, respectively.

The complex-to-real transform is implicitly inverse. For in-place complex-to-real FFTs where FFTW compatible output is selected (default padding mode), the input size is assumed to be $\lfloor \frac{N}{2}\rfloor +1$cufftComplex elements. Note that in-place complex-to-real FFTs may overwrite arbitrary imaginary input point values when non-unit input and output strides are chosen. For out-of-place transforms, input and output sizes match the logical transform non-redundant size $\lfloor \frac{N}{2}\rfloor +1$ and size $\mathrm{N}$, respectively.

Multidimensional DFT map a $d$-dimensional array $x_{\mathbf{n}}$, where $\mathbf{n}=(n_1,n_2,\dots ,n_d)$ into its frequency domain array given by:

cuFFT supports one-dimensional, two-dimensional and three-dimensional transforms, which can all be called by the same cufftExec* functions (see Fourier Transform Types).

Similar to the one-dimensional case, the frequency domain representation of real-valued input data satisfies Hermitian symmetry, defined as: $x_{(n_1,n_2,\dots ,n_d)}=x_{(N_1-n_1,N_2-n_2,\dots ,N_d-n_d)}^{\ast }$.

C2R and R2C algorithms take advantage of this fact by operating only on half of the elements of signal array, namely on: $x_{\mathbf{n}}$ for $\mathbf{n}\in \{1,\dots ,N_1\}\times \dots \times \{1,\dots ,N_{d-1}\}\times \{1,\dots ,\lfloor \frac{N_d}{2}\rfloor +1\}$.

The general rules of data alignment described in Data Layout apply to higher-dimensional transforms. The following table summarizes input and output data sizes for multidimensional DFTs:

For example, static declaration of a three-dimensional array for the output of an out-of-place real-to-complex transform will look like this:

cufftComplex odata[N1][N2][N3/2+1];

The advanced data layout feature allows transforming only a subset of an input array, or outputting to only a portion of a larger data structure. It can be set by calling function:

cufftResult cufftPlanMany(cufftHandle *plan, int rank, int *n, int *inembed,
    int istride, int idist, int *onembed, int ostride,
    int odist, cufftType type, int batch);

Passing inembed or onembed set to NULL is a special case and is equivalent to passing n for each. This is same as the basic data layout and other advanced parameters such as istride are ignored.

If the advanced parameters are to be used, then all of the advanced interface parameters must be specified correctly. Advanced parameters are defined in units of the relevant data type (cufftReal, cufftDoubleReal, cufftComplex, or cufftDoubleComplex).

Advanced layout can be perceived as an additional layer of abstraction above the access to input/output data arrays. An element of coordinates [z][y][x] in signal number b in the batch will be associated with the following addresses in the memory:

• 1D

  input[ b * idist + x * istride]

  output[ b * odist + x * ostride]

• 2D

  input[b * idist + (x * inembed[1] + y) * istride]

  output[b * odist + (x * onembed[1] + y) * ostride]

• 3D

  input[b * idist + ((x * inembed[1] + y) * inembed[2] + z) * istride]

  output[b * odist + ((x * onembed[1] + y) * onembed[2] + z) * ostride]

The istride and ostride parameters denote the distance between two successive input and output elements in the least significant (that is, the innermost) dimension respectively. In a single 1D transform, if every input element is to be used in the transform, istride should be set to $1$; if every other input element is to be used in the transform, then istride should be set to $2$. Similarly, in a single 1D transform, if it is desired to output final elements one after another compactly, ostride should be set to $1$; if spacing is desired between the least significant dimension output data, ostride should be set to the distance between the elements.

The inembed and onembed parameters define the number of elements in each dimension in the input array and the output array respectively. The inembed[rank-1] contains the number of elements in the least significant (innermost) dimension of the input data excluding the istride elements; the number of total elements in the least significant dimension of the input array is then istride*inembed[rank-1]. The inembed[0] or onembed[0] corresponds to the most significant (that is, the outermost) dimension and is effectively ignored since the idist or odist parameter provides this information instead. Note that the size of each dimension of the transform should be less than or equal to the inembed and onembed values for the corresponding dimension, that is n[i] ≤ inembed[i], n[i] ≤ onembed[i], where $i\in \{0,\dots ,rank-1\}$.

The idist and odist parameters indicate the distance between the first element of two consecutive batches in the input and output data.

Every cuFFT plan may be associated with a CUDA stream. Once so associated, all launches of the internal stages of that plan take place through the specified stream. Streaming of cuFFT execution allows for potential overlap between transforms and memory copies. (See the NVIDIA CUDA Programming Guide for more information on streams.) If no stream is associated with a plan, launches take place in stream(0), the default CUDA stream. Note that many plan executions require multiple kernel launches.

cufftSetStream() returns an error in the multiple GPU case as multiple GPU plans perform operations in their own streams.

cuFFT supports using two GPUs connected to a CPU to perform Fourier Transforms whose calculations are distributed across the GPUs. An API has been defined to allow users to write new code or modify existing code to use this functionality.

Some existing functions such as the creation of a plan using cufftCreate() also apply in the two GPU case. Multiple GPU unique routines contain Xt in their name.

The memory on the GPUs is managed by helper functions cufftXtMalloc()/cufftXtFree() and cufftXtMemcpy() using the cudaLibXtDesc descriptor.

Performance is a function of the bandwidth between the GPUs, the computational ability of the individual GPUs, and the type and number of FFT to be performed. The fastest performance is obtained using PCI Express 3.0 between the GPUs and ensuring that both GPUs are on the same switch. Note that two GPU execution is not guaranteed to solve a given size problem in a shorter time than single GPU execution.

The multiple GPU extensions to cuFFT are built on the extensible cuFFT API. The general steps in defining and executing a transform with this API are:

• cufftCreate() - create an empty plan, as in the single GPU case
• cufftXtSetGPUs() - define which GPUs are to be used
• Optional: cufftEstimate{1d,2d,3d,Many}() - estimate the sizes of the work areas required. These are the same functions used in the single GPU case although the definition of the argument workSize reflects the number of GPUs used.
• cufftMakePlan{1d,2d,3d,Many}() - create the plan. These are the same functions used in the single GPU case although the definition of the argument workSize reflects the number of GPUs used.
• Optional: cufftGetSize{1d,2d,3d,Many}() - refined estimate of the sizes of the work areas required. These are the same functions used in the single GPU case although the definition of the argument workSize reflects the number of GPUs used.
• Optional: cufftGetSize() - check workspace size. This is the same function used in the single GPU case although the definition of the argument workSize reflects the number of GPUs used.
• Optional: cufftXtSetWorkArea() - do your own workspace allocation.
• cufftXtMalloc() - allocate descriptor and data on the GPUs
• cufftXtMemcpy() - copy data to the GPUs
• cufftXtExecDescriptorC2C()/cufftXtExecDescriptorZ2Z() - execute the plan
• cufftXtMemcpy() - copy data from the GPUs
• cufftXtFree() - free any memory allocated with cufftXtMalloc()
• cufftDestroy() - free cuFFT plan resources

cuFFT APIs are thread safe as long as different host threads execute FFTs using different plans and the output data are disjoint.

A DFT can be implemented as a matrix vector multiplication that requires $\mathrm{O}({\mathrm{N}}^2)$ operations. However, the cuFFT Library employs the Cooley-Tukey algorithm to reduce the number of required operations to optimize the performance of particular transform sizes. This algorithm expresses the DFT matrix as a product of sparse building block matrices. The cuFFT Library implements the following building blocks: radix-2, radix-3, radix-5, and radix-7. Hence the performance of any transform size that can be factored as $2^a\times 3^b\times 5^c\times 7^d$ (where a, b, c, and d are non-negative integers) is optimized in the cuFFT library. There are also radix-m building blocks for other primes, m, whose value is < 128. When the length cannot be decomposed as multiples of powers of primes from 2 to 127, Bluestein's algorithm is used. Since the Bluestein implementation requires more computations per output point than the Cooley-Tukey implementation, the accuracy of the Cooley-Tukey algorithm is better. The pure Cooley-Tukey implementation has excellent accuracy, with the relative error growing proportionally to ${\log }_2(\mathrm{N})$ , where $\mathrm{N}$ is the transform size in points.

For sizes handled by the Cooley-Tukey code path, the most efficient implementation is obtained by applying the following constraints (listed in order from the most generic to the most specialized constraint, with each subsequent constraint providing the potential of an additional performance improvement).

All cuFFT Library return values except for CUFFT_SUCCESS indicate that the current API call failed and the user should reconfigure to correct the problem. The possible return values are defined as follows:

typedef enum cufftResult_t {
    CUFFT_SUCCESS        = 0,  //  The cuFFT operation was successful
    CUFFT_INVALID_PLAN   = 1,  //  cuFFT was passed an invalid plan handle
    CUFFT_ALLOC_FAILED   = 2,  //  cuFFT failed to allocate GPU or CPU memory
    CUFFT_INVALID_TYPE   = 3,  //  No longer used
    CUFFT_INVALID_VALUE  = 4,  //  User specified an invalid pointer or parameter
    CUFFT_INTERNAL_ERROR = 5,  //  Driver or internal cuFFT library error
    CUFFT_EXEC_FAILED    = 6,  //  Failed to execute an FFT on the GPU
    CUFFT_SETUP_FAILED   = 7,  //  The cuFFT library failed to initialize
    CUFFT_INVALID_SIZE   = 8,  //  User specified an invalid transform size
    CUFFT_UNALIGNED_DATA = 9,  //  No longer used
    CUFFT_INCOMPLETE_PARAMETER_LIST = 10, //  Missing parameters in call
    CUFFT_INVALID_DEVICE = 11, //  Execution of a plan was on different GPU than plan creation
    CUFFT_PARSE_ERROR    = 12, //  Internal plan database error 
    CUFFT_NO_WORKSPACE   = 13  //  No workspace has been provided prior to plan execution
} cufftResult;

Users are encouraged to check return values from cuFFT functions for errors as shown in cuFFT Code Examples.

This API separates handle creation from plan generation. This makes it possible to change plan settings, which may alter the outcome of the plan generation phase, before the plan is actually generated.

During plan execution, cuFFT requires a work area for temporary storage of intermediate results. The cufftEstimate*() calls return an estimate for the size of the work area required, given the specified parameters, and assuming default plan settings. Some problem sizes require much more storage than others. In particular powers of 2 are very efficient in terms of temporary storage. Large prime numbers, however, use different algorithms and may need up to the eight times that of a similarly sized power of 2. These routines return estimated workSize values which may still be smaller than the actual values needed especially for values of n that are not multiples of powers of 2, 3, 5 and 7. More refined values are given by the cufftGetSize*() routines, but these values may still be conservative.

The cufftGetSize*() routines give a more accurate estimate of the work area size required for a plan than the cufftEstimate*() routines as they take into account any plan settings that may have been made. As discussed in the section cuFFT Estimated Size of Work Area, the workSize value(s) returned may be conservative especially for values of n that are not multiples of powers of 2, 3, 5 and 7.

cufftResult 
    cufftGetSize(cufftHandle plan, size_t *workSize);

Once plan generation has been done, either with the original API or the extensible API, this call returns the actual size of the work area required to support the plan. Callers who choose to manage work area allocation within their application must use this call after plan generation, and after any cufftSet*() calls subsequent to plan generation, if those calls might alter the required work space size.

cufftResult 
    cufftDestroy(cufftHandle plan);

Frees all GPU resources associated with a cuFFT plan and destroys the internal plan data structure. This function should be called once a plan is no longer needed, to avoid wasting GPU memory.

cufftResult 
    cufftSetStream(cufftHandle plan, cudaStream_t stream);

Associates a CUDA stream with a cuFFT plan. All kernel launches made during plan execution are now done through the associated stream, enabling overlap with activity in other streams (e.g. data copying). The association remains until the plan is destroyed or the stream is changed with another call to cufftSetStream(). This call will return an error for multiple GPU plans.

cufftResult 
    cufftGetVersion(int *version);

Returns the version number of cuFFT.

cufftResult 
    cufftSetCompatibilityMode(cufftHandle plan, cufftCompatibility mode);

Configures the layout of cuFFT output in FFTW-compatible modes. When desired, FFTW compatibility can be configured for padding only, for asymmetric complex inputs only, or for full compatibility. If the cufftSetCompatibilityMode() API fails, later cufftExecute*() calls are not guaranteed to work.

cuFFT Library defines FFTW compatible data layouts using the following enumeration of values. See FFTW Compatibility Mode for more details.

typedef enum cufftCompatibility_t {
    CUFFT_COMPATIBILITY_NATIVE       = 0, // Compact data in native format 
                                          // (deprecated. Use
                                          // CUFFT_COMPATIBILITY_FFTW_PADDING)
    CUFFT_COMPATIBILITY_FFTW_PADDING = 1, // FFTW-compatible alignment 
                                          // (default value)
    CUFFT_COMPATIBILITY_FFTW_ASYMMETRIC = 2, // Waives the C2R symmetry 
                                             //requirement
    CUFFT_COMPATIBILITY_FFTW_ALL        = 3
} cufftCompatibility;

In this example a one-dimensional complex-to-complex transform is applied to the input data. Afterwards an inverse transform is performed on the computed frequency domain representation.

#define NX 256
#define BATCH 1

cufftHandle plan;
cufftComplex *data;
cudaMalloc((void**)&data, sizeof(cufftComplex)*NX*BATCH);
if (cudaGetLastError() != cudaSuccess){
	fprintf(stderr, "Cuda error: Failed to allocate\n");
	return;	
}

if (cufftPlan1d(&plan, NX, CUFFT_C2C, BATCH) != CUFFT_SUCCESS){
	fprintf(stderr, "CUFFT error: Plan creation failed");
	return;	
}	

...

/* Note:
 *  Identical pointers to input and output arrays implies in-place transformation
 */

if (cufftExecC2C(plan, data, data, CUFFT_FORWARD) != CUFFT_SUCCESS){
	fprintf(stderr, "CUFFT error: ExecC2C Forward failed");
	return;	
}

if (cufftExecC2C(plan, data, data, CUFFT_INVERSE) != CUFFT_SUCCESS){
	fprintf(stderr, "CUFFT error: ExecC2C Inverse failed");
	return;	
}

/* 
 *  Results may not be immediately available so block device until all 
 *  tasks have completed
 */

if (cudaDeviceSynchronize() != cudaSuccess){
	fprintf(stderr, "Cuda error: Failed to synchronize\n");
	return;	
}	

/* 
 *  Divide by number of elements in data set to get back original data
 */

...
	
cufftDestroy(plan);
cudaFree(data);

In this example a one-dimensional real-to-complex transform is applied to the input data.

#define NX 256
#define BATCH 1

cufftHandle plan;
cufftComplex *data;
cudaMalloc((void**)&data, sizeof(cufftComplex)*(NX/2+1)*BATCH);
if (cudaGetLastError() != cudaSuccess){
	fprintf(stderr, "Cuda error: Failed to allocate\n");
	return;	
}

if (cufftPlan1d(&plan, NX, CUFFT_R2C, BATCH) != CUFFT_SUCCESS){
	fprintf(stderr, "CUFFT error: Plan creation failed");
	return;	
}	

...

/* Use the CUFFT plan to transform the signal in place. */
if (cufftExecR2C(plan, (cufftReal*)data, data) != CUFFT_SUCCESS){
	fprintf(stderr, "CUFFT error: ExecC2C Forward failed");
	return;	
}

if (cudaDeviceSynchronize() != cudaSuccess){
	fprintf(stderr, "Cuda error: Failed to synchronize\n");
	return;	
}

...

cufftDestroy(plan);
cudaFree(data);

In this example a two-dimensional complex-to-real transform is applied to the input data arranged according to the requirements of the default FFTW padding mode.

#define NX 256
#define NY 128
#define NRANK 2
#define BATCH 1

cufftHandle plan;
cufftComplex *data;
int n[NRANK] = {NX, NY};

cudaMalloc((void**)&data, sizeof(cufftComplex)*NX*(NY/2+1));
if (cudaGetLastError() != cudaSuccess){
	fprintf(stderr, "Cuda error: Failed to allocate\n");
	return;	
}

/* Create a 2D FFT plan. */
if (cufftPlanMany(&plan, NRANK, n,
				  NULL, 1, 0,
				  NULL, 1, 0,
				  CUFFT_C2R,BATCH) != CUFFT_SUCCESS){
	fprintf(stderr, "CUFFT Error: Unable to create plan\n");
	return;	
}

...

if (cufftExecC2R(plan, data, data) != CUFFT_SUCCESS){
	fprintf(stderr, "CUFFT Error: Unable to execute plan\n");
	return;		
}

if (cudaDeviceSynchronize() != cudaSuccess){
  	fprintf(stderr, "Cuda error: Failed to synchronize\n");
   	return;
}	

...

cufftDestroy(plan);
cudaFree(data);

In this example a three-dimensional complex-to-complex transform is applied to the input data.

#define NX 64
#define NY 128
#define NZ 128
#define BATCH 10
#define NRANK 3

cufftHandle plan;
cufftComplex *data;
int n[NRANK] = {NX, NY, NZ};

cudaMalloc((void**)&data, sizeof(cufftComplex)*NX*NY*NZ*BATCH);
if (cudaGetLastError() != cudaSuccess){
	fprintf(stderr, "Cuda error: Failed to allocate\n");
	return;	
}

/* Create a 3D FFT plan. */
if (cufftPlanMany(&plan, NRANK, n, 
				  NULL, 1, NX*NY*NZ, // *inembed, istride, idist 
				  NULL, 1, NX*NY*NZ, // *onembed, ostride, odist
				  CUFFT_C2C, BATCH) != CUFFT_SUCCESS){
	fprintf(stderr, "CUFFT error: Plan creation failed");
	return;	
}	

/* Use the CUFFT plan to transform the signal in place. */
if (cufftExecC2C(plan, data, data, CUFFT_FORWARD) != CUFFT_SUCCESS){
	fprintf(stderr, "CUFFT error: ExecC2C Forward failed");
	return;	
}

if (cudaDeviceSynchronize() != cudaSuccess){
	fprintf(stderr, "Cuda error: Failed to synchronize\n");
	return;	
}	

...
	
cufftDestroy(plan);
cudaFree(data);

In this example a two-dimensional complex-to-complex transform is applied to the input data arranged according to the requirements the advanced layout.

#define NX 128
#define NY 256
#define BATCH 10
#define NRANK 2

/* Advanced interface parameters, arbitrary strides */
#define ISTRIDE 2   // distance between successive input elements in innermost dimension
#define OSTRIDE 1   // distance between successive output elements in innermost dimension
#define IX (NX+2)
#define IY (NY+1)
#define OX (NX+3)
#define OY (NY+4)
#define IDIST (IX*IY*ISTRIDE+3) // distance between first element of two consecutive signals in a batch of input data
#define ODIST (OX*OY*OSTRIDE+5) // distance between first element of two consecutive signals in a batch of output data

cufftHandle plan;
cufftComplex *idata, *odata;
int isize = IDIST * BATCH;
int osize = ODIST * BATCH;
int n[NRANK] = {NX, NY};
int inembed[NRANK] = {IX, IY}; // pointer that indicates storage dimensions of input data
int onembed[NRANK] = {OX, OY}; // pointer that indicates storage dimensions of output data

cudaMalloc((void **)&idata, sizeof(cufftComplex)*isize);
cudaMalloc((void **)&odata, sizeof(cufftComplex)*osize);
if (cudaGetLastError() != cudaSuccess){
	fprintf(stderr, "Cuda error: Failed to allocate\n");
	return;	
}

/* Create a batched 2D plan */
if (cufftPlanMany(&plan, NRANK, n,
				  inembed,ISTRIDE,IDIST,
				  onembed,OSTRIDE,ODIST,
				  CUFFT_C2C,BATCH) != CUFFT_SUCCESS){
	fprintf(stderr, "CUFFT Error: Unable to create plan\n");
	return;	
}

...

/* Execute the transform out-of-place */
if (cufftExecC2C(plan, idata, odata, CUFFT_FORWARD) != CUFFT_SUCCESS){
	fprintf(stderr, "CUFFT Error: Failed to execute plan\n");
	return;		
}

if (cudaDeviceSynchronize() != cudaSuccess){
  	fprintf(stderr, "Cuda error: Failed to synchronize\n");
   	return;
}

...

cufftDestroy(plan);
cudaFree(idata);
cudaFree(odata);

In this example a three-dimensional complex-to-complex transform is applied to the input data using two GPUs.

// Demonstrate how to use CUFFT to perform 3-d FFTs using 2 GPUs
//
// cufftCreate() - Create an empty plan
    cufftHandle plan_input; cufftResult result;
    result = cufftCreate(&plan_input);
    if (result != CUFFT_SUCCESS) { printf ("*Create failed\n"); return; }
//
// cufftXtSetGPUs() - Define which GPUs to use
    int nGPUs = 2, whichGPUs[2];
    whichGPUs[0] = 0; whichGPUs[1] = 1;
    result = cufftXtSetGPUs (plan_input, nGPUs, whichGPUs);
    if (result != CUFFT_SUCCESS) { printf ("*XtSetGPUs failed\n"); return; }
//
// Initialize FFT input data
    size_t worksize[2];
    cufftComplex *host_data_input, *host_data_output;
    int nx = 64, ny = 128, nz = 32;
    int size_of_data = sizeof(cufftComplex) * nx * ny * nz;
    host_data_input = malloc(size_of_data);
    if (host_data_input == NULL) { printf ("malloc failed\n"); return; }
    host_data_output = malloc(size_of_data);
    if (host_data_output == NULL) { printf ("malloc failed\n"); return; }
    initialize_3d_data (nx, ny, nz, host_data_input, host_data_output);
//
// cufftMakePlan3d() - Create the plan
    result = cufftMakePlan3d (plan_input, nz, ny, nx, CUFFT_C2C, worksize);
    if (result != CUFFT_SUCCESS) { printf ("*MakePlan* failed\n"); return; }
//
// cufftXtMalloc() - Malloc data on multiple GPUs
    cudaLibXtDesc *device_data_input;
    result = cufftXtMalloc (plan_input, (void*)&device_data_input,
        CUFFT_XT_FORMAT_INPLACE);
    if (result != CUFFT_SUCCESS) { printf ("*XtMalloc failed\n"); return; }
//
// cufftXtMemcpy() - Copy data from host to multiple GPUs
    result = cufftXtMemcpy (plan_input, device_data_input,
        host_data_input, CUFFT_COPY_HOST_TO_DEVICE);
    if (result != CUFFT_SUCCESS) { printf ("*XtMemcpy failed\n"); return; }
//
// cufftXtExecDescriptorC2C() - Execute FFT on multiple GPUs
    result = cufftXtExecDescriptorC2C (plan_input, device_data_input,
        device_data_input, CUFFT_FORWARD);
    if (result != CUFFT_SUCCESS) { printf ("*XtExec* failed\n"); return; }
//
// cufftXtMemcpy() - Copy data from multiple GPUs to host
    result = cufftXtMemcpy (plan_input, host_data_output,
        device_data_input, CUFFT_COPY_DEVICE_TO_HOST);
    if (result != CUFFT_SUCCESS) { printf ("*XtMemcpy failed\n"); return; }
//
// Print output and check results
    int output_return = output_3d_results (nx, ny, nz,
        host_data_input, host_data_output);
    if (output_return != 0) { return; }
//
// cufftXtFree() - Free GPU memory
    result = cufftXtFree(device_data_input);
    if (result != CUFFT_SUCCESS) { printf ("*XtFree failed\n"); return; }
//
// cufftDestroy() - Destroy FFT plan
    result = cufftDestroy(plan_input);
    if (result != CUFFT_SUCCESS) { printf ("*Destroy failed: code\n"); return; }
    free(host_data_input); free(host_data_output);

In this example a one-dimensional complex-to-complex transform is applied to the input data using two GPUs. The output data is in natural order in GPU memory.

// Demonstrate how to use CUFFT to perform 1-d FFTs using 2 GPUs 
// Output on the GPUs is in natural output
// Function return codes should be checked for errors in actual code
//
// cufftCreate() - Create an empty plan
    cufftHandle plan_input; cufftResult result;
    result = cufftCreate(&plan_input);
//
// cufftXtSetGPUs() - Define which GPUs to use
    int nGPUs = 2, whichGPUs[2];
    whichGPUs[0] = 0; whichGPUs[1] = 1;
    result = cufftXtSetGPUs (plan_input, nGPUs, whichGPUs);
//
// Initialize FFT input data
    size_t worksize[2];
    cufftComplex *host_data_input, *host_data_output;
    int nx = 1024, batch = 1, rank = 1, n[1];
    int inembed[1], istride, idist, onembed[1], ostride, odist;
    n[0] = nx;
    int size_of_data = sizeof(cufftComplex) * nx * batch;
    host_data_input = malloc(size_of_data);
    host_data_output = malloc(size_of_data);
    initialize_1d_data (nx, batch, rank, n, inembed, &istride, &idist,
        onembed, &ostride, &odist, host_data_input, host_data_output);
//
// cufftMakePlanMany() - Create the plan
    result = cufftMakePlanMany (plan_input, rank, n, inembed, istride, idist,
        onembed, ostride, odist, CUFFT_C2C, batch, worksize);
//
// cufftXtMalloc() - Malloc data on multiple GPUs
    cudaLibXtDesc *device_data_input, *device_data_output;
    result = cufftXtMalloc (plan_input, (void*)&device_data_input,
        CUFFT_XT_FORMAT_INPLACE);
    result = cufftXtMalloc (plan_input, (void*)&device_data_output,
        CUFFT_XT_FORMAT_INPLACE);
//
// cufftXtMemcpy() - Copy data from host to multiple GPUs
    result = cufftXtMemcpy (plan_input, device_data_input,
        host_data_input, CUFFT_COPY_HOST_TO_DEVICE);
//
// cufftXtExecDescriptorC2C() - Execute FFT on multiple GPUs
    result = cufftXtExecDescriptorC2C (plan_input, device_data_input,
        device_data_input, CUFFT_FORWARD);
//
// cufftXtMemcpy() - Copy the data to natural order on GPUs
    result = cufftXtMemcpy (plan_input, device_data_output,
        device_data_input, CUFFT_COPY_DEVICE_TO_DEVICE);
//
// cufftXtMemcpy() - Copy natural order data from multiple GPUs to host
    result = cufftXtMemcpy (plan_input, host_data_output,
        device_data_output, CUFFT_COPY_DEVICE_TO_HOST);
//
// Print output and check results
    int output_return = output_1d_results (nx, batch,
        host_data_input, host_data_output);
//
// cufftXtFree() - Free GPU memory
    result = cufftXtFree(device_data_input);
    result = cufftXtFree(device_data_output);
//
// cufftDestroy() - Destroy FFT plan
    result = cufftDestroy(plan_input);
    free(host_data_input); free(host_data_output);

In this example a one-dimensional convolution is calculated using complex-to-complex transforms.

//
// Demonstrate how to use CUFFT to perform a convolution using 1-d FFTs and
// 2 GPUs. The forward FFTs use both GPUs, while the inverse FFT uses one.
// Function return codes should be checked for errors in actual code.
//
// cufftCreate() - Create an empty plan
    cufftResult result; cudaError_t cuda_status;
    cufftHandle plan_forward_2_gpus, plan_inverse_1_gpu;
    result = cufftCreate(&plan_forward_2_gpus);
    result = cufftCreate(&plan_inverse_1_gpu);
//
// cufftXtSetGPUs() - Define which GPUs to use
    int nGPUs = 2, whichGPUs[2];
    whichGPUs[0] = 0; whichGPUs[1] = 1;
    result = cufftXtSetGPUs (plan_forward_2_gpus, nGPUs, whichGPUs);
//
// Initialize FFT input data
    size_t worksize[2];
    cufftComplex *host_data_input, *host_data_output;
    int nx = 1048576, batch = 2, rank = 1, n[1];
    int inembed[1], istride, idist, onembed[1], ostride, odist;
    n[0] = nx;
    int size_of_one_set = sizeof(cufftComplex) * nx;
    int size_of_data = size_of_one_set * batch;
    host_data_input = (cufftComplex*)malloc(size_of_data);
    host_data_output = (cufftComplex*)malloc(size_of_one_set);
    initialize_1d_data (nx, batch, rank, n, inembed, &istride, &idist,
        onembed, &ostride, &odist, host_data_input, host_data_output);
//
// cufftMakePlanMany(), cufftPlan1d - Create the plans
    result = cufftMakePlanMany (plan_forward_2_gpus, rank, n, inembed,
        istride, idist, onembed, ostride, odist, CUFFT_C2C, batch, worksize);
    result = cufftPlan1d (&plan_inverse_1_gpu, nx, CUFFT_C2C, 1);
//
// cufftXtMalloc(), cudaMallocHost - Allocate data for GPUs
    cudaLibXtDesc *device_data_input; cufftComplex *GPU0_data_from_GPU1;
    result = cufftXtMalloc (plan_forward_2_gpus, &device_data_input,
        CUFFT_XT_FORMAT_INPLACE);
    int device0 = device_data_input->descriptor->GPUs[0];
    cudaSetDevice(device0) ;
    cuda_status = cudaMallocHost ((void**)&GPU0_data_from_GPU1,size_of_one_set);
//
// cufftXtMemcpy() - Copy data from host to multiple GPUs
    result = cufftXtMemcpy (plan_forward_2_gpus, device_data_input,
        host_data_input, CUFFT_COPY_HOST_TO_DEVICE);
//
// cufftXtExecDescriptorC2C() - Execute forward FFTs on multiple GPUs
    result = cufftXtExecDescriptorC2C (plan_forward_2_gpus, device_data_input,
        device_data_input, CUFFT_FORWARD);
//
// cudaMemcpy result from GPU1 to GPU0
    cufftComplex *device_data_on_GPU1;
    device_data_on_GPU1 = (cufftComplex*)
        (device_data_input->descriptor->data[1]);
    cuda_status = cudaMemcpy (GPU0_data_from_GPU1, device_data_on_GPU1,
        size_of_one_set, cudaMemcpyDeviceToDevice);
//
// Continued on next page
//

//
// Demonstrate how to use CUFFT to perform a convolution using 1-d FFTs and
// 2 GPUs. The forward FFTs use both GPUs, while the inverse FFT uses one.
// Function return codes should be checked for errors in actual code.
//
// Part 2
//
// Multiply results and scale output
    cufftComplex *device_data_on_GPU0;
    device_data_on_GPU0 = (cufftComplex*)
        (device_data_input->descriptor->data[0]);
    cudaSetDevice(device0) ;
    ComplexPointwiseMulAndScale<<<32, 256>>>((cufftComplex*)device_data_on_GPU0,
        (cufftComplex*) GPU0_data_from_GPU1, nx);
//   
// cufftExecC2C() - Execute inverse FFT on one GPU
    result = cufftExecC2C (plan_inverse_1_gpu, GPU0_data_from_GPU1,
        GPU0_data_from_GPU1, CUFFT_INVERSE);
//
// cudaMemcpy() - Copy results from GPU0 to host
    cuda_status = cudaMemcpy(host_data_output, GPU0_data_from_GPU1,
        size_of_one_set, cudaMemcpyDeviceToHost);
//
// Print output and check results
    int output_return = output_1d_results (nx, batch,
        host_data_input, host_data_output);
//
// cufftDestroy() - Destroy FFT plans
    result = cufftDestroy(plan_forward_2_gpus);
    result = cufftDestroy(plan_inverse_1_gpu);
//
// cufftXtFree(), cudaFreeHost(), free() - Free GPU and host memory
    result = cufftXtFree(device_data_input);
    cuda_status = cudaFreeHost (GPU0_data_from_GPU1);
    free(host_data_input); free(host_data_output);

// 
// Utility routine to perform complex pointwise multiplication with scaling
__global__ void ComplexPointwiseMulAndScale
    (cufftComplex *a, cufftComplex *b, int size)
{
    const int numThreads = blockDim.x * gridDim.x;
    const int threadID = blockIdx.x * blockDim.x + threadIdx.x;
    float scale = 1.0f / (float)size;
    cufftComplex c;
    for (int i = threadID; i < size; i += numThreads)
    {
        c = cuCmulf(a[i], b[i]);
        b[i] = make_cuFloatComplex(scale*cuCrealf(c), scale*cuCimagf(c));
    }
    return;
}

For batches of transforms, each individual transform is executed on a single GPU. The first $⌊\frac{\mathrm{batch\; size\; +\; 1}}{2}⌋$ transforms are executed on GPU 0, and the remainder are executed on GPU 1. This approach removes the need for data exchange between the GPUs, and results in nearly perfect scaling for even batch sizes.

Single transforms performed on multiple GPUs require the data to be divided between the GPUs. Then execution takes place in phases. First, for 2D and 3D transforms, each GPU does half of the transform in (rank - 1) dimensions. Then data are exchanged between the GPUs so that the final dimension can be processed.

For 2D transforms using an array declared in C as data[x][y], the surface is distributed prior to the transform such that each GPU receives half the surface, with dimensions [x/2][y]. GPU 0 receives all y elements for x indices ranging from 0...(x/2-1). GPU 1 receives all y elements for x indices ranging from (x/2)...(x-1). After the transform, each GPU again has half the surface, but with dimensions [x][y/2]. GPU 0 has all x elements for y indices ranging from 0...(y/2-1). GPU 1 has all x elements for y indices ranging from (y/2)...(y-1).

For 3D transforms using an array declared in C as data[x][y][z], the volume is distributed prior to the transform such that each GPU receives half the volume, with dimensions [x/2][y][z]. GPU 0 receives the [y][z] planes with x indices ranging from 0...(x/2-1). GPU 1 receives the [y][z] planes with x indices ranging from (x/2)...(x-1). After the transform, each GPU again has half the volume, but with dimensions [x][y/2][z]. GPU 0 has the [x][z] planes with y indices ranging from 0...(y/2-1). GPU 1 has the [x][z] planes with y indices ranging from (y/2)...(y-1).

For 1D transforms, the initial distribution of data to the GPUs is similar to the 2D and 3D cases. For a transform of dimension x, GPU 0 receives data ranging from 0...(x/2-1). GPU 1 receives data ranging from (x/2)...(x-1). The partitioning of the work between the GPUs results in a more complex output arrangement than in the 2D and 3D cases.

cuFFT performs multiple GPU 1D transforms by decomposing the transform size into factors Factor1 and Factor2, and treating the data as a grid of size Factor1 x Factor2. The four steps done to calculate the 1D FFT are: Factor1 transforms of size Factor2, data exchange between the GPUs, a pointwise twiddle multiplication, and Factor2 transforms of size Factor1. To gain efficiency by overlapping computation with data exchange, cuFFT breaks the whole transform into independent segments or strings, which can be processed while others are in flight. A side effect of this algorithm is that the output of the transform is not in linear order. The output in GPU memory is in strings, each of which is composed of Factor2 substrings of equal size. Each substring contains contiguous results starting Factor1 elements subsequent to start of the previous substring. Each string starts substring size elements after the start of the previous string. The strings appear in order, the first half on GPU 0, and the second half on GPU 1. See the example below:

transform size = 1024
number of strings = 8
Factor1 = 64
Factor2 = 16
substrings per string for output layout is Factor2 (16)
string size = 1024/8 = 128
substring size = 128/16 = 8
stride between substrings = 1024/16 = Factor1 (64)

On GPU 0:
string 0 has substrings with indices 0...7   64...71   128...135 ... 960...967
string 1 has substrings with indices 8...15  72...79   136...143 ... 968...975
...
On GPU 1:
string 4 has substrings with indices 32...39  96...103 160...167 ... 992...999
...
string 7 has substrings with indices 56...63 120...127 184...191 ... 1016...1023
    

The cufftXtQueryPlan API allows the caller to retrieve a structure containing the number of strings, the decomposition factors, and (in the case of power of 2 size) some useful mask and shift elements. The example below shows how cufftXtQueryPlan is invoked. It also shows how to translate from an index in the host input array to the corresponding index on the device, and vice versa.

/* 
 * These routines demonstrate the use of cufftXtQueryPlan to get the 1D 
 * factorization and convert between permuted and linear indexes.
 */
/*
 * Set up a 1D plan that will execute on GPU 0 and GPU1, and query
 * the decomposition factors
 */
int main(int argc, char **argv){
    cufftHandle plan;
    cufftResult stat;
    int whichGPUs[2] = { 0, 1 };
    cufftXt1dFactors factors;
    stat = cufftCreate( &plan );
    if (stat != CUFFT_SUCCESS) {
        printf("Create error %d\n",stat);
        return 1;
    }
    stat = cufftXtSetGPUs( plan, 2, whichGPUs );
    if (stat != CUFFT_SUCCESS) {
        printf("SetGPU error %d\n",stat);
        return 1;
    }
    stat = cufftMakePlan1d( plan, size, CUFFT_C2C, 1, workSizes );
    if (stat != CUFFT_SUCCESS) {
        printf("MakePlan error %d\n",stat);
        return 1;
    }
    stat = cufftXtQueryPlan( plan, (void *) &factors, CUFFT_QUERY_1D_FACTORS );
    if (stat != CUFFT_SUCCESS) {
        printf("QueryPlan error %d\n",stat);
        return 1;
    }
    printf("Factor 1 %zd, Factor2 %zd\n",factors.factor1,factors.factor2);
    cufftDestroy(plan);
    return 0;
}

    

/* 
 * Given an index into a permuted array, and the GPU index return the 
 * corresponding linear index from the beginning of the input buffer.
 * 
 * Parameters:
 *      factors     input:  pointer to cufftXt1dFactors as returned by 
 *                          cufftXtQueryPlan
 *      permutedIx  input:  index of the desired element in the device output
 *                          array
 *      linearIx    output: index of the corresponding input element in the  
 *                          host array
 *      GPUix       input:  index of the GPU containing the desired element
 */
cufftResult permuted2Linear( cufftXt1dFactors * factors,
                             size_t permutedIx,
                             size_t *linearIx,
                             int GPUIx ) {
    size_t indexInSubstring;
    size_t whichString;
    size_t whichSubstring;
    // the low order bits of the permuted index match those of the linear index
    indexInSubstring = permutedIx & factors->substringMask;
    // the next higher bits are the substring index
    whichSubstring = (permutedIx >> factors->substringShift) &
                      factors->factor2Mask;
    // the next higher bits are the string index on this GPU
    whichString = (permutedIx >> factors->stringShift) & factors->stringMask;
    // now adjust the index for the second GPU
    if (GPUIx) {
        whichString += factors->stringCount/2;
    }
    // linear index low order bits are the same
    // next higher linear index bits are the string index
    *linearIx = indexInSubstring + ( whichString << factors->substringShift );
    // next higher bits of linear address are the substring index
    *linearIx += whichSubstring << factors->factor1Shift;
    return CUFFT_SUCCESS;
}

    

/* 
 * Given a linear index into a 1D array, return the GPU containing the permuted 
 * result, and index from the start of the data buffer for that element.
 * 
 * Parameters:
 *      factors     input:  pointer to cufftXt1dFactors as returned by 
 *                          cufftXtQueryPlan
 *      linearIx    input:  index of the desired element in the host input 
 *                          array
 *      permutedIx  output: index of the corresponding result in the device 
 *                          output array
 *      GPUix       output: index of the GPU containing the result
 */
cufftResult linear2Permuted( cufftXt1dFactors * factors,
                             size_t linearIx,
                             size_t *permutedIx,
                             int *GPUIx ) {
    size_t indexInSubstring;
    size_t whichString;
    size_t whichSubstring;
    size_t whichStringMask;
    int whichStringShift;
    if (linearIx >= factors->size) {
        return CUFFT_INVALID_VALUE;
    }
    // get a useful additional mask and shift count
    whichStringMask = factors->stringCount -1;
    whichStringShift = (factors->factor1Shift + factors->factor2Shift) -
                        factors->stringShift ;
    // the low order bits identify the index within the substring
    indexInSubstring = linearIx & factors->substringMask;
    // first determine which string has our linear index.
    // the low order bits indentify the index within the substring.
    // the next higher order bits identify which string.
    whichString = (linearIx >> factors->substringShift) & whichStringMask;
    // the first stringCount/2 strings are in the first GPU, 
    // the rest are in the second.
    *GPUIx = whichString/(factors->stringCount/2);
    // next determine which substring within the string has our index
    // the substring index is in the next higher order bits of the index
    whichSubstring = (linearIx >>(factors->substringShift + whichStringShift)) &
                      factors->factor2Mask;
    // now we can re-assemble the index 
    *permutedIx = indexInSubstring;
    *permutedIx += whichSubstring << factors->substringShift;
    if ( !*GPUIx ) {
        *permutedIx += whichString << factors->stringShift;
    } else {
        *permutedIx += (whichString - (factors->stringCount/2) ) <<
                        factors->stringShift;
    }
    return CUFFT_SUCCESS;
}