!-------------------------------------------------------------------------!
!                                                                         !
!        N  A  S     P A R A L L E L     B E N C H M A R K S  3.3         !
!                                                                         !
!                       O p e n M P     V E R S I O N                     !
!                                                                         !
!                                   F T                                   !
!                                                                         !
!-------------------------------------------------------------------------!
!                                                                         !
!    This benchmark is an OpenMP version of the NPB FT code.              !
!    It is described in NAS Technical Report 99-011.                      !
!                                                                         !
!    Permission to use, copy, distribute and modify this software         !
!    for any purpose with or without fee is hereby granted.  We           !
!    request, however, that all derived work reference the NAS            !
!    Parallel Benchmarks 3.3. This software is provided "as is"           !
!    without express or implied warranty.                                 !
!                                                                         !
!    Information on NPB 3.3, including the technical report, the          !
!    original specifications, source code, results and information        !
!    on how to submit new results, is available at:                       !
!                                                                         !
!           http://www.nas.nasa.gov/Software/NPB/                         !
!                                                                         !
!    Send comments or suggestions to  npb@nas.nasa.gov                    !
!                                                                         !
!          NAS Parallel Benchmarks Group                                  !
!          NASA Ames Research Center                                      !
!          Mail Stop: T27A-1                                              !
!          Moffett Field, CA   94035-1000                                 !
!                                                                         !
!          E-mail:  npb@nas.nasa.gov                                      !
!          Fax:     (650) 604-3957                                        !
!                                                                         !
!-------------------------------------------------------------------------!

c---------------------------------------------------------------------
c
c Authors: D. Bailey
c          W. Saphir
c          H. Jin
c
c---------------------------------------------------------------------

c---------------------------------------------------------------------

c---------------------------------------------------------------------
c FT benchmark
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      program ft

c---------------------------------------------------------------------
c---------------------------------------------------------------------


      implicit none

      include 'global.h'
      integer i
      
c---------------------------------------------------------------------
c u0, u1, u2 are the main arrays in the problem. 
c Depending on the decomposition, these arrays will have different 
c dimensions. To accomodate all possibilities, we allocate them as 
c one-dimensional arrays and pass them to subroutines for different 
c views
c  - u0 contains the initial (transformed) initial condition
c  - u1 and u2 are working arrays
c  - twiddle contains exponents for the time evolution operator. 
c---------------------------------------------------------------------

      double complex   u0(ntotalp), 
     >                 u1(ntotalp)
c     >                 u2(ntotalp)
      double precision twiddle(ntotalp)
c---------------------------------------------------------------------
c Large arrays are in common so that they are allocated on the
c heap rather than the stack. This common block is not
c referenced directly anywhere else. Padding is to avoid accidental 
c cache problems, since all array sizes are powers of two.
c---------------------------------------------------------------------

c      double complex pad1(3), pad2(3), pad3(3)
c      common /bigarrays/ u0, pad1, u1, pad2, u2, pad3, twiddle
      double complex pad1(3), pad2(3)
      common /bigarrays/ u0, pad1, u1, pad2, twiddle

      integer iter
      double precision total_time, mflops
      logical verified
      character class

c---------------------------------------------------------------------
c Run the entire problem once to make sure all data is touched. 
c This reduces variable startup costs, which is important for such a 
c short benchmark. The other NPB 2 implementations are similar. 
c---------------------------------------------------------------------
      do i = 1, t_max
         call timer_clear(i)
      end do
      call setup()
      call init_ui(u0, u1, twiddle, dims(1), dims(2), dims(3))
      call compute_indexmap(twiddle, dims(1), dims(2), dims(3))
      call compute_initial_conditions(u1, dims(1), dims(2), dims(3))
      call fft_init (dims(1))
      call fft(1, u1, u0)

c---------------------------------------------------------------------
c Start over from the beginning. Note that all operations must
c be timed, in contrast to other benchmarks. 
c---------------------------------------------------------------------
      do i = 1, t_max
         call timer_clear(i)
      end do
#ifdef HOOKS
      call roi_begin
#endif

      call timer_start(T_total)
      if (timers_enabled) call timer_start(T_setup)

      call compute_indexmap(twiddle, dims(1), dims(2), dims(3))

      call compute_initial_conditions(u1, dims(1), dims(2), dims(3))

      call fft_init (dims(1))

      if (timers_enabled) call timer_stop(T_setup)
      if (timers_enabled) call timer_start(T_fft)
      call fft(1, u1, u0)
      if (timers_enabled) call timer_stop(T_fft)

      do iter = 1, niter
         if (timers_enabled) call timer_start(T_evolve)
         call evolve(u0, u1, twiddle, dims(1), dims(2), dims(3))
         if (timers_enabled) call timer_stop(T_evolve)
         if (timers_enabled) call timer_start(T_fft)
