
\documentstyle[11pt]{article}      

% Stuff for incorporating postscript figures in Latex
\input epsf

\begin{titlepage}

\baselineskip=24pt

\begin{center}
{\LARGE \bf
Parallel Wolff Cluster Algorithms \\}
\bigskip
\bigskip
{\large 
SUNG-HOON KO  SUNG-YONG PARK\\}
\smallskip
{\it
School of Computer Science, Syracuse University,\\
111 College Place, Syracuse, NY 13244, U.S.A.\\}
%
\bigskip


\bigskip
\bigskip
November 17, 1994 \\
\end{center}
\bigskip
\bigskip

\begin{abstract}
\noindent
The Wolff single-cluster algorithm is the most efficient method known
for Monte Carlo simulation of many spin models. 
Due to the irregular size, shape and position of the Wolff clusters,
this method does not easily lend itself to efficient parallel
implementation, so that simulations using this method have thus far
been confined to workstations and vector machines.
Here we present parallel implementations of this algorithm, 
\end{abstract}

\end{titlepage}

\textwidth=6.0in
\textheight=8.5in
\topmargin=-0.5in
\oddsidemargin=0.25in

\begin{document}
\baselineskip=24pt

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Monte Carlo simulation is widely used in studying the statistical
mechanics of spin models of magnetism~\cite{MONTECARLO}.
Unfortunately, traditional Monte Carlo algorithms for these models,
such as the commonly used Metropolis algorithm
\cite{MONTECARLO,METROPOLIS}, suffer from critical slowing down
\cite{MONTECARLO,SOKALREV}, which means the computational
efficiency of these methods decreases rapidly as the size of
the system is increased.
Algorithms have been developed by Swendsen and Wang~\cite{SW}
and Wolff~\cite{WOLFF} that make large-scale, 
non-local changes in the spin configurations, and greatly reduce
critical slowing down for certain types of spin models.
In these so-called cluster algorithms, clusters of spins
(rather than single spins) are collectively updated at each step
of the Monte Carlo procedure. The clusters are formed by generating
bonds connecting neighboring sites, using a probabilistic
procedure that varies between different spin models and algorithms
(for reviews of cluster algorithms, see
Refs.~\cite{SOKALREV,LAT89REV,SWREV,CLIVEREV}).

Here we will consider the simplest of spin models, the two 
dimensional Ising model~\cite{MONTECARLO}, for which the spin at
any site of the lattice can take only two values, $+1$ or $-1$.
The parallel algorithms we will present are, however, easily 
generalizable to any dimension and any spin model for which a 
Wolff cluster algorithm can be defined. 

In the Swendsen-Wang algorithm for the Ising model, 
bonds are placed with probability ${1-e^{-2\beta}}$ between 
all neighboring sites of the lattice that have the same spin value,
where $\beta$ is the inverse temperature parameter.
This procedure creates clusters of bonded sites of the same spin.
The spins are updated by assigning a random spin value to each
cluster.

In the Wolff algorithm, only a single cluster is created.
A site is chosen at random and a cluster is constructed around it,
using the same bond probabilities as for the Swendsen-Wang algorithm.
All the spins in this cluster are then flipped, that is, collectively
changed to the alternate spin value. 

The performance of a Monte Carlo algorithm is best measured in terms
of statistically independent configurations produced per unit 
(computer) time. Using this measure, cluster algorithms can easily 
outperform the Metropolis algorithm on sequential computers,
with the Wolff algorithm generally better than Swendsen-Wang
\cite{WOLFF,SWREV,ITO,WOLFFCRAY,KLEIN,POTTS,ISING}. 
However this may not be the case on supercomputers, which have 
vector and/or parallel architectures. Very efficient vector and
parallel programs can be written for the Metropolis algorithm, 
since it is regular and local. Cluster algorithms, which are highly
irregular and non-local, are much more difficult to implement
efficiently on vector and parallel machines.
Unless efficient vector or parallel cluster algorithms can be
developed, cluster algorithms may give inferior performance
to the Metropolis algorithm on supercomputers~\cite{ITO}.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Computational algorithms}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The major computational task for the Swendsen-Wang algorithm is the
identification and labeling of the clusters of connected sites, 
given the configuration of bonds. 
This is an instance of a connected component labeling problem
for an undirected graph~\cite{COMPSCI}, where the vertices are the 
lattice sites and the edges are the bonds between connected sites.
The goal of the component labeling algorithm is to end up with the 
same label on all connected sites, and different labels for all 
disconnected clusters.

A number of methods have been developed for implementing the 
Swendsen-Wang algorithm on parallel
\cite{BURKITT,DMCC5,CPE,BROWER,TAMAYO,SIMDLETTER,SIMDPAPER,ROSSI}
and vector~\cite{MINO} computers.
Most of these methods are variations on parallel algorithms that
were developed for the general problem of connected component 
labeling of a graph~\cite{SHIVISH,EMBRECHTS,WOO}.  
Parallel Swendsen-Wang algorithms can be quite efficient,
especially for large lattice sizes on coarse-grain MIMD machines.
The Swendsen-Wang algorithm seems inherently better suited to a
parallel implementation than the Wolff algorithm, since clusters
are constructed for all the sites in the lattice, so a data 
parallel algorithm using standard domain 
decomposition~\cite{THEBOOK} can be used. 
This is not the case for the Wolff algorithm,
for which only a (random) subset of the sites are updated at 
every iteration, so a standard domain decomposition would give
extremely poor load balance between processors.

