percolate/C-MPI/percolate.c
percolate/C-MPI/unirand.c
percolate/C-MPI/percolate.h
percolate/C-MPI/Makefile
percolate/C-MPI/percio.c
percolate/F-MPI/percmod.f90
percolate/F-MPI/percolate.f90
percolate/F-MPI/percio.f90
percolate/F-MPI/Makefile
percolate/F-MPI/unirand.f90
/*
* Simple parallel program to test for percolation of a cluster
*/
#include
#include
#include
#include “percolate.h”
/*
* Simple parallel program to test for percolation of a cluster.
*/
int main(int argc, char *argv[])
{
/*
* Define the main arrays for the simulation
*/
int old[M+2][N+2], new[M+2][N+2];
/*
* Additional array WITHOUT halos for initialisation and IO. This
* is of size LxL because, even in our parallel program, we do
* these two steps in serial
*/
int map[L][L];
/*
* Array to store local part of map
*/
int smallmap[M][N];
/*
* Variables that define the simulation
*/
int seed;
double rho;
/*
* Local variables
*/
int i, j, nhole, step, maxstep, oldval, newval;
int nchangelocal, nchange, printfreq;
int itop, ibot, perc;
double r;
/*
* MPI variables
*/
MPI_Comm comm = MPI_COMM_WORLD;
MPI_Status status;
int size, rank, prev, next;
int tag = 1;
MPI_Init(&argc, &argv);
MPI_Comm_size(comm, &size);
MPI_Comm_rank(comm, &rank);
next = rank + 1;
prev = rank – 1;
/*
* Non-periodic boundary conditions
*
* Note that the special rank of MPI_PROC_NULL is a “black hole” for
* communications. Using this value for processes off the edges of the
* image means there is no additional logic needed to ensure processes
* at the edges do not attempt to send to or receive from invalid
* ranks (i.e. rank = -1 and rank = NPROC).
*
* Proper solution would compute neighbours with a Cartesian topology
* and MPI_Cart_shift, where MPI_PROC_NULL is assigned automatically.
*/
if (next >= size)
{
next = MPI_PROC_NULL;
}
if (prev < 0)
{
prev = MPI_PROC_NULL;
}
if (NPROC != size)
{
if (rank == 0)
{
printf("percolate: ERROR, NPROC = %d but running on %d\n",
NPROC, size);
}
MPI_Finalize();
return 0;
}
if (argc != 2)
{
if (rank == 0)
{
printf("Usage: percolate
}
MPI_Finalize();
return 0;
}
/*
* Update for a fixed number of steps, periodically report progress
*/
maxstep = 5*L;
printfreq = 100;
if (rank == 0)
{
printf(“percolate: running on %d process(es)\n”, size);
/*
* Set most important value: the rock density rho (between 0 and 1)
*/
rho = 0.4064;
/*
* Set the randum number seed and initialise the generator
*/
seed = atoi(argv[1]);
printf(“percolate: L = %d, rho = %f, seed = %d, maxstep = %d\n”,
L, rho, seed, maxstep);
rinit(seed);
/*
* Initialise map with density rho. Zero indicates rock, a positive
* value indicates a hole. For the algorithm to work, all the holes
* must be initialised with a unique integer
*/
nhole = 0;
for (i=0; i < L; i++)
{
for (j=0; j < L; j++)
{
r=uni();
if(r < rho)
{
map[i][j] = 0;
}
else
{
nhole++;
map[i][j] = nhole;
}
}
}
printf("percolate: rho = %f, actual density = %f\n",
rho, 1.0 - ((double) nhole)/((double) L*L) );
}
/*
* Now scatter map to smallmap
*/
MPI_Scatter(map, M*N, MPI_INT, smallmap, M*N, MPI_INT, 0, comm);
/*
* Initialise the old array: copy the LxL array smallmap to the centre of
* old, and set the halo values to zero.
*/
for (i=1; i <= M; i++)
{
for (j=1; j <= N; j++)
{
old[i][j] = smallmap[i-1][j-1];
}
}
for (i=0; i <= M+1; i++) // zero the bottom and top halos
{
old[i][0] = 0;
old[i][N+1] = 0;
}
for (j=0; j <= N+1; j++) // zero the left and right halos
{
old[0][j] = 0;
old[M+1][j] = 0;
}
step = 1;
nchange = 1;
while (step <= maxstep)
{
/*
* Swap halos up and down
*/
/*
* Communications is done using the sendrecv routine; a proper
* solution would use non-blocking communications (e.g. some
* combination of issend/recv or ssend/irecv)
*/
MPI_Sendrecv(&old[M][1], N, MPI_INT, next, tag,
&old[0][1], N, MPI_INT, prev, tag,
comm, &status);
MPI_Sendrecv(&old[1][1], N, MPI_INT, prev, tag,
&old[M+1][1], N, MPI_INT, next, tag,
comm, &status);
nchangelocal = 0;
for (i=1; i<=M; i++)
{
for (j=1; j<=N; j++)
{
oldval = old[i][j];
newval = oldval;
/*
* Set new[i][j] to be the maximum value of old[i][j]
* and its four nearest neighbours
*/
if (oldval != 0)
{
if (old[i][j-1] > newval) newval = old[i][j-1];
if (old[i][j+1] > newval) newval = old[i][j+1];
if (old[i-1][j] > newval) newval = old[i-1][j];
if (old[i+1][j] > newval) newval = old[i+1][j];
if (newval != oldval)
{
++nchangelocal;
}
}
new[i][j] = newval;
}
}
/*
* Compute global number of changes on rank 0
*/
MPI_Reduce(&nchangelocal, &nchange, 1, MPI_INT, MPI_SUM, 0, comm);
/*
* Report progress every now and then
*/
if (step % printfreq == 0)
{
if (rank == 0)
{
printf(“percolate: changes on step %d is %d\n”,
step, nchange);
}
}
/*
* Copy back in preparation for next step, omitting halos
*/
for (i=1; i<=M; i++)
{
for (j=1; j<=N; j++)
{
old[i][j] = new[i][j];
}
}
step++;
}
/*
* We set a maximum number of steps to ensure the algorithm always
* terminates. However, if we hit this limit before the algorithm
* has finished then there must have been a problem (e.g. the value
* of maxstep is too small)
*/
if (rank == 0)
{
if (nchange != 0)
{
printf("percolate: WARNING max steps = %d reached but nchange != 0\n",
maxstep);
}
}
/*
* Copy the centre of old, excluding the halos, into smallmap
*/
for (i=1; i<=M; i++)
{
for (j=1; j<=N; j++)
{
smallmap[i-1][j-1] = old[i][j];
}
}
/*
* Now gather smallmap back to map
*/
MPI_Gather(smallmap, M*N, MPI_INT, map, M*N, MPI_INT, 0, comm);
/*
* Test to see if percolation occurred by looking for positive numbers
* that appear on both the top and bottom edges
*/
if (rank == 0)
{
perc = 0;
for (itop=0; itop < L; itop++)
{
if (map[itop][L-1] > 0)
{
for (ibot=0; ibot < L; ibot++)
{
if (map[ibot][0] == map[itop][L-1])
{
perc = 1;
}
}
}
}
if (perc != 0)
{
printf("percolate: cluster DOES percolate\n");
}
else
{
printf("percolate: cluster DOES NOT percolate\n");
}
/*
* Write the map to the file "map.pgm", displaying the two
* largest clusters. If the last argument here was 3, it would
* display the three largest clusters etc. The picture looks
* cleanest with only a single cluster, but multiple clusters
* are useful for debugging.
*/
mapwrite("map.pgm", map, 2);
}
MPI_Finalize();
return 0;
}
/*
* 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)
*/
#include
#include
/*
* Global variables for rstart & uni
*/
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);
}
/*
* Main header file for percolation code.
*/
/*
* System size L
*/
#define L 480
/*
* Use 1D decomposition over NPROC processes across first dimension
* For an LxL simulation, the local arrays are of size MxN
*/
#define NPROC 4
#define M L/NPROC
#define N L
/*
* Prototypes for supplied functions
*/
/*
* Visualisation
*/
void mapwrite(char *percfile, int map[L][L], int ncluster);
void mapwritedynamic(char *percfile, int **map, int l, int ncluster);
/*
* Random numbers
*/
void rinit(int ijkl);
float uni(void);
MF= Makefile
CC= mpicc
CFLAGS= -cc=icc -O3 -Wall
LFLAGS= $(CFLAGS)
EXE= percolate
INC= \
percolate.h
SRC= \
percolate.c \
percio.c \
unirand.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) $(LFLAGS) -o $@ $(OBJ)
$(OBJ): $(MF)
clean:
rm -f $(EXE) $(OBJ) core
#include
#include
#include “percolate.h”
/*
* Function to write a percolation map in greyscale Portable Grey Map
* (PGM) format. The largest “ncluster” clusters are identified and
* shown as shades of grey against a black background, with the
* largest cluster shown in white.
*/
#define MAXNCLUSTER 9 // Must be able to identify by a single digit
int foundcluster[MAXNCLUSTER];
#define MAXCLUSTERID (L*L) // Extreme case of zero density
int clustersize[MAXCLUSTERID+1];
void mapwrite(char *percfile, int map[L][L], int ncluster)
{
FILE *fp;
int i, j, colour, npix;
int clusterid, icluster, maxcluster, prevcluster;
static int pixperline = 32; // PGM format limits to 70 characters per line
if (ncluster > MAXNCLUSTER)
{
printf(“mapwrite: WARNING ncluster too large, resetting to %d\n”,
MAXNCLUSTER);
ncluster = MAXNCLUSTER;
}
if (ncluster > 1)
{
printf(“mapwrite: visualising the largest %d clusters\n”, ncluster);
}
else
{
printf(“mapwrite: only visualising the largest cluster\n”);
}
/*
* Count up the size of each cluster
*/
for (i=0; i<=MAXCLUSTERID; i++)
{
clustersize[i] = 0;
}
for (i=0; i < L; i++)
{
for (j=0; j < L; j++)
{
clusterid = map[i][j];
if (clusterid > 0)
{
clustersize[clusterid]++;
}
}
}
/*
* Find the size of the “ncluster” largest clusters (by brute force!)
*/
prevcluster = MAXCLUSTERID+1; // Larger than the largest possible cluster id
for (icluster = 0; icluster < ncluster; icluster++)
{
maxcluster = 0;
for (i=0; i<=MAXCLUSTERID; i++)
{
if (clustersize[i] > maxcluster && clustersize[i] < prevcluster)
{
maxcluster = clustersize[i];
}
}
foundcluster[icluster] = maxcluster;
prevcluster = maxcluster;
}
if (ncluster > 1)
{
printf(“mapwrite: cluster sizes are “);
}
else
{
printf(“mapwrite: maximum cluster size is “);
}
for (icluster = 0; icluster < ncluster-1; icluster++) { printf("%d, ", foundcluster[icluster]); } printf("%d\n", foundcluster[ncluster-1]); /* * Write the file */ printf("mapwrite: opening file <%s>\n”, percfile);
fp = fopen(percfile, “w”);
printf(“mapwrite: writing data …\n”);
/*
* Start with the PGM header
*/
fprintf(fp, “P2\n”);
fprintf(fp, “%d %d\n%d\n”, L, L, ncluster);
/*
* Now write the cells to file so that map[0][0] is in the
* bottom-left-hand corner and map[l-1][l-1] is in the
* top-right-hand corner
*/
npix = 0;
for (j=L-1; j >= 0; j–)
{
for (i=0; i < L; i++)
{
clusterid = map[i][j];
/*
* Write out the largest cluster(s), shading appropriately
*/
colour = 0;
if (clusterid > 0)
{
for (icluster = 0; icluster < ncluster; icluster++)
{
if (clustersize[clusterid] == foundcluster[icluster])
{
// Largest (first) cluster is white
colour = ncluster - icluster;
}
}
}
npix++;
// Make sure lines wrap after "npix" pixels
if (npix == 1)
{
fprintf(fp, "%1d", colour);
}
else if (npix < pixperline)
{
fprintf(fp, " %1d", colour);
}
else
{
fprintf(fp, " %1d\n", colour);
npix = 0;
}
}
}
if (npix != 0) fprintf(fp, "\n");
printf("mapwrite: ... done\n");
fclose(fp);
printf("mapwrite: file closed\n");
}
/*
* Function to write a percolation map in greyscale Portable Grey Map
* (PGM) format. The largest "ncluster" clusters are identified and
* shown as shades of grey against a black background, with the
* largest cluster shown in white.
*
* Note that this version expects the map array to have been
* dynamically allocated, e.g. using the arralloc() routine:
*
* int **map;
* map = (int **) arralloc(sizeof(int), 2, L, L);
* ...
* mapwritedynamic("map.pgm", map, L, 1);
*/
#define MAXNCLUSTER 9 // Must be able to identify by a single digit
int foundcluster[MAXNCLUSTER];
void mapwritedynamic(char *percfile, int **map, int l, int ncluster)
{
FILE *fp;
int i, j, colour, npix;
int clusterid, icluster, maxcluster, prevcluster, maxclusterid;
int *clustersize;
static int pixperline = 32; // PGM format limits to 70 characters per line
if (ncluster > MAXNCLUSTER)
{
printf(“mapwrite: WARNING ncluster too large, resetting to %d\n”,
MAXNCLUSTER);
ncluster = MAXNCLUSTER;
}
if (ncluster > 1)
{
printf(“mapwrite: visualising the largest %d clusters\n”, ncluster);
}
else
{
printf(“mapwrite: only visualising the largest cluster\n”);
}
/*
* Allocate the local clustersize array
*/
maxclusterid = l*l;
if ((clustersize = (int *) malloc((maxclusterid+1)*sizeof(int))) == NULL)
{
printf(“mapwrite: allocation of clustersize failed\n”);
exit(1);
}
/*
* Count up the size of each cluster
*/
for (i=0; i<=maxclusterid; i++)
{
clustersize[i] = 0;
}
for (i=0; i < l; i++)
{
for (j=0; j < l; j++)
{
clusterid = map[i][j];
if (clusterid > 0)
{
clustersize[clusterid]++;
}
}
}
/*
* Find the size of the “ncluster” largest clusters (by brute force!)
*/
prevcluster = maxclusterid+1; // Larger than the largest possible cluster id
for (icluster = 0; icluster < ncluster; icluster++)
{
maxcluster = 0;
for (i=0; i<=maxclusterid; i++)
{
if (clustersize[i] > maxcluster && clustersize[i] < prevcluster)
{
maxcluster = clustersize[i];
}
}
foundcluster[icluster] = maxcluster;
prevcluster = maxcluster;
}
if (ncluster > 1)
{
printf(“mapwrite: cluster sizes are “);
}
else
{
printf(“mapwrite: maximum cluster size is “);
}
for (icluster = 0; icluster < ncluster-1; icluster++) { printf("%d, ", foundcluster[icluster]); } printf("%d\n", foundcluster[ncluster-1]); /* * Write the file */ printf("mapwrite: opening file <%s>\n”, percfile);
fp = fopen(percfile, “w”);
printf(“mapwrite: writing data …\n”);
/*
* Start with the PGM header
*/
fprintf(fp, “P2\n”);
fprintf(fp, “%d %d\n%d\n”, l, l, ncluster);
/*
* Now write the cells to file so that map[0][0] is in the
* bottom-left-hand corner and map[l-1][l-1] is in the
* top-right-hand corner
*/
npix = 0;
for (j=l-1; j >= 0; j–)
{
for (i=0; i < l; i++)
{
clusterid = map[i][j];
/*
* Write out the largest cluster(s), shading appropriately
*/
colour = 0;
if (clusterid > 0)
{
for (icluster = 0; icluster < ncluster; icluster++)
{
if (clustersize[clusterid] == foundcluster[icluster])
{
// Largest (first) cluster is white
colour = ncluster - icluster;
}
}
}
npix++;
// Make sure lines wrap after "npix" pixels
if (npix == 1)
{
fprintf(fp, "%1d", colour);
}
else if (npix < pixperline)
{
fprintf(fp, " %1d", colour);
}
else
{
fprintf(fp, " %1d\n", colour);
npix = 0;
}
}
}
if (npix != 0) fprintf(fp, "\n");
printf("mapwrite: ... done\n");
fclose(fp);
printf("mapwrite: file closed\n");
/*
* De-allocate the local clustersize array
*/
free(clustersize);
}
!
! Main module file for percolation code.
!
module percmod
implicit none
!
! System size L
!
integer, parameter :: L = 480
!
! Use 1D decomposition over nproc processes across second dimension
! For an LxL simulation, the local arrays are of size MxN
!
integer, parameter :: nproc = 4
integer, parameter :: M = L
integer, parameter :: N = L/nproc
end module percmod
program percolate
!
! Simple parallel program to test for percolation of a cluster
!
use mpi
use percmod
use unirand
use percio
!
! Define the main arrays for the simulation
!
integer, dimension(0:M+1,0:N+1) :: old, new
!
! Additional array WITHOUT halos for initialisation and IO. This
! is of size LxL because, even in our parallel program, we do
! these two steps in serial
!
integer, dimension(L,L) :: map
!
! Array to store local part of map
!
integer, dimension(M,N) :: smallmap
!
! Variables that define the simulation
!
double precision :: rho
integer :: seed
!
! Local variables
!
integer :: argc, i, j, nhole, step, maxstep
integer :: nchange, nchangelocal, printfreq
integer :: itop, ibot
logical :: perc
character*(16) :: tmparg
double precision :: r
!
! MPI variables
!
integer :: comm = MPI_COMM_WORLD
integer, dimension(MPI_STATUS_SIZE) :: status
integer :: size, rank, prev, next, ierr
integer :: tag = 1
call MPI_Init(ierr)
call MPI_Comm_size(comm, size, ierr)
call MPI_Comm_rank(comm, rank, ierr)
next = rank + 1
prev = rank - 1
!
! Non-periodic boundary conditions
!
! Note that the special rank of MPI_PROC_NULL is a "black hole" for
! communications. Using this value for processes off the edges of the
! image means there is no additional logic needed to ensure processes
! at the edges do not attempt to send to or receive from invalid
! ranks (i.e. rank = -1 and rank = NPROC).
!
! Proper solution would compute neighbours with a Cartesian topology
! and MPI_Cart_shift, where MPI_PROC_NULL is assigned automatically.
!
if (next >= size) then
next = MPI_PROC_NULL
end if
if (prev < 0) then
prev = MPI_PROC_NULL
end if
if (nproc /= size) then
if (rank == 0) then
write(*,*) "percolate: ERROR, NPROC = ",nproc," but running on ",size
end if
call MPI_Finalize(ierr)
stop
end if
argc = command_argument_count()
if (argc /= 1) then
if (rank == 0) then
write(*,*) "Usage: percolate
end if
call MPI_Finalize(ierr)
stop
end if
!
! Update for a fixed number of iterations
!
maxstep = 5*L
printfreq = 100
if (rank == 0) then
write(*,*) “percolate: running on “, size, ” process(es)”
!
! Set most important value: the rock density rho (between 0 and 1)
!
rho = 0.4064
!
! Set the randum number seed and initialise the generator
!
call get_command_argument(1, tmparg)
read(tmparg,*) seed
write(*,fmt=”(“” percolate: params are L = “”,i5,””, rho = “”,f6.4, &
&””, seed = “”, i5,””, maxstep = “”,i5)”) L, rho, seed, maxstep
call rinit(seed)
!
! Initialise map with density rho. Zero indicates rock, a positive
! value indicates a hole. For the algorithm to work, all the holes
! must be initialised with a unique integer.
!
! Difficult to do with array syntax due to non-pure nature of random
! number generation routine uni().
!
nhole = 0
do i = 1, L
do j = 1, L
r = uni()
if (r < rho) then map(i,j) = 0 else nhole = nhole + 1 map(i,j) = nhole end if end do end do write(*,fmt="("" percolate: rho = "",f8.6,"", actual density = "",f8.6)")& rho, 1.0-float(nhole)/float(L*L) end if ! ! Now scatter map to smallmap ! call MPI_Scatter(map, M*N, MPI_INTEGER, & smallmap, M*N, MPI_INTEGER, 0, comm, ierr); ! ! Initialise the old array: copy the LxL array smallmap to the centre of ! old, and set the halo values to zero. ! old(1:M,1:N) = smallmap(:,:) ! Zero the bottom and top halos old(:,0) = 0 old(:,N+1) = 0 ! Zero the left and right halos old(0,:) = 0 old(M+1,0) = 0 step = 1 nchange = 1 do while (step <= maxstep) ! ! Swap halos up and down ! ! ! Communications is done using the sendrecv routine; a proper ! solution would use non-blocking communications (e.g. some ! combination of issend/recv or ssend/irecv) ! call MPI_Sendrecv(old(1,N), M, MPI_INTEGER, next, tag, & old(1,0), M, MPI_INTEGER, prev, tag, & comm, status, ierr) call MPI_Sendrecv(old(1,1), M, MPI_INTEGER, prev, tag, & old(1,N+1), M, MPI_INTEGER, next, tag, & comm, status, ierr) ! ! Set new(i,j) to be the maximum value of old(i,j) ! and its four nearest neighbours ! where(old(1:M,1:N) /= 0) & new(1:M,1:N) = max(old(1:M,1:N), old(1:M,0:N-1), & old(1:M,2:N+1), & old(0:M-1,1:N), & old(2:M+1,1:N) ) nchangelocal = count(new(1:M,1:N) /= old(1:M,1:N)) ! ! Compute global number of changes on rank 0 ! call MPI_Reduce(nchangelocal, nchange, 1, MPI_INTEGER, MPI_SUM, & 0, comm, ierr); ! ! Report progress every now and then ! if (mod(step, printfreq) == 0 ) then if (rank == 0) then write(*,*) "percolate: changes on step ", step, " is ", nchange end if end if ! ! Copy back in preparation for next step, omitting halos ! old(1:M,1:N) = new(1:M,1:N) step = step + 1 end do ! ! We set a maximum number of steps to ensure the algorithm always ! terminates. However, if we hit this limit before the algorithm ! has finished then there must have been a problem (e.g. the value ! of maxstep is too small) ! if (rank == 0) then if (nchange /= 0) then write(*,*) "percolate: WARNING max steps = ", maxstep, & " reached but nchange /= 0" end if end if ! ! Copy the centre of old, excluding the halos, into smallmap ! smallmap(:,:) = old(1:M, 1:N) ! ! Now gather smallmap back to map ! call MPI_Gather(smallmap, M*N, MPI_INTEGER, & map, M*N, MPI_INTEGER, 0, comm, ierr); ! ! Test to see if percolation occurred by looking for positive numbers ! that appear on both the top and bottom edges ! if (rank == 0) then perc = .FALSE. do itop = 1, L if (map(itop,L) > 0) then
do ibot = 1, L
if (map(ibot,1) == map(itop,L)) then
perc = .TRUE.
end if
end do
end if
end do
if (perc) then
write(*,*) “percolate: cluster DOES percolate”
else
write(*,*) “percolate: cluster DOES NOT percolate”
end if
!
! Write the map to the file “map.pgm”, displaying the two
! largest clusters. If the last argument here was 3, it would
! display the three largest clusters etc. The picture looks
! cleanest with only a single cluster, but multiple clusters
! are useful for debugging.
!
call mapwrite(“map.pgm”, map, 2)
end if
call MPI_Finalize(ierr)
end program percolate
module percio
implicit none
public :: mapwrite
contains
! Function to write a percolation map in greyscale Portable Grey
! Map (PGM) format. The largest “ncluster” clusters are identified
! and shown as shades of grey against a black background, with the
! largest cluster shown in white.
subroutine mapwrite(percfile, map, ncluster)
character (len = *), intent(in) :: percfile
integer, dimension(:,:), intent(in) :: map
integer, intent(in) :: ncluster
integer, parameter :: maxncluster = 9 ! Must be a single digit
integer, parameter :: pixperline = 32 ! PGM limit 70 chars per line
integer, dimension(pixperline) :: pgmline
integer, dimension(maxncluster) :: foundcluster
integer, allocatable, dimension(:) :: clustersize
integer, parameter :: fmtlen = 64
character (len = fmtlen) :: fmtstring
integer, parameter :: iounit = 12
integer :: m, n
integer :: i, j, npix, colour
integer :: clusterid, icluster, lncluster, prevcluster, maxcluster
m = size(map, 1)
n = size(map, 2)
allocate(clustersize(m*n))
lncluster = ncluster
if (lncluster > maxncluster) then
write(*,*) “mapwrite: WARNING ncluster too large, resetting to “, &
maxncluster
lncluster = maxncluster
end if
if (lncluster > 1) then
write(*,*) “mapwrite: visualising the largest “, lncluster, “clusters”
else
write(*,*) “mapwrite: only visualising the largest cluster”
end if
! Count up the size of each cluster
clustersize(:) = 0
do i = 1, m
do j = 1, n
clusterid = map(i,j)
if (clusterid > 0) then
clustersize(clusterid) = clustersize(clusterid) + 1
end if
end do
end do
! Find the size of the “lncluster” largest clusters (by brute force!)
prevcluster = m*n+1 ! Larger than the largest possible cluster id
do icluster = 1, lncluster
maxcluster = 0
do i = 1, m*n
if (clustersize(i) > maxcluster .and. &
clustersize(i) < prevcluster ) then
maxcluster = clustersize(i)
end if
end do
foundcluster(icluster) = maxcluster
prevcluster = maxcluster
end do
if (lncluster > 1) then
write(*,*) “mapwrite: cluster sizes are “, foundcluster(1:lncluster)
else
write(*,*) “mapwrite: maximum cluster size is “, foundcluster(1)
end if
write(*,*) ‘mapwrite: opening file ‘, percfile
open(unit=iounit, file=percfile)
write(*,*) ‘mapwrite: writing data …’
write(fmtstring, fmt='(”(”, ”i1,”,i5,”(1x, i1))”)’) pixperline-1
write(iounit,fmt='(”P2”)’)
write(iounit,*) m, n
write(iounit,*) lncluster
npix = 0
do j = n, 1, -1
do i = 1, m
clusterid = map(i,j)
! Write out the largest cluster(s), shading appropriately
colour = 0
if (clusterid > 0) then
do icluster = 1, lncluster
if (clustersize(clusterid) == foundcluster(icluster)) then
! Largest (first) cluster is white
colour = lncluster – icluster + 1
end if
end do
end if
npix = npix + 1
pgmline(npix) = colour
if (npix == pixperline) then
write(iounit,fmt=fmtstring) pgmline(1:npix)
npix = 0
end if
end do
end do
if (npix /= 0) then
write(iounit,fmt=fmtstring) pgmline(1:npix)
end if
write(*,*) ‘mapwrite: … done’
close(iounit)
write(*,*) ‘mapwrite: file closed’
return
end subroutine mapwrite
end module percio
MF= Makefile
FC= mpif90
FFLAGS= -fc=ifort -O3 -Warn all
LFLAGS= $(FFLAGS)
EXE= percolate
MOD= \
percmod.f90 \
percio.f90 \
unirand.f90
SRC= percolate.f90
#
# No need to edit below this line
#
.SUFFIXES:
.SUFFIXES: .f90 .o
OBJ= $(MOD:.f90=.o) $(SRC:.f90=.o)
TMP= $(MOD:.f90=.mod)
.f90.o:
$(FC) $(FFLAGS) -c $<
all: $(EXE)
$(OBJ): $(MOD)
$(EXE): $(OBJ)
$(FC) $(FFLAGS) -o $@ $(OBJ)
$(OBJ): $(MF)
clean:
rm -f $(EXE) $(OBJ) $(TMP) core
module 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