1! test of fftsg2d.f
2!
3      program main
4      integer nmax, nmaxsqrt
5      parameter (nmax = 1024)
6      parameter (nmaxsqrt = 32)
7      integer ip(0 : nmaxsqrt + 1), n1, n2, i
8      real*8 a(0 : nmax - 1, 0 : nmax - 1), t(0 : 8 * nmax - 1),
9     &    w(0 : nmax * 3 / 2 - 1), err, errorcheck2d
10!
11      write (*, *) 'data length n1=? (n1 = power of 2) '
12      read (*, *) n1
13      write (*, *) 'data length n2=? (n2 = power of 2) '
14      read (*, *) n2
15      ip(0) = 0
16!
17!   check of CDFT
18      call putdata2d(nmax, n1, n2, a)
19      call cdft2d(nmax, n1, n2, 1, a, t, ip, w)
20      call cdft2d(nmax, n1, n2, -1, a, t, ip, w)
21      err = errorcheck2d(nmax, n1, n2, 2.0d0 / n1 / n2, a)
22      write (*, *) 'cdft2d err= ', err
23!
24!   check of RDFT
25      call putdata2d(nmax, n1, n2, a)
26      call rdft2d(nmax, n1, n2, 1, a, t, ip, w)
27      call rdft2d(nmax, n1, n2, -1, a, t, ip, w)
28      err = errorcheck2d(nmax, n1, n2, 2.0d0 / n1 / n2, a)
29      write (*, *) 'rdft2d err= ', err
30!
31!   check of DDCT
32      call putdata2d(nmax, n1, n2, a)
33      call ddct2d(nmax, n1, n2, 1, a, t, ip, w)
34      call ddct2d(nmax, n1, n2, -1, a, t, ip, w)
35      do i = 0, n1 - 1
36          a(i, 0) = a(i, 0) * 0.5d0
37      end do
38      do i = 0, n2 - 1
39          a(0, i) = a(0, i) * 0.5d0
40      end do
41      err = errorcheck2d(nmax, n1, n2, 4.0d0 / n1 / n2, a)
42      write (*, *) 'ddct2d err= ', err
43!
44!   check of DDST
45      call putdata2d(nmax, n1, n2, a)
46      call ddst2d(nmax, n1, n2, 1, a, t, ip, w)
47      call ddst2d(nmax, n1, n2, -1, a, t, ip, w)
48      do i = 0, n1 - 1
49          a(i, 0) = a(i, 0) * 0.5d0
50      end do
51      do i = 0, n2 - 1
52          a(0, i) = a(0, i) * 0.5d0
53      end do
54      err = errorcheck2d(nmax, n1, n2, 4.0d0 / n1 / n2, a)
55      write (*, *) 'ddst2d err= ', err
56!
57      end
58!
59!
60      subroutine putdata2d(n1max, n1, n2, a)
61      integer n1max, n1, n2, j1, j2, seed
62      real*8 a(0 : n1max - 1, 0 : *), drnd
63      seed = 0
64      do j2 = 0, n2 - 1
65          do j1 = 0, n1 - 1
66              a(j1, j2) = drnd(seed)
67          end do
68      end do
69      end
70!
71!
72      function errorcheck2d(n1max, n1, n2, scale, a)
73      integer n1max, n1, n2, j1, j2, seed
74      real*8 scale, a(0 : n1max - 1, 0 : *), drnd, err, e,
75     &    errorcheck2d
76      err = 0
77      seed = 0
78      do j2 = 0, n2 - 1
79          do j1 = 0, n1 - 1
80              e = drnd(seed) - a(j1, j2) * scale
81              err = max(err, abs(e))
82          end do
83      end do
84      errorcheck2d = err
85      end
86!
87!
88! random number generator, 0 <= drnd < 1
89      real*8 function drnd(seed)
90      integer seed
91      seed = mod(seed * 7141 + 54773, 259200)
92      drnd = seed * (1.0d0 / 259200)
93      end
94!
95