sll_m_fft_fftw3.F90

Go to the documentation of this file.

!**************************************************************
!  Copyright INRIA
!  Authors :
!     CALVI project team
!
!  This code SeLaLib (for Semi-Lagrangian-Library)
!  is a parallel library for simulating the plasma turbulence
!  in a tokamak.
!
!  This software is governed by the CeCILL-B license
!  under French law and abiding by the rules of distribution
!  of free software.  You can  use, modify and redistribute
!  the software under the terms of the CeCILL-B license as
!  circulated by CEA, CNRS and INRIA at the following URL
!  "http://www.cecill.info".
!**************************************************************
 
module sll_m_fft
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#include "sll_working_precision.h"
#include "sll_assert.h"
#include "sll_memory.h"
#include "sll_errors.h"
#include "sll_fftw.h"
 
   use sll_m_fftw3
   use, intrinsic :: iso_c_binding, only: &
      c_f_pointer, &
      c_loc
 
   implicit none
 
   public :: &
      sll_t_fft, &
      sll_p_fft_forward, &
      sll_p_fft_backward, &
      sll_p_fft_measure, &
      sll_p_fft_patient, &
      sll_p_fft_estimate, &
      sll_p_fft_exhaustive, &
      sll_p_fft_wisdom_only, &
      sll_s_fft_print_defaultlib, &
      sll_f_fft_allocate_aligned_complex, &
      sll_f_fft_allocate_aligned_real, &
      sll_s_fft_deallocate_aligned_complex, &
      sll_s_fft_deallocate_aligned_real, &
      sll_s_fft_init_r2r_1d, &
      sll_s_fft_init_c2r_1d, &
      sll_s_fft_init_r2c_1d, &
      sll_s_fft_init_c2c_1d, &
      sll_s_fft_init_r2c_2d, &
      sll_s_fft_init_c2r_2d, &
      sll_s_fft_init_c2c_2d, &
      sll_s_fft_init_c2c_3d, &
      sll_s_fft_exec_r2r_1d, &
      sll_s_fft_exec_c2r_1d, &
      sll_s_fft_exec_r2c_1d, &
      sll_s_fft_exec_c2c_1d, &
      sll_s_fft_exec_r2c_2d, &
      sll_s_fft_exec_c2r_2d, &
      sll_s_fft_exec_c2c_2d, &
      sll_s_fft_exec_c2c_3d, &
      sll_s_fft_set_mode_c2r_1d, &
      sll_f_fft_get_mode_r2c_1d, &
      sll_s_fft_get_k_list_c2c_1d, &
      sll_s_fft_get_k_list_r2r_1d, &
      sll_s_fft_get_k_list_r2c_1d, &
      sll_s_fft_free
 
   private
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   type sll_t_fft
 
      fftw_plan, private               :: fftw       
      logical                          :: normalized
      sll_int32                        :: direction  
      sll_int32                        :: problem_rank 
      sll_int32, allocatable           :: problem_shape(:) 
 
      sll_comp64, allocatable, private :: scratch(:) 
      sll_int32, private               :: transform_type 
 
   end type sll_t_fft
 
   ! Flags for direction (values as in fftw3.f03)
   integer, parameter :: sll_p_fft_forward = fftw_forward 
   integer, parameter :: sll_p_fft_backward = fftw_backward 
 
   ! Flags for initialization of the plan (values as in fftw3.f03)
   integer(C_INT), parameter :: sll_p_fft_measure = fftw_measure 
   integer(C_INT), parameter :: sll_p_fft_patient = fftw_patient 
   integer(C_INT), parameter :: sll_p_fft_estimate = fftw_estimate 
   integer(C_INT), parameter :: sll_p_fft_exhaustive = fftw_exhaustive 
   integer(C_INT), parameter :: sll_p_fft_wisdom_only = fftw_wisdom_only ! 2097152 !< FFTW planning-rigor flag FFTW_WISDOM_ONLY (planner only initialized if wisdom is available)  [value 2097152]
 
