# FFT Code Examples

This section presents code examples of using the FFT interface functions described in “Fourier Transform Functions”. Here are the examples of two one-dimensional computations. These examples use the default settings for all of the configuration parameters, which are specified in “Configuration Settings”.

## One-dimensional In-place FFT (Fortran Interface)

! Fortran example.

! 1D complex to complex, and real to conjugate even

Use MKL_DFTI

Complex :: X(32)

Real :: Y(34)

type(DFTI_DESCRIPTOR), POINTER :: My_Desc1_Handle, My_Desc2_Handle

Integer :: Status

!...put input data into X(1),...,X(32); Y(1),...,Y(32)

! Perform a complex to complex transform

Status = DftiCreateDescriptor( My_Desc1_Handle, DFTI_SINGLE,&

DFTI_COMPLEX, 1, 32 )

Status = DftiCommitDescriptor( My_Desc1_Handle )

Status = DftiComputeForward( My_Desc1_Handle, X )

Status = DftiFreeDescriptor(My_Desc1_Handle)

! result is given by {X(1),X(2),...,X(32)}

! Perform a real to complex conjugate even transform

Status = DftiCreateDescriptor(My_Desc2_Handle, DFTI_SINGLE,&

DFTI_REAL, 1, 32)

Status = DftiCommitDescriptor(My_Desc2_Handle)

Status = DftiComputeForward(My_Desc2_Handle, Y)

Status = DftiFreeDescriptor(My_Desc2_Handle)

! result is given in CCS format.

## One-dimensional Out-of-place FFT (Fortran Interface)

! Fortran example.

! 1D complex to complex, and real to conjugate even

Use MKL_DFTI

Complex :: X_in(32)

Complex :: X_out(32)

Real :: Y_in(32)

Real :: Y_out(34)

type(DFTI_DESCRIPTOR), POINTER :: My_Desc1_Handle, My_Desc2_Handle

Integer :: Status

...put input data into X_in(1),...,X_in(32); Y_in(1),...,Y_in(32)

! Perform a complex to complex transform

Status = DftiCreateDescriptor( My_Desc1_Handle, DFTI_SINGLE,

DFTI_COMPLEX, 1, 32 )

Status = DftiSetValue( My_Desc1_Handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE)

Status = DftiCommitDescriptor( My_Desc1_Handle )

Status = DftiComputeForward( My_Desc1_Handle, X_in, X_out )

Status = DftiFreeDescriptor(My_Desc1_Handle)

! result is given by {X_out(1),X_out(2),...,X_out(32)}

! Perform a real to complex conjugate even transform

Status = DftiCreateDescriptor(My_Desc2_Handle, DFTI_SINGLE,

DFTI_REAL, 1, 32)

Status = DftiSetValue( My_Desc2_Handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE)

Status = DftiCommitDescriptor(My_Desc2_Handle)

Status = DftiComputeForward(My_Desc2_Handle, Y_in, Y_out)

Status = DftiFreeDescriptor(My_Desc2_Handle)

! result is given by Y_out in CCS format.

## One-dimensional In-place FFT (C Interface)

/* C example, float _Complex is defined in C9X */

#include "mkl_dfti.h"

float _Complex x[32];

float y[34];

DFTI_DESCRIPTOR *my_desc1_handle,*my_desc2_handle;

/* .... or alternatively

dfti_descriptor_handle my_desc1_handle, my_desc2_handle; */

long status;

//...put input data into x[0],...,x[31]; y[0],...,y[31]

status = DftiCreateDescriptor( &my_desc1_handle, DFTI_SINGLE,

DFTI_COMPLEX, 1, 32);

status = DftiCommitDescriptor( my_desc1_handle );

status = DftiComputeForward( my_desc1_handle, x);

status = DftiFreeDescriptor(&my_desc1_handle);

/* result is x[0], ..., x[31]*/

status = DftiCreateDescriptor( &my_desc2_handle, DFTI_SINGLE,

DFTI_REAL, 1, 32);

status = DftiCommitDescriptor( my_desc2_handle);

status = DftiComputeForward( my_desc2_handle, y);

status = DftiFreeDescriptor(&my_desc2_handle);

/* result is given in CCS format*/

## One-dimensional Out-of-place FFT (C Interface)

/* C example, float _Complex is defined in C9X */

#include "mkl_dfti.h"

float _Complex x_in[32];

float _Complex x_out[32];

float y_in[32];

float y_out[34];

DFTI_DESCRIPTOR *my_desc1_handle,*my_desc2_handle;