c         call fft(-1, u1, u2)
         call fft(-1, u1, u1)
         if (timers_enabled) call timer_stop(T_fft)
         if (timers_enabled) call timer_start(T_checksum)
c         call checksum(iter, u2, dims(1), dims(2), dims(3))
         call checksum(iter, u1, dims(1), dims(2), dims(3))
         if (timers_enabled) call timer_stop(T_checksum)
      end do

      call verify(nx, ny, nz, niter, verified, class)

      call timer_stop(t_total)
      total_time = timer_read(t_total)

#ifdef HOOKS
      call roi_end
#endif

      if( total_time .ne. 0. ) then
         mflops = 1.0d-6*float(ntotal) *
     >             (14.8157+7.19641*log(float(ntotal))
     >          +  (5.23518+7.21113*log(float(ntotal)))*niter)
     >                 /total_time
      else
         mflops = 0.0
      endif
      call print_results('FT', class, nx, ny, nz, niter,
     >  total_time, mflops, '          floating point', verified, 
     >  npbversion, compiletime, cs1, cs2, cs3, cs4, cs5, cs6, cs7)
      if (timers_enabled) call print_timers()

      end

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine init_ui(u0, u1, twiddle, d1, d2, d3)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c touch all the big data
c---------------------------------------------------------------------

      implicit none
      integer d1, d2, d3
      double complex   u0(d1+1,d2,d3)
      double complex   u1(d1+1,d2,d3)
      double precision twiddle(d1+1,d2,d3)
      integer i, j, k

!$omp parallel do default(shared) private(i,j,k)
      do k = 1, d3
         do j = 1, d2
            do i = 1, d1
               u0(i,j,k) = 0.d0
               u1(i,j,k) = 0.d0
               twiddle(i,j,k) = 0.d0
            end do
         end do
      end do

      return
      end

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine evolve(u0, u1, twiddle, d1, d2, d3)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c evolve u0 -> u1 (t time steps) in fourier space
c---------------------------------------------------------------------

      implicit none
      include 'global.h'
      integer d1, d2, d3
      double complex   u0(d1+1,d2,d3)
      double complex   u1(d1+1,d2,d3)
      double precision twiddle(d1+1,d2,d3)
      integer i, j, k

!$omp parallel do default(shared) private(i,j,k)
      do k = 1, d3
         do j = 1, d2
            do i = 1, d1
               u0(i,j,k) = u0(i,j,k) * twiddle(i,j,k)
               u1(i,j,k) = u0(i,j,k)
            end do
         end do
      end do

      return
      end


c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine compute_initial_conditions(u0, d1, d2, d3)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c Fill in array u0 with initial conditions from 
c random number generator 
c---------------------------------------------------------------------
      implicit none
      include 'global.h'
      integer d1, d2, d3
      double complex u0(d1+1, d2, d3)
      integer k, j
      double precision x0, start, an, dummy, starts(nz)
      

      start = seed
c---------------------------------------------------------------------
c Jump to the starting element for our first plane.
c---------------------------------------------------------------------
      call ipow46(a, 0, an)
      dummy = randlc(start, an)
      call ipow46(a, 2*nx*ny, an)

      starts(1) = start
      do k = 2, dims(3)
         dummy = randlc(start, an)
         starts(k) = start
      end do
      
c---------------------------------------------------------------------
c Go through by z planes filling in one square at a time.
c---------------------------------------------------------------------
!$omp parallel do default(shared) private(k,j,x0)
      do k = 1, dims(3) 
         x0 = starts(k)
         do j = 1, dims(2) 
            call vranlc(2*nx, x0, a, u0(1, j, k))
         end do
      end do

      return
      end


c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine ipow46(a, exponent, result)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c compute a^exponent mod 2^46
c---------------------------------------------------------------------

      implicit none
      double precision a, result, dummy, q, r
      integer exponent, n, n2
      external randlc
      double precision randlc
c---------------------------------------------------------------------
c Use
c   a^n = a^(n/2)*a^(n/2) if n even else
c   a^n = a*a^(n-1)       if n odd
c---------------------------------------------------------------------
      result = 1
      if (exponent .eq. 0) return
      q = a
      r = 1
      n = exponent


      do while (n .gt. 1)
         n2 = n/2
         if (n2 * 2 .eq. n) then
            dummy = randlc(q, q) 
            n = n2
         else
            dummy = randlc(r, q)
            n = n-1
         endif
      end do
      dummy = randlc(r, q)
      result = r
      return
      end


