
% Figure for Burkitt-Heermann efficiencies ??

% Figure for breadth-first search?

% Figure for local label propagation??

% Figure for multigrid ??

% Compare block and cyclic performance for parallel Wolff?? 
%  Row cyclic instead of column cyclic for Fortran instead of C?

% Performance plot ??


\section{Parallel Cluster Algorithms}


\begin{slide}

\subsection{\bf Parallel Cluster Algorithms}

For the Swendsen-Wang cluster algorithm, all sites need to calculate their
cluster labels, so can use standard domain decomposition.

The measurements can be done in parallel, just as for Metropolis. 

Updating the spins can easily be done in parallel, 
as can setting up the bonds --- all the sites just look at their 
nearest neighbors in the positive directions, and put a bond between them
with probability 
$p = \left(1-e^{-\beta}\right)\,\delta \left( s_i,s_j \right)$ for 
the Potts model. 

The only difficult problem in parallel is the connected component 
labeling, that is, assigning to all the sites a cluster label that 
is unique for each cluster. So implementing an efficient parallel
Swendsen-Wang cluster algorithm is basically equivalent to the problem of
implementing an efficient parallel connected component labeling 
algorithm.

\end{slide}


\begin{slide}

\subsection{\bf The Need for a Parallel Algorithm}

The reason we want to use cluster algorithms is that they offer such a large 
reduction in computation time over the traditional Metropolis algorithm.  
For example, for the Ising model on a $512^2$ lattice, the number of 
Metropolis iterations required to generate an independent configuration is of 
order 1000, whereas for the cluster algorithms it is of order 10, i.e., a 
100-fold improvement in speed.  

However, early attempts at running cluster algorithms on parallel 
supercomputers gave speed-ups only of order 10, 
whereas the Metropolis algorithm can easily have speed-ups of 
100--1000 on vector and parallel supercomputers.

So on modern supercomputers, the improvement of the cluster algorithm is 
lost unless an efficient parallel algorithm can be found.

\end{slide}


\begin{slide}

% Burkitt-Heermann efficiencies -- scan the figure??

% \centerline{\psfig{figure=SIMD-iterations.eps,height=5in}}
% \vspace{-1cm}
{\Smallsize
The first attempts at parallelizing irregular cluster algorithms got 
poor speedups compared to regular Metropolis algorithms.
Figure taken from A.N.~Burkitt and D.W.~Heermann,
``Parallelization of a Cluster Algorithm'',
{\it Comp. Phys. Comm.} {\bf 54}, 210 (1989).
}

\end{slide}


\begin{slide}

\subsection{\bf Regular and Local vs. Irregular and Non-Local}

The same type of problem occurs in many different guises in
computational physics.  
In many applications, there are standard algorithms that are regular, simple, 
often local, and thus parallelize very efficiently.  However, it is this 
very regularity and locality that makes them poor algorithms.  

Complex physical problems and physical systems are usually irregular, 
often non-local, and change rapidly with time, so simple regular and 
local algorithms tend to have problems, such as critical slowing down. 

More complex algorithms, such as 
non-local, irregular cluster algorithms\@; multiscale and multigrid methods 
for PDEs and spin models\@; adaptive, irregular grids for finite element 
calculations\@; hierarchical, adaptive N-body solvers\@; etc., all work much 
better, but are much more difficult to parallelize!

\end{slide}


\begin{slide}

\subsection{\bf Labeling or Growing a Single Cluster}

First look at how sequential component labeling algorithms work.  
Start with the Wolff algorithm, which ``grows'' a single cluster.  
Can liken this to the growth of an ant colony, in the 
following manner:

\begin{minipage}[h]{5in}
{\parskip=0pt
\Smallsize
\begin{enumerate}
\item
Pick a site at random and place an ant on it. 
\item
The ant reproduces by placing a child on the four neighboring sites 
(in 2-$d$) with probability $p$.  
\item
This ``generation'' of ants also looks at each of its four neighbors in 
turn, and if it is not already occupied, places a child of its own there with 
probability $p$.  
\item
The ant colony continues to grow until
the final generation produces no children.
\end{enumerate}
}
\end{minipage} \ \
%
\begin{minipage}[h]{3in}
\begin{center}
\begin{picture}(173,280)
%%BoundingBox: 0 0 173 266
\psfig{figure=wolff.eps}
\end{picture}
\end{center}
\end{minipage} 


