MPP-traffic.tar

traffic/C-SER/0000775000175000017500000000000012620656170011432 5ustar  dshdshtraffic/C-SER/uni.h0000664000175000017500000000004712534300645012374 0ustar  dshdshfloat uni(void);
void rinit(int seed);
traffic/C-SER/traffic.c0000664000175000017500000000302612551223420013204 0ustar  dshdsh#include <stdlib.h>
#include <stdio.h>

#include "traffic.h"

#define NCELL 100000

int main(int argc, char **argv)
{
  int *oldroad, *newroad;

  int i, iter, nmove, ncars;
  int maxiter, printfreq; 

  float density; 

  double tstart, tstop;

  oldroad = (int *) malloc((NCELL+2)*sizeof(int));
  newroad = (int *) malloc((NCELL+2)*sizeof(int));

  maxiter = 200000000/NCELL; 
  printfreq = maxiter/10; 

  // Set target density of cars

  density = 0.52;

  printf("Length of road is %d\n", NCELL);
  printf("Number of iterations is %d \n", maxiter);
  printf("Target density of cars is %f \n", density);

  // Initialise road accordingly using random number generator

  printf("Initialising road ...\n");
  
  ncars = initroad(&oldroad[1], NCELL, density, SEED);

  printf("...done\n");

  printf("Actual density of cars is %f\n\n", (float) ncars / (float) NCELL);

  tstart = gettime();

  for (iter=1; iter<=maxiter; iter++)
    { 
      updatebcs(oldroad, NCELL);

      // Apply CA rules to all cells

      nmove = updateroad(newroad, oldroad, NCELL);
      
      // Copy new to old array

      for (i=1; i<=NCELL; i++)
	{
	  oldroad[i] = newroad[i]; 
	}

      if (iter%printfreq == 0)
	{
	  printf("At iteration %d average velocity is %f \n",
		 iter, (float) nmove / (float) ncars);
	} 
    }

  tstop = gettime();

  free(oldroad);
  free(newroad);
  
  printf("\nFinished\n");
  printf("\nTime taken was  %f seconds\n", tstop-tstart);
  printf("Update rate was %f MCOPs\n\n", \
	 1.e-6*((double) NCELL)*((double) maxiter)/(tstop-tstart));
}
traffic/C-SER/trafficlib.c0000664000175000017500000000177712540550517013716 0ustar  dshdsh#include "uni.h"

int initroad(int *road, int n, float density, int seed)
{
  int i, ncar;
  float rng;

  // seed random number generator

  rinit(seed);

  ncar = 0;

  for (i=0; i < n; i++)
    {
      rng = uni();

      if (rng < density)
	{
	  road[i] = 1;
	}
      else
	{
	  road[i] = 0;
	}

      ncar += road[i];
    }

  return ncar;
}


int updateroad(int *newroad, int *oldroad, int n)
{
  int i, nmove;

  nmove = 0;

  for (i=1; i<=n; i++)
    {
      if (oldroad[i] == 1)
	{
	  if (oldroad[i+1] == 1)
	    {
	      newroad[i] = 1;  
	    }
	  else
	    {
	      newroad[i] = 0; 
	      nmove++; 
	    }
	}
      else
	{
	  if (oldroad[i-1] == 1)
	    { 
	      newroad[i] = 1;  
	    }
	  else
	    {
	      newroad[i] = 0; 
	    }
	}
    }

  return nmove;
}


void updatebcs(int *road, int n)
{
  road[0]   = road[n];
  road[n+1] = road[1];
}


#include <stdlib.h>
#include <sys/time.h>

double gettime()
{
  struct timeval tp;
  gettimeofday (&tp, NULL);
  return tp.tv_sec + tp.tv_usec/(double)1.0e6;
}
traffic/C-SER/uni.c0000664000175000017500000000650612534300645012375 0ustar  dshdsh/*
 *	C version of Marsaglia's UNI random number generator
 *	More or less transliterated from the Fortran -- with 1 bug fix
 *	Hence horrible style
 *
 *	Features:
 *		ANSI C
 *		not callable from Fortran (yet)
 */


/*
 *	Global variables for rstart & uni
 */
#include <stdio.h>
#include <stdlib.h>

float uni_u[98];	/* Was U(97) in Fortran version -- too lazy to fix */
float uni_c, uni_cd, uni_cm;
int uni_ui, uni_uj;

float uni(void)
{
	float luni;			/* local variable for uni */

	luni = uni_u[uni_ui] - uni_u[uni_uj];
	if (luni < 0.0)
		luni += 1.0;
	uni_u[uni_ui] = luni;
	if (--uni_ui == 0)
		uni_ui = 97;
	if (--uni_uj == 0)
		uni_uj = 97;
	if ((uni_c -= uni_cd) < 0.0)
		uni_c += uni_cm;
	if ((luni -= uni_c) < 0.0)
		luni += 1.0;
	return (float) luni;
}