c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine setup

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      implicit none
      include 'global.h'

      integer fstatus
!$    integer  omp_get_max_threads
!$    external omp_get_max_threads
      debug = .FALSE.

      open (unit=2,file='timer.flag',status='old',iostat=fstatus)
      if (fstatus .eq. 0) then
         timers_enabled = .true.
         close(2)
      else
         timers_enabled = .false.
      endif

      write(*, 1000)

      niter = niter_default

      write(*, 1001) nx, ny, nz
      write(*, 1002) niter
!$    write(*, 1003) omp_get_max_threads()
      write(*, *)


 1000 format(//,' NAS Parallel Benchmarks (NPB3.3-OMP)',
     >          ' - FT Benchmark', /)
 1001 format(' Size                : ', i4, 'x', i4, 'x', i4)
 1002 format(' Iterations                  :', i7)
 1003 format(' Number of available threads :', i7)

      dims(1) = nx
      dims(2) = ny
      dims(3) = nz


c---------------------------------------------------------------------
c Set up info for blocking of ffts and transposes.  This improves
c performance on cache-based systems. Blocking involves
c working on a chunk of the problem at a time, taking chunks
c along the first, second, or third dimension. 
c
c - In cffts1 blocking is on 2nd dimension (with fft on 1st dim)
c - In cffts2/3 blocking is on 1st dimension (with fft on 2nd and 3rd dims)

c Since 1st dim is always in processor, we'll assume it's long enough 
c (default blocking factor is 16 so min size for 1st dim is 16)
c The only case we have to worry about is cffts1 in a 2d decomposition. 
c so the blocking factor should not be larger than the 2nd dimension. 
c---------------------------------------------------------------------

      fftblock = fftblock_default
      fftblockpad = fftblockpad_default

      if (fftblock .ne. fftblock_default) fftblockpad = fftblock+3

      return
      end

      
c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine compute_indexmap(twiddle, d1, d2, d3)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c compute function from local (i,j,k) to ibar^2+jbar^2+kbar^2 
c for time evolution exponent. 
c---------------------------------------------------------------------

      implicit none
      include 'global.h'
      integer d1, d2, d3
      double precision twiddle(d1+1, d2, d3)
      integer i, j, k, kk, kk2, jj, kj2, ii
      double precision ap

c---------------------------------------------------------------------
c basically we want to convert the fortran indices 
c   1 2 3 4 5 6 7 8 
c to 
c   0 1 2 3 -4 -3 -2 -1
c The following magic formula does the trick:
c mod(i-1+n/2, n) - n/2
c---------------------------------------------------------------------

      ap = - 4.d0 * alpha * pi *pi

!$omp parallel do default(shared) private(i,j,k,kk,kk2,jj,kj2,ii)
      do k = 1, dims(3)
         kk =  mod(k-1+nz/2, nz) - nz/2
         kk2 = kk*kk
         do j = 1, dims(2)
            jj = mod(j-1+ny/2, ny) - ny/2
            kj2 = jj*jj+kk2
            do i = 1, dims(1)
               ii = mod(i-1+nx/2, nx) - nx/2
               twiddle(i,j,k) = dexp(ap*dble(ii*ii+kj2))
            end do
         end do
      end do

      return
      end


c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine print_timers()

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      implicit none
      integer i
      include 'global.h'
      double precision t, t_m
      character*25 tstrings(T_max)
      data tstrings / '          total ', 
     >                '          setup ', 
     >                '            fft ', 
     >                '         evolve ', 
     >                '       checksum ', 
     >                '           fftx ', 
     >                '           ffty ', 
     >                '           fftz ' /

      t_m = timer_read(T_total)
      if (t_m .le. 0.0d0) t_m = 1.0d0
      do i = 1, t_max
         t = timer_read(i)
         write(*, 100) i, tstrings(i), t, t*100.0/t_m
      end do
 100  format(' timer ', i2, '(', A16,  ') :', F9.4, ' (',F6.2,'%)')
      return
      end


c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine fft(dir, x1, x2)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      implicit none
      include 'global.h'
      integer dir
      double complex x1(ntotalp), x2(ntotalp)

      double complex y1(fftblockpad_default*maxdim),
     >               y2(fftblockpad_default*maxdim)