! Flags to pass when we create a new plan
! We can define 31 different flags.
! The value assigned to the flag can only be a power of two.
! See section "How-to manipulate flags ?" for more information.
 
   ! Flags for the various types of transform (to make sure same type of init and execute functions are used)
   integer, parameter :: p_fftw_c2c_1d = 0
   integer, parameter :: p_fftw_r2r_1d = 1
   integer, parameter :: p_fftw_r2c_1d = 2
   integer, parameter :: p_fftw_c2r_1d = 3
   integer, parameter :: p_fftw_c2c_2d = 4
   integer, parameter :: p_fftw_r2r_2d = 5
   integer, parameter :: p_fftw_r2c_2d = 6
   integer, parameter :: p_fftw_c2r_2d = 7
   integer, parameter :: p_fftw_c2c_3d = 8
 
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
contains
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_print_defaultlib()
      print *, 'The library used is FFTW'
   end subroutine
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   function sll_f_fft_allocate_aligned_complex(n) result(data)
      sll_int32, intent(in)  :: n       
      complex(C_DOUBLE_COMPLEX), pointer     :: data(:)
 
      type(c_ptr) :: ptr_data         ! C pointer needed for allocation
 
#ifdef FFTW_F2003
      ptr_data = fftw_alloc_complex(int(n, kind=c_size_t))
      call c_f_pointer(ptr_data, data, [n])
#else
      allocate (data(n))
      sll_warning('sll_f_fft_allocate_aligned_complex', 'Aligned allocation only supported by F2003. Usual allocation.')
#endif
 
   end function sll_f_fft_allocate_aligned_complex
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   function sll_f_fft_allocate_aligned_real(n) result(data)
      sll_int32, intent(in) :: n       
      real(c_double), pointer    :: data(:)
 
      type(c_ptr) :: ptr_data         ! C pointer needed for allocation
 
#ifdef FFTW_F2003
      ptr_data = fftw_alloc_real(int(n, kind=c_size_t))
      call c_f_pointer(ptr_data, data, [n])
#else
      allocate (data(n))
      sll_warning('sll_f_fft_allocate_aligned_real', 'Aligned allocation only supported by F2003. Usual allocation.')
#endif
 
   end function sll_f_fft_allocate_aligned_real
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_deallocate_aligned_complex(data)
      complex(C_DOUBLE_COMPLEX), pointer :: data(:)
 
#ifdef FFTW_F2003
      type(c_ptr) :: ptr_data ! C pointer needed for deallocation
      ptr_data = c_loc(data(1))
      call fftw_free(ptr_data)
#else
      deallocate (data)
      sll_warning('sll_s_fft_deallocate_aligned_complex', 'Aligned deallocation only supported by F2003. Usual deallocation.')
#endif
 
   end subroutine sll_s_fft_deallocate_aligned_complex
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_deallocate_aligned_real(data)
      real(c_double), pointer :: data(:)
 
#ifdef FFTW_F2003
      type(c_ptr) :: ptr_data ! C pointer needed for deallocation
      ptr_data = c_loc(data(1))
      call fftw_free(ptr_data)
#else
      deallocate (data)
      sll_warning('sll_s_fft_deallocate_aligned_real', 'Aligned deallocation only supported by F2003. Usual deallocation.')