Constructing (or labeling) a single cluster of connected sites is 
done sequentially by a simple breadth-first-search 
algorithm~\cite{WOLFF,COMPSCI,CPE,ANTS,EVERTZ}.
Two data structures are used:
a list {\tt C} of sites that belong to the 
cluster;\footnote{It is actually more efficient to store this
information as a logical array with a TRUE or FALSE value
for each site in the lattice, to save searching the list, but
here we will consider it as a list of sites in the cluster.}
and a queue {\tt Q} of sites that have been added to the cluster
but not yet used in the search.
The algorithm is described by the pseudo-code in 
Fig.~\ref{fig:SeqAnts}.
The sites in the queue can be visualized as a ``wavefront''
expanding outwards from the initial site. 
It is possible to rewrite the loop over sites in the queue as
a double loop: one loop over all sites in a particular wavefront,
and another loop over successive wavefronts. 
% CHECK THAT SEQ CODE DOES THIS, AND GENS WE HAVE ARE RIGHT ??
We will refer to these successive wavefronts as {\it generations} of
sites, with the first site being the parent (the first generation), 
the set of sites connected by a bond to this site being its
children (the second generation), 
the set of sites connected to these children (and not already
part of the cluster) being the next generation, 
and so on~\cite{CPE,ANTS}.
In this algorithm, the loop over all sites in a particular 
generation can be done in parallel~\cite{CPE,ANTS,EVERTZ}. 
We implement this parallel breadth-first-search by splitting the
queue {\tt Q} into two data structures, the current generation
{\tt G} and the next generation {\tt G$^\prime$}.
The algorithm is described by the pseudo-code in 
Fig.~\ref{fig:ParAnts}.


For the Wolff algorithm only nearest-neighbor information is 
required, so the communication overhead for the parallel 
algorithm is not great. However this method performs poorly on the 
type of fine-grain, massively parallel SIMD machine for which it 
was originally proposed~\cite{ANTS}. 
This is because at any point in the algorithm, 
the current generation of sites will occupy at most O($L$)
sites of a lattice of size $L \! \times \! L$, 
and much less than that in the early
generations when the cluster is small (see Fig.~\ref{fig:domain}). 
If there is a single processor per lattice site, 
the load balance will be extremely poor.
The situation is little better for a coarse-grain parallel machine
with many sites per processor, as can be seen in 
Fig.~\ref{fig:domain}. This parallel breadth-first-search
algorithm is well suited to a vector machine, for which the 
degree of parallelism (the vector length of the machine) is much
smaller~\cite{EVERTZ}.
However on a Cray YMP this method achieves speed-ups of less than 10
over non-vectorized algorithms that realize only the scalar
performance of the machine~\cite{WOLFFCRAY,EVERTZ}.

Methods for implementing the Wolff algorithm in parallel on MIMD 
machines have so far been confined to the trivial parallelism of 
running independent simulations, with different parameters or
random number streams, on different processors~\cite{DMCC5,CPE}.
In this case the size of the system that can be simulated is
limited by the memory and speed of a single node of the parallel
machine. In the next two sections we introduce two methods for
parallelizing a single simulation over many processors using 
domain decomposition. The second method performs quite well
on coarse-grain MIMD parallel computers.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The single-cluster Swendsen-Wang algorithm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The Swendsen-Wang algorithm requires the identification of clusters
of connected sites.
This is done using a connected component labeling algorithm to give
each site a number that labels the cluster to which it belongs.
All of the parallel connected component labeling algorithms use
an iterative procedure for updating the cluster label of each site.
At the end of each iteration, a test is performed to see if the 
label of any site has been changed. 
If no change has been made, then all sites have received the label
of their cluster, and the procedure terminates. 

A parallel Wolff algorithm can be implemented using a simple 
change to a parallel Swendsen-Wang algorithm, which creates
a single-cluster version of the algorithm.
As in the sequential Wolff algorithm, a random site is chosen, 
around which the cluster will grow.
The iterative procedure mentioned above remains the same, 
however now the termination condition is not that the labels of 
{\it all sites} are unchanged,
but rather the labels of {\it all sites in the single cluster 
containing the initial random site} are unchanged.
This allows all sites of the Wolff cluster to be worked on in
parallel, rather than just the expanding wavefront. 
However in this method the computation is being done for all clusters
at the same time, not just the single Wolff cluster,
so there is a lot of wasted effort. 