c---------------------------------------------------------------------
c note: args x1, x2 must be different arrays
c note: args for cfftsx are (direction, layout, xin, xout, scratch)
c       xin/xout may be the same and it can be somewhat faster
c       if they are
c---------------------------------------------------------------------

      if (dir .eq. 1) then
         call cffts1(1, dims(1), dims(2), dims(3), x1, x1, y1, y2)
         call cffts2(1, dims(1), dims(2), dims(3), x1, x1, y1, y2)
         call cffts3(1, dims(1), dims(2), dims(3), x1, x2, y1, y2)
      else
         call cffts3(-1, dims(1), dims(2), dims(3), x1, x1, y1, y2)
         call cffts2(-1, dims(1), dims(2), dims(3), x1, x1, y1, y2)
         call cffts1(-1, dims(1), dims(2), dims(3), x1, x2, y1, y2)
      endif
      return
      end


c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine cffts1(is, d1, d2, d3, x, xout, y1, y2)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      implicit none

      include 'global.h'
      integer is, d1, d2, d3, logd1
      double complex x(d1+1,d2,d3)
      double complex xout(d1+1,d2,d3)
      double complex y1(fftblockpad, d1), y2(fftblockpad, d1)
      integer i, j, k, jj

      logd1 = ilog2(d1)

      if (timers_enabled) call timer_start(T_fftx)
!$omp parallel do default(shared) private(i,j,k,jj,y1,y2)
!$omp&  shared(is,logd1,d1)
      do k = 1, d3
         do jj = 0, d2 - fftblock, fftblock
            do j = 1, fftblock
               do i = 1, d1
                  y1(j,i) = x(i,j+jj,k)
               enddo
            enddo
            
            call cfftz (is, logd1, d1, y1, y2)


            do j = 1, fftblock
               do i = 1, d1
                  xout(i,j+jj,k) = y1(j,i)
               enddo
            enddo
         enddo
      enddo
      if (timers_enabled) call timer_stop(T_fftx)

      return
      end


c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine cffts2(is, d1, d2, d3, x, xout, y1, y2)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      implicit none

      include 'global.h'
      integer is, d1, d2, d3, logd2
      double complex x(d1+1,d2,d3)
      double complex xout(d1+1,d2,d3)
      double complex y1(fftblockpad, d2), y2(fftblockpad, d2)
      integer i, j, k, ii

      logd2 = ilog2(d2)

      if (timers_enabled) call timer_start(T_ffty)
!$omp parallel do default(shared) private(i,j,k,ii,y1,y2)
!$omp&  shared(is,logd2,d2)
      do k = 1, d3
        do ii = 0, d1 - fftblock, fftblock
           do j = 1, d2
              do i = 1, fftblock
                 y1(i,j) = x(i+ii,j,k)
              enddo
           enddo

           call cfftz (is, logd2, d2, y1, y2)
           
           do j = 1, d2
              do i = 1, fftblock
                 xout(i+ii,j,k) = y1(i,j)
              enddo
           enddo
        enddo
      enddo
      if (timers_enabled) call timer_stop(T_ffty)

      return
      end


c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine cffts3(is, d1, d2, d3, x, xout, y1, y2)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      implicit none

      include 'global.h'
      integer is, d1, d2, d3, logd3
      double complex x(d1+1,d2,d3)
      double complex xout(d1+1,d2,d3)
      double complex y1(fftblockpad, d3), y2(fftblockpad, d3)
      integer i, j, k, ii

      logd3 = ilog2(d3)

      if (timers_enabled) call timer_start(T_fftz)
!$omp parallel do default(shared) private(i,j,k,ii,y1,y2)
!$omp&  shared(is)
      do j = 1, d2
        do ii = 0, d1 - fftblock, fftblock
           do k = 1, d3
              do i = 1, fftblock
                 y1(i,k) = x(i+ii,j,k)
              enddo
           enddo

           call cfftz (is, logd3, d3, y1, y2)

           do k = 1, d3
              do i = 1, fftblock
                 xout(i+ii,j,k) = y1(i,k)
              enddo
           enddo
        enddo
      enddo
      if (timers_enabled) call timer_stop(T_fftz)

      return
      end


c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine fft_init (n)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c compute the roots-of-unity array that will be used for subsequent FFTs. 
c---------------------------------------------------------------------

      implicit none
      include 'global.h'

      integer m,n,nu,ku,i,j,ln
      double precision t, ti


