sll_fft

Interface around fftpack, fftw and the selalib fft.

Author

Edwin Chacon-Golcher, Samuel De Santis, Pierre Navaro, Katharina Kormann

How to use sll_m_fft module?

The first thing is to add the line

use sll_m_fft 

The sll_m_fft module can use internal or external librarie s

1. Declare a fft plan

   type(sll_t_fft) :: p 

2. Initialize the plan

   call sll_s_fft_init_c2c_1d( p, size, in, out, direction, normalized, aligned, optimization ) 

The arrays in and out can be real and/or complex, 1d or 2d. Change c2c_1d accordingly.

Warning

For complex to real and real to complex transform, there is no direction flag.

call sll_s_fft_init_r2c_1d( p, size, in, out, normalized, aligned, optimization ) 

direction can take two values : sll_p_fft_forward and sll_p_fft_backward

normalized is a boolean optional argument: If true, the transformed data is normalized by the (product) of the element length. [default: false] aligned is a boolean optional argument (only used by FFTW): If true, FFTW assumes the the data is aligned (use fft_alloc for in/out data allocation in this case). [default: false] optimization is an optional argument (only used by FFTW): With this argument, you can set how rigorous the initialization of the planner will be. It can take the values sll_p_fft_patient, sll_p_fft_estimate, sll_p_fft_exhaustive, sll_p_fft_wisdom_only. Note that in and out are overwritten for the last three options. sll_p_fft_wisdom_only only works if wisdom is available. [default: sll_p_fft_estimate]

1. Execute the plan

   call sll_s_fft_exec_c2c_1d( p, in, out ) 

2. Delete the plan

   call sll_s_fft_free( p ) 

Summary:

1D

size problem	type of in	type of out	size of in	size of out	direction	normalized	aligned	optimization
n	real	real	n	n	sll_p_fft_forward sll_p_fft_backward	.TRUE. .FALSE.	.TRUE. .FALSE.	sll_p_fft_estimate sll_p_wisdom_only sll_p_fft_measure sll_p_fft_patient sll_p_fft_exhaustive
n	complex	complex	n	n	sll_p_fft_forward sll_p_fft_backward	.TRUE. .FALSE.	.TRUE. .FALSE.	sll_p_fft_estimate sll_p_wisdom_only sll_p_fft_measure sll_p_fft_patient sll_p_fft_exhaustive
n	real	complex	n	n/2+1	--—	.TRUE. .FALSE.	.TRUE. .FALSE.	sll_p_fft_estimate sll_p_wisdom_only sll_p_fft_measure sll_p_fft_patient sll_p_fft_exhaustive
n	complex	real	n/2+1	n	--—	.TRUE. .FALSE.	.TRUE. .FALSE.	sll_p_fft_estimate sll_p_wisdom_only sll_p_fft_measure sll_p_fft_patient sll_p_fft_exhaustive

2D

size problem	type of in	type of out	size of in	size of out	direction present	normalized	aligned	optimization
n,m	real	real	n,m	n,m	sll_p_fft_forward sll_p_fft_backward	.TRUE. .FALSE.	.TRUE. .FALSE.	sll_p_fft_estimate sll_p_wisdom_only sll_p_fft_measure sll_p_fft_patient sll_p_fft_exhaustive
n,m	complex	complex	n,m	n,m	sll_p_fft_forward sll_p_fft_backward	.TRUE. .FALSE.	.TRUE. .FALSE.	sll_p_fft_estimate sll_p_wisdom_only sll_p_fft_measure sll_p_fft_patient sll_p_fft_exhaustive
n,m	real	complex	n,m	n/2,m	--—	.TRUE. .FALSE.	.TRUE. .FALSE.	sll_p_fft_estimate sll_p_wisdom_only sll_p_fft_measure sll_p_fft_patient sll_p_fft_exhaustive
n,m	complex	real	n/2,m	n,m	--—	.TRUE. .FALSE.	.TRUE. .FALSE.	sll_p_fft_estimate sll_p_wisdom_only sll_p_fft_measure sll_p_fft_patient sll_p_fft_exhaustive

Examples:

In-place transform

integer(kind=i32), parameter :: n = 2**5
complex(kind=f64), dimension(0,n-1) :: in
type(sll_t_fft)  :: p
 
!** INIT DATA **
 
call sll_s_fft_init_c2c_1d( p, n, in, in, sll_p_fft_forward, normalized = .true. )
call sll_s_fft_exec_c2c_1d( p, in, in )
call sll_s_fft_free( p )

Two-dimensional transform

integer(kind=i32), parameter :: n = 2**5
integer(kind=i32), parameter :: m = 2**3
complex(kind=f64), dimension(n/2,m) :: in
real(kind=f64), dimension(n,m) :: out
type(sll_t_fft)  :: p
 
!** INIT DATA **
 
call sll_s_fft_init_c2r_2d( p, n, m, in, out )
call sll_s_fft_exec_c2r_2d( p, in, out )
call sll_s_fft_free( p )

What sll_m_fft really computes

The forward (sll_p_fft_forward) DFT of a 1d complex array x of size n computes an array X, where:

\[ X_k = \sum_{i=0}^{n-1} x_i e^{-2\pi i j k/n}. \]

The backward (sll_p_fft_backward) DFT computes:

\[ x_i = \sum_{k=0}^{n-1} X_k e^{2\pi k j i/n}. \]

For the real transform, we have \( (x_0,x_1,\dots,x_{n-1}) \rightarrow (r_0,r_1,\dots,r_{n/2-1},r_{n/2},i_{n/2-1},\dots,i_1)\) which must be interpreted as the complex array

\[ \begin{pmatrix} r_0 &,& 0 \\ r_1 &,& i_1 \\ \vdots & & \vdots \\ r_{n/2-1} &,& i_{n/2-1} \\ r_{n/2} &,& 0 \\ r_{n/2-1} &,& -i_{n/2-1} \\ \vdots & & \vdots \\ r_1 &,& -i_1 \end{pmatrix}\]

Warning

Note that ffw uses \((r_0,r_1,\dots,r_{n/2-1},r_{n/2},i_{n/2-1},\dots,i_1)\) convention whereas fftpack uses \((r_0,r_1,i_1,\dots,r_{n/2-1},i_{n/2-1},r_{n/2})\). Through the interface, FFTW ordering is enforced also for FFTPACK. For r2c and r2c, the lower half (plus one element) is stored (non-negative frequencies). In higher dimensions, the first dimension is reduced to n/2+1 whereas the others remain n.