#endif
 
   end subroutine sll_s_fft_deallocate_aligned_real
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_get_k_list_c2c_1d(plan, k_list)
      type(sll_t_fft), intent(in) :: plan
      sll_int32, intent(out) :: k_list(0:) 
 
      sll_int32 :: k  ! single mode
      sll_int32 :: n  ! problem size
      n = plan%problem_shape(1)
 
      if (size(k_list) /= n) then
         sll_error("sll_s_fft_get_k_list_c2c_1d", "k_list must have n elements")
      end if
 
      k_list(0:n/2) = [(k, k=0, n/2)]
      k_list(n/2 + 1:n - 1) = [(k - n, k=n/2 + 1, n - 1)]
 
   end subroutine sll_s_fft_get_k_list_c2c_1d
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_get_k_list_r2r_1d(plan, k_list)
      type(sll_t_fft), intent(in) :: plan
      sll_int32, intent(out) :: k_list(0:) 
 
      sll_int32 :: k  ! single mode
      sll_int32 :: n  ! problem size
      n = plan%problem_shape(1)
 
      if (size(k_list) /= n) then
         sll_error("sll_s_fft_get_k_list_r2r_1d", "k_list must have n elements")
      end if
 
      k_list(0:n/2) = [(k, k=0, n/2)]
      k_list(n/2 + 1:n - 1) = [(n - k, k=n/2 + 1, n - 1)]
 
   end subroutine sll_s_fft_get_k_list_r2r_1d
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_get_k_list_r2c_1d(plan, k_list)
      type(sll_t_fft), intent(in) :: plan
      sll_int32, intent(out) :: k_list(0:) 
 
      sll_int32 :: k    ! single mode
      sll_int32 :: n_2  ! problem size / 2 (rounded down)
      n_2 = plan%problem_shape(1)/2
 
      if (size(k_list) /= n_2 + 1) then
         sll_error("sll_s_fft_get_k_list_r2c_1d", "k_list must have n/2+1 elements")
      end if
 
      k_list(:) = [(k, k=0, n_2)]
 
   end subroutine sll_s_fft_get_k_list_r2c_1d
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   function sll_f_fft_get_mode_r2c_1d(plan, data, k) result(ck)
      type(sll_t_fft), intent(in) :: plan
      sll_real64, intent(in) :: data(0:) 
      sll_int32, intent(in) :: k        
      sll_comp64                  :: ck       
 
      sll_int32 :: n  ! problem size
      n = plan%problem_shape(1)
 
      if (k == 0) then; ck = cmplx(data(k), 0.0_f64, kind=f64)
      else if (k <= (n + 1)/2 - 1) then; ck = cmplx(data(k), data(n - k), kind=f64)
      else if (k == n/2) then; ck = cmplx(data(k), 0.0_f64, kind=f64)
      else if (k <= n - 1) then; ck = cmplx(data(n - k), -data(k), kind=f64)
      else
         sll_error("sll_f_fft_get_mode_r2c_1d", "k must be between 0 and n-1")
      end if
 
   end function sll_f_fft_get_mode_r2c_1d
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_set_mode_c2r_1d(plan, data, ck, k)
      type(sll_t_fft), intent(in) :: plan
      sll_real64, intent(inout) :: data(0:)  
      sll_comp64, intent(in) :: ck        
      sll_int32, intent(in) :: k         
 
      sll_int32 :: n  ! problem size
      n = plan%problem_shape(1)
 
      if (k == 0) then; data(0) = real(ck) ! imaginary part ignored
      else if (k <= (n + 1)/2 - 1) then; data(k) = real(ck); data(n - k) = aimag(ck)
      else if (k == n/2) then; data(k) = real(ck) ! imaginary part ignored
      else if (k <= n - 1) then; data(n - k) = real(ck); data(k) = -aimag(ck)
      else
         sll_error("sll_s_fft_set_mode_c2r_1d", "k must be between 0 and n-1")
      end if
 
   end subroutine sll_s_fft_set_mode_c2r_1d
 
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!                             COMPLEX to COMPLEX (1D)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_init_c2c_1d(plan, nx, array_in, array_out, direction, normalized, aligned, optimization)
      type(sll_t_fft)                         :: plan
      sll_int32, intent(in)         :: nx 
      sll_comp64, intent(inout)      :: array_in(:) 
      sll_comp64, intent(inout)      :: array_out(:) 
      sll_int32, intent(in)         :: direction  
      logical, optional, intent(in)         :: normalized
      logical, optional, intent(in)         :: aligned
      sll_int32, optional, intent(in)         :: optimization
 
      ! local variables
      sll_int32                                :: ierr
      sll_int32                                :: flag_fftw
 
      plan%transform_type = p_fftw_c2c_1d
      plan%direction = direction
      if (present(normalized)) then
         plan%normalized = normalized
      else
         plan%normalized = .false.
      end if
      ! Set the information about the algorithm to compute the plan. The default is FFTW_ESTIMATE
      if (present(optimization)) then
         flag_fftw = optimization
      else
         flag_fftw = fftw_estimate
      end if
      if (present(aligned)) then
         if (aligned .EQV. .false.) then
            flag_fftw = flag_fftw + fftw_unaligned
         end if
      else
         flag_fftw = flag_fftw + fftw_unaligned
      end if
      plan%problem_rank = 1
      sll_allocate(plan%problem_shape(1), ierr)
      plan%problem_shape = (/nx/)
 
#ifdef FFTW_F2003
      plan%fftw = fftw_plan_dft_1d(nx, array_in, array_out, direction, flag_fftw)
#else
      call dfftw_plan_dft_1d(plan%fftw, nx, array_in, array_out, direction, flag_fftw)
