


\subsection{ Parallel Monte Carlo Algorithms}

{\bf Metropolis} algorithm for spin models is easy to parallelize efficiently,
since it is regular and local.  
{\parskip=0pt
\begin{itemize}
\item
Regular data structures $\Rightarrow$ simple domain decomposition 
\item
Regular algorithm $\Rightarrow$ lots of parallelism and perfect load balance 
\item
Local interactions $\Rightarrow$ local communications  
\end{itemize}
}
\vspace*{1mm}

{\bf Swendsen-Wang} algorithm is difficult to parallelize efficiently,
since it is irregular and non-local.  
{\parskip=0pt
\begin{itemize}
\item
Clusters are irregular (fractal!) in size and shape $\Rightarrow$ 
      hard to load balance, lots of non-local communication
\end{itemize}
}
\vspace*{1mm}

{\bf Wolff} algorithm is almost impossible to parallelize efficiently. 
{\parskip=0pt
\begin{itemize}
\item
Grow a single, irregular cluster starting at a random lattice site
      $\Rightarrow$ little parallelism, very hard to load balance 
\end{itemize}
}




\section{Parallel Metropolis Algorithms}




\subsection{ Parallel Metropolis Algorithm}

The Metropolis algorithm for a spin model is well suited to parallelism, 
since it is 
{\parskip=0pt
\begin{itemize}
\item {\bf regular} --- so can use standard data decomposition of the lattice 
      onto the processors and get good load balance
\item {\bf local} --- the update of a site depends only on its nearest 
neighbors, so only local communication is required
\end{itemize}
}

However, cannot update all sites at the same time.
This would violate detailed balance, 
since a site and its nearest neighbors are dependent, 
and thus must be updated separately.
% EXPLAIN ?? 






\subsection{ Red/Black or Checkerboard Update}

Instead, the lattice is partitioned into a checkerboard of alternating 
red and black sites (note this can also be done for $d>2$ dimensions).
Can then update all the black sites in parallel (since they are independent),
followed by all the red sites in parallel.

\begin{figure}[h]
% \vspace*{0.10in}
\begin{center}
\begin{picture}(144,144)
%%BoundingBox: 0 0 144 144
\psfig{figure=checkerboard.eps}
\end{picture}
\end{center}
\end{figure}

This is known as a {\em red/black} or {\em checkerboard} update.

Note that unless there is more than one site per processor,  
the efficiency is at most 50\% since half the processors will be idle.







\subsection{ Parallel Metropolis -- Data Parallel}

Parallel Metropolis using a checkerboard update is a data parallel
SIMD algorithm.  
It is very simple to code in Fortran 90 and High Performance Fortran (HPF).
{\parskip=0pt
\begin{itemize}
\item 
Use standard domain decomposition and distribute the 2-$d$ lattice 
among processors using {\tt (BLOCK,BLOCK)} form.
\item 
Use a {\em logical mask} ({\tt black}) which only activates 
the appropriate (abstract) processors which deal with black sites. 
Then change the mask {\tt (black = .NOT.black)} and update the other sites.
\item 
To get data from neighboring sites, 
use periodic shift operations {\tt (CSHIFT)} for 
periodic boundary conditions.
\end{itemize}
}






\subsection{ Metropolis Code in a Data Parallel Language}







\subsection{ Parallel Metropolis -- Message Passing}

\begin{minipage}[h]{6in}
If the number of lattice sites per processor is $>2^d$,  
need not use the checkerboard scheme, since cannot 
update neighboring sites at once.
\end{minipage} \ \
%
\begin{minipage}[h]{2.0in}
\begin{center}
\begin{picture}(108,73)
%%BoundingBox: 0 0 108 73
\psfig{figure=coarse-nbors.eps}
\end{picture}
\end{center}
\end{minipage} 

