      PROGRAM FFTPDE
C
C   This is the standard Fortran-77 version of the APP Benchmark 4, the
C   3-D FFT PDE benchmark.
C   On 64 bit systems, double precision should be disabled.
C   Computer specific and tuning notes may be located by searching for C>>.
C>>
C   David H. Bailey     January 8, 1991
C
C   In the following parameter statement, M1, M2 and M3 are the Log_2 of the
C   three dimensions of the 3-D input array.  Set MX = MAX (M1, M2, M3).
C   A is the multiplier of the random number generator (here set to 5^13),
C   and S is the initial seed.  AL is the value of alpha.  NT is the number
C   of iterations.  timer is a double precision function that returns elapsed
C   CPU time in timers.
C
      IMPLICIT real*8 (A-H, O-Z)
c     PARAMETER (M1 = 4, M2 = 4, M3 = 4, MM = M1 + M2 + M3, MP = M3,     
c    $  MX = 4, N1 = 2 ** M1, N2 = 2 ** M2, N3 = 2 ** M3, NN = 2 ** MM, 
c    $  NP = 2 ** MP, NT = 6, NX = 2 ** MX, A = 1220703125.D0,
c    $  AL = 1D-6, PI = 3.141592653589793238D0, S = 314159265.D0)
      PARAMETER (M1 = 8, M2 = 8, M3 = 7, MM = M1 + M2 + M3,
     $  MX = 8, N1 = 2 ** M1, N2 = 2 ** M2, N3 = 2 ** M3, NN = 2 ** MM, 
     $  NT = 6, NX = 2 ** MX, A = 1220703125.D0,
     $  AL = 1D-6, PI = 3.141592653589793238D0, S = 314159265.D0)
c     DIMENSION U(2*NX), X0(2*N1*N2*N3), X1(N1,N2,2*N3),                
c    $  X2(N1,N2,2*N3), X3(N1,N2,N3), Y(4*NX)
      DIMENSION X1real(N1,N2,N3), X2real(N1,N2,N3), X3(N1,N2,N3)
      DIMENSION X1imag(N1,N2,N3), X2imag(N1,N2,N3)
      real*8 tm1,tm0,timer
c
c for some reason SPARC doesnt like passing double prec const to vranlc
c make it a variable.
CAPR$ PARTITION X3 (*,*,SHRUNKBLOCK)
CAPR$ PARTITION X1REAL (*,*,SHRUNKBLOCK)
CAPR$ PARTITION X1IMAG (*,*,SHRUNKBLOCK)
CAPR$ PARTITION X2REAL (*,*,SHRUNKBLOCK)
CAPR$ PARTITION X2IMAG (*,*,SHRUNKBLOCK)
      open (unit=6,file='fft.out',status='unknown')
      aa=a
C
C   Initialize.
C
      WRITE (6, 1) N1, N2, N3
 1    FORMAT ('3-D FFT PDE TEST'/'DIMENSIONS =',3I5)
      TM0 = timer ()
c     CALL CFFTZ (0, MX, U, X0, Y)
      CALL VRANLC (0, T1, Aa, X0)
      RN = 1.D0 / NN
      AP = - 4.D0 * AL * PI ** 2
      N12 = N1 / 2
      N22 = N2 / 2
      N32 = N3 / 2
C
C   Compute AN = A ^ (2 * NQ) (mod 2^46).
C
      T1 = A
C
      DO 100 I = 1, M1+M2+1
        T1 = RANDLC (T1, T1)
 100  CONTINUE
C
      AN = T1
      TT = S
C
C   Each instance of this loop may be performed independently.
C
cgene since x1real and x1imag are partitioned on the third dimension at
cgene the top of this program, and there is only one route to this loop,
cgene no preloop or postloop comm is required for <:,1,1~1>
c
CAPR$ DO PAR ON X1REAL<:, :, 1~1>
CAPR$ IGNORE ALL INHIBITORS
CAPR$ IGNORE ALL COM
      DO 130 K = 1, N3
        KK = K - 1
        KL = KK
        T1 = S
        T2 = AN
C
C   Find starting seed T1 for this KK using the binary rule for exponentiation.
C
        DO 110 I = 1, 100
          IK = KK / 2
          IF (2 * IK .NE. KK) T2 = RANDLC (T1, T2)
          IF (IK .EQ. 0) GOTO 120
          T2 = RANDLC (T2, T2)
          KK = IK
 110    CONTINUE
