%\documentstyle[11pt]{article}
\documentstyle[11pt,psfig]{article}


%\input epsf

\begin{titlepage}

\topmargin=-10mm
\baselineskip=18pt


\begin{center}

{\Huge 
\underline{CIS 712  Term Project} \\}

\bigskip
\bigskip
\bigskip
{\LARGE \bf
(Parallel Wolff Cluster Algorithms for Two-Dimension Spin Model)\\}
\bigskip
\bigskip
\bigskip

{\Large \bf
SUNG-HOON KO \hspace{2cm} SUNG-YONG PARK\\  }
{\large \sl(shko@npac.syr.edu) \hspace{3.5cm} (sypark@npac.syr.edu)}


\smallskip
\bigskip
{\large \it
School of Computer Science, Syracuse University\\}
%

\bigskip
\bigskip
{\large December 1994} \\
\end{center}
\bigskip


\vspace{0.5in}

%#########################################################################
\begin{abstract}
%\noindent
Monte Carlo simulation is a very important numerical method for 
studying spin models in physics. Unfortunately, traditional Monte Carlo 
algorithms like Metropolis algorithm suffer from critical slowing down 
near the regions of interest - the critical points which separate different 
phase of the system.

In traditional Monte Carlo algorithms like Metropolis algorithm, update 
is local, that is, update only one site at a time, and this update depends 
only on the spin values at the site and its neighbors. It is simple to 
program, but the performance is poor. To overcome this poor performance 
and accelerate the dynamics of Monte Carlo updates, a number of different 
methods have been attempted. One of these acceleration algorithms is Cluster 
algorithm.

The Wolff single-cluster algorithm is the most efficient method known
for Monte Carlo simulation of many spin models. 
Due to the irregular size, shape and position of the Wolff clusters,
this method does not easily lend itself to efficient parallel
implementation, so that simulations using this method have thus far
been confined to workstations and vector machines.
Here two parallel implementations of this algorithm are presented.
\end{abstract}

\end{titlepage}

\vspace{1.0in}

\textwidth=5.5in
\textheight=8.5in
\topmargin=-0.5in
\oddsidemargin=0.25in

\begin{document}

\baselineskip=17pt

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Monte Carlo simulation is a very important numerical method for 
studying spin models in physics~\cite{MONTECARLO}. However, 
the performance of traditional Monte Carlo algorithms, 
like Metropolis algorithm~\cite{MONTECARLO,METROPOLIS}, is poor. 
One of the acceleration 
algorithms to overcome poor performance and accelerate the 
dynamics of Monte Carlo update is Cluster algorithm
~\cite{SOKALREV,LAT89REV,SWREV}.

Cluster algorithms are non-local Monte Carlo update schemes 
which can greatly increase the efficiency of computer simulation 
of spin models. The most well known sequential algorithms in 
Monte Carlo cluster update algorithms are Swendsen-Wang algorithm 
~\cite{SW} and Wolff algorithm~\cite{WOLFF}.

The Wolff single cluster algorithm is a variation of the Swendsen-Wang 
cluster algorithm for the Ising model. 

In the Swendsen-Wang algorithm for the Ising model, 
bonds are placed with probability ${1-e^{-2\beta}}$ between 
all neighboring sites of the lattice that have the same spin value,
where $\beta$ is the inverse temperature parameter.
This procedure creates clusters of bonded sites of the same spin.
The spins are updated by assigning a random spin value to each
cluster.

In the Wolff algorithm, only a single cluster is created.
A site is chosen at random and a cluster is constructed around it,
using the same bond probabilities as for the Swendsen-Wang algorithm.
All the spins in this cluster are then flipped, that is, collectively
changed to the alternate spin value. 