#endif
 
   end subroutine sll_s_fft_init_c2c_1d
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_exec_c2c_1d(plan, array_in, array_out)
      type(sll_t_fft), intent(in)            :: plan
      sll_comp64, intent(inout)         :: array_in(:) 
      sll_comp64, intent(inout)         :: array_out(:) 
 
      sll_real64 :: factor
 
      sll_assert(plan%transform_type == p_fftw_c2c_1d)
 
      call fftw_execute_dft(plan%fftw, array_in, array_out)
 
      if (plan%normalized .EQV. .true.) then
         factor = 1.0_f64/real(plan%problem_shape(1), kind=f64)
         array_out = cmplx(factor, 0.0, f64)*array_out
      end if
   end subroutine
 
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!                             COMPLEX to COMPLEX (2D)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_init_c2c_2d(plan, nx, ny, array_in, array_out, direction, normalized, aligned, optimization)
      type(sll_t_fft), intent(out)    :: plan
      sll_int32, intent(in)     :: nx 
      sll_int32, intent(in)     :: ny 
      sll_comp64, intent(inout)  :: array_in(0:, 0:) 
      sll_comp64, intent(inout)  :: array_out(0:, 0:) 
      sll_int32, intent(in)     :: direction  
      logical, optional, intent(in)     :: normalized
      logical, optional, intent(in)     :: aligned
      sll_int32, optional, intent(in)      :: optimization
      
 
      sll_int32                                     :: ierr
      sll_int32                                     :: flag_fftw
 
      plan%transform_type = p_fftw_c2c_2d
      plan%direction = direction
      if (present(normalized)) then
         plan%normalized = normalized
      else
         plan%normalized = .false.
      end if
      ! Set the information about the algorithm to compute the plan. The default is FFTW_ESTIMATE
      if (present(optimization)) then
         flag_fftw = optimization
      else
         flag_fftw = fftw_estimate
      end if
      if (present(aligned)) then
         if (aligned .EQV. .false.) then
            flag_fftw = flag_fftw + fftw_unaligned
         end if
      else
         flag_fftw = flag_fftw + fftw_unaligned
      end if
      plan%problem_rank = 2
      sll_allocate(plan%problem_shape(2), ierr)
      plan%problem_shape = (/nx, ny/)
 
      !We must switch the dimension. It's a fftw convention.
 
#ifdef FFTW_F2003
      plan%fftw = fftw_plan_dft_2d(ny, nx, array_in, array_out, direction, flag_fftw)!FFTW_ESTIMATE + FFTW_UNALIGNED)
#else
      call dfftw_plan_dft_2d(plan%fftw, ny, nx, array_in, array_out, direction, flag_fftw)!FFTW_ESTIMATE + FFTW_UNALIGNED)
#endif
 
   end subroutine sll_s_fft_init_c2c_2d
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_exec_c2c_2d(plan, array_in, array_out)
      type(sll_t_fft), intent(in)     :: plan
      sll_comp64, intent(inout)  :: array_in(0:, 0:) 
      sll_comp64, intent(inout)  :: array_out(0:, 0:) 
 
      sll_real64                                   :: factor
 
      call fftw_execute_dft(plan%fftw, array_in, array_out)
 
      if (plan%normalized .EQV. .true.) then
         factor = 1.0_f64/real(plan%problem_shape(1)*plan%problem_shape(2), kind=f64)
         array_out = cmplx(factor, 0.0, f64)*array_out
      end if
 
   end subroutine
 
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!                             COMPLEX to COMPLEX (3D)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_init_c2c_3d(plan, nx, ny, nz, array_in, array_out, direction, normalized, aligned, optimization)
      type(sll_t_fft), intent(out)    :: plan
      sll_int32, intent(in)     :: nx 
      sll_int32, intent(in)     :: ny 
      sll_int32, intent(in)     :: nz 
      sll_comp64, intent(inout)  :: array_in(0:, 0:, 0:) 
      sll_comp64, intent(inout)  :: array_out(0:, 0:, 0:) 
      sll_int32, intent(in)     :: direction  
      logical, optional, intent(in)     :: normalized
      logical, optional, intent(in)     :: aligned
      sll_int32, optional, intent(in)      :: optimization
      
 
      sll_int32                                     :: ierr
      sll_int32                                     :: flag_fftw
 
      plan%transform_type = p_fftw_c2c_3d
      plan%direction = direction
      if (present(normalized)) then
         plan%normalized = normalized
      else
         plan%normalized = .false.
      end if
      ! Set the information about the algorithm to compute the plan. The default is FFTW_ESTIMATE
      if (present(optimization)) then
         flag_fftw = optimization
      else
         flag_fftw = fftw_estimate
      end if
      if (present(aligned)) then
         if (aligned .EQV. .false.) then
            flag_fftw = flag_fftw + fftw_unaligned
         end if
      else
         flag_fftw = flag_fftw + fftw_unaligned
      end if
      plan%problem_rank = 3
      sll_allocate(plan%problem_shape(3), ierr)
      plan%problem_shape = [nx, ny, nz]
 
      !We must switch the dimension. It's a fftw convention.
 
