Odd(sine) Fourier transform with FFTW package

mrudko · June 3, 2021, 9:02pm

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

certik · June 3, 2021, 10:23pm

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.

milancurcic · June 3, 2021, 10:30pm

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

mrudko · June 3, 2021, 10:46pm

Thanks, Milan!

mrudko · June 3, 2021, 10:56pm

Hello @certik ,

I installed fftw using standard procedure

./configure
make 
make install

The package was downloaded from here http://www.fftw.org/fftw-3.3.9.tar.gz

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

pcosta · June 4, 2021, 9:03pm

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

mrudko · June 6, 2021, 9:25pm

Thanks @pcosta

mrudko · June 13, 2021, 9:16am

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

pcosta · June 14, 2021, 9:49am

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

mrudko · June 15, 2021, 1:06am

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.

Thanks for your help anyways,

Mike

themos · June 15, 2021, 12:27pm

Lesson learned: Always quote the exact thing you have tried when asking for help.

Topic		Replies	Views
Spectral derivative in 2-d Help	10	453	December 13, 2023
Old Code - Need help Help	9	375	November 23, 2023
Problems with coarrays	20	1518	February 4, 2023
DGBMV blas issue with 2d array Help	7	373	May 23, 2023
Curious result in a simple program Help	10	413	April 12, 2022

Odd(sine) Fourier transform with FFTW package

Related Topics