However, cluster algorithms which are non-local algorithms are more 
difficult to implement efficiently on parallel machines than local 
algorithms such as Metropolis algorithm. Specially, Wolff algorithm is
a classic example of a problem which is virtually impossible to 
parallelize efficiently~\cite{PL}. Only one cluster is grown and the cluster
is irregular in size and shape. These properties make efficient 
parallelization to be difficult. In this project, we have implemented two
different SIMD algorithms, and compared their performance.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{The Wolff Algorithm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The Wolff algorithm~\cite{WOLFF} is a single cluster Monte Carlo algorithm 
as a variation on the Swendsen-Wang cluster algorithm for the Ising model. 
The main difference for the Ising model is that only single cluster is 
formed and flipped with probability 1. The single cluster is grown 
starting from a site selected at random. On the otherhand, in the S-W 
algorithm, all clusters are formed and flipped with probability 1/2. 
But two algorithms use the same bond activation probabilities~\cite{DA,PRW}.
The algorithm is described by the pseudo-code in Fig.~\ref{fig:wolff}.

Since the Wolff algorithm is a variation on the SW algorithm, we can 
make a less efficient but equivalent implementation of Wolff algorithm. 
We first activate all bonds of the lattice with the appropriate 
probabilities, and decompose all sites into clusters. These are the SW 
clusters. Then we pick a random site $x$, and flip all spins within the 
cluster to which $x$ belongs and ignore the remaining ones~\cite{WOLFF}. 
One sweep of Wolff algorithm has a computational cost of the cluster size, 
while the S-W algorithm has the cost of the entire lattice size. 

\begin{figure}[tbh]
\vspace{10mm}
\baselineskip=17pt
\begin{tabbing}
xxx\=xxx\=xxx\=XXXXXX\= \kill

\>  1) Pick a site $i_0$ at random as the first site in a cluster. \\
\>  2) Grow a cluster from a site $i$(initially, $i_0$) in the cluster by \\
\>\>   connecting bonds to nearest neighbor $j$ with probability $p(i,j)$\\
\>\>\>  where \> $p(i,j)={1-e^{-\beta(1+S_iS_j)}}$, \\
\>\>\>\>         {\it $S_i,S_j$ = spin value,}      \\
\>\>\>\>         {\it $\beta$ = J/KT,}             \\
\>\>\>\>         {\it J =  interaction strength between two spins,}  \\
\>\>\>\>         {\it K =  Boltzmann's constant,}   \\
\>\>\>\>         {\it T =  temperature}            \\
\>  3) Flip the spins in the cluster just formed. \\
\>  4) Continue step 2 and 3 iteratively until no new sites are generated.\\

\end{tabbing}
%\vspace*{2cm}
\protect\caption { 
\protect\label {fig:wolff}
The Wolff algorithm
}
\vspace{1cm}
\end{figure}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Sequential Ants-in-the-labyrinth}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In the Wolff algorithm, the key point is how to grow a single cluster.
One of the simple methods is an {\it Ants-in-the-labyrinth} algorithm which 
is the exactly same as the breadth-first-search in an undirected graph
~\cite{CP,PL}. 

First we pick a site at random and place an ant on it. The ant reproduces 
by placing a child on the four neighboring sites with probability $p$. 
This first generation of ants also looks at each of 4 neighbors in turn, 
and if it is not already occupied, places a child of its own there with 
probability $p$. The ant colony continues to grow until the final generation
produces no children~\cite{PL}. We can consider a given configuration as an 
undirected graph, so that this ants algorithm is the same as a 
breadth-first-search. A breadth-first-traversal of a cluster associates 
a spanning tree to the cluster. The randomly selected site at first 
becomes a root of the spanning tree, and each generation of ants 
corresponds to the each level of the tree.
The algorithm is described by the pseudo-code in Fig.~\ref{fig:SeqAnts}.



\subsection{Algorithm}
For the sequential Ants algorithm, we need a type queue that allows two 
operations enque and deque. The following pseudo code is exactly same as
a breadth-first-search. The only difference is that there is no global 
adjacent matrix(or list). When a site is dequed from a queue, the bond 
is generated dynamically between the site and its neighbors.

\begin{figure}[tbh]
\vspace{6mm}
\baselineskip=4pt
\begin{tabbing}
xxx\=xxx\=xxx\=xxx\=xxx\= \kill
\> PROCEDURE \it Sequential-Ants; \\
\>\>    {\tt Q $\leftarrow$ \{\} } \\
\>\>    select a site {\tt V} at random \\
\>\>    {\tt cluster $\leftarrow$ mark[V]}\\
\>\>	{\tt spin[V] $\leftarrow$ -spin[V]}\\
\>\>	{\tt {\it enque} V {\it into} Q}\\
\>\>	{\tt WHILE Q $\neq$ \{\} DO}\\
\>\>\>	  {\tt U $\leftarrow$  first(Q)}\\
\>\>\>	  {\tt {\it deque} U {\it from} Q}\\
\>\>\>	  make a bond for each of four nearest neighbors using probability p\\
\>\>\>	  {\tt FOR each neighbor W of U DO}\\
\>\>\>\>    {\tt IF mark[W] $\notin$ cluster THEN}\\
\>\>\>\>\>    {\tt cluster $\leftarrow$ mark[W]}\\
\>\>\>\>\>    {\tt spin[W] $\leftarrow$ -spin[W]}\\
\>\>\>\>\>    {\tt {\it enque} W {\it into} Q}\\
\>\>\>\>    {\tt END-IF}\\
\>\>\>	  {\tt END-FOR}\\
\>\>	{\tt END-WHILE}\\
\> END PROCEDURE \it Sequential-Ants\\

