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
|
c-----------------------------------------------------------------------
c
c-----------------------------------------------------------------------
c
parameter (SIZE = 1024, ILOOP = 100)
complex a, w
real breal(SIZE), bimag(SIZE), qbreal(SIZE), qbimag(SIZE)
c
ioutd = i1mach(2)
nn = SIZE
tpi = 8.*atan(1.)
tpion = tpi/float(nn)
w = cmplx(cos(tpion),-sin(tpion))
c
c generate a**k as test function
c result to b(i) for modification by dft and idft subroutines and
c a copy to qb(i) to compare final result with for error difference.
c
a = (.9,.3)
breal(1) = 1.0
bimag(1) = 0.0
qbreal(1) = 1.0
qbimag(1) = 0.0
do 10 k=2,nn
w = a**(k-1)
breal(k) = real(w)
bimag(k) = aimag(w)
qbreal(k) = breal(k)
qbimag(k) = bimag(k)
10 continue
c
c now compute dft, idft, dft, idft, ...
c first dft is computed specially, in case subroutine needs to be started.
c
call fft842(0, SIZE, breal, bimag)
do 25 icount = 1, ILOOP
call fft842(1, SIZE, breal, bimag)
call fft842(0, SIZE, breal, bimag)
25 continue
call fft842(1, SIZE, breal, bimag)
c
c calculate rms error between b(i) and qb(i).
c
err = 0.
do 30 i=1,nn
err = err + (breal(i)-qbreal(i))**2
* + (bimag(i)-qbimag(i))**2
30 continue
err = sqrt(err / float(nn))
write (ioutd,9994) ILOOP, err
9994 format(' rms error, after ', i5, ' loops = ', e15.8)
stop
end
|
