!-------------------------------------------------------------------------!
!                                                                         !
!        N  A  S     P A R A L L E L     B E N C H M A R K S  3.3         !
!                                                                         !
!                                   F T                                   !
!                                                                         !
!-------------------------------------------------------------------------!
!                                                                         !
!    This benchmark is part of the NAS Parallel Benchmark 3.3 suite.      !
!    It is described in NAS Technical Reports 95-020 and 02-007           !
!                                                                         !
!    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                                        !
!                                                                         !
!-------------------------------------------------------------------------!

!TO REDUCE THE AMOUNT OF MEMORY REQUIRED BY THE BENCHMARK WE NO LONGER
!STORE THE ENTIRE TIME EVOLUTION ARRAY "EX" FOR ALL TIME STEPS, BUT
!JUST FOR THE FIRST. ALSO, IT IS STORED ONLY FOR THE PART OF THE GRID
!FOR WHICH THE CALLING PROCESSOR IS RESPONSIBLE, SO THAT THE MEMORY 
!USAGE BECOMES SCALABLE. THIS NEW ARRAY IS CALLED "TWIDDLE" (SEE
!NPB3.0-SER)

!TO AVOID PROBLEMS WITH VERY LARGE ARRAY SIZES THAT ARE COMPUTED BY
!MULTIPLYING GRID DIMENSIONS (CAUSING INTEGER OVERFLOW IN THE VARIABLE
!NTOTAL) AND SUBSEQUENTLY DIVIDING BY THE NUMBER OF PROCESSORS, WE
!COMPUTE THE SIZE OF ARRAY PARTITIONS MORE CONSERVATIVELY AS
!((NX*NY)/NP)*NZ, WHERE NX, NY, AND NZ ARE GRID DIMENSIONS AND NP IS
!THE NUMBER OF PROCESSORS, THE RESULT IS STORED IN "NTDIVNP". FOR THE 
!PERFORMANCE CALCULATION WE STORE THE TOTAL NUMBER OF GRID POINTS IN A 
!FLOATING POINT NUMBER "NTOTAL_F" INSTEAD OF AN INTEGER.
!THIS FIX WILL FAIL IF THE NUMBER OF PROCESSORS IS SMALL.

!UGLY HACK OF SUBROUTINE IPOW46: FOR VERY LARGE GRIDS THE SINGLE EXPONENT
!FROM NPB2.3 MAY NOT FIT IN A 32-BIT INTEGER. HOWEVER, WE KNOW THAT THE
!"EXPONENT" ARGUMENT OF THIS ROUTINE CAN ALWAYS BE FACTORED INTO A TERM 
!DIVISIBLE BY NX (EXP_1) AND ANOTHER TERM (EXP_2). NX IS USUALLY A POWER
!OF TWO, SO WE CAN KEEP HALVING IT UNTIL THE PRODUCT OF EXP_1
!AND EXP_2 IS SMALL ENOUGH (NAMELY EXP_2 ITSELF). THIS UPDATED VERSION
!OF IPWO46, WHICH NOW TAKES THE TWO FACTORS OF "EXPONENT" AS SEPARATE
!ARGUMENTS, MAY BREAK DOWN IF EXP_1 DOES NOT CONTAIN A LARGE POWER OF TWO.

c---------------------------------------------------------------------
c
c Authors: D. Bailey
c          W. Saphir
c          R. F. Van der Wijngaart
c
c---------------------------------------------------------------------

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

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

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

      program ft

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

      implicit none

      include 'mpif.h'
      include 'global.h'
      integer i, ierr
      
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---------------------------------------------------------------------

      double complex   u0(ntdivnp), 
     >                 u1(ntdivnp), 
     >                 u2(ntdivnp)
      double precision twiddle(ntdivnp)
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---------------------------------------------------------------------

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

      integer iter
      double precision total_time, mflops
      logical verified
      character class

      call MPI_Init(ierr)

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 timer_start(T_init)
      call setup()
      call compute_indexmap(twiddle, dims(1,3), dims(2,3), dims(3,3))
      call compute_initial_conditions(u1, dims(1,1), dims(2,1), 
     >                                dims(3,1))
      call fft_init (dims(1,1))
      call fft(1, u1, u0)
      call timer_stop(T_init)
      if (me .eq. 0) then
         print *,'Initialization time =',timer_read(T_init)
      endif

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
      call MPI_Barrier(MPI_COMM_WORLD, ierr)

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

      call compute_indexmap(twiddle, dims(1,3), dims(2,3), dims(3,3))
      call compute_initial_conditions(u1, dims(1,1), dims(2,1), 
     >                                dims(3,1))
      call fft_init (dims(1,1))

!      if (timers_enabled) call synchup()
      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,1), dims(2,1), dims(3,1))
         if (timers_enabled) call timer_stop(T_evolve)
         if (timers_enabled) call timer_start(T_fft)
         call fft(-1, u1, u2)
         if (timers_enabled) call timer_stop(T_fft)
