1! test of fftsg3d.f
2!
3      program main
4      integer nmax, nmaxsqrt
5      parameter (nmax = 128)
6      parameter (nmaxsqrt = 16)
7      integer ip(0 : nmaxsqrt + 1), n1, n2, n3, i, j
8      real*8 a(0 : nmax - 1, 0 : nmax - 1, 0 : nmax - 1),
9     &    t(0 : 8 * nmax - 1),
10     &    w(0 : nmax * 3 / 2 - 1), err, errorcheck3d
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      write (*, *) 'data length n3=? (n3 = power of 2) '
17      read (*, *) n3
18      ip(0) = 0
19!
20!   check of CDFT
21      call putdata3d(nmax, nmax, n1, n2, n3, a)
22      call cdft3d(nmax, nmax, n1, n2, n3, 1, a, t, ip, w)
23      call cdft3d(nmax, nmax, n1, n2, n3, -1, a, t, ip, w)
24      err = errorcheck3d(nmax, nmax, n1, n2, n3,
25     &    2.0d0 / n1 / n2 / n3, a)
26      write (*, *) 'cdft3d err= ', err
27!
28!   check of RDFT
29      call putdata3d(nmax, nmax, n1, n2, n3, a)
30      call rdft3d(nmax, nmax, n1, n2, n3, 1, a, t, ip, w)
31      call rdft3d(nmax, nmax, n1, n2, n3, -1, a, t, ip, w)
32      err = errorcheck3d(nmax, nmax, n1, n2, n3,
33     &    2.0d0 / n1 / n2 / n3, a)
34      write (*, *) 'rdft3d err= ', err
35!
36!   check of DDCT
37      call putdata3d(nmax, nmax, n1, n2, n3, a)
38      call ddct3d(nmax, nmax, n1, n2, n3, 1, a, t, ip, w)
39      call ddct3d(nmax, nmax, n1, n2, n3, -1, a, t, ip, w)
40      do j = 0, n2 - 1
41          do i = 0, n1 - 1
42              a(i, j, 0) = a(i, j, 0) * 0.5d0
43          end do
44      end do
45      do j = 0, n3 - 1
46          do i = 0, n1 - 1
47              a(i, 0, j) = a(i, 0, j) * 0.5d0
48          end do
49          do i = 0, n2 - 1
50              a(0, i, j) = a(0, i, j) * 0.5d0
51          end do
52      end do
53      err = errorcheck3d(nmax, nmax, n1, n2, n3,
54     &    8.0d0 / n1 / n2 / n3, a)
55      write (*, *) 'ddct3d err= ', err
56!
57!   check of DDST
58      call putdata3d(nmax, nmax, n1, n2, n3, a)
59      call ddst3d(nmax, nmax, n1, n2, n3, 1, a, t, ip, w)
60      call ddst3d(nmax, nmax, n1, n2, n3, -1, a, t, ip, w)
61      do j = 0, n2 - 1
62          do i = 0, n1 - 1
63              a(i, j, 0) = a(i, j, 0) * 0.5d0
64          end do
65      end do
66      do j = 0, n3 - 1
67          do i = 0, n1 - 1
68              a(i, 0, j) = a(i, 0, j) * 0.5d0
69          end do
70          do i = 0, n2 - 1
71              a(0, i, j) = a(0, i, j) * 0.5d0
72          end do
73      end do
74      err = errorcheck3d(nmax, nmax, n1, n2, n3,
75     &    8.0d0 / n1 / n2 / n3, a)
76      write (*, *) 'ddst3d err= ', err
77!
78      end
79!
80!
81      subroutine putdata3d(n1max, n2max, n1, n2, n3, a)
82      integer n1max, n2max, n1, n2, n3, j1, j2, j3, seed
83      real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : *), drnd
84      seed = 0
85      do j3 = 0, n3 - 1
86          do j2 = 0, n2 - 1
87              do j1 = 0, n1 - 1
88                  a(j1, j2, j3) = drnd(seed)
89              end do
90          end do
91      end do
92      end
93!
94!
95      function errorcheck3d(n1max, n2max, n1, n2, n3, scale, a)
96      integer n1max, n2max, n1, n2, n3, j1, j2, j3, seed
97      real*8 scale, a(0 : n1max - 1, 0 : n2max - 1, 0 : *),
98     &    drnd, err, e, errorcheck3d
99      err = 0
100      seed = 0
101      do j3 = 0, n3 - 1
102          do j2 = 0, n2 - 1
103              do j1 = 0, n1 - 1
104                  e = drnd(seed) - a(j1, j2, j3) * scale
105                  err = max(err, abs(e))
106              end do
107          end do
108      end do
109      errorcheck3d = err
110      end
111!
112!
113! random number generator, 0 <= drnd < 1
114      real*8 function drnd(seed)
115      integer seed
116      seed = mod(seed * 7141 + 54773, 259200)
117      drnd = seed * (1.0d0 / 259200)
118      end
119!
120