| ! test of fftsg3d.f |
| ! |
| program main |
| integer nmax, nmaxsqrt |
| parameter (nmax = 128) |
| parameter (nmaxsqrt = 16) |
| integer ip(0 : nmaxsqrt + 1), n1, n2, n3, i, j |
| real*8 a(0 : nmax - 1, 0 : nmax - 1, 0 : nmax - 1), |
| & t(0 : 8 * nmax - 1), |
| & w(0 : nmax * 3 / 2 - 1), err, errorcheck3d |
| ! |
| write (*, *) 'data length n1=? (n1 = power of 2) ' |
| read (*, *) n1 |
| write (*, *) 'data length n2=? (n2 = power of 2) ' |
| read (*, *) n2 |
| write (*, *) 'data length n3=? (n3 = power of 2) ' |
| read (*, *) n3 |
| ip(0) = 0 |
| ! |
| ! check of CDFT |
| call putdata3d(nmax, nmax, n1, n2, n3, a) |
| call cdft3d(nmax, nmax, n1, n2, n3, 1, a, t, ip, w) |
| call cdft3d(nmax, nmax, n1, n2, n3, -1, a, t, ip, w) |
| err = errorcheck3d(nmax, nmax, n1, n2, n3, |
| & 2.0d0 / n1 / n2 / n3, a) |
| write (*, *) 'cdft3d err= ', err |
| ! |
| ! check of RDFT |
| call putdata3d(nmax, nmax, n1, n2, n3, a) |
| call rdft3d(nmax, nmax, n1, n2, n3, 1, a, t, ip, w) |
| call rdft3d(nmax, nmax, n1, n2, n3, -1, a, t, ip, w) |
| err = errorcheck3d(nmax, nmax, n1, n2, n3, |
| & 2.0d0 / n1 / n2 / n3, a) |
| write (*, *) 'rdft3d err= ', err |
| ! |
| ! check of DDCT |
| call putdata3d(nmax, nmax, n1, n2, n3, a) |
| call ddct3d(nmax, nmax, n1, n2, n3, 1, a, t, ip, w) |
| call ddct3d(nmax, nmax, n1, n2, n3, -1, a, t, ip, w) |
| do j = 0, n2 - 1 |
| do i = 0, n1 - 1 |
| a(i, j, 0) = a(i, j, 0) * 0.5d0 |
| end do |
| end do |
| do j = 0, n3 - 1 |
| do i = 0, n1 - 1 |
| a(i, 0, j) = a(i, 0, j) * 0.5d0 |
| end do |
| do i = 0, n2 - 1 |
| a(0, i, j) = a(0, i, j) * 0.5d0 |
| end do |
| end do |
| err = errorcheck3d(nmax, nmax, n1, n2, n3, |
| & 8.0d0 / n1 / n2 / n3, a) |
| write (*, *) 'ddct3d err= ', err |
| ! |
| ! check of DDST |
| call putdata3d(nmax, nmax, n1, n2, n3, a) |
| call ddst3d(nmax, nmax, n1, n2, n3, 1, a, t, ip, w) |
| call ddst3d(nmax, nmax, n1, n2, n3, -1, a, t, ip, w) |
| do j = 0, n2 - 1 |
| do i = 0, n1 - 1 |
| a(i, j, 0) = a(i, j, 0) * 0.5d0 |
| end do |
| end do |
| do j = 0, n3 - 1 |
| do i = 0, n1 - 1 |
| a(i, 0, j) = a(i, 0, j) * 0.5d0 |
| end do |
| do i = 0, n2 - 1 |
| a(0, i, j) = a(0, i, j) * 0.5d0 |
| end do |
| end do |
| err = errorcheck3d(nmax, nmax, n1, n2, n3, |
| & 8.0d0 / n1 / n2 / n3, a) |
| write (*, *) 'ddst3d err= ', err |
| ! |
| end |
| ! |
| ! |
| subroutine putdata3d(n1max, n2max, n1, n2, n3, a) |
| integer n1max, n2max, n1, n2, n3, j1, j2, j3, seed |
| real*8 a(0 : n1max - 1, 0 : n2max - 1, 0 : *), drnd |
| seed = 0 |
| do j3 = 0, n3 - 1 |
| do j2 = 0, n2 - 1 |
| do j1 = 0, n1 - 1 |
| a(j1, j2, j3) = drnd(seed) |
| end do |
| end do |
| end do |
| end |
| ! |
| ! |
| function errorcheck3d(n1max, n2max, n1, n2, n3, scale, a) |
| integer n1max, n2max, n1, n2, n3, j1, j2, j3, seed |
| real*8 scale, a(0 : n1max - 1, 0 : n2max - 1, 0 : *), |
| & drnd, err, e, errorcheck3d |
| err = 0 |
| seed = 0 |
| do j3 = 0, n3 - 1 |
| do j2 = 0, n2 - 1 |
| do j1 = 0, n1 - 1 |
| e = drnd(seed) - a(j1, j2, j3) * scale |
| err = max(err, abs(e)) |
| end do |
| end do |
| end do |
| errorcheck3d = 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 |
| ! |