\documentstyle[12pt]{article}       %\documentstyle[12pt,upcite]{article}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Change citations to appear as superscripts
%   ie.  upcite included here to make it self-contained.  % \uppercite
%
\catcode`\@=11 \def\@cite#1{[#1]}                  % \uppercite
\def\@cite#1{\raisebox{1ex}{{\footnotesize #1}}} % Produces the output of \cite
\catcode`\@=12                                     % \uppercite
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{titlepage}

\begin{flushright}
{\bf SCCS-279} \\[2.5em]
\end{flushright}

% double spacing
% \baselineskip=24pt

\begin{center}
{\LARGE
New SIMD Algorithms for Cluster Labeling on Parallel Computers\\}
\bigskip
{\Large
John Apostolakis, Paul Coddington \\
and Enzo Marinari$^{(*)}$\\[0.5em]}
{\large
Physics Department, \\
Syracuse University,\\
Syracuse, NY 13244, U.S.A.\\[1.0em]}
%
\bigskip
% single spacing
\bigskip
\today \\
\end{center}
\bigskip
% single spacing
\bigskip


\begin{abstract}
\noindent
Cluster algorithms are non-local Monte Carlo update sche\-mes which
can greatly increase the efficiency of computer simulations of spin models
of magnets.
The major computational task in these algorithms is connected component
labeling, to identify clusters of connected sites on a lattice. We have
devised some new SIMD component labeling algorithms, and implemented them on
the Connection Machine. We investigate their performance when applied to
the cluster update of the two dimensional Ising spin model.
These algorithms could also be applied to other problems which use connected
component labeling, such as percolation and image analysis.
\end{abstract}

