!--------------------------
! The original code was written by Tom Haupt and a student and later, in Fall 1996,
! as part of class assignment, Ulf Dittmer modified the random generator part of the code.
! 1996 NPAC.
!--------------------------
PROGRAM nbody

INTERFACE
	EXTRINSIC (hpf_local) SUBROUTINE timer (return_time, initial_time)
		REAL, DIMENSION(:), INTENT(in)  :: initial_time
		REAL, DIMENSION(:), INTENT(out) :: return_time
!HPF$ DISTRIBUTE *(BLOCK) :: initial_time
!HPF$ ALIGN return_time WITH initial_time
	END SUBROUTINE timer

	SUBROUTINE rng_seed(seed)
		INTEGER, INTENT(out) :: seed
	END SUBROUTINE rng_seed

	SUBROUTINE rng_number (x, N, seed)
		INTEGER :: N, seed
		REAL, DIMENSION(0:N-1), INTENT(out) :: x
	END SUBROUTINE rng_number

	SUBROUTINE p_rng_seed(seed)
		INTEGER, DIMENSION(:), INTENT(out) :: seed
	END SUBROUTINE p_rng_seed

	EXTRINSIC (hpf_local) SUBROUTINE p_rng_number (x, N, seed)
		INTEGER :: N
		INTEGER, DIMENSION(:), INTENT(inout) :: seed
		REAL, DIMENSION(:), INTENT(out) :: x
!HPF$ DISTRIBUTE *(BLOCK) :: x, seed
	END SUBROUTINE p_rng_number
END INTERFACE

	INTEGER, PARAMETER :: NPROCS = 1
! change NPROCS to the number of processors used
	INTEGER, PARAMETER :: NN = 65536
	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, DIMENSION(0:NN-1) :: x
!	REAL, DIMENSION(0:NN-1) :: t,dx,dy,dz,tx,ty,tz,sq,dist,fac

	REAL, DIMENSION(0:NN-1,0:6) :: p
!	REAL, DIMENSION(0:NN-1,0:6) :: q
	INTEGER :: i, j, k, npts
	REAL, DIMENSION(NPROCS) :: ti, tf
	INTEGER, DIMENSION(NPROCS) :: rnd
! Strange: When compiling from the kestrel command line, in the two
! preceding lines you could write "..., DIMENSION(NUMBER_OF_PROCS())..."
! In VPL, the compiler complains about that.

!HPF$ DISTRIBUTE (BLOCK) :: ti, tf, rnd, x
!HPF$ ALIGN (:,*) WITH x(:) :: p

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

!	print *, 'Enter the number of data points'
!	read  *, npts
	npts = 1024


!	Start Timer
	call timer(tf, ti)

!	call RANDOM_SEED
	call p_rng_seed(rnd)

	print *, 'seeds: ', rnd

!	call RANDOM_NUMBER(x)
	call p_rng_number(x, npts, rnd)
	p(:,MM) = x(:)

!	call RANDOM_NUMBER(x)
	call p_rng_number(x, npts, rnd)
	p(:,PX) = x(:)

!	call RANDOM_NUMBER(x)
	call p_rng_number(x, npts, rnd)
	p(:,PY) = x(:)

!	call RANDOM_NUMBER(x)
	call p_rng_number(x, npts, rnd)
	p(:,PZ) = x(:)

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

! Stop timer
	call timer(ti, tf)


!	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)

	print *, 'running time ', ti, 'seconds'

!	print *, 'FX', p(0:6,FX)
!	print *, 'FY', p(0:6,FY)
!	print *, 'FZ', p(0:6,FZ)

	END

!======================================================================

EXTRINSIC (hpf_local) SUBROUTINE timer (return_time, initial_time)
	REAL, DIMENSION(:), INTENT(in)  :: initial_time
	REAL, DIMENSION(:), INTENT(out) :: return_time

!HPF$ DISTRIBUTE *(BLOCK) :: initial_time
!HPF$ ALIGN return_time WITH initial_time
	
	return_time(lbound(return_time, 1)) = secnds(initial_time(lbound(initial_time, 1)))

END SUBROUTINE timer

!======================================================================

SUBROUTINE rng_seed(seed)
	IMPLICIT none

	INTEGER, PARAMETER :: A = 1103515245
	INTEGER, PARAMETER :: C = 12345
	INTEGER, PARAMETER :: MASK = Z"7FFFFFFF"

	INTEGER, INTENT(out) :: seed
	CHARACTER(LEN=10) :: date, time, zone
	INTEGER, DIMENSION(8) :: values

	call DATE_AND_TIME(date, time, zone, values)
	seed = IAND(A * (values(8) * C + values(7)) + C, MASK)
END SUBROUTINE rng_seed

!======================================================================

SUBROUTINE rng_number (x, N, seed)
	IMPLICIT none

	INTEGER, PARAMETER :: A = 1103515245
	INTEGER, PARAMETER :: C = 12345
	INTEGER, PARAMETER :: MASK = Z"7FFFFFFF"	! lower 31 bits

	INTEGER :: N, i, seed
	REAL, DIMENSION(0:N-1), INTENT(out) :: x

	DO i=0, N-1
		seed = IAND(A * seed + C, MASK)
		x(i) = seed / 2147483648.0
	END DO
END SUBROUTINE rng_number

!======================================================================

SUBROUTINE p_rng_seed(seed)
	IMPLICIT none

	INTEGER, PARAMETER :: A = 1103515245
	INTEGER, PARAMETER :: C = 12345
	INTEGER, PARAMETER :: MASK = Z"7FFFFFFF"	! lower 31 bits

	INTEGER, DIMENSION(:), INTENT(out) :: seed
	CHARACTER(LEN=10) :: date, time, zone
	INTEGER, DIMENSION(8) :: values
	INTEGER :: i

	call DATE_AND_TIME(date, time, zone, values)
	seed(1) = IAND(A * (values(8) * C + values(7)) + C, MASK)
	DO i=2, ubound(seed, 1)
		seed(i) = IAND(ISHFT(A * seed(i-1) + C, -3), MASK)
	END DO
END SUBROUTINE p_rng_seed

!======================================================================

EXTRINSIC (hpf_local) SUBROUTINE p_rng_number (x, N, seed)
	IMPLICIT none

	INTEGER, PARAMETER :: A = 1103515245
	INTEGER, PARAMETER :: C = 12345
	INTEGER, PARAMETER :: MASK = Z"7FFFFFFF"	! lower 31 bits

	INTEGER :: N, i, j
	INTEGER, DIMENSION(:), INTENT(inout) :: seed
	REAL, DIMENSION(:), INTENT(out) :: x

!HPF$ DISTRIBUTE *(BLOCK) :: x, seed

	j = lbound(seed, 1)

	DO i = lbound(x, 1), ubound(x, 1)
		seed(j) = IAND(A * seed(j) + C, MASK)
		x(i) = seed(j) / 2147483648.0
	END DO
END SUBROUTINE p_rng_number