c---------------------------------------------------------------------
c   Initialize the U array with sines and cosines in a manner that permits
c   stride one access at each FFT iteration.
c---------------------------------------------------------------------
      nu = n
      m = ilog2(n)
      u(1) = m
      ku = 2
      ln = 1

      do j = 1, m
         t = pi / ln
         
         do i = 0, ln - 1
            ti = i * t
            u(i+ku) = dcmplx (cos (ti), sin(ti))
         enddo
         
         ku = ku + ln
         ln = 2 * ln
      enddo
      
      return
      end

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine cfftz (is, m, n, x, y)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c   Computes NY N-point complex-to-complex FFTs of X using an algorithm due
c   to Swarztrauber.  X is both the input and the output array, while Y is a 
c   scratch array.  It is assumed that N = 2^M.  Before calling CFFTZ to 
c   perform FFTs, the array U must be initialized by calling CFFTZ with IS 
c   set to 0 and M set to MX, where MX is the maximum value of M for any 
c   subsequent call.
c---------------------------------------------------------------------

      implicit none
      include 'global.h'

      integer is,m,n,i,j,l,mx
      double complex x, y

      dimension x(fftblockpad,n), y(fftblockpad,n)

c---------------------------------------------------------------------
c   Check if input parameters are invalid.
c---------------------------------------------------------------------
      mx = u(1)
      if ((is .ne. 1 .and. is .ne. -1) .or. m .lt. 1 .or. m .gt. mx)    
     >  then
        write (*, 1)  is, m, mx
 1      format ('CFFTZ: Either U has not been initialized, or else'/    
     >    'one of the input parameters is invalid', 3I5)
        stop
      endif

c---------------------------------------------------------------------
c   Perform one variant of the Stockham FFT.
c---------------------------------------------------------------------
      do l = 1, m, 2
        call fftz2 (is, l, m, n, fftblock, fftblockpad, u, x, y)
        if (l .eq. m) goto 160
        call fftz2 (is, l + 1, m, n, fftblock, fftblockpad, u, y, x)
      enddo

      goto 180

c---------------------------------------------------------------------
c   Copy Y to X.
c---------------------------------------------------------------------
 160  do j = 1, n
        do i = 1, fftblock
          x(i,j) = y(i,j)
        enddo
      enddo

 180  continue

      return
      end

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine fftz2 (is, l, m, n, ny, ny1, u, x, y)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

c---------------------------------------------------------------------
c   Performs the L-th iteration of the second variant of the Stockham FFT.
c---------------------------------------------------------------------

      implicit none

      integer is,k,l,m,n,ny,ny1,n1,li,lj,lk,ku,i,j,i11,i12,i21,i22
      double complex u,x,y,u1,x11,x21
      dimension u(n), x(ny1,n), y(ny1,n)


c---------------------------------------------------------------------
c   Set initial parameters.
c---------------------------------------------------------------------

      n1 = n / 2
      lk = 2 ** (l - 1)
      li = 2 ** (m - l)
      lj = 2 * lk
      ku = li + 1

      do i = 0, li - 1
        i11 = i * lk + 1
        i12 = i11 + n1
        i21 = i * lj + 1
        i22 = i21 + lk
        if (is .ge. 1) then
          u1 = u(ku+i)
        else
          u1 = dconjg (u(ku+i))
        endif

c---------------------------------------------------------------------
c   This loop is vectorizable.
c---------------------------------------------------------------------
        do k = 0, lk - 1
          do j = 1, ny
            x11 = x(j,i11+k)
            x21 = x(j,i12+k)
            y(j,i21+k) = x11 + x21
            y(j,i22+k) = u1 * (x11 - x21)
          enddo
        enddo
      enddo

      return
      end

c---------------------------------------------------------------------


c---------------------------------------------------------------------
c---------------------------------------------------------------------

      integer function ilog2(n)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      implicit none
      integer n, nn, lg
      if (n .eq. 1) then
         ilog2=0
         return
      endif
      lg = 1
      nn = 2
      do while (nn .lt. n)
         nn = nn*2
         lg = lg+1
      end do
      ilog2 = lg
      return
      end


c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine checksum(i, u1, d1, d2, d3)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      implicit none
      include 'global.h'
      integer i, d1, d2, d3
      double complex u1(d1+1,d2,d3)
      integer j, q,r,s
      double complex chk
      chk = (0.0,0.0)

!$omp parallel do default(shared) private(i,q,r,s) reduction(+:chk)
      do j=1,1024
         q = mod(j, nx)+1
         r = mod(3*j,ny)+1
         s = mod(5*j,nz)+1
         chk=chk+u1(q,r,s)
      end do

      chk = chk/dble(ntotal)
      
      write (*, 30) i, chk
 30   format (' T =',I5,5X,'Checksum =',1P2D22.12)
      sums(i) = chk
      return
      end