\end{tabbing}
%\vspace*{2cm}
\protect\caption { 
\protect\label {fig:SeqAnts}
The Sequential Ants-in-the-labyrinth Algorithm
}
%\vspace{1cm}
\end{figure}



\subsection{Analysis of algorithm}

Let the size of a cluster is $N$(that is, there are $N$ sites in a cluster). 
Each site in the cluster is visited exactly once, and we look at all 
four nearest neighbors to make a bond when we visit a site. Therefore 
the sequential ants algorithm takes a time $\Theta(4N)=\Theta(N)$. Note that 
the size of a cluster completely depends on the bond activation probability
$p$, and $p$ depends on the value of $\beta$. 

At large $\beta, p$ is large, so that the size $N$ is also large. 
In the sequential ants algorithm, the execution 
time only depends on the size of a cluster. But, in parallel ants 
algorithm, the execution time depends on not only the size but also 
the shape of a cluster. 
We shall explain this below.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Parallel(SIMD) Ants-in-the-labyrinth}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Unfortunately the ants scheme requires an unavoidable sequential 
process - each ant generation has to follow the previous one. 
The only parallelism is within a single generation - that is, over 
the wavefront of an expanding cluster. In the spanning tree, all sites 
within a single level can search their children in parallel. Therefore 
there is some parallelism, but it is very sparse, especially if the 
cluster size is very small when $\beta$ is small, so that the probability $p$ 
is small. Therefore, on parallel or vector machine, the Wolff algorithm 
is very inefficient. Also, we have serious load imbalance problem. 
If the number of virtual processor is $N^2$, then the maximum number of 
active virtual processor is $O(N)$. This ants works better at large $\beta$, 
where the average cluster size is large, than at small $\beta$, where the 
small cluster sizes severely limit the amount of parallelism which 
can be obtained~\cite{CP,PL}. An example is shown in Figure~\ref{fig1}.


%\vspace{5mm}
\pssilent
%**********************************************************
\begin{figure}[h]
%  \leavevmode
%  \epsfxsize=6in
\psfig{figure=pal.eps,height=2in,width=5.5in}
%  \epsfbox{pal.eps}
  \caption{Parallel ants algorithm}
  \label{fig1}
%\vspace{10mm}
\end{figure}
%*********************************************************


\begin{figure}[tbh]
\vspace{11mm}
\baselineskip=10pt
\begin{tabbing}
xxx\=xxx\=xxx\=xxx\=xxxxxxxx\= \kill
\> PROCEDURE \it Parallel-Ants; \\
\>\>    generate random number \it R \\
\>\>    set up right-neighbors {\it Nr} and left-neighbors {\it Nl}
         as using \it R\\
\>\>    generate random number \it R\\
\>\>    set up down-neighbors {\it Nd} up-neighbors {\it Nu}
         as using \it R\\
\>\>    {\tt exit $\leftarrow$ FALSE} \\
\>\>    {\tt {\it Cluster} $\leftarrow$ FALSE} \\
\>\>    generate random number {\it i} and {\it j}\\
\>\>    {\tt {\it Children[i,j]} $\leftarrow$ TRUE} \\
\>\>	{\tt WHILE ($\neg$exit) DO}\\
\>\>\>	  {\tt \it Parent $\leftarrow$ Children}\\
\>\>\>	  {\tt \it Cluster $\leftarrow$ Cluster $\vee$ Parent}\\
\>\>\>	  {\tt {\it Children} $\leftarrow$ 
                  ({\it Nr} $\wedge$ SHIFT ({\it Parent},right)) $\vee$}\\
\>\>\>\>\>   {\tt ({\it Nl} $\wedge$ SHIFT ({\it Parent}, left)) $\vee$}\\
\>\>\>\>\>   {\tt ({\it Nd} $\wedge$ SHIFT ({\it Parent}, down)) $\vee$}\\
\>\>\>\>\>   {\tt ({\it Nu} $\wedge$ SHIFT ({\it Parent}, up)) }\\
\>\>\>	  {\tt \it Children $\leftarrow$ Children $\wedge$ $\neg$Cluster}\\
\>\>\>	  {\tt IF} there are no more new children {\tt THEN} \\
\>\>\>\>	  {\tt exit $\leftarrow$ TRUE}\\
\>\>\>	  {\tt END IF}\\
\>\>	  {\tt END WHILE }\\
\>\>	  {\tt WHERE (\it Cluster)}\\
\>\>\>	  {\tt \it Spin $\leftarrow$ - Spin}\\
\>\>	  {\tt END WHERE }\\
\> END PROCEDURE \it Parallel-Ants 

