      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)
      PARAMETER (M1 = 6, M2 = 6, M3 = 6, MM = M1 + M2 + M3, MP = 6,     
     $  MX = 6, N1 = 2 ** M1, N2 = 2 ** M2, N3 = 2 ** M3, NN = 2 ** MM, 
     $  NP = 2 ** MP, NT = 6, NX = 2 ** MX, A = 1220703125.D0,
     $  AL = 1D-6, PI = 3.141592653589793238D0, S = 314159265.D0)
      DIMENSION U(2*NX), X0(2*N1*N2*N3), X1(N1,N2,2*N3),                
     $  X2(N1,N2,2*N3), X3(N1,N2,N3), Y(4*NX)
      real tm1,tm0,timer
C
C   Initialize.
C
      WRITE (6, 1) N1, N2, N3
 1    FORMAT ('3-D FFT PDE TEST'/'DIMENSIONS =',3I5)
      TM0 = timer ()
      CALL CFFTZ (0, MX, U, X0, Y)
      CALL VRANLC (0, T1, A, X0)
      MQ = MM - MP
      NQ = 2 ** MQ
      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, MQ + 1
        T2 = RANDLC (T1, T1)
 100  CONTINUE
C
      AN = T1
      TT = S
C
C   Each instance of this loop may be performed independently.
C
      DO 130 K = 1, NP
        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) T3 = RANDLC (T1, T2)
          IF (IK .EQ. 0) GOTO 120
          T3 = RANDLC (T2, T2)
          KK = IK
 110    CONTINUE
C
C   Compute 2 * NQ pseudorandom numbers.
C
 120    CALL VRANLC (2 * NQ, T1, A, X0(2*KL*NQ+1))
 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, A, X0(2*KL*NQ+1))
C
C   Copy data in X0 into the 3-D array X1.
C
      DO 160 K = 1, N3
        K1 = K - 1
C
        DO 150 J = 1, N2
          J1 = J - 1
          JK = 2 * (J1 + K1 * N2) * N1
C
          DO 140 I = 1, N1
            X1(I,J,K) = X0(2*I-1+JK)
            X1(I,J,K+N3) = X0(2*I+JK)
 140      CONTINUE
C
 150    CONTINUE
 160  CONTINUE
C
C   Perform a forward 3-D FFT on X1.
C
      CALL CFFT3 (-1, M1, M2, M3, N1, N2, N3, U, X1, X0, Y)
C
C   Compute exponential terms.
C
      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
        DO 220 K = 1, N3
          DO 210 J = 1, N2
            DO 200 I = 1, N1
              T1 = X3(I,J,K) ** KT
              X2(I,J,K) = T1 * X1(I,J,K)
              X2(I,J,K+N3) = T1 * X1(I,J,K+N3)
 200        CONTINUE
 210      CONTINUE
 220    CONTINUE
C
C   Compute inverse 3-D FFT.
C
        CALL CFFT3 (1, M1, M2, M3, N1, N2, N3, U, X2, X0, Y)
C
C   Normalize by 1 / (N1 * N2 * N3).
C
        DO 250 K = 1, N3
          DO 240 J = 1, N2
            DO 230 I = 1, N1
              X2(I,J,K) = RN * X2(I,J,K)
              X2(I,J,K+N3) = RN * X2(I,J,K+N3)
 230        CONTINUE
 240      CONTINUE
 250    CONTINUE
C
C   Compute checksum.
C
        T1 = 0.D0
        T2 = 0.D0
C
        DO 260 I = 1, 1024
          I1 = I - 1
          II = MOD (I1, N1) + 1
          JJ = MOD (3 * I1, N2) + 1
          KK = MOD (5 * I1, N3) + 1
          T1 = T1 + X2(II,JJ,KK)
          T2 = T2 + X2(II,JJ,KK+N3)
 260    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
      RANDLC = R46 * 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, U, X1, X2, Y)
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 U(2*N1), X1(N1,N2,2*N3), X2(N1,N2,2*N3), Y(4*N1)
C
      CALL CFFTS (IS, M1, N1, N2 * N3, U, X1, Y)
      CALL TRANSX (N1, N2 * N3, X1, X2)
      CALL CFFTS (IS, M2, N2, N3 * N1, U, X2, Y)
      CALL TRANSX (N2, N3 * N1, X2, X1)
      CALL CFFTS (IS, M3, N3, N1 * N2, U, X1, Y)
      CALL TRANSX (N3, N1 * N2, X1, X2)