c---------------------------------------------------------------------
c---------------------------------------------------------------------

      subroutine verify (d1, d2, d3, nt, verified, class)

c---------------------------------------------------------------------
c---------------------------------------------------------------------

      implicit none
      include 'global.h'
      integer d1, d2, d3, nt
      character class
      logical verified
      integer i
      double precision err, epsilon

c---------------------------------------------------------------------
c   Reference checksums
c---------------------------------------------------------------------
      double complex csum_ref(25)


      class = 'U'

      epsilon = 1.0d-12
      verified = .FALSE.

      if (d1 .eq. 64 .and.
     >    d2 .eq. 64 .and.
     >    d3 .eq. 64 .and.
     >    nt .eq. 6) then
c---------------------------------------------------------------------
c   Sample size reference checksums
c---------------------------------------------------------------------
         class = 'S'
         csum_ref(1) = dcmplx(5.546087004964D+02, 4.845363331978D+02)
         csum_ref(2) = dcmplx(5.546385409189D+02, 4.865304269511D+02)
         csum_ref(3) = dcmplx(5.546148406171D+02, 4.883910722336D+02)
         csum_ref(4) = dcmplx(5.545423607415D+02, 4.901273169046D+02)
         csum_ref(5) = dcmplx(5.544255039624D+02, 4.917475857993D+02)
         csum_ref(6) = dcmplx(5.542683411902D+02, 4.932597244941D+02)

      else if (d1 .eq. 128 .and.
     >    d2 .eq. 128 .and.
     >    d3 .eq. 32 .and.
     >    nt .eq. 6) then
c---------------------------------------------------------------------
c   Class W size reference checksums
c---------------------------------------------------------------------
         class = 'W'
         csum_ref(1) = dcmplx(5.673612178944D+02, 5.293246849175D+02)
         csum_ref(2) = dcmplx(5.631436885271D+02, 5.282149986629D+02)
         csum_ref(3) = dcmplx(5.594024089970D+02, 5.270996558037D+02)
         csum_ref(4) = dcmplx(5.560698047020D+02, 5.260027904925D+02)
         csum_ref(5) = dcmplx(5.530898991250D+02, 5.249400845633D+02)
         csum_ref(6) = dcmplx(5.504159734538D+02, 5.239212247086D+02)

      else if (d1 .eq. 256 .and.
     >    d2 .eq. 256 .and.
     >    d3 .eq. 128 .and.
     >    nt .eq. 6) then
c---------------------------------------------------------------------
c   Class A size reference checksums
c---------------------------------------------------------------------
         class = 'A'
         csum_ref(1) = dcmplx(5.046735008193D+02, 5.114047905510D+02)
         csum_ref(2) = dcmplx(5.059412319734D+02, 5.098809666433D+02)
         csum_ref(3) = dcmplx(5.069376896287D+02, 5.098144042213D+02)
         csum_ref(4) = dcmplx(5.077892868474D+02, 5.101336130759D+02)
         csum_ref(5) = dcmplx(5.085233095391D+02, 5.104914655194D+02)
         csum_ref(6) = dcmplx(5.091487099959D+02, 5.107917842803D+02)
      
      else if (d1 .eq. 512 .and.
     >    d2 .eq. 256 .and.
     >    d3 .eq. 256 .and.
     >    nt .eq. 20) then
c---------------------------------------------------------------------
c   Class B size reference checksums
c---------------------------------------------------------------------
         class = 'B'
         csum_ref(1)  = dcmplx(5.177643571579D+02, 5.077803458597D+02)
         csum_ref(2)  = dcmplx(5.154521291263D+02, 5.088249431599D+02)
         csum_ref(3)  = dcmplx(5.146409228649D+02, 5.096208912659D+02)
         csum_ref(4)  = dcmplx(5.142378756213D+02, 5.101023387619D+02)
         csum_ref(5)  = dcmplx(5.139626667737D+02, 5.103976610617D+02)
         csum_ref(6)  = dcmplx(5.137423460082D+02, 5.105948019802D+02)
         csum_ref(7)  = dcmplx(5.135547056878D+02, 5.107404165783D+02)
         csum_ref(8)  = dcmplx(5.133910925466D+02, 5.108576573661D+02)
         csum_ref(9)  = dcmplx(5.132470705390D+02, 5.109577278523D+02)
         csum_ref(10) = dcmplx(5.131197729984D+02, 5.110460304483D+02)
         csum_ref(11) = dcmplx(5.130070319283D+02, 5.111252433800D+02)
         csum_ref(12) = dcmplx(5.129070537032D+02, 5.111968077718D+02)
         csum_ref(13) = dcmplx(5.128182883502D+02, 5.112616233064D+02)
         csum_ref(14) = dcmplx(5.127393733383D+02, 5.113203605551D+02)
         csum_ref(15) = dcmplx(5.126691062020D+02, 5.113735928093D+02)
         csum_ref(16) = dcmplx(5.126064276004D+02, 5.114218460548D+02)
         csum_ref(17) = dcmplx(5.125504076570D+02, 5.114656139760D+02)
         csum_ref(18) = dcmplx(5.125002331720D+02, 5.115053595966D+02)
         csum_ref(19) = dcmplx(5.124551951846D+02, 5.115415130407D+02)
         csum_ref(20) = dcmplx(5.124146770029D+02, 5.115744692211D+02)

      else if (d1 .eq. 512 .and.
     >    d2 .eq. 512 .and.
     >    d3 .eq. 512 .and.
     >    nt .eq. 20) then
