! This code calculates the full-range interactions amongst a
! number of particles, with no cut-off. The code also avoids
! superfluous computations by exploiting Newton's third law to
! halve the possible computations. The timing information output lists time spent
! executing statements with communication required, without communication, and
! others.  The no-communciation times should scale nicely, but those with
! communication obviously will not.


      PROGRAM nbody
      IMPLICIT NONE
C _____________________________________________________________________
C|                                                                     |
C| Direct N-body simulation code.                                      |
C| Computes the forces between N-bodies after one time-step.           |
C|                                                                     |
C| MM: mass.                                                           |
C| PX, PY, PZ: position in x, y, z coordinates.                        |
C| FX, FY, FZ: force in x, y, z coordinates.                           |
C|                                                                     |
C| H W Yau.  25th of August, 1996.                                     |
C| Northeast Parallel Architectures Center.                            |
C| Syracuse University.                                                |
C|_____________________________________________________________________|
C|                                                                     |
C| Edit record:                                                        |
C| 06/Jun/96 Cleaned up by M. McMahon.  Modified all the array         |
C|           statements to dimension.                                  |
C|           Modified by M.McMahon 6/6/96 Changed EXTRINSIC and ALIGN  |
C|           statements.                                               |
C| 17/Jun/96 (MM) fixed up timer.                                      |
C| 22/Aug/1996: (HWY) Made problem smaller from NN=4096.               |
C|              Fixed code to remove hard-wiring of `4096'.            |
C| 25/Aug/1996: Tidying up.                                            |
C|              Removed superfluous `NP' parameter.                    |
C|              Fortran 90 version of the same code.                   |
C| 08/Apr/1997: (HWY) Bug fixes.                                       |
C|              Rewrote initialisation of states, to be predictable.   |
C| 09/Apr/1997: Use of allocatable arrays.                             |
C|_____________________________________________________________________|
      INTERFACE
         SUBROUTINE timer(return_time, initial_time)
         REAL, INTENT(IN)  :: initial_time
         REAL, INTENT(OUT) :: return_time
         END SUBROUTINE
      END INTERFACE
      INTEGER, PARAMETER :: NN = 3200
      INTEGER, PARAMETER :: MM=0, PX=1, PY=2, PZ=3, FX=4, FY=5, FZ=6
      REAL, PARAMETER :: G = 1.0
      REAL,DIMENSION(:),ALLOCATABLE :: dx,dy,dz,tx,ty,tz,
     1                                 sq,dist,fac
      REAL,DIMENSION(:,:),ALLOCATABLE :: p,q
      INTEGER :: i,j,k,npts
      REAL    :: ti, tf,tcalc,tcom,start,a