/* .... or alternatively

dfti_descriptor_handle my_desc1_handle, my_desc2_handle; */

long status;

//...put input data into x_in[0],...,x_in[31]; y_in[0],...,y_in[31]

status = DftiCreateDescriptor( &my_desc1_handle, DFTI_SINGLE,

DFTI_COMPLEX, 1, 32);

status = DftiSetValue( my_desc1_handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);

status = DftiCommitDescriptor( my_desc1_handle );

status = DftiComputeForward( my_desc1_handle, x_in, x_out);

status = DftiFreeDescriptor(&my_desc1_handle);

/* result is x_out[0], ..., x_out[31]*/

status = DftiCreateDescriptor( &my_desc2_handle, DFTI_SINGLE,

DFTI_REAL, 1, 32);

Status = DftiSetValue( My_Desc2_Handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);

status = DftiCommitDescriptor( my_desc2_handle);

status = DftiComputeForward( my_desc2_handle, y_in, y_out);

status = DftiFreeDescriptor(&my_desc2_handle);

/* result is given by y_out in CCS format*

## Two-dimensional FFT (Fortran Interface)

The following is an example of two simple two-dimensional transforms. Notice that the data and result parameters in computation functions are all declared as assumed-size rank-1 array DIMENSION(0:*). Therefore two-dimensional array must be transformed to one-dimensional array by EQUIVALENCE statement or other facilities of Fortran.

! Fortran example.

! 2D complex to complex, and real to conjugate even

Use MKL_DFTI

Complex :: X_2D(32,100)

Real :: Y_2D(34, 102)

Complex :: X(3200)

Real :: Y(3468)

Equivalence (X_2D, X)

Equivalence (Y_2D, Y)

type(DFTI_DESCRIPTOR), POINTER :: My_Desc1_Handle, My_Desc2_Handle

Integer :: Status, L(2)

!...put input data into X_2D(j,k), Y_2D(j,k), 1<=j=32,1<=k<=100

!...set L(1) = 32, L(2) = 100

!...the transform is a 32-by-100

! Perform a complex to complex transform

Status = DftiCreateDescriptor( My_Desc1_Handle, DFTI_SINGLE,&

DFTI_COMPLEX, 2, L)

Status = DftiCommitDescriptor( My_Desc1_Handle)

Status = DftiComputeForward( My_Desc1_Handle, X)

Status = DftiFreeDescriptor(My_Desc1_Handle)

! result is given by X_2D(j,k), 1<=j<=32, 1<=k<=100

! Perform a real to complex conjugate even transform

Status = DftiCreateDescriptor( My_Desc2_Handle, DFTI_SINGLE,&

DFTI_REAL, 2, L)

Status = DftiCommitDescriptor( My_Desc2_Handle)

Status = DftiComputeForward( My_Desc2_Handle, Y)

Status = DftiFreeDescriptor(My_Desc2_Handle)

! result is given by the complex value z(j,k) 1<=j<=32; 1<=k<=100

! and is stored in CCS format

## Two-dimensional FFT (C Interface)

/* C99 example */

#include "mkl_dfti.h"

float _Complex x[32][100];

float y[34][102];

DFTI_DESCRIPTOR *my_desc1_handle, my_desc2_handle;

long status, l[2];

//...put input data into x[j][k] 0<=j<=31, 0<=k<=99

//...put input data into y[j][k] 0<=j<=31, 0<=k<=99

l[0] = 32; l[1] = 100;

status = DftiCreateDescriptor( &my_desc1_handle, DFTI_SINGLE,

DFTI_COMPLEX, 2, l);

status = DftiCommitDescriptor( my_desc1_handle);

status = DftiComputeForward( my_desc1_handle, x);

status = DftiFreeDescriptor(&my_desc1_handle);

/* result is the complex value x[j][k], 0<=j<=31, 0<=k<=99 */

status = DftiCreateDescriptor( &my_desc2_handle, DFTI_SINGLE,

DFTI_REAL, 2, l);

status = DftiCommitDescriptor( my_desc2_handle);

status = DftiComputeForward( my_desc2_handle, y);

status = DftiFreeDescriptor(&my_desc2_handle);

/* result is the complex value z(j,k) 0<=j<=31; 0<=k<=99

/* and is stored in CCS format*/

## Changing Default Settings (Fortran)

The following examples demonstrate how you can change the default configuration settings by using the DftiSetValue function.

