
\input vanilla.sty
\font\tenbf=cmbx10
\font\tenrm=cmr10
\font\tenit=cmti10
\font\ninebf=cmbx9
\font\ninerm=cmr9
\font\nineit=cmti9
\font\eightbf=cmbx8
\font\eightrm=cmr8
\font\eightit=cmti8
\font\sevenrm=cmr7
\TagsOnRight
\hsize=5.0truein
\vsize=7.7truein
\parindent=15pt
\nopagenumbers
\baselineskip=10pt
\centerline{\eightrm International Journal of Modern Physics C\hfil}
\centerline{\eightrm $\copyright$\, World
Scientific Publishing Company \hfil}
\vglue 5pc
\baselineskip=13pt
\centerline{\tenbf A PARALLEL CLUSTER LABELING METHOD}
\centerline{\tenbf FOR MONTE CARLO DYNAMICS}
\vglue 1cm
\centerline{\tenrm Mike Flanigan and Pablo Tamayo}
\baselineskip=12pt
\centerline{\eightit Thinking Machines Corp.}
\baselineskip=10pt
\centerline{\eightit 245 First. St. Cambridge, MA 02142 U.S.A}
\vglue 1cm
\centerline{\eightrm Received (........)}
\baselineskip=10pt
\centerline{\eightrm Revised (........)}
\vglue 0.5cm

\centerline{\eightrm ABSTRACT}
{\rightskip=1.5pc \leftskip=1.5pc \eightrm\baselineskip=10pt\noindent We
present a parallel algorithm for cluster dynamics to efficiently simulate
large systems on MIMD computers with large numbers of processing nodes. The
method divides physical space into cells which are assigned to processing
nodes and combines a local (serial) cluster finding procedure with a global
relaxation process: first clusters are labeled locally inside each cell
(processing node) and then the nodes exchange labels with neighbors until a
fixed point is reached. By controlling overhead and reducing
inter-processor communication this method attains good performance and
speed-up. We analyze the complexity and scaling properties of the algorithm
and present timings for the two-dimensional Swendsen-Wang Ising dynamics.
Typical updating times for large systems are 93 nanosecs/site using 256
processing nodes and 0.36 $\mu$secs/site using 64 processing nodes on
CM-5 parallel supercomputers.  Lattices up to 12288$^2$ were simulated
using 256 processing nodes.

\vglue 0.3cm
\noindent
{\eightit Keywords}\/: Ising Model, Cluster Labeling, Percolation, Monte Carlo Simulations, Accelerated Dynamics.
\vglue 0.5cm}
\baselineskip=13pt

\vglue 12pt
\line{\tenbf 1. Introduction\hfil}
\vglue 5pt

Over the last five years, since the introduction of the Swendsen-Wang (SW)
algorithm,$^1$ cluster dynamics have been extended to a variety of models
in Statistical Mechanics and Field Theory$^2$. These new ``accelerated
dynamics'' update percolation clusters instead of single spins producing
faster decorrelation times and reduced critical slowing down$^3$. In fact,
the dynamic critical exponent $z$ is reduced significantly or even
eliminated in some cases. This acceleration does not come for free; the
price to pay is an increased computational complexity compared with
standard Metropolis, heat bath or other local algorithms. The cluster
labeling procedure requires a fair amount of non-local (unstructured) data
movements that makes cluster algorithms intrinsically hard to vectorize or
paralellize. It is important to develop efficient cluster labeling methods
because they are used every time step inside the core of the simulations.
One has to worry not only about their complexity and scaling properties but
also about absolute execution times. One of the modern challenges of
computational science is to find efficient algorithms to exploit the
unprecedented computational power of today's massively parallel
supercomputers.

Cluster labeling on a lattice is also relevant for the analysis of
structures obtained in many computer simulations of statistical systems:
Ising and percolation clusters,$^4$ nucleation droplets, polymers,
crystals, fractal structures, and particle tracks in experimental high
energy Physics. Cluster labeling is a special case of the more general
problem of finding the connected components of a graph,$^5$ which has
applications in computer vision, image processing and network analysis,
among others.

\headline={\ifodd\pageno\rightheadline \else\leftheadline\fi}
\def\rightheadline{\eightit\hfil A Parallel Cluster Labeling 
Method for Monte Carlo Dynamics\quad
\eightrm\folio}
\def\leftheadline{\eightrm\folio\quad\eightit M. Flanigan \&
 P. Tamayo \hfil}
