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
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