This method was implemented using a SIMD Swendsen-Wang 
program~\cite{SIMDPAPER}
written in CMFortran~\cite{CMF} for the Thinking Machines CM-5 
computer.
We had planned to also implement this method in a message-passing 
language using a MIMD programming model, which can better
handle this type of irregular problem and is known to provide
much better performance than SIMD algorithms~\cite{CPE,SIMDPAPER}.
However we instead developed and implemented an improved MIMD 
algorithm.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The scattered strip partitioning algorithm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Parallel algorithms with local data dependencies generally use a 
regular domain decomposition in order to minimize communications.
However for some dynamic irregular problems, the computational
load tends to be concentrated in different regions of the domain
at different times. In order to balance the load, some kind of 
dynamic data partitioning is often required. However in some cases
a much simpler static approach can be adopted, in which the data
in neighboring regions is scattered among the processors. 
This scattered domain decomposition has a higher communications
cost since data locality is sacrificed,
however this may be more than compensated for by
the improvement in load balance when the computational load is
concentrated in a certain region of the domain, since now all 
processors contain part of that domain. This method has been very
effective in many parallel algorithms, including certain types of
matrix solvers, irregular finite element problems, and problems
with dynamic data distributions~\cite{THEBOOK}.

Scattered domain decomposition has also been effective in a parallel
implementation of Lee's maze routing algorithm~\cite{FANG}.
Maze routing is an algorithm used for routing wires in VLSI 
circuits~\cite{LEE}. The model is basically that of a regular grid
where some paths are blocked and some are allowed. The problem is
to find the shortest allowed path to place a wire between two 
points. The maze routing algorithm performs a breadth-first-search
of the space, starting from the initial point, and following all
allowed paths until one path reaches the final point. It is 
therefore very similar to the growth of a cluster in the Wolff
cluster algorithm. 

Fang et al.~\cite{FANG} found that the most effective domain 
decomposition for a parallel maze routing algorithm was scattered
strip partitioning (this is often referred to as a cyclic or
column-cyclic data distribution, for example in the 
High Performance Fortran language specification~\cite{HPF}). 
In this scheme, the lattice is split into $N$ partitions, where
$N$ is chosen to be an integer multiple of the number of processors
$P$ to ensure load balancing. 
Each partition contains $W$ columns of the lattice 
(we shall refer to $W$ as the {\it partition width}). 
Thus for a lattice of size $L \! \times \! L$, $L = W\!N$ and 
$W$ must be in the range $1 \leq W \leq L/P$.
Partition $M$ is assigned to processor $M \bmod P$
(see Fig.~\ref{fig:cyclic}).
The mapping of data to processors is:
%
\begin{itemize}
\item  Global lattice coordinate $(i,j)$ is assigned to processor
       $(j \bmod W\!P)$ div $W$
\item  Local lattice coordinate $(a,b)$ is given by  
\begin{list}{ }
\item  $a = i$ 
\item  $b = (j \: \bmod \: W\!P) \, \bmod \, W  \: + \:
	    W \times \, ( j\,$  div $W\!P)$ 
\end{list}
\end{itemize}
%
An example of the distribution of data after a scattered strip 
partitioning is shown in Fig.~\ref{fig:partitioning}. 
The load balance is still far from perfect, but it is much 
better than for the standard domain decomposition in 
Fig.~\ref{fig:domain}.
As $W$ becomes smaller, we can expect better load balancing, but
a larger communications overhead due to non-locality of data.
It is therefore important to find the optimal partition width $W$.

We have implemented this algorithm in a MIMD message-passing
paradigm. The following is an outline of the steps required:
\begin{enumerate}
%
\item {\bf Lattice partitioning} \vspace*{1mm} \\
Do a scattered strip (column-cyclic) partitioning of the lattice,
so that the portion of the lattice allocated to each processor 
has $L / (W\!P)$ partitions of width $W$, i.e. $L/P$ columns and
$L^2 / P$ sites in total.
%
\item {\bf Selection of an initial random site} \vspace*{1mm} \\
One of the processors (processor 0) selects the initial random site 
and broadcasts the information to the other processors. 
%
\item {\bf Local expansion of the cluster}  \vspace*{1mm} \\
Each processor maintains a queue that stores those sites in the
local lattice that are in the current generation, and also
maintains a communication buffer to collect all communication 
requests. \vspace*{2mm} \\
Each local site in the current generation is fetched from the queue,
and used as a parent site to expand the cluster. \vspace*{2mm} \\
Bonds are made with the appropriate probability between the parent 
site and all its neighboring sites, unless the neighboring site
is already an element of the cluster. \vspace*{2mm} \\ 
At this time, if the parent site is on the boundary 
(the left or right edge) of a partition,
communication with the left or right processor may be needed. 
Such communication requests are saved into the communication buffer
for collective communication. \vspace*{2mm} \\
All local connected sites are added to the cluster and the queue for
the next generation of sites, and the spin values at these sites are 
updated. 
%
\item {\bf Collective communication} \vspace*{1mm} \\
Any communication requests are collected from the communication
buffer and sent to the destination processors.  
%
\item {\bf Remote expansion of the cluster} \vspace*{1mm}\\ 
Each site in the received buffer is checked to see whether
it is already an element of the cluster.
If not, the site is added to the local cluster list and the queue
for the next generation, and its spin value is updated. 
%
\item {\bf Check for Termination} \vspace*{1mm} \\
The final step checks whether there are any sites in the new
generation. This can be done easily by using a global reduction 
function such as logical AND. 
If every processor has the empty queue, the algorithm halts.
If not, steps 3, 4 and 5 are repeated.
%
\end{enumerate}


