forked from patrickbwarren/SunlightHNC
-
Notifications
You must be signed in to change notification settings - Fork 0
/
fftw_test.f90
82 lines (61 loc) · 2.3 KB
/
fftw_test.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
! This file is part of SunlightDPD - a home for open source software
! related to the dissipative particle dynamics (DPD) simulation
! method.
! Copyright (c) 2009-2017 Unilever UK Central Resources Ltd
! (Registered in England & Wales, Company No 29140; Registered Office:
! Unilever House, Blackfriars, London, EC4P 4BQ, UK).
! SunlightDPD is free software: you can redistribute it and/or
! modify it under the terms of the GNU General Public License as
! published by the Free Software Foundation, either version 3 of the
! License, or (at your option) any later version.
! SunlightDPD is distributed in the hope that it will be useful, but
! WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
! General Public License for more details.
! You should have received a copy of the GNU General Public License
! along with SunlightDPD. If not, see <http://www.gnu.org/licenses/>.
! Test code to comprehend how FFTW works.
program fftw_test
implicit none
include "fftw3.f"
integer, parameter :: dp = kind(1.0d0)
integer, parameter :: n = 47
integer :: i, j, k
integer*8 :: planA, planB
real(kind=dp) :: a(n), b(n), c(n), d(n)
real(kind=dp), parameter :: pi = 4.0_dp * atan(1.0_dp)
do i=1,n-1
a(i) = 1.0_dp / sqrt(dble(i))
end do
call dfftw_plan_r2r_1d(planA, n-1, a, b, FFTW_RODFT00, FFTW_ESTIMATE)
call dfftw_plan_r2r_1d(planB, n-1, b, c, FFTW_RODFT00, FFTW_ESTIMATE)
call dfftw_execute(planA)
b = b / sqrt(2.0_dp*n)
call dfftw_execute(planB)
c = c / sqrt(2.0_dp*n)
do k = 1, n-1
d(k) = 0.0_dp
do j = 1, n-1
d(k) = d(k) + a(j) * sin(j*k*pi/n)
end do
d(k) = d(k) * sqrt(2.0_dp/n)
end do
print *, "n = ", n
print *, "fftw_rodft00 = ", fftw_rodft00
print *, "fftw_estimate = ",fftw_estimate
print *, "address of planA = ", planA
print *, "address of planB = ", planB
print *, ' a b c d'
print *, ' A FFT(FFT(A)) FFT(A) FT(A)'
do i = 1, n-1
print *, i, a(i), c(i), b(i), d(i)
end do
print *, "||a|| = ", sqrt(sum(a(:)**2))
print *, "||b|| = ", sqrt(sum(b(:)**2))
print *, "||c|| = ", sqrt(sum(c(:)**2))
print *, "||d|| = ", sqrt(sum(d(:)**2))
print *, "||a-c|| = ", sqrt(sum((a(:)-c(:))**2))
print *, "||b-d|| = ", sqrt(sum((b(:)-d(:))**2))
call dfftw_destroy_plan(planA)
call dfftw_destroy_plan(planB)
end