1 | ! |
---|
2 | ! $Id: mod_fft_fftw.F90 1910 2013-11-29 08:40:25Z fairhead $ |
---|
3 | ! |
---|
4 | |
---|
5 | MODULE mod_fft_fftw |
---|
6 | |
---|
7 | #ifdef FFT_FFTW |
---|
8 | |
---|
9 | REAL, SAVE :: scale_factor |
---|
10 | INTEGER, SAVE :: vsize |
---|
11 | INTEGER, PARAMETER :: inc=1 |
---|
12 | |
---|
13 | INTEGER*8, ALLOCATABLE, DIMENSION(:), SAVE :: plan_forward |
---|
14 | INTEGER*8, ALLOCATABLE, DIMENSION(:), SAVE :: plan_backward |
---|
15 | |
---|
16 | CONTAINS |
---|
17 | |
---|
18 | SUBROUTINE Init_fft(iim,nvectmax) |
---|
19 | IMPLICIT NONE |
---|
20 | #include <fftw3.f> |
---|
21 | INTEGER :: iim |
---|
22 | INTEGER :: nvectmax |
---|
23 | |
---|
24 | INTEGER :: itmp |
---|
25 | |
---|
26 | INTEGER :: rank |
---|
27 | INTEGER :: howmany |
---|
28 | INTEGER :: istride, idist |
---|
29 | INTEGER :: ostride, odist |
---|
30 | INTEGER, DIMENSION(1) :: n_array, inembed, onembed |
---|
31 | |
---|
32 | REAL, DIMENSION(iim+1,nvectmax) :: dbidon |
---|
33 | COMPLEX, DIMENSION(iim/2+1,nvectmax) :: cbidon |
---|
34 | |
---|
35 | vsize = iim |
---|
36 | scale_factor = 1./SQRT(1.*vsize) |
---|
37 | |
---|
38 | dbidon = 0 |
---|
39 | cbidon = 0 |
---|
40 | |
---|
41 | ALLOCATE(plan_forward(nvectmax)) |
---|
42 | ALLOCATE(plan_backward(nvectmax)) |
---|
43 | |
---|
44 | WRITE(*,*)"!---------------------!" |
---|
45 | WRITE(*,*)"! !" |
---|
46 | WRITE(*,*)"! INITIALISATION FFTW !" |
---|
47 | WRITE(*,*)"! !" |
---|
48 | WRITE(*,*)"!---------------------!" |
---|
49 | |
---|
50 | ! On initialise tous les plans |
---|
51 | DO itmp = 1, nvectmax |
---|
52 | rank = 1 |
---|
53 | n_array(1) = iim |
---|
54 | howmany = itmp |
---|
55 | inembed(1) = iim + 1 ; onembed(1) = iim/2 + 1 |
---|
56 | istride = 1 ; ostride = 1 |
---|
57 | idist = iim + 1 ; odist = iim/2 + 1 |
---|
58 | |
---|
59 | CALL dfftw_plan_many_dft_r2c(plan_forward(itmp), rank, n_array, howmany, & |
---|
60 | & dbidon, inembed, istride, idist, & |
---|
61 | & cbidon, onembed, ostride, odist, & |
---|
62 | & FFTW_ESTIMATE) |
---|
63 | |
---|
64 | rank = 1 |
---|
65 | n_array(1) = iim |
---|
66 | howmany = itmp |
---|
67 | inembed(1) = iim/2 + 1 ; onembed(1) = iim + 1 |
---|
68 | istride = 1 ; ostride = 1 |
---|
69 | idist = iim/2 + 1 ; odist = iim + 1 |
---|
70 | CALL dfftw_plan_many_dft_c2r(plan_backward(itmp), rank, n_array, howmany, & |
---|
71 | & cbidon, inembed, istride, idist, & |
---|
72 | & dbidon, onembed, ostride, odist, & |
---|
73 | & FFTW_ESTIMATE) |
---|
74 | |
---|
75 | ENDDO |
---|
76 | |
---|
77 | WRITE(*,*)"!-------------------------!" |
---|
78 | WRITE(*,*)"! !" |
---|
79 | WRITE(*,*)"! FIN INITIALISATION FFTW !" |
---|
80 | WRITE(*,*)"! !" |
---|
81 | WRITE(*,*)"!-------------------------!" |
---|
82 | |
---|
83 | END SUBROUTINE Init_fft |
---|
84 | |
---|
85 | |
---|
86 | SUBROUTINE fft_forward(vect,TF_vect,nb_vect) |
---|
87 | IMPLICIT NONE |
---|
88 | #include <fftw3.f> |
---|
89 | INTEGER,INTENT(IN) :: nb_vect |
---|
90 | REAL,INTENT(IN) :: vect(vsize+inc,nb_vect) |
---|
91 | COMPLEX,INTENT(OUT) :: TF_vect(vsize/2+1,nb_vect) |
---|
92 | |
---|
93 | CALL dfftw_execute_dft_r2c(plan_forward(nb_vect),vect,TF_vect) |
---|
94 | |
---|
95 | TF_vect = scale_factor * TF_vect |
---|
96 | |
---|
97 | END SUBROUTINE fft_forward |
---|
98 | |
---|
99 | SUBROUTINE fft_backward(TF_vect,vect,nb_vect) |
---|
100 | IMPLICIT NONE |
---|
101 | #include <fftw3.f> |
---|
102 | INTEGER,INTENT(IN) :: nb_vect |
---|
103 | REAL,INTENT(OUT) :: vect(vsize+inc,nb_vect) |
---|
104 | COMPLEX,INTENT(IN ) :: TF_vect(vsize/2+1,nb_vect) |
---|
105 | |
---|
106 | CALL dfftw_execute_dft_c2r(plan_backward(nb_vect),TF_vect,vect) |
---|
107 | |
---|
108 | vect = scale_factor * vect |
---|
109 | |
---|
110 | END SUBROUTINE fft_backward |
---|
111 | |
---|
112 | #endif |
---|
113 | |
---|
114 | END MODULE mod_fft_fftw |
---|