
% Better explanation of why parallel updated violate detailed balance 

% Figures for parallel Metropolis Code in message passing and parallel.
%   Use HPF instead of CMF.

% {\tt MPI_COMBINE} in MPI, ??  instead of using Express as an example


\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{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 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 ?? 

\end{slide}


\begin{slide}

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


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

% Figure ??

\end{slide}


\begin{slide}

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

% Figure ??

\end{slide}


\begin{slide}

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

\end{slide}


\begin{slide}

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

\end{slide}


\begin{slide}

\begin{figure}[h]
\begin{center}
\begin{picture}(250,250)
\psfig{figure=metrop-effic.ps}
\end{picture}
\end{center}
\vspace{1in}
\centerline{\Smallsize
Efficiency of coarse and fine grain 
Metropolis algorithms on 16 nodes of the nCUBE/2.
}
\end{figure}

\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 \qquad energy = SUM (energy-per-site) / volume}

\end{slide}


\begin{slide}

\subsection{\bf Non-local Measurements}

Some measurements are \Underline{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. 

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

\end{slide}