\voffset=2\baselineskip

In this paper we propose a simple method for cluster labeling that can be
used at the core of Monte Carlo simulations on
MIMD\footnote"$^i$"{\eightrm\baselineskip=10pt Multiple Instruction
Multiple Data} parallel computers and in particular the CM-5 (Thinking
Machines Corp.).  The algorithm we will present is a general method for
cluster finding in $n-$dimensional Euclidean lattices, but we will
concentrate on the specific problem of cluster labeling for the Ising Model
with Swendsen-Wang dynamics.  Percolation bonds are defined between aligned
spins with bond probability $p_{bond} = 1 - e^{-2\beta}$, and the clusters
of connected spins, the Coniglio-Klein$^{6,7,8}$ percolation clusters, are
flipped with 50\% probability. At the critical point, where most of the
interesting behavior occurs, the clusters span the system and information
has to be propagated across the entire computational domain.  

\midinsert\vskip 8.5cm
{\rightskip=1.5pc \leftskip=1.5pc \eightrm\baselineskip=10pt\noindent
Fig. 1.\,\,\, Typical critical clusters in the
two-dimensional Ising model. The lines show the partitioning of the system
into 16 cells (processing nodes).}
\endinsert

Several methods to find connected components have been introduced in the
Computer Science literature. There are general methods based on
union-find,$^{9,10}$ transitive closure,$^{11}$ vertex
collapse$^{11,12,13,14}$ and vector models$^{15}$. We realized that most of
these methods are not well suited to practical lattice Monte Carlo
simulations on MIMD machines where one has to worry about absolute
execution times over average configurations instead of worse cases. Most of
these methods are designed for idealized
PRAM\footnote"$^{ii}$"{\eightrm\baselineskip=10pt Parallel Random-Access
Machine.} models and require extensive global communications, shared memory
or data format transformations. Some other methods are more appropriate for
SIMD\footnote"$^{iii}$"{\eightrm\baselineskip=10pt Single Instruction
Multiple Data.} fine-grained machines$^{15,16,17}$.  The choice of data
partitioning is very important$^{18,9}$. A number of MIMD cluster labeling
algorithms, specifically designed for cluster Monte Carlo simulations have
appeared in the literature,$^{19,21}$ but have very limited speed-up and
efficiency which limits their application to the simulation of large
systems using large numbers of processing nodes.  Recently, Kertesz and
Stauffer$^{20}$ have introduced a new geometric parallelization algorithm
for large systems ($6400^2$) on the Intel iPSC/860.

The method we will describe here is based on a rectangular domain
decomposition strategy\footnote"$^{iv}$"{\eightrm\baselineskip=10pt Similar
partitioning techniques have been used by Tucker,$^{18}$ Embrechts {\it
\eightrm et al},$^{22}$ and Baillie and Coddington$^{21}$}: the clusters
are first labeled locally inside each processing node and then labels are
propagated across processing nodes by a relaxation process.  In section two we
describe the algorithm and discuss some of its properties.  In section
three we will analyze its time complexity and scaling properties.
Numerical results are discussed in section four.  Finally, section five
contains conclusions.

\midinsert\vskip 8cm
{\rightskip=1.5pc \leftskip=1.5pc \eightrm\baselineskip=10pt\noindent
Fig. 2.\,\,\, (a) shows the initial state 
and connectivity of a sub-system on one processing node. (b) shows the
local roots (in black) and the sites that point to them. (c) illustrates
the data structures setup for the relaxation cycle.}
\endinsert

\vglue 12pt
\line{\tenbf 2. Description of the Algorithm\hfil}
\vglue 5pt

Physical space is divided into rectangular cells in such way that each cell
is assigned to one processing node (see Figure 1). The algorithm labels the
clusters in two stages: first it finds all the clusters inside each
processing node using a serial local algorithm, and then it performs a global
relaxation process in which processing nodes exchange clusters labels
with nearest neighbors until a fixed point is reached. The operations of
the algorithm are shown in Figures 2 and 3. The procedure can be described
as follows:

\vglue 0.35cm
{\tenbf Procedure Cluster-Labeling:}\quad
\halign{\indent#\hfil&#\hfil \cr

(i) &Define connectivity for the sites: for Swendsen-Wang dynamics connectivity\cr
    &bonds are defined with probability $p_{bond} = 1 - e^{-2\beta}$
between aligned spins.\cr

(ii) &In each processing node independently apply the serial algorithm to find\cr
     &clusters (Procedure Local). At boundary sites the off-node bonds are ignored.\cr 
     &At the end all sites are labeled with their corresponding ``local roots'',\cr 
     & which are then globalized\footnote"$^{v}$"{\eightrm\baselineskip=10pt i.e. made unique
over the whole system. This is done by adding an offset equal to the
processing node number times the size of the local system ($n$). For the
Swendsen-Wang  dynamics we also multiply the result by two and add a random
bit -- in this way we can ``piggy-back'' flipping information in the parity
of the cluster label.}.\cr

(iii) &Iterate relaxation cycles, exchanging local root labels with neighboring\cr
      &processing nodes until there is no change in any node (Procedure
Relax).\cr
      & At the end all sites get their final global label from their local roots.\cr

(iv) &Clusters are flipped with 50\% probability and measurements of
relevant\cr
     &quantities are accumulated (energy, magnetization etc.).\cr
}


\midinsert\vskip 8cm
{\rightskip=1.5pc \leftskip=1.5pc \eightrm\baselineskip=10pt\noindent
Fig. 3.\,\,\, Processing nodes exchange
local root labels with nearest neighbors. Minimum labels are found and
local root labels are updated. The process is repeated until there are no
more changes in local root labels.}
\endinsert

The serial local algorithm we used has similarities to both
Hoshen-Kopelman$^{23}$ and union-find$^{10}$ algorithms. This particular
serial algorithm proved efficient for our approach but other serial
algorithms can be used.  It can be briefly described as follows:

\vglue 0.35cm
{\tenbf Procedure Local:}\quad

\halign{\indent#\hfil&#\hfil \cr
(i) &For each site make a table of connected neighbors for directions east and\cr
    &south, and assign a unique label to each site (e.g. the array index of the\cr
    &site itself).\cr

(ii) &For each site, follow the label paths for the site and each of its connected\cr
     &neighbors until their roots (sites pointing to themselves) are found. Then\cr
     &obtain the minimum label and set the site label, its root, and the roots of\cr
     &the connected neighbors to the minimum value.\cr

(iii) &To finish the labeling process a final collapse of trees is done by a pass\cr
      &over all sites setting each site to point to its corresponding root.\cr
}

\vglue 0.35cm

After the local labeling is completed (see Fig. 2b) a number of relaxation
cycles are executed to label the clusters globally. The relaxation
procedure consists of the following:


\vglue 0.35cm
{\tenbf Procedure Relax:}\quad

\halign{\indent#\hfil&#\hfil \cr
(i) &Execute a preparation step to set data structures: for each node boundary\cr
    &define a list of pointers to local roots (see Fig. 2c), one pointer per each bond\cr
    &crossing node boundaries\footnote"$^{vi}$"{\eightrm\baselineskip=10pt 
To improve efficiency we also eliminate redundant connections across
processing node boundaries.}.\cr

(ii) &Execute a sequence of relaxation cycles until all nodes detect no change\cr
     &in the labels (see Fig. 3): \cr

{} {} {} {} (a) &{} {} Each processing node interchanges boundary labels with the neighboring\cr
       &{} {} node in each direction. The labels are sent in a single block of data using\cr  
       &{} {} synchronous message-passing calls.\cr

{} {} {} {} (b) &{} {} The local root labels are compared with the ones obtained from the\cr
       &{} {} neighbors and then set to the minimum values.\cr
}
\vglue 0.35cm


We find that even for large numbers of processing nodes the execution time is
dominated by the local part. Originally we intended to implement a
multi-grid scheme for the relaxation part, as was used in ref. 24.  This
extension is straightforward and makes the algorithm more scalable but we
have found that in practice it is not necessary given the sizes of
coarse-grained MIMD machines available today. Even a large parallel
computer with 1024 processing nodes forms relatively small lattices in 2
and 3 dimensions which would require only a few multigrid levels.
Nearest-neighbor relaxation is very efficient for this problem at the scale
of the processing nodes grid.

\vglue 12pt
\line{\tenbf 3. Analysis: Complexity and Scaling\hfil}
\vglue 5pt

Considering the fact that the cluster labeling algorithm operates at the
core of equilibrium Monte Carlo simulations the analysis will focus on
average case instead of worst case performance. The probability of
obtaining a given configuration of clusters is determined by the Boltzmann
weight of that particular configuration. The probability of observing the
worst case is negligible. The execution times we will consider correspond
to averages over equilibrium configurations generated in the Monte Carlo
process. Furthermore, we will analyze the properties of the algorithm at
the critical point for the Ising model.