\end{slide}


\begin{slide}

\subsection{\bf Depth-First Search (``Ants-in-the-Labyrinth'')}

A Wolff cluster is grown in the same way as the ant colony analogy. 
Placing a baby ant on a new site corresponds to putting a bond between the 
two sites (with probability $p$), and an ant colony corresponds to a cluster.  

This is straightforward 
to program --- just need to keep track of which sites are in the cluster, 
and which are the newly labeled sites (latest generation of ants) 
to be looped over in the next iteration.

This algorithm can also be used to do the connected component labeling for all 
sites in the Swendsen-Wang algorithm, given the configuration of bonds, by
just adding an extra loop over all sites.  If a site has already been labeled 
(assigned to a cluster), go to the next site.  If not, increment the cluster
label counter and give it to a new parent ant on that site, which then 
propagates the label through the labyrinth of bonds connecting sites in the 
cluster.

\end{slide}


\begin{slide}

\subsection{\bf Breadth-First Search (FGHK)} 

First proposed by Fisher and Galler, implemented in a more efficient manner 
by Hoshen and Kopelmann for percolation studies.

Uses the idea of a {\em rooted tree}.  Each site has a pointer to 
another site in the cluster, and so on down a tree-like structure.  The site 
that is the root of the tree has the label for the cluster.  Sites find 
their labels by traversal of the tree. Periodically pruning the tree 
(making shortcuts) makes the algorithm efficient.

The breadth-first search is better suited to a parallelism than the 
depth-first search, since it involves
an outer loop over all sites, which can be parallelized using standard 
domain decomposition.

% Figure ??

\end{slide}


\begin{slide}

\subsection{\bf Local Label Propagation}

The simplest and most obvious parallel component labeling algorithm proceeds 
as follows, assuming one lattice site per (abstract) processor.  
{\parskip=0pt
% \Smallsize
\begin{enumerate}
\item
Assign a unique label to each site (e.g., the processor number). 
\item
Every site looks at its neighbor in the positive $x$ direction. \\ 
If there is a bond between them, and the neighbor's cluster label 
is smaller than the label at the site in question, 
then it sets its label to be the same as its neighbor's.\\ 
If the neighbor's label is bigger, then it gets the smaller label.\\ 
This is repeated in the $y$ direction (and other directions in higher 
dimensions).
\item
Repeat the previous step until none of the labels has changed,
in which case the clusters have been correctly labeled.
\end{enumerate}
}

% Figure??

\end{slide}


\begin{slide}

\subsection{\bf SIMD Local Propagation}

Notice that this algorithm is purely SIMD, since every processor 
does the same thing 
at the same time. Since it propagates the label only one lattice spacing at a 
step, the number of steps required to label the clusters is the maximum path 
length (along bonded sites) between any two sites in the same cluster.  For 
cluster algorithms this will be order $L$ (the length of the lattice), since 
at $T_c$ there is one large cluster.  

This method will work better for 
labeling small connected components, and worse for labeling objects such as 
spirals.  The method is much too inefficient for cluster algorithms, 
basically because the parallel algorithm suffers from the same kind of 
``critical slowing down'' (trying to label a large cluster using local steps) 
as does the sequential Metropolis algorithm!  

Need a \Underline{non-local} parallel algorithm.

\end{slide}


\begin{slide}

\subsection{\bf MIMD Local Propagation}

On a coarse grain MIMD machine, this algorithm is much more effective.  
In this case, 
we can label the clusters on the sub-lattice on each processor using a fast, 
efficient sequential algorithm.  The local label propagation is now done only 
on the \Underline{edges} of the sub-lattices, 
to match up labels of clusters that cross processor boundaries.

Now have good load balancing, and the label of the largest cluster is 
propagated across the number of processors, not the number of sites. 
For an $L^2$ lattice on a $P^2$ processor array:
%
\begin{displaymath}
\mbox{communication} \sim P\left( \frac{L}{P} \right) = L \qquad
\mbox{calculation}   \sim \frac{L^2}{P^2} 
\end{displaymath}
%
\begin{displaymath}
\mbox{comm/calc} \sim \frac{P^2}{L} 
\end{displaymath}
 
As long as $L \gg P$, will get good efficiencies for this method.

\end{slide}


\begin{slide}