C HPF directives
! First a 1D abstract processor array is declared, with its size (NUMBER_OF_PROCESSORS()
! given at run time.  A template is created and distributed in a block fashion onto
! the processor array.  The 1D arrays in the code are aligned with the template,
! thereby insuring that they are distributed in the same BLOCK fashion.  
! The 2d arrays have their 1st dimension aligned with the template, lso it is distributed
! in the same manner as the template.  The second dimension of these arrays is
! collapsed, as the * indicates in the distribute directive .

!HPF$ PROCESSORS :: procs(NUMBER_OF_PROCESSORS())
!HPF$ TEMPLATE,DIMENSION(0:NN-1) :: tplate
!HPF$ DISTRIBUTE (BLOCK) ONTO procs :: tplate
!HPF$ ALIGN dx(:)  WITH tplate(:)
!HPF$ ALIGN dy(:)  WITH tplate(:)
!HPF$ ALIGN dz(:)  WITH tplate(:)
!HPF$ ALIGN tx(:)  WITH tplate(:)
!HPF$ ALIGN ty(:)  WITH tplate(:)
!HPF$ ALIGN tz(:)  WITH tplate(:)
!HPF$ ALIGN sq(:)  WITH tplate(:)
!HPF$ ALIGN dist(:)WITH tplate(:)
!HPF$ ALIGN fac(:) WITH tplate(:)
!HPF$ ALIGN p(:,*) WITH tplate(:)
!HPF$ ALIGN q(:,*) WITH tplate(:)
C Executable Code.
      npts = NN
      WRITE(UNIT=*,FMT=101) NN
      READ  *, npts
      IF( npts .GT. NN )THEN
         WRITE(UNIT=*,FMT=102) NN
         STOP
      END IF
C Allocate arrays.
      ALLOCATE( dx(0:npts-1) )
      ALLOCATE( dy(0:npts-1) )
      ALLOCATE( dz(0:npts-1) )
      ALLOCATE( tx(0:npts-1) )
      ALLOCATE( ty(0:npts-1) )
      ALLOCATE( tz(0:npts-1) )
      ALLOCATE( sq(0:npts-1) )
      ALLOCATE( dist(0:npts-1) )
      ALLOCATE( fac(0:npts-1) )
      ALLOCATE( p(0:npts-1,0:6) )
      ALLOCATE( q(0:npts-1,0:6) )
      call timer(start,0.0)
      call timer(ti,0.0)
      tcalc = 0.0
      tcom  = 0.0
C Define inital state of particles.
      p = 0
      FORALL( i=0:npts-1 )
         p(i,MM) = 10.05 + i
         p(i,PX) = 30.0 * i
         p(i,PY) = 20.0 * i
         p(i,PZ) = 10.0 * i
      END FORALL
C  Initialize q to be the same as p.  As the force between particle i and j
C  is calculated and added into p, q will be used to accumulate the
C  symmetrical opposite force from j to i.
      q = p
      k = ((npts - 1) / 2) - 1
      call timer(tcalc,ti)
      do i = 0, k
         call timer(ti,0.0)
         q = CSHIFT(q,DIM=1,SHIFT=-1)
         call timer(a,ti)
         tcom = tcom + a
         call timer(ti,0.0)
         dx(:) = p(:,PX) - q(:,PX)
         dy(:) = p(:,PY) - q(:,PY)
         dz(:) = p(:,PZ) - q(:,PZ)
         sq(:) = dx(:)**2 + dy(:)**2 + dz(:)**2
         dist(:) = SQRT(sq(:))
         fac(:) = p(:,MM) * q(:,MM) / (dist(:) * sq(:))
         tx(:) = fac(:) * dx(:)
         ty(:) = fac(:) * dy(:)
         tz(:) = fac(:) * dz(:)
         p(:,FX) = p(:,FX) - tx(:)
         p(:,FY) = p(:,FY) - ty(:)
         p(:,FZ) = p(:,FZ) - tz(:)
         q(:,FX) = q(:,FX) + tx(:)
         q(:,FY) = q(:,FY) + ty(:)
         q(:,FZ) = q(:,FZ) + tz(:)
         call timer(a,ti)
         tcalc = tcalc + a
      end do
C Final one-way forces.
      IF( MOD(npts,2) .EQ. 0 )THEN
         CALL timer(ti,0.0)
         q = CSHIFT(q,DIM=1,SHIFT=-1)
         call timer(a,ti)
         tcom = tcom + a
         call timer(ti,0.0)
         dx(:) = p(:,PX) - q(:,PX)
         dy(:) = p(:,PY) - q(:,PY)
         dz(:) = p(:,PZ) - q(:,PZ)
         sq(:) = dx(:)**2 + dy(:)**2 + dz(:)**2
         dist(:) = SQRT(sq(:))
         fac(:) = p(:,MM) * q(:,MM) / (dist(:) * sq(:))
         tx(:) = fac(:) * dx(:)
         ty(:) = fac(:) * dy(:)
         tz(:) = fac(:) * dz(:)
         p(:,FX) = p(:,FX) - tx(:)
         p(:,FY) = p(:,FY) - ty(:)
         p(:,FZ) = p(:,FZ) - tz(:)
         call timer(a,ti)
         tcalc = tcalc + a
      END IF
C Combine final forces for the final result
C CSHIFT operation to bring q() back to original indexing state.
      call timer(ti,0.0)
      q = CSHIFT( q,DIM=1,SHIFT=(npts/2) )
      call timer(a,ti)
      tcom = tcom + a
      call timer(ti,0.0)
      p(:,FX) = p(:,FX) + q(:,FX)
      p(:,FY) = p(:,FY) + q(:,FY)
      p(:,FZ) = p(:,FZ) + q(:,FZ)
      call timer(a,ti)
      tcalc = tcalc + a
C Stop timer
      call timer(tf,start)
      j = MIN(65,npts-1)
      print *, 'State After='
      print *, 'p(FX)=',(p(i,FX),i=0,j)
      print *, 'p(FY)=',(p(i,FY),i=0,j)
      print *, 'p(FZ)=',(p(i,FZ),i=0,j)
      WRITE(6,6001) NUMBER_OF_PROCESSORS(),npts,
     1    tcom,tcalc,(tf-tcom-tcalc),tf
 101  FORMAT('Enter the number of data points - <= ',I8)
 102  FORMAT('The choice of npts=',I8,' will probably cause a crash.')
 6001 FORMAT('Number of Processors = ',I4/
     1     'Problem size = ',I5/
     2     'Communications = ',F9.3/
     3     'Compute        = ',F9.3/
     4     'Others         = ',F9.3/
     5     'Total time     = ',F9.3)
      STOP
      END
C======================================================================
      SUBROUTINE timer(return_time, initial_time)
      real, intent(in)  :: initial_time
      real,  intent(out) :: return_time
      integer finish,rate
      CALL system_clock( COUNT=finish,COUNT_RATE=rate)
      return_time = FLOAT(finish) / FLOAT(rate) - initial_time
      RETURN
      END
C======================================================================