!         if (timers_enabled) call synchup()
         if (timers_enabled) call timer_start(T_checksum)
         call checksum(iter, u2, dims(1,1), dims(2,1), dims(3,1))
         if (timers_enabled) call timer_stop(T_checksum)
      end do

      call verify(nx, ny, nz, niter, verified, class)
      call timer_stop(t_total)
      if (np .ne. np_min) verified = .false.
      total_time = timer_read(t_total)

      if( total_time .ne. 0. ) then
         mflops = 1.0d-6*ntotal_f *
     >             (14.8157+7.19641*log(ntotal_f)
     >          +  (5.23518+7.21113*log(ntotal_f))*niter)
     >                 /total_time
      else
         mflops = 0.0
      endif
      if (me .eq. 0) then
         call print_results('FT', class, nx, ny, nz, niter, np_min, np,
     >     total_time, mflops, '          floating point', verified, 
     >     npbversion, compiletime, cs1, cs2, cs3, cs4, cs5, cs6, cs7)
      endif
      if (timers_enabled) call print_timers()
      call MPI_Finalize(ierr)
      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 precision exi
      double complex u0(d1,d2,d3)
      double complex u1(d1,d2,d3)
      double precision twiddle(d1,d2,d3)
      integer 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, d2, d3)
      integer k
      double precision x0, start, an, dummy
      
c---------------------------------------------------------------------
c 0-D and 1-D layouts are easy because each processor gets a contiguous
c chunk of the array, in the Fortran ordering sense. 
c For a 2-D layout, it's a bit more complicated. We always
c have entire x-lines (contiguous) in processor. 
c We can do ny/np1 of them at a time since we have
c ny/np1 contiguous in y-direction. But then we jump
c by z-planes (nz/np2 of them, total). 
c For the 0-D and 1-D layouts we could do larger chunks, but
c this turns out to have no measurable impact on performance. 
c---------------------------------------------------------------------


      start = seed                                    
c---------------------------------------------------------------------
c Jump to the starting element for our first plane.
c---------------------------------------------------------------------
      call ipow46(a, 2*nx, (zstart(1)-1)*ny + (ystart(1)-1), an)
      dummy = randlc(start, an)
      call ipow46(a, 2*nx, ny, an)
      
c---------------------------------------------------------------------
c Go through by z planes filling in one square at a time.
c---------------------------------------------------------------------
      do k = 1, dims(3, 1) ! nz/np2
         x0 = start
         call vranlc(2*nx*dims(2, 1), x0, a, u0(1, 1, k))
         if (k .ne. dims(3, 1)) dummy = randlc(start, an)
      end do

      return
      end


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

      subroutine ipow46(a, exp_1, exp_2, result)

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

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

      implicit none
      double precision a, result, dummy, q, r
      integer exp_1, exp_2, n, n2, ierr
      external randlc
      double precision randlc
      logical  two_pow
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 (exp_2 .eq. 0 .or. exp_1 .eq. 0) return
      q = a
      r = 1
      n = exp_1
      two_pow = .true.

      do while (two_pow)
         n2 = n/2
         if (n2 * 2 .eq. n) then
            dummy = randlc(q, q)
            n = n2
         else
            n = n * exp_2
            two_pow = .false.
         endif
      end do

      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 'mpinpb.h'
      include 'global.h'

      integer ierr, i, j, fstatus
      debug = .FALSE.
      
      call MPI_Comm_size(MPI_COMM_WORLD, np, ierr)
      call MPI_Comm_rank(MPI_COMM_WORLD, me, ierr)

      if (.not. convertdouble) then
         dc_type = MPI_DOUBLE_COMPLEX
      else
         dc_type = MPI_COMPLEX
      endif


      if (me .eq. 0) then
         write(*, 1000)

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

         open (unit=2,file='inputft.data',status='old', iostat=fstatus)

         if (fstatus .eq. 0) then
            write(*,233) 
 233        format(' Reading from input file inputft.data')
            read (2,*) niter
            read (2,*) layout_type
            read (2,*) np1, np2
            close(2)

c---------------------------------------------------------------------
c check to make sure input data is consistent
c---------------------------------------------------------------------

    
c---------------------------------------------------------------------
c 1. product of processor grid dims must equal number of processors
c---------------------------------------------------------------------

            if (np1 * np2 .ne. np) then
               write(*, 238)
 238           format(' np1 and np2 given in input file are not valid.')
               write(*, 239) np1*np2, np
 239           format(' Product is ', i5, ' and should be ', i5)
               call MPI_Abort(MPI_COMM_WORLD, 1, ierr)
            endif

