
\begin{slide}

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

\end{slide}


\section{\bf Parallel Metropolis Algorithms}


\begin{slide}

\subsection{\bf 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
\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 a site and its neighbors at the same time, 
since they are dependent, and this would violate detailed balance.  % WHY??

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.

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

\end{slide}


\begin{slide}

\subsection{\bf Illustration of Parallel Metropolis Algorithm}

Figure??

\end{slide}


\begin{slide}

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

\end{slide}


\begin{slide}

\subsection{\bf Metropolis Code in a Data Parallel Language}


\end{slide}


\begin{slide}

\subsection{\bf Parallel Metropolis -- Message Passing}

If the number of lattice sites per processor is $>2^d$,  
need not bother with the checkboard scheme, 
since there is no longer the 
possibility of updating neighboring sites at once.

Notice that this algorithm, when written using a message passing language
such as Express or MPI, looks $\Underline{\rm 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.

\end{slide}


\begin{slide}

\subsection{\bf Metropolis Code in a Message Passing Language}


\end{slide}


\begin{slide}

\subsection{\bf Efficiency of the Parallel Algorithm}

Since the Metropolis algorithm for spin models is local and regular,
it parallelizes very efficiently 
(even for the Ising model, with very little computation), 
$\Underline{\rm provided}$ there is more than one site per processor. 
Otherwise the efficiency is at most 50\% since half the processors are idle
for the checkerboard update.  

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)$.  

Communication time $\propto l$ 
(\# edge sites of sub-lattice, i.e. perimeter).

Calculation time $\propto l^2$ 
(\# sites of sub-lattice, i.e. volume).  

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

\end{slide}


\begin{slide}

\subsection{\bf Coarse Grained Algorithm}

On a coarse grained message passing machine, there is a more efficient 
parallel Metropolis algorithm.  

This requires a checkerboard partitioning of the lattice. 
But now the data communication is done as large blocks of spins 
rather than single spins, by passing all the edge data to neighboring 
processors after each red or black update.  

This {\em blocked communication} reduces the latency time for 
communication enormously, and greatly improves efficiency.

It is also possible to overlap communications with computation, 
by computing edge values first, and then sending them while the
interior points are being computed.

\end{slide}


\begin{slide}

\subsection{\bf 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 energy = SUM (energy-per-site) / volume}.

\end{slide}


\begin{slide}

\subsection{\bf Non-local Measurements}

Some measurements are $\Underline{\rm non-local}$, e.g., the correlation 
function $\Gamma (n) = \sum_i S_i S_{i+n}$.
This can also be easily calculated in parallel.  
{\parskip=0pt
\begin{enumerate}
\item
Take a copy of the spin variables and keep shifting them one site 
at a time. 
\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. $\Gamma (n) = \mbox{SUM}\ \left( S_i S_{i+1} \right)$.
\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 in turn, in order to reduce latency. 

\end{slide}


\begin{slide}

\subsection{\bf 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 parallelizeed just like coarse grain parallel Metroplis algorithm.
{\parskip=0pt
\begin{itemize}
\item
standard {\tt (BLOCK,BLOCK)} data decomposition
\item
red/black or checkerboard update
\item
local (blocked) communications 
\end{itemize}
}

\end{slide}


\section{\bf Parallel Cluster Algorithms}


\begin{slide}

\subsection{\bf Parallel Cluster Algorithms}

For the Swendsen-Wang cluster algorithm, the standard 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 direction, and put a bond there with 
probability $\left(1-d^{-\beta}\right)\,\delta \left( s_i, s_j \right)$ for 
the Potts model.  Note that $\beta=\frac{J}{kT}$.

The only difficult problem in parallel is the $\Underline{\rm connected\ 
comonent\ labeling}$, that is, assigning to all the sites a cluster label that 
is unique for each cluster.

\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 computational 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 configuratio is of 
order 1000, whereas for the cluster algorithms it is of order 10, that is, a 
100 fold improvement in speed.  

Howver, we have seen that 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 the improvement of the 
cluster algorithm is lost unless we can find an efficient parallel algorithm.

\end{slide}


\begin{slide}

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

This type of problem occurs in many different guises in computational physics.  
In many applications, there are standard algorithms that are regular, simple, 
possibly 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, 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 DEs and spins 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}

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:

Pick a site at random and place an ant on it.  The ant reproduces by placing a 
child on the four neighboring sites (in two dimensions) with probability $p$.  
This ``first generation'' of ants also looks at each of its four neighbors in 
turn, and if it is not already accupied, places a child of its own there with 
probability $p$.  The ant colony continues to grow unti the final generation 
produces no children.

\end{slide}


\begin{slide}

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

Clearly, we grow a Wolff cluster in exactly the same fashion.  Placing a baby 
ant on a new site corresponds to putting a bond between the two sites (with 
probability $p$), and the ant colony is the cluster.  This is straightforward 
to program---we merely 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.  We 
just add an extra loop over all sites.  If a site has already been labeld 
(assigned to a cluster), we go to the next site.  If not, we increment the 
label counter and give it to our ant, which propagates it through the 
labyrinth of bonds connecting sites in the cluster.

\end{slide}


\begin{slide}

\subsection{\bf Fisher-Galler-Hoshen-Kopelmann (Breadth-First Search)}

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.

\end{slide}


\begin{slide}

\subsection{\bf Parallel Ants}

If we want to parallelize the Wolff algorithm, we need to use a parallel 
``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{\rm within}$ 
a generation, that is, over the ``wavefront'' of the expanding cluster.  There 
is some parallelism, but it is very sparse, especially if the cluster being 
grown is very small (which is often the case for the Wolff algorithm).  So the 
Wolff algorithm is very inefficient on a vector or parallel machine.

Clearly, this is not what we want to do for a parallel Swendsen-Wang algorithm, 
either.

\end{slide}


\begin{slide}

\subsection{\bf Local Label Propagation}

The simplest and most obvious parallel component labeling algorithm proceeds 
as follows.  Suppose we have one site per processor.  Assign a unique label to 
each site (e.g., the processor number).  Every site looks at its neighbor in 
the positive $x$ direction.  If they have 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 receives the smaller label.  This is repreated in the 
$y$ direction.

\end{slide}


\begin{slide}

\subsection{\bf SIMD Local Propagation}

Notice that this algorithm is purely SIMD, since everyone 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.  This 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{\rm 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.  Here, 
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 edges of the sub-lattices, to match up labels of clusters that cross 
processor boundaries.

Now we have good load balancing, and label of largest cluster is propagated 
across number of processors, not number of sites.  For an $L^2$ lattice on a 
$p^2$ processor array:
\begin{eqnarray*}
\mbox{communication} & \sim & p\left( \frac{L}{p} \right) = L \\
\mbox{calculation}   & \sim & \frac{L^2}{p^2} \\
\end{eqnarray*}

As long as $L>>p$, will get good efficiencies for this method.
\begin{displaymath}
\mbox{comm/calc} \sim \frac{p^2}{L} \hspace{1em} \mbox{\# proc} << L
\end{displaymath}

\end{slide}


\begin{slide}

\subsection{\bf Parallel Labeling Algorithms 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 our application.
\begin{enumerate}
\item Pixel connectivity is a special case of component labeling, and some 
algorithms cannot be generalized.

\item Most computer science analsysis of algorithms concnetrates on 
$\Underline{\rm worse\ case}$ complexity.  We are interested only in 
$\Underline{\rm average}$ time complexity to label a large number of 
configurations.

\item Optimal asymptotic efficiencies are $\Underline{\rm not}$ necessarily 
optimal on real machines.  May need $>> 10^6$ processors!

\item 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:
\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}

\end{slide}


\begin{slide}

\subsection{\bf Hierarchical Method}

Problems with global equivalencing:
\begin{enumerate}
\item Equivalence step sequential bottleneck
\item One node needs enough memory for all cluster labels
\end{enumerate}

Solved (partially) by hierarchically merging the labels.

Same algorithm used by Embrechts, et al., for component labeling problem. Take 
$\log P$ steps rather than $P$.  Complicated to implement---could be improved.

\end{slide}


\begin{slide}

\subsection{\bf ``Multi-grid'' Algorithm (SIMD)}

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}{\rightarrow} B 
\stackrel{2^m}{\rightarrow} C$, make a connection $A 
\stackrel{2^{m+1}}{\rightarrow} 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$ multi-grid steps, each taking time $o(1)$ on hypercube.  So time 
complexity of $O(\log L)^2$.

\end{slide}


\begin{slide}

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

Another non-local SIMD algorithm that we implemented was 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).  

Again, we start with every processor having a unique label (i.e., the root of 
its own tree).  Each iteration consists of the following steps:
\begin{enumerate}
\item Local label propagation

\item Send

If a site's label changes during step (1), it sends the new label to the 
processor labeled by this site's old label (i.e., the root of the current sub-
cluster), which picks the smallest of all the labels it receives as its new 
label.

\item Scan

Propagate labels fast along rows and columns.

\item Get

Every site gets the (possibly new) label from the root processor for its 
current label.  The get/send algorithm again has $\sim \log L$ iterations 
(rather than $L$), and $\nu (\log L)^2$ time complexity since non-local 
communications take time $O(\log L)$.

Get/send faster in practice than multigrid, both much faster than local label 
propagation.

Only implemented on SIMD machines, where performance is poor due to usual load 
balance problems.  Again, should be much better on MIMD machines (not 
implemented yet).
\end{enumerate}

\end{slide}


\begin{slide}

\subsection{\bf Independent Parallelism}

All is not lost for the parallelization of highly irregular Monte Carlo 
methods, such as the Wolff algorithm.  This is because domain decomposition is 
not the only level of parallelism available for MC algorithms.

Remember that the basic idea of MC simulation is that we are approximating 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 MC Simulations}

This is very easy to do.  We just set up different random number streams on 
each processor (which we would do 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.  So now we get 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, which may be a problem for parallel machines, but not too bad for 
workstation networks.  Note that this only works for MIMD machines, since 
cluster algorithm is irregular.

\end{slide}


\begin{slide}

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

We have seen that parallel S-W 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.  We 
can use a hybrid of domain (data) parallelism with independent parallelism 
(over different MC runs) by, for example, running eight independent 
simulations of 64 nodes each on a 512-node machine.

We did this for the old nCUBE-1 at Caltech, and got performances equivalent to 
an eight-processor CRAY~Y-MP for this problem.

\end{slide}