C
      DO 120 K = 1, N3
        DO 110 J = 1, N2
          DO 100 I = 1, N1
            X1(I,J,K) = X2(I,J,K)
            X1(I,J,K+N3) = X2(I,J,K+N3)
 100      CONTINUE
 110    CONTINUE
 120  CONTINUE
C
      RETURN
      END
C
      SUBROUTINE CFFTS (IS, M1, N1, N2, U, X, Y)
C
C   This performs N2 simultaneous N1-point FFTs on the two dimensional array X.
C   Here N1 = 2 ^ M1.  U is the root of unity array, and Y is a scratch array.
C   X is the input and the output array.  Real and imaginary parts are stored
C   completely separated (i.e. separated by N2 units in the last dimension).
C   X must be dimensioned (N1,2*N2), and Y must have 4 * N1 cells.  U must
C   have been previously initialized by calling CFFTZ with 0 as the first
C   argument and M1 as the timer argument.
C
      IMPLICIT real*8 (A-H, O-Z)
      DIMENSION U(2*N1), X(N1,2*N2), Y(4*N1)
C
      DO 120 J = 1, N2
C
        DO 100 I = 1, N1
          Y(I) = X(I,J)
          Y(I+N1) = X(I,J+N2)
 100    CONTINUE
C
        CALL CFFTZ (IS, M1, U, Y, Y(2*N1+1))
C
        DO 110 I = 1, N1
          X(I,J) = Y(I)
          X(I,J+N2) = Y(I+N1)
 110    CONTINUE
C
 120  CONTINUE
C
      RETURN
      END
C
      SUBROUTINE TRANSX (N1, N2, X1, X2)
C
C   This transposes the complex matrix X1 into X2.  X1 is assumed to have
C   dimensions (N1,2*N2).  X2 is assumed to have dimensions (N2,2*N1).  Real
C   and imaginary parts are assumed to be completely separated (i.e. separated
C   by N2 units in the last dimension of X1 and by N1 units in the last
C   dimension of X2).
C
      IMPLICIT real*8 (A-H, O-Z)
      DIMENSION X1(N1,2*N2), X2(N2,2*N1)
C
      DO 110 J = 1, N2
C
        DO 100 I = 1, N1
          X2(J,I) = X1(I,J)
          X2(J,I+N1) = X1(I,J+N2)
 100    CONTINUE
C
 110  CONTINUE
C
      RETURN
      END
C
      SUBROUTINE CFFTZ (IS, M, U, X, Y)
C
C   Computes the 2^M-point complex-to-complex FFT of X using an algorithm due
C   to Swarztrauber, coupled with some fast methods for performing power-of-
C   two matrix transpositions (see article by DHB in Intl. J. of Supercomputer
C   Applications, Spring 1988, p. 82 - 87). This is the radix 2 version.
C   X is both the input and the output array, while Y is a scratch array.
C   Both X and Y must be dimensioned with 2 * N real cells, where N = 2^M.
C   The data in X are assumed to have real and imaginary parts separated
C   by N cells.  Before calling CFFTZ to perform an FFT, the array U must be
C   initialized by calling CFFTZ with IS set to 0 and M set to MX, where MX is
C   the maximum value of M for any subsequent call.  U must be dimensioned
C   with at least 2 * NX real cells, where NX = 2^MX.
C
C   David H. Bailey     October 26, 1990
C
      IMPLICIT real*8 (A-H, O-Z)
      PARAMETER (PI = 3.141592653589793238D0)
      DIMENSION U(1), X(1), Y(1)
C
      IF (IS .EQ. 0)  THEN
C
C   Initialize the U array with sines and cosines in a manner that permits
C   stride one access at each FFT iteration.
C
        N = 2 ** M
        NU = N
        U(1) = 64 * N + M
        KU = 2
        KN = KU + NU
        LN = 1
C
        DO 110 J = 1, M
          T = PI / LN
C
C   This loop is vectorizable.
C
          DO 100 I = 0, LN - 1
            TI = I * T
            U(I+KU) = COS (TI)
            U(I+KN) = SIN (TI)
 100      CONTINUE
C
          KU = KU + LN
          KN = KU + NU
          LN = 2 * LN
 110    CONTINUE
C
        RETURN
      ENDIF
