Does real dft work with sizes different from power of 2?

[ShatrovOA]

Hi all!

I saw in documentation that real dft works only with even numbers, but my test shows that it fails when size of transform /= 2**n. Example can be found below. Program has to be linked with kfr_capi and fftw3

module kfr
use iso_c_binding
implicit none

integer(c_int), parameter :: KFR_PACK_PERM = 0
integer(c_int), parameter :: KFR_PACK_CCS = 1
! C_INT8_T
  interface 
    type(c_ptr) function kfr_dft_create_plan_f32(size) bind(C)
    import
    integer(c_size_t), value :: size
    end function kfr_dft_create_plan_f32

    integer(c_size_t) function kfr_dft_get_temp_size_f32(plan) bind(C)
      import
      type(c_ptr), value :: plan
    end function

    subroutine kfr_dft_execute_f32(plan, out, in, temp) bind(C)
      import
      type(c_ptr), value :: plan
      type(c_ptr), value :: out
      type(c_ptr), value :: in
      type(c_ptr), value :: temp
    end subroutine kfr_dft_execute_f32

    subroutine kfr_dft_execute_inverse_f32(plan, out, in, temp) bind(C)
      import
      type(c_ptr), value :: plan
      type(c_ptr), value :: out
      type(c_ptr), value :: in
      type(c_ptr), value :: temp
    end subroutine kfr_dft_execute_inverse_f32

    subroutine kfr_dft_delete_plan_f32(plan) bind(C)
      import
      type(c_ptr), value :: plan
    end subroutine kfr_dft_delete_plan_f32

    type(c_ptr) function kfr_dft_create_plan_f64(size) bind(C)
    import
    integer(c_size_t), value :: size
    end function kfr_dft_create_plan_f64

    integer(c_size_t) function kfr_dft_get_temp_size_f64(plan) bind(C)
      import
      type(c_ptr), value :: plan
    end function

    subroutine kfr_dft_execute_f64(plan, out, in, temp) bind(C)
      import
      type(c_ptr), value :: plan
      type(c_ptr), value :: out
      type(c_ptr), value :: in
      type(c_ptr), value :: temp
    end subroutine kfr_dft_execute_f64

    subroutine kfr_dft_execute_inverse_f64(plan, out, in, temp) bind(C)
      import
      type(c_ptr), value :: plan
      type(c_ptr), value :: out
      type(c_ptr), value :: in
      type(c_ptr), value :: temp
    end subroutine kfr_dft_execute_inverse_f64

    subroutine kfr_dft_delete_plan_f64(plan) bind(C)
      import
      type(c_ptr), value :: plan
    end subroutine kfr_dft_delete_plan_f64

    type(c_ptr) function kfr_dft_real_create_plan_f32(size, pack_format) bind(C)
      import
      integer(c_size_t),  value :: size
      integer(c_int),        value :: pack_format
    end function kfr_dft_real_create_plan_f32

    type(c_ptr) function kfr_dft_real_create_plan_f64(size, pack_format) bind(C)
      import
      integer(c_size_t),  value :: size
      integer(c_int),        value :: pack_format
    end function kfr_dft_real_create_plan_f64

    integer(c_size_t) function kfr_dft_real_get_temp_size_f32(plan) bind(C)
      import
      type(c_ptr), value :: plan
    end function kfr_dft_real_get_temp_size_f32

    integer(c_size_t) function kfr_dft_real_get_temp_size_f64(plan) bind(C)
      import
      type(c_ptr), value :: plan
    end function kfr_dft_real_get_temp_size_f64

    ! subroutine kfr_dft_execute_f64(plan, out, in, temp) bind(C)
    !   import
    !   type(c_ptr), value :: plan
    !   type(c_ptr), value :: out
    !   type(c_ptr), value :: in
    !   type(c_ptr), value :: temp
    ! end subroutine kfr_dft_execute_f64

    ! subroutine kfr_dft_execute_inverse_f64(plan, out, in, temp) bind(C)
    !   import
    !   type(c_ptr), value :: plan
    !   type(c_ptr), value :: out
    !   type(c_ptr), value :: in
    !   type(c_ptr), value :: temp
    ! end subroutine kfr_dft_execute_inverse_f64

    subroutine kfr_dft_real_execute_f32(plan, out, in, temp) bind(C)
      import
      type(c_ptr), value :: plan
      type(c_ptr), value :: out
      type(c_ptr), value :: in
      type(c_ptr), value :: temp
    end subroutine

    subroutine kfr_dft_real_execute_inverse_f32(plan, out, in, temp) bind(C)
      import
      type(c_ptr), value :: plan
      type(c_ptr), value :: out
      type(c_ptr), value :: in
      type(c_ptr), value :: temp
    end subroutine

    
    subroutine kfr_dft_real_execute_f64(plan, out, in, temp) bind(C)
      import
      type(c_ptr), value :: plan
      type(c_ptr), value :: out
      type(c_ptr), value :: in
      type(c_ptr), value :: temp
    end subroutine

        
    subroutine kfr_dft_real_execute_inverse_f64(plan, out, in, temp) bind(C)
      import
      type(c_ptr), value :: plan
      type(c_ptr), value :: out
      type(c_ptr), value :: in
      type(c_ptr), value :: temp
    end subroutine
  endinterface