The parallel code to implement this message-passing program
was written in C for the CM-5, using CMMD~\cite{CMMD} 
reduction functions for global communications (steps 2 and 6)
and Active Messages~\cite{CMMD,ACTIVE} for point-to-point 
communications (step 4). 

% OVERLAP COMPUTATION AND COMMUNICATION?? IS THIS DONE?? no

This algorithm should be easily implementable in High Performance 
Fortran (HPF) using a column-cyclic data distribution~\cite{HPF}.
The HPF program would be much simpler than the message-passing
program we have implemented, since there would be no need to 
set up communications buffers and the reading and writing
of messages to access non-local data  
-- this would be handled by the compiler.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Results}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The two algorithms have been run on a 32-node Thinking Machines
CM-5 using lattices of size $128^2$ to $2048^2$. 
The results were compared with a sequential program written in C 
and run on a single node of the CM-5.
The inverse temperature $\beta$ was taken to be the critical
inverse temperature 
$\beta_{c} = \log (1 + \sqrt 2) / 2 \approx 0.4406868$ 
for the 2-$d$ Ising model, 
where the clusters are most irregular in size and shape. 

To measure the average execution time, we performed 100 iterations
to thermalize the system, and then measured the execution time 
for at least 100 (and usually 1000) iterations. 
The time for each iteration depends on the size of the cluster 
that is grown, which can be any size from a single site to many
thousands of sites. In order to make a sensible comparison to other
Monte Carlo algorithms such as Metropolis and Swendsen-Wang, 
for which every site is updated at each iteration, 
the execution time for the Wolff algorithm is given as the time
per spin update, rather than the time per iteration.

We measured the performance using both the scalar processor and
the vector processors on the CM-5. 
As expected, using the vector units greatly improves the performance 
of the SIMD algorithms, so the results shown for these are
using the vector units. 
However, for the sequential and MIMD algorithms, the vector units are not
effective since the algorithm is so irregular,
and in fact the program runs slower if they are used, 
presumably due to the overhead of loading the data onto the vector
registers. For these algorithms the times quoted are without vector units.

For the scattered strip partitioning algorithm, we tested various
partition widths $W$ to find the optimal width for each lattice size.
Fig.~\ref{fig:width} shows the scaled time for different values of
$W$. The performance is very sensitive to the partition width. 
Larger widths mean smaller communications overhead, but greater
load imbalance. Since load imbalance is less of a problem
for larger generation sizes (and hence larger lattice sizes), 
the optimal width starts at $W=1$ for small lattices and slowly 
increases with lattice size.

The timings and efficiencies (taken relative to the sequential Wolff
algorithm run on a single processor) for the single-cluster Swendsen-Wang
algorithm, the full Swendsen-Wang algorithm from which it is derived,
and the scattered-strip partitioning algorithm with $W=2$,
are shown in Table~\ref{tab:efficiencies} for a 32 node CM-5.
Reasonable speed-ups of around 10 are obtained for the MIMD algorithm
for the larger lattice sizes. The corresponding efficiencies of around
35\% are not particularly good, but perhaps better than expected for
this highly irregular and dynamic problem.
Even using the vector units, the efficiencies for the SIMD algorithms
are very low. This is partly because this irregular problem is not 
well-suited to SIMD methods, but partly because the CMFortran compiler
is not very
efficient at implementing this type of algorithm, for which the ratio of
computation to communication is low. The same SIMD algorithms programmed
in a message-passing paradigm using CMMD or Active Messages
would be much more efficient.

The amount of parallelism that can be extracted from the
scattered-strip partitioning algorithm 
depends on the average number of sites in each generation.
Fig.~\ref{fig:gensize} shows the average generation (or wavefront) 
size $G$ as a function of the linear lattice size $L$ for the 2-d 
Ising model. 
We would expect $G$ to scale slower than $\sqrt C$, where $C$
is the average cluster size. Since $C$ for the Wolff algorithm is an
estimator of the susceptibility~\cite{KLEIN,CSIZE}, we have 
${C \sim L^{\gamma / \nu}}$, where ${\gamma / \nu = 1.75}$ 
is the finite-size scaling result for the susceptibility 
of the 2-d Ising model. 
Thus $G$ should scale slower than $L^{0.875}$, and we measured the
actual exponent to be 0.66(1).
% IS THERE ANY THEORETICAL (PERCOLATION??) REASON FOR IT TO BE 2/3??

For this algorithm, the number of processors that can be
effectively used is approximately $G/2$, so on a parallel machine
with 32 nodes (as we have used) reasonable efficiencies should be
obtained for lattices of size greater that $512^2$, which is
indeed what we found. In fact these are the kinds of lattice sizes
one would like to be able to run in parallel -- anything smaller
than this could more easily be run on a single processor, with 
multiple simulations on the different processors of a parallel 
machine~\cite{DMCC5,CPE}.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

