| 1 | MODULE mod_fft_fftw |
|---|
| 2 | |
|---|
| 3 | #ifdef FFT_FFTW |
|---|
| 4 | |
|---|
| 5 | REAL,SAVE,ALLOCATABLE :: Table_forward(:) |
|---|
| 6 | REAL,SAVE,ALLOCATABLE :: Table_backward(:) |
|---|
| 7 | REAL,SAVE :: scale_factor |
|---|
| 8 | INTEGER,SAVE :: vsize |
|---|
| 9 | INTEGER,PARAMETER :: inc=1 |
|---|
| 10 | |
|---|
| 11 | INTEGER,SAVE :: plan_forward |
|---|
| 12 | INTEGER,SAVE :: plan_backward |
|---|
| 13 | |
|---|
| 14 | CONTAINS |
|---|
| 15 | |
|---|
| 16 | SUBROUTINE Init_fft(iim) |
|---|
| 17 | IMPLICIT NONE |
|---|
| 18 | #include <rfftw.h> |
|---|
| 19 | INTEGER :: iim |
|---|
| 20 | REAL :: rtmp=1. |
|---|
| 21 | COMPLEX*16 :: ctmp |
|---|
| 22 | INTEGER :: itmp=1 |
|---|
| 23 | INTEGER :: isign=0 |
|---|
| 24 | INTEGER :: ierr |
|---|
| 25 | |
|---|
| 26 | vsize=iim |
|---|
| 27 | scale_factor=1./SQRT(1.*vsize) |
|---|
| 28 | ALLOCATE(Table_forward(2*vsize+64)) |
|---|
| 29 | ALLOCATE(Table_backward(2*vsize+64)) |
|---|
| 30 | |
|---|
| 31 | ! CALL DZFFTM(isign,vsize,itmp,scale_factor,rtmp,vsize+inc,ctmp,vsize/2+1,table_forward,rtmp,ierr) |
|---|
| 32 | |
|---|
| 33 | ! CALL ZDFFTM(isign,vsize,itmp,scale_factor,ctmp,vsize/2+1,rtmp,vsize+inc,table_backward,rtmp,ierr) |
|---|
| 34 | |
|---|
| 35 | CALL rfftw_f77_create_plan(plan_forward,iim,FFTW_REAL_TO_COMPLEX,FFTW_ESTIMATE) |
|---|
| 36 | CALL rfftw_f77_create_plan(plan_backward,iim,FFTW_COMPLEX_TO_REAL,FFTW_ESTIMATE) |
|---|
| 37 | |
|---|
| 38 | END SUBROUTINE Init_fft |
|---|
| 39 | |
|---|
| 40 | |
|---|
| 41 | SUBROUTINE fft_forward(vect,TF_vect,nb_vect) |
|---|
| 42 | IMPLICIT NONE |
|---|
| 43 | #include <rfftw.h> |
|---|
| 44 | INTEGER,INTENT(IN) :: nb_vect |
|---|
| 45 | REAL,INTENT(IN) :: vect(vsize+inc,nb_vect) |
|---|
| 46 | COMPLEX*16,INTENT(OUT) :: TF_vect(vsize/2+1,nb_vect) |
|---|
| 47 | REAL :: work(4*vsize*nb_vect) |
|---|
| 48 | INTEGER :: ierr |
|---|
| 49 | INTEGER, PARAMETER :: isign=-1 |
|---|
| 50 | |
|---|
| 51 | ! CALL DZFFTM(isign,vsize,nb_vect,scale_factor,vect,vsize+inc,TF_vect,vsize/2+1,table_forward,work,ierr) |
|---|
| 52 | CALL rfftwnd_f77_real_to_complex(plan_forward,nb_vect,vect, 1, vsize+inc , TF_vect, 1, vsize/2+1); |
|---|
| 53 | |
|---|
| 54 | END SUBROUTINE fft_forward |
|---|
| 55 | |
|---|
| 56 | SUBROUTINE fft_backward(TF_vect,vect,nb_vect) |
|---|
| 57 | IMPLICIT NONE |
|---|
| 58 | #include <rfftw.h> |
|---|
| 59 | INTEGER,INTENT(IN) :: nb_vect |
|---|
| 60 | REAL,INTENT(OUT) :: vect(vsize+inc,nb_vect) |
|---|
| 61 | COMPLEX*16,INTENT(IN ) :: TF_vect(vsize/2+1,nb_vect) |
|---|
| 62 | REAL :: work(4*vsize*nb_vect) |
|---|
| 63 | INTEGER :: ierr |
|---|
| 64 | INTEGER, PARAMETER :: isign=1 |
|---|
| 65 | |
|---|
| 66 | ! CALL ZDFFTM(isign,vsize,nb_vect,scale_factor,TF_vect,vsize/2+1,vect,vsize+inc,table_backward,work,ierr) |
|---|
| 67 | CALL rfftwnd_f77_complex_to_real(plan_forward,nb_vect,TF_vect, 1, vsize/2+1 , vect, 1, vsize+inc); |
|---|
| 68 | |
|---|
| 69 | END SUBROUTINE fft_backward |
|---|
| 70 | |
|---|
| 71 | #endif |
|---|
| 72 | |
|---|
| 73 | END MODULE mod_fft_fftw |
|---|
| 74 | |
|---|