For instance, to preserve the input data after the FFT computation, the configuration of the DFTI_PLACEMENT should be changed to "not in place" from the default choice of "in place."

The code below illustrates how this can be done:

! Fortran example

! 1D complex to complex, not in place

Use MKL_DFTI

Complex :: X_in(32), X_out(32)

type(DFTI_DESCRIPTOR), POINTER :: My_Desc_Handle

Integer :: Status

!...put input data into X_in(j), 1<=j<=32

Status = DftiCreateDescriptor( My_Desc_Handle,& DFTI_SINGLE, DFTI_COMPLEX, 1, 32)

Status = DftiSetValue( My_Desc_Handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE)

Status = DftiCommitDescriptor( My_Desc_Handle)

Status = DftiComputeForward( My_Desc_Handle, X_in, X_out)

Status = DftiFreeDescriptor (My_Desc_Handle)

! result is X_out(1),X_out(2),...,X_out(32)

## Changing Default Settings (C)

/* C99 example */

#include "mkl_dfti.h"

float _Complex x_in[32], x_out[32];

DFTI_DESCRIPTOR_HANDLE my_desc_handle;

/* or alternatively

DFTI_DESCRIPTOR *my_desc_handle;*/

long status;

//...put input data into x_in[j], 0 <= j < 32

status = DftiCreateDescriptor( &my_desc_handle, DFTI_SINGLE,

DFTI_COMPLEX, 1, 32);

status = DftiSetValue( my_desc_handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE);

status = DftiCommitDescriptor( my_desc_handle);

status = DftiComputeForward( my_desc_handle, x_in, x_out);

status = DftiFreeDescriptor(&my_desc_handle);

/* result is x_out[0], x_out[1], ..., x_out[31] */

## Using Status Checking Function

The Example “Using Status Checking Function” illustrates the use of the status checking functions described in “Fourier Transform Functions”.

/* C */

DFTI_DESCRIPTOR_HANDLE desc;

long status, class_error, value;

char* error_message;

// . . . descriptor creation and other code

status = DftiGetValue( desc, DFTI_PRECISION, &value); //

class_error = DftiErrorClass(status, DFTI_NO_ERROR);