#ifdef FFTW_F2003
      plan%fftw = fftw_plan_dft_3d(nz, ny, nx, array_in, array_out, direction, flag_fftw)!FFTW_ESTIMATE + FFTW_UNALIGNED)
#else
      call dfftw_plan_dft_3d(plan%fftw, nz, ny, nx, array_in, array_out, direction, flag_fftw)!FFTW_ESTIMATE + FFTW_UNALIGNED)
#endif
 
   end subroutine sll_s_fft_init_c2c_3d
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_exec_c2c_3d(plan, array_in, array_out)
      type(sll_t_fft), intent(in)     :: plan
      sll_comp64, intent(inout)  :: array_in(0:, 0:, 0:) 
      sll_comp64, intent(inout)  :: array_out(0:, 0:, 0:) 
 
      sll_real64                                   :: factor
 
      call fftw_execute_dft(plan%fftw, array_in, array_out)
 
      if (plan%normalized) then
         factor = 1.0_f64/real(plan%problem_shape(1)*plan%problem_shape(2)*plan%problem_shape(3), kind=f64)
         array_out = cmplx(factor, 0.0, f64)*array_out
      end if
 
   end subroutine
 
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!                                 REAL to REAL (1D)
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_init_r2r_1d(plan, nx, array_in, array_out, direction, normalized, aligned, optimization)
      type(sll_t_fft), intent(out)   :: plan
      sll_int32, intent(in)    :: nx 
      sll_real64, intent(inout) :: array_in(:) 
      sll_real64, intent(inout) :: array_out(:) 
      sll_int32, intent(in)    :: direction  
      logical, optional, intent(in)    :: normalized
      logical, optional, intent(in)    :: aligned
      sll_int32, optional, intent(in)    :: optimization
      
 
      sll_int32 :: ierr
      sll_int32 :: flag_fftw
 
      plan%transform_type = p_fftw_r2r_1d
      plan%direction = direction
      if (present(normalized)) then
         plan%normalized = normalized
      else
         plan%normalized = .false.
      end if
      ! Set the information about the algorithm to compute the plan. The default is FFTW_ESTIMATE
      if (present(optimization)) then
         flag_fftw = optimization
      else
         flag_fftw = fftw_estimate
      end if
      if (present(aligned)) then
         if (aligned .EQV. .false.) then
            flag_fftw = flag_fftw + fftw_unaligned
         end if
      else
         flag_fftw = flag_fftw + fftw_unaligned
      end if
      plan%problem_rank = 1
      sll_allocate(plan%problem_shape(1), ierr)
      plan%problem_shape = (/nx/)
 
      if (direction .eq. sll_p_fft_forward) then
#ifdef FFTW_F2003
         plan%fftw = fftw_plan_r2r_1d(nx, array_in, array_out, fftw_r2hc, flag_fftw)
#else
         ! This would be the call to actually create a r2r plan
         !call dfftw_plan_r2r_1d(plan%fftw,nx,array_in,array_out,FFTW_R2HC,flag_fftw)
 
         allocate (plan%scratch(nx/2 + 1))
         call dfftw_plan_dft_r2c_1d(plan%fftw, nx, array_in, plan%scratch, flag_fftw)
       sll_warning('sll_s_fft_init_r2r_1d', 'R2HC not supported without FFTW_F2003. Use implementation based on r2c/c2r transform.')
#endif
      else if (direction .eq. sll_p_fft_backward) then
#ifdef FFTW_F2003
         plan%fftw = fftw_plan_r2r_1d(nx, array_in, array_out, fftw_hc2r, flag_fftw)