\begin{figure}[h]
\begin{center}
\begin{picture}(250,250)
\psfig{figure=self-label.ps}
\end{picture}
\end{center}
\vspace{1in}
\centerline{\Smallsize
Efficiency of MIMD local propagation algorithm on the nCUBE/2.
}
\end{figure}


\end{slide}


\begin{slide}

\subsection{\bf Parallel Labeling for Image Processing}

There is a large literature of parallel algorithms for pixel connectivity 
and connected component labeling applications in image processing. 
But many of these are useless for spin models:
{\parskip=0pt
\begin{enumerate}
\item Pixel connectivity is a special case of component labeling, and some 
of the algorithms cannot be generalized.

\item Most computer science analysis of algorithms concentrates on 
{\it worse case} complexity.  We are interested only in 
{\it average} time complexity to label a large number of 
configurations.

\item Algorithms with optimal asymptotic efficiencies are 
\Underline{not} necessarily optimal on real machines.  
For one ``asymptotically optimal'' parallel component labeling algorithm, 
need $\sim 10^6$ processors to get a speedup of 1!

\item Some algorithms may work well on small, regular objects,
but not on large, irregular spin model clusters.
\end{enumerate}
}

\end{slide}


\begin{slide}

\subsection{\bf Global Equivalencing}

Cluster labeling can be looked on as a special case of an equivalence 
problem\@: given a list of equivalences (site $i \equiv$ site $j$), put all 
the elements into equivalence classes (cluster labels).  There are many 
algorithms known for solving this problem.

Parallel implementation:
{\parskip=0pt
\begin{enumerate}
\item Label sub-clusters sequentially for sublattice on each node.

\item Each processor looks at the edges of its neighbors and creates a list 
of equivalences.

\item Lists are sent to one of the nodes, 
which creates the global equivalence classes 
and returns the results to all nodes.
\end{enumerate}
}

\end{slide}


\begin{slide}

\subsection{\bf Hierarchical Method}

Problems with global equivalencing:
{\parskip=0pt
\begin{itemize}
\item Equivalence step is a sequential bottleneck
\item One node needs enough memory for all cluster labels
\end{itemize}
}
These problems are alleviated by hierarchically merging the labels,
using a quad-tree approach as shown below.
Takes $\log P$ steps rather than $P$.  More complicated to implement,
a number of optimizations are required to get good performance.

\vspace{2mm}

\begin{minipage}[h]{2.5in}
\begin{center}
\begin{picture}(148,148)
%%BoundingBox: 0 0 148 148
\psfig{figure=block64.eps,height=148pt,width=148pt}
\end{picture}
\end{center}
\end{minipage} \ \
%
\begin{minipage}[h]{2.5in}
\begin{center}
\begin{picture}(148,148)
%%BoundingBox: 0 0 219 219
\psfig{figure=block16.eps,height=148pt,width=148pt}
\end{picture}
\end{center}
\end{minipage} \ \
%
\begin{minipage}[h]{2.5in}
\begin{center}
\begin{picture}(148,148)
%%BoundingBox: 0 0 216 216
\psfig{figure=block4.eps,height=148pt,width=148pt}
\end{picture}
\end{center}
\end{minipage} \ \

\end{slide}


\begin{slide}

\subsection{\bf Multigrid Algorithm}

Instead of just looking at neighboring sites, look at sites a distance $2^m$ 
away (this is fast on a hypercube).

If the label is the same at any point, put a connection between the sites.  If 
there are connections $A \stackrel{2^m}{\longrightarrow} B 
\stackrel{2^m}{\longrightarrow} C$, make a connection $A 
\stackrel{2^{m+1}}{\longrightarrow} C$.

Now changes in the labels can be propagated long distances in one iteration.  
Expect labeling in $\sim \log_2 L$ iterations, each of which takes $\sim 
\log_2 L$ multigrid steps, each taking time $O(1)$ on hypercube.  So time 
complexity of $O(\log L)^2$.

% Figure??

\end{slide}


\begin{slide}

\subsection{\bf Get/Send Algorithm (Shiloach-Vishkin)}

Another non-local SIMD algorithm that has been implemented is a method 
similar to an algorithm due to Shiloach and Vishkin. 
It is based on the rooted tree 
approach of the FGHK algorithm (breadth-first search).  

This ``Get/Send'' algorithm again has $\sim \log L$ iterations 
(rather than $L$), and $O(\log L)^2$ time complexity, since non-local 
communications take time $O(\log L)$.