\bigskip
\bigskip
% single spacing
\bigskip
\bigskip
\bigskip
\bigskip
\noindent
$^{(*)}$: and Dipartimento di Fisica,
Universit\`a di Roma {\it Tor Vergata}, Via E. Carnevale, 00173 Roma, Italy.

\end{titlepage}
% _____________________________________________________________________

% DRAFT:
\oddsidemargin  0pt
\evensidemargin 0pt
\textwidth    461pt

\begin{document}
\baselineskip=18pt

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Introduction}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Monte Carlo simulation is a very important numerical technique for
studying a wide range of problems in the physical sciences, and in
particular, the statistical mechanics of spin models of
magnets\cite{PARISI,BINDER}.
Unfortunately, traditional Monte Carlo algorithms for these models,
such as the commonly used Metropolis
algorithm\cite{METROPOLIS},
suffer from {\em critical slowing down} near the regions of interest~--~the
critical points separating different phases of the system.
This means that the autocorrelation time (the number of
iterations needed to generate a new, statistically independent, data point)
increases as $L^z$ at the critical point, where $L$ is the linear size of
the system, and $z$ is the dynamic critical
exponent\cite{SOKALREV,LAT90REV}.
$z$ is at least 2 for most local algorithms, such as Metropolis,
so the efficiency of these methods decreases rapidly as the size of
the system is increased.

The reason for this poor performance is that standard Monte Carlo algorithms
are local. In the lattice of spins which represents, for example, a
magnetic material, only a single spin at a time is changed, and this
change is influenced only by the spins on neighboring sites.
Information undergoes a random walk on the lattice, and thus takes a time
of order $L^2$ to propagate throughout the lattice.
In the last few years, algorithms have been invented for certain types of
spin models which make large-scale, non-local changes, and greatly reduce
critical slowing down\cite{SW,WOLFF}.
In these so-called {\em cluster} algorithms,
clusters of spins (rather than single spins) are changed at each step
of the Monte Carlo procedure. The clusters are formed by generating bonds
connecting neighboring sites, using a probabilistic
procedure which varies between different models and algorithms
(for reviews of cluster algorithms, see
Refs.~%
{\catcode`\@=11 \def\@cite#1{[#1]}            % \inlinecite
 \cite{SOKALREV,LAT90REV,LAT89REV,SWREV}       % \inlinecite
 \catcode`\@=12                                % \inlinecite
} ).                                           % \inlinecite

%\def\@cite#1{\raisebox{1ex}{{\footnotesize #1}}}     %\uppercite  

The major computational task of these cluster algorithms 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,IMAGE},
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.

Sequentially, this can be done in time of order $V$ (the number of vertices,
which in our case is the volume, or number of sites, in the lattice), and
consequently the cluster algorithms run about as fast as the local algorithms
%(see   Refs.~\cite{COMPSCI,PAUL,LIEDSOK} 
(see   Refs.~%
{\catcode`\@=11  \def\@cite#1{[#1]}            % \inlinecite
 \cite{COMPSCI,PAUL,LIEDSOK} 
 \catcode`\@=12                                % \inlinecite
} ).                                           % \inlinecite

for a discussion of sequential labeling algorithms).
However this may not be the case if our computer is a distributed memory
parallel machine. Local algorithms perform very efficiently on parallel
machines, whereas efficient component labeling on a parallel machine is a very
difficult problem\cite{PAUL}. Here the information concerning
the connectivity of a given physical part of the lattice is only
contained in a single processor, and obtaining information from distant
regions of the lattice (and hence also of the computer) can be very slow
if the clusters are large, and thus contain sites which are distributed over
many processors.

Let us assume that the time taken to label the clusters scales asymptotically
as $L^{d+y}$ for a lattice of $L^d$ sites, where $y$ is an exponent indicating
{\em computational slowing down}, in analogy with the dynamical exponent $z$
expressing the critical slowing down of a Monte Carlo simulation.
This means that the overall computational cost of a Monte Carlo cluster
algorithm simulation at the critical point will scale as $L^{d+y+z}$.
If we cannot find a parallel labeling algorithm for which $y$ is zero, the
advantages of cluster update algorithms over traditional local algorithms
may be eliminated on a parallel machine by the computational complexity of
labeling the clusters.

Our aim is to find a parallel component labeling algorithm with no
computational slowing down (i.e. $y=0$).
We will consider here the case where the parallel computer is a Single
Instruction Multiple Data (SIMD) machine, although the ideas described here
could also be applied to Multiple Instruction Multiple Data (MIMD) machines.
We have implemented all the algorithms on the CM-2 Connection Machine, which
is a typical massively parallel SIMD computer\cite{HILLIS}.

In order to test the algorithms we have studied the clusters formed in the
physically interesting case of the Swendsen-Wang cluster algorithm applied
to the Ising spin model at its critical point\cite{SW}.
These clusters are very similar to those created by the simple procedure
of randomly connecting neighboring sites on a two dimensional lattice with
probability $\frac{1}{2}$.
Clusters created in this way are very difficult to label efficiently on a
parallel machine, since the clusters in a particular configuration of the
connections come in many different sizes,
have extremely irregular shapes, with small clusters embedded in larger ones,
and typically including a very large cluster which will span
the lattice.
This type of problem is consequently an excellent test of parallel component
labeling algorithms.

We also note that
the worst case behavior of the labeling algorithms is not relevant for this
problem -- what we are really interested in is the {\em average}
time to label physically realistic configurations of clusters which occur in
the cluster update of the spin model.
We have therefore obtained all our data by averaging over a large number
(typically 400) of different realizations of the site connections,
taken from different Swendsen-Wang bond configurations for the Ising model
at its critical point. This is
in order to get statistically significant results from which we can obtain
the scaling behavior of our algorithms, and timings for our implementations
of these algorithms on the CM-2.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Simple Parallel Algorithms}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Local label propagation. **********************************************

The simplest and most obvious SIMD component labeling algorithm is local
label pro\-pa\-ga\-tion\cite{PAUL,BROWER}. We start with a different label
on each site, and with a list of nearest neighbor connections
(these will be Boolean variables in the following: {\em off} means no
connection is present and {\em on} means that there is a connection).
Each site then looks to each of its neighbors in turn.
If it is connected to this neighbor, and if its neighbor's label is smaller
than its own label, then it replaces its label with that of its neighbor.
This procedure is repeated until there is no change to the labels, at
which time each cluster will be labeled by the minimum initial label of all
the sites in that cluster.

% Computational slowing down. ******************************************

This local algorithm suffers from computational slowing down, and for many
problems of interest (such as spin models at their critical point)
its performance degrades very fast with increasing volume.
This is because there is typically a large cluster whose graph-theoretic
diameter or {\em chemical distance},
defined as the maximum value of the shortest path length between
two points in the cluster,
scales as $L^f$, where the exponent $f$ is approximately 1 for the
two dimensional Ising model\cite{MIRANDA}, i.e. the diameter of the
largest cluster scales approximately linearly with $L$.
For any local labeling algorithm, the minimum label has to diffuse across
this large cluster, so we expect that $y = f \approx 1$.

% Scan algorithm. *********************************************************

This algorithm can be improved by making the propagation step non-local.
One way of doing this is, instead of propagating the labels only to
neighboring connected sites, to propagate them as far as we can along a
given direction, until we come to a site with no connection in that
direction. On the CM-2 this can be done very quickly by using the intrinsic
{\bf scan\_with\_minimum} function\cite{CMDOCS}.
This routine operates on a row of numbers (labels in our case), each of
which has an associated Boolean flag (the connections).
It runs along each connected section of the row and deposits at each site
the minimum of the numbers in the section up to that point
(this is done in a distance doubling way, taking $\log_2 L$ steps).
One step of this labeling method consists of a {\it scan} in each of the
forward and backward directions of every axis. If periodic boundary
conditions are used, this must be supplemented by a local label propagation
step, since the {\it scan} routine does not wrap around the lattice.

% Bond improvement. *********************************************************

Another way of improving the above algorithms is the notion of
{\em connection improvement}. So far we have considered the bonds between
sites to be static, in other words they are set up at the beginning
and remain unchanged throughout the labeling procedure.
However it is actually very useful to change, or improve, the
connections as the labeling procedure progresses, and we learn more about the
connectivity of the sites. It will often happen that neighboring sites
will not have a bond between them, but will still be part of the same
cluster, as shown in Fig.~\ref{fig:connection_imp}.
If we compare the labels of neighboring sites
at each step of the labeling algorithm, then at some point we will find
that these two neighboring sites have the same label. We could then place
a connection between these sites, since we now know that they are in the
same cluster. Improving the bonds in this way means that new labeling
information can now flow directly between these two sites, rather than by
an indirect route via the original bonds.

Connection improvement is especially useful when applied to the {\it scan}
algorithm, since in that case the addition of extra connections means that
labels may be propagated much further in a single {\it scan} operation.
This can be seen in Fig.~\ref{fig:local_iters}, which shows a log-log plot of
the average number of iterations required to complete the cluster labeling
for the local algorithm and the {\it scan} algorithm, both with and without
connection improvement.
In Table~\ref{tab:exponents} we show the exponent $y$ for computational
slowing down for each of these algorithms, which are obtained from the
straight line fits shown in Fig.~\ref{fig:local_iters}.
As expected, the exponents are all near 1, except for {\it scan} with
connection improvement, which is substantially smaller, although still far
from zero. However we should note
that these results are very dependent on the type of bond configurations
used. Note for example that configurations for which the clusters are fully
connected, smooth, regular shapes, such as may occur in labeling objects in
image processing applications, would be labeled in a very small number of
{\it scan} operations. We might expect that $y$ would be zero for the
{\it scan} algorithm for those particular types of configurations.

%---------------------------------------------------------------------------
% \input tab1.tex
%---------------------------------------------------------------------------


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{A Multi-Scale Algorithm}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In this section we describe a regular, synchronous, multi-scale
algorithm for cluster la\-be\-ling\cite{LETTER}.
We present numerical evidence that the average number of iterations and
the average time taken
do not undergo any power-law computational slowing down (i.e. $y=0$)
for our application of labeling Ising model clusters.

The algorithm is effective on a general SIMD machine provided that the
switching network has some very basic non-local connections.
In the following we will assume that
the machine allows very fast communication between sites which
are a distance of $2^m$ sites away in any direction of the physical lattice.
These are the only non-local connections we need in order to build an
algorithm which is not affected by power-law slowing down.
Such connections would be provided, for example, by a machine with a
hypercube topology.

On the CM-2 the mapping of the physical structure of the lattice
to the (almost) hypercube processor communication network
provides specific communication to nodes of the lattice that
are at a distance of any power of two away, known as {\bf power\_of\_two}
operations.
This involves the transfer of information over not more than
two links of the hypercube, and should thus be executed at
not less than half the speed of local communications;
however the relative timing compared to a local communication depends on
the virtual processor to physical processor (VP) ratio (i.e. the number
of lattice sites per processor), as we will see later.

In common with the method proposed by Brower, Tamayo and
York\cite{BROWER},
this method uses a {\em multi-scale} approach in propagating cluster labels,
in order to overcome the slowing down inherent in local labeling algorithms.
However this algorithm is much simpler, and seems to have better scaling
properties.

Our algorithm works for a lattice of any dimensionality $d$, but for ease of
description we will consider a two dimensional problem.
In this case the key variables used are
Boolean connections that are set up in the $x$ and the $y$ axis
at a distance $2^m$, for $m=1,\ ... ,\ l-1$
(where the lattice size $L \equiv 2^l$),
by a logical {\small AND} of connections at level $m-1$.
For example, the distance 2 connection between sites $i$ and
$i + 2 \hat{x}$ is set (i.e turned {\em on}) if both the connections between
sites $i$ and $i + \hat{x}$ and sites $i + \hat{x}$ and $i + 2 \hat{x}$
are already {\em on}.
These connections are rebuilt in this way at each iteration.
In addition to building up the long distance connections in this manner,
at each iteration we also use improve the connections,
thus a connection between two sites at a generic distance $M$ which was
originally {\em off} can be set (i.e. declared to be {\em on})
if the two sites are found to have the same label.
Using connection improvement greatly reduces the number of iterations
needed to converge to the final values of the labels.

Thus, during one multi-grid label updating cycle each site will look in turn
at each of its $2d$ neighbors at each level $m$ of the multi-scale connections.
It will update, when possible, its label and also update its
connection by merging the level $m-1$ connections and by using connection
improvement.
A full cycle of the algorithm sweeps all $l$ connection levels, and a
single such cycle solves the trivial case where all connections are {\em on}.
As the labeling progresses, what happens is that an increasing fraction
of ever longer distance connections are set as sites are
recognized as belonging to the same cluster, and these connections
become fast long distance communication channels.

In Fig.~\ref{fig:global_iters} we show the average number of iterations needed
to label the Ising clusters as a function of $\log_2{L}$.
The logarithmic slowing down is very clear. We do not see any sign of
power-law behavior, or of a higher power of the logarithm.
Each iteration of the algorithm
involves a multi-grid cycle of $\log_2 L$ steps, with each step taking
approximately the same amount of time, which is proportional to $L^d / N$,
where $N$ is the number of processors ($N \le L^d$).
Thus the total CPU time goes as $L^d (\log L)^2 / N$,
or $(\log L)^2$ for a machine with $L^d$ processors.
Hence this algorithm adds only a $(\log L)^2$ term
to the overall slowing down of a spin model cluster algorithm.
The average labeling time per site as a function of $\log_2 L$
is shown in Fig.~\ref{fig:timings}.

The effect of different VP ratios means that the times for the multi-scale
algorithm on a fixed number of processors for different lattice sizes
do not scale simply as $L^d (\log L)^2 / N$.
% as we would expect for a constant VP ratio.
%
%%  eg. scan w/ improvement's time decreases linearly with log L
%% on 8K CM-2 even though the number of iterations taken to converge increases
%% as a power of L. (REALLY?
%
Firstly, we note that local operations are more efficient at higher VP ratios,
since a greater proportion of the neighboring sites will be on the same
processor, so less inter-processor communication is required.
This effect decreases as the depth of the multi-scale procedure is increased,
since more sites at a distance $2^m$ are going to be on different processors
as $m$ increases.
Eventually the communication distance will be greater than the the size of
the subdomain on each processor, so that all data must be communicated between
processors. The communication time will therefore be roughly constant at
this level and higher, and at the highest levels it is roughly independent
of the VP ratio.
(The situation is actually slightly more complicated than this, since
on the CM-2 there are 16 processors per chip, and it is inter-chip, rather
than inter-processor, communication which is costly.)
Thus the ratio of the time taken to do a step at the largest depth to
the time to do a local labeling step increases from about 2 at $L=128$
to about 6 at $L=2048$ on a $16K$ CM-2.

There is a way however to combat the higher cost of deep iterations
and significantly reduce the running time of our algorithm with only
a simple modification.
Clearly at the beginning of the labeling procedure the long distance
connections are all {\em off}, and due to the fractal structure of the
connections, it takes several iterations before a significant number of
long distance connections are generated. It is thus very useful to tailor
the number of multi-grid levels as a function of the cycle number:
a lower depth is  useful at the beginning, while using longer distance
connections is more useful towards the end of the procedure.
Fig.~\ref{fig:global_iters} shows the average iterations for the simplest case,
where the depth is constant, while Fig.~\ref{fig:timings} gives the timings
for both the full depth version and the optimized method.

In order to investigate the effect of varying the depth of the multi-grid
procedure, we have measured the number of connections at every multi-grid
level after each iteration of the full depth algorithm. We show this for a
typical configuration in Fig.~\ref{fig:multiscale_conns}. It can be seen that
the points where the different levels become useful, i.e. where there are a
reasonable number of connections (of the order of 10\%, for instance),
increases roughly linearly with the number of iterations.
We therefore chose in our modified algorithm to make the depth
a linear function of the iteration number.
Since the usefulness of a connection at a certain level depends on the
relative timings of different operations on a specific machine,
this relation must be determined empirically.

Although Fig.~\ref{fig:multiscale_conns} shows that
a large number of connections at higher levels exist,
we found that steps at the highest levels cost too much and were thus of
comparatively little use, so a maximum depth of $N_{max}$ steps was used.
We thus parametrize the depth at each iteration by
%
\begin{equation}
depth = \min \{ slope * iteration, N_{max} \} .
\end{equation}
%
and seek to find the optimum value for the $slope$ and the maximum
depth $N_{max}$.

The behavior of the average labeling time versus the slope for some different
values of $L$ and $N_{max}$ can be seen in Fig.~\ref{fig:multiscale_opt}.
A minimum exists between $0.3$ and $0.5$ in all cases.
The minimum is fairly broad and its breadth tends to increase as $L$ increases.
Fig.~\ref{fig:timings} shows the average labeling time per site for the
optimized algorithm as well as the full depth multi-scale procedure.
Note that the time for the optimized procedure is significantly smaller.
The optimal value of $N_{max}$ is $\log_2 L - 3$ (i.e. the 3 highest levels
are not used) for most lattice sizes, although for large $L$ it is
slightly more efficient to exclude the fourth highest level as well, since
the ratio of the time taken at the higher levels to the time for a
local iteration is greater at large $L$.

   We also tried a very simple telescoping scheme in order to determine
 whether any benefit could be derived from reducing the size of the problem
 in such a way. For this we used only the top level and one lower level,
 consisting of the even-even coordinate sites,
 and attempted to solve this partial problem by iterating until the labels
 of this sublattice did not change. Neither this nor an attempt to do a
 fixed number of multi-scale iterations on the lower level
 managed to reduce the amount of time required for the full labeling
 procedure. We note that only a single reduction and expansion was tried.
 This poor result was not due to the overhead of communicating between the
 different lattices, which cost very little time, but must have been due
 to the lack of information at the lower level
 about enough of the important connections.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{ Get/Send : An algorithm using general communications }

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

We now present a different component labeling algorithm\cite{JAPHCOMP}
%{\marginpar{\scriptsize JAPHCOMP}} % as a primitive DRAFT mode
that was inspired by the efficiency of the
SIMD algorithm of Hillis and Steele\cite{HILSTEEL} for finding the
end of a linked list.
We treat the labels as pointers in a dynamic tree-like structure,
making our method similar to the sequential algorithms
 of Galler and Fisher\cite{GALLFISH}
and of Hoshen and Kopelman\cite{HOSHKOPL}
and the paralell algorithm of Shiloach and Vishkin\cite{SHIVISH}.

All the cluster labeling algorithms discussed to this point start with each
site being given a unique label.
It is convenient to set the original label to be the site number.
We will again label a cluster by the smallest starting label of all
the sites in that cluster.
For a two dimensional lattice we could assign the original
cluster label $C$ of the site $(x,y)$ to be ${(y * L) + x}$, for example.
During the cluster labeling procedure we can consider the current label
as a pointer to the site where it originated (i.e.
    $x_{orig} = {\it C \ \tt mod \ } L $ ,
    $y_{orig} = C / L $ ).
The site ($x_{orig}$,$y_{orig}$) started with the current label of the
site ($x,y$), and as we shall see, 
at any later time in the labeling procedure it must have a
label which is less than or equal to this value.
Thus at any time each site can get the label of the origin
% , or root,
($x_{orig}$,$y_{orig}$) of its current label
and use it as its new label. On the Connection Machine this is
done using the general communication routine {\bf get}.
To ensure that the algorithm eventually gives the correct result,
at each stage a local label propagation step must be performed as well.

By itself this method will correctly label any lattice but performs
very badly, because in many places labels propagate
only with the local step, on paths that can be very long.
To overcome this problem we have added an important supplement --
an inverse step which propagates information large distances in the
opposite direction and proceeds as follows.
Each site saves its label, and then performs a local iteration. It then
compares its current label with that old label and, if they are different,
sends the current label to the originating site of its old label.
For this step we use the Connection Machine routine {\bf send\_with\_minimum},
for which any site that is sent more than one value keeps only the smallest.
Each site then takes the minimum of the labels it is sent, if that is
smaller than its current label.

Doing a {\it send} step before each {\it get} means that if the label of
any site is changed, the new label is then propagated immediately
(by the {\it get}) to all sites with the old label.
Thus we can wholly relabel a large area, or subcluster, in one step
as soon as it contacts another large subcluster with a smaller label.
In Fig.~\ref{fig:global_iters} we show the average number of iterations
required to label the Ising model clusters. The data fit perfectly
to a logarithmic increase with the lattice size.

The costliest parts of this algorithm are the {\it get} and {\it send}
steps, which require general communication routines which take about
ten times longer than local grid communications on the CM-2.
This is compensated of course by the value of the information which is
passed over large distances. However a way of doing some of the work by
a less costly method will reduce the total running time.
Since in the first few iterations the {\it get/send} step is able to do
very little useful work, at a large cost,
it is always faster to wait for $i$ initial iterations before using it.
Also, depending on the machine specifics and the efficiency with which
{\it get} and {\it send} are implemented,
it may be more efficient to do the {\it get/send} step only at every $p$th
iteration of the algorithm, with just local label propagation for the other
iterations.
The parameters $i$ and $p$ can be tuned to optimize the algorithm for any
particular application, system size, and parallel computer.
We found that for our application on the CM-2,
${p = 1}$ or $2$ and $i$ between 4 and 8
generally gave the best results.

We also note that for large lattices it is not necessary to {\it get}
the label of the originating site for every point since,
in any sizable cluster, most neighboring points belong to the same cluster.
We have experimented with having only a portion (we use a quarter)
of the sites perform a {\it get}. The new labels received are then
transferred to neighboring sites by the subsequent local iteration.
Whether this modification proves to be the fastest option depends of
course on the communications hardware and software of the particular machine.
For the CM-2, the time for each {\it get} is roughly halved,
although the number of iterations required is increased slightly,
so there is a trade-off. For an earlier implementation
of the {\it get} routine on the CM-2 this method was substantially faster;
however, an improved {\it get} now means that it is about $10\%$ faster to
do the {\it get} step at every site.

Another way of improving the performance of this algorithm is to combine
it with the {\it scan} operation.
Using {\em scan} in some of the initial iterations
pushes the label more efficiently over moderate distances.
On the CM-2, doing a {\em get} and a {\em send} is about as expensive as
doing the $4$ {\it scan}s (one in each direction) of a {\it scan} step,
however it provides a much better way to propagate labels
by making large changes at ever increasing length scales, and
of handling the large, irregular and labyrinthine
clusters for which the {\it scan} algorithm fares poorly.
Including a few initial {\it scan} steps makes the algorithm slightly more
efficient at large VP ratios (i.e. larger lattice sizes), but again this
will be highly dependent on the specifics of the problem and the machine.

Of course we continue to use connection improvement for the local steps,
and see a benefit for all the variations of the algorithm.
Fig.~\ref{fig:getsend_opt} shows the average time per site to label lattices of
different sizes.
It is evident that between sizes of 128 and 512 the efficiency of
large VP ratios reduces the average time, while for larger lattices this
behavior subsides and it is dominated by the increase in the number of
iterations needed to converge.

The average labeling time per site for the basic {\it get/send} algorithm
($i=0$, $p=0$, and no {\it scan}s), as well as the optimized algorithm,
is shown in Fig.~\ref{fig:timings}.
We can see that this method is substantially faster than the multi-scale
algorithm.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Discussion and Conclusions}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Our labeling algorithms are very general, and can be applied to any
application where component labeling is necessary, such as
percolation\cite{STAUFFER}, image analysis\cite{IMAGE},
and for the various cluster Monte Carlo algorithms which have been
proposed for many different spin models\cite{SOKALREV,LAT90REV,SWREV}.

We have presented numerical evidence
that the average number of iterations required by our
algorithms to label percolation-like Swendsen-Wang clusters in the
Ising model, which are highly irregular in both shape and size,
scales with the logarithm of the lattice size.
Up to corrections caused by differing VP rations on the Connection Machine,
the times required for the labeling using these algorithms scale as
$(\log_2 L)^2$ per lattice site.

There is only a subtle difference between the multi-scale and {\it scan}
algorithms: both methods look at connections at distances 1, 2, 4, etc.,
but for multi-scale we also do a comparison at each distance and set
the connection accordingly. This connection improvement is enough to
give the multi-scale algorithm significantly better scaling behavior.

Our multi-scale algorithm is simpler than that of Brower
{\it et al.}\cite{BROWER}, and appears to scale better with increasing
lattice size. For the lattice sizes of interest (of order
${1024 \times 1024}$), our optimized algorithm gives an average
labeling time for the Ising problem of 6.5 microseconds per site on a
16K CM-2 running at 7 MHz, which is comparable with the time of
6.0 $\mu$s per site obtained by Brower {\it et al.} for the same size machine.
However this time for the optimized {\it get/send} algorithm is substantially
better, at 2.6 $\mu$s per site.

These kind of SIMD algorithms work quite well on massively parallel fine
grained SIMD machines like the CM-2, as long as the
objects to be labeled are fairly small, for example in image processing
applications such as analyzing images on a radar screen.
However fine grained SIMD parallelism does not usually work well for problems
which are very irregular and require a lot of non-local communication.
Unfortunately the clusters to be labeled in spin model and percolation
applications are very large and irregularly shaped, and we would therefore
expect that it would be very hard to get good performance for labeling
algorithms on these problems using fine grained SIMD machines.
This is reflected in our results, since cluster labeling for the Ising
model can be done at a rate of about 5 $\mu$s per site on a single
IBM RS/6000-550 workstation, compared to $2.6$ $\mu$s per site
with our best algorithm on a 16K CM-2.

However we have previously obtained quite good efficiencies on coarse
grained MIMD machines for parallel component labeling algorithms which
use only local propagation of labels\cite{PAUL}, and thus do not
scale well for very large numbers of processors.
Incorporating the above multi-scale and
general communication ({\it get/send}) ideas into these
MIMD algorithms promises to allow us to greatly improve their efficiency
and scalability, and thus exploit the power of large MIMD parallel
supercomputers such as the nCUBE, the Intel machines, and the CM-5.

\bigskip
\noindent
{\it Note added:} \quad After this work had been substantially completed,
and preliminary results reported at a conference\cite{JAPHCOMP},
we found that the {\it get/send} algorithm had been independently
proposed by P.~Rossi and G.P.~Tecchiolli\cite{ROSSI}.

\bigskip
{\bf Acknowledgements}
\medskip

This work was done using Connection Machines at
the Northeast Parallel Architecture Center at Syracuse University,
Sandia National Laboratory, and Rice University.
Work supported in part by the
Center for Research on Parallel Computation with NSF cooperative
agreement No. CCR-9120008, and a grant from the IBM Corporation.

\newpage

\begin{thebibliography}{99}
%
\bibitem{PARISI}
 G.~Parisi, {\em Statistical Field Theory,} (Addison Wesley 1988).
\bibitem{BINDER}
 {\it Monte Carlo Methods in Statistical Physics, Topics in Current
 Physics 7, Second Edition}, Ed. K. Binder (Springer-Verlag, Berlin, 1986);
 {\it Applications of the Monte Carlo Method in Statistical Physics, Topics in
 Current Physics 36, Second Edition}, Ed. K. Binder (Springer-Verlag, Berlin,
 1987).
\bibitem{METROPOLIS}
 N. Metropolis {\it et al.}, {\it J. Chem. Phys.} {\bf 21}, 1087 (1953).
\bibitem{SOKALREV}
 A.~D.~Sokal, in {\em Computer Simulation Studies in Condensed Matter
 Physics:  Recent Developments}, eds.\ D.~P.~Landau {\em et al.}
 (Springer-Verlag, Berlin-Heidelberg, 1988).
\bibitem{LAT90REV}
 A.~D.~Sokal, % {\em How to Beat Critical Slowing Down: 1990 Update},
 in Proc. of the International Conference on Lattice Field Theory,
 Tallahassee, October 1990, {\it Nucl. Phys. B (Proc. Suppl.)} {\bf 20}, 55
(1991).
\bibitem{SW}
 R. H. Swendsen and J.-S. Wang, {\it Phys. Rev. Lett.} {\bf 58}, 86 (1987).
\bibitem{WOLFF}
 U. Wolff, {\it Phys. Rev. Lett.} {\bf 62}, 361 (1989).
\bibitem{LAT89REV}
 U.~Wolff,
 in Proc. of the Symposium on Lattice Field Theory, Capri, September 1989,
 {\it Nucl. Phys. B (Proc. Suppl.)} {\bf 17}, 93 (1990).
\bibitem{SWREV}
 J.-S.~Wang and R.~H.~Swendsen, {\it Physica A} {\bf 167}, 565 (1990).
\bibitem{COMPSCI}
 E.~M.~Reingold, J.~Nievergelt and N.~Deo,
 {\em Combinatorial Algorithms: Theory and Practice}
 (Prentice-Hall, Englewood Cliffs, N.J., 1977);
 E.~Horowitz and S.~Sahni, {\em Fundamentals of Computer Algorithms},
 (Computer Science Press, Rockville, Maryland, 1978).
\bibitem{IMAGE}
 A.~Rosenfeld and A.~C.~Kak, {\em Digital Picture Processing},
 (Academic Press, New York, 1982).
\bibitem{PAUL}
 C.~F.~Baillie and P.~D.~Coddington,
 {\it Concurrency: Practice and Experience} {\bf 3}, 129 (1991).
\bibitem{LIEDSOK}
 R.~G.~Edwards, X.-J.~Li and A.~D.~Sokal, {\em Sequential and Vectorized
 Algorithms for Computing the Connected Components of an Undirected
 Graph}, in preparation.
\bibitem{HILLIS}
 D.~Hillis, {\em The Connection Machine}, (MIT Press, Cambridge, Mass., 1985).
\bibitem{BROWER}
 R.~C.~Brower, P.~Tamayo and B.~York, {\it J. Stat. Phys.} {\bf 63}, 73 (1991).
\bibitem{MIRANDA}
 E.~N.~Miranda, {\it Physica A} {\bf 175}, 229 (1991).
\bibitem{CMDOCS}
 {\it Programming in Paris}, (Thinking Machines Corporation,
 Cambridge, Mass., 1989).
\bibitem{LETTER}
 J. Apostolakis, P. Coddington and E. Marinari,
 {\it Europhys. Lett.} {\bf 17}, 189 (1992).
\bibitem{JAPHCOMP}
 Talk presented at the conference ``Physics Computing '91'', San Jose,
 California, June 1991.
\bibitem{HILSTEEL}
 W.~D.~Hillis and Guy~L.~Steele~Jr., {\it Comm. of the ACM}
 {\bf 29}, 1170 (1986).
\bibitem{GALLFISH}
 B.~A.~Galler and M.~J.~Fisher, {\it Commun. ACM} {\bf 7}, 301 (1964);
 D.~E. Knuth, {\it Fundamental Algorithms}, vol. 1 of
 {\it The Art of Computer Programming} (Addison-Wesley, Reading, 1968).
\bibitem{HOSHKOPL}
 J.~Hoshen and R.~Kopelman, {\it Phys. Rev. B} {\bf 14}, 3438 (1976).
\bibitem{SHIVISH}
 Y.~Shiloach and U.~Vishkin, {\it J. Algorithms} {\bf 3}, 57 (1982).
\bibitem{STAUFFER}
 D.~Stauffer, {\it Introduction to Percolation Theory},
 (Taylor and Francis, Philadelphia, 1985).
\bibitem{ROSSI}
 P.~Rossi and G.~P.~Tecchiolli, {\it Finding Clusters in a Parallel
 Environment}, unpublished.
\end{thebibliography}

\vfill
\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{center}

{\large \bf Tables}
\bigskip
\bigskip

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{table}[ph]
\renewcommand{\arraystretch}{1.5}
\centering
\begin{tabular}{cc} \hline \hline
\quad algorithm \quad           & \quad {\it y} \quad \\ \hline
\quad local \quad               & \quad 1.08(2) \quad \\
\quad {\it scan} \quad          & \quad 1.09(3) \quad \\
\quad local improved \quad      & \quad 1.01(2) \quad \\
\quad {\it scan} improved \quad & \quad 0.84(3) \quad \\ \hline \hline
\end{tabular}
\medskip
\caption{Exponents $y$ of computational slowing down for some simple
component labeling algorithms applied to clusters of Swendsen-Wang
bonds for the Ising model.}
\label{tab:exponents}
\end{table}

\vfill
\newpage

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

{\large \bf Figure Captions}
\bigskip
\bigskip

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% THESE NEED TO BE CLOSER TOGETHER, BUT LATEX SCREWS THE SPACINGS UP.
% HOW CAN WE FIX THIS???
% THE FOLLOWING COMMANDS DON'T DO ANY GOOD:
% \textfloatsep 0.5cm
% \floatsep 0.5cm

\begin{figure}[ph]
\caption{An illustration of connection improvement. The original bonds
are shown as the thick lines. If at some point in the labeling algorithm
it is found that sites $i$ and $j$ (denoted by the filled circles)
have the same label, then a new bond (the dashed line) is introduced,
so that changes in the labels are now
propagated faster between these two points. For the multi-scale algorithm
the same idea is used, except that $i$ and $j$ do not have to be neighboring
points.}
\label{fig:connection_imp}
\end{figure}
%
\begin{figure}[ph]
\caption{A log-log plot of the average number of iterations versus lattice
size for labeling Ising model clusters using the local propagation
and {\it scan} algorithms, with and without connection improvement. The
errors are not shown, but are smaller than the points.}
\label{fig:local_iters}
\end{figure}
%
\begin{figure}[ph]
\vspace{ 0 in}
\caption{Average number of iterations versus $\log_2 L$  for labeling
Ising model clusters using the full depth multi-scale and
the {\it get/send} algorithms.}
\label{fig:global_iters}
\end{figure}
%
\begin{figure}[ph]
\caption{Average times per iteration per site versus $\log_2 L$ for 
labeling Ising model clusters
on the CM-2 using the {\it scan} with connection improvement,
multi-scale (optimized and full depth) and {\it get/send}
(with and without optimization) algorithms.}
\label{fig:timings}
\end{figure}
%
\begin{figure}[ph]
\caption{The number of connections at different levels (shown at right)
as a function of iteration number for the multi-scale algorithm, for a 
lattice of size $L=1024$.}
\label{fig:multiscale_conns}
\end{figure}
%
\begin{figure}[ph]
\caption{Average times per iteration per site on the CM-2
for labeling Ising model clusters using the multi-scale algorithm
as a function of the slope parameter,
for different lattice sizes and maximum depths. }
\label{fig:multiscale_opt}
\end{figure}
%
\begin{figure}[ph]
\caption{Average times per iteration per site as a function of lattice size
for labeling Ising model clusters on the CM-2,
for different values of the parameters for the {\it get/send} algorithm.
$i_{send}$ is the number of steps without a {\it get} or {\it send}.
Here the interval between {\it get}s is $p_{get}=1$, while 
$p_{send}=p_{scan}=2$. $l_{scan}$ denotes the number of steps after which 
{\it scan}s are not performed.}
\label{fig:getsend_opt}
\end{figure}

\end{center}

\end{document}