#else
         ! This would be the code to actually create an r2r plan
         !call dfftw_plan_r2r_1d(plan%fftw,nx,array_in,array_out,FFTW_HC2R,flag_fftw)
 
         allocate (plan%scratch(nx/2 + 1))
         call dfftw_plan_dft_c2r_1d(plan%fftw, nx, plan%scratch, array_out, flag_fftw)
       sll_warning('sll_s_fft_init_r2r_1d', 'HC2R not supported without FFTW_F2003. Use implementation based on r2c/c2r transform.')
#endif
      end if
   end subroutine sll_s_fft_init_r2r_1d
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_exec_r2r_1d(plan, array_in, array_out)
      type(sll_t_fft), intent(inout) :: plan
      sll_real64, intent(inout) :: array_in(:) 
      sll_real64, intent(inout) :: array_out(:) 
 
      sll_real64                              :: factor
#ifndef FFTW_F2003
      sll_int32                               :: j
#endif
 
      sll_assert(plan%transform_type == p_fftw_r2r_1d)
 
#ifdef FFTW_F2003
      call fftw_execute_r2r(plan%fftw, array_in, array_out)
 
      if (plan%normalized .EQV. .true.) then
         factor = 1.0_f64/real(plan%problem_shape(1), kind=f64)
         array_out = factor*array_out
      end if
 
#else
 
      if (plan%direction == sll_p_fft_forward) then
         call sll_s_fft_exec_r2c_1d(plan, array_in, plan%scratch)
         array_out(1) = real(plan%scratch(1), f64)
         array_out(plan%problem_shape(1)/2 + 1) = real(plan%scratch((plan%problem_shape(1)/2 + 1)), f64)
         do j = 1, plan%problem_shape(1)/2 - 1
            call sll_s_fft_set_mode_c2r_1d(plan, array_out, plan%scratch(j + 1), j)
         end do
      elseif (plan%direction == sll_p_fft_backward) then
         do j = 0, plan%problem_shape(1)/2
            plan%scratch(j + 1) = sll_f_fft_get_mode_r2c_1d(plan, array_out, j)
         end do
         call sll_s_fft_exec_c2r_1d(plan, plan%scratch, array_out)
      end if
 
#endif
 
   end subroutine
 
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!                               REAL to COMPLEX
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_init_r2c_1d(plan, nx, array_in, array_out, normalized, aligned, optimization)
      type(sll_t_fft), intent(out)   :: plan
      sll_int32, intent(in)    :: nx 
      sll_real64, intent(inout) :: array_in(:) 
      sll_comp64, intent(out)   :: array_out(:) 
      logical, optional, intent(in)    :: normalized
      logical, optional, intent(in)    :: aligned
      sll_int32, optional, intent(in)    :: optimization
      
 
      sll_int32 :: ierr
      sll_int32 :: flag_fftw
 
      plan%transform_type = p_fftw_r2c_1d
      plan%direction = 0
      if (present(normalized)) then
         plan%normalized = normalized
      else
         plan%normalized = .false.
      end if
      ! Set the information about the algorithm to compute the plan. The default is FFTW_ESTIMATE
      if (present(optimization)) then
         flag_fftw = optimization
      else
         flag_fftw = fftw_estimate
      end if
      if (present(aligned)) then
         if (aligned .EQV. .false.) then
            flag_fftw = flag_fftw + fftw_unaligned
         end if
      else
         flag_fftw = flag_fftw + fftw_unaligned
      end if
      plan%problem_rank = 1
      sll_allocate(plan%problem_shape(1), ierr)
      plan%problem_shape = (/nx/)
 
#ifdef FFTW_F2003
      plan%fftw = fftw_plan_dft_r2c_1d(nx, array_in, array_out, flag_fftw)
#else
      call dfftw_plan_dft_r2c_1d(plan%fftw, nx, array_in, array_out, flag_fftw)