\end{tabbing}
%\vspace*{2cm}
\protect\caption { 
\protect\label {fig:ParAnts}
The Parallel Ants-in-the-labyrinth Algorithm
}
%\vspace{10mm}
\end{figure}


%=====================================================================
\subsection{Algorithm}

The algorithm is described by the pseudo-code in Fig.~\ref{fig:ParAnts}.
In the sequential ants algorithm, each bond is dynamically made when
 we look at the nearest neighbors according to the probability $p$. 
Note that we use random number for the sampling method with probability $p$. 
In the sequential program, this can be done easily and efficiently. 
But in the parallel program, this is not good method. At each iteration, 
we should know how many sites are within a generation, and then we should 
use logical mask array to generate random numbers. This will decrease 
efficiency. Instead we make all bonds through the entire lattice at first, 
as generating random numbers in parallel. Then the next step(growing a
cluster) is performed in parallel. 
At each generation level(iteration), each site within a generation 
search its nearest neighbors by using already made bonds information, 
and all sites within a generation perform this search in parallel. 
This is only parallelism which we can expect.


%=======================================================================
\subsection{Analysis of algorithm}

The time cost of the sequential ants algorithm is dependent on only the 
size of a cluster. It is completely independent on the shape of a cluster.
But, in the parallel ants, the factor-the shape of cluster is  very 
important to decide the time cost. 

As mentioned before, the only parallelism is within a single 
generation(i.e., within a single level of a spanning tree). Therefore, 
if the number of generations is $H$(i.e. the height of a spanning tree 
= the longest length from the root to a leaf), this parallel ants 
algorithm takes a time $\Theta(H)$. The problem is that the height of a 
spanning tree $H$ depends on two factors - the size of a cluster and
the shape of a cluster. 
The size and shape of a cluster is dependent on the bond activation 
probability $p$. Note that the probability $p$ is dependent on $\beta$ 
and we can think $\beta$ as the inverse temperature. 

Therefore, the speedup $S$ is 

\hspace{4cm}   $S = \frac{\Theta(N)}{\Theta(H)}$ 
               where $N$ is the size of a cluster.


Let $N$ be fixed. Then the best case is when the spanning tree is almost 
balanced - i.e., the shape of a cluster is almost regular . For any $H$, 
the maximum number of sites are $2H(H+1)$. On the other hand, 
the worst case is when the spanning tree is in the form of a chain. 
Hence, the algorithm takes a time

\hspace{4cm} $\Omega(H) = \Omega(\sqrt{\frac{N}{2}})$
                       and $O(H) = O(N)$.



Unfortunately the size and shape of a cluster in the Wolff algorithm is 
highly irregular, so that we can not expect the best case. In the best case, 
the speedup is

\hspace{4cm}   $S = \frac{\Theta(N)}{\Omega(\sqrt{N})} = O(\sqrt{2N})$.

But, in the worst case, we can not expect any speedup.



%=====================================================================
\subsection{Performance}

We have executed the program on the CM-5 with vector unit to update various 
lattice sizes($64^2$, $128^2$, $256^2$, $512^2$) at various values 
of the inverse temperature $\beta=$(0.3, 0.4, $\beta_c$, 0.5, 0.6) to 
understand how the algorithm behave as the cluster size distribution changes. 

We have executed at least 1000 sweeps to have some statistics. The critical 
inverse temperature $\beta_c$ of the phase transition for the two-dimension 
spin model is known exactly $\beta_c = \ln(\sqrt{2}+1$)=0.4406868...
~\cite{CP}. To compare the time taken at various $\beta$, we used the time 
per site for one sweep(time for one sweep divided by average cluster size) 
rather than the time for one sweep, since average cluster size is very 
different according to the value of $\beta$. The results of lattice $256^2$
are in Figure~\ref{fig:pal3},Figure~\ref{fig:pal4},Figure~\ref{fig:pal5}
and Figure~\ref{fig:pal6}.

From the result, we can see this algorithm works better at large 
$\beta(>\beta_c)$ - i.e., average cluster is large, since the only parallelism
is within a single generation. Therefore the time per site for one sweep is 
inversely proportional to the cluster size, and we have more amount of 
parallelism at large $\beta$.

