C     ==========================================================
C     DEBUG VARIABLES
C     ==========================================================

      MODULE DEBUG_VARS

       INTEGER, PARAMETER                     :: LEVEL0 = 0
       INTEGER, PARAMETER                     :: LEVEL1 = 1
       INTEGER, PARAMETER                     :: LEVEL2 = 2 
       INTEGER, PARAMETER                     :: DEBUG  = LEVEL0

      END MODULE DEBUG_VARS

C     ==========================================================
C     LINEAR CONGRUENTIAL GENERATOR MODULE : INTERNAL DATA
C     ==========================================================
      MODULE LCG_PRIVATE_VARIABLES

C     **********************************************************
C     THESE ARE THE PARAMETERS(MULTIPLIER AND MODULUS) OF THE
C     GENERATOR. I USED THE VALUES PROVIDED BY L'ECUYER PAPER 
C     FOR A SIMPLE 32-BITS GENERATOR
C     **********************************************************

      INTEGER, PARAMETER     :: M = 2147483563 
      INTEGER, PARAMETER     :: A = 40014

C     **********************************************************
C     THESE ARE THE PARAMETERS FOR SPLITTING THE SEQUENCE AMONG
C     THE PROCESSORS. SUBSEQUENCES ARE 2**LAG VALUES APART
C     **********************************************************

      INTEGER                :: LAG     
      INTEGER                :: A_2_LAG

      END MODULE LCG_PRIVATE_VARIABLES


C     ==========================================================
C     NBODY PROBLEM
C     ==========================================================
      PROGRAM NBODY

      USE DEBUG_VARS
      IMPLICIT NONE

      INTERFACE

        SUBROUTINE TIMER(T)
	      DOUBLE PRECISION, INTENT(OUT) :: T
        END SUBROUTINE TIMER

        EXTRINSIC(HPF_LOCAL) SUBROUTINE L_TIMER(R_TIME,I_TIME) 
              REAL, DIMENSION(:), INTENT(IN)  :: I_TIME
              REAL, DIMENSION(:), INTENT(OUT) :: R_TIME
!HPF$         DISTRIBUTE *(BLOCK)             :: I_TIME
!HPF$         ALIGN R_TIME WITH I_TIME
        END SUBROUTINE L_TIMER

        SUBROUTINE LCG_INIT(S,SEED,GAP)
              INTEGER, DIMENSION(:), INTENT(OUT)  :: S
              INTEGER, OPTIONAL, INTENT(IN)       :: SEED
              INTEGER, OPTIONAL, INTENT(IN)       :: GAP
        END SUBROUTINE LCG_INIT

        EXTRINSIC(HPF_LOCAL) SUBROUTINE LCG_IRANDOM(S,V)
              INTEGER, DIMENSION(:), INTENT(OUT)  :: S
!HPF$         DISTRIBUTE *(BLOCK)                 :: S
              INTEGER, DIMENSION(:), INTENT(OUT)  :: V
!HPF$         DISTRIBUTE *(BLOCK)                 :: V
        END SUBROUTINE LCG_IRANDOM

        EXTRINSIC(HPF_LOCAL) SUBROUTINE LCG_RRANDOM(S,V)
              INTEGER, DIMENSION(:), INTENT(OUT)  :: S
!HPF$         DISTRIBUTE *(BLOCK)                 :: S
              REAL, DIMENSION(:), INTENT(OUT)     :: V
!HPF$         DISTRIBUTE *(BLOCK)                 :: V
        END SUBROUTINE LCG_RRANDOM

      END INTERFACE   
   
C     ==========================================================
C     VARIABLES , ALIGNED WRT THE PROBLEM'S SPACE
C     ==========================================================

      integer, parameter :: MM=0, PX=1, PY=2, PZ=3, FX=4, FY=5, FZ=6
!     MM:mass, PX:position in x coordinate, PY:position in y coordinate
!     FX: force in x coordinate, FY: force in y coordinate

      real, parameter    :: G = 1.0

      integer, parameter          :: NN = 1024 
      real, dimension(0:NN-1)     :: x,t,dx,dy,dz,tx,ty,tz,sq,dist,fac
      real, dimension(0:NN-1,0:6) :: p,q

!   Specify the HPf directives

!HPF$ DISTRIBUTE x(block)
!HPF$ ALIGN (:)   WITH x(:) :: t,dx,dy,dz,tx,ty,tz,sq,dist,fac
!HPF$ ALIGN (:,*) WITH x(:) :: p,q

      integer           :: i,j,k,npts
      logical           :: odd, IsOdd

      INTEGER, ALLOCATABLE, DIMENSION(:) :: RAND_STATE
