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