void rstart(int i, int j, int k, int l)
{
	int ii, jj, m;
	float s, t;

	for (ii = 1; ii <= 97; ii++) {
		s = 0.0;
		t = 0.5;
		for (jj = 1; jj <= 24; jj++) {
			m = ((i*j % 179) * k) % 179;
			i = j;
			j = k;
			k = m;
			l = (53*l+1) % 169;
			if (l*m % 64 >= 32)
				s += t;
			t *= 0.5;
		}
		uni_u[ii] = s;
	}
	uni_c  = 362436.0   / 16777216.0;
	uni_cd = 7654321.0  / 16777216.0;
	uni_cm = 16777213.0 / 16777216.0;
	uni_ui = 97;	/*  There is a bug in the original Fortran version */
	uni_uj = 33;	/*  of UNI -- i and j should be SAVEd in UNI()     */
}


/* ~rinit: this takes a single integer in the range
		0 <= ijkl <= 900 000 000
	and produces the four smaller integers needed for rstart. It is
 *	based on the ideas contained in the RMARIN subroutine in
 *		F. James, "A Review of Pseudorandom Number Generators",
 *			Comp. Phys. Commun. Oct 1990, p.340
 *	To reduce the modifications to the existing code, rinit now
 *	takes the role of a preprocessor for rstart.
 *
 *	This is useful for the parallel version of the code as James
 *	states that any integer ijkl will produce a statistically
 *	independent sequence of random numbers.
 *
 *     Very funny. If that statement was worth anything he would have provided
 *     a proof to go with it. spb 12/12/90 
 */

void rinit(int ijkl)
{
	int i, j, k, l, ij, kl;

	/* check ijkl is within range */
	if( (ijkl < 0) || (ijkl > 900000000) )
		{
		printf("rinit: ijkl = %d -- out of range\n\n", ijkl);
		exit(3);
               	}

/*        printf("rinit: seed_ijkl = %d\n", ijkl); */

	/* decompose the long integer into the the equivalent four
	 * integers for rstart. This should be a 1-1 mapping
 	 *	ijkl <--> (i, j, k, l)
	 * though not quite all of the possible sets of (i, j, k, l)
	 * can be produced.
	 */

	ij = ijkl/30082;
	kl = ijkl - (30082 * ij);

	i = ((ij/177) % 177) + 2;
	j = (ij % 177) + 2;
	k = ((kl/169) % 178) + 1;
	l = kl % 169;

	if( (i <= 0) || (i > 178) )
		{
		printf("rinit: i = %d -- out of range\n\n", i);
		exit(3);
               	}

	if( (j <= 0) || (j > 178) )
		{
		printf("rinit: j = %d -- out of range\n\n", j);
		exit(3);
               	}

	if( (k <= 0) || (k > 178) )
		{
		printf("rinit: k = %d -- out of range\n\n", k);
		exit(3);
               	}

	if( (l < 0) || (l > 168) )
		{
		printf("rinit: l = %d -- out of range\n\n", l);
		exit(3);
               	}

	if (i == 1 && j == 1 && k == 1)
		{
                printf("rinit: 1 1 1 not allowed for 1st 3 seeds\n\n");
		exit(4);
                }

/*        printf("rinit: initialising RNG via rstart(%d, %d, %d, %d)\n",
				i, j, k, l); */

        rstart(i, j, k, l);

}

traffic/C-SER/Makefile0000664000175000017500000000063712620656170013100 0ustar  dshdshMF=	Makefile

# For morar
CC=	pgcc

# For ARCHER
#CC=	cc

CFLAGS=
LFLAGS=	-lm

EXE=	traffic

INC= \
	traffic.h \
	uni.h

SRC= \
	traffic.c \
	trafficlib.c \
	uni.c

#
# No need to edit below this line
#

.SUFFIXES:
.SUFFIXES: .c .o

OBJ=	$(SRC:.c=.o)

.c.o:
	$(CC) $(CFLAGS) -c $<

all:	$(EXE)

$(OBJ):	$(INC)

$(EXE):	$(OBJ)
	$(CC) $(CFLAGS) -o $@ $(OBJ) $(LFLAGS)

$(OBJ):	$(MF)

clean:
	rm -f $(OBJ) $(EXE) core
traffic/C-SER/traffic.h0000664000175000017500000000026412540550566013226 0ustar  dshdsh#define SEED  5743

int initroad(int *road, int n, float density, int seed);
void updatebcs(int *road, int n);
int updateroad(int *newroad, int *oldroad, int n);
double gettime();
traffic/C-MPI/0000775000175000017500000000000012633773606011436 5ustar  dshdshtraffic/C-MPI/uni.h0000664000175000017500000000004712534300645012370 0ustar  dshdshfloat uni(void);
void rinit(int seed);
traffic/C-MPI/traffic.c0000664000175000017500000000423712540772577013231 0ustar  dshdsh#include <stdlib.h>
#include <stdio.h>