In this section we will make a simple model for the time complexity of the
algorithm based on simple scaling arguments. The predictions of the model
will be compared with experimental results in the next section.

In the following discussion, $N = L \times L$ is the total size of the
system, $n = l \times l$ is the size for the subsystem in each processing
node, and $p$ is the total number of processing nodes. Clearly the total number of
sites is $N = np$.  In addition, $a$ and $b$ are constants. We start by
considering the total time to perform the cluster labeling as consisting
of two contributions: a ``local time'', the time spent by the serial
algorithm inside each processing node, and a ``relax time'' which is the
time spent in the global relaxation procedure until completion.  The
scaling of the local part is basically $O(n\; log^*n)$ where $log^*n$ is a
very slow growing function\footnote"$^{vii}$"{\eightrm\baselineskip=10pt
$\log^*n$ is equal to the number of times the $\log$
function has to be applied recursively to the argument until it converges
to one.}. In practice, as it is discussed in ref, 10, the value of $log^*n$
can be safely considered a constant and the serial time is basically
$O(n)$,

$$ T_{local} \; = \; a l^2 \; = \; a n. \,\, \tag 1$$

To compute the relaxation time we need to compute two contributions:
$t_{relax}$, the time to do one relaxation cycle and $n_{relax}$, the number of
relaxation cycles needed to complete the process. $t_{relax}$ is proportional to
the size of the boundary between processing nodes, therefore, assuming the
communication times grow linearly with message size, we have,

$$ t_{relax} \; \sim \; l \; \sim \; n^{1/2}. \,\, \tag 2$$

If we consider relaxation as a nearest-neighbor propagation of labels, then
$n_{relax}$ should be proportional to the maximum depth of the clusters
embedded in the lattice. The maximum depth of the clusters is equivalent to
the maximum shortest path joining two sites over the set of clusters:
$\lambda_{min}$; this length, also known as the chemical distance$^{25}$,
characterizes the way information is transmitted inside a cluster by
step-by-step processes such as nearest neighbor label propagation.
Typically for percolation clusters, $\lambda_{min}$ scales with cluster
linear size $r$ with a characteristic exponent $d_{min}$,

 
$$ \lambda_{min} \; \sim \; r^{d_{min}}, \,\, \tag 3$$

\noindent
which implies that the number of relaxation cycles to label a given cluster
will scale with $d_{min}$. Since the clusters can be as large as the entire
system $(r = L)$, then $n_{relax}$ will be proportional to $L^{d_{min}}$.
However, in our case we have relaxation only at the scale of the processing
nodes, not at the scale of single sites, and then the correct length to use
in the scaling expression for $n_{relax}$ is the renormalized length $L/l =
p^{1/2}$,

$$ n_{relax} \; \sim \; (L/l)^{d_{min}} \; \sim \; p^{d_{min}/2}. \,\, \tag
4$$
 

The connectivity at lower scales is ``integrated out'' by the local procedure
but we are assuming that the chemical distance of the resulting
renormalized lattice still scales with the same exponent $d_{min}$. The
value of $d_{min}$ for 2d Coniglio-Klein clusters is reported to be$^{26}$
$1.08 \pm 0.01$; consequently we expect $n_{relax}$ to grow slightly faster
than the square root of the number of
processing nodes\footnote"$^{viii}$"{\eightrm\baselineskip=10pt The use of a
multigrid or hierarchical scheme can make $n_{relax}$ scale as $\log^2 p$,
see for example refs. 24 and 27.}. Now we can express the total relaxation
time, $T_{relax}$, as

$$ T_{relax} \; = \; n_{relax}\;  t_{relax} \; = \; b p^{d_{min}/2} n^{1/2}, \,\, \tag
5$$
\noindent
where we can see that $d_{min}$ plays the role of a ``computational''
critical slowing down exponent for the algorithm. This is one of the cases
in which a physical parameter of the simulated system, such as the chemical
length exponent, determines its computational properties.

The total time for the parallel algorithm is then the sum of local plus relax times,


$$ T_{parallel} \; = \; T_{local} + T_{relax}\; = \; a n + b p^{d_{min}/2}
n^{1/2}. \,\, \tag 6$$