This algorithm is faster in practice than multigrid, 
both of which are much faster than local label propagation.
Currently these are only implemented on SIMD machines, where performance 
is quite poor due to load balance problems. 
Should be much better on MIMD machines.

\end{slide}


\begin{slide}

\subsection{\bf Get/Send Algorithm (cont.)}

Start with every site (abstract processor) having the processor number
as its unique label (i.e., the root of its own tree). 
Each iteration consists of the following steps:
%
{\parskip=0pt
\begin{enumerate}
\item {\bf Local label propagation}

\item {\bf Send} ---
if a processor's label changes during step (1), it sends the new label to the 
root of its current sub-cluster 
(i.e., the processor number given by its old label),  
which picks the smallest of all the labels it receives as its new label.

\item {\bf Scan} (optional) ---
propagate labels fast along rows and columns, using {\em scan} or
{\em parallel prefix} operations.

\item {\bf Get} ---
every site gets the (possibly new) label from the root processor for its 
current label.  

\item {\bf Check for completion} --- continue if any labels have changed.
\end{enumerate}
}

\end{slide}


\begin{slide}

\centerline{\psfig{figure=SIMD-iterations.eps,height=5in}}
\vspace{-1cm}
{\Smallsize 
The Get/Send and MultiGrid algorithms are logarithmic
in the lattice size $L$, unlike the simple local label propagation 
algorithm which usually takes time of order $L$, but can be order $L^2$
for complex clusters like spirals.
}

\end{slide}



\begin{slide}

\begin{figure}[h]
\begin{center}
\begin{picture}(380,380)
%%BoundingBox: 116.000000 221.120003 496.000000 601.120003
\psfig{figure=new1.ps}
\end{picture}
\end{center}
{\Smallsize
1. Get/Send algorithm after local label propagation and a scan operation. 
The colors represent different cluster labels, and the algorithm starts
with random labels at each site.
}
\end{figure}

\end{slide}


\begin{slide}

\begin{center}
\begin{picture}(380,380)
\psfig{figure=new2.ps}
\end{picture}

{\Smallsize
2. Get/Send algorithm after a send.
}
\end{center}

\end{slide}


\begin{slide}

\begin{center}
\begin{picture}(380,380)
%%BoundingBox: 116.000000 221.120003 496.000000 601.120003
\psfig{figure=new3.ps}
\end{picture}
\end{center}
\centerline{\Smallsize 3. Get/Send algorithm after a few iterations.}


\end{slide}


\begin{slide}

\begin{center}
\begin{picture}(380,380)
%%BoundingBox: 116.000000 221.120003 496.000000 601.120003
\psfig{figure=new11.ps}
\end{picture}
\end{center}
\centerline{\Smallsize 4. The final cluster labels.}

\end{slide}


\begin{slide}

\subsection{\bf Parallel Wolff Algorithm}

To parallelize the Wolff algorithm, need to use a parallel 
depth-first search (``ants-in-the-labyrinth'') algorithm, 
since that is how the Wolff cluster is grown. 
Unfortunately, this is a sequential process --- each ant generation has 
to follow the previous one.  

The only parallelism is \Underline{within} 
a generation, i.e. over the ``wavefront'' of the expanding cluster.  
The amount of parallelism is very limited, especially if the cluster being 
grown is small. 

A major problem in parallelizing depth-first search methods
is that the load balance can be very poor, since the computation is usually
confined to a particular (contiguous) section of the data.
Can sometimes alleviate this problem by using {\em scattered} or {\em cyclic}
data distribution, rather than the standard block distribution.
Howver there is a tradeoff --- the data is no longer local so 
extra communication is required.

\end{slide}


\begin{slide}

\subsection{\bf Problems with Parallel Ants}

% \begin{figure}[h]
\begin{center}
\begin{picture}(256,256)
%%BoundingBox: 0 0 432 432
\psfig{figure=block-wolff.eps,height=256pt,width=256pt}
\end{picture}
\end{center}
% \end{figure}
{\Smallsize
The Wolff cluster is shown here in gray, with the initial site in black.
A particular generation of ants is shown in red. The black squares 
represent processor boundaries for a standard {\tt (BLOCK,BLOCK)}
data distribution.
Note that the load balance (number of red sites per processor)
is very poor.
}
\end{slide}


\begin{slide}
\subsection{\bf Parallel Ants With Cyclic Data Distribution}

