Parallel Cluster Identification for Multidimensional Lattices.
AbstractThe cluster identification problem is a variant of connected component labeling that arises in cluster algorithms for spin models in statistical physics. We present a multidimensional version of Belkhale and Banerjee's Quad algorithm for connected component labeling on distributed memory parallel computers. Our extension abstracts away extraneous spatial connectivity information in more than two dimensions, simplifying implementation for higher dimensionality. We identify two types of locality present in cluster configurations, and present optimizations to exploit locality for better performance. Performance results from 2D, 3D, and 4D Ising model simulations with Swendson-Wang dynamics show that the optimizations improve performance by 20-80 percent.
Introduction
The cluster identification problem is a variant of connected component labeling that arises
in cluster algorithms for spin models in statistical mechanics. In these applications, the
graph to be labeled is a d-dimensional hypercubic lattice of variables called spins, with edges
(bonds) that may exist between nearest-neighbor spins. A cluster of spins is a set of spins
defined by the transitive closure of the relation "is a bond between". Cluster algorithms
require the lattice to be labeled such that any two spins have the same label if and only if
they belong to the same cluster.
Since the cluster identification step is often the bottleneck in cluster spin model applica-
tions, it is a candidate for parallelization. However, implementation on a distributed memory
parallel computer is problematic since clusters may span the entire spatial domain, requiring
global information propagation. Furthermore, cluster configurations may be highly irregular,
preventing a priori analysis of communication and computation patterns. Parallel algorithms
for cluster identification must overcome these difficulties to achieve good performance.
We have developed a multidimensional extension of Belkhale and Banerjee's Quad algorithm
[1, 2], a 2D connected component labeling algorithm developed for VLSI circuit
extraction on a hypercube multiprocessor. This paper presents performance results from
applying the algorithm to Ising model simulations with Swendson-Wang dynamics [3] in
2D, 3D, and 4D. Our extension abstracts away extraneous spatial information so that distributed
data structures are managed in a dimension-independent manner. This strategy
considerably simplifies implementation in more than two dimensions. To our knowledge,
this implementation is the first parallelization of cluster identification in 4D.
To improve performance, we identify two types of locality present in Swendson-Wang
cluster configurations and present optimizations to exploit each type of locality. The optimizations
work with an abstract representation of the spatial connectivity information, so
they are no more complicated to implement in d ? 2 dimensions than in 2D. Performance
results show that the optimizations effectively exploit cluster locality, and can improve performance
by 20-80% for the multidimensional Quad algorithm.
The remainder of his paper proceeds as follows. Section 2 discusses previous approaches
to the cluster identification problem on parallel computers. Section 3 describes the Ising
model and the Swendson-Wang dynamics. Section 4 reviews Belkhale and Banerjee's Quad
algorithm and presents extensions for more than two dimensions. Section 5 presents two
optimizations to exploit cluster locality, and section 6 gives performance results in 2D, 3D,
and 4D.
Related Work
Several algorithms for 2D cluster identification on distributed memory MIMD computers
have been presented in recent years.
Flanigan and Tamayo present a relaxation algorithm for a block domain decomposition
[4]. In this method, neighboring processors compare cluster labels and iterate until a
steady state is reached. Baillie and Coddington consider a similar approach in their self-
labeling algorithm [5]. Both relaxation methods demonstrate reasonable scaleup for 2D
problems, but for critical cluster configurations the number of relaxation iterations grows
as the distance between the furthest two processors (2
P for block decompositions on P
processors). Other approaches similar to relaxation have been presented with strip decompositions
[6, 7]. Strip decompositions result in only two external surfaces per processor.
However, the distance between two processors can be as large as P , which increases the
number of stages to reach a steady state. Multigrid methods to accelerate the relaxation
algorithm for large clusters have been presented for SIMD architectures [8, 9].
Host-node algorithms involve communicating global connectivity information to a single
processor. This host processor labels the global graph and then communicates the results to
other processors. Host-node algorithms [10, 11, 5] do not scale to more than a few processors
since the serialized host process becomes a bottleneck.
Hierarchical methods for connected component labeling are characterized by a spatial
domain decomposition and propagation of global information in log P stages. Our approach
is based on the hierarchical Quad algorithm for VLSI circuit extraction on a hypercube
multiprocessor [1]. Other hierarchical methods for distributed memory computers have been
used for image component labeling [12, 13]. Baillie and Coddington consider a MIMD
hierarchical algorithm for the Ising model, but do not achieve good parallel efficiency [5].
Mino presents a hierarchical labeling algorithm for vector architectures [14].
There has been comparably little work evaluating MIMD cluster identification algorithms
in more than two dimensions. Bauernfeind et al. consider both relaxation and a host-
node approaches to the 3D problem [15]. They introduce the channel reduction and net list
optimizations to reduce communication and computation requirements in 3D. They conclude
that the host-node approach is inappropriate for 3D due to increased memory requirements
on the host node.
Fink et al. present 2D and 3D results from a preliminary implementation of the multidimensional
Quad algorithm [2]. This paper includes 4D results and introduces issues
pertaining to a dimension-independent implementation.
Ising Model
Many physical systems such as binary fluids, liquid and gas systems, and magnets exhibit
phase transitions. In order to understand these "critical phenomena," simple effective models
have been constructed in statistical mechanics. The simplest such model, the Ising model,
gives qualitative insights into the properties of phase transitions and sometimes can even
provide quantitative predictions for measurable physical quantities [16].
The Ising model can be solved exactly in 2D [17]. In more than two dimensions, exact
solutions are not known and numerical simulations are often used to obtain approximate
results. For example, numerical simulations of the 3D Ising model can be used to determine
properties of phase transitions in systems like binary liquids [18]. The 4D Ising model is
a prototype of a relativistic field theory and can be used to learn about non-perturbative
aspects, in particular phase transitions, of such theories [19].
In d dimensions, the Ising model consists of a d-dimensional lattice of variables (called
spins) that take discrete values of \Sigma1. Neighboring spins are coupled, with a coupling
strength - which is inversely proportional to the temperature T .
Monte Carlo simulations of the Ising model generate a sequence of spin configurations.
In traditional local-update Monte Carlo Ising model simulations, a spin's value may or may
not change depending on the values of its neighbors and a random variable [5]. Since each
spin update depends solely on local information, these algorithms map naturally onto a
distributed memory architecture.
The interesting physics arises from spin configurations in the critical region, where phase
transitions occur. In these configurations, neighboring spins form large clusters in which all
spins have the same value. Unfortunately, if - is the length over which spins are correlated
(the correlation length), then the number of iterations required to reach a statistically independent
configuration grows as - z . For local update schemes the value z (the dynamical
critical exponent) is z - 2. Thus, even for correlation lengths - as small as 10 to 100, critical
slowing-down severely limits the effectiveness of local-update algorithms for the Ising
model [20].
In order to avoid critical slowing-down, Swendson and Wang's cluster algorithm updates
whole regions of spins simultaneously [3]. This non-local update scheme generates
independent configurations in fewer iterations that the conventional algorithms. The cluster
algorithm has a much smaller value of z, often approaching 0. Therefore, it eliminates critical
slowing-down completely. The Swendson-Wang cluster algorithm proceeds as follows:
1. Compute bonds between spins. A bond exists with probability
adjacent spins with the same value.
2. Label clusters of spins, where clusters are defined by the transitive closure of the
relation "is a bond between".
3. Randomly assign all spins in each cluster a common spin value, \Sigma1.
These steps are repeated in each iteration.
On a distributed memory computer, a very large spin lattice must be partitioned spatially
across processors. With a block decomposition, step 1 is simple to parallelize, since we only
compute bonds between neighboring spins. Each processor must only communicate spins on
the boundaries to neighboring processors. The work in step 3 is proportional to the number
of clusters, which is typically much less than the number of lattice sites.
Step 2 is the bottleneck in the computation. A single cluster may span the entire lattice,
and thus the entire processor array. To label such a cluster requires global propagation
of information. Thus the labeling step is not ideally matched to a distributed memory
architecture, and requires an efficient parallel algorithm.
4.1 2D Quad Algorithm
Our cluster identification method is based on Belkhale and Banerjee's Quad algorithm for
geometric connected component labeling [1], which was developed to label connected sets
of rectangles that represent VLSI circuits in a plane. It is straightforward to apply the
same algorithm to label clusters in a 2D lattice of spin values. A brief description of the
Quad algorithm as applied to a 2D lattice of spins is presented here. For a more complete
description of the Quad algorithm, see [1].
The cluster labeling algorithm consists of a local labeling phase and a global combining
phase. First, the global 2D lattice is partitioned blockwise across the processor array. Each
processor labels the clusters in its local partition of the plane with some sequential labeling
algorithm. The Quad algorithm merges the local labels across processors to assign the correct
global label to each spin site on each processor.
On P processors, the Quad algorithm takes log P stages, during which each processor
determines the correct global labels for spins in its partition of the plane. Before each stage,
each processor has knowledge of a rectangular information region that spans an ever-growing
section of the plane. Intuitively, processor Q's information region represents the portion of
the global domain from which Q has already collected the information necessary to label Q's
local spins. The data associated with an information region consists of
ffl A list of labels of clusters that touch at least one border of the information region.
These clusters are called CCOMP sets.
ffl For each of the four borders of the information region, a list representing the off-
processor bonds that touch the border.
Each bond in a border list connects a spin site in the current information region with a
spin site that is outside the region. Each bond is associated with the CCOMP set containing
the local spin site. The border list data structure is a list of offsets into the list of CCOMP
set labels, where each offset represents one bond on the border. This indirect representation
facilitates Union-Find cluster mergers, which are described below (see figure 1).
l 1
l
l
l
l
l
CCOMP Labels
Border Bond Lists