If $n$ is large compared with $p$, and the communication to computation
ratio $a/b$ is small, the scaling will be dominated by the local part. This
can be seen more explicitly by calculating the speed-up function $S$, which
is the ratio between the serial and parallel times,


$$ S \; = \; { T_{serial} \over T_{parallel}}. \,\, \tag 7$$


We will use as $T_{serial}$ the time of the local algorithm running on only
one processing node.  We can express $S$ as a function of $n$ and $p$ in
the following way,


$$ S(n, p)  \; = \; { anp \over an + b p^{d_{min}/2} n^{1/2}}. \,\, \tag 8$$


As expected the speed-up improves with large $n$ and gets worse as $p$
increases. Usually, the speed-up is computed as a function of $p$ for fixed
system size $N$,

$$ S_N(p) \; = \; { aNp \over aN + b p^{(d_{min} + 1)/2}
N^{1/2}}. \,\, \tag 9$$

The corresponding efficiency $E_N(p) = S_N(p)/p$ is

$$ E_N(p) \; = \; { 1 \over 1 + {b \over a} [p^{(d_{min} + 1)} N^{-1}]^{1/2}}.
\,\, \tag 10$$


The efficiency decreases as the number of processing nodes increases because
inter-processor communication times will eventually dominate. In practice it is
important to make $b/a$ as small as possible. As we will see in
the next section, this can be done effectively on the CM-5 where $b/a
\simeq 3/4$. It is interesting to notice that the efficiency is a universal
function of $p^{(d_{min}+1)} N^{-1}$. In terms of efficiency the only
important parameter in the simulation is the value of
$p^{(d_{min}+1)}N^{-1}$: the larger $N$ is in relation to $p$ the better the
algorithm will perform.

\midinsert\vskip 8cm
{\rightskip=1.5pc \leftskip=1.5pc \eightrm\baselineskip=10pt\noindent
Fig. 4.\,\,\, Local, relax and total times for
2d SW dynamics as a function of $N^{1/2}$ for fixed $p = 256$.}
\endinsert

\vglue 12pt
\line{\tenbf 4. Numerical Results\hfil}
\vglue 5pt

We have implemented the algorithm in standard ANSI C plus message-passing
calls using the CM-5 CMMD library. We have performed simulations of the 2d
Ising model with Swendsen-Wang dynamics at the critical point and measured
execution times for different values of $n$ and $p$, and total system sizes
$N$ from $128^2$ to $12288^2$, using up to $256$ processing nodes on CM-5
machines without vector units.  Tables 1 and 2, and Figures 4 and 5 show
times (local, relax, and total) for different system sizes and numbers of
processing nodes.  Measurements times (energy, magnetization, etc.) were
not included in the timings. As we can see in Figures 4 and 5 the local
times dominate for large systems and we obtain very good performance: $0.36
\mu$secs/site for a $64$ node CM-5 and $93$ nanosecs/site on a $256$ node
CM-5.  


\bigskip
\centerline{\eightrm Table 1. Timings for 2d SW dynamics, $p$ = 64
processing nodes.}
$$\vbox{\eightrm\halign{\hfil#\quad\hfil &#\hfil\quad
 &\hfil#\hfil\quad
 &\hfil#\hfil\quad &#\hfil\quad &#\hfil\cr
\multispan5\hrulefill\cr
$N$ &\hfil local [secs] \hfil &\hfil relax [secs] \hfil &\hfil 
total [secs] \hfil &\hfil $\mu$secs/site\hfil\cr
\multispan5\hrulefill\cr
128$^2$ & 0.006 & 0.008 & 0.014 & 0.854\cr
256$^2$ & 0.020 & 0.009 & 0.029 & 0.442\cr
512$^2$ & 0.081 & 0.012 & 0.093 & 0.355\cr
1024$^2$ & 0.359 & 0.019 & 0.378 & 0.360\cr
2048$^2$ & 1.545 & 0.038 & 1.583 & 0.377\cr
4096$^2$ & 6.170 & 0.072 & 6.242 & 0.372\cr
\multispan5\hrulefill\cr}}$$