We have introduced two new parallel implementations of the Wolff
single-cluster Monte Carlo algorithm.
%
The first method involves a simple change to a parallel
Swendsen-Wang cluster algorithm, so that the iterative procedure
for identifying the clusters is halted as soon as the single Wolff
cluster is identified.
This method performs about as well as the parallel Swendsen-Wang
algorithms on which it is based,
however neither of these parallel cluster algorithms are well-suited
to implementation on massively parallel SIMD machines.

The scattered strip partitioning algorithm successfully alleviates
the load balance problem that is the main difficulty in constructing
a parallel Wolff cluster algorithm. This is done by partitioning the
lattice using a scattered type of decomposition, so that the 
expanding wavefront is distributed over almost all processors.
This method of distributing the data produces reasonable speed-ups 
for the parallel breadth-first-search algorithm for growing
a single cluster, as long as the linear dimension of the lattice
is much greater than the number of processors, and the number of
processors is not too great. 

On a massively parallel machine, parallelism could be extracted
by a combination of the strip partitioning algorithm and the
independent (or job-level) parallelism of running independent 
Monte Carlo simulations with different random number streams on
different groups of processors. For example, one might run 16 
independent simulations, each of which use 16 processors.
In this scenario, the main advantage of using the parallel 
algorithm is that it avoids the memory limitations of a single 
processor, and allows the use of larger lattice sizes.

Note that the results for the speed-ups and efficiencies of our
parallel algorithms are for the simplest of spin models, 
the Ising model. For more complex models, for example continuous 
spin models such as the O(N) model, the computation required to 
calculate the bond probabilities is much greater, 
and we would therefore expect greater efficiencies due to the
higher ratio of computation to communication.
The speed-up also depends heavily on the average generation size,
which will vary between different spin models and lattice sizes.

It is possible that this technique could also be used to improve
the performance of parallel Swendsen-Wang algorithms. These 
algorithms also suffer from load imbalance, since at the critical
point there is generally one large cluster that takes the longest
time to update. With a standard domain decomposition, there are
generally processors that contain little or no sites belonging to
the largest cluster. These processors will be mainly idle during the
latter stages of the cluster labeling. If a scattered strip
partitioning of the lattice were used, this will no longer be the 
case, and the load will be much more evenly divided. 
However for the MIMD algorithm the strips will have to be much wider,
in order for the sequential labeling algorithm on each processor 
to work effectively, 
so there are tradeoffs that need to be investigated to find the
most effective implementation for a given lattice size and number 
of processors.

Simulation of the 2-d Ising spin model using the Wolff algorithm has been
shown to be a very sensitive test of the randomness properties of random
number generators~\cite{FERRENBERG,CODDINGTON}.
One of the motivations for this work was to create a parallel Wolff 
algorithm that can be used as a test for parallel random number generators.
This work is currently in progress.

% MULTIPLE CLUSTER METHOD AS WORK IN PROGRESS??


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section*{Acknowledgements}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This work was done using the Thinking Machines CM-5 at the
Northeast Parallel Architectures Center (NPAC) at Syracuse University.
Work supported in part by the Center for Research on Parallel
Computation with NSF cooperative agreement No. CCR-9120008
and ARPA under contract No. DABT63-91-K-0005.



\begin{thebibliography}{00}
%
\bibitem{MONTECARLO}
K.~Binder ed., 
 Monte Carlo Methods in Statistical Physics, 
(Springer-Verlag, Berlin, 1986);
K.~Binder and D.W.~Heermann,  Monte Carlo Simulation in 
Statistical Physics, (Springer-Verlag, Berlin, 1988);
H.~Gould and J.~Tobochnik,
 An Introduction to Computer Simulation Methods, Vol. 2,
(Addison-Wesley, Reading, Mass., 1988).
%
\bibitem{METROPOLIS}
N.~Metropolis  et al., 
 J. Chem. Phys. 21 (1953) 1087.
%
\bibitem{SOKALREV}
A.~D.~Sokal, in Computer Simulation Studies in Condensed Matter
Physics:  Recent Developments, eds.\ D.~P.~Landau et al.
(Springer-Verlag, Berlin, 1988);
A.~D.~Sokal, % How to Beat Critical Slowing Down: 1990 Update,
in Proc. of the International Conference on Lattice Field Theory,
Tallahassee, October 1990,  Nucl. Phys. B (Proc. Suppl.) 
20 (1991) 55.
%
\bibitem{SW}
R.H.~Swendsen and J.-S.~Wang, 
% Nonuniversal critical dynamics in Monte Carlo simulations,
 Phys. Rev. Lett. 58 (1987) 86.
%
\bibitem{WOLFF}
U. Wolff, 
 Phys. Rev. Lett. 62 (1989) 361.
%
\bibitem{LAT89REV}
U.~Wolff,
in Proc. of the Symposium on Lattice Field Theory, Capri, September
1989,  Nucl. Phys. B (Proc. Suppl.) 17 (1990) 93. 
%
\bibitem{SWREV}
J.-S.~Wang and R.~H.~Swendsen,  Physica A 167 (1990) 565.
%
\bibitem{CLIVEREV}
C. F. Baillie,
 Int. J. Mod. Phys. C 1 (1990) 91.
