-
Notifications
You must be signed in to change notification settings - Fork 0
/
fft_fftw_serial.f90
162 lines (109 loc) · 3.91 KB
/
fft_fftw_serial.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
module fft_pak
! Arrow of time adjusted......8 Nov 04
!
! This one is for the FFTW 3.0.1 fma package, uniprocessor, double precision
! This one tries to use the same plan_f, plan_b all the time.
!
! Routine: Getfft
!
! Oliver's notes: Here it is. There is one explicit kind selection
! for the Fortran integer plan that is used as a C pointer by fftw; this
! must be 32 bit on my G4. Otherwise this is almost exactly as fftw2,
! so I don't know why they made such a fuss about not making it
! backwards compatible.
include 'fftw3.f'
integer*8 :: plan_f,plan_b
! FFTW needs plan to be pointer on machine. 32bits for G4? Yes.
!integer (kind=3) :: plan_f,plan_b
integer :: hres
integer,dimension(2) :: nsize
private
save
public :: Init_fft, fft
interface fft
module procedure fft2, fft3
end interface
contains
!********************************************************************
subroutine Init_fft(kmax)
! Initialize fft routine
! Does little for version 3
integer,intent(in) :: kmax
complex,dimension(:,:), allocatable :: f,fr
hres = 2*kmax+2
nsize(1) = hres
nsize(2) = hres
allocate(f(nsize(1),nsize(2)))
allocate(fr(nsize(1),nsize(2)))
f=0.
fr=0.
call dfftw_plan_dft_2d(plan_b,nsize(1),nsize(2),f,fr,&
FFTW_BACKWARD,FFTW_MEASURE)
call dfftw_plan_dft_2d(plan_f,nsize(1),nsize(2),f,fr,&
FFTW_FORWARD,FFTW_MEASURE)
deallocate(f)
deallocate(fr)
end subroutine Init_fft
!********************************************************************
function fft2(f,dirn) result(fr)
! Calculate 2d complex to complex fft. dirn = +1 or -1.
! these values correspond to sign of exponent in spectral
! sum - arrow of time?
complex,dimension(:,:) :: f
complex,dimension(size(f,1),size(f,2)) :: fr
integer,intent(in) :: dirn
real :: scale=1.0
fr=0.
if (dirn==-1) then
scale=1.0
call dfftw_execute_dft(plan_b,f,fr)
elseif (dirn==1) then
scale=1.0/float(hres)**2
call dfftw_execute_dft(plan_f,f,fr)
endif
fr = scale*fr
end function fft2
function fft3(f,dirn) result(fr)
! Calculate 2d complex to complex fft. dirn = +1 or -1.
! these values correspond to sign of exponent in spectral
! sum (see man page on zzfft2d).
complex,dimension(:,:,:) :: f
complex,dimension(size(f,1),size(f,2),size(f,3)) :: fr
complex,dimension(size(f,1),size(f,2)) :: g,gr
integer,intent(in) :: dirn
real :: scale
fr=0.
if (dirn==1) scale=1.0/float(nx*ny)
if (dirn==-1) scale=1.0
do n = 1,size(f,3)
! gr=0.
! g=0.
! g(1:nx,1:ny)=f(1:nx,1:ny,n)
! call ccfft2d(dirn,nx,ny,scale,g,nx1,gr,nx1,table,work,0)
g=f(:,:,n)
gr=fft2(g,dirn)
fr(:,:,n)=gr(:,:)
enddo
end function fft3
!********************************************************************
function fft3bis(f,dirn) result(fr)
! Calculate 3d complex to complex fft. dirn = +1 or -1.
! these values correspond to sign of exponent in spectral
! sum (see man page on ccfft2d).
complex,dimension(:,:,:) :: f
complex,dimension(size(f,1),size(f,2),size(f,3)) :: fr
integer,intent(in) :: dirn
real :: scale=1.0
do n=1,size(f,3)
if (dirn==-1) then
scale=1.0
call dfftw_execute_dft(plan_b,f(:,:,n),fr(:,:,n))
elseif (dirn==1) then
scale=1.0/float(hres)**2
call dfftw_execute_dft(plan_f,f(:,:,n),fr(:,:,n))
endif
enddo
fr = scale*fr
end function fft3bis
!********************************************************************
end module fft_pak