Intel® oneAPI Math Kernel Library Developer Reference - Fortran
This section presents 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”.
In the Fortran examples, the use mkl_dfti statement assumes that:
! 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.
! 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.
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
The following example demonstrates 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 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)
This example illustrates the use of status checking functions described in “Fourier Transform Functions”.
! Fortran
type(DFTI_DESCRIPTOR), POINTER :: desc
integer status
! ...descriptor creation and other code
status = DftiCommitDescriptor(desc)
if (status .ne. 0) then
if (.not. DftiErrorClass(status,DFTI_NO_ERROR) then
print *, 'Error: ‘, DftiErrorMessage(status)
endif
endif
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 )
The following code illustrates 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. 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.
! 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.