C
C   Compute 2 * NQ pseudorandom numbers.
C
 120   continue
        CALL VRANLC (N1*N2, T2, aa, x1real(1,1,k))
        CALL VRANLC (N1*N2, T2, aa, x1imag(1,1,k))
 130  CONTINUE
C
C   On a single processor system, the 110 loop and line 120 can be replaced
C   by the following single line.
C
C120     CALL VRANLC (2 * NQ, TT, Aa, X0(2*KL*NQ+1))
C
C   Copy data in X0 into the 3-D array X1.
C
c     DO 160 K = 1, N3
c       K1 = K - 1
c
c       DO 150 J = 1, N2
c         J1 = J - 1
c         JK = 2 * (J1 + K1 * N2) * N1
c
c         DO 140 I = 1, N1
c           X1(I,J,K) = X0(2*I-1+JK)
c           X1(I,J,K+N3) = X0(2*I+JK)
c140      CONTINUE
C
c150    CONTINUE
c160  CONTINUE
C
C   Perform a forward 3-D FFT on X1.
C
      CALL CFFT3 (-1, M1, M2, M3, N1, N2, N3, X1real,x1imag)
C
C   Compute exponential terms.
C
CAPR$ DO PAR ON X3<:, :, 1~1>
      DO 190 K = 1, N3
        K1 = K - 1
        IF (K .GT. N32) K1 = K1 - N3
C
        DO 180 J = 1, N2
          J1 = J - 1
          IF (J .GT. N22) J1 = J1 - N2
          JK = J1 ** 2 + K1 ** 2
C
          DO 170 I = 1, N1
            I1 = I - 1
            IF (I .GT. N12) I1 = I1 - N1
            X3(I,J,K) = EXP (AP * (I1 ** 2 + JK))
 170      CONTINUE
C
 180    CONTINUE
 190  CONTINUE
C
C   Perform the following for KT = 1, ..., NT.
C
      DO 270 KT = 1, NT
C
C   Multiply by the exponential term raised to the KT power.
C
CAPR$ PARTITION X1REAL (*,*,SHRUNKBLOCK)
CAPR$ PARTITION X1IMAG (*,*,SHRUNKBLOCK)
CAPR$ PARTITION X2REAL (*,*,SHRUNKBLOCK)
CAPR$ PARTITION X2IMAG (*,*,SHRUNKBLOCK)
CAPR$ DO PAR ON X3<:, :, 1~1>
CAPR$ IGNORE ALL COM
        DO 220 K = 1, N3
          DO 210 J = 1, N2
            DO 200 I = 1, N1
              T1 = X3(I,J,K) ** KT
              X2real(I,J,K) = T1 * X1real(I,J,K)
              X2imag(I,J,K) = T1 * X1imag(I,J,K)
 200        CONTINUE
 210      CONTINUE
 220    CONTINUE
C
C   Compute inverse 3-D FFT.
C
        CALL CFFT3 (1, M1, M2, M3, N1, N2, N3, x2real,x2imag)
C
C   Normalize by 1 / (N1 * N2 * N3).
C
CAPR$ PARTITION X2REAL (*,*,SHRUNKBLOCK)
CAPR$ PARTITION X2IMAG (*,*,SHRUNKBLOCK)
CAPR$ DO PAR ON X2REAL<:, :, 1~1>
CAPR$ IGNORE ALL COM
        DO 250 K = 1, N3
          DO 240 J = 1, N2
            DO 230 I = 1, N1
              X2real(I,J,K) = RN * X2real(I,J,K)
              X2imag(I,J,K) = RN * X2imag(I,J,K)
 230        CONTINUE
 240      CONTINUE
 250    CONTINUE
C
C   Compute checksum.
C
        T1 = 0.D0
        T2 = 0.D0
C
CAPR$ DO PAR ON X2REAL<:, :, 1~1>
CAPR$ IGNORE ALL COM
        DO 260 I = 1, N3
        DO 260 J = 1, 1024/N3
		  I1 = J - 1 + (I-1) * (1024/N3)
          JJ = MOD (3 * I1, N2) + 1
          KK = MOD (5 * I1, N3) + 1
          T1 = T1 + X2real(JJ,KK,I)
          T2 = T2 + X2imag(JJ,KK,I)
 260    CONTINUE