#include "traffic.h"

#define NCELL 100000

int main(int argc, char **argv)
{
  int *oldroad, *newroad, *bigroad;

  int i, iter, nmove, ncars;
  int maxiter, printfreq; 

  float density; 

  double tstart, tstop;

  int procid, nproc, nlocal;

  maxiter = 200000000/NCELL; 
  printfreq = maxiter/10; 

  // Set target density of cars

  density = 0.52;

  // Start message passing system

  mpstart(&nproc, &procid);

  if (procid == 0)
    {
      printf("Length of road is %d\n", NCELL);
      printf("Number of iterations is %d \n", maxiter);
      printf("Target density of cars is %f \n", density);
      printf("Running on %d processes\n", nproc);
    }

  nlocal = NCELL/nproc;

  oldroad = (int *) malloc((nlocal+2)*sizeof(int));
  newroad = (int *) malloc((nlocal+2)*sizeof(int));

  if (procid == 0)
    {
      bigroad = (int *) malloc(NCELL*sizeof(int));

      // Initialise road accordingly using random number generator

      printf("Initialising road ...\n");
  
      ncars = initroad(bigroad, NCELL, density, SEED);

      printf("...done\n");
      printf("Actual density is %f\n", (float) ncars / (float) NCELL);
      printf("Scattering data ...\n");
    }

  mpscatter(bigroad, &oldroad[1], nlocal);

  if (procid == 0)
    {
      printf("... done\n\n");
    }

  tstart = gettime();

  for (iter=1; iter<=maxiter; iter++)
    { 
      mpupdatebcs(oldroad, nlocal, procid, nproc);

      // Apply CA rules to all cells

      nmove = updateroad(newroad, oldroad, nlocal);

      // Globally sum the value

      mpgsum(&nmove);

      // Copy new to old array

      for (i=1; i<=nlocal; i++)
	{
	  oldroad[i] = newroad[i]; 
	}

      if (iter%printfreq == 0)
	{
	  if (procid == 0)
	    {
	      printf("At iteration %d average velocity is %f \n",
		     iter, (float) nmove / (float) ncars);
	    }
	} 
    }

  tstop = gettime();

  free(oldroad);
  free(newroad);

  if (procid == 0)
    {
      free(bigroad);

      printf("\nFinished\n");
      printf("\nTime taken was  %f seconds\n", tstop-tstart);
      printf("Update rate was %f MCOPs\n\n", \
	     1.e-6*((double) NCELL)*((double) maxiter)/(tstop-tstart));
    }

  mpstop();
}
traffic/C-MPI/trafficlib.c0000664000175000017500000000221012534300645013667 0ustar  dshdsh#include "uni.h"
#include "mplib.h"

int initroad(int *road, int n, float density, int seed)
{
  int i, ncar;
  float rng;

  // seed random number generator

  rinit(seed);

  ncar = 0;

  for (i=0; i < n; i++)
    {
      rng = uni();

      if (rng < density)
	{
	  road[i] = 1;
	}
      else
	{
	  road[i] = 0;
	}

      ncar += road[i];
    }

  return ncar;
}


int updateroad(int *newroad, int *oldroad, int n)
{
  int i, nmove;

  nmove = 0;

  for (i=1; i<=n; i++)
    {
      if (oldroad[i] == 1)
	{
	  if (oldroad[i+1] == 1)
	    {
	      newroad[i] = 1;  
	    }
	  else
	    {
	      newroad[i] = 0; 
	      nmove++; 
	    }
	}
      else
	{
	  if (oldroad[i-1] == 1)
	    { 
	      newroad[i] = 1;  
	    }
	  else
	    {
	      newroad[i] = 0; 
	    }
	}
    }

  return nmove;
}


void updatebcs(int *road, int n)
{
  road[0]   = road[n];
  road[n+1] = road[1];
}


void mpupdatebcs(int *road, int n, int procid, int nproc)
{
  int procup, procdown;

  procup   = (procid + 1) % nproc;
  procdown = (procid + nproc - 1) % nproc;

  mpsendrecv(&road[n], procup,   &road[0],   procdown);
  mpsendrecv(&road[1], procdown, &road[n+1], procup  );
}
traffic/C-MPI/mplib.c0000664000175000017500000000145212534300645012674 0ustar  dshdsh#include <stdio.h>
#include <mpi.h> 

void mpstart(int *numproc, int *procid)
{ 
  MPI_Init(NULL, NULL);
  MPI_Comm_rank(MPI_COMM_WORLD, procid);
  MPI_Comm_size(MPI_COMM_WORLD, numproc);
}
  
void mpstop(void)
{ 
  MPI_Finalize();
}
 