As mentioned before, this algorithm takes a time $\Theta(H)$, and speedup 
$S=\Theta(N)/\Theta(H)$ where $N$ is the size of a cluster and $H$ is the 
number of heights of the spanning tree(i.e., the number of iterations in 
the program). From some statistical result, we can estimate speedup roughly 
by dividing an average cluster size by an average number of iterations. 
The result is in Figure~\ref{fig:pal7}. At $\beta$=0.6, this estimated 
speedup is almost $N^{1/2}(N=cluster size)$. 
Since the best speedup is $O(\sqrt{2N})$ , 
we have very good performance(almost the best case). The reason is that a 
cluster at large $\beta$ like 0.6 is more regular in shape, and its size 
is almost equal to the whole lattice size.
We have also measured the time taken as a function of lattice size at 
$\beta=\beta_c$, the result is in Table~\ref{pal1}..





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Parallel Multiple Ants}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In the Wolff algorithm, we make bonds simultaneously as growing a single 
cluster. That means we do not need to make all bonds through all lattice, 
but need to make bonds within the single cluster. Note that we make a bond 
between two sites using a random number. In the previous parallel ants 
program, we make all bonds before growing a cluster. If the lattice size 
is $L^2$, we need $2L^2$ random numbers, since each site needs two random 
numbers(one is for right-left direction, and the other is up-down direction). 
During one sweep of parallel ants program, we just use some part of all 
bonds(that is, some parts of random numbers) which belong to a single 
cluster.
 
Hence we can use the unused bonds(unused random numbers) which do not belong 
to the single cluster during the next sweep if the next cluster is not 
connected to the previous cluster(that is, two clusters are independent 
each other). Update of a single cluster is completely independent to the 
other sites which do not belong to the cluster. 

Therefore we can perform two sweeps of parallel ant algorithm in parallel 
if two selected clusters are not connected. If two clusters are connected, 
we have some problem. In such case, we need to make new bonds(generate 
new random numbers) again, and then perform the ants algorithm again 
starting from the same initial site. 

In the same way, we can do more than two sweeps in parallel, checking any two 
clusters are connected. If any two clusters are connected, we merge two 
clusters into one cluster, and then save the initial random site of the 
merged cluster for the next step. After one step is finished, we do again 
for the merged clusters starting from the same initial site. The previous 
ants works better at large $\beta$ (i.e., large cluster), but does not works  
well at small $\beta$ (i.e., small cluster). On the contrary, this new 
algorithm works very well at small $\beta$, since the possibilities that any 
two clusters are connected is very low. Some examples are shown in 
Figure~\ref{fig2}.


%=======================================================================
\subsection{Algorithm}

The algorithm is described by the pseudo-code in Fig.~\ref{fig:MulAnts}.

This algorithm consists of three phase:

\begin{itemize}
\item {\it Phase} 1: Select $n$ initial random sites, and label them.

\item {\it Phase} 2: Grow $m$ multiple clusters from the corresponding 
                     initial random sites.

\item {\it Phase} 3: Merge any connected clusters, and decide which 
                     clusters should be grown again in the next step.
\end{itemize}

In {\it phase} 1, we assign a label to each of $n$ initial random sites. 
During growing multiple clusters, each cluster corresponds to one of $n$
initial sites, and all sites in a cluster should be assigned the same 
label as the corresponding initial site.

In {\it phase} 2, growing all $m$ clusters is same as in the previous parallel 
ants. Only difference is that we need to assign the same label to a new 
child as its parent's label, because we must distinguish $m$ clusters. 
Hence, if any site $i$ belongs to a cluster and this cluster is grown form 
the initial site $j$ whose label is $l$, the label of $i$ is same as $j$.


In {\it phase} 3, we need some sequential step. First, we find all sites which 
belong to more than one cluster, and then make merge table by using labels 
of those sites. We then optimize this table to decide which clusters are 
merged into and should be grown again in the next step(This step is a 
sequential step). Note that since we have already selected n random sites, 
any merged cluster to be grown again in the next step should be grown from 
the corresponding initial random sites.

\vspace{20mm}
%**********************************************************
%\clearpage
\begin{figure}[bh]
%  \begin{center}
%\begin{boxit}
%  \leavevmode
%  \epsfxsize=5in
%  \epsfxsize=5in
%  \epsfbox{mal.eps}
\psfig{figure=mal.eps,height=6.1in,width=6.1in}
%  \end{center}
  \caption{Parallel Multiple ants algorithm}
  \label{fig2}
%\vspace{10mm}
%\end{boxit}
\end{figure}
%*********************************************************


\begin{figure}[tbh]
%\vspace{6mm}
\baselineskip=5pt
\begin{tabbing}
xxx\=xxx\=xxx\= \kill