Figure

1: Fields of an information region data structure.
The initial information region for a processor consists of the CCOMP set labels and
border lists for its local partition of the plane. At each stage, each processor Q 1 exchanges
messages with a processor Q 2 such that Q 1 and Q 2 's information regions are adjacent. The
messages contain the CCOMP set labels and border lists of the current information region.
Processor merges the CCOMP sets on the common border of the two information regions
using a Union-Find data structure [21]. The other border lists of the two information regions
are concatenated to form the information region for processor Q 1 in the next stage. In this
manner, the size of a processor's information region doubles at each stage so after log P
stages each processor's information region spans the entire plane. Figure 2 illustrates how
the information region grows to span the entire global domain.
For a planar topology, a processor's global combining is complete when its information
region spans the entire plane. If the global domain has a toroidal topology, clusters on
opposite borders of the last information region are merged in a post-processing step.
current information region
stage 1 stage 2
stage 4
stage 3
partner information region
done

Figure

2: Information regions in each stage of the Quad algorithm. There are sixteen
processors. At each stage, the information region of the processor in the top left corner is
the current information region. The partner processor and its information region are also
shown. In each stage, the two information regions are merged, forming the information
region for the subsequent stage.
4.2 Extending the Quad Algorithm to Higher Dimensions
A straightforward extension of the Quad algorithm to more than two dimensions results in
fairly complex multidimensional information region data structures. To simplify implemen-
tation, we present a multidimensional extension using an abstract dimension-independent
information region representation.
The divide-and-conquer Quad algorithm strategy can be naturally extended to d ? 2
dimensions by partitioning the global domain into d-dimensional blocks, and assigning them
one to a processor. Each processor performs a sequential labeling method on its local domain,
and then the domains are translated into information regions for the global combining step.
An information region represents a d-dimensional subset of the d-dimensional global domain.
These d-dimensional information regions are merged at each stage of the algorithm, so after
log P stages the information region spans the entire global domain.
In two dimensions, the list of bonds on each border is just a 1D list, corresponding to the
1D border between two 2D information regions. Since bonds do not exist at every lattice
site, the border lists are sparse. For a 3D lattice, the border lists must represent sparse 2D
borders. In general, the border between two d-dimensional information regions is a d \Gamma 1-
dimensional hyperplane. Thus a straightforward 3D or 4D implementation would be much
more complex than in two dimensions, because sparse multidimensional hyperplanes must
be communicated and traversed in order to merge clusters.
To avoid this complication, note that if we impose an order on the bonds touching an
information region border, the actual spatial location of each bond within the border is not
needed to merge sets across processors. As long as each processor stores the border bonds
in the same order, we can store the bonds in a 1D list and merge clusters from different
processors by traversing corresponding border lists in order. Figure 3 illustrates this for 3D
lattices. This concept was first applied by Fink et al. to the 3D Quad algorithm [2], and a