void mpgsum(int *ival)
{
  int itmp;
 
  MPI_Allreduce(ival, &itmp, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);

  *ival = itmp;
} 

void mpscatter(int *inbuff, int *outbuff, int ndata)
{
  MPI_Scatter(inbuff,  ndata, MPI_INT,
	      outbuff, ndata, MPI_INT,
	      0, MPI_COMM_WORLD);
}

void mpsendrecv(int *senddata, int sendproc, int *recvdata, int recvproc)
{
  MPI_Status status;
    
  MPI_Sendrecv(senddata, 1, MPI_INT, sendproc, 1,
	       recvdata, 1, MPI_INT, recvproc, 1,
               MPI_COMM_WORLD, &status);
}
  
double gettime(void)
{ 
  return MPI_Wtime(); 
} 
traffic/C-MPI/uni.c0000664000175000017500000000650612534300645012371 0ustar  dshdsh/*
 *	C version of Marsaglia's UNI random number generator
 *	More or less transliterated from the Fortran -- with 1 bug fix
 *	Hence horrible style
 *
 *	Features:
 *		ANSI C
 *		not callable from Fortran (yet)
 */


/*
 *	Global variables for rstart & uni
 */
#include <stdio.h>
#include <stdlib.h>

float uni_u[98];	/* Was U(97) in Fortran version -- too lazy to fix */
float uni_c, uni_cd, uni_cm;
int uni_ui, uni_uj;

float uni(void)
{
	float luni;			/* local variable for uni */

	luni = uni_u[uni_ui] - uni_u[uni_uj];
	if (luni < 0.0)
		luni += 1.0;
	uni_u[uni_ui] = luni;
	if (--uni_ui == 0)
		uni_ui = 97;
	if (--uni_uj == 0)
		uni_uj = 97;
	if ((uni_c -= uni_cd) < 0.0)
		uni_c += uni_cm;
	if ((luni -= uni_c) < 0.0)
		luni += 1.0;
	return (float) luni;
}

void rstart(int i, int j, int k, int l)
{
	int ii, jj, m;
	float s, t;

	for (ii = 1; ii <= 97; ii++) {
		s = 0.0;
		t = 0.5;
		for (jj = 1; jj <= 24; jj++) {
			m = ((i*j % 179) * k) % 179;
			i = j;
			j = k;
			k = m;
			l = (53*l+1) % 169;
			if (l*m % 64 >= 32)
				s += t;
			t *= 0.5;
		}
		uni_u[ii] = s;
	}
	uni_c  = 362436.0   / 16777216.0;
	uni_cd = 7654321.0  / 16777216.0;
	uni_cm = 16777213.0 / 16777216.0;
	uni_ui = 97;	/*  There is a bug in the original Fortran version */
	uni_uj = 33;	/*  of UNI -- i and j should be SAVEd in UNI()     */
}


/* ~rinit: this takes a single integer in the range
		0 <= ijkl <= 900 000 000
	and produces the four smaller integers needed for rstart. It is
 *	based on the ideas contained in the RMARIN subroutine in
 *		F. James, "A Review of Pseudorandom Number Generators",
 *			Comp. Phys. Commun. Oct 1990, p.340
 *	To reduce the modifications to the existing code, rinit now
 *	takes the role of a preprocessor for rstart.
 *
 *	This is useful for the parallel version of the code as James
 *	states that any integer ijkl will produce a statistically
 *	independent sequence of random numbers.
 *
 *     Very funny. If that statement was worth anything he would have provided
 *     a proof to go with it. spb 12/12/90 
 */

void rinit(int ijkl)
{
	int i, j, k, l, ij, kl;

	/* check ijkl is within range */
	if( (ijkl < 0) || (ijkl > 900000000) )
		{
		printf("rinit: ijkl = %d -- out of range\n\n", ijkl);
		exit(3);
               	}

/*        printf("rinit: seed_ijkl = %d\n", ijkl); */

	/* decompose the long integer into the the equivalent four
	 * integers for rstart. This should be a 1-1 mapping
 	 *	ijkl <--> (i, j, k, l)
	 * though not quite all of the possible sets of (i, j, k, l)
	 * can be produced.
	 */

	ij = ijkl/30082;
	kl = ijkl - (30082 * ij);

	i = ((ij/177) % 177) + 2;
	j = (ij % 177) + 2;
	k = ((kl/169) % 178) + 1;
	l = kl % 169;

	if( (i <= 0) || (i > 178) )
		{
		printf("rinit: i = %d -- out of range\n\n", i);
		exit(3);
               	}

	if( (j <= 0) || (j > 178) )
		{
		printf("rinit: j = %d -- out of range\n\n", j);
		exit(3);
               	}

	if( (k <= 0) || (k > 178) )
		{
		printf("rinit: k = %d -- out of range\n\n", k);
		exit(3);
               	}

	if( (l < 0) || (l > 168) )
		{
		printf("rinit: l = %d -- out of range\n\n", l);
		exit(3);
               	}

	if (i == 1 && j == 1 && k == 1)
		{
                printf("rinit: 1 1 1 not allowed for 1st 3 seeds\n\n");
		exit(4);
                }

/*        printf("rinit: initialising RNG via rstart(%d, %d, %d, %d)\n",
				i, j, k, l); */

        rstart(i, j, k, l);

}