c       DO 260 I = 1, 1024
c         I1 = I - 1
c         II = MOD (I1, N1) + 1
c         JJ = MOD (3 * I1, N2) + 1
c         KK = MOD (5 * I1, N3) + 1
c         T1 = T1 + X2real(II,JJ,KK)
c         T2 = T2 + X2imag(II,JJ,KK)
c260    CONTINUE
C
        WRITE (6, 2) KT, T1, T2
 2      FORMAT ('T =',I5,5X,'CHECKSUM =',1P2D22.12)
 270  CONTINUE
C
      TM1 = timer ()
      TM = TM1 - TM0
      WRITE (6, 3) TM
 3    FORMAT ('TEST COMPLETED'/'CPU TIME =',F12.6)
      STOP
      END
C
      FUNCTION RANDLC (X, A)
C
C   This routine returns a uniform pseudorandom double precision number in the
C   range (0, 1) by using the linear congruential generator
C
C   x_{k+1} = a x_k  (mod 2^46)
C
C   where 0 < x_k < 2^46 and 0 < a < 2^46.  This scheme generates 2^44 numbers
C   before repeating.  The argument A is the same as 'a' in the above formula,
C   and X is the same as x_0.  A and X must be odd double precision integers
C   in the range (1, 2^46).  The returned value RANDLC is normalized to be
C   between 0 and 1, i.e. RANDLC = 2^(-46) * x_1.  X is updated to contain
C   the new seed x_1, so that subsequent calls to RANDLC using the same
C   arguments will generate a continuous sequence.
C
C   This routine should produce the same results on any computer with at least
C   48 mantissa bits in double precision floating point data.  On 64 bit
C   systems, double precision should be disabled.
C
C   David H. Bailey     October 26, 1990
C
      IMPLICIT real*8 (A-H, O-Z)
      SAVE KS, R23, R46, T23, T46
      DATA KS/0/
C
C   If this is the first call to RANDLC, compute R23 = 2 ^ -23, R46 = 2 ^ -46,
C   T23 = 2 ^ 23, and T46 = 2 ^ 46.  These are computed in loops, rather than
C   by merely using the ** operator, in order to insure that the results are
C   exact on all systems.  This code assumes that 0.5D0 is represented exactly.
C
      IF (KS .EQ. 0) THEN
        R23 = 1.D0
        R46 = 1.D0
        T23 = 1.D0
        T46 = 1.D0
C
        DO 100 I = 1, 23
          R23 = 0.5D0 * R23
          T23 = 2.D0 * T23
 100    CONTINUE
C
        DO 110 I = 1, 46
          R46 = 0.5D0 * R46
          T46 = 2.D0 * T46
 110    CONTINUE
C
        KS = 1
      ENDIF
C
C   Break A into two parts such that A = 2^23 * A1 + A2.
C
      T1 = R23 * A
      A1 = AINT (T1)
      A2 = A - T23 * A1
C
C   Break X into two parts such that X = 2^23 * X1 + X2, compute
C   Z = A1 * X2 + A2 * X1  (mod 2^23), and then
C   X = 2^23 * Z + A2 * X2  (mod 2^46).
C
      T1 = R23 * X
      X1 = AINT (T1)
      X2 = X - T23 * X1
      T1 = A1 * X2 + A2 * X1
      T2 = AINT (R23 * T1)
      Z = T1 - T23 * T2
      T3 = T23 * Z + A2 * X2
      T4 = AINT (R46 * T3)
      X = T3 - T46 * T4
c     RANDLC = R46 * X
      RANDLC = X
C
      RETURN
      END
C
      SUBROUTINE VRANLC (N, X, A, Y)
C
C   This routine generates N uniform pseudorandom double precision numbers in
C   the range (0, 1) by using the linear congruential generator
C
C   x_{k+1} = a x_k  (mod 2^46)
C
C   where 0 < x_k < 2^46 and 0 < a < 2^46.  This scheme generates 2^44 numbers
C   before repeating.  The argument A is the same as 'a' in the above formula,
C   and X is the same as x_0.  A and X must be odd double precision integers
C   in the range (1, 2^46).  The N results are placed in Y and are normalized
C   to be between 0 and 1.  X is updated to contain the new seed, so that
C   subsequent calls to VRANLC using the same arguments will generate a
C   continuous sequence.  If N is zero, only initialization is performed, and
C   the variables X, A and Y are ignored.
C
C   This routine is the standard version designed for scalar or RISC systems.
C   However, it should produce the same results on any single processor
C   computer with at least 48 mantissa bits in double precision floating point
C   data.  On 64 bit systems, double precision should be disabled.
C
C   David H. Bailey     October 26, 1990
C
      IMPLICIT real*8 (A-H, O-Z)
      DIMENSION Y(N)
      SAVE KS, R23, R46, T23, T46
      DATA KS/0/