%
\bibitem{ITO}
N.~Ito and G.~A.~Koring,
Cluster vs single-spin algorithms -- which are more efficient?,
to be published in Int. J. Mod. Phys. C.
% ????
%
\bibitem{COMPSCI}
E.~Horowitz and S.~Sahni, Fundamentals of Computer Algorithms,
(Computer Science Press, Rockville, Maryland, 1978);
G.~Brassard and P.~Bratley, Algorithmics: Theory and Practice,
(Prentice Hall, Englewood Cliffs, N.J., 1988). 
% MORE RECENT REFS CONTAINING BOTH CCL and BFS??
%
\bibitem{WOLFFCRAY}
U. Wolff,  Phys. Lett. B228 (1989) 379. 
%
\bibitem{KLEIN}
P.~Tamayo, R.~C.~Brower and W.~Klein,
 J. Stat. Phys. 58 (1990) 1083.
%
\bibitem{POTTS}
C.~F.~Baillie and P.~D.~Coddington,
 Phys. Rev. B 43 (1991) 10617.
%
\bibitem{ISING}
P.~D.~Coddington and C.~F.~Baillie,
 Phys. Rev. Lett. 68 (1992) 962.
%
\bibitem{BURKITT}
A.~N.~Burkitt and D.~W.~Heermann,  Comp. Phys. Comm. 54 
(1989) 210.
%
\bibitem{DMCC5}
P.D.~Coddington and C.F.~Baillie, in
Proc. of the 5th Annual Distributed Memory Computing Conference,
Charleston, April 1990, eds. D.W. Walker and Q.F. Stout
(IEEE Computer Society Press, Los Alamitos, California, 1990).
%
\bibitem{CPE}
C.~F.~Baillie and P.~D.~Coddington,
 Concurrency: Practice and Experience 3 (1991) 129.
%
\bibitem{BROWER}
R.~C.~Brower, P.~Tamayo and B.~York,  J. Stat. Phys. 63
(1991) 73.
%
\bibitem{TAMAYO}
M.~Flanigan and P.~Tamayo,  Int. J. Mod. Phys. C 3 (1992)
1235.
%
\bibitem{SIMDLETTER}
J. Apostolakis, P. Coddington and E. Marinari,
 Europhys. Lett. 17 (1992) 189.
%
\bibitem{SIMDPAPER}
J. Apostolakis, P. Coddington and E. Marinari,
 Int. J. Mod. Phys. C 4 (1993) 749.
%
\bibitem{ROSSI}
P.~Rossi and G.~P.~Tecchiolli,  Finding Clusters in a Parallel
Environment, unpublished.
% (REF???)
%
\bibitem{MINO}
H. Mino,
 Comp. Phys. Comm. 66 (1991) 25.
%
\bibitem{SHIVISH}
Y.~Shiloach and U.~Vishkin,  J. Algorithms 3 (1982) 57. 
%
\bibitem{EMBRECHTS}
H.~Embrechts, D.~Roose, and P.~Wambacq, 
in  Proc. First European Workshop on
Hypercube and Distributed Computers, F. Andre and J.P. Verjus eds., 
(North-Holland, Amsterdam, 1989); 
H.~Embrechts, D.~Roose, and P.~Wambacq, 
 Computer Vision Graphics and Image Processing:
Image Understanding, 57 (1993) 155.
%
\bibitem{WOO}
J. Woo and S. Sahni,  J. of Supercomputing 3 (1989) 209.
%
\bibitem{ANTS}
R.~Dewar and C.~K.~Harris,  J. Phys. A 20 (1987) 985.
%
\bibitem{EVERTZ}
H.~G.~Evertz,
 J. Stat. Phys. 70 (1993) 1075.
%
\bibitem{THEBOOK}
G.~C.~Fox  et al.,
 Solving Problems on Concurrent Processors,
(Prentice-Hall, Englewood Cliffs, New Jersey, 1988).
%
\bibitem{CMF}
CMFortran Reference Manual, (Thinking Machines Corporation,
Cambridge, Mass., 1993).
%
\bibitem{FANG}
Y.-Y.~Fang, I.-L.~Yen and R.~Dubash,
Improving the performance of Lee's maze routing algorithm
on parallel computers via semi-dynamic mapping strategies,
Technical Report CPS-93-35, Michigan State University, December, 1993.
%
\bibitem{LEE}
C.~Y.~Lee, IRE Trans. Electronic Computers EC-10 (1961) 346.
%
\bibitem{HPF}
The High Performance Fortran Forum,
High Performance Fortran Language Specification,
Center for Research in Parallel Computing Technical Report 
CRPC-TR92225.
Available via anonymous ftp from
titan.cs.rice.edu in the directory public/HPFF/draft.
%
\bibitem{CSIZE}
U.~Wolff,
Nucl. Phys. B 334 (1990) 581.
%
\bibitem{CMMD}
CMMD Reference Manual, (Thinking Machines Corporation,
Cambridge, Mass., 1993).
%
\bibitem{ACTIVE}
T.~von Eicken et al.,
Active Messages: a Mechanism for Integrated Communication and Computation,
in Proc. of the 19th Int. Symposium on Computer Architecture,
Gold Coast, Australia, May, 1992.
% Report No. UCD/CSD 92/675, March 1992.
%
\bibitem{FERRENBERG}
A.M.~Ferrenberg, D.P.~Landau and Y.J.~Wong,
% Monte Carlo simulations: Hidden errors from ``good'' random number generators,
Phys. Rev. Lett. 69 (1992) 3382.
%
\bibitem{CODDINGTON}
P.D. Coddington,
Analysis of Random Number Generators Using Monte Carlo Simulation,
NPAC technical report SCCS-526, to be published in Int. J. Mod. Phys. C.
%
% \bibitem{STAUFFER}
% D.~Stauffer,  Introduction to Percolation Theory,
% (Taylor and Francis, Philadelphia, 1985).
%
\end{thebibliography}