#endif
   end subroutine sll_s_fft_init_r2c_1d
 
   subroutine sll_s_fft_exec_r2c_1d(plan, array_in, array_out)
      type(sll_t_fft), intent(in)            :: plan
      sll_real64, intent(inout)         :: array_in(:) 
      sll_comp64, intent(out)           :: array_out(:) 
 
      sll_real64                                      :: factor
 
      sll_assert(plan%transform_type == p_fftw_r2c_1d)
 
      call fftw_execute_dft_r2c(plan%fftw, array_in, array_out)
 
      if (plan%normalized .EQV. .true.) then
         factor = 1.0_f64/real(plan%problem_shape(1), kind=f64)
         array_out = cmplx(factor, 0.0, kind=f64)*array_out
      end if
   end subroutine
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_init_r2c_2d(plan, nx, ny, array_in, array_out, normalized, aligned, optimization)
      type(sll_t_fft), intent(out)   :: plan
      sll_int32, intent(in)    :: nx 
      sll_int32, intent(in)    :: ny 
      sll_real64, intent(inout) :: array_in(:, :) 
      sll_comp64, intent(out)   :: array_out(:, :) 
      logical, optional, intent(in)    :: normalized
      logical, optional, intent(in)    :: aligned
      sll_int32, optional, intent(in)    :: optimization
      
 
      sll_int32 :: ierr
      sll_int32 :: flag_fftw
 
      plan%transform_type = p_fftw_r2c_2d
      plan%direction = 0
      if (present(normalized)) then
         plan%normalized = normalized
      else
         plan%normalized = .false.
      end if
      ! Set the information about the algorithm to compute the plan. The default is FFTW_ESTIMATE
      if (present(optimization)) then
         flag_fftw = optimization
      else
         flag_fftw = fftw_estimate
      end if
      if (present(aligned)) then
         if (aligned .EQV. .false.) then
            flag_fftw = flag_fftw + fftw_unaligned
         end if
      else
         flag_fftw = flag_fftw + fftw_unaligned
      end if
 
      plan%problem_rank = 2
      sll_allocate(plan%problem_shape(2), ierr)
      plan%problem_shape = (/nx, ny/)
 
#ifdef FFTW_F2003
      plan%fftw = fftw_plan_dft_r2c_2d(ny, nx, array_in, array_out, fftw_estimate + fftw_unaligned)
#else
      call dfftw_plan_dft_r2c_2d(plan%fftw, nx, ny, array_in, array_out, fftw_estimate)
#endif
 
   end subroutine sll_s_fft_init_r2c_2d
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_exec_r2c_2d(plan, array_in, array_out)
      type(sll_t_fft), intent(in)           :: plan
      sll_real64, intent(inout)        :: array_in(:, :)  
      sll_comp64, intent(out)          :: array_out(:, :) 
 
      sll_int32     :: nx, ny
      sll_real64    :: factor
 
      sll_assert(plan%transform_type == p_fftw_r2c_2d)
 
      nx = plan%problem_shape(1)
      ny = plan%problem_shape(2)
 
      if ((size(array_in, dim=1) .ne. nx) .and. (size(array_in, dim=2) .ne. ny)) then
         print *, 'Error in file sll_m_fft.F90'
         print *, '      in subroutine fft_apply_r2c_2d'
         print *, '      array_in size problem'
         stop ''
      else if (size(array_in, dim=1) .ne. nx/2 + 1 .and. size(array_in, dim=2) .ne. ny) then
         print *, 'Error in file sll_m_fft.F90'
         print *, '      in subroutine fft_apply_r2c_2d'
         print *, '      array_out size problem'
         stop ''
      end if
 
      call fftw_execute_dft_r2c(plan%fftw, array_in, array_out)
 
      if (plan%normalized .EQV. .true.) then
         factor = 1.0_f64/real(nx*ny, kind=f64)
         array_out = cmplx(factor, 0.0, kind=f64)*array_out
      end if
   end subroutine
 
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!                               COMPLEX to REAL
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_init_c2r_1d(plan, nx, array_in, array_out, normalized, aligned, optimization)
      type(sll_t_fft), intent(out)      :: plan
      sll_int32, intent(in)       :: nx 
      sll_comp64, intent(inout)    :: array_in(:)  
      sll_real64, intent(out)      :: array_out(:) 
      logical, optional, intent(in)       :: normalized
      logical, optional, intent(in)       :: aligned
      sll_int32, optional, intent(in)       :: optimization
      
 
      sll_int32 :: ierr
      sll_int32 :: flag_fftw
 
      plan%transform_type = p_fftw_c2r_1d
      plan%direction = 0
      if (present(normalized)) then
         plan%normalized = normalized
      else
         plan%normalized = .false.
      end if
      ! Set the information about the algorithm to compute the plan. The default is FFTW_ESTIMATE
      if (present(optimization)) then
         flag_fftw = optimization
      else
         flag_fftw = fftw_estimate
      end if
      if (present(aligned)) then
         if (aligned .EQV. .false.) then
            flag_fftw = flag_fftw + fftw_unaligned
         end if
      else
         flag_fftw = flag_fftw + fftw_unaligned
      end if
      plan%problem_rank = 1
      sll_allocate(plan%problem_shape(1), ierr)
      plan%problem_shape = (/nx/)
 