\>  Procedure {\it Multiple-Ants;} \\
\>\>  Initially select $n(>1)$ random sites.\\
\>\>  Record these $n$ initial sites for the next step.\\
\>\>  {\tt REPEAT} \\
\>\>\>   Decide the number of clusters, m to be grown in parallel,\\
\>\>\>   Initially $m = n$, and after the first step,\\
\>\>\>   $m$ = number of the merged clusters.\\
\>\>\>   Make all bonds.\\
\>\>\>   Grow all $m$ clusters simultaneously starting from \\
\>\>\>   the corresponding initial random site.\\
\>\>\>   After all $m$ clusters are made, flip all spins in the clusters.\\
\>\>\>   Find which clusters are connected.\\
\>\>\>   Merge all connected clusters.\\
\>\>  {\tt UNTIL} there is no merged cluster.\\
\>  End-Procedure\\

\end{tabbing}
%\vspace*{2cm}
\protect\caption { 
\protect\label {fig:MulAnts}
The Parallel Multiple Ants algorithm
}
\vspace{1.5cm}
\end{figure}




%======================================================================
\subsection{Comparison with parallel ants}

We define speedup $S$(parallel ants to parallel multiple ants):

\hspace{2cm}  $S = N + A - B - C - D$

  where  $N$ = sweep unit(number of sweeps to be executed in parallel)

\hspace{1cm} $A$ = speed up for growing a cluster in case that some clusters

\hspace{2cm}  are connected in the cluster.

\hspace{1cm} $B$ = number of steps(iterations), $1 \leq B \leq N$ 

\hspace{1cm} $C$ = overhead for labeling clusters

\hspace{1cm} $D$ = overhead for merging clusters


In this algorithm, we expect better load balance than parallel ants, 
since we grow multiple clusters in parallel. The speedup $S$ depends on the 
size of cluster(i.e., $\beta$) like parallel ants. While the previous 
parallel ants works better at large $\beta$(large cluster size), this 
multiple ants works better at small $\beta$(small cluster size) since we 
have small number of steps($B$) and less overhead for labeling ($C$) and 
merging ($D$) clusters. 

Suppose that $N$ is fixed. Then, as $\beta$ increases, $A$ increases, 
but $B$, $C$, and $D$ also increase. On the other hand, as $\beta$ decrease, 
$A$ decrease, but $B$, $C$, and $D$ also decrease. 
Therefore, we expect good speedup $S$ at small $\beta$.



%=======================================================================
\subsection{Performance}

We have executed the program on the CM-5 with vector unit to update various 
lattice sizes($64^2$, $128^2$, $256^2$, $512^2$) at various values 
of the inverse temperature $\beta$(0.3, 0.4, $\beta_c$, 0.5, 0.6) to 
understand how the algorithm behave as the cluster size distribution changes. 
Also, we have executed the program with various {\it sweep unit N}(=5, 10, 20) 
to compare speedup. 
The time per site for one sweep(time for one sweep divided by average 
cluster size) taken on the CM-5 to update a $256^2$ lattice are shown in 
Figure~\ref{fig:mal6}. 
This figure shows that, at small $\beta$, as the {\it sweep unit(N)}
increases, average time per site for one sweep decreases. 
But, at large $\beta$, average time per site is almost same although $N$
increases.

The performances are much better at small $\beta(<\beta_c)$ than the previous 
parallel ants algorithm. The reason is that we have small number of 
steps(factor $B$) and less overhead for labeling($C$) and merging($D$) 
clusters.
But the performances are poor at large $\beta(>\beta_c)$ and $\beta_c$.
The speedup factor($A$) does not seems to play an important role. 
Note that on the contrary, the parallel ants worked well at large $\beta$. 
The average number of steps(factor $B$)are shown in Figure~\ref{fig:mal7}. 
 At large $\beta$, the number of steps($B$) is almost equal to the sweep 
unit($N$) since the average cluster size is almost same as the whole 
lattice size. 
In this case, we have serious speed-down factor $B$.

We measured the time taken as a function of lattice size at $\beta=\beta_c$,
the result is shown in Table~\ref{mal5}. 




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion and Discussion}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In this project, we have implemented two different parallel(SIMD) algorithms 
for the Wolff algorithm with CM Fortran, and compared their performances 
on the CM-5. 
For each of two algorithms, the time per site for one sweep(time for one 
sweep divided by average cluster size) taken on the CM-5 to update a $64^2$, 
$128^2$, $256^2$, $512^2$ lattice and the time per site for one sweep of 
these two algorithms is shown in Figure~\ref{fig:mal6} for $256^2$ and in 
Table~\ref{con1} for $512^2$.