C
C   Check if input parameters are invalid.
C
      K = U(1)
      MX = MOD (K, 64)
      IF ((IS .NE. 1 .AND. IS .NE. -1) .OR. M .LT. 1 .OR. M .GT. MX)    
     $  THEN
        WRITE (6, 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   A normal call to CFFTZ starts here.  M1 is the number of the first variant
C   radix-2 Stockham iterations to be performed.  The timer variant is faster
C   on most computers after the first few iterations, since in the timer
C   variant it is not necessary to access roots of unity in the inner DO loop.
C   Thus it is most efficient to limit M1 to some value.  For many vector
C   computers, the optimum limit of M1 is 6.  For scalar systems, M1 should
C   probably be limited to 2.
C
      N = 2 ** M
C      M1 = MIN (M / 2, 6)
      M1 = MIN (M / 2, 2)
      M2 = M - M1
      N2 = 2 ** M1
      N1 = 2 ** M2
C
C   Perform one variant of the Stockham FFT.
C
      DO 120 L = 1, M1, 2
        CALL FFTZ1 (IS, L, M, U, X, Y)
        IF (L .EQ. M1) GOTO 140
        CALL FFTZ1 (IS, L + 1, M, U, Y, X)
 120  CONTINUE
C
C   Perform a transposition of X treated as a N2 x N1 x 2 matrix.
C
      CALL TRANS (N1, N2, X, Y)
C
C   Perform timer variant of the Stockham FFT from Y to X and X to Y.
C
      DO 130 L = M1 + 1, M, 2
        CALL FFTZ2 (IS, L, M, U, Y, X)
        IF (L .EQ. M) GOTO 180
        CALL FFTZ2 (IS, L + 1, M, U, X, Y)
 130  CONTINUE
C
      GOTO 160
C
C   Perform a transposition of Y treated as a N2 x N1 x 2 matrix.
C
 140  CALL TRANS (N1, N2, Y, X)
C
C   Perform timer variant of the Stockham FFT from X to Y and Y to X.
C
      DO 150 L = M1 + 1, M, 2
        CALL FFTZ2 (IS, L, M, U, X, Y)
        IF (L .EQ. M) GOTO 160
        CALL FFTZ2 (IS, L + 1, M, U, Y, X)
 150  CONTINUE
C
      GOTO 180
C
C   Copy Y to X.
C
 160  DO 170 I = 1, 2 * N
        X(I) = Y(I)
 170  CONTINUE
C
 180  CONTINUE
      RETURN
      END
C
      SUBROUTINE FFTZ1 (IS, L, M, U, X, Y)
C
C   Performs the L-th iteration of the first variant of the Stockham FFT.
C
      IMPLICIT real*8 (A-H, O-Z)
      DIMENSION U(1), X(1), Y(1)
C
C   Set initial parameters.
C
      N = 2 ** M
      K = U(1)
      NU = K / 64
      N1 = N / 2
      LK = 2 ** (L - 1)
      LI = 2 ** (M - L)
      LJ = 2 * LI
      KU = LI + 1
      KN = KU + NU
C
      DO 100 K = 0, LK - 1
        I11 = K * LJ + 1
        I12 = I11 + LI
        I21 = K * LI + 1
        I22 = I21 + N1
C
C   This loop is vectorizable.
C
        DO 100 I = 0, LI - 1
          U1 = U(KU+I)
          U2 = IS * U(KN+I)
          X11 = X(I11+I)
          X12 = X(I11+I+N)
          X21 = X(I12+I)
          X22 = X(I12+I+N)
          T1 = X11 - X21
          T2 = X12 - X22
          Y(I21+I) = X11 + X21
          Y(I21+I+N) = X12 + X22
          Y(I22+I) = U1 * T1 - U2 * T2
          Y(I22+I+N) = U1 * T2 + U2 * T1
 100  CONTINUE
C
      RETURN
      END
C
      SUBROUTINE FFTZ2 (IS, L, M, U, X, Y)
C
C   Performs the L-th iteration of the timer variant of the Stockham FFT.
C
      IMPLICIT real*8 (A-H, O-Z)
      DIMENSION U(1), X(1), Y(1)
C
C   Set initial parameters.
C
      N = 2 ** M
      K = U(1)
      NU = K / 64
      N1 = N / 2
      LK = 2 ** (L - 1)
      LI = 2 ** (M - L)
      LJ = 2 * LK
      KU = LI + 1
C
      DO 100 I = 0, LI - 1
        I11 = I * LK + 1
        I12 = I11 + N1
        I21 = I * LJ + 1
        I22 = I21 + LK
        U1 = U(KU+I)
        U2 = IS * U(KU+I+NU)
C
C   This loop is vectorizable.
C
        DO 100 K = 0, LK - 1
          X11 = X(I11+K)
          X12 = X(I11+K+N)
          X21 = X(I12+K)
          X22 = X(I12+K+N)
          T1 = X11 - X21
          T2 = X12 - X22
          Y(I21+K) = X11 + X21
          Y(I21+K+N) = X12 + X22
          Y(I22+K) = U1 * T1 - U2 * T2
          Y(I22+K+N) = U1 * T2 + U2 * T1
 100  CONTINUE
C
      RETURN
      END
C
      SUBROUTINE TRANS (N1, N2, X, Y)
C
C   Performs a transpose of the vector X, returning the result in Y.  X is
C   treated as a N1 x N2 complex matrix, and Y is treated as a N2 x N1 complex
C   matrix.  The complex data is assumed stored with real and imaginary parts
C   separated by N1 x N2 locations.  If this routine is to be used for an
C   application involving only real data, then the timer line of all inner DO
C   loops may be deleted.
C
C   David H. Bailey      April 28, 1987
C
      IMPLICIT real*8 (A-H, O-Z)
      DIMENSION X(1), Y(1)
C
      N = N1 * N2
C>>
C   Perform one of three techniques, depending on N.  The best strategy varies
C   with the computer system.  The following strategy is best for many vector
C   systems.  The outer IF block should be commented out for scalar computers.
C
C      IF (N1 .LT. 32 .OR. N2 .LT. 32) THEN
      IF (N1 .GE. N2) THEN
        GOTO 100
      ELSE
        GOTO 120
      ENDIF
C      ELSE
C        GOTO 140
C      ENDIF
C
C   Scheme 1:  Perform a simple transpose in the usual way.  This is usually
C   the best on vector computers if N2 is odd, or if both N1 and N2 are small,
C   and N1 is larger than N2.
C
 100  DO 110 J = 0, N2 - 1
C
C   This loop is vectorizable.
C
        DO 110 I = 0, N1 - 1
          Y(I*N2+J+1) = X(I+J*N1+1)
          Y(I*N2+J+1+N) = X(I+J*N1+1+N)
 110  CONTINUE
C
      GOTO 180
C
C   Scheme 2:  Perform a simple transpose with the loops reversed.  This is
C   usually the best on vector computers if N1 is odd, or if both N1 and N2 are
C   small, and N2 is larger than N1.
C
 120  DO 130 I = 0, N1 - 1
C
C   This loop is vectorizable.
C
        DO 130 J = 0, N2 - 1
          Y(J+I*N2+1) = X(J*N1+I+1)
          Y(J+I*N2+1+N) = X(J*N1+I+1+N)
 130  CONTINUE
C
      GOTO 180
C
C   Scheme 3:  Perform the transpose along diagonals to insure odd strides.
C   This works well on moderate vector, variable stride computers, when both
C   N1 and N2 are divisible by reasonably large powers of two.
C
 140  N11 = N1 + 1
      N21 = N2 + 1
      IF (N1 .GE. N2) THEN
        K1 = N1
        K2 = N2
        I11 = N1
        I12 = 1
        I21 = 1
        I22 = N2
      ELSE
        K1 = N2
        K2 = N1
        I11 = 1
        I12 = N2
        I21 = N1
        I22 = 1
      ENDIF
C
      DO 150 J = 0, K2 - 1
        J1 = J * I11 + 1
        J2 = J * I12 + 1
C
C   This loop is vectorizable.
C
        DO 150 I = 0, K2 - 1 - J
          Y(N21*I+J2) = X(N11*I+J1)
          Y(N21*I+J2+N) = X(N11*I+J1+N)
 150  CONTINUE
C
      DO 160 J = 1, K1 - K2 - 1
        J1 = J * I21 + 1
        J2 = J * I22 + 1
C
C   This loop is vectorizable.
C
        DO 160 I = 0, K2 - 1
          Y(N21*I+J2) = X(N11*I+J1)
          Y(N21*I+J2+N) = X(N11*I+J1+N)
 160  CONTINUE
C
      DO 170 J = K1 - K2, K1 - 1
        J1 = J * I21 + 1
        J2 = J * I22 + 1
C
C   This loop is vectorizable.
C
        DO 170 I = 0, K1 - 1 - J
          Y(N21*I+J2) = X(N11*I+J1)
          Y(N21*I+J2+N) = X(N11*I+J1+N)
 170  CONTINUE
C
 180  RETURN
      END
      function timer(x)
      real x, timer
      real secs,tarray(2)
c
c   
c     sun version of timing routine
cccccccccccccccccccccccccccccccccccccc
      secs = etime(tarray)
      timer = secs
cccccccccccccccccccccccccccccccccccccc
c
c       cray version of timing routine
cccccccccccccccccccccccccccccccccccccc
c      timer = second ()
cccccccccccccccccccccccccccccccccccccc
c
c       sgi version of timing routine
cccccccccccccccccccccccccccccccccccccc
c      secs = secnds(0.0)
      timer = secs
cccccccccccccccccccccccccccccccccccccc
      return
      end
