%% TERRI - I have added a couple of things to the Latex file, and
% fixed a couple of grammatical errors. The main thing I have 
% done is added a lot more references. The references are all at
% the end of this file. I have marked the additional citations in
% the text with a comment like
%
% % TERRI - check the numbering on these new refs.
% \cite{Coddington:90?}, \cite{Baillie:91?}.  
%
% so you can insert the appropriate letter for the ?
%

\documentstyle{book}
\tolerance=10000
\def\lea{\mathrel{\raise .4ex\hbox{\rlap{$<$}\lower 1.2ex\hbox{$\sim$}}}}
\def\gea{\mathrel{\raise .4ex\hbox{\rlap{$>$}\lower 1.2ex\hbox{$\sim$}}}}
\newcommand{\medStrut}{{\rule{0cm}{4.0ex}}}
%
\newdimen\digitwidth
\setbox1=\hbox{0}
\digitwidth=\wd1
\catcode`"=\active
\def"{\kern\digitwidth}
%

\begin{document}

\setcounter{chapter}{12}
\setcounter{section}{5}

\section{Cluster Algorithms for Spin Models}
\label{cluster}

\subsection{Monte Carlo Calculations of Spin Models}

The goal of computer simulations of spin models is to generate configurations 
of spins typical of statistical equilibrium and measure physical observables 
on this ensemble of configurations.  The generation of configurations is 
traditionally performed by Monte Carlo methods such as the Metropolis 
algorithm \cite{Metropolis:53a}, which produce configurations with a 
probability given by the Boltzmann distribution $e^{- \beta S(\phi)}$, where 
$S(\phi)$ is the action, or energy, of the system in configuration $\phi$, 
and $\beta$ is the inverse temperature.  One of the main problems with these 
methods in practice is that successive configurations are not statistically 
independent, but rather are correlated with some autocorrelation time,
$\tau$, between effectively independent configurations. 

A key feature about traditional (Metropolis-like) Monte Carlo algorithms is
that the updates are {\em local}, that is, one spin at a time is updated, and 
its new value depends only on the values of spins which affect its 
contribution to the action, {\em i.e.,} only on local (usually nearest 
neighbor) spins.  Thus, in a single step of the algorithm, information about 
the state of a spin is transmitted only to its nearest neighbors.  In order 
for the system to reach a new effectively independent configuration, this 
information must travel a distance of order the (static or spatial) 
correlation length $\xi$.  As the information executes a random walk around 
the lattice, one would suppose that the autocorrelation time, $\tau \sim 
\xi^2$.  However, in general $\tau \sim \xi^z$, where $z$ is called the 
dynamical critical exponent.  Almost all numerical simulations of spin models
have measured $z \approx 2$ for local update algorithms.  (See also
Sections~\ref{qcd}, \ref{spins}, \ref{tc}, \ref{cluster}, and \ref{melting}).

For a spin model with a phase transition, as the inverse temperature $\beta$ 
approaches the critical value, $\xi$ diverges to infinity so that the 
computational efficiency rapidly goes to zero!  This behavior is called 
critical slowing down (CSD), and until very recently it has plagued Monte 
Carlo simulations of statistical mechanical systems, in particular spin 
models, at or near their phase transitions.  Recently, however, several new 
``cluster algorithms'' have been introduced which decrease $z$ dramatically 
by performing {\em non-local\/} spin updates, thus reducing (or even 
eliminating) CSD and facilitating much more efficient computer simulations.

\subsection{Cluster Algorithms}

The aim of the cluster update algorithms is to find a suitable collection of 
spins which can be flipped with relatively little cost in energy.  We could 
obtain non-local updating very simply by using the standard Metropolis Monte 
Carlo algorithm to flip randomly selected bunches of spins, but then the 
acceptance would be tiny.  Therefore, we need a method which picks sensible 
bunches or clusters of spins to be updated.  The first such algorithm was 
proposed by Swendsen and Wang \cite{Swendsen:87a}, and was based on an 
equivalence between a Potts spin model \cite{Potts:52a}, \cite{Wu:82a} and 
percolation models \cite{Stauffer:78a}, \cite{Essam:80a}, for which cluster 
properties play a fundamental role.

The Potts model is a very simple spin model of a ferromagnet, in which the 
spins can take $q$ different values.  The case $q=2$ is just the well-known 
Ising model.  In the Swendsen and Wang algorithm, clusters of spins are 
created by introducing bonds between neighboring spins with probability $
1-e^{-\beta}$ if the two spins are the same, and zero if they are not.
All such clusters are created and then updated by choosing a random new 
spin value for each cluster and assigning it to all the spins in that cluster.

A variant of this algorithm, for which only a single cluster is constructed
and updated at each sweep, has been proposed by Wolff \cite{Wolff:89a}. The 
implementation of this algorithm is shown in Figure~\ref{fig1.cluster} to
Figure~\ref{fig3.cluster}, which shows a $q=3$ Potts model at its critical 
temperature, with different colors representing the three different spin 
values.  From the starting configuration (Figure~\ref{fig1.cluster}), we 
choose a site at random, and construct a cluster around it by bonding 
together neighboring sites with the appropriate probabilities 
(Figure~\ref{fig2.cluster}).  All sites in this cluster are then given the 
same new spin value, producing the new configuration shown in 
Figure~\ref{fig3.cluster}, which is obviously far less correlated with the 
initial configuration than the result of a single Metropolis update 
(Figure~\ref{fig4.cluster}).  Although Wolff's method is probably the best 
{\em sequential} cluster algorithm, the Swendsen and Wang algorithm seems to 
be better suited for parallelization, since it involves the entire lattice 
rather than just a single cluster.  We have, therefore, concentrated our 
attention on parallelizing the method of Swendsen and Wang, where {\em all} 
the clusters must be identified and labeled.

\begin{figure}
\vspace*{2.0in}
\caption{Initial Configuration of 3-state Potts Spins on which Wolff
Algorithm is to be Applied}
\label{fig1.cluster}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{Bonds of Cluster Constructed by Wolff Algorithm Indicated in
Yellow}
\label{fig2.cluster}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{Result of Wolff Algorithm Applied to Spin Configuration in
Figure~\protect\ref{fig1.cluster}---all spins in cluster flipped to same
new value (in this case from blue to red)}
\label{fig3.cluster}
\end{figure}

\begin{figure}
\vspace*{2.0in}
\caption{Result of Metropolis Algorithm Applied to Spin Configuration in
Figure~\protect\ref{fig1.cluster}---only a few single spins flipped}
\label{fig4.cluster}
\end{figure}

First we outline a sequential method for labeling clusters, the so-called 
``ants in the labyrinth'' algorithm.  The reason for its name is that we can 
visualize the algorithm as follows \cite{Dewar:87a}.  An ant is put somewhere 
on the lattice, and notes which of the neighboring sites are connected to the 
site it is on.  At the next time-step, this ant places children on each of 
these connected sites which are not already occupied.  The children then 
proceed to reproduce likewise until the entire cluster is populated.  In 
order to label all the clusters, we start by giving every site a negative
label, set the initial cluster label to be zero, and then loop through all 
the sites in turn.  If a site's label is negative, then the site has not 
already been assigned to a cluster so we place an ant on this site, give it 
the current cluster label, and let it reproduce passing the label on to all 
its offspring.  When this cluster is identified, we increment the cluster 
label and carry on repeating the ant-colony birth, growth and death cycle 
until all the clusters have been identified.

\subsection{Parallel Cluster Algorithms}

As with the percolation models upon which the cluster algorithms are based, 
the phase transition in a spin model occurs when the clusters of bonded 
spins become large enough to span the entire lattice.  Thus, near criticality 
(which in most cases is where we want to perform the simulation), clusters 
come in all sizes, from order $N$ (where $N$ is the number of sites in the 
lattice) right down to a single site.  The highly irregular and non-local 
nature of the clusters means that cluster update algorithms do not vectorize 
well, and hence give poor performance on vector machines.  On this problem, a 
CRAY~X-MP is only about ten times faster than a Sun~4 workstation.  The 
irregularity of the clusters also means that SIMD machines are not well 
suited to this problem 
% TERRI - check the numbering on these new refs.
\cite{Baillie:91?}, \cite{Brower:91?}, \cite{Apostolakis:92?},
\cite{Apostolakis:93?},
whereas for the Metropolis type algorithms, they are perhaps 
the best machines available.  It therefore appears that the optimum 
performance for this type of problem will come from MIMD parallel computers.

A parallel cluster algorithm involves distributing the lattice onto an array 
of processors using the usual domain decomposition.  Clearly, a sequential 
algorithm can be used to label the clusters on each processor, but we need a 
procedure for converting these labels to their correct {\em global} values. 
We need to be able to tell many processors, which may be any distance apart, 
that some of their clusters are actually the same, to agree on which of the 
many different local labels for a given cluster should be assigned to be the 
global cluster label, and to pass this label to all the processors containing 
a part of that cluster.  We have implemented two such algorithms, 
``self-labeling'' and ``global equivalencing''
% TERRI - check the numbering on these new refs.
\cite{Coddington:90?}, \cite{Baillie:91?}.  

\subsection{Self-labeling}

We shall refer to this algorithm as ``self-labeling'', since each site 
figures out which cluster it is in by itself from local information.  
This method has also been referred to as ``local label propagation''
% TERRI - check the numbering on these new refs.
\cite{Brower:91?}, \cite{Tamayo:93?}.
We begin by assigning each site, $i$, a unique cluster label, $S_i$.  In 
practice, this is simply chosen as the position of that site in the lattice.
At each step of the algorithm in parallel, every site looks in turn at each 
of its neighbors in the positive directions.  If it is bonded to a 
neighboring site, $n$, which has a different cluster label, $S_n$, then both 
$S_i$ and $S_n$ are set to the minimum of the two.  This is continued until 
nothing changes, by which time all the clusters will have been labeled with 
the minimum initial label of all the sites in the cluster.  Note that to 
check termination of the algorithm involves each processor sending a 
termination flag (finished or not finished) to every other processor after 
each step, which can become very costly for a large processor array.  This is 
an SIMD algorithm and can, therefore, be run on machines like the AMT DAP and 
TMC Connection Machine.  However, the SIMD nature of these computers leads to 
very poor load balancing.  Most processors end up waiting for the few in the 
largest cluster which are the last to finish.  We implemented this on the AMT 
DAP and obtained only about 20\% efficiency.

We can improve this method on a MIMD machine by using a faster sequential 
algorithm, such as ``ants in the labyrinth'', to label the clusters in the 
sub-lattice on each processor, and then just use self-labeling on the sites 
at the edges of each processor to eventually arrive at the global cluster 
labels  
% TERRI - check the numbering on these new refs.
\cite{Coddington:90?}, \cite{Baillie:91?}, \cite{Tamayo:93?}.
The number of steps required to do the self-labeling will depend on 
the largest cluster which, at the phase transition, will generally span the 
entire lattice.  The number of self-labeling steps will, therefore, be of the 
order of the maximum distance between processors, which for a square array
of $P$ processors is just $2\sqrt{P}$.  Hence, the amount of communication 
(and calculation) involved in doing the self-labeling, which is proportional 
to the number of iterations times the perimeter of the sub-lattice, goes like 
$L$ for an $L\times L$ lattice, whereas the time taken on each processor to do 
the local cluster labeling goes like the area of the sub-lattice, which is 
$L^2/P$.  Therefore, as long as $L$ is substantially greater than the number 
of processors, we can expect to obtain a reasonable speedup.  Of course, this 
algorithm suffers from the same type of load imbalance as does the SIMD 
version.  However, in this case, it is much less severe since most of the 
work is done with ``ants in the labyrinth'' which is well load balanced.  The 
speedups obtained on the Symult~2010, for a variety of lattice sizes, are 
shown in Figure~\ref{fig5.cluster}.  The dashed line indicates perfect speedup 
({\em i.e.,} 100\% efficiency).  The lattice sizes for which we actually need 
large numbers of processors are of the order of $512^2$ or greater, and we 
can see that running on 64 nodes (or running multiple simulations of 64 nodes 
each) gives us quite acceptable efficiencies of about 70\% for $512^2$ and 
80\% for $1024^2$.

\begin{figure}
\vspace*{2.0in}
\caption{Speedups for Self-Labeling Algorithm}
\label{fig5.cluster}
\end{figure}

\subsection{Global Equivalencing}

In this method, we again use the fastest sequential algorithm to identify the 
clusters in the sub-lattice on every processor.  Each processor then looks at 
the labels of sites along the edges of the neighboring processors in the 
positive directions, and works out which ones are connected and should be 
matched up.  These lists of ``equivalences'' are all passed to one of the 
processors, which uses an algorithm for finding equivalence classes 
\cite{Knuth:68a}, \cite{Press:86a} (which, in this case, are the global 
cluster labels) to match up the connected clusters.  This processor then 
broadcasts the results back to all the other processors. 

This part of the algorithm is purely sequential, and is thus a potentially
disastrous bottleneck for large numbers of processors.   It also requires 
this processor to have a large amount of memory in which to store all the 
labels from every other processor.  The amount of work involved in doing the 
global matchup goes like $P$ times the perimeter of the sub-lattice on each
processor, i.e., like $L\sqrt{P}$ so that the efficiency should be less than 
for self-labeling, although we might still expect reasonable speedups if the 
number of processors is not extremely large.  The speedups obtained on the 
Symult~2010 for a variety of lattice sizes are shown in 
Figure~\ref{fig6.cluster}.  The figure for $512^2$ on 128 processors is 
missing due to memory constraints.  Global equivalencing gives about the same 
speedups as self-labeling for small numbers of processors, but as expected 
self-labeling does much better as the number of processors increases.

\begin{figure}
\vspace*{2.0in}
\caption{Speedups for Global Equivalencing Algorithm}
\label{fig6.cluster}
\end{figure}

\subsection{Other Algorithms}

The problem of labeling clusters of spins is an example of a standard graph
problem known as {\it connected component labeling}
% TERRI - check the numbering of this ref.
\cite{Horowitz:78?}.
Another important instance occurs in image analysis, in identifying and 
labeling the connected components of a binary or multi-colored image 
composed of an array of pixels 
% TERRI - check the numbering of this new ref.
\cite{Rosenfeld:82?}. 
There have been a number of parallel algorithms implemented for this 
problem \cite{Cypher:89a}, \cite{Embrechts:89a},  
% TERRI - check the numbering of these new refs.
\cite{Woo:89?}, \cite{Alnuweiri:92?}.
The most promising of these parallel algorithms for spin models has a 
hierarchical divide-and-conquer approach \cite{Embrechts:89a},  
% TERRI - check the numbering of this new ref.
\cite{Baillie:91?}.  
The processor array is 
divided up into smaller sub-arrays of, for example, $2\times2$ processors.  
In each subarray, the processors look at the edges of their neighbors for 
clusters which are connected across processor boundaries.  As in global 
equivalencing, these equivalences are all passed to one of the nodes of the 
sub-array, which places them in equivalence classes.  The results of these 
partial matchings are similarly combined on each $4\times4$ sub-array, and 
this process is continued until finally all the partial results are merged 
together on a single processor to give the global cluster values.

Finally, we should mention the trivial parallelization technique of running 
independent Monte Carlo simulations on different processors.  This method 
works well until the lattice size gets too big to fit into the memory of each 
node.  In the case of the Potts model, for example, only lattices of size 
less than about $300^2$ or $50^3$ will fit into 1~Mbyte, though most other 
spin models are more complicated and more memory intensive.  The smaller 
lattices which are seen to give poor speedups in Figure~\ref{fig5.cluster} 
and Figure~\ref{fig6.cluster} can be run with 100\% efficiency in this way.  
Note, of course, that this requires an MIMD computer.  In fact, we have used 
this method to calculate the dynamical critical exponents of various cluster 
algorithms for Potts models \cite{Baillie:90m}, 
% TERRI - check the numbering of these new refs.
\cite{Baillie:91?}, \cite{Coddington:92?}
(see Section~\ref{potts}).

%TERRI - should this really go here? I haven't seen any Credits
% in any other sections.
\subsection{Credits}

This research was performed by C.~F.~Baillie, P.~D.~Coddington,
J.~Apostolakis and E.~Marinari.

\newpage
\bibliography{/home/sashimi/kea/bib/doc,/home/sashimi/kea/bib/outside}
\bibliographystyle{gcf_unsrt}
\end{document}



% TERRI - here are the new refs, check the numbers (I've just put ?).
% Each have a unique name and year so there shouldn't be any confusion. 


[Coddington:90?] 
P.~D.~Coddington and C.~F.~Baillie,
``Cluster Algorithms for Spin Models on MIMD Parallel Computers'',
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).
Caltech Report C3P-862.

[Baillie:91?]
C.F.~Baillie and P.D.~Coddington,
``Cluster Identification Algorithms for Spin Models - Sequential and Parallel'',
{\it Concurrency: Practice and Experience} {\bf 3}, 129-144 (1991).
Caltech Report C3P-855.


[Brower:91?]
R.~C.~Brower, P.~Tamayo and B.~York,
``A Parallel Multigrid Algorithm for Percolation Clusters'',
{\it J. Stat. Phys.} {\bf 63}, 73 (1991).

[Flanigan:92?]
M.~Flanigan and P.~Tamayo,
``A Parallel Cluster Labeling Method for Monte Carlo Dynamics'',
To be published in {\it Int. J. Mod. Phys. C}.


[Apostolakis:92?]
J.~Apostolakis, P.~Coddington and E.~Marinari,
``A Multi-Grid Cluster Labeling Scheme'',
{\it Europhys. Lett.} {\bf 17}, 189-194 (1992).
NPAC Report SCCS-103.

[Apostolakis:93?]
J.~Apostolakis, P.~Coddington and E.~Marinari,
``New SIMD Algorithms for Cluster Labeling on Parallel Computers''
To be published in {\it Int. J. Mod. Phys. C}.
NPAC Report SCCS-279.


[Horowitz:78?]
E.~Horowitz and S.~Sahni, {\em Fundamentals of Computer Algorithms},
(Computer Science Press, Rockville, Maryland, 1978).

[Rosenfeld:82?]
A.~Rosenfeld and A.~C.~Kak, {\em Digital Picture Processing},
(Academic Press, New York, 1982).


[Baillie:91?]
C.~F.~Baillie and P.~D.~Coddington,
``Comparison of Cluster Algorithms for 2-D Potts Models'',
{\it Phys. Rev.} {\bf B 43}, 10617 (1991).
Caltech Report C3P-945, NPAC Report SCCS-12.

[Coddington:92?]
P.D.~Coddington and C.F.~Baillie,
``Empirical Relations Between Static and Dynamic Exponents for
Ising Model Cluster Algorithms'',
{\it Phys. Rev. Lett.} {\it 68}, 962-965 (1992).
NPAC Report SCCS-141.


[Woo:89?]
J.~Woo and S.~Sahni,
``Hypercube Computing: Connected Components'',
{\it J. Supercomput.} {\bf 3}, 209-234 (1989).

[Alnuweiri:92?]
H.~M.~Alnuweiri and V.~K.~Prasanna,
``Parallel Architectures and Algorithms for Image Component Labeling'',
{\it IEEE Trans. Patt. Anal. Machine Intell.} {\bf 10}, 1014-1034 (1992).