\smallskip
\centerline{\eightrm Table 2. Timings for 2d SW dynamics, $p$ = 256
processing nodes.}
$$\vbox{\eightrm\halign{\hfil#\quad\hfil &#\hfil\quad
 &\hfil#\hfil\quad
 &\hfil#\hfil\quad &#\hfil\quad &#\hfil\cr
\multispan5\hrulefill\cr
$N$ &\hfil local [secs] \hfil &\hfil relax [secs] \hfil &\hfil 
total [secs] \hfil &\hfil $\mu$secs/site\hfil\cr
\multispan5\hrulefill\cr
256$^2$ & 0.006 & 0.014 & 0.020 & 0.305\cr
512$^2$ & 0.021 & 0.017 & 0.038 & 0.145\cr
1024$^2$ & 0.082 & 0.019 & 0.101 & 0.096\cr
2048$^2$ & 0.370 & 0.034 & 0.404 & 0.096\cr
4096$^2$ & 1.573 & 0.062 & 1.635 & 0.097\cr
8192$^2$ & 6.265 & 0.077 & 6.342 & 0.095\cr
12288$^2$ & 13.927 & 0.1075 & 14.0345 & 0.093\cr
\multispan5\hrulefill\cr}}$$

\midinsert\vskip 8cm
{\rightskip=1.5pc \leftskip=1.5pc \eightrm\baselineskip=10pt\noindent
Fig. 5.\,\,\, Local, relax and total times for
2d SW dynamics as a function of $p$ for fixed $n$.}
\endinsert

Simulating a $12288^2$ system requires about $15$ Mbytes of local memory on
a $256$ node machine. Using $32$ Mbytes of memory systems as large as
$16384^2$ can be simulated. Bit packing techniques will make possible to
store larger systems but the overhead will increase due to the additional
complexity and bookkeeping. We have not explored this possibility. The use
of the vector units\footnote"$^i$"{\eightrm\baselineskip=10pt The CM-5
vector units also support integer operations.} will increase the speed of
the local procedure significantly.

\midinsert\vskip 8cm
{\rightskip=1.5pc \leftskip=1.5pc \eightrm\baselineskip=10pt\noindent
Fig. 6.\,\,\, The number of relaxation
cycles as a function of the number of processing nodes for fixed $n=128^2$. The
dot-dashed line gives the asymptotic slope predicted by $d_{min}/2$.}
\endinsert

 The scaling behavior of the measured times agrees well with the simple
scaling model of the previous section (Equations 1-5): local times scale
linearly with $n$ (see Fig. 4), and relaxation times with $n^{1/2}$ (see
Fig. 5). The number of relaxation cycles as a function of $p$ for fixed $n$
approaches the asymptotic slope $d_{min}/2$ as can be seen in Fig. 6. From
these data we obtain $b/a \simeq 3/4$.  Figure 7 shows the total time
as a function of temperature where the peak corresponds to the critical
point. 


\midinsert\vskip 8cm
\centerline{\eightrm\baselineskip=8pt Fig. 7.\,\,\, Total time as a function of inverse temperature $\beta$.}
\endinsert

The speed-up $S_N(p)$ as a function of the number of processing nodes for
different system sizes $N$ is shown in Figure 8. Notice that for large
system sizes the speed-up keeps increasing even for our maximum $p$ equal
to $256$ processing nodes. The algorithm has much better speed-up than
previous methods\footnote"$^i$"{\eightrm\baselineskip=10pt For example
compare our Fig. 8 with Fig. 5 of ref. 19 (first paper) or Figures 4, 5 and
6 of ref. 21.}  The functional form of these curves is given by Eq. 9.


\midinsert\vskip 9cm
\centerline{\eightrm\baselineskip=8pt Fig. 8.\,\,\, Speed-up data for up to 256 processing nodes for different system sizes.}
\endinsert

\midinsert\vskip 9cm
\centerline{\eightrm\baselineskip=8pt Fig. 9.\,\,\, Efficiency data for up to 256 processing nodes for different system sizes.}
\endinsert

The efficiency $E_N(p) = S_N(p)/p$ as a function of $p$ is shown in Fig. 9.
Simulations of large systems ($ > 1024^2$) attain efficiencies higher than
80\% even for $256$ processing nodes.  When these points are plotted as a
function of $[p^{(d_{min}+1)}N^{-1}]^{1/2}$ (see Fig. 10) they show the
scaling described by Eq. 10: data for different system sizes
collapses reasonably well to a universal curve. In this plot, as well as in
the others, there are fluctuations in the timings that cause deviations
from the expected behavior. This is due to cache memory effects and
non-linear behavior in the message-passing communications. We have not
tried to take these effects into account in the scaling model because they
are hard to estimate and will make the analysis unnecessarily complicated.