traffic/C-MPI/Makefile0000664000175000017500000000066612620655617013103 0ustar  dshdshMF=	Makefile

# For morar
CC=	mpicc

# For ARCHER
#CC=	cc

CFLAGS=
LFLAGS=	-lm

EXE=	traffic

INC= \
	traffic.h \
	uni.h \
	mplib.h

SRC= \
	traffic.c \
	trafficlib.c \
	uni.c \
	mplib.c

#
# No need to edit below this line
#

.SUFFIXES:
.SUFFIXES: .c .o

OBJ=	$(SRC:.c=.o)

.c.o:
	$(CC) $(CFLAGS) -c $<

all:	$(EXE)

$(OBJ):	$(INC)

$(EXE):	$(OBJ)
	$(CC) $(CFLAGS) -o $@ $(OBJ) $(LFLAGS)

$(OBJ):	$(MF)

clean:
	rm -f $(OBJ) $(EXE) core
traffic/C-MPI/traffic.h0000664000175000017500000000036112540576420013215 0ustar  dshdsh#define SEED  5743

#include "mplib.h"

int initroad(int *road, int n, float density, int seed);
void updatebcs(int *road, int n);
void mpupdatebcs(int *road, int n, int procid, int nproc);
int updateroad(int *newroad, int *oldroad, int n);
traffic/C-MPI/mplib.h0000664000175000017500000000035312534300645012700 0ustar  dshdshvoid mpstart(int *numproc, int *procid);
void mpstop(void);
void mpgsum(int *ival);
void mpscatter(int *inbuff, int *outbuff, int ndata);
void mpsendrecv(int *senddata, int sendproc, int *recvdata, int recvproc);
double gettime(void);
traffic/F-SER/0000775000175000017500000000000012620656243011436 5ustar  dshdshtraffic/F-SER/unirand.f900000664000175000017500000000335612534300645013421 0ustar  dshdshmodule unirand

  implicit none

  ! Some magic constants

  real, private, parameter :: cd = ( 7654321.0/16777216.0), &
                              cm = (16777213.0/16777216.0) 

  ! The state variables

  integer, private, parameter :: nstate = 97

  real,    private, dimension(nstate) :: u
  integer, private                    :: iu, ju
  real,    private                    :: c
  
  ! The routines

  public  :: uni, rinit
  private :: rstart

contains

  real function uni()

    ! First call rstart(i,j,k,l), with i,j,k,l integers from 1...168 NOT all 1

    uni = u(iu) - u(ju)

    if (uni < 0.0) uni = uni + 1.0

    u(iu) = uni

    iu = iu - 1
    if (iu == 0) iu = nstate

    ju = ju - 1
    if (ju == 0) ju = nstate

    c = c - cd
    if (c < 0.0) c = c + cm

    uni = uni - c
    if (uni < 0.0) uni = uni + 1.0

  end function uni

  !
  !  The algorithm from James RMARIN to generate the 4 seeds from an
  !  integer in the range 0 <= seed <= 900,000,000
  !

  subroutine rinit(ijkl)

    integer :: ijkl, ij, kl, i, j, k, l

    ij = ijkl/30082
    kl = ijkl - 30082*ij

    i  = mod(ij/177,177) + 2
    j  = mod(ij,    177) + 2
    k  = mod(kl/169,178) + 1
    l  = mod(kl,    169)

    call rstart(i, j, k, l)

  end subroutine rinit


  subroutine rstart(i, j, k, l)

    integer :: i, j, k, l
    integer :: ii, jj, m
    real    :: s, t

    do ii = 1, nstate

       s = 0.0
       t = 0.5

       do jj = 1, 24

          m = mod(mod(i*j,179)*k,179)
          i = j
          j = k
          k = m
          l = mod(53*L+1,169)
          if (mod(l*m,64) >= 32) s = s + t
          t = 0.5*t

       end do

       u(ii) = s

    end do

    iu = nstate
    ju = 33
    c  = (362436.0/16777216.0)


  end subroutine rstart

end module unirand
traffic/F-SER/traffic.f900000664000175000017500000000306412540772577013411 0ustar  dshdshprogram traffic

  use trafficlib

  implicit none

  integer, parameter :: NCELL = 100000
  integer, parameter :: maxiter = 200000000/NCELL
  integer, parameter :: printfreq = maxiter/10

  integer :: i, iter, nmove, ncars
  real    :: density

  integer, allocatable, dimension(:) :: newroad, oldroad

  double precision :: tstart, tstop