C
C   If this is the first call to VRANLC, compute R23 = 2 ^ -23, R46 = 2 ^ -46,
C   T23 = 2 ^ 23, and T46 = 2 ^ 46.  See comments in RANDLC.
C
      IF (KS .EQ. 0 .OR. N .EQ. 0) THEN
        R23 = 1.D0
        R46 = 1.D0
        T23 = 1.D0
        T46 = 1.D0
C
        DO 100 I = 1, 23
          R23 = 0.5D0 * R23
          T23 = 2.D0 * T23
 100    CONTINUE
C
        DO 110 I = 1, 46
          R46 = 0.5D0 * R46
          T46 = 2.D0 * T46
 110    CONTINUE
C
        KS = 1
        IF (N .EQ. 0) RETURN
      ENDIF
C
C   Break A into two parts such that A = 2^23 * A1 + A2.
C
      T1 = R23 * A
      A1 = AINT (T1)
      A2 = A - T23 * A1
C
C   Generate N results.   This loop is not vectorizable.
C
      DO 120 I = 1, N
C
C   Break X into two parts such that X = 2^23 * X1 + X2, compute
C   Z = A1 * X2 + A2 * X1  (mod 2^23), and then
C   X = 2^23 * Z + A2 * X2  (mod 2^46).
C
        T1 = R23 * X
        X1 = AINT (T1)
        X2 = X - T23 * X1
        T1 = A1 * X2 + A2 * X1
        T2 = AINT (R23 * T1)
        Z = T1 - T23 * T2
        T3 = T23 * Z + A2 * X2
        T4 = AINT (R46 * T3)
        X = T3 - T46 * T4
        Y(I) = R46 * X
 120  CONTINUE
C
      RETURN
      END
C
      SUBROUTINE CFFT3 (IS, M1, M2, M3, N1, N2, N3, X1real,X1imag)
C
C   This performs a 3-D complex-to-complex FFT on the array X1, which is
C   assumed to have dimensions (N1,N2,2*N3).  It is assumed that N1 = 2 ^ M1,
C   N2 = 2 ^ M2 and N3 = 2 ^ M3.  Real and imaginary parts are stored
C   completely separated (i.e. separated by N3 units in the last dimension of
C   X1).  IS is the sign of the transform, either 1 or -1.  U is the root of
C   unity array, which must have been previously initialized by calling CFFTZ
C   with 0 as the first argument and MM = MAX (M1,M2,M3) as the timer
C   argument.  X2 is used as a scratch array and must be the same size as X1.
C   Y is a scratch array of size 4 * 2 ^ MM.
C
C   David H. Bailey     October 26, 1990
C
      IMPLICIT real*8 (A-H, O-Z)
      DIMENSION X1real(N1,N2,N3),X1imag(N1,N2,N3)
	  parameter(mxsz=8192)
	  complex z(mxsz),w(mxsz/2)
	  logical inverse
	  real timer
	  if(n1.gt.mxsz.or.n2.gt.mxsz.or.n3.gt.mxsz)then
		print *,'Error in CFFT3. Maximum size allowed=',mxsz
		stop 998
	  endif
c	  print *,'cfft3'
c	  print *,'is=',is,' m1=',m1,' m2=',m2,' m3=',m3
c	  print *,' n1=',n1,' n2=',n2,' n3=',n3
c	  print *,'time=',timer()
c	  print *,'input x1 real ='
c	  do k=1,n3
c		do j=1,n2
c			print '(4e12.5)',(x1(i,j,k),i=1,n1)
c		end do
c	  end do
c	  print *,'input x1 imaj ='
c	  do k=n3+1,2*n3
c		do j=1,n2
c			print '(4e12.5)',(x1(i,j,k),i=1,n1)
c		end do
c	  end do
	  if(is.eq.1)then 
		inverse=.false.
	  else
		inverse=.true.
	  endif
	  call tabfft(w,n1)
c	  print *,'first w:'
c	  print '(4e12.5)',(w(i),i=1,n1/2)
CAPR$ PARTITION X1IMAG (*,*,SHRUNKBLOCK)
CAPR$ PARTITION X1REAL (*,*,SHRUNKBLOCK)
CAPR$ DO PAR ON X1REAL<:, :, 1~1>
CAPR$ IGNORE ALL INHIBITORS
CAPR$ IGNORE ALL COM
	  do k=1,n3
		do j=1,n2
			do i=1,n1
				z(i)=cmplx(x1real(i,j,k),x1imag(i,j,k))
			end do
			call fft(z,inverse,w,n1,m1)
			do i=1,n1
				x1real(i,j,k)=real(z(i))
				x1imag(i,j,k)=aimag(z(i))
			end do
		end do
	  end do
