# Odd(sine) Fourier transform with FFTW package

Hello everyone,

I am having problem performing forward and backward sine tranforms with fftw package. The problem is that the repeated application of the sine transform does not return the original array (even after rescaling). The same package in python works perfectly fine. Please see the code below

`````` program test

use, intrinsic :: iso_c_binding
include 'fftw3.f03'

!      implicit none

type(C_PTR) :: plan5

integer(C_FFTW_R2R_KIND) :: RODFT00

integer(C_INT), parameter :: kk = 6
real(C_DOUBLE) :: cff
real(C_DOUBLE), dimension(kk) :: in1, out1
real(C_DOUBLE), dimension(kk) :: test_in, test_out, test_out1

plan5 = fftw_plan_r2r_1d(kk,in1,out1,RODFT00,FFTW_ESTIMATE)

cff = dble(2*(kk+1))

print*, 'cff = ', cff

test_in(1) = 1d0
test_in(2) = 2d0
test_in(3) = 3d0
test_in(4) = 4d0
test_in(5) = 5d0
test_in(6) = 6d0

print*, 'test_in', test_in
call fftw_execute_r2r(plan5,test_in,test_out)

print*, 'test_out', test_out
call fftw_execute_r2r(plan5,test_out,test_out1)

print*, 'test_out1', test_out1/cff

stop

end program test
``````

The output of the program is

`````` cff =    14.0000000000000
test_in   1.00000000000000        2.00000000000000        3.00000000000000
4.00000000000000        5.00000000000000        6.00000000000000
test_out   21.0000000000000       -3.00000000000000
-3.00000000000000       -3.00000000000000        1.73205080756888
5.19615242270663
test_out1   1.35201451644825        1.83800362911206
1.25256417034730        1.46684988463302       0.214285714285714
0.799725173050474
``````

Anyone can help?

Thanks,

Mike

3 Likes

Hi @mrudko, thanks for the question and welcome to the forum!

How did you install fftw? What version? What operating system? What Fortran compiler and version? If you can help us reproduce your setup, we can help debug this.

1 Like

Welcome, @mrudko! (Mike is my colleague at U. Miami)

1 Like

Thanks, Milan!

Hello @certik ,

I installed fftw using standard procedure

``````./configure
make
make install
``````

I used Intel’s compiler (ifort)

Below is Makefile to compile the code

``````%:%.f90
\$(FC) -o \$@ \$(FFLAGS) \$<
%.o:%.f90
\$(FC) -c \$(FFLAGS) \$<

FC = ifort
FFLAGS = -free -O3 -convert big_endian
OBJS = fftw3.o qgparams.o aux.o diagnostics.o forcing.o elliptic.o correct.o rhs.o rk.o
LIBS = libfftw3.a

test: \$(OBJS) \$(LIBS) test.f90
\$(FC) \$(FFLAGS) \$(OBJS) \$(LIBS) test.f90 -o \$@
clean:
rm *.o *.mod
depend:
./sfmakedepend *.f90
# DO NOT DELETE THIS LINE - used by makedepend
# DO NOT DELETE THIS LINE - used by make depend
aux.o: fftw3.o qgparams.o

diagnostics.o: aux.o qgparams.o

elliptic.o: fftw3.o qgparams.o

fftw3.o: fftw3.f03

forcing.o: aux.o fftw3.o qgparams.o

qgparams.o: params.h
qgparams.o: fftw3.o

test.o: aux.o diagnostics.o elliptic.o FFTW3.o forcing.o qgparams.o rk.o

correct.o: qgparams.o elliptic.o

rhs.o: aux.o elliptic.o forcing.o qgparams.o

rk.o: forcing.o qgparams.o rhs.o
``````

Most of the object files are not needed for this test code except for fftw3.o. The fftw3.f90 module is

``````      module fftw3

use, intrinsic :: iso_c_binding
include 'fftw3.f03'

end module fftw3
``````

The file fftw3.f03 is generated upon installation. It defines the interfaces and constants used by fftw subroutines. libfftw3.a is a library created upon installation. Please let me know if any other information is needed.

Thanks,

Mike

1 Like

Hi @mrudko ,

`integer(C_FFTW_R2R_KIND) :: RODFT00`

This variable is undefined. I guess you could import the constant `FFTW_RODFT00` from the `fftw3.f03` file: `integer(C_FFTW_R2R_KIND), parameter :: RODFT00 = FFTW_RODFT00`.

But you can also pass directly that parameter to the planner. Here is the revised code (essentially you can just substitute `RODFT00` with `FFTW_RODFT00`), but it is also very important practice to use `implicit none` statement in the Fortran program, to make sure all variables are explicitly declared in the program.