\midinsert\vskip 8cm
{\eightrm\baselineskip=8pt Fig. 10.\,\,\, Universal scaling of 
efficiency data plotted as a function of $[p^{(d_{min}+1)}N^{-1}]^{1/2}$.
The solid line is the predicted behavior given by Equation 10 (with $a/b =
3/4$).}
\endinsert

\vglue 12pt
\line{\tenbf 5. Conclusions\hfil}
\vglue 5pt

We have found that rectangular domain decomposition and a combination of a
local method with global relaxation produces very good results in terms of
performance and speed-up. The method shows much better scaling properties
than others$^{19,21}$. A careful choice of data
structures, partitioning and communication strategies is fundamental to
produce practical and efficient algorithms for Monte Carlo simulations. A
simple scaling model for the complexity of the algorithm gives a reasonable
prediction of observed times, speed-up and efficiency.


It is interesting to notice that the most efficient cluster algorithms, SIMD
and MIMD, employ relaxation or relaxation plus multigrid methods: Brower
{\it et al}$^{24}$ (1.5 $\mu$sec/site, 64K CM-2), Baillie and
Coddington$^{21}$ (1 $\mu$sec/site, Symult 192 nodes), Apostolakis {\it et
al}$^{27,28}$ (2.5 $\mu$secs/site, 16K CM-2) and this work (93
nanosecs/site, CM-5 256 nodes). The fastest algorithm using geometric
parallelization with strips we are aware of is the one introduced by Kertesz and
Stauffer$^{20}$.

Our algorithm running on a 256 node CM-5 is about 68 times faster than a
Swendsen-Wang program$^{29}$ on a Cray-XMP (6.5 $\mu$secs/ site), which clearly
indicates that MIMD parallel computers can be used very effectively for
this kind of computational problems.

We hope our cluster labeling method will provide a useful tool for the
study of large systems and the collection of statistics over long runs. It
can be applied to Swendsen-Wang simulations of Ising and Potts
models,$^{1,2}$ embedded dynamics simulations of Landau-Ginzburg and
Heisenberg models,$^2$ virtual bond percolation dynamics$^{30}$ and the
study of static and dynamic properties of percolation clusters$^{4,6}$.


\vglue 12pt
\line{\tenbf Acknowledgments\hfil}
\vglue 5pt

We want to thank L. Tucker, C. Feynmann, D. Stauffer, G. Blelloch, A.
Greenberg, R. Giles, R.  Brower, W. Klein, N. Copty, B. Boghosian, G.
Batrouni, R. Jones, M. Drumheller, G. Drescher, R. Frye, M. Best, P.
Coddington, J. Apostolakis and E. Marinari, for interesting discussions
about cluster methods and their applications. We also thank the members of
the CM-5 CMMD team for their technical advice and useful suggestions. We
express our appreciation to C. Madsen, J. Mesirov, L. Tucker and J. Mucci
for supporting this project. Some of the simulations were done on CM-5
machines at the Pittsburgh Supercomputer Center (PSC) and the National
Center for Supercomputer Applications (NCSA).

\vglue 12pt
\line{\tenbf References\hfil}
\vglue 5pt

\medskip
\ninerm
\baselineskip=11pt
\frenchspacing

\item{1.} R. Swendsen and J. S. Wang,  {\nineit Phys. Rev. Lett.}, {\ninebf
58,}  (1987) 86.

\item{2.} R. Brower and P. Tamayo, {\nineit Phys. Rev. Lett.}, {\ninebf
62},  (1989) 1087; R. Brower and S. Huang, {\nineit Phys. Rev. D } {\ninebf
41},  (1990) 708; U. Wolff, {\nineit Phys. Rev. Lett.} {\ninebf 60}, 
(1988) 1461; R.  Edwards and A. Sokal, {\nineit Phys.  Rev. D} {\ninebf 38},
 (1988) 2009; W.  Klein, T. Ray and P. Tamayo, {\nineit Phys.  Rev. Lett.},
{\ninebf 62},  (1989) 163; P.  Tamayo, R. C.  Brower and W. Klein, {\nineit J. of
Stat. Phys.}, {\ninebf 58},  (1990) 1083;

\item{3.} P. C. Hohenberg and B. Halperin, {\nineit Rev. of Mod. Phys.},
{\ninebf 49},  (1977) 435.

\item{4.} D. Stauffer and A. Aharony, {\nineit Introduction to Percolation 
Theory 2nd. ed.,} Taylor and Francis, London (1992).