! Set target density of cars

  density = 0.52

  write(*,*) 'Length of road is ', NCELL
  write(*,*) 'Number of iterations is ', maxiter
  write(*,*) 'Target density of cars is ', density

! Allocate arrays

  allocate(newroad(0:NCELL+1))
  allocate(oldroad(0:NCELL+1))

! Initialise road accordingly using random number generator

  write(*,*) 'Initialising ...'

  ncars = initroad(oldroad(1), NCELL, density, seed)

  write(*,*) '... done'

  write(*,*) 'Actual density of cars is ', float(ncars)/float(NCELL)
  write(*,*)

  tstart = gettime()

  do iter = 1, maxiter

     call updatebcs(oldroad, NCELL)

     nmove = updateroad(newroad, oldroad, NCELL)

! Copy new to old array

     do i = 1, NCELL

        oldroad(i) = newroad(i)

     end do

     if (mod(iter, printfreq) .eq. 0) then

        write(*,*) 'At iteration ', iter, ' average velocity is ', &
             float(nmove)/float(ncars)
     end if

  end do

  tstop = gettime()

  deallocate(newroad)
  deallocate(oldroad)

  write(*,*)
  write(*,*) 'Finished'
  write(*,*)
  write(*,*) 'Time taken was  ', tstop - tstart, ' seconds'
  write(*,*) 'Update rate was ', &
              1.d-6*float(NCELL)*float(maxiter)/(tstop-tstart), ' MCOPs'
  write(*,*)

end program traffic
traffic/F-SER/Makefile0000664000175000017500000000071712620656162013103 0ustar  dshdshMF=	Makefile

# For morar
FC=	pgf90

# For ARCHER
#FC=	ftn

FFLAGS=
LFLAGS=

EXE=	traffic

SRC= \
	traffic.f90 \
	trafficlib.f90 \
	unirand.f90

#
# No need to edit below this line
#

.SUFFIXES:
.SUFFIXES: .f90 .o

OBJ=	$(SRC:.f90=.o)
MOD=	$(SRC:.f90=.mod)

.f90.o:
	$(FC) $(FFLAGS) -c $<

all:	$(EXE)

$(EXE):	$(OBJ)
	$(FC) $(FFLAGS) -o $@ $(OBJ) $(LFLAGS)

$(OBJ):	$(MF)

trafficlib.o:	unirand.o
traffic.o:	trafficlib.o

clean:
	rm -f $(OBJ) $(MOD) $(EXE) core
traffic/F-SER/trafficlib.f900000664000175000017500000000362312540772664014076 0ustar  dshdshmodule trafficlib

  use unirand

  implicit none

  integer, parameter :: seed = 5743

contains

integer function initroad(road, n, density, seedval)

  integer :: n, seedval
  integer :: road(n)
  real    :: density
  
  integer :: i, ncar
  real    :: rng

  call rinit(seedval)

  ncar = 0

  do i = 1, n

     rng = uni()

     if (rng .lt. density) then

        road(i) = 1

     else

        road(i) = 0

     end if

     ncar = ncar + road(i)
    
  end do

  initroad = ncar

end function initroad


integer function updateroad(newroad, oldroad, n)

  integer :: n
  integer :: newroad(0:n+1), oldroad(0:n+1)

  integer :: i, nmove

  nmove = 0

  do i = 1, n

     if (oldroad(i) .eq. 1) then

        if (oldroad(i+1) .eq. 1) then

           newroad(i) = 1
           
        else

           newroad(i) = 0
           nmove = nmove + 1

        end if

     else

        if (oldroad(i-1) .eq. 1) then

           newroad(i) = 1

        else

           newroad(i) = 0
           
        end if

     end if

  end do

  updateroad = nmove

end function updateroad


subroutine updatebcs(road, n)

  integer :: n
  integer :: road(0:n+1)

  road(0)   = road(n)
  road(n+1) = road(1)

end subroutine updatebcs


double precision function gettime()

    logical, save :: firstcall = .true.

    integer, parameter :: int32kind = selected_int_kind( 9)
    integer, parameter :: int64kind = selected_int_kind(18)

    integer, parameter :: intkind = int64kind

    integer(kind = intkind) :: count,rate

    double precision, save :: ticktime

    if (firstcall) then

       firstcall = .false.

       call system_clock(count, rate)

       ticktime = 1.0d0/dble(rate)
       gettime  = dble(count)*ticktime

!          write(*,*) 'Clock resolution is ', ticktime*1.0e6, ', usecs'

    else

       call system_clock(count)

       gettime = dble(count)*ticktime

    end if

  end function gettime