``````program test
use, intrinsic :: iso_c_binding
implicit none
include 'fftw3.f03'
type(C_PTR) :: plan5
integer(C_INT), parameter :: kk = 6
real(C_DOUBLE) :: cff
real(C_DOUBLE), dimension(kk) :: in1, out1
real(C_DOUBLE), dimension(kk) :: test_in, test_out, test_out1
!
plan5 = fftw_plan_r2r_1d(kk,in1,out1,FFTW_RODFT00,FFTW_ESTIMATE)
cff = dble(2*(kk+1))
print*, 'cff = ', cff
!
test_in(1) = 1d0
test_in(2) = 2d0
test_in(3) = 3d0
test_in(4) = 4d0
test_in(5) = 5d0
test_in(6) = 6d0
!
print*, 'test_in', test_in
call fftw_execute_r2r(plan5,test_in,test_out)
!
print*, 'test_out', test_out
call fftw_execute_r2r(plan5,test_out,test_out1)
!
print*, 'test_out1', test_out1/cff
stop
end program test

``````

Edit: here the expected results:

``````\$ gfortran test.f90 -lfftw3 && ./a.out
cff =    14.000000000000000
test_in   1.0000000000000000        2.0000000000000000        3.0000000000000000        4.0000000000000000        5.0000000000000000        6.0000000000000000
test_out   30.669003872743762       -14.535649776006355        8.7777223636389259       -5.5823137221768286        3.3710223316526999       -1.5977043207310500
test_out1  0.99999999999999967        2.0000000000000004        3.0000000000000000        4.0000000000000009        5.0000000000000000        6.0000000000000000
``````
3 Likes

Thanks @pcosta

2 Likes

Hey guys,

I just got another weird question. Now, I am trying to take 2d sine transform with FFTW library. The code is as follows:

``````  program test

use, intrinsic :: iso_c_binding
implicit none

include 'fftw3.f03'

type(C_PTR) :: plan5, plan6

integer(C_FFTW_R2R_KIND) :: RODFT00 = FFTW_RODFT00

integer(C_INT), parameter :: kk = 6
integer(C_INT), parameter :: ll = kk*kk
real(C_DOUBLE) :: cff, cff1
real(C_DOUBLE), dimension(kk) :: in1, out1
real(C_DOUBLE), dimension(kk,kk) :: in2, out2
real(C_DOUBLE), dimension(kk) :: test_in, test_out, test_out1
real(C_DOUBLE), dimension(kk,kk) :: test_2d_in, test_2d_out, test_2d_out1

plan5 = fftw_plan_r2r_1d(kk,in1,out1,RODFT00,FFTW_ESTIMATE)
plan6 = fftw_plan_r2r_2d(kk,kk,in2,out2,RODFT00,RODFT00,FFTW_ESTIMATE)

cff = dble(2*(kk+1))
cff1 = cff**2

print*, 'cff = ', cff
print*, 'cff1 = ', cff1

call RANDOM_NUMBER(test_2d_in(:,:))

print*, 'test_2d_in', test_2d_in
call fftw_execute_r2r(plan6,test_2d_in,test_2d_out)

stop
print*, 'test_2d_out', test_2d_out
call fftw_execute_r2r(plan6,test_2d_out,test_2d_out1)

print*, 'test_2d_out1', test_2d_out1/cff1

stop

test_in(1) = 1d0
test_in(2) = 2d0
test_in(3) = 3d0
test_in(4) = 4d0
test_in(5) = 5d0
test_in(6) = 6d0

print*, 'test_in', test_in
call fftw_execute_r2r(plan5,test_in,test_out)

print*, 'test_out', test_out
call fftw_execute_r2r(plan5,test_out,test_out1)

print*, 'test_out1', test_out1/cff

stop

end program test
``````

When I try to compile this code with fort compiler it gives me the following error message

/tmp/ifortmuCvfo.o: In function `MAIN__’:

test1.f90:(.text+0x86): undefined reference to `fftw_plan_r2r_2d’

make: *** [Makefile:13: test1] Error 1

Can anybody help?

Thanks,

Mike

I was able to compile and run your file with the FFTW library installed in my system, as I wrote above.

It seems the FFTW library is not being linked properly. Are you using the same makefile as above, and is the `libfftw3.a` file still in the same folder?

PS:

best replacing

`integer(C_FFTW_R2R_KIND) :: RODFT00 = FFTW_RODFT00`

with

`integer(C_FFTW_R2R_KIND), parameter :: RODFT00 = FFTW_RODFT00`

See here why.

Best,
Pedro

Hello @pcosta.

The problem was that I was putting libfftw3.a before test1.f90 when compiling. The correct version is ifort -free -03 -convert big_endian fftw3.o test1.f90 libfftw3.a -o test1.