c---------------------------------------------------------------------
c 2. layout type must be valid
c---------------------------------------------------------------------

            if (layout_type .ne. layout_0D .and.
     >          layout_type .ne. layout_1D .and.
     >          layout_type .ne. layout_2D) then
               write(*, 240)
 240           format(' Layout type specified in inputft.data is 
     >                  invalid ')
               call MPI_Abort(MPI_COMM_WORLD, 1, ierr)
            endif

c---------------------------------------------------------------------
c 3. 0D layout must be 1x1 grid
c---------------------------------------------------------------------

            if (layout_type .eq. layout_0D .and.
     >            (np1 .ne.1 .or. np2 .ne. 1)) then
               write(*, 241)
 241           format(' For 0D layout, both np1 and np2 must be 1 ')
               call MPI_Abort(MPI_COMM_WORLD, 1, ierr)
            endif
c---------------------------------------------------------------------
c 4. 1D layout must be 1xN grid
c---------------------------------------------------------------------

            if (layout_type .eq. layout_1D .and. np1 .ne. 1) then
               write(*, 242)
 242           format(' For 1D layout, np1 must be 1 ')
               call MPI_Abort(MPI_COMM_WORLD, 1, ierr)
            endif

         else
            write(*,234) 
            niter = niter_default
            if (np .eq. 1) then
               np1 = 1
               np2 = 1
               layout_type = layout_0D
            else if (np .le. nz) then
               np1 = 1
               np2 = np
               layout_type = layout_1D
            else
               np1 = nz
               np2 = np/nz
               layout_type = layout_2D
            endif
         endif

         if (np .lt. np_min) then
            write(*, 10) np_min
 10         format(' Error: Compiled for ', I5, ' processors. ')
            write(*, 11) np
 11         format(' Only ',  i5, ' processors found ')
            call MPI_Abort(MPI_COMM_WORLD, 1, ierr)
         endif

 234     format(' No input file inputft.data. Using compiled defaults')
         write(*, 1001) nx, ny, nz
         write(*, 1002) niter
         write(*, 1004) np
         write(*, 1005) np1, np2
         if (np .ne. np_min) write(*, 1006) np_min

         if (layout_type .eq. layout_0D) then
            write(*, 1010) '0D'
         else if (layout_type .eq. layout_1D) then
            write(*, 1010) '1D'
         else
            write(*, 1010) '2D'
         endif

 1000 format(//,' NAS Parallel Benchmarks 3.3 -- FT Benchmark',/)
 1001    format(' Size                : ', i4, 'x', i4, 'x', i4)
 1002    format(' Iterations          : ', 7x, i7)
 1004    format(' Number of processes : ', 7x, i7)
 1005    format(' Processor array     : ', 5x, i4, 'x', i4)
 1006    format(' WARNING: compiled for ', i5, ' processes. ',
     >          ' Will not verify. ')
 1010    format(' Layout type         : ', 9x, A5)
      endif


c---------------------------------------------------------------------
c Broadcast parameters 
c---------------------------------------------------------------------
      call MPI_BCAST(np1, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
      call MPI_BCAST(np2, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
      call MPI_BCAST(layout_type, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, 
     &               ierr)
      call MPI_BCAST(niter, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
      call MPI_BCAST(timers_enabled, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD,
     &               ierr)

      if (np1 .eq. 1 .and. np2 .eq. 1) then
        layout_type = layout_0D
      else if (np1 .eq. 1) then
         layout_type = layout_1D
      else
         layout_type = layout_2D
      endif

      if (layout_type .eq. layout_0D) then
         do i = 1, 3
            dims(1, i) = nx
            dims(2, i) = ny
            dims(3, i) = nz
         end do
      else if (layout_type .eq. layout_1D) then
         dims(1, 1) = nx
         dims(2, 1) = ny
         dims(3, 1) = nz

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

         dims(1, 3) = nz
         dims(2, 3) = nx
         dims(3, 3) = ny
      else if (layout_type .eq. layout_2D) then
         dims(1, 1) = nx
         dims(2, 1) = ny
         dims(3, 1) = nz

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

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

      endif
      do i = 1, 3
         dims(2, i) = dims(2, i) / np1
         dims(3, i) = dims(3, i) / np2
      end do


c---------------------------------------------------------------------
c Determine processor coordinates of this processor
c Processor grid is np1xnp2. 
c Arrays are always (n1, n2/np1, n3/np2)
c Processor coords are zero-based. 
c---------------------------------------------------------------------
      me2 = mod(me, np2)  ! goes from 0...np2-1
      me1 = me/np2        ! goes from 0...np1-1
c---------------------------------------------------------------------
c Communicators for rows/columns of processor grid. 
c commslice1 is communicator of all procs with same me1, ranked as me2
c commslice2 is communicator of all procs with same me2, ranked as me1
c mpi_comm_split(comm, color, key, ...)
c---------------------------------------------------------------------
      call MPI_Comm_split(MPI_COMM_WORLD, me1, me2, commslice1, ierr)
      call MPI_Comm_split(MPI_COMM_WORLD, me2, me1, commslice2, ierr)
!      if (timers_enabled) call synchup()

      if (debug) print *, 'proc coords: ', me, me1, me2

c---------------------------------------------------------------------
c Determine which section of the grid is owned by this
c processor. 
c---------------------------------------------------------------------
      if (layout_type .eq. layout_0d) then

         do i = 1, 3
            xstart(i) = 1
            xend(i)   = nx
            ystart(i) = 1
            yend(i)   = ny
            zstart(i) = 1
            zend(i)   = nz
         end do

      else if (layout_type .eq. layout_1d) then

         xstart(1) = 1
         xend(1)   = nx
         ystart(1) = 1
         yend(1)   = ny
         zstart(1) = 1 + me2 * nz/np2
         zend(1)   = (me2+1) * nz/np2

         xstart(2) = 1
         xend(2)   = nx
         ystart(2) = 1
         yend(2)   = ny
         zstart(2) = 1 + me2 * nz/np2
         zend(2)   = (me2+1) * nz/np2

         xstart(3) = 1
         xend(3)   = nx
         ystart(3) = 1 + me2 * ny/np2
         yend(3)   = (me2+1) * ny/np2
         zstart(3) = 1
         zend(3)   = nz

      else if (layout_type .eq. layout_2d) then

         xstart(1) = 1
         xend(1)   = nx
         ystart(1) = 1 + me1 * ny/np1
         yend(1)   = (me1+1) * ny/np1
         zstart(1) = 1 + me2 * nz/np2
         zend(1)   = (me2+1) * nz/np2

         xstart(2) = 1 + me1 * nx/np1
         xend(2)   = (me1+1)*nx/np1
         ystart(2) = 1
         yend(2)   = ny
         zstart(2) = zstart(1)
         zend(2)   = zend(1)

         xstart(3) = xstart(2)
         xend(3)   = xend(2)
         ystart(3) = 1 + me2 *ny/np2
         yend(3)   = (me2+1)*ny/np2
         zstart(3) = 1
         zend(3)   = nz
      endif

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 (layout_type .eq. layout_2d) then
         if (dims(2, 1) .lt. fftblock) fftblock = dims(2, 1)
         if (dims(2, 2) .lt. fftblock) fftblock = dims(2, 2)
         if (dims(2, 3) .lt. fftblock) fftblock = dims(2, 3)
      endif
      
      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 'mpinpb.h'
      include 'global.h'
      integer d1, d2, d3
      integer i, j, k, ii, ii2, jj, ij2, kk
      double precision ap, twiddle(d1, d2, d3)

c---------------------------------------------------------------------
c this function is very different depending on whether 
c we are in the 0d, 1d or 2d layout. Compute separately. 
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

      if (layout_type .eq. layout_0d) then ! xyz layout
         do i = 1, dims(1,3)
            ii =  mod(i+xstart(3)-2+nx/2, nx) - nx/2
            ii2 = ii*ii
            do j = 1, dims(2,3)
               jj = mod(j+ystart(3)-2+ny/2, ny) - ny/2
               ij2 = jj*jj+ii2
               do k = 1, dims(3,3)
                  kk = mod(k+zstart(3)-2+nz/2, nz) - nz/2
                  twiddle(i,j,k) = dexp(ap*dfloat(kk*kk+ij2))
               end do
            end do
         end do
      else if (layout_type .eq. layout_1d) then ! zxy layout 
         do i = 1,dims(2,3)
            ii =  mod(i+xstart(3)-2+nx/2, nx) - nx/2
            ii2 = ii*ii
            do j = 1,dims(3,3)
               jj = mod(j+ystart(3)-2+ny/2, ny) - ny/2
               ij2 = jj*jj+ii2
               do k = 1,dims(1,3)
                  kk = mod(k+zstart(3)-2+nz/2, nz) - nz/2
                  twiddle(k,i,j) = dexp(ap*dfloat(kk*kk+ij2))
               end do
            end do
         end do
      else if (layout_type .eq. layout_2d) then ! zxy layout
         do i = 1,dims(2,3)
            ii =  mod(i+xstart(3)-2+nx/2, nx) - nx/2
            ii2 = ii*ii
            do j = 1, dims(3,3)
               jj = mod(j+ystart(3)-2+ny/2, ny) - ny/2
               ij2 = jj*jj+ii2
               do k =1,dims(1,3)
                  kk = mod(k+zstart(3)-2+nz/2, nz) - nz/2
                  twiddle(k,i,j) = dexp(ap*dfloat(kk*kk+ij2))
               end do
            end do
         end do
      else
         print *, ' Unknown layout type ', layout_type
         stop
      endif

      return
      end


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

      subroutine print_timers()

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

      implicit none
      integer i, ierr
      include 'global.h'
      include 'mpinpb.h'
      character*25 tstrings(T_max+2)
      double precision t1(T_max+2), tsum(T_max+2),
     >                 tming(T_max+2), tmaxg(T_max+2)
      data tstrings / '          total ', 
     >                '          setup ', 
     >                '            fft ', 
     >                '         evolve ', 
     >                '       checksum ', 
     >                '         fftlow ', 
     >                '        fftcopy ', 
     >                '      transpose ', 
     >                ' transpose1_loc ', 
     >                ' transpose1_glo ', 
     >                ' transpose1_fin ', 
     >                ' transpose2_loc ', 
     >                ' transpose2_glo ', 
     >                ' transpose2_fin ', 
     >                '           sync ',
     >                '           init ',
     >                '        totcomp ',
     >                '        totcomm ' /

      do i = 1, t_max
         t1(i) = timer_read(i)
      end do
      t1(t_max+2) = t1(t_transxzglo) + t1(t_transxyglo) + t1(t_synch)
      t1(t_max+1) = t1(t_total) - t1(t_max+2)

      call MPI_Reduce(t1, tsum,  t_max+2, MPI_DOUBLE_PRECISION, 
     >                MPI_SUM, 0, MPI_COMM_WORLD, ierr)
      call MPI_Reduce(t1, tming, t_max+2, MPI_DOUBLE_PRECISION, 
     >                MPI_MIN, 0, MPI_COMM_WORLD, ierr)
      call MPI_Reduce(t1, tmaxg, t_max+2, MPI_DOUBLE_PRECISION, 
     >                MPI_MAX, 0, MPI_COMM_WORLD, ierr)

      if (me .ne. 0) return
      write(*, 800) np
      do i = 1, t_max+2
         if (tsum(i) .ne. 0.0d0) then
            write(*, 810) i, tstrings(i), tming(i), tmaxg(i), tsum(i)/np
         endif
      end do
 800  format(' nprocs =', i6, 19x, 'minimum', 5x, 'maximum', 
     >       5x, 'average')
 810  format(' timer ', i2, '(', A16, ') :', 3(2X,F10.4))
      return
      end


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

      subroutine fft(dir, x1, x2)

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

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

      double complex scratch(fftblockpad_default*maxdim*2)

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 note: args for transpose are (layout1, layout2, xin, xout)
c       xin/xout must be different
c---------------------------------------------------------------------

      if (dir .eq. 1) then
         if (layout_type .eq. layout_0d) then
            call cffts1(1, dims(1,1), dims(2,1), dims(3,1), 
     >                  x1, x1, scratch)
            call cffts2(1, dims(1,2), dims(2,2), dims(3,2), 
     >                  x1, x1, scratch)
            call cffts3(1, dims(1,3), dims(2,3), dims(3,3), 
     >                  x1, x2, scratch)
         else if (layout_type .eq. layout_1d) then
            call cffts1(1, dims(1,1), dims(2,1), dims(3,1), 
     >                  x1, x1, scratch)
            call cffts2(1, dims(1,2), dims(2,2), dims(3,2), 
     >                  x1, x1, scratch)
            if (timers_enabled) call timer_start(T_transpose)
            call transpose_xy_z(2, 3, x1, x2)
            if (timers_enabled) call timer_stop(T_transpose)
            call cffts1(1, dims(1,3), dims(2,3), dims(3,3), 
     >                  x2, x2, scratch)
         else if (layout_type .eq. layout_2d) then
            call cffts1(1, dims(1,1), dims(2,1), dims(3,1), 
     >                  x1, x1, scratch)
            if (timers_enabled) call timer_start(T_transpose)
            call transpose_x_y(1, 2, x1, x2)
            if (timers_enabled) call timer_stop(T_transpose)
            call cffts1(1, dims(1,2), dims(2,2), dims(3,2), 
     >                  x2, x2, scratch)
            if (timers_enabled) call timer_start(T_transpose)
            call transpose_x_z(2, 3, x2, x1)
            if (timers_enabled) call timer_stop(T_transpose)
            call cffts1(1, dims(1,3), dims(2,3), dims(3,3), 
     >                  x1, x2, scratch)
         endif
      else
         if (layout_type .eq. layout_0d) then
            call cffts3(-1, dims(1,3), dims(2,3), dims(3,3), 
     >                  x1, x1, scratch)
            call cffts2(-1, dims(1,2), dims(2,2), dims(3,2), 
     >                  x1, x1, scratch)
            call cffts1(-1, dims(1,1), dims(2,1), dims(3,1), 
     >                  x1, x2, scratch)
         else if (layout_type .eq. layout_1d) then
            call cffts1(-1, dims(1,3), dims(2,3), dims(3,3), 
     >                  x1, x1, scratch)
            if (timers_enabled) call timer_start(T_transpose)
            call transpose_x_yz(3, 2, x1, x2)
            if (timers_enabled) call timer_stop(T_transpose)
            call cffts2(-1, dims(1,2), dims(2,2), dims(3,2), 
     >                  x2, x2, scratch)
            call cffts1(-1, dims(1,1), dims(2,1), dims(3,1), 
     >                  x2, x2, scratch)
         else if (layout_type .eq. layout_2d) then
            call cffts1(-1, dims(1,3), dims(2,3), dims(3,3), 
     >                  x1, x1, scratch)
            if (timers_enabled) call timer_start(T_transpose)
            call transpose_x_z(3, 2, x1, x2)
            if (timers_enabled) call timer_stop(T_transpose)
            call cffts1(-1, dims(1,2), dims(2,2), dims(3,2), 
     >                  x2, x2, scratch)
            if (timers_enabled) call timer_start(T_transpose)
            call transpose_x_y(2, 1, x2, x1)
            if (timers_enabled) call timer_stop(T_transpose)
            call cffts1(-1, dims(1,1), dims(2,1), dims(3,1), 
     >                  x1, x2, scratch)
         endif
      endif
      return
      end


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

      subroutine cffts1(is, d1, d2, d3, x, xout, y)

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

      implicit none

      include 'global.h'
      integer is, d1, d2, d3, logd1
      double complex x(d1,d2,d3)
      double complex xout(d1,d2,d3)
      double complex y(fftblockpad, d1, 2) 
      integer i, j, k, jj

      logd1 = ilog2(d1)

      do k = 1, d3
         do jj = 0, d2 - fftblock, fftblock
            if (timers_enabled) call timer_start(T_fftcopy)
            do j = 1, fftblock
               do i = 1, d1
                  y(j,i,1) = x(i,j+jj,k)
               enddo
            enddo
            if (timers_enabled) call timer_stop(T_fftcopy)
            
            if (timers_enabled) call timer_start(T_fftlow)
            call cfftz (is, logd1, d1, y, y(1,1,2))
            if (timers_enabled) call timer_stop(T_fftlow)

            if (timers_enabled) call timer_start(T_fftcopy)
            do j = 1, fftblock
               do i = 1, d1
                  xout(i,j+jj,k) = y(j,i,1)
               enddo
            enddo
            if (timers_enabled) call timer_stop(T_fftcopy)
         enddo
      enddo

      return
      end


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

      subroutine cffts2(is, d1, d2, d3, x, xout, y)

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

      implicit none

      include 'global.h'
      integer is, d1, d2, d3, logd2
      double complex x(d1,d2,d3)
      double complex xout(d1,d2,d3)
      double complex y(fftblockpad, d2, 2) 
      integer i, j, k, ii

      logd2 = ilog2(d2)

      do k = 1, d3
        do ii = 0, d1 - fftblock, fftblock
           if (timers_enabled) call timer_start(T_fftcopy)
           do j = 1, d2
              do i = 1, fftblock
                 y(i,j,1) = x(i+ii,j,k)
              enddo
           enddo
           if (timers_enabled) call timer_stop(T_fftcopy)

           if (timers_enabled) call timer_start(T_fftlow)
           call cfftz (is, logd2, d2, y, y(1, 1, 2))
           if (timers_enabled) call timer_stop(T_fftlow)

           if (timers_enabled) call timer_start(T_fftcopy)
           do j = 1, d2
              do i = 1, fftblock
                 xout(i+ii,j,k) = y(i,j,1)
              enddo
           enddo
           if (timers_enabled) call timer_stop(T_fftcopy)
        enddo
      enddo

      return
      end


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

      subroutine cffts3(is, d1, d2, d3, x, xout, y)

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

      implicit none

      include 'global.h'
      integer is, d1, d2, d3, logd3
      double complex x(d1,d2,d3)
      double complex xout(d1,d2,d3)
      double complex y(fftblockpad, d3, 2) 
      integer i, j, k, ii

      logd3 = ilog2(d3)

      do j = 1, d2
        do ii = 0, d1 - fftblock, fftblock
           if (timers_enabled) call timer_start(T_fftcopy)
           do k = 1, d3
              do i = 1, fftblock
                 y(i,k,1) = x(i+ii,j,k)
              enddo
           enddo
           if (timers_enabled) call timer_stop(T_fftcopy)

           if (timers_enabled) call timer_start(T_fftlow)
           call cfftz (is, logd3, d3, y, y(1, 1, 2))
           if (timers_enabled) call timer_stop(T_fftlow)

           if (timers_enabled) call timer_start(T_fftcopy)
           do k = 1, d3
              do i = 1, fftblock
                 xout(i+ii,j,k) = y(i,k,1)
              enddo
           enddo
           if (timers_enabled) call timer_stop(T_fftcopy)
        enddo
      enddo

      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 transpose_x_yz(l1, l2, xin, xout)

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

      implicit none
      include 'global.h'
      integer l1, l2
      double complex xin(ntdivnp), xout(ntdivnp)

      call transpose2_local(dims(1,l1),dims(2, l1)*dims(3, l1),
     >                          xin, xout)

      call transpose2_global(xout, xin)

      call transpose2_finish(dims(1,l1),dims(2, l1)*dims(3, l1),
     >                          xin, xout)

      return
      end


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

      subroutine transpose_xy_z(l1, l2, xin, xout)

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

      implicit none
      include 'global.h'
      integer l1, l2
      double complex xin(ntdivnp), xout(ntdivnp)

      call transpose2_local(dims(1,l1)*dims(2, l1),dims(3, l1),
     >                          xin, xout)
      call transpose2_global(xout, xin)
      call transpose2_finish(dims(1,l1)*dims(2, l1),dims(3, l1),
     >                          xin, xout)

      return
      end


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

      subroutine transpose2_local(n1, n2, xin, xout)

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

      implicit none
      include 'mpinpb.h'
      include 'global.h'
      integer n1, n2
      double complex xin(n1, n2), xout(n2, n1)
      
      double complex z(transblockpad, transblock)

      integer i, j, ii, jj

      if (timers_enabled) call timer_start(T_transxzloc)

c---------------------------------------------------------------------
c If possible, block the transpose for cache memory systems. 
c How much does this help? Example: R8000 Power Challenge (90 MHz)
c Blocked version decreases time spend in this routine 
c from 14 seconds to 5.2 seconds on 8 nodes class A.
c---------------------------------------------------------------------

      if (n1 .lt. transblock .or. n2 .lt. transblock) then
         if (n1 .ge. n2) then 
            do j = 1, n2
               do i = 1, n1
                  xout(j, i) = xin(i, j)
               end do
            end do
         else
            do i = 1, n1
               do j = 1, n2
                  xout(j, i) = xin(i, j)
               end do
            end do
         endif
      else
         do j = 0, n2-1, transblock
            do i = 0, n1-1, transblock
               
c---------------------------------------------------------------------
c Note: compiler should be able to take j+jj out of inner loop
c---------------------------------------------------------------------
               do jj = 1, transblock
                  do ii = 1, transblock
                     z(jj,ii) = xin(i+ii, j+jj)
                  end do
               end do
               
               do ii = 1, transblock
                  do jj = 1, transblock
                     xout(j+jj, i+ii) = z(jj,ii)
                  end do
               end do
               
            end do
         end do
      endif
      if (timers_enabled) call timer_stop(T_transxzloc)

      return
      end


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

      subroutine transpose2_global(xin, xout)

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

      implicit none
      include 'global.h'
      include 'mpinpb.h'
      double complex xin(ntdivnp)
      double complex xout(ntdivnp) 
      integer ierr

!      if (timers_enabled) call synchup()

      if (timers_enabled) call timer_start(T_transxzglo)
      call mpi_alltoall(xin, ntdivnp/np, dc_type,
     >                  xout, ntdivnp/np, dc_type,
     >                  commslice1, ierr)
      if (timers_enabled) call timer_stop(T_transxzglo)

      return
      end


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

      subroutine transpose2_finish(n1, n2, xin, xout)

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

      implicit none
      include 'global.h'
      integer n1, n2, ioff
      double complex xin(n2, n1/np2, 0:np2-1), xout(n2*np2, n1/np2)
      
      integer i, j, p

      if (timers_enabled) call timer_start(T_transxzfin)
      do p = 0, np2-1
         ioff = p*n2
         do j = 1, n1/np2
            do i = 1, n2
               xout(i+ioff, j) = xin(i, j, p)
            end do
         end do
      end do
      if (timers_enabled) call timer_stop(T_transxzfin)

      return
      end


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

      subroutine transpose_x_z(l1, l2, xin, xout)

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

      implicit none
      include 'global.h'
      integer l1, l2
      double complex xin(ntdivnp), xout(ntdivnp)

      call transpose_x_z_local(dims(1,l1),dims(2,l1),dims(3,l1),
     >                         xin, xout)
      call transpose_x_z_global(dims(1,l1),dims(2,l1),dims(3,l1), 
     >                          xout, xin)
      call transpose_x_z_finish(dims(1,l2),dims(2,l2),dims(3,l2), 
     >                          xin, xout)
      return
      end


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

      subroutine transpose_x_z_local(d1, d2, d3, xin, xout)

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

      implicit none
      include 'global.h'
      integer d1, d2, d3
      double complex xin(d1,d2,d3)
      double complex xout(d3,d2,d1)
      integer block1, block3
      integer i, j, k, kk, ii, i1, k1

      double complex buf(transblockpad, maxdim)
      if (timers_enabled) call timer_start(T_transxzloc)
      if (d1 .lt. 32) goto 100
      block3 = d3
      if (block3 .eq. 1)  goto 100
      if (block3 .gt. transblock) block3 = transblock
      block1 = d1
      if (block1*block3 .gt. transblock*transblock) 
     >          block1 = transblock*transblock/block3
c---------------------------------------------------------------------
c blocked transpose
c---------------------------------------------------------------------
      do j = 1, d2
         do kk = 0, d3-block3, block3
            do ii = 0, d1-block1, block1
               
               do k = 1, block3
                  k1 = k + kk
                  do i = 1, block1
                     buf(k, i) = xin(i+ii, j, k1)
                  end do
               end do

               do i = 1, block1
                  i1 = i + ii
                  do k = 1, block3
                     xout(k+kk, j, i1) = buf(k, i)
                  end do
               end do

            end do
         end do
      end do
      goto 200
      

c---------------------------------------------------------------------
c basic transpose
c---------------------------------------------------------------------
 100  continue
      
      do j = 1, d2
         do k = 1, d3
            do i = 1, d1
               xout(k, j, i) = xin(i, j, k)
            end do
         end do
      end do

c---------------------------------------------------------------------
c all done
c---------------------------------------------------------------------
 200  continue

      if (timers_enabled) call timer_stop(T_transxzloc)
      return 
      end


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

      subroutine transpose_x_z_global(d1, d2, d3, xin, xout)

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

      implicit none
      include 'global.h'
      include 'mpinpb.h'
      integer d1, d2, d3
      double complex xin(d3,d2,d1)
      double complex xout(d3,d2,d1) ! not real layout, but right size
      integer ierr

!      if (timers_enabled) call synchup()

c---------------------------------------------------------------------
c do transpose among all  processes with same 1-coord (me1)
c---------------------------------------------------------------------
      if (timers_enabled)call timer_start(T_transxzglo)
      call mpi_alltoall(xin, d1*d2*d3/np2, dc_type,
     >                  xout, d1*d2*d3/np2, dc_type,
     >                  commslice1, ierr)
      if (timers_enabled) call timer_stop(T_transxzglo)
      return
      end
      

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

      subroutine transpose_x_z_finish(d1, d2, d3, xin, xout)

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

      implicit none
      include 'global.h'
      integer d1, d2, d3
      double complex xin(d1/np2, d2, d3, 0:np2-1)
      double complex xout(d1,d2,d3)
      integer i, j, k, p, ioff
      if (timers_enabled) call timer_start(T_transxzfin)
c---------------------------------------------------------------------
c this is the most straightforward way of doing it. the
c calculation in the inner loop doesn't help. 
c      do i = 1, d1/np2
c         do j = 1, d2
c            do k = 1, d3
c               do p = 0, np2-1
c                  ii = i + p*d1/np2
c                  xout(ii, j, k) = xin(i, j, k, p)
c               end do
c            end do
c         end do
c      end do
c---------------------------------------------------------------------

      do p = 0, np2-1
         ioff = p*d1/np2
         do k = 1, d3
            do j = 1, d2
               do i = 1, d1/np2
                  xout(i+ioff, j, k) = xin(i, j, k, p)
               end do
            end do
         end do
      end do
      if (timers_enabled) call timer_stop(T_transxzfin)
      return
      end


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

      subroutine transpose_x_y(l1, l2, xin, xout)

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

      implicit none
      include 'global.h'
      integer l1, l2
      double complex xin(ntdivnp), xout(ntdivnp)

c---------------------------------------------------------------------
c xy transpose is a little tricky, since we don't want
c to touch 3rd axis. But alltoall must involve 3rd axis (most 
c slowly varying) to be efficient. So we do
c (nx, ny/np1, nz/np2) -> (ny/np1, nz/np2, nx) (local)
c (ny/np1, nz/np2, nx) -> ((ny/np1*nz/np2)*np1, nx/np1) (global)
c then local finish. 
c---------------------------------------------------------------------


      call transpose_x_y_local(dims(1,l1),dims(2,l1),dims(3,l1),
     >                         xin, xout)
      call transpose_x_y_global(dims(1,l1),dims(2,l1),dims(3,l1), 
     >                          xout, xin)
      call transpose_x_y_finish(dims(1,l2),dims(2,l2),dims(3,l2), 
     >                          xin, xout)

      return
      end


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

      subroutine transpose_x_y_local(d1, d2, d3, xin, xout)

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

      implicit none
      include 'global.h'
      integer d1, d2, d3
      double complex xin(d1, d2, d3)
      double complex xout(d2, d3, d1)
      integer i, j, k
      if (timers_enabled) call timer_start(T_transxyloc)

      do k = 1, d3
         do i = 1, d1
            do j = 1, d2
               xout(j,k,i)=xin(i,j,k)
            end do
         end do
      end do
      if (timers_enabled) call timer_stop(T_transxyloc)
      return 
      end


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

      subroutine transpose_x_y_global(d1, d2, d3, xin, xout)

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

      implicit none
      include 'global.h'
      include 'mpinpb.h'
      integer d1, d2, d3
c---------------------------------------------------------------------
c array is in form (ny/np1, nz/np2, nx)
c---------------------------------------------------------------------
      double complex xin(d2,d3,d1)
      double complex xout(d2,d3,d1) ! not real layout but right size
      integer ierr

!      if (timers_enabled) call synchup()

c---------------------------------------------------------------------
c do transpose among all processes with same 1-coord (me1)
c---------------------------------------------------------------------
      if (timers_enabled) call timer_start(T_transxyglo)
      call mpi_alltoall(xin, d1*d2*d3/np1, dc_type,
     >                  xout, d1*d2*d3/np1, dc_type,
     >                  commslice2, ierr)
      if (timers_enabled) call timer_stop(T_transxyglo)

      return
      end
      

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

      subroutine transpose_x_y_finish(d1, d2, d3, xin, xout)

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

      implicit none
      include 'global.h'
      integer d1, d2, d3
      double complex xin(d1/np1, d3, d2, 0:np1-1)
      double complex xout(d1,d2,d3)
      integer i, j, k, p, ioff
      if (timers_enabled) call timer_start(T_transxyfin)
c---------------------------------------------------------------------
c this is the most straightforward way of doing it. the
c calculation in the inner loop doesn't help. 
c      do i = 1, d1/np1
c         do j = 1, d2
c            do k = 1, d3
c               do p = 0, np1-1
c                  ii = i + p*d1/np1
c note order is screwy bcz we have (ny/np1, nz/np2, nx) -> (ny, nx/np1, nz/np2)
c                  xout(ii, j, k) = xin(i, k, j, p)
c               end do
c            end do
c         end do
c      end do
c---------------------------------------------------------------------

      do p = 0, np1-1
         ioff = p*d1/np1
         do k = 1, d3
            do j = 1, d2
               do i = 1, d1/np1
                  xout(i+ioff, j, k) = xin(i, k, j, p)
               end do
            end do
         end do
      end do
      if (timers_enabled) call timer_stop(T_transxyfin)
      return
      end


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

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

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

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

      do j=1,1024
         q = mod(j, nx)+1
         if (q .ge. xstart(1) .and. q .le. xend(1)) then
            r = mod(3*j,ny)+1
            if (r .ge. ystart(1) .and. r .le. yend(1)) then
               s = mod(5*j,nz)+1
               if (s .ge. zstart(1) .and. s .le. zend(1)) then
                  chk=chk+u1(q-xstart(1)+1,r-ystart(1)+1,s-zstart(1)+1)
               end if
            end if
         end if
      end do
      chk = chk/ntotal_f

      if (timers_enabled) call timer_start(T_synch)
      call MPI_Reduce(chk, allchk, 1, dc_type, MPI_SUM, 
     >                0, MPI_COMM_WORLD, ierr)      
      if (timers_enabled) call timer_stop(T_synch)
      if (me .eq. 0) then
            write (*, 30) i, allchk
 30         format (' T =',I5,5X,'Checksum =',1P2D22.12)
      endif

c      sums(i) = allchk
c     If we compute the checksum for diagnostic purposes, we let i be
c     negative, so the result will not be stored in an array
      if (i .gt. 0) sums(i) = allchk

      return
      end

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

      subroutine synchup

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

      implicit none
      include 'global.h'
      include 'mpinpb.h'
      integer ierr
      call timer_start(T_synch)
      call mpi_barrier(MPI_COMM_WORLD, ierr)
      call timer_stop(T_synch)
      return
      end


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

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

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

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

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


      class = 'U'

      if (me .ne. 0) return

      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

      call MPI_COMM_SIZE(MPI_COMM_WORLD, size, ierr)
      if (size .ne. np) then
         write(*, 4010) np
         write(*, 4011)
         write(*, 4012)
c---------------------------------------------------------------------
c multiple statements because some Fortran compilers have
c problems with long strings. 
c---------------------------------------------------------------------
 4010    format( ' Warning: benchmark was compiled for ', i5, 
     >           'processors')
 4011    format( ' Must be run on this many processors for official',
     >           ' verification')
 4012    format( ' so memory access is repeatable')
         verified = .false.
      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