At small $\beta(<\beta_c)$, the parallel multiple ants labeling algorithm have 
the best performance. 
On the other hand, at large $\beta(>\beta_c))$, the parallel ants has the best 
performance.

 As mentioned before, the Wolff algorithm unfortunately requires an 
unavoidable sequential process - i.e., each ant generation has to follow the 
previous one. The only parallelism is within a single generation - that is, 
over the wavefront of an expanding cluster. To overcome this problem, 
we developed the parallel multiple ants , so that we have improved the 
performances at small $\beta(<\beta_c)$ and $\beta_c$. But at large 
$\beta(>\beta_c)$, unfortunately their performance are poor, 
so that the parallel ants still has the best performance.

In this project, we have implemented SIMD algorithms. However, even in MIMD 
machine, the Wolff algorithm will suffer from more serious load imbalance 
problem and highly irregular communication. The dynamic domain 
decomposition might be one of the solutions.





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{thebibliography}{00}
%-----------------------
\bibitem{MONTECARLO}
K.~Binder ed., 
 Monte Carlo Methods in Statistical Physics, 
(Springer-Verlag, Berlin, 1986);
K.~Binder and D.W.~Heermann,  Monte Carlo Simulation in 
Statistical Physics, (Springer-Verlag, Berlin, 1988);
H.~Gould and J.~Tobochnik,
 An Introduction to Computer Simulation Methods, Vol. 2,
(Addison-Wesley, Reading, Mass., 1988).
%-----------------------
\bibitem{METROPOLIS}
N.~Metropolis  et al., 
 J. Chem. Phys. 21 (1953) 1087.
%-----------------------
\bibitem{SOKALREV}
A.~D.~Sokal, in Computer Simulation Studies in Condensed Matter
Physics:  Recent Developments, eds.\ D.~P.~Landau et al.
(Springer-Verlag, Berlin, 1988);
A.~D.~Sokal, % How to Beat Critical Slowing Down: 1990 Update,
in Proc. of the International Conference on Lattice Field Theory,
Tallahassee, October 1990,  Nucl. Phys. B (Proc. Suppl.) 
20 (1991) 55.
%-----------------------
\bibitem{SW}
R.H.~Swendsen and J.-S.~Wang, 
% Nonuniversal critical dynamics in Monte Carlo simulations,
 Phys. Rev. Lett. 58 (1987) 86.
%-----------------------
\bibitem{WOLFF}
U. Wolff, 
 Phys. Rev. Lett. 62 (1989) 361.
%-----------------------
\bibitem{LAT89REV}
U.~Wolff,
in Proc. of the Symposium on Lattice Field Theory, Capri, September
1989,  Nucl. Phys. B (Proc. Suppl.) 17 (1990) 93. 
%-----------------------
\bibitem{SWREV}
J.-S.~Wang and R.~H.~Swendsen,  Physica A 167 (1990) 565.
%-----------------------
\bibitem{CP}
Clive F.Baillie and Paul D. Coddington; "Cluster Identification 
Algorithms for Spin Models-Sequential and Parallel", June, 1990.
%-----------------------
\bibitem{PL}
Paul Coddington, Lecture Note of CPS713.
%-----------------------
\bibitem{DA}
D. W. Heermann and A. N. Burki, "Parallel Algorithms in Computational 
Science" (Springer-Verlag).
%-----------------------
\bibitem{PRW}
P. Tamayo, R. C. Brower, and W. Klein, J. Statis. Phys. 58, 1083(1990).
%-----------------------
\bibitem{GP}
G. Brassard and P. Bratley, "Algorithmics: Theory \& Practice"
(Prentice Hall, 1987).
%-----------------------
\bibitem{COMPSCI}
E.~Horowitz and S.~Sahni, Fundamentals of Computer Algorithms,
(Computer Science Press, Rockville, Maryland, 1978);
G.~Brassard and P.~Bratley, Algorithmics: Theory and Practice,
%-----------------------
\bibitem{CMF}
CMFortran Reference Manual, (Thinking Machines Corporation,
Cambridge, Mass., 1993).

%

\end{thebibliography}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\newpage

\vspace{0.5in}

\begin{table}

\begin{tabular}{|c||c||c|}  \hline
 Lattice Size&Time per One Sweep($msec$)&Time for Site per One Sweep($msec$)\\ \hline\hline
 64          & 64.84                    & 0.0352     \\ \hline
 128         & 202.19                   & 0.0397   \\ \hline
 256         & 879.34                   & 0.0440   \\ \hline
 512         & 4547.09                  & 0.0696  \\ \hline