c---------------------------------------------------------------------
c   Class C size reference checksums
c---------------------------------------------------------------------
         class = 'C'
         csum_ref(1)  = dcmplx(5.195078707457D+02, 5.149019699238D+02)
         csum_ref(2)  = dcmplx(5.155422171134D+02, 5.127578201997D+02)
         csum_ref(3)  = dcmplx(5.144678022222D+02, 5.122251847514D+02)
         csum_ref(4)  = dcmplx(5.140150594328D+02, 5.121090289018D+02)
         csum_ref(5)  = dcmplx(5.137550426810D+02, 5.121143685824D+02)
         csum_ref(6)  = dcmplx(5.135811056728D+02, 5.121496764568D+02)
         csum_ref(7)  = dcmplx(5.134569343165D+02, 5.121870921893D+02)
         csum_ref(8)  = dcmplx(5.133651975661D+02, 5.122193250322D+02)
         csum_ref(9)  = dcmplx(5.132955192805D+02, 5.122454735794D+02)
         csum_ref(10) = dcmplx(5.132410471738D+02, 5.122663649603D+02)
         csum_ref(11) = dcmplx(5.131971141679D+02, 5.122830879827D+02)
         csum_ref(12) = dcmplx(5.131605205716D+02, 5.122965869718D+02)
         csum_ref(13) = dcmplx(5.131290734194D+02, 5.123075927445D+02)
         csum_ref(14) = dcmplx(5.131012720314D+02, 5.123166486553D+02)
         csum_ref(15) = dcmplx(5.130760908195D+02, 5.123241541685D+02)
         csum_ref(16) = dcmplx(5.130528295923D+02, 5.123304037599D+02)
         csum_ref(17) = dcmplx(5.130310107773D+02, 5.123356167976D+02)
         csum_ref(18) = dcmplx(5.130103090133D+02, 5.123399592211D+02)
         csum_ref(19) = dcmplx(5.129905029333D+02, 5.123435588985D+02)
         csum_ref(20) = dcmplx(5.129714421109D+02, 5.123465164008D+02)

      else if (d1 .eq. 2048 .and.
     >    d2 .eq. 1024 .and.
     >    d3 .eq. 1024 .and.
     >    nt .eq. 25) then