!HPF$ DISTRIBUTE RAND_STATE(BLOCK)

      INTEGER, ALLOCATABLE, DIMENSION(:) :: TEST_RAND
!HPF$ DISTRIBUTE TEST_RAND(BLOCK)
      REAL, ALLOCATABLE, DIMENSION(:) :: RTEST_RAND
!HPF$ DISTRIBUTE RTEST_RAND(BLOCK)

C     ==========================================================
C     TIMING VARIABLES
C     ==========================================================

      DOUBLE PRECISION                       :: START_TIME
      DOUBLE PRECISION                       :: END_TIME

      REAL, ALLOCATABLE,DIMENSION(:)         :: ti, tf
!HPF$ DISTRIBUTE ti(BLOCK) 
!HPF$ DISTRIBUTE tf(BLOCK) 

C     ==========================================================
C     INITIALIZATION
C     ==========================================================

        print *,'********************************************'
        print *,'              NBODY PROBLEM'
        print *,'NUMBER OF POINTS ' , NN
        print *,'NUMBER OF PROCESSORS ', NUMBER_OF_PROCESSORS()
 
        ALLOCATE(ti(NUMBER_OF_PROCESSORS()))
        ALLOCATE(tf(NUMBER_OF_PROCESSORS()))

!       Start Timer

        CALL TIMER(START_TIME)
        call l_timer(tf, ti)

        ALLOCATE(RAND_STATE(NUMBER_OF_PROCESSORS()))
        CALL LCG_INIT(RAND_STATE)

     
        CALL LCG_RRANDOM(RAND_STATE,x)
        p(:,MM) = x
        CALL LCG_RRANDOM(RAND_STATE,x)
        p(:,PX) = x
        CALL LCG_RRANDOM(RAND_STATE,x)
        p(:,PY) = x
        CALL LCG_RRANDOM(RAND_STATE,x)
        p(:,PZ) = x

        p(:,FX) = 0.0
        p(:,FY) = 0.0
        p(:,FZ) = 0.0

! Stop timer
        call l_timer(ti, tf)
        CALL TIMER(END_TIME)

        print *,'time taken for initialization is ',ti,'seconds'
        print *,'INITIALIZATION ELAPSED TIME ', END_TIME-START_TIME

C     ==========================================================
C     PROCESSING
C     ==========================================================

!       Start Timer
        ti = 0.0
        tf = 0.0
        CALL TIMER(START_TIME)
        call l_timer(tf, ti)

        npts = NN

        q = p
        k = npts/2 - 1
        do i = 0, k - 1
          q = CSHIFT(q,DIM=1,SHIFT=1)
          dx(0:npts-1) = p(0:npts-1,PX) - q(0:npts-1,PX)
          dy(0:npts-1) = p(0:npts-1,PY) - q(0:npts-1,PY)
          dz(0:npts-1) = p(0:npts-1,PZ) - q(0:npts-1,PZ)
          sq(0:npts-1) = dx(0:npts-1)**2+dy(0:npts-1)**2+
     &        dz(0:npts-1)**2
          dist(0:npts-1) = SQRT(sq(0:npts-1))
          fac(0:npts-1) = p(0:npts-1,MM) * q(0:npts-1,MM) 
     &                   / (dist(0:npts-1) * sq(0:npts-1))
          tx(0:npts-1) = fac(0:npts-1) * dx(0:npts-1)
          ty(0:npts-1) = fac(0:npts-1) * dy(0:npts-1)
          tz(0:npts-1) = fac(0:npts-1) * dz(0:npts-1)
          p(0:npts-1,FX) = p(0:npts-1,FX)-tx(0:npts-1)
          q(0:npts-1,FX) = q(0:npts-1,FX)+tx(0:npts-1)
          p(0:npts-1,FY) = p(0:npts-1,FY)-ty(0:npts-1)
          q(0:npts-1,FY) = q(0:npts-1,FY)+ty(0:npts-1)
          p(0:npts-1,FZ) = p(0:npts-1,FZ)-tz(0:npts-1)
          q(0:npts-1,FZ) = q(0:npts-1,FZ)+tz(0:npts-1)
        end do

        q = CSHIFT(q,DIM=1,SHIFT=1)
        k = npts/2
        dx(0:k-1) = p(0:k-1,PX) - q(0:k-1,PX)
        dy(0:k-1) = p(0:k-1,PY) - q(0:k-1,PY)
        dz(0:k-1) = p(0:k-1,PZ) - q(0:k-1,PZ)
        sq(0:k-1) = dx(0:k-1)**2+dy(0:k-1)**2+dz(0:k-1)**2
        dist(0:k-1) = SQRT(sq(0:k-1))
        fac(0:k-1) = p(0:k-1,MM) * q(0:k-1,MM) 
     &              / (dist(0:k-1) * sq(0:k-1))
        tx(0:k-1) = fac(0:k-1) * dx(0:k-1)
        ty(0:k-1) = fac(0:k-1) * dy(0:k-1)
        tz(0:k-1) = fac(0:k-1) * dz(0:k-1)
        p(0:k-1,FX) = p(0:k-1,FX)-tx(0:k-1)
        q(0:k-1,FX) = q(0:k-1,FX)+tx(0:k-1)
        p(0:k-1,FY) = p(0:k-1,FY)-ty(0:k-1)
        q(0:k-1,FY) = q(0:k-1,FY)+ty(0:k-1)
        p(0:k-1,FZ) = p(0:k-1,FZ)-tz(0:k-1)
        q(0:k-1,FZ) = q(0:k-1,FZ)+tz(0:k-1)

        q = CSHIFT(q,DIM=1,SHIFT=npts/2)
        p(:,FX) = p(:,FX) + q(:,FX)
        p(:,FY) = p(:,FY) + q(:,FY)
        p(:,FZ) = p(:,FZ) + q(:,FZ)
        