end module trafficlib
traffic/F-MPI/0000775000175000017500000000000012633773613011437 5ustar  dshdshtraffic/F-MPI/unirand.f900000664000175000017500000000335612536040135013412 0ustar  dshdshmodule unirand

  implicit none

  ! Some magic constants

  real, private, parameter :: cd = ( 7654321.0/16777216.0), &
                              cm = (16777213.0/16777216.0) 

  ! The state variables

  integer, private, parameter :: nstate = 97

  real,    private, dimension(nstate) :: u
  integer, private                    :: iu, ju
  real,    private                    :: c
  
  ! The routines

  public  :: uni, rinit
  private :: rstart

contains

  real function uni()

    ! First call rstart(i,j,k,l), with i,j,k,l integers from 1...168 NOT all 1

    uni = u(iu) - u(ju)

    if (uni < 0.0) uni = uni + 1.0

    u(iu) = uni

    iu = iu - 1
    if (iu == 0) iu = nstate

    ju = ju - 1
    if (ju == 0) ju = nstate

    c = c - cd
    if (c < 0.0) c = c + cm

    uni = uni - c
    if (uni < 0.0) uni = uni + 1.0

  end function uni

  !
  !  The algorithm from James RMARIN to generate the 4 seeds from an
  !  integer in the range 0 <= seed <= 900,000,000
  !

  subroutine rinit(ijkl)

    integer :: ijkl, ij, kl, i, j, k, l

    ij = ijkl/30082
    kl = ijkl - 30082*ij

    i  = mod(ij/177,177) + 2
    j  = mod(ij,    177) + 2
    k  = mod(kl/169,178) + 1
    l  = mod(kl,    169)

    call rstart(i, j, k, l)

  end subroutine rinit


  subroutine rstart(i, j, k, l)

    integer :: i, j, k, l
    integer :: ii, jj, m
    real    :: s, t

    do ii = 1, nstate

       s = 0.0
       t = 0.5

       do jj = 1, 24

          m = mod(mod(i*j,179)*k,179)
          i = j
          j = k
          k = m
          l = mod(53*L+1,169)
          if (mod(l*m,64) >= 32) s = s + t
          t = 0.5*t

       end do

       u(ii) = s

    end do

    iu = nstate
    ju = 33
    c  = (362436.0/16777216.0)


  end subroutine rstart

end module unirand
traffic/F-MPI/traffic.f900000664000175000017500000000427512540772577013412 0ustar  dshdshprogram traffic

  use trafficlib

  implicit none

  integer, parameter :: NCELL = 100000
  integer, parameter :: maxiter = 200000000/NCELL
  integer, parameter :: printfreq = maxiter/10

  integer procid, nproc, nlocal

  integer :: i, iter, nmove, ncars
  real    :: density

  integer, allocatable, dimension(:) :: newroad, oldroad, bigroad

  double precision :: tstart, tstop

! Set target density of cars

  density = 0.52

! Start message passing system

  call mpstart(nproc, procid)

  nlocal = NCELL/nproc

  if (procid == 0) then

     write(*,*) 'Length of road is ', NCELL
     write(*,*) 'Number of iterations is ', maxiter
     write(*,*) 'Target density of cars is ', density
     write(*,*) 'Running on ', nproc, ' processes'

  end if

! Allocate arrays

  allocate(newroad(0:nlocal+1))
  allocate(oldroad(0:nlocal+1))

  if (procid == 0) then

     allocate(bigroad(NCELL))

     ! Initialise road accordingly using random number generator

     write(*,*) 'Initialising ...'

     ncars = initroad(bigroad, NCELL, density, seed)

     write(*,*) '... done'

     write(*,*) 'Actual density of cars is ', float(ncars)/float(NCELL)
     write(*,*) 'Scattering data ...'

  end if

  call mpscatter(bigroad, oldroad(1:nlocal), nlocal)

  if (procid == 0) then

     write(*,*) '... done'
     write(*,*)

  end if

  tstart = gettime()

  do iter = 1, maxiter

     call mpupdatebcs(oldroad, nlocal, procid, nproc)

     nmove = updateroad(newroad, oldroad, nlocal)

! Globally sum the value

     call mpgsum(nmove)

! Copy new to old array

     do i = 1, nlocal

        oldroad(i) = newroad(i)

     end do

     if (mod(iter, printfreq) == 0) then

        if (procid == 0) write(*,*) 'At iteration ', iter, &
             ' average velocity is ', float(nmove)/float(ncars)
     end if

  end do

  tstop = gettime()

  deallocate(oldroad)
  deallocate(newroad)

  if (procid == 0) then

     deallocate(bigroad)

     write(*,*)
     write(*,*) 'Finished'
     write(*,*)
     write(*,*) 'Time taken was  ', tstop - tstart, ' seconds'
     write(*,*) 'Update rate was ', &
                 1.d-6*float(NCELL)*float(maxiter)/(tstop-tstart), ' MCOPs'
     write(*,*)

  end if

  call mpstop()