Notice that this algorithm, when written using a message passing language
such as Express or MPI, looks {\bf exactly} the same as 
the sequential algorithm!
%
The only difference is how the {\tt shift} subroutine, which finds 
neighboring spin values, is implemented:
{\parskip=0pt
\begin{itemize}
\item
sequential algorithm --- {\tt shift} handles the boundary conditions.
\item
parallel algorithm --- {\tt shift} handles boundary conditions and does 
message passing if the neighboring spin is on a different processor.
\end{itemize}
}
Also, for the parallel algorithm, the {\tt size} variable refers to
the size of the sub-lattice on each processor, not the total lattice size.






\subsection{ Metropolis Code in a Message Passing Language}







\subsection{ Efficiency of the Parallel Algorithm}

Since the Metropolis algorithm for spin models is local and regular,
it should parallelize very efficiently, even for the Ising model, 
which has very little computation. 

For a 2-$d$ grid of $N=P \! \times \! P$ 
processors, use a {\tt (BLOCK,BLOCK)} distribution of the 
$V=L \! \times \! L$ sites of a 2-$d$ lattice over the processor grid so 
that every processor has an $l \! \times \! l$ sub-lattice $(l=L/P)$.  

\begin{minipage}[h]{3.4in}
{\parskip=5pt
Communication time $\propto l$ \\
(\# edge sites of sub-lattice, \\
i.e. perimeter). 
  
Calculation time $\propto l^2$ \\
(\# sites of sub-lattice, \\
i.e. volume).  
}
\end{minipage} \ \
%
\begin{minipage}[h]{4.0in}
\begin{center}
\begin{picture}(140,138)
%%BoundingBox: 0 0 140 138
\psfig{figure=metrop-mesh.eps}
\end{picture}
\end{center}
\end{minipage} 

Thus, as long as $l$ is large enough, the communication/calculation ratio 
will be small, and the efficiency will be near 1.






\subsection{ Coarse Grained Algorithm}

On a coarse grained MIMD machine, there is a more efficient 
parallel Metropolis algorithm.  
%
This requires a checkerboard partitioning of the lattice. 
The data communication is done as large blocks rather than single values,
by passing a block of edge data to neighboring 
processors after each red or black update.  

This {\em blocked communication} greatly reduces the latency time for 
communication, thereby improving efficiency.
%
It is also possible to overlap communications with computation, 
by computing edge values first, and sending them while the
interior points are being computed.

\vspace*{-5mm}

\begin{center}
\begin{picture}(210,115)
%%BoundingBox: 0 0 210 115
\psfig{figure=metrop-coarse.eps}
\end{picture}
\end{center}



% FIGURE metrop-effic.ps  ??



\subsection{ Measurements}

Measurements of standard quantities such as energy and magnetization are 
easily done in parallel.  
{\parskip=0pt
\begin{enumerate}
\item
Calculate the quantities locally (e.g., energy per site), 
which can be done completely in parallel (no checkerboarding required)
\item
Do a partial sum over sites in each processor 
\item
Combine the partial sums on each processor, using for example 
% {\tt MPI_COMBINE} in MPI, ??
and {\tt EXCOMBINE} in Express C or {\tt KXCOMB} in Express Fortran.
\end{enumerate}
}
In data parallel languages such as HPF, the last two steps can be
implemented using the {\tt SUM} primitive. For example, \\
{\tt \qquad energy = SUM (energy-per-site) / volume}






\subsection{ Non-local Measurements}

Some measurements are {\bf non-local}, e.g., the correlation 
function $\Gamma (n) = \sum_i S_i S_{i+n}$.
This can also be easily calculated in parallel, by continually shifting
a copy of the spin variables one site at a time.  
{\parskip=0pt
\begin{enumerate}
\item
After the first shift, every (abstract) processor calculates the local value 
of $\Gamma (1)$, i.e., $S_i S_{i+1}$, at each site, and the sum over sites is 
taken (e.g. using {\tt SUM}).
\item
After the second shift every processor calculates the local value of 
$\Gamma (2)$, and the sum over sites is done.
\item
This is continued until the spin has been shifted halfway around the lattice 
(only halfway due to periodic boundary conditions).
\end{enumerate}
}
For message passing languages it is more efficient to only do the partial 
sums on each processor after every shift.
The global sums are done at the end on the vector $\Gamma (n)$ rather than
each of the individual values $\Gamma (1), \Gamma (2)$, \ldots in turn, 
so as to reduce latency. 






