Context Navigation

mod_fft_fftw.F90 @ 5192

Last change on this file since 5192 was 1601, checked in by Ehouarn Millour, 13 years ago

Updates and upgrades of filter:

add minor corrections (r1591 of trunk) in filtreg_mod.F90: some arrays were oversized and the computation of the indexes from which the filter is applied could go wrong in extreme cases.
adapt filtreg_p.F (r1597 of trunk) so that BLAS routine DGEMM is only used if 'BLAS' preprocessing flag is set.
update FFT filter routines to match current 'trunk' versions (mostly cosmetic changes, exept for the implementation of use of FFTW in mod_fft_fftw.F90).
update "arch-PW6_VARGAS.fcm" to enable use of FFTW

NB: implemented FFTs assume working precision to be double precision ; trying to use them when working precision is single precision will clearly end in tragedy.
EM

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 3.0 KB

Line
1	!
2	! $Id: mod_fft_fftw.F90 1601 2011-12-07 07:05:43Z abarral $
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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format