if (! class_error) {

printf ("DftiGetValue() fixes the wrong situation and

returns the corresponding value n");

error_message = DftiErrorMessage(status);

printf("error_message = %s \n", error_message);

}

// . . .

! Fortran

type(DFTI_DESCRIPTOR), POINTER :: desc

integer value, status

character(DFTI_MAX_MESSAGE_LENGTH) error_message

logical class_error

! . . . descriptor creation and other code

status = DftiGetValue( desc, DFTI_PRECISION, value)

class_error = DftiErrorClass(status, DFTI_NO_ERROR)

if (.not. class_error) then

print *, ' DftiGetValue() fixes the wrong situation and ',&

'returns the corresponding value '

error_message = DftiErrorMessage(status)

print *, 'error_message = ', error_message

endif

## Computing 2D FFT by One-Dimensional Transforms

Below is an example where a 20-by-40 two-dimensional FFT is computed explicitly using one-dimensional transforms. Notice that the data and result parameters in computation functions are all declared as assumed-size rank-1 array DIMENSION(0:*). Therefore two-dimensional array must be transformed to one-dimensional array by EQUIVALENCE statement or other facilities of Fortran.

! Fortran

use mkl_dfti

Complex :: X_2D(20,40)

Complex :: X(800)

Equivalence (X_2D, X)

INTEGER :: STRIDE(2)

type(DFTI_DESCRIPTOR), POINTER :: Desc_Handle_Dim1

type(DFTI_DESCRIPTOR), POINTER :: Desc_Handle_Dim2

! ...

Status = DftiCreateDescriptor(Desc_Handle_Dim1, DFTI_SINGLE,&

DFTI_COMPLEX, 1, 20 )

Status = DftiCreateDescriptor(Desc_Handle_Dim2, DFTI_SINGLE,&

DFTI_COMPLEX, 1, 40 )

! perform 40 one-dimensional transforms along 1st dimension

Status = DftiSetValue( Desc_Handle_Dim1, DFTI_NUMBER_OF_TRANSFORMS, 40 )

Status = DftiSetValue( Desc_Handle_Dim1, DFTI_INPUT_DISTANCE, 20 )

Status = DftiSetValue( Desc_Handle_Dim1, DFTI_OUTPUT_DISTANCE, 20 )

Status = DftiCommitDescriptor( Desc_Handle_Dim1 )

Status = DftiComputeForward( Desc_Handle_Dim1, X )

! perform 20 one-dimensional transforms along 2nd dimension

Stride(1) = 0; Stride(2) = 20

Status = DftiSetValue( Desc_Handle_Dim2, DFTI_NUMBER_OF_TRANSFORMS, 20 )

Status = DftiSetValue( Desc_Handle_Dim2, DFTI_INPUT_DISTANCE, 1 )

Status = DftiSetValue( Desc_Handle_Dim2, DFTI_OUTPUT_DISTANCE, 1 )

Status = DftiSetValue( Desc_Handle_Dim2, DFTI_INPUT_STRIDES, Stride )

Status = DftiSetValue( Desc_Handle_Dim2, DFTI_OUTPUT_STRIDES, Stride )

Status = DftiCommitDescriptor( Desc_Handle_Dim2 )

Status = DftiComputeForward( Desc_Handle_Dim2, X )

Status = DftiFreeDescriptor( Desc_Handle_Dim1 )

Status = DftiFreeDescriptor( Desc_Handle_Dim2 )

/* C */

#include "mkl_dfti.h"

float _Complex x[20][40];

long stride[2];

long status;

DFTI_DESCRIPTOR_HANDLE desc_handle_dim1;

DFTI_DESCRIPTOR_HANDLE desc_handle_dim2;

//...

status = DftiCreateDescriptor( &desc_handle_dim1, DFTI_SINGLE,

DFTI_COMPLEX, 1, 20 );

status = DftiCreateDescriptor( &desc_handle_dim2, DFTI_SINGLE,

DFTI_COMPLEX, 1, 40 );

/* perform 40 one-dimensional transforms along 1st dimension */

/* note that the 1st dimension data are not unit-stride */

stride[0] = 0; stride[1] = 40;

status = DftiSetValue( desc_handle_dim1, DFTI_NUMBER_OF_TRANSFORMS, 40 );

status = DftiSetValue( desc_handle_dim1, DFTI_INPUT_DISTANCE, 1 );

status = DftiSetValue( desc_handle_dim1, DFTI_OUTPUT_DISTANCE, 1 );

status = DftiSetValue( desc_handle_dim1, DFTI_INPUT_STRIDES, stride );

status = DftiSetValue( desc_handle_dim1, DFTI_OUTPUT_STRIDES, stride );

status = DftiCommitDescriptor( desc_handle_dim1 );

status = DftiComputeForward( desc_handle_dim1, x );

/* perform 20 one-dimensional transforms along 2nd dimension */

/* note that the 2nd dimension is unit stride */

status = DftiSetValue( desc_handle_dim2, DFTI_NUMBER_OF_TRANSFORMS, 20 );

status = DftiSetValue( desc_handle_dim2,

DFTI_INPUT_DISTANCE, 40 );

status = DftiSetValue( desc_handle_dim2,

DFTI_OUTPUT_DISTANCE, 40 );

status = DftiCommitDescriptor( desc_handle_dim2 );

status = DftiComputeForward( desc_handle_dim2, x );

status = DftiFreeDescriptor( &Desc_Handle_dim1 );

status = DftiFreeDescriptor( &Desc_Handle_dim2 );

## Two-Dimensional REAL In-place FFT (Fortran Interface)

The following are examples of real multi-dimensional transforms with CCE format storage of conjugate-even complex matrix. Example “Two-Dimensional REAL In-place FFT (Fortran Interface)” is two-dimensional in-place transform and Example “Two-Dimensional REAL Out-of-place FFT (Fortran Interface)” is two-dimensional out-of-place transform in Fortran interface. Example “Three-Dimensional REAL FFT (C Interface)” is three-dimensional out-of-place transform in C interface. Note that the data and result parameters in computation functions are all declared as assumed-size rank-1 array DIMENSION(0:*). Therefore two-dimensional array must be transformed to one-dimensional array by EQUIVALENCE statement or other facilities of Fortran.

`! Fortran example.`
`! 2D and real to conjugate even`
`Use MKL_DFTI`
`Real :: X_2D(34,100) ! 34 = (32/2 + 1)*2`
`Real :: X(3400)`
`Equivalence (X_2D, X)`
`type(DFTI_DESCRIPTOR), POINTER :: My_Desc_Handle`
`Integer :: Status, L(2)`
`Integer :: strides_in(3)`
`Integer :: strides_out(3)`
`! ...put input data into X_2D(j,k), 1<=j=32,1<=k<=100`
`! ...set L(1) = 32, L(2) = 100`
`! ...set strides_in(1) = 0, strides_in(2) = 1, strides_in(3) = 34`
`! ...set strides_out(1) = 0, strides_out(2) = 1, strides_out(3) = 17`
`! ...the transform is a 32-by-100`
`! Perform a real to complex conjugate even transform`
`Status = DftiCreateDescriptor( My_Desc_Handle, DFTI_SINGLE,&`
`DFTI_REAL, 2, L )`
`Status = DftiSetValue(My_Desc_Handle, DFTI_CONJUGATE_EVEN_STORAGE,&`
`DFTI_COMPLEX_COMPLEX)`
`Status = DftiSetValue(My_Desc_Handle, DFTI_INPUT_STRIDES, strides_in)`
`Status = DftiSetValue(My_Desc_Handle, DFTI_OUTPUT_STRIDES, strides_out)`
`Status = DftiCommitDescriptor( My_Desc_Handle)`
`Status = DftiComputeForward( My_Desc_Handle, X )`
`Status = DftiFreeDescriptor(My_Desc_Handle)`
`! result is given by the complex value z(j,k) 1<=j<=17; 1<=k<=100 and`
`! is stored in real matrix X_2D in CCE format.`

## Two-Dimensional REAL Out-of-place FFT (Fortran Interface)

! Fortran example.

! 2D and real to conjugate even

Use MKL_DFTI

Real :: X_2D(32,100)

Complex :: Y_2D(17, 100) ! 17 = 32/2 + 1

Real :: X(3200)

Complex :: Y(1700)

Equivalence (X_2D, X)

Equivalence (Y_2D, Y)

type(DFTI_DESCRIPTOR), POINTER :: My_Desc_Handle

Integer :: Status, L(2)

Integer :: strides_out(3)

! ...put input data into X_2D(j,k), 1<=j=32,1<=k<=100

! ...set L(1) = 32, L(2) = 100

! ...set strides_out(1) = 0, strides_out(2) = 1, strides_out(3) = 17

! ...the transform is a 32-by-100

! Perform a real to complex conjugate even transform

Status = DftiCreateDescriptor( My_Desc_Handle, DFTI_SINGLE,&

DFTI_REAL, 2, L )

Status = DftiSetValue(My_Desc_Handle,&

DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX)

Status = DftiSetValue( My_Desc_Handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE )

Status = DftiSetValue(My_Desc_Handle,&

DFTI_OUTPUT_STRIDES, strides_out)

Status = DftiCommitDescriptor(My_Desc_Handle)

Status = DftiComputeForward(My_Desc_Handle, X, Y)

Status = DftiFreeDescriptor(My_Desc_Handle)

! result is given by the complex value z(j,k) 1<=j<=17; 1<=k<=100 and

! is stored in complex matrix Y_2D in CCE format.

## Three-Dimensional REAL FFT (C Interface)

/* C99 example */

#include "mkl_dfti.h"

float x[32][100][19];

float _Complex y[32][100][10]; /* 10 = 19/2 + 1 */

DFTI_DESCRIPTOR_HANDLE my_desc_handle;

/* or alternatively

DFTI_DESCRIPTOR *my_desc_handle*/

long status, l[3];

long strides_out[4];

//...put input data into x[j][k][s] 0<=j<=31, 0<=k<=99, 0<=s<=18

l[0] = 32; l[1] = 100; l[2] = 19;

strides_out[0] = 0; strides_out[1] = 1000;

strides_out[2] = 10; strides_out[3] = 1;

status = DftiCreateDescriptor( &my_desc_handle, DFTI_SINGLE,

DFTI_REAL, 3, l );

status = DftiSetValue(my_desc_handle,

DFTI_CONJUGATE_EVEN_STORAGE, DFTI_COMPLEX_COMPLEX);

status = DftiSetValue( my_desc_handle, DFTI_PLACEMENT, DFTI_NOT_INPLACE );

status = DftiSetValue(my_desc_handle,

DFTI_OUTPUT_STRIDES, strides_out);

status = DftiCommitDescriptor(my_desc_handle);

status = DftiComputeForward(my_desc_handle, x, y);

status = DftiFreeDescriptor(&my_desc_handle);

/* result is the complex value z(j,k,s) 0<=j<=31; 0<=k<=99, 0<=s<=9

and is stored in complex matrix y in CCE format. */

- Examples of Using Multi-Threading for FFT Computation

Parent topic: Fourier Transform Functions Code Examples