percolate/C-SER/percolate.c
percolate/C-SER/arralloc.c
percolate/C-SER/percolatedynamic.c
percolate/C-SER/unirand.c
percolate/C-SER/percolate.h
percolate/C-SER/Makefile
percolate/C-SER/percio.c
percolate/C-SER/arralloc.h
percolate/F-SER/percmod.f90
percolate/F-SER/percolate.f90
percolate/F-SER/percio.f90
percolate/F-SER/Makefile
percolate/F-SER/unirand.f90
#include
#include
#include “percolate.h”
/*
* Simple serial 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];
/*
* Variables that define the simulation
*/
int seed;
double rho;
/*
* Local variables
*/
int i, j, nhole, step, maxstep, oldval, newval, nchange, printfreq;
int itop, ibot, perc;
double r;
if (argc != 2)
{
printf(“Usage: percolate
return 0;
}
/*
* Update for a fixed number of steps, periodically report progress
*/
maxstep = 5*L;
printfreq = 100;
/*
* 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) );
/*
* Initialise the old array: copy the LxL array map 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] = map[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)
{
nchange = 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)
{
++nchange;
}
}
new[i][j] = newval;
}
}
/*
* Report progress every now and then
*/
if (step % printfreq == 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 (nchange != 0)
{
printf("percolate: WARNING max steps = %d reached but nchange != 0\n",
maxstep);
}
/*
* Copy the centre of old, excluding the halos, into map
*/
for (i=1; i<=M; i++)
{
for (j=1; j<=N; j++)
{
map[i-1][j-1] = old[i][j];
}
}
/*
* Test to see if percolation occurred by looking for positive numbers
* that appear on both the top and bottom edges
*/
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);
return 0;
}
/******************************************************************************
* Alloc Interface functions to dynamic store allocators. *
* arralloc() Allocate rectangular dope-vector (ie using pointers) array *
******************************************************************************/
/*========================== Library include files ===========================*/
#include
#include
#include
/*========================== Library declarations ============================*/
/* char *calloc(); */
/* char *malloc(); */
/*========================== External function declarations ==================*/
#ifdef DEBUG
int malloc_verify();
int malloc_debug();
#endif
/******************************************************************************
* ~arralloc. Allocate a psuedo array of any dimensionality and type with *
* specified size for each dimension. Each dimension is *
* an array of pointers, and the actual data is laid out in standard ‘c’ *
* fashion ie last index varies most rapidly. All storage is got in one *
* block, so to free whole array, just free the pointer array. *
* array = (double***) arralloc(sizeof(double), 3, 10, 12, 5); *
******************************************************************************/
/* ALIGN returns the next b byte aligned address after a */
#define ALIGN(a,b) (int*)( (((long)(a) + (b) – 1)/(b))*(b) )
/* If on an I860 align arrays to cache line boundaries */
#ifdef I860
#define MIN_ALIGN 32
#else
#define MIN_ALIGN 1
#endif
/*———————————————————————-*/
void subarray(align_size, size, ndim, prdim, pp, qq, dimp, index)
size_t align_size; /* size of object to align the data on */
size_t size; /* actual size of objects in the array */
int ndim, prdim; /* ndim- number of dim left to do */
/* prdim – no of pointers in previous iteration */
void ***pp, **qq; /* pp – pointer to previous level of the array */
/* qq – pointer to start of this level */
int *dimp, index;
{
int *dd = ALIGN(qq,align_size); /*aligned pointer only used in last recursion*/
int **dpp = (int**)pp;
int i, dim = dimp[index];
if(ndim > 0) /* General case – set up pointers to pointers */
{
for( i = 0; i < prdim; i++)
pp[i] = qq + i*dim; /* previous level points to us */
subarray(align_size, size, ndim-1, prdim*dim,
(void***)qq, /* my level filled in next */
qq+prdim*dim, /* next level starts later */
dimp, (index+1) );
}
else /* Last recursion - set up pointers to data */
for( i = 0; i < prdim; i++)
dpp[i] = dd + (i*dim)*size/sizeof(int);
}
/*-----------------------------------------------------------------------*/
/*
* if REFS is defined the va macros are dummied out. This is because the
* GreenHills va_arg macro will not get past the cref utility.
* This way the call tree can still be constructed. Do NOT under
* any circumstance define REFS when compiling the code.
*/
#if REFS
#undef va_start
#undef va_arg
#undef va_end
#undef va_list
#define va_list int
#define va_start( A , B) ( A = (int) B)
#define va_end( A ) ( A = 0 )
#define va_arg( A , T ) ( A = (T) 0)
#endif
void *arralloc(size_t size, int ndim, ...)
{
va_list ap;
void **p, **start;
int idim;
long n_ptr = 0, n_data = 1;
int *dimp;
size_t align_size;
va_start(ap, ndim);
if( size % sizeof(int) != 0 ) /* Code only works for 'word' objects */
return 0;
/* we want to align on both size and MIN_ALIGN */
if( size > MIN_ALIGN )
{
align_size = size;
}
else
{
align_size = MIN_ALIGN;
}
while( (align_size % size) || (align_size % MIN_ALIGN) )
{
align_size++;
}
/*
* Cycle over dims, accumulate # pointers & data items.
*/
if ( NULL == (dimp=(int *)malloc( ndim * sizeof(int) )))
return 0;
for(idim = 0; idim < ndim; idim++)
{
dimp[idim] = va_arg(ap, int);
n_data *= dimp[idim];
if( idim < ndim-1 )
n_ptr += n_data;
}
va_end(ap);
/*
* Allocate space for pointers and data.
*/
if( (start = (void**)malloc(
(size_t)((n_data*size)+align_size+(n_ptr*sizeof(void**))))) == 0)
return 0;
/*
* Set up pointers to form dope-vector array.
*/
subarray(align_size, size, ndim-1, 1, &p, start, dimp, 0);
free( dimp );
return (void*)p;
}
#include
#include
#include “percolate.h”
#include “arralloc.h”
/*
* Simple serial program to test for percolation of a cluster.
*
* This version illustrates the use of dynamic array allocation via the
* arralloc routine.
*
* In a real program, the array dimensions would be variables
* (e.g. taken from the command line arguments). For simplicty, Here
* we just use constant values as in the static allocation version.
*/
int main(int argc, char *argv[])
{
/*
* Define the main arrays for the simulation
*/
// int old[M+2][N+2], new[M+2][N+2];
int **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
*/
// int map[L][L];
int **map;
/*
* Variables that define the simulation
*/
int seed;
double rho;
/*
* Local variables
*/
int i, j, nhole, step, maxstep, oldval, newval, nchange, printfreq;
int itop, ibot, perc;
double r;
if (argc != 2)
{
printf(“Usage: percolate
return 0;
}
/*
* Update for a fixed number of steps, periodically report progress
*/
maxstep = 5*L;
printfreq = 100;
/*
* 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);
/*
* Allocate 2D LxL integer arrays dynamically
*/
old = (int **) arralloc(sizeof(int), 2, M+2, N+2);
new = (int **) arralloc(sizeof(int), 2, M+2, N+2);
map = (int **) arralloc(sizeof(int), 2, L, L);
if (NULL == old || NULL == new || NULL == map)
{
printf(“percolate: array allocation failed\n”);
return 1;
}
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) );
/*
* Initialise the old array: copy the LxL array map 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] = map[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)
{
nchange = 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)
{
++nchange;
}
}
new[i][j] = newval;
}
}
/*
* Report progress every now and then
*/
if (step % printfreq == 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 (nchange != 0)
{
printf("percolate: WARNING max steps = %d reached but nchange != 0\n",
maxstep);
}
/*
* Copy the centre of old, excluding the halos, into map
*/
for (i=1; i<=M; i++)
{
for (j=1; j<=N; j++)
{
map[i-1][j-1] = old[i][j];
}
}
/*
* Test to see if percolation occurred by looking for positive numbers
* that appear on both the top and bottom edges
*/
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.
*/
mapwritedynamic("map.pgm", map, L, 2);
/*
* Free the arrays
*/
free(old);
free(new);
free(map);
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 exercise.
*/
/*
* System size L
*/
#define L 480
/*
* Although overall system is square, i.e. size L x L, we will define
* different variables for the first and second dimensions. This is
* because, in the parallel code, the local arrays will not be
* square. For example, using a simple 1D decomposition over NPROC
* processes, then M = L/NPROC and N = L
*/
#define M L
#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= icc
CFLAGS= -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);
}
void *arralloc(size_t size, int ndim, ...);
!
! Main module file for percolation code.
!
module percmod
implicit none
!
! System size L
!
integer, parameter :: L = 480
!
! Although overall system is square, i.e. size L x L, we will
! define different variables for the first and second
! dimensions. This is because, in the parallel code, the local
! arrays will not be square. For example, using a simple 1D
! decomposition over nproc processes, then M = L and N = L/nproc
!
integer, parameter :: M = L
integer, parameter :: N = L
end module percmod
!
! Simple serial program to test for percolation of a cluster
!
program percolate
!
! Main simulation parameters
!
use percmod
!
! Modules containing supplied functions
!
use percio
use unirand
!
! 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
!
! Variables that define the simulation
!
double precision :: rho
integer :: seed
!
! Local variables
!
integer :: argc, i, j, nhole, step, maxstep
integer :: nchange, printfreq
integer :: itop, ibot
logical :: perc
character*(16) :: tmparg
double precision :: r
argc = command_argument_count()
if (argc /= 1) then
write(*,*) "Usage: percolate
stop
end if
!
! Update for a fixed number of iterations
!
maxstep = 5*L
printfreq = 100
!
! 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) ! ! Initialise the old array: copy the LxL array map to the centre of ! old, and set the halo values to zero. ! old(1:M,1:N) = map(:,:) ! 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) ! ! 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) ) nchange = count(new(1:M,1:N) /= old(1:M,1:N)) ! ! Report progress every now and then ! if (mod(step, printfreq) == 0 ) then write(*,*) "percolate: changes on step ", step, " is ", nchange 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 (nchange /= 0) then write(*,*) "percolate: WARNING max steps = ", maxstep, & " reached but nchange /= 0" end if ! ! Copy the centre of old, excluding the halos, into map ! map(:,:) = old(1:L, 1:L) ! ! Test to see if percolation occurred by looking for positive numbers ! that appear on both the top and bottom edges ! 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
map(:,:) = old(1:L, 1:L)
!
! 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 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= ifort
FFLAGS= -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