Figure

3: The 2D borders of a 3D information region are linearized by enumerating the
border bonds in the same order on each processor.
similar optimization was applied to 3D lattices by Bauernfeind et al.[15].
We define an order on the border bonds by considering each (d \Gamma 1)-dimensional border
as a subset of the d-dimensional global lattice. Enumerate the bonds touching a
dimensional border in column-major order relative to the d-dimensional global lattice. Since
each processor enumerates the sites relative to the same global indices, each processor stores
the sets on a border in the same order, without regard to the orientation of the border in
space.
This ordering linearizes (d \Gamma 1)-dimensional borders, resulting in an abstract information
region whose border representations are independent of the problem dimensionality. When
two of these information regions are merged, the order of bonds in the new border lists is
consistent on different processors. Therefore, the logic of merging clusters on a border of
two information regions does not change for multidimensional lattices. No sparse hyperplane
data structures are required, and a 2D cluster implementation can be extended to 3D and
4D with few modifications.
4.3 Performance analysis
Belkhale and Banerjee show that the 2D Quad algorithm for VLSI circuit extraction runs in
is the number of processors, ff() is the inverse of
Ackerman's function, t s is the message startup time, t b is the communication time per byte,
and B is is the number of border rectangles along a cross section of the global domain [1].
The number of border rectangles in VLSI circuit extraction applications corresponds to the
number of border bonds in cluster identification applications. For cluster identification on a
lattice, let N be the lattice size and p be the probability that there is a bond between
two adjacent lattice points. Then
giving a running time of O(log P t s
N)).
For a d-dimensional problem, define N and p as above. Assume the global domain is
a d-dimensional hypercube with sides N 1=d , which is partitioned onto the d-dimensional
logical hypercube of processors with sides P 1=d . Suppose at stage i a processor's information
region is a hypercube with sides of length a. Then at stage i + d the information region is
a hypercube with sides of length 2a. Thus, the surface area of information region increases
by a factor of 2 d\Gamma1 every d stages. Let b(i) be the surface area of the information region at
stage i. Then b(i) is at most
d
e (1)
It is easy to see that
d . The total number of bonds on the border of an information
region is proportional to the surface area. Summing over log P stages, we find the total
number of bytes that a processor communicates during the algorithm is O(2dpN
d ). There
are log P message starts, so the total time spent in communication is O(log
The total number of Union-Find operations performed by a processor at each stage is equal
to the number of bonds on a border of the information region. Using the path compression
and union by rank optimizations of Union-Find operations [21], the total work spent merging
clusters is O(pdN
d ff(pdN
d )). (Our implementation uses the path compression heuristic
but not union-by-rank.) Adding together communication and computation, the running
time for global combining is O(log
d ff(pdN
d )).
Breadth-First Search(BFS) has been shown to be an efficient algorithm to perform the
sequential local labeling step [5]. Since BFS runs in O(jV j+jEj) [21], the local labeling phase
runs in O(( N
)). Thus, for any dimension lattice, the time for the local phase will
dominate the time for the global phase as long as N is large. However, as d increases, the
global time increases relative to the local time for a fixed problem size. We must therefore
scale the problem size along with the problem dimensionality in order to realize equivalent
parallel efficiency for higher dimensional lattices.
Optimizations
One limitation of the Quad algorithm is that the surface area of the information region
grows in each stage. By the last stage, each processor must handle a cross-section of the
entire global domain. With many processors and large problem sizes, this can degrade the
algorithm's performance [1]. To mitigate this effect, we have developed optimizations that
exploit properties of the cluster configuration for better performance.
In Monte Carlo Ising model simulations, the cluster configuration structure depends
heavily on the coupling constant -. Recall that the probability that a bond exists between
two adjacent same-valued spins is . For subcritical (low) -, bonds are relative
sparse and most clusters are small. For supercritical (high) -, bonds are relatively dense
and one large cluster tends to permeate the entire lattice. At criticality, the system is in
transition between these two cases, and the cluster configurations are combinations of small
and large clusters.
How any particular spin affects the labels of other spins depends on the cluster configuration
properties. We identify the following two types of locality that may exist in a cluster
configuration:
clusters only affect cluster labels in a limited area.
ffl Type 2: Adjacent lattice points are likely to belong the same cluster.
Subcritical configurations exhibit Type 1 locality, and supercritical configurations exhibit
Type 2 locality. Configurations at criticality show some aspects of both types.
Belkhale and Banerjee exploit Type 1 locality in two dimensions with the Overlap Quad
algorithm [1]. In this algorithm, information regions overlap and only clusters that span the
overlap region must be merged. Intuitively, small clusters are eliminated in early stages of the
algorithm, leaving only large clusters to merge in later stages. The Overlap Quad algorithm
requires that the positions of bonds within borders be maintained, precluding use of an
abstract dimension-independent information region data structure. Instead, we present two
simpler optimizations, Bubble Elimination and Border Compression. These optimizations
work with the abstract border representations, so they are no more complicated to implement
in d ? 2 dimensions than in 2D.
5.1 Bubble Elimination
Bubble Elimination exploits Type 1 locality by eliminating small clusters in a preprocessing
phase to the Quad algorithm. A local cluster that touches only one border of the information
region is called a bubble. Immediately after initializing its information region, each processor
identifies the bubbles along each border. This information is exchanged with each neighbor,
and clusters marked as bubbles on both sides of a border are merged and deleted from the
borders. Thus, small clusters are eliminated from the information regions before performing
the basic Quad algorithm. During the course of the Quad algorithm, communication and
computation is reduced since the bubble clusters are not considered.
Bubble elimination incurs a communication overhead of 3 d \Gamma1 messages for a d-dimensional
problem. If we communicate with Manhattan neighbors only, the communication overhead
drops to 2d messages. Although bubbles on the corners and edges of an information region
are not eliminated, this effect is insignificant if the granularity of the problem is sufficiently
large.
5.2 Border Compression
Border Compression exploits Type 2 locality by changing the representation of the border
lists. We compress the representation of each list using run-length encoding [22]. That is,
a border list of set labels is replaced by a sequence of pairs ((l
where s(l i ) is the number of times value l i appears in succession in a border list.
If Type 2 locality is prevalent, border compression aids performance in two ways: it
reduces the length of messages, and we can exploit the compressed representation to reduce
the number of Union-Find operations that are performed. Before two compressed borders are
merged, they are decompressed to form two corresponding lists of cluster labels to combine.
From the compressed representation, it is simple to determine when two clusters are merged
together several times in succession. During decompression, it is simple to filter redundant
mergers out of the lists, reducing the number of Union-Find mergers to be performed. Thus,
border compression reduces both communication volume and computation.
For some cluster configurations, bubble elimination increases the effectiveness of border
compression. Suppose the global cluster configuration resembles swiss cheese, in that there
are many small clusters interspersed with one large cluster. This phenomenon occurs in Ising
model cluster configurations with - at or above criticality. Bubble elimination removes most
small clusters during preprocessing, leaving most active clusters belonging to the one large
cluster. In this case, there will be long runs of identical labels along a border of an information
region. Border compression collapses these runs, leaving small effective information region
borders.
6 Performance Results
6.1 Implementation
We have implemented the cluster algorithm and Ising model simulation in 2D, 3D, and
4D with C++. The global lattice has a toroidal topology in all directions. When using
bubble elimination, only Manhattan neighbors are considered. The local labeling method is
Breadth-First Search [21].
According to the Swendson-Wang algorithm, clusters must be flipped randomly after
the cluster identification step. For a spatially decomposed parallel implementation, it is
necessary that all processors obtain consistent pseudorandom numbers when generating the
new spins for clusters that span more than one processor. In our implementation, we generate
the new random spins for each local cluster prior to global cluster merging, and store the
new spin as the high-order bit in the cluster label. Thus, after cluster merging, all spins in
a cluster are guaranteed to be consistent.
To simplify implementation in more than 2 dimensions, we use the LPARX programming
library, version 1.1 [23]. LPARX manages data distribution and communication between
Cartesian lattices, and greatly simplifies the sections of the code that manages the regular
spin lattice. The kernel of the cluster algorithm is written in a message-passing style using the
message-passing layer of LPARX[24], a generic message-passing system resembling the
Message Passing Interface[25]. Since the cluster algorithm is largely dimension-independent,
the message-passing code is almost identical for each problem dimensionality. In fact, the
same code generates the 2D,3D, and 4D versions; compile-time macros determine the problem
dimensionality.
The code was developed and debugged on a Sun workstation using LPARX mechanisms
to simulate parallel processes and message-passing. All performance results were obtained
from runs on an Intel Paragon under OSF/1 R1.2, with 32MB/node. The code was compiled
with gcc v.2.5.7 with optimization level 2.
6.2 Performance
The total cluster identification time consists of a local stage, to perform local labeling with
a sequential algorithm, and a global stage, to combine labels across processors using the
multidimensional Quad algorithm. All times reported are wall clock times and were obtained
with the Paragon dclock() system call, which has an accuracy of 100 nanoseconds.
Intuitively, we expect the benefits from bubble elimination and border compression to
vary with -, the coupling constant. Figures 4, 5, and 6 show the global stage running times
at varying values of -. For the problem sizes shown, the critical region occurs at - c - 0:221
in 2D, - c - 0:111 in 3D, and - c - 0:08 in 4D.
Since the surface area-to-volume ratio is larger in 3D and 4D than in 2D, the optimizations
are more important for these problems. As expected, figures 5 and 6 show that bubble
elimination is effective in the subcritical region, and border compression is effective in the
supercritical region. In the critical region, the two optimizations together are more effective
than either optimization alone. Presumably this is due to the "swiss cheese" effect discussed
in Section 5. Together the optimizations improve performance by 35-85%, depending on -,
in both 3D and 4D.
Kappa1.03.05.07.0
Global
Combining
Time
Per
4D Global Combining Time
Nodes, 68x68x34x34 Lattice
Bubble Elimination
Border Compression
Both
In 2D, the optimizations improve performance by 20-70%, but do not show the intuitive
dependence on - as in 3D and 4D. We suspect this is due to cache effects. As - increases,
the number of global clusters decreases. Thus, during cluster merging, more union-find
data structure accesses will be cache hits at higher - since a greater portion of the union-
find data structure fits in cache. In 2D, the surface area-to-volume ratio is low, so these
union-find accesses become the dominant factor in the algorithm's performance. In 3D and
4D, information region borders are much larger, overflowing the cache and causing many
more cache misses traversing the borders of the information region. Since these borders are
larger than the union-find data structures, union-find data structure memory accesses are
less critical.

Figure

7 shows the relative costs of the local stage and global stages with and without
optimizations. The breakdown shows that in 2D, the local labeling time dominates the global
time, so the benefit from optimizations is limited by Amdahl's law [26]. However, in 3D and
4D, the global stage is the bottleneck in the computation, so the two optimizations have a
significant impact on performance.
Timing results are instructive, but depend on implementation details and machine ar-
chitecture. To evaluate the optimizations with a more objective measure, Table 1 shows the
total number of bytes transmitted in the course of the Quad algorithm. Since the amount
of work to merge clusters is directly proportional to the length of messages, these numbers
give a good indication of how successfully the optimizations exploit various cluster configu-
rations. The communication volume reduction varies depending on the cluster configuration
structure, ranging up to a factor of twenty to one.
Since physicists are interested in using parallel processing capability to run larger problems
than possible on a single workstation, it is appropriate to examine the algorithm's
performance as the problem size and number of processors grow. For an ideal linear parallel
Time
per
Site
(ns)
Normalized Labeling Perfomance
Local Labeling
Global Labeling
Optimization
Both
Optimizations
subcritical critical supercritical

Figure

7: Breakdown of algorithm costs, normalized per spin site. All runs here are with 64
processors of an Intel paragon. The lattice sizes are 4680x4680 in 2D, 280x280x280 in 3D,
and 68x68x34x34 in 4D. For subcritical runs, in 3D, and 0.04 in 4D. For
critical runs, in 2D, 0.111 in 3D, and 0.08 in 4D. For supercritical runs,
in 2D, 0.2 in 3D, and 0.2 in 4D.
Opt. Elim. Compress
2D 4680x4680 lattice
3D 280x280x280 lattice
4D 68x68x34x34 lattice

Table

1: Total number of bytes transmitted during global combining. All runs are with 64
processors.
algorithm, if the problem size and number of processors are scaled together asymptotically,
the running time remains constant. Due to the global nature of the cluster identification
problem, the basic Quad algorithm cannot achieve ideal scaled speedup in practice. Since
the Quad algorithm takes log P stages, the global work should increase by at least log P . A
further difficulty is that in d dimensions, the work in the last stage of the algorithm doubles
every d stages.
However, the bubble elimination and border compression optimizations vastly reduce the
work in later stages of the algorithm. Thus, with the optimizations, we can get closer to
achieving ideal scaled speedup. Table 2 shows these scaled speedup results for a fixed number
of spins sites per processor for critical cluster configurations. The results show that as the
number of processors and problem size are scaled together, the performance benefit from the
optimizations increases. In 2D, the scaled speedup with the optimizations is nearly ideal.
The 3D and especially 4D versions do not scale as well, although figure 7 shows that better
performance is achieved away from criticality.
Although the optimizations were developed with the multidimensional Quad algorithm
in mind, we conjecture that they would also be effective for other cluster identification
algorithms, such as relaxation methods [4, 6, 7]. The multidimensional Quad algorithm and
optimizations may be also be appropriate for other variants of connected component labeling.
One open question is whether the border compression and bubble elimination optimizations
would effectively exploit the graph structure of other applications, such as image component
labeling applications.
7 Conclusion
We have presented an efficient multidimensional extension to Belkhale and Banerjee's Quad
algorithm for connected component labeling. Since the extension deals with abstract spatial
Number of No Optimizations Both Optimizations
Processors 2D 3D 4D 2D 3D 4D

Table

2: Global combining time, in seconds, when the lattice size and number of processors
is scaled together. Each processor's partition is 585x585 in 2D, 70x70x70 in 3D, and
17x17x17x17 in 4D. All runs are at - c .
connectivity information, distributed data structures are managed in a dimension-independent
manner. This technique considerably simplifies implementations in more than two dimen-
sions. We introduced two optimizations to the basic algorithm that effectively exploit locality
in Ising model cluster configurations. Depending on the structure of cluster configurations,
the optimizations improve performance by 20-80% on the Intel Paragon. With the opti-
mizations, large lattices can be labeled on many processors with good parallel efficiency.
The optimizations are especially important in more than two dimensions, where the surface
area-to-volume ratio is high.



--R

"Parallel algorithms for geometric connected component labeling on a hypercube multiprocessor,"
"Cluster identification on a distributed memory multiprocessor,"
"Nonuniversal critical dynamics in monte carlo simu- lations,"
"A parallel cluster labeling method for monte carlo dy- namics,"
"Cluster identification algorithms for spin models - sequential and parallel,"
"Parallelization of the Ising model and its performance evaluation,"
"Swendson-wang dynamics on large 2d critical Ising models,"
"A multi-grid cluster labeling scheme,"
"A parallel multigrid algorithm for percolation clusters,"
"Parallel simulation of the Ising model,"
"Paralleliza- tion of the 2d swendson-wang algorithm,"
"Evaluation of connected component labeling algorithms on shared and distributed memory multiprocessors,"
"Component labeling algorithms on an intel ipsc/2 hypercube,"
"A vectorized algorithm for cluster formation in the Swendson-Wang dynam- ics,"
"3D Ising model with swendson-wang dynamics: A parallel approach,"
Statistical Field Theory.
"Crystal statistics. i. a two-dimensional model with an order-disorder tran- sition,"
"Numerical investigation of the interface tension in the three-dimensional Ising model,"
"Broken phase of the 4-dimensional Ising model in a finite volume,"
Computer simulation methods in theoretical physics.
Introduction to Algorithms.

"A robust parallel programming model for dynamic, non-uniform scientific computation,"
"The LPARX user's guide v2.0,"
Message Passing Interface Forum
Computer Architecture A Quantitative Approach.
--TR

--CTR
Scott B. Baden, Software infrastructure for non-uniform scientific computations on parallel processors, ACM SIGAPP Applied Computing Review, v.4 n.1, p.7-10, Spring 1996
swendson-wang dynamics;parallel algorithm;cluster identification;ising model;connected component labeling