end module kfr

module fftw
use iso_c_binding
include "fftw3.f03"
end module fftw

program main
  use kfr
  use fftw
  use iso_fortran_env
  implicit none
  integer, parameter :: n = 14
  real(real64), pointer :: kfr_in(:), check(:), fftw_in(:)
  complex(real64), pointer :: kfr_out(:), fftw_out(:)
  integer(c_size_t) :: s, i
  type(c_ptr) :: plan, ckfr_in, ckfr_out, ctemp, plan_fftwr, plan_fftwc
  integer(C_INT8_T), pointer :: temp(:)
  real(real64) :: rnd


  allocate(kfr_in(n), check(n), kfr_out(int(n / 2) + 1))
  allocate(fftw_in(n), fftw_out(int(n/2) + 1))
  plan = kfr_dft_real_create_plan_f64(int(n, c_size_t), KFR_PACK_CCS)

  ! plan = kfr_dft_create_plan_f64(int(n, c_size_t))

  s = kfr_dft_real_get_temp_size_f64(plan)
  ! s = kfr_dft_get_temp_size_f64(plan)
  print*,s
  allocate(temp(s))
  
  plan_fftwr = fftw_plan_dft_r2c_1d(n, fftw_in, fftw_out, FFTW_ESTIMATE)
  plan_fftwc = fftw_plan_dft_c2r_1d(n, fftw_out, fftw_in, FFTW_ESTIMATE)

  do i = 1, n
    call random_number(rnd)
    kfr_in(i) = rnd
  enddo
  check = kfr_in
  fftw_in = kfr_in

  ! in = cmplx(1.0, 1.0)

  ckfr_in = c_loc(kfr_in)
  ckfr_out = c_loc(kfr_out)
  ctemp = c_loc(temp)

  print*,kfr_in

  call kfr_dft_real_execute_f64(plan, ckfr_out, ckfr_in, ctemp)
  ! call kfr_dft_execute_f64(plan, cin, cin, ctemp)
  ! print*,'temp = ',temp(1:20)
  print*,'kfr complex = ', kfr_out

  call fftw_execute_dft_r2c(plan_fftwr, fftw_in, fftw_out)

  print*, 'fftw_complex = ', fftw_out

  call kfr_dft_real_execute_inverse_f64(plan, ckfr_in, ckfr_out, ctemp)
  ! call kfr_dft_execute_inverse_f64(plan, cin, cin, ctemp)

  kfr_in = kfr_in / real(n, real64)

  print*,'kfr real = ',kfr_in

  call fftw_execute_dft_c2r(plan_fftwc, fftw_out, fftw_in)

  fftw_in = fftw_in / real(n, real64)
  print*,'fftw real = ', fftw_in

  print*,'kfr error = ',maxval(abs(check-kfr_in))
  print*,'fftw error = ',maxval(abs(check-fftw_in))

  ! call kfr_dft_delete_plan_f64(plan)
end program main

[gopher2008]

According to my test, realdft ONLY works when input size equals to 4*N. I think it is a bug.