! Stop timer
        call l_timer(ti, tf)
        CALL TIMER(END_TIME)

        print *,'time taken is ',ti,'seconds'
        print *,'ELAPSED TIME PROCESSING ', END_TIME-START_TIME
 
        print *, 'FX', p(0:6,FX)
        print *, 'FY', p(0:6,FY)
        print *, 'FZ', p(0:6,FZ)

      STOP
      END PROGRAM NBODY

C     ==========================================================
C     TIMER SUBROUTINE
C        RETURN NUMBER OF ELPASED SECONDS
C     ==========================================================

      SUBROUTINE TIMER(T)
        DOUBLE PRECISION, INTENT(OUT) :: T
      
        INTEGER, DIMENSION(8) :: VAL

        CALL DATE_AND_TIME(VALUES=VAL) 
        T = (DBLE(VAL(8)) * 1D-3) + (DBLE(VAL(7))) +
     &      (DBLE(VAL(6)) * 60D0) + (DBLE(VAL(5)) * 3600D0)

      RETURN
      END SUBROUTINE TIMER

C     ==========================================================
C     L_TIMER SUBROUTINE
C        RETURN NUMBER OF SECONDS IN EACH PROCESSOR
C     ==========================================================

      extrinsic(hpf_local) subroutine l_timer(r_time,i_time)
        real, dimension(:), intent(in)  :: i_time
        real, dimension(:), intent(out) :: r_time
        
        r_time(lbound(r_time, 1)) =
     &  secnds(i_time(lbound(i_time, 1)))

      end

C     ==========================================================
C     MULTIMOD
C       THIS CALCULATES A*S MOD M FOR ARBITRARY A,S,M  WITHOUT
C       INTEGER OVERFLOW (E.G. ALL INTERMEDIATE RESULTS ARE WELL
C       REPRESENTED INTEGERS). THIS FUNCTION USES THE RUSSIAN
C       PEASANT METHOD DESCRIBED IN L'ECUYER PAPER
C     ==========================================================

      INTEGER FUNCTION MULTIMOD (A,S,M)
        INTEGER,INTENT(IN) :: A,S,M

        REAL    :: LARGE_POSITIVE
        INTEGER :: H
        INTEGER :: AA,SS
        INTEGER :: Q,R,K,P

        LARGE_POSITIVE =REAL(HUGE(M))+1
        H = INT(SQRT(LARGE_POSITIVE/2))

        AA = A ; SS = S
        P = -M
        RUSSIAN_PEASANT: DO


          IF (NOT(AA .GT. H)) EXIT RUSSIAN_PEASANT

          IF (MOD(AA,2).EQ.1) THEN
             P = P+SS
             IF (P .GT. 0) THEN
                P = P-M
             END IF
          END IF

          AA = AA/2
          SS = (SS-M)+SS
          IF (SS .LT. 0) THEN
              SS = SS+M
          END IF
        END DO RUSSIAN_PEASANT
  
        Q = M / AA
        K = SS / Q
        SS = (AA*(SS-(K*Q))) - (K*(M-(Q*AA)))
        DO
          IF (NOT(SS .LT. 0)) EXIT
          SS = SS+M
        END DO     

        P = P+SS
        IF (P .LT. 0) THEN
            P = P+M 
        END IF

        MULTIMOD = P
        RETURN  
      END FUNCTION MULTIMOD

