      PROGRAM EMBAR
C
C   This is the standard Fortran-77 version of the APP Benchmark 1,
C   the "embarassingly parallel" benchmark.
C   On 64 bit systems, real*8 should be disabled.
C
C   David H. Bailey     January 8, 1991
C
C   M is the Log_2 of the number of complex pairs of uniform (0, 1) random
C   numbers.  MK is the Log_2 of the size of each batch of uniform random
C   numbers.  MK can be set for convenience on a given system, since it does
C   not affect the results.  TIMER is a real*8 function that
C   returns elapsed CPU time in TIMERs.
C
      IMPLICIT real*8 (A-H, O-Z)
      PARAMETER (M = 24, MK = 16, MM = M - MK, NN = 2 ** MM,            
     $  NK = 2 ** MK, NQ = 10, A = 1220703125.d0, S = 271828183.d0)
      DIMENSION X(2*NK), Q(0:NQ-1)
C
      TM0 = TIMER ()
      CALL VRANLC (0, T1, A, X)
C
C   Compute AN = A ^ (2 * NK) (mod 2^46).
C
      T1 = A
C
      DO 100 I = 1, MK + 1
        T2 = RANDLC (T1, T1)
 100  CONTINUE
C
      AN = T1
      TT = S
      GC = 0.d0
C
      DO 110 I = 0, NQ - 1
        Q(I) = 0.d0
 110  CONTINUE
C
C   Each instance of this loop may be performed independently.
C
      DO 150 K = 1, NN
        KK = K - 1
        T1 = S
        T2 = AN
C
C   Find starting seed T1 for this KK.
C
        DO 120 I = 1, 100
          IK = KK / 2
          IF (2 * IK .NE. KK) T3 = RANDLC (T1, T2)
          IF (IK .EQ. 0) GOTO 130
          T3 = RANDLC (T2, T2)
          KK = IK
 120    CONTINUE
C
C   Compute uniform pseudorandom numbers.
C
 130    CALL VRANLC (2 * NK, T1, A, X)
C
C   On a single processor system, the 120 loop and line 130 can be replaced
C   by the following single line.
C
C        CALL VRANLC (2 * NK, TT, A, X)
C
C   Compute Gaussian deviates by acceptance-rejection method and tally counts
C   in concentric square annuli.   This loop is not vectorizable.
C
        DO 140 I = 1, NK
          X1 = 2.d0 * X(2*I-1) - 1.d0
          X2 = 2.d0 * X(2*I) - 1.d0
          T1 = X1 ** 2 + X2 ** 2
          IF (T1 .LE. 1.d0) THEN
            GC = GC + 1.d0
            T2 = SQRT (-2.d0 * LOG (T1) / T1)
            T3 = ABS (X1 * T2)
            T4 = ABS (X2 * T2)
            L = MAX (T3, T4)
            Q(L) = Q(L) + 1.d0
          ENDIF
 140    CONTINUE
C
 150  CONTINUE
C
      TM1 = TIMER ()
      WRITE (6, 1) TM1 - TM0, M, GC, (I, Q(I), I = 0, NQ - 1)
 1    FORMAT ('BENCHMARK 1 RESULTS:'//'CPU TIME =',F10.4/'N = 2^',I5/   
     $  'NO. GAUSSIAN PAIRS =',F15.0/'COUNTS:'/(I3,F15.0))
      STOP
      END
C
      function timer(x)
      real*8 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
      FUNCTION RANDLC (X, A)
C
C   This routine returns a uniform pseudorandom real*8 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 real*8 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 real*8 floating point data.  On 64 bit
C   systems, real*8 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 real*8 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 real*8 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 real*8 floating point
C   data.  On 64 bit systems, real*8 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
