! test of fft4f2d.f!      program main      integer nmax, nmaxsqrt      parameter (nmax = 1024)      parameter (nmaxsqrt = 32)      integer ip(0 : nmaxsqrt + 1), n1, n2, i      real*8 a(0 : nmax - 1, 0 : nmax - 1),      &    t(0 : nmax - 1, 0 : nmax - 1), w(0 : nmax * 3 / 2 - 1),      &    err, errorcheck2d!      write (*, *) 'data length n1=? (n1 = power of 2) '      read (*, *) n1      write (*, *) 'data length n2=? (n2 = power of 2) '      read (*, *) n2      ip(0) = 0!!   check of CDFT      call putdata2d(nmax, n1, n2, a)      call cdft2d(nmax, n1, n2, 1, a, ip, w)      call cdft2d(nmax, n1, n2, -1, a, ip, w)      err = errorcheck2d(nmax, n1, n2, 2.0d0 / n1 / n2, a)      write (*, *) 'cdft2d err= ', err!!   check of RDFT      call putdata2d(nmax, n1, n2, a)      call rdft2d(nmax, n1, n2, 1, a, ip, w)      call rdft2d(nmax, n1, n2, -1, a, ip, w)      err = errorcheck2d(nmax, n1, n2, 2.0d0 / n1 / n2, a)      write (*, *) 'rdft2d err= ', err!!   check of DDCT      call putdata2d(nmax, n1, n2, a)      call ddct2d(nmax, n1, n2, 1, a, t, ip, w)      call ddct2d(nmax, n1, n2, -1, a, t, ip, w)      do i = 0, n1 - 1          a(i, 0) = a(i, 0) * 0.5d0      end do      do i = 0, n2 - 1          a(0, i) = a(0, i) * 0.5d0      end do      err = errorcheck2d(nmax, n1, n2, 4.0d0 / n1 / n2, a)      write (*, *) 'ddct2d err= ', err!!   check of DDST      call putdata2d(nmax, n1, n2, a)      call ddst2d(nmax, n1, n2, 1, a, t, ip, w)      call ddst2d(nmax, n1, n2, -1, a, t, ip, w)      do i = 0, n1 - 1          a(i, 0) = a(i, 0) * 0.5d0      end do      do i = 0, n2 - 1          a(0, i) = a(0, i) * 0.5d0      end do      err = errorcheck2d(nmax, n1, n2, 4.0d0 / n1 / n2, a)      write (*, *) 'ddst2d err= ', err!      end!!      subroutine putdata2d(n1max, n1, n2, a)      integer n1max, n1, n2, j1, j2, seed      real*8 a(0 : n1max - 1, 0 : *), drnd      seed = 0      do j2 = 0, n2 - 1          do j1 = 0, n1 - 1              a(j1, j2) = drnd(seed)          end do      end do      end!!      function errorcheck2d(n1max, n1, n2, scale, a)      integer n1max, n1, n2, j1, j2, seed      real*8 scale, a(0 : n1max - 1, 0 : *), drnd, err, e,      &    errorcheck2d      err = 0      seed = 0      do j2 = 0, n2 - 1          do j1 = 0, n1 - 1              e = drnd(seed) - a(j1, j2) * scale              err = max(err, abs(e))          end do      end do      errorcheck2d = err      end!!! random number generator, 0 <= drnd < 1      real*8 function drnd(seed)      integer seed      seed = mod(seed * 7141 + 54773, 259200)      drnd = seed * (1.0d0 / 259200)      end!