\subsection{ Relation to PDE Solvers}

The Metropolis algorithm for a 2-$d$ spin model is similar in many ways
to numerical methods for solving differential equations, 
such as Laplace's equation $\nabla^2 \phi = 0$. 
This can be discretized onto a 2-$d$ grid, with the update depending only 
on nearest neighbor points, e.g., an iterative scheme to solve this equation 
would replace $\phi_i$ by the average value of its four neighbors.  

It is possible to iteratively update all sites at once (Jacobi algorithm). 
However, using a red/black update gives better convergence (like Gauss-Seidel).
This is parallelized just like coarse grain parallel Metropolis algorithm.
{\parskip=0pt
\begin{itemize}
\item
standard block data decomposition
\item
red/black or checkerboard update
\item
local (blocked) communications 
\end{itemize}
}




\section{Parallel Cluster Algorithms}




\subsection{ 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.






\subsection{ 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.






\subsection{ Early Parallel Algorithms}

Burkitt-Heermann efficiencies






\subsection{ 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!






\subsection{ Labeling or Growing a Single Cluster}

Let us first look at how sequential component labeling algorithms work.  
We will start with the Wolff algorithm, which ``grows'' a single cluster.  
Amusingly enough, we can liken this to the growth of an ant colony, in the 
following manner:
{\parskip=0pt
\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 ``first 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}
}






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

A Wolff cluster is grown in the same fashion 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 (baby ants) that need 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 ant, which propagates it through the 
labyrinth of bonds connecting sites in the cluster.

% Figure ??






\subsection{ 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, 
which is the root of the tree, has the label for the cluster.  Sites find 
their labels by traversal of the tree. Prune the tree during the algorithm 
(make shortcuts) to make 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.






\subsection{ Local Label Propagation}

The simplest and most obvious parallel component labeling algorithm proceeds 
as follows, assuming one lattice site per (abstract) processor.  
{\parskip=0pt
\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}
}






\subsection{ 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 {\bf non-local} parallel algorithm.






\subsection{ 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 {\bf 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.






\subsection{ Parallel Labeling for Image Processing}

There is a very large literature of parallel algorithms for pixel connectivity 
and connected component labeling applications in image processing. 
But many of these are useless for our application:
{\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 
{\bf worse case} complexity.  We are interested only in 
{\bf average} time complexity to label a large number of 
configurations.

\item Algorithms with optimal asymptotic efficiencies are 
{\bf 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}
}






\subsection{ 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 clusters sequentially on each node.

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

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






\subsection{ 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 (mostly) solved by hierarchically merging the labels.
Takes $\log P$ steps rather than $P$.  More complicated to implement,
a number of optimizations are required to get good performance.






\subsection{ 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$.






\subsection{ 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.






\subsection{ 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}
}






\subsection{ Parallel Wolff Algorithm}

To parallelize the Wolff algorithm, need to use a parallel 
``ants-in-the-labyrinth'' (depth-first search) 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 {\bf within} 
a generation, i.e. over the ``wavefront'' of the expanding cluster.  
This parallelism is very sparse, especially if the cluster being 
grown is small (often the case for the Wolff algorithm). 
So a naive parallel Wolff algorithm is very inefficient on a vector or 
parallel machine.

Another major problem in parallelizing depth-first search algorithms
is that the load balance is 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.






\subsection{ Parallel Ants}

Figures??






\subsection{ Independent Parallelism}

In spite of what we have just said,  
it {\em is} actually possible to efficiently parallelize highly 
irregular Monte Carlo methods, such as the Wolff algorithm. 
This is because domain decomposition is 
not the only level of parallelism available for Monte Carlo 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.






\subsection{ 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 not
a worry for workstation networks and machines like the IBM SP2. 
Note that this only works for MIMD machines, since 
cluster algorithms are irregular.






\subsection{ ``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.