\newpage


\begin{table}[htbp]
\vspace*{2.0cm}
\renewcommand{\arraystretch}{1.3}
\begin{center}
\begin{tabular} {|c|ccccc|} 
\hline
Lattice Size   & $128^2$ & $256^2$ & $512^2$ & $1024^2$ & $2048^2$ \\
\hline
Time ($\mu s$) Sequential & 31.42 & 32.36 & 32.81 & 33.74 & --
\\ \hline
Time ($\mu s$) MIMD   & 
\phantom{0}6.50  & \phantom{0}5.10  & \phantom{0}3.96  & 
\phantom{0}3.26  & \phantom{0}2.93 \\
Speed-up MIMD         & 
\phantom{0}4.8\phantom{0}  & \phantom{0}6.3\phantom{0} & 
\phantom{0}8.3\phantom{0}  & 10.4\phantom{0}  &  11.5\phantom{0} \\
Efficiency MIMD       &
\phantom{0}0.15  & \phantom{0}0.20  & \phantom{0}0.26  &
\phantom{0}0.32  & \phantom{0}0.36 \\
Time ($\mu s$) Wolff SIMD   & 
40.70 & 44.90 & 46.56 & 63.09 & -- \\
Time ($\mu s$) S-W SIMD   & 
20.22 & 19.23 & 17.38 & 18.34 & -- \\
\hline
\end{tabular}
\end{center}
\vspace*{2.0cm}
\protect\caption{\protect\label{tab:efficiencies}
Times per spin update, speed-ups, and 
efficiencies as a function of lattice size for 
the MIMD scattered strip partitioning algorithm,
the SIMD Wolff algorithm (the single-cluster Swendsen-Wang algorithm),
and the SIMD Swendsen-Wang (S-W) algorithm from which it is derived,
all for a 32-node CM-5.
Speed-ups and efficiencies for $2048^2$ are based on sequential
$1024^2$ times, so are probably underestimated.
}
\end{table}


% \newpage

\begin{figure}[tb]
\vspace*{2cm}
\begin{tabbing}
xx\=xx\=xx\=xx\=XXXXXXXXXXXXXXXXXXXXXXXXXXX\= \kill
 \>  Chose a site {\tt i} at random 
\>\>\>\>  {\it First site in the cluster} \\
 \>  {\tt C $\leftarrow$ \{i\} }
\>\>\>\>  {\it Initialize the cluster list} \\
 \>  {\tt Q $\leftarrow$ \{i\} }
\>\>\>\>  {\it Initialize the queue} \\
 \>  {\tt WHILE Q $\neq$ \{\} DO }
\>\>\>\>  {\it Search over sites still in queue} \\
 \>  \>  {\tt j $\leftarrow\:$} next element in {\tt Q}
\>\>\>    {\it Get the next site from the queue} \\
 \>  \>  {\tt Q $\leftarrow$ Q - j}
\>\>\>    {\it Remove this site from the queue} \\
 \>  \>  Find the set {\tt\{n\}} of connected neighbors of {\tt j} 
\>\>\>    {\it Identify bonds to neighbors} \\
 \>  \>  {\tt FOR} each {\tt n $\!\epsilon\!$ \{n\} DO}
\>\>\>        {\it For all connected neighbors} \\
 \>  \>  \>  {\tt IF n $\!\not\!\epsilon\!$ C THEN} 
\>\>          {\it If neighbor is not in cluster list} \\
 \>  \>  \>  \>  {\tt C $\leftarrow$ C + n}
\>              {\it Add it to the cluster list} \\
 \>  \>  \>  \>  {\tt Q $\leftarrow$ Q + n}
\>              {\it Add it to the queue} \\
 \>  \>  \>  {\tt END IF} \\
 \>  \>  {\tt END DO} \\
 \>  {\tt END WHILE} \\
\end{tabbing}
\vspace*{2cm}
\protect\caption{\protect\label{fig:SeqAnts}
The breadth-first-search algorithm for creating a single cluster.
For the Wolff algorithm, the spin is also changed to the new value 
($s_n \leftarrow - s_n$ for the Ising model) 
when a site $n$ is added to the cluster.
}
\vspace*{2cm}
\end{figure}
 