end program traffic
traffic/F-MPI/Makefile0000664000175000017500000000076312620655641013101 0ustar  dshdshMF=	Makefile

# For morar
FC=	mpif90

# For ARCHER
#FC=	ftn

FFLAGS=
LFLAGS=

EXE=	traffic

SRC= \
	traffic.f90 \
	trafficlib.f90 \
	unirand.f90 \
	mplib.f90

#
# No need to edit below this line
#

.SUFFIXES:
.SUFFIXES: .f90 .o

OBJ=	$(SRC:.f90=.o)
MOD=	$(SRC:.f90=.mod)

.f90.o:
	$(FC) $(FFLAGS) -c $<

all:	$(EXE)

$(EXE):	$(OBJ)
	$(FC) $(FFLAGS) -o $@ $(OBJ) $(LFLAGS)

$(OBJ):	$(MF)

traffic.o:	trafficlib.o
trafficlib.o:	mplib.o
trafficlib.o:	unirand.o

clean:
	rm -f $(OBJ) $(MOD) $(EXE) core
traffic/F-MPI/mplib.f900000664000175000017500000000247512536062312013057 0ustar  dshdshmodule mplib

  use mpi

  implicit none

  integer :: ierror

contains

subroutine mpstart(numproc, procid)

  integer :: numproc, procid

  call MPI_Init(ierror)
  call MPI_Comm_rank(MPI_COMM_WORLD, procid, ierror)
  call MPI_Comm_size(MPI_COMM_WORLD, numproc, ierror)

end subroutine mpstart

subroutine mpstop

  call MPI_Finalize(ierror)

end subroutine mpstop
 
subroutine mpgsum(ival)

  integer :: ival, itmp
 
  call MPI_Allreduce(ival, itmp, 1, MPI_INTEGER, MPI_SUM, &
                     MPI_COMM_WORLD, ierror)

  ival = itmp

end subroutine mpgsum

subroutine mpscatter(inbuff, outbuff, ndata)

  integer :: ndata

  integer, dimension(ndata) :: outbuff
  integer, dimension(:)     :: inbuff

  call MPI_Scatter(inbuff,  ndata, MPI_INTEGER, &
	           outbuff, ndata, MPI_INTEGER, &
                   0, MPI_COMM_WORLD, ierror)

end subroutine mpscatter

subroutine mpsendrecv(senddata, sendproc, recvdata, recvproc)

  integer :: senddata, sendproc, recvdata, recvproc
  integer, dimension(MPI_STATUS_SIZE) :: status

  call MPI_Sendrecv(senddata, 1, MPI_INTEGER, sendproc, 1, &
                    recvdata, 1, MPI_INTEGER, recvproc, 1, &
                    MPI_COMM_WORLD, status, ierror)
 
end subroutine mpsendrecv

double precision function gettime()

  gettime = MPI_Wtime()

end function gettime

end module mplib
traffic/F-MPI/trafficlib.f900000664000175000017500000000312712540772676014074 0ustar  dshdshmodule trafficlib

  use unirand
  use mplib

  implicit none

  integer, parameter :: seed = 5743

contains

integer function initroad(road, n, density, seedval)

  integer :: n, seedval
  integer :: road(n)
  real    :: density
  
  integer :: i, ncar
  real    :: rng

  call rinit(seedval)

  ncar = 0

  do i = 1, n

     rng = uni()

     if (rng .lt. density) then

        road(i) = 1

     else

        road(i) = 0

     end if

     ncar = ncar + road(i)
    
  end do

  initroad = ncar

end function initroad


integer function updateroad(newroad, oldroad, n)

  integer :: n
  integer :: newroad(0:n+1), oldroad(0:n+1)

  integer :: i, nmove

  nmove = 0

  do i = 1, n

     if (oldroad(i) .eq. 1) then

        if (oldroad(i+1) .eq. 1) then

           newroad(i) = 1
           
        else

           newroad(i) = 0
           nmove = nmove + 1

        end if

     else

        if (oldroad(i-1) .eq. 1) then

           newroad(i) = 1

        else

           newroad(i) = 0
           
        end if

     end if

  end do

  updateroad = nmove

end function updateroad


subroutine updatebcs(road, n)

  integer :: n
  integer :: road(0:n+1)

  road(0)   = road(n)
  road(n+1) = road(1)

end subroutine updatebcs


subroutine mpupdatebcs(road, n, procid, nproc)

  integer :: n, procid, nproc
  integer :: road(0:n+1)

  integer procup, procdown

  procup   = mod(procid + 1, nproc)
  procdown = mod(procid + nproc - 1, nproc)

  call mpsendrecv(road(n), procup,   road(0),   procdown)
  call mpsendrecv(road(1), procdown, road(n+1), procup  )

end subroutine mpupdatebcs

end module trafficlib