\item{5.} T. H. Cormen, C. E. Leiserson and R. L. Rivest, {\nineit Introduction to
Algorithms}, MIT Press, (1990); D. Knuth, {\nineit The Art of Computer
Programming,} vol 3, Addison-Wesley (1973).

\item{6.} A. Coniglio and W. Klein, {\nineit J. Phys. A,} {\ninebf 13},
(1980)   2775.

\item{7.} D. Stauffer, {\nineit Physica A} {\ninebf 171},  (1991) 471; D. Stauffer and J. Kertesz,
{\nineit Physica A} {\ninebf 177}, (1991) 381; E. N. Miranda, {\nineit
Physica A} {\ninebf 175}, (1991) 235 , {\ninebf 179}, (1991) 340. E. N.
Miranda, {\nineit Physica A} {\ninebf 175}, (1991) 229.

\item{8.} A. Coniglio, {\nineit Phys. Rev. Lett.}, {\ninebf 62},  (1989) 3054.


\item{9.} J. Woo and S. Sahni, {\nineit Jour. of Supercomputing} {\ninebf
3},  (1989) 209.

\item{10.} E. Horowitz and S. Sahni, {\nineit Fundamentals of Computers
Algorithms, Potomac, Md. Computer Science Press, 1978.}  

\item{11.} M. J. Quinn and N. Deo, {\nineit Computing Surveys,} Vol. 16, No 3, September 1984.

\item{12.} Y. Shiloach and U. Vishkin, {\nineit Jour. of Algorithms}
{\ninebf 3}, (1982)  57.

\item{13.} U. Vishkin, {\nineit Discrete Applied Mathematics} {\ninebf 9},
(1984) 197.

\item{14.} P. S. Gopalakrishnan, I. V. Ramakrishnan and L. N. Kanal, 1985
{\nineit IEEE, Int. Conf. on Parallel Processing.}

\item{15.} G. H. Blelloch, {\nineit Vector Models for Data-Parallel
Computing},  MIT Press, 1990.

\item{16.} R. Cypher, J. L. C. Sanz and L. Snyder, {\nineit Journal of
Algorithms} {\ninebf 10}, (1989) 140.

\item{17.} W. Lim, A. Agrawal and L. Nekludova, {\nineit Thinking Machines Tech. Report} NA86-2.

\item{18.} L. W. Tucker, Proc. {\nineit IEEE Conference on Computer Vision and Pattern
Recognition,} June 1986, Miami, Florida.

\item{19.} A. N. Burkitt and D. W. Heermann, {\nineit Comp. Phys. Comm.}
{\ninebf 54}, (1989) 210; D. W. Heermann and A. N. Burkitt, {\it Parallel
Algorithms in Computational Science}, Springer Verlag, Heidelberg 1991.

\item{20.} J. Kertesz and D. Stauffer, private communication.

\item{21.}  C. F. Baillie and P. D. Coddington, {\nineit Concurrency: Practice and
Experience} {\ninebf 3}(2), (1991)  129.

\item{22.} H. Embrechts, D. Roose and P. Wambacq, {\nineit Hypercube and Distributed
Computers,} F. Andre and J. P. Verjus eds.  Elsevier Science Pubs. B. V.
(North-Holland) 1989.

\item{23.}  J. Hoshen and R. Kopelman, {\nineit Phys. Rev. B,} {\ninebf
14}, (1976)  3438. 

\item{24.} R. C. Brower, P. Tamayo and B. York, {\nineit Jour. of Stat.
Phys.} {\ninebf 63}, (1991)  73.

\item{25.} H. J. Herrmann and H. E. Stanley, {\nineit J. Phys. A} {\ninebf
21}  (1988) L829.

\item{26.} E. N. Miranda, {\nineit Physica A} {\ninebf 175},  (1991) 229.

\item{27.}  J. Apostolakis, P. Coddington and E. Marinari, {\nineit
Europhys. Lett.} {\ninebf 17}(3)  (1992) 198;

\item{28.}  J. Apostolakis, P. Coddington and E. Marinari, private communication.

\item{29.}  U. Wolff, {\nineit Phys. Lett. B} {\ninebf 228}, (1989)  379 .

\item{30.}  R. Brower and P. Tamayo, Virtual Bond Percolation for Ising Cluster
Dynamics, {\nineit Thinking Machines Corp. tech.  rep.} TMC 231 (1992).



\vfil\supereject
\bye

.