\end{tabular}
\caption{Parallel Ants algorithm at $\beta_c$}
\label{pal1}
\end{table}


\begin{table}
\begin{tabular}{|c||c|c|c||c|c|c|}  \hline
  &  \multicolumn{3}{|c||}{ Time per One Sweep($msec$)} 
  & \multicolumn{3}{  c|}{Time for Site per One Sweep($msec$)}\\ \cline{2-7}
 Lattice Size &5       &10      &20      &5      &10     &20\\ \hline\hline
 64           &75.67   &61.16   &41.99   &0.0618 &0.0584 &0.0578 \\ \hline
 128          &247.84  &165.81  &128.96  &0.0588 &0.0552 &0.0488 \\ \hline
 256          &1020.06 &785.87  &515.80  &0.0861 &0.0743 &0.0644 \\ \hline
 512          &6529.00 &5276.51 &3187.25 &0.1438 &0.1229 &0.1072 \\ \hline
\end{tabular}
\caption{Parallel Multiple Ants algorithm(sweep=5,10,20) at $\beta_c$}
\label{mal5}
\end{table}



\begin{table}
\begin{tabular}{|c||c|c|c||c|c|c|}  \hline
&\multicolumn{2}{|c|}{Time for Site per One Sweep($msec$)} \\ \cline{2-3}
 $\beta$     &parallel ants  &parallel multiple ants($sweep=20$)\\ \hline\hline
 0.3          &20.641          &1.686    \\ \hline
 0.4          &9.612           &0.927   \\ \hline
 $\beta_c$    &0.0696          &0.107   \\ \hline
 0.5          &0.0185          &0.026  \\ \hline
 0.6          &0.0160          &0.023  \\ \hline
\end{tabular}
\caption{Parallel Ants vs Parallel Multiple Ants algorithm with size $512^2$}
\label{con1}
\end{table}





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%\vspace{15mm}
%**********************************************************
\clearpage


\pssilent
%**********************************************************

\begin{figure}
%  \leavevmode
%  \epsfxsize=6in
\centerline{
\psfig{figure=cluster.eps,height=3in,width=4.5in}
}
\protect\caption{\protect\label{fig:pal3}
Average cluster size vs beta value for 
the parallel ants algorithm on CM5 with vector unit.}
\end{figure}

%\vspace{15mm}

\begin{figure}
\centerline{
\psfig{figure=iter.eps,height=3in,width=4.5in}
}
\protect\caption{\protect\label{fig:pal4}
Average iteration  vs beta value for the parallel ants 
algorithm on CM5 with vector unit.
}
\end{figure}
%*********************************************************



%**********************************************************
\begin{figure}[th]
%  \leavevmode
%  \epsfxsize=6in
\centerline{
\psfig{figure=tos.eps,height=3in,width=4.5in}
}
\protect\caption{\protect\label{fig:pal5}
Average time for one sweep vs beta value for the parallel ants 
algorithm on CM5 with vector unit.
}
\end{figure}

\vspace{15mm}

\begin{figure}[th]
\centerline{
\psfig{figure=tss.eps,height=3in,width=4.5in}
}
\protect\caption{\protect\label{fig:pal6}
Average time per site for one sweep vs beta value for the parallel ants 
algorithm on CM5 with vector unit.
}
\end{figure}
%*********************************************************


%**********************************************************
\begin{figure}[th]
%  \leavevmode
%  \epsfxsize=6in
\centerline{
\psfig{figure=sit.eps,height=4in,width=5.3in}
}
\protect\caption{\protect\label{fig:pal7}
Average site divided by average iteration for one sweep vs beta value 
for the parallel ants algorithm on CM5 with vector unit.
}
\end{figure}
%*********************************************************

\newpage
\clearpage
%\vspace{30mm}

%**********************************************************
\begin{figure}[t]
%  \leavevmode
%  \epsfxsize=6in
\centerline{
\psfig{figure=w33tss.eps,height=3in,width=4.5in}
}
\protect\caption{\protect\label{fig:mal6}
Average time per site for one sweep vs beta value for the parallel
multiple ants  algorithm with different sweep-unit(5,10,20)
on CM5 with vector unit.
}
\end{figure}

%\vspace{15mm}

\begin{figure}
\centerline{
\psfig{figure=w33nos.eps,height=3in,width=4.5in}
}
\protect\caption{\protect\label{fig:mal7}
Average number of steps vs beta value for the parallel
multiple ants  algorithm with different sweep-unit(5,10,20)
on CM5 with vector unit.
}
\end{figure}
%*********************************************************




\end{document}