\begin{minipage}[h]{3.9in}
\begin{center}
\begin{picture}(240,240)
%%BoundingBox: 0 0 432 432
\psfig{figure=column-cyclic.eps,height=240pt,width=240pt}
\end{picture}
\end{center}
{\Smallsize 
In a \mbox{column-cyclic} \nolinebreak[4] distribution \linebreak[4]
\mbox{\tt (*,CYCLIC(2))}
data is dealt out to processors (indicated with different colors)
like cards to card players, 2 columns at a time. 
}
\end{minipage} \ \quad \
%
\begin{minipage}[h]{3.9in}
\begin{center}
\begin{picture}(240,240)
%%BoundingBox: 0 0 432 432
\psfig{figure=cyclic-wolff.eps,height=240pt,width=240pt}
\end{picture}
\end{center}
{\Smallsize
The Wolff cluster shown in the previous figure but now with a  
column-cyclic distribution of data to processors. 
Note that the load balance is much improved by this non-local 
data distribution.
}
\end{minipage} 
\end{slide}


% Compare block and cyclic performance?? (do we do this in paper??)

% Performance plot ??

\begin{slide}

\subsection{\bf Independent Parallelism}

For a data parallel Wolff algorithm, it is not possible to get the 
high efficiences that can be obtained for regular (e.g. Metropolis)
algorithms, even using a cyclic data distribution. 
However highly irregular Monte Carlo methods, such as the Wolff algorithm,
{\em can} be parallelized efficiently, since domain decomposition is 
not the only level of parallelism available for MC algorithms.

Remember that the basic idea of MC simulation is to 
approximate an extremely large or infinite sum 
(or large or infinite dimensional integral) by 
a finite sum over important configurations. The error is predominantly 
statistical, that is, proportional to $1/\sqrt{N}$, where $N$ is the number of 
independent configurations. So here is another way the simulation can be done 
in parallel --- if different processors are generating independent 
configurations to add to the sum.

\end{slide}


\begin{slide}

\subsection{\bf Independent Monte Carlo Simulations}

This is very easy to do.  Just set up different random number streams on 
each processor (which is done in any case on a parallel machine, but here 
different processors could be different workstations in a network), set up 
different random initial configurations on each processor, and then run 
completely independent simulations on each processor.  This gives a
parallel efficiency of 100\%!

This simple idea works very well for cluster algorithms, as long as the 
lattice is small enough that all the data will fit into the memory of a single 
processor. This was a problem for early parallel machines, but is usually no
problem for workstation networks and machines like the IBM SP2. 
Note that this only works for MIMD machines, since 
cluster algorithms are irregular.

\end{slide}


\begin{slide}

\subsection{\bf ``Hybrid'' Parallelism}

Parallel Swendsen-Wang cluster algorithms for moderate sized lattices 
$(\sim 512^2)$ work fairly efficiently for moderate numbers of nodes $(\sim 
64)$ on a MIMD machine.  However, this does not prevent us from using 
massively parallel MIMD machines with $\sim 1000$ nodes for this problem.  

Can use a hybrid of domain (data) parallelism and independent parallelism 
(over different MC runs). For example, can run 8 independent 
simulations of 64 nodes each on a 512-node machine.

This gives much better performance than running parallel cluster 
algorithms on vector supercomputers, since these applications do not
vectorize well.

\end{slide}


\begin{slide}
\subsection{\bf Message Passing vs Data Parallel}

Coarse grain parallel algorithms using message passing on MIMD machines
give excellent performance for regular algorithms such as Metropolis, 
and fairly good performance for irregular cluster algorithms. This is because 
the efficient sequential algorithms can be used on the sub-domain on 
each processor.

Data parallel languages such as High Performance Fortran handle regular
algorithms very efficiently, and are 
much easier to program that explicit message passing languages.
However they do not perform as well for irregular cluster algorithms.
This is mainly because the slower SIMD algorithms must be used
even for the domain within each processor. There are also load balancing 
problems. 

\end{slide}


\begin{slide}
\subsection{\bf Irregular Problems in HPF}

These types of irregular, non-local problems provide a real challenge to
efficient implementation (for algorithms and compilers) in data parallel 
languages such as HPF. 

It may be that good performance
can only come by using library routines, e.g. for connected component
labeling in this case, which would be implemented using an efficient 
message passing routine.

HPF can handle independent and hybrid parallelism using the 
{\tt INDEPENDENT} construct.

\end{slide}


