CS计算机代考程序代写 Fortran algorithm percolate/C-MPI/percolate.c

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 \n”);
}

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