#ifdef FFTW_F2003
      plan%fftw = fftw_plan_dft_c2r_1d(nx, array_in, array_out, flag_fftw)
#else
      call dfftw_plan_dft_c2r_1d(plan%fftw, nx, array_in, array_out, flag_fftw)
#endif
 
   end subroutine sll_s_fft_init_c2r_1d
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_exec_c2r_1d(plan, array_in, array_out)
      type(sll_t_fft), intent(in)    :: plan
      sll_comp64, intent(inout) :: array_in(:) 
      sll_real64, intent(inout) :: array_out(:) 
 
      sll_real64                  :: factor
 
      sll_assert(plan%transform_type == p_fftw_c2r_1d)
 
      call fftw_execute_dft_c2r(plan%fftw, array_in, array_out)
 
      if (plan%normalized .EQV. .true.) then
         factor = 1.0_f64/real(plan%problem_shape(1), kind=f64)
         array_out = factor*array_out
      end if
   end subroutine
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_init_c2r_2d(plan, nx, ny, array_in, array_out, normalized, aligned, optimization)
      type(sll_t_fft), intent(out)     :: plan
      sll_int32, intent(in)      :: nx 
      sll_int32, intent(in)      :: ny 
      sll_comp64, intent(inout)   :: array_in(:, :) 
      sll_real64, intent(out)     :: array_out(:, :) 
      logical, optional, intent(in)      :: normalized
      logical, optional, intent(in)      :: aligned
      sll_int32, optional, intent(in)      :: optimization
      
 
      sll_int32 :: ierr
      sll_int32 :: flag_fftw
 
      plan%transform_type = p_fftw_c2r_2d
      plan%direction = 0
      if (present(normalized)) then
         plan%normalized = normalized
      else
         plan%normalized = .false.
      end if
      ! Set the information about the algorithm to compute the plan. The default is FFTW_ESTIMATE
      if (present(optimization)) then
         flag_fftw = optimization
      else
         flag_fftw = fftw_estimate
      end if
      if (present(aligned)) then
         if (aligned .EQV. .false.) then
            flag_fftw = flag_fftw + fftw_unaligned
         end if
      else
         flag_fftw = flag_fftw + fftw_unaligned
      end if
 
      plan%problem_rank = 2
      sll_allocate(plan%problem_shape(2), ierr)
      plan%problem_shape = (/nx, ny/)
#ifdef FFTW_F2003
      plan%fftw = fftw_plan_dft_c2r_2d(ny, nx, array_in, array_out, fftw_estimate + fftw_unaligned)
#else
      call dfftw_plan_dft_c2r_2d(plan%fftw, nx, ny, array_in, array_out, fftw_estimate)
#endif
   end subroutine sll_s_fft_init_c2r_2d
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_exec_c2r_2d(plan, array_in, array_out)
      type(sll_t_fft), intent(in)          :: plan
      sll_comp64, intent(inout)       :: array_in(1:, 1:)  
      sll_real64, intent(out)         :: array_out(1:, 1:) 
 
      sll_int32                                         :: nx, ny
      sll_real64                                        :: factor
 
      sll_assert(plan%transform_type == p_fftw_c2r_2d)
 
      nx = plan%problem_shape(1)
      ny = plan%problem_shape(2)
      call fftw_execute_dft_c2r(plan%fftw, array_in(1:nx/2 + 1, 1:ny), array_out(1:nx, 1:ny))
 
      if (plan%normalized .EQV. .true.) then
         factor = 1.0_f64/real(nx*ny, kind=f64)
         array_out = factor*array_out
      end if
   end subroutine
!END C2R
 
   !-----------------------------------------------------------------------------
   !-----------------------------------------------------------------------------
   subroutine sll_s_fft_free(plan)
      type(sll_t_fft), intent(inout) :: plan
 
      call fftw_destroy_plan(plan%fftw)
      if (allocated(plan%problem_shape)) then
         deallocate (plan%problem_shape)
      end if
 
   end subroutine
 
end module sll_m_fft