c	  print *,'Time after first stage=',timer()
	  call tabfft(w,n2)
c	  print *,'second w:'
c	  print '(4e12.5)',(w(i),i=1,n2/2)
CAPR$ DO PAR ON X1REAL<:, :, 1~1>
CAPR$ IGNORE ALL INHIBITORS
CAPR$ IGNORE ALL COM
	  do k=1,n3
		do j=1,n1
			do i=1,n2
				z(i)=cmplx(x1real(j,i,k),x1imag(j,i,k))
			end do
			call fft(z,inverse,w,n2,m2)
			do i=1,n2
				x1real(j,i,k)=real(z(i))
				x1imag(j,i,k)=aimag(z(i))
			end do
		end do
	  end do
c	  print *,'Time after second stage=',timer()
	  call tabfft(w,n3)
c	  print *,'second w:'
c	  print '(4e12.5)',(w(i),i=1,n3/2)
CAPR$ PARTITION X1IMAG (*,SHRUNKBLOCK,*)
CAPR$ PARTITION X1REAL (*,SHRUNKBLOCK,*)
CAPR$ DO PAR ON X1REAL<:, 1~1, :>
CAPR$ IGNORE ALL INHIBITORS
CAPR$ IGNORE ALL COM
	  do k=1,n2
		do j=1,n1
			do i=1,n3
				z(i)=cmplx(x1real(j,k,i),x1imag(j,k,i))
			end do
			call fft(z,inverse,w,n3,m3)
			do i=1,n3
				x1real(j,k,i)=real(z(i))
				x1imag(j,k,i)=aimag(z(i))
			end do
		end do
	  end do
c	  print *,'Time after third stage=',timer()
c	  print *,'output x1 real ='
c	  do k=1,n3
c		do j=1,n2
c			print '(4e12.5)',(x1(i,j,k),i=1,n1)
c		end do
c	  end do
c	  print *,'output x1 imaj ='
c	  do k=n3+1,2*n3
c		do j=1,n2
c			print '(4e12.5)',(x1(i,j,k),i=1,n1)
c		end do
c	  end do
C
      RETURN
      END

      subroutine fft( x,inverse,w,nx,lnx )
      complex x(nx),
     $        w(nx/2)
      logical  inverse

c   fft
c   (note that variable iter = log(base 2) of nx

      real    rnx
      complex t,
     $        wk

      integer itab(12)
      integer iter

      data itab / 1, 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024, 2048 /
      iter=lnx

      nxp2 = nx
      rnx  = 1./float(nx)

      do 30 it = 1,iter
         n = 1
         nxp = nxp2
         nxp2 = nxp/2
         do 20 m = 1,nxp2
			if( inverse ) then
				wk = conjg(w(n))
			else
				wk = w(n)
			end if
            do 10 mxp = nxp,nx,nxp
               j1 = mxp - nxp + m
               j2 = j1 + nxp2
                  t = x(j1) - x(j2)
                  x(j1) = x(j1) + x(j2)
                  x(j2) = t*wk
c	  print *,'it=',it,' m=',m,' mxp=',mxp,' j1=',j1,' j2=',j2
  10        continue
            n = n +itab(it)
  20     continue
  30  continue

      n2 = nx/2
      n1 = nx - 1
      j = 1

      do i = 1,n1
		if( i .lt. j ) then
			t = x(j)
			x(j) = x(i)
			x(i) = t
		end if
        k = n2
  42    if( k .lt. j ) then
            j = j - k
            k = k/2
            go to 42
		end if
         j = j + k
      end do

      if( inverse ) then
         do k = 1,nx
            x(k) = rnx*x(k)
         end do
      end if

      end
      subroutine tabfft(w,nx)
c----------------------------------------------------------------------
c     calculate some arrays needed for the fft filter
c----------------------------------------------------------------------
      complex w(nx/2)

      do 1 i=1,nx/2
      arg = (i-1)*3.141592653589793/float(nx/2)
      w(i) = cmplx(cos(arg),sin(arg))
   1  continue
      return
      end
	  real*8 function timer()
	  real*8 second
	  timer = second()
	  end