\begin{figure}[tb]
\vspace*{2cm}
\begin{tabbing}
xx\=xx\=xx\=xx\=xx\=XXXXXXXXXXXXXXXXXXXXXXXXXX\= \kill
 \>  {\tt G$\;$ $\leftarrow$ \{\} }
\>\>\>\>\>  {\it Initialize the current generation} \\
 \>  Chose a site {\tt i} at random 
\>\>\>\>\>  {\it First site in the cluster} \\
 \>  {\tt C$\;$ $\leftarrow$ \{i\} }
\>\>\>\>\>  {\it Initialize the cluster list} \\
 \>  {\tt G$^\prime$ $\leftarrow$ \{i\} }
\>\>\>\>\>  {\it Initialize the next generation} \\
 \>  {\tt WHILE G$^\prime$ $\neq$ \{\} DO }
\>\>\>\>\>  {\it Loop over generations} \\
 \>  \>  {\tt G$\;$ $\leftarrow$ G$^\prime$ }
\>\>\>\>  {\it Go to new generation} \\
 \>  \>   {\tt G$^\prime$ $\leftarrow$ \{\} }
\>\>\>\>  {\it Initialize the next generation} \\
 \>  \>   {\tt WHILE G $\neq$ \{\} DO IN PARALLEL}
\>\>\>\>  {\it Loop over this generation of sites} \\
 \>  \>  \>   {\tt j $\leftarrow\:$} next element in {\tt G}
\>\>\>    {\it Get the next site from the queue} \\
 \>  \>  \>   {\tt G $\leftarrow$ G - j}
\>\>\>    {\it Remove this site from the queue} \\
 \>  \>  \> Find the set {\tt\{n\}} of connected neighbors of {\tt j} 
\>\>\>    {\it Identify bonds to neighbors} \\
 \>  \>  \>   {\tt FOR} each {\tt n $\!\epsilon\!$ \{n\} DO}
\>\>\>        {\it For all connected neighbors} \\
 \>  \>  \>  \>   {\tt IF n $\!\not\!\epsilon\!$ C THEN} 
\>\>          {\it If neighbor is not in cluster list} \\
 \>  \>  \>  \>  \>   {\tt C$\;$ $\leftarrow$ C$\;$ + n}
\>              {\it Add it to the cluster list} \\
 \>  \>  \>  \>  \>   {\tt G$^\prime$ $\leftarrow$ G$^\prime$ + n}
\>              {\it Add it to the next generation} \\
 \>  \>  \>  \>  {\tt END IF} \\
 \>  \>  \>  {\tt END DO} \\
 \>  \>  {\tt END WHILE} \\
 \>  {\tt END WHILE} \\
\end{tabbing}
\vspace*{2cm}
\protect\caption{\protect\label{fig:ParAnts}
The parallel breadth-first-search algorithm for creating a single
cluster. 
}
\vspace*{2cm}
\end{figure}



\begin{figure}
\vspace*{2.0cm}
\epsfxsize=3.0in
\centerline{\epsfbox{fig3.ps}}
\vspace*{4.0cm}
\protect\caption{\protect\label{fig:domain}
An example of a Wolff cluster for a configuration of the 2-$d$ 
Ising spin model.
The light gray sites indicate the cluster, 
the sites in dark gray show a particular wavefront or generation,
and the black site is the initial randomly chosen site.
The black lines represent processor boundaries for a standard
domain decomposition of the lattice onto 16 processors.
}
\end{figure}


\vspace*{2.0cm}

\begin{figure}
\epsfxsize=3.0in
\centerline{\epsffile{fig4.ps}}
\vspace*{4.0cm}
\protect\caption{\protect\label{fig:cyclic}
The mapping of 16 columns of data onto 4 processors for 
scattered strip (or column-cyclic) partitioning with a column
width of 2.
}
\end{figure}


\vspace*{2.0cm}

\begin{figure}
\epsfxsize=3.0in
\centerline{\epsfbox{fig5.ps}}
\vspace*{4.0cm}
\protect\caption{\protect\label{fig:partitioning}
A scattered strip (or column-cyclic) partitioning onto 16
processors for the data in Fig.~3. 
}
\end{figure}


\vspace*{2.0cm}

\begin{figure}
\epsfxsize=4.0in
\centerline{\epsfbox{fig6.ps}}
\vspace*{4.0cm}
\protect\caption{\protect\label{fig:width}
Average update time per spin 
as a function of partition width for different lattice sizes,
for the scattered strip partitioning algorithm on a 32-node CM-5.
}
\end{figure}


\vspace*{2.0cm}

\begin{figure}
\epsfxsize=4.0in
\centerline{\epsfbox{fig7.ps}}
\vspace*{4.0cm}
\protect\caption{\protect\label{fig:gensize}
The average number of sites per generation in the Wolff 
cluster update as a function of lattice size
for the 2-$d$ Ising model. 
The line is a fit to a power with exponent 0.66(1).
}
\end{figure}



\end{document}