c---------------------------------------------------------------------
c   Class D size reference checksums
c---------------------------------------------------------------------
         class = 'D'
         csum_ref(1)  = dcmplx(5.122230065252D+02, 5.118534037109D+02)
         csum_ref(2)  = dcmplx(5.120463975765D+02, 5.117061181082D+02)
         csum_ref(3)  = dcmplx(5.119865766760D+02, 5.117096364601D+02)
         csum_ref(4)  = dcmplx(5.119518799488D+02, 5.117373863950D+02)
         csum_ref(5)  = dcmplx(5.119269088223D+02, 5.117680347632D+02)
         csum_ref(6)  = dcmplx(5.119082416858D+02, 5.117967875532D+02)
         csum_ref(7)  = dcmplx(5.118943814638D+02, 5.118225281841D+02)
         csum_ref(8)  = dcmplx(5.118842385057D+02, 5.118451629348D+02)
         csum_ref(9)  = dcmplx(5.118769435632D+02, 5.118649119387D+02)
         csum_ref(10) = dcmplx(5.118718203448D+02, 5.118820803844D+02)
         csum_ref(11) = dcmplx(5.118683569061D+02, 5.118969781011D+02)
         csum_ref(12) = dcmplx(5.118661708593D+02, 5.119098918835D+02)
         csum_ref(13) = dcmplx(5.118649768950D+02, 5.119210777066D+02)
         csum_ref(14) = dcmplx(5.118645605626D+02, 5.119307604484D+02)
         csum_ref(15) = dcmplx(5.118647586618D+02, 5.119391362671D+02)
         csum_ref(16) = dcmplx(5.118654451572D+02, 5.119463757241D+02)
         csum_ref(17) = dcmplx(5.118665212451D+02, 5.119526269238D+02)
         csum_ref(18) = dcmplx(5.118679083821D+02, 5.119580184108D+02)
         csum_ref(19) = dcmplx(5.118695433664D+02, 5.119626617538D+02)
         csum_ref(20) = dcmplx(5.118713748264D+02, 5.119666538138D+02)
         csum_ref(21) = dcmplx(5.118733606701D+02, 5.119700787219D+02)
         csum_ref(22) = dcmplx(5.118754661974D+02, 5.119730095953D+02)
         csum_ref(23) = dcmplx(5.118776626738D+02, 5.119755100241D+02)
         csum_ref(24) = dcmplx(5.118799262314D+02, 5.119776353561D+02)
         csum_ref(25) = dcmplx(5.118822370068D+02, 5.119794338060D+02)

      else if (d1 .eq. 4096 .and.
     >    d2 .eq. 2048 .and.
     >    d3 .eq. 2048 .and.
     >    nt .eq. 25) then
c---------------------------------------------------------------------
c   Class E size reference checksums
c---------------------------------------------------------------------
         class = 'E'
         csum_ref(1)  = dcmplx(5.121601045346D+02, 5.117395998266D+02)
         csum_ref(2)  = dcmplx(5.120905403678D+02, 5.118614716182D+02)
         csum_ref(3)  = dcmplx(5.120623229306D+02, 5.119074203747D+02)
         csum_ref(4)  = dcmplx(5.120438418997D+02, 5.119345900733D+02)
         csum_ref(5)  = dcmplx(5.120311521872D+02, 5.119551325550D+02)
         csum_ref(6)  = dcmplx(5.120226088809D+02, 5.119720179919D+02)
         csum_ref(7)  = dcmplx(5.120169296534D+02, 5.119861371665D+02)
         csum_ref(8)  = dcmplx(5.120131225172D+02, 5.119979364402D+02)
         csum_ref(9)  = dcmplx(5.120104767108D+02, 5.120077674092D+02)
         csum_ref(10) = dcmplx(5.120085127969D+02, 5.120159443121D+02)
         csum_ref(11) = dcmplx(5.120069224127D+02, 5.120227453670D+02)
         csum_ref(12) = dcmplx(5.120055158164D+02, 5.120284096041D+02)
         csum_ref(13) = dcmplx(5.120041820159D+02, 5.120331373793D+02)
         csum_ref(14) = dcmplx(5.120028605402D+02, 5.120370938679D+02)
         csum_ref(15) = dcmplx(5.120015223011D+02, 5.120404138831D+02)
         csum_ref(16) = dcmplx(5.120001570022D+02, 5.120432068837D+02)
         csum_ref(17) = dcmplx(5.119987650555D+02, 5.120455615860D+02)
         csum_ref(18) = dcmplx(5.119973525091D+02, 5.120475499442D+02)
         csum_ref(19) = dcmplx(5.119959279472D+02, 5.120492304629D+02)
         csum_ref(20) = dcmplx(5.119945006558D+02, 5.120506508902D+02)
         csum_ref(21) = dcmplx(5.119930795911D+02, 5.120518503782D+02)
         csum_ref(22) = dcmplx(5.119916728462D+02, 5.120528612016D+02)
         csum_ref(23) = dcmplx(5.119902874185D+02, 5.120537101195D+02)
         csum_ref(24) = dcmplx(5.119889291565D+02, 5.120544194514D+02)
         csum_ref(25) = dcmplx(5.119876028049D+02, 5.120550079284D+02)

      endif


      if (class .ne. 'U') then

         do i = 1, nt
            err = abs( (sums(i) - csum_ref(i)) / csum_ref(i) )
            if (.not.(err .le. epsilon)) goto 100
         end do
         verified = .TRUE.
 100     continue

      endif

         
      if (class .ne. 'U') then
         if (verified) then
            write(*,2000)
 2000       format(' Result verification successful')
         else
            write(*,2001)
 2001       format(' Result verification failed')
         endif
      endif
      print *, 'class = ', class

      return
      end