C     ==========================================================
C     LCG_INIT
C       THIS PROCEDURE MUST BE CALLED BEFORE ANY OF THE OTHER
C       RANDOM ROUTINES PROVIDED. 
C       YOU CAN OPTIONALLY CHOOSE THE INITIAL SEED AND THE SIZE OF
C       THE SEQUENCE ASSIGNED TO EACH PROCESSOR (E.G. THE INITIAL
C       VALUES OF THE  SUB-SEQUENCE IN TWO DIFFERENT PROCESSORS 
C       WILL BE 2**GAP VALUES APART IN THE ORIGINAL SEQUENCE). 
C     ==========================================================

      SUBROUTINE LCG_INIT(S,SEED,GAP)

        USE LCG_PRIVATE_VARIABLES
        USE DEBUG_VARS

        INTEGER, DIMENSION(:), INTENT(OUT) :: S
        INTEGER, OPTIONAL, INTENT(IN)      :: SEED
        INTEGER, OPTIONAL, INTENT(IN)      :: GAP

        INTEGER  :: I,NPROC

        LAG = 50
        IF (PRESENT(GAP)) THEN LAG = GAP

C       COMPUTES A**(2**LAG) MOD M  
C       ITERATIVILY SQUARE A MODULUS M, LAG TIMES

        A_2_LAG = A
        DO I = 1,LAG
           A_2_LAG = MULTIMOD(A_2_LAG,A_2_LAG,M)
        END DO

        IF (NOT(DEBUG .LT. LEVEL2)) THEN
           print *,'LAG = ', LAG
           PRINT *,'A**(2**LAG) = ', A_2_LAG
        END IF 

C       COMPUTES THE SEED OF THE SEQUENCE OF EACH PROCESSOR

        IF (PRESENT(SEED)) THEN
          S(1) = SEED
        ELSE
          S(1) = 1234567890
        END IF

        DO NPROC = 2,NUMBER_OF_PROCESSORS()
          S(NPROC) = MULTIMOD(S(NPROC-1),A_2_LAG,M)
        END DO

        IF (NOT(DEBUG .LT. LEVEL2)) THEN
            print *,'INITIAL VALUES ',S
        END IF

      END SUBROUTINE LCG_INIT


C     ==========================================================
C     INITIALIZES THE SECTION OF THE VECTOR MAPPED TO THE LOCAL
C     PROCESSOR USING THE LCG. ADVANCES THE STATE OF THE GENERATOR
C     PRODUCES INTEGER PSEUDO-RANDOM NUMBERS
C     ==========================================================

      EXTRINSIC(HPF_LOCAL) SUBROUTINE LCG_IRANDOM(S,V)
        USE LCG_PRIVATE_VARIABLES

        INTEGER, DIMENSION(:), INTENT(OUT)  :: S
        INTEGER, DIMENSION(:), INTENT(OUT)  :: V

        INTEGER  :: I

        DO I=LBOUND(V,1),UBOUND(V,1)     
           S(LBOUND(S,1)) = MULTIMOD(A,S(LBOUND(S,1)),M)
           V(I)           = S(LBOUND(S,1))
        END DO
      END SUBROUTINE LCG_IRANDOM


C     ==========================================================
C     INITIALIZES THE SECTION OF THE VECTOR MAPPED TO THE LOCAL
C     PROCESSOR USING THE LCG. ADVANCES THE STATE OF THE GENERATOR
C     PRODUCES REAL PSEUDO-RANDOM NUMBERS IN THE RANGE (0,1)
C     ==========================================================

      EXTRINSIC(HPF_LOCAL) SUBROUTINE LCG_RRANDOM(S,V)
        USE LCG_PRIVATE_VARIABLES

        INTEGER, DIMENSION(:), INTENT(OUT)  :: S
        REAL, DIMENSION(:), INTENT(OUT)     :: V

        INTEGER  :: I

        DO I=LBOUND(V,1),UBOUND(V,1)     
           S(LBOUND(S,1)) = MULTIMOD(A,S(LBOUND(S,1)),M)
           V(I)           = REAL(S(LBOUND(S,1))) / REAL(M)
        END DO
      END SUBROUTINE LCG_RRANDOM



