\documentstyle[12pt]{book}

\begin{document}

\section{Quantum Chromodynamics}
\label{qcd}

\subsection{Introduction}

Quantum Chromodynamics\index{Quantum Chromodynamics} is the proposed theory 
of the so called strong interactions that mediate the forces between protons, 
neutrons, and mesons and so control and formation of nuclei.  Quantum 
Chromodynamics is always abbreviated QCD.  The fundamental properties of QCD 
cannot be directly tested, but a wealth of indirect evidence supports this 
theory.

Gauge theories\index{Gauge theories} are ubiquitous in elementary particle 
physics\@: the electromagnetic interaction between electrons and photons is 
described by quantum electrodynamics (QED) based on the gauge group U(1), the 
strong force between quarks and gluons is believed to be explained by quantum
chromodynamics (QCD) based on SU(3), and there is a unified 
description of the weak and electromagnetic interactions in terms of the 
gauge group $\rm SU(2)\times U(1)$.  The strength of these interactions is 
measured by a coupling constant.  This coupling constant is small for QED, so 
very accurate analytical calculations can be performed using perturbation 
theory.  However, for QCD, the coupling constant appears to increase with 
distance so that perturbative calculations are only possible at short 
distances.  In order to solve QCD at longer distances, Wilson 
\cite{Wilson:74a} introduced lattice gauge theory\index{lattice gauge
theory} in which the space-time continuum is discretized to provide a cut-off 
that regulates ultraviolet divergences non-perturbatively.  This 
discretization onto a lattice, which is typically hypercubic, also makes the 
gauge theory amenable to numerical simulation by computer.

{\small
\begin{table}
\begin{center}
\begin{tabular}{|l|l|c|l|}
\hline
\multicolumn{1}{|c|}{Authors}&\multicolumn{1}{c|}{Project}&
\multicolumn{1}{c|}{Hypercube}&\multicolumn{1}{c|}{Reference}
\\
\hline
\medStrut
Brooks, Fox, Otto, & Glueball mass (SU(2)) & Cosmic Cube &
\protect\cite{Brooks:83a} \\
Randeria, Athas, & $4^3 \times 8$ lattice & 4 node & \\
De Benedictis, & & & \\
Newton, Seitz & & & \\
\hline
Otto, Randeria & Glueball mass, modified action & Cosmic Cube &
\protect\cite{Otto:83a} \\
& $4^3 \times 8$ lattice & 4 node & \\
\hline
Otto, Stack & Heavy quark potential & Cosmic Cube & \protect\cite{Otto:84a}
\\
& $12^3 \times 16$ lattice & 64 node & \\
\hline
Otto, Stolorz & Glueball mass, enhanced statistics & Cosmic Cube &
\protect\cite{Otto:85b} \\
& $12^3 \times 16$ lattice & 64 node & \\
\hline
Patel, Otto, Gupta & Monte Carlo renormalization group, & Cosmic Cube &
\protect\cite{Patel:85a} \\
& nonperturbative $\beta$-function & 64 node & \\
\hline
Fucito, Soloman & Chiral symmetry breaking, mass & Mark II &
\protect\cite{Fucito:84a}, \protect\cite{Fucito:85b}, \\
& spectrum, deconfinement transition, & 64 node &
\protect\cite{Fucito:85c}, \protect\cite{Fucito:85d} \\
& pseudo-fermions & & \protect\cite{Fucito:86a} \\
\hline
Flower, Otto & Field distribution & Mark II & \protect\cite{Flower:85a}
\\
& $10^2 \times 12^2$ lattice & 32 node & \\
\hline
Flower, Otto & Heavy quark potential & Mark II &
\protect\cite{Flower:86b} \\
& $20^4$ lattice & 128 node & \\
\hline
Kolawa, & Glueball mass (SU(2)), & Mark II &
\protect\cite{Furmanski:87c} \\
Furmanski & Hamiltonian ``loop'' formalism & 32 node & \\
\hline
Stolorz, Otto & Microcanonical renormalization & Mark I &
\protect\cite{Stolorz:86b} \\
& group (SU(2)) & 64 node & \\
\hline
Flower & Baryons on the lattice $20^4$ & Mark II &
\protect\cite{Flower:87e} \\
& & 128 node & \\
\hline
Chiu & Random block lattice & nCUBE-1 & \protect\cite{Chiu:89b},
\protect\cite{Chiu:88c}, \\
& & 1024 node & \protect\cite{Chiu:88e}, \protect\cite{Chiu:88f} \\
& & & \protect\cite{Chiu:89a} \\
\hline
Ding, Baillie, & Heavy quark potential & Mark IIIfp &
\protect\cite{Ding:90b} \\
Fox & $24^3 \times 10$, $24^3 \times 12$, $12^3 \times 24$, & 32 node & \\
& $16^3 \times 24$ & and 64 node & \\
\hline
Baillie, Brickner, & QCD with dynamical Wilson & TMC CM-2 &
\protect\cite{Baillie:89e} \\
Gupta & fermions $16^4$ lattice & 16K & \\
\hline
\end{tabular}
\caption{Quantum Chromodynamic (QCD) Calculations within C$^3$P}
\label{tab1.qcd}
\end{center}
\end{table}
}

Most of the work on lattice gauge theory has been directed towards solving
lattice QCD and, thus, deriving the hadron mass spectrum from first 
principles.  This would confirm QCD as the theory of the strong force.  Other 
calculations have also been performed, in particular, the properties of QCD 
at finite temperature and/or finite baryon density have been determined.  
Unfortunately, in order to simulate lattice QCD on a computer, one must 
integrate out the quark variables (because they are fermions, 
{\em i.e.,}\@ anticommuting elements of a Grassmann algebra rather than 
numbers) leaving a highly non-local fermion determinant for each flavor of 
quark.  Physically, this determinant arises from closed quark loops.  The 
simplest way to proceed is to ignore these quark loops and work in the 
so-called quenched approximation with no flavors of quark present.  (This may 
be valid for heavy quarks.)  Current state-of-the-art quenched QCD 
calculations are performed on lattices of size $24^3\times 48$.

However, to investigate the physically more realistic fully interacting
QCD, with the inclusion of dynamical quarks, we must go beyond the quenched 
approximation and tackle the problem of the fermion determinant.  We shall 
return to this later.

\subsection{Monte Carlo}

In order to explain the computations for QCD, we use the Feynman path 
integral formalism.  For any field theory described by a 
Lagrangian density ${\cal L}$, the dynamics of the fields $\phi(x)$ are 
determined through the action functional
\begin{equation}
S(\phi) = \int dx {\cal L}(\phi).
\end{equation}
In this language, the measurement of a physical observable represented by an 
operator ${\cal O}$ is given as the expectation value
\begin{equation}
\langle {\cal O}(\phi)\rangle = \frac{1}{Z} \int D\phi {\cal O}(\phi) 
e^{-S(\phi)},
\label{eqn2.qcd}
\end{equation}
where the partition function\index{partition function} $Z$ is
\begin{equation}
Z = \int D\phi e^{-S(\phi)}.
\end{equation}
In these expressions, the integral $D\phi$ indicates a sum over all possible
configurations of the field $\phi$.  A typical observable would be a product 
of fields ${\cal O}=\phi(x)\phi(y)$, which says how the fluctuations in the 
field are correlated, and in turn, tells us something about the particles 
that can propagate from point $x$ to point $y$.  The appropriate correlation 
functions give us, for example, the masses of the various particles in the 
theory.  Thus, to evaluate most any quantity in field theories like QCD,
one must simply evaluate the corresponding path integral.  The catch is that 
the integrals are over an infinite-dimensional space.

To put the field theory onto a computer, we begin by discretizing space and 
time into a lattice of points.  Then the functional
integral is simply defined as the product of the 
integrals over the fields at every site of the lattice $\phi(n)$\@:
\begin{equation}
\int D\phi \rightarrow \int \prod_n \lbrack d\phi(n)\rbrack.
\end{equation}
Restricting space and time to a finite box, we end up with a finite (but 
large) number of ordinary integrals, something we might imagine doing 
directly on a computer.  However, the high multi-dimensionality of these 
integrals renders conventional mesh techniques impractical.  Fortunately, the 
presence of the exponential $e^{-S}$ means that the integrand is sharply 
peaked in one region of configuration space.  Hence, we resort to a 
statistical treatment and use Monte Carlo type algorithms to sample the 
important parts of the integration region.

Monte Carlo algorithms typically begin with some initial 
configuration of fields, and then make pseudo-random changes on the fields 
such that the ultimate probability $P$ of generating a particular field 
configuration $\phi$ is proportional to the Boltzmann factor,
\begin{equation}
P(\phi) = e^{-S(\phi)},
\end{equation}
where $S(\phi)$ is the action associated with the given configuration.
There are several ways to implement such a scheme, but for many theories the 
simple Metropolis algorithm\index{Metropolis algorithm} \cite{Metropolis:53a} 
is effective.  In this algorithm, a new configuration $\phi'$ is generated by 
updating a single variable in the old configuration $\phi$ and calculating 
the change in energy (action)
\begin{equation}
\Delta S = S(\phi') - S(\phi).
\end{equation}
If $\Delta S \le 0$, the change is accepted; if $\Delta S > 0$ the change
is accepted with probability $\exp(- \Delta S)$.  In practice, this is done 
by generating a pseudo-random number $r$ in the interval [0,1] with uniform 
probability distribution and accepting the change if $r < \exp(- \Delta S)$.
This is guaranteed to generate the correct (Boltzmann) distribution of 
configurations, provided ``detailed balance'' is satisfied.  That condition 
means that the probability of proposing the change $\phi\rightarrow\phi'$ is 
the same as that of proposing the reverse process $\phi'\rightarrow\phi$.
In practice, this is true if we never simultaneously update two fields which 
interact directly via the action.  Note that this constraint has important 
ramifications for parallel computers as we shall see below.

Whichever method one chooses to generate field configurations, one updates 
the fields for some equilibration time of $E$ steps, and then calculates the 
expectation value of ${\cal O}$ in equation~\ref{eqn2.qcd} from the next $T$ 
configurations as
\begin{equation}
\langle {\cal O}\rangle = \frac{1}{T} \sum_{i=E+1}^{E+T} {\cal O}(\phi_i).
\end{equation}
The statistical error in Monte Carlo behaves as $1/\sqrt{N}$, where $N$ is 
the number of effectively independent configurations.  $N = T/2\tau$, where 
$\tau$ is the autocorrelation time.  This autocorrelation time can easily be 
large, and most of the computer time spent in generating statistically 
independent configurations.  The operator measurements then become a small 
overhead on the whole calculation.

\subsection{QCD}

Quantum Chromo-dynamics (QCD) is the theory describing the strong interaction 
in high energy physics. This interaction binds quarks together to form 
hadrons---the constituents of nuclear matter---for example, the proton and 
neutron.  Currently, we know of five types (referred to as ``flavors'') of 
quark: up, down, strange, charm and bottom; and expect at least one more 
(top) to show up soon.  In addition to having a ``flavor,'' quarks can carry 
one of three possible charges known as ``color'' (this has nothing to do with 
color in the macroscopic world!); hence, Quantum {\em Chromo\/}-dynamics.
The strong color force is mediated by particles called gluons, just as 
photons mediate light in electromagnetism.  Unlike photons, though, gluons 
themselves carry a color charge and, therefore, interact with one another.
This makes QCD a nonlinear theory which is impossible to solve
analytically.  Therefore, we turn to the computer for numerical
solutions.

\subsection{Lattice QCD}

To put QCD on a computer, we proceed as follows\@:  Four-dimensional 
space-time continuum is replaced by a four-dimensional hypercubic periodic 
lattice, of size $N = N_s\!\times\! N_s\!\times\! N_s\!\times\! N_t$ with the 
quarks living on the sites and the gluons living on the links of the lattice.
$N_s$ is the spatial and $N_t$ is the temporal extent of the lattice.  The 
lattice has a finite spacing $a$.  The gluons are represented by 
$3\!\times\! 3$ complex $SU(3)$ matrices associated with each link in the 
lattice.  The $3$ in $SU(3)$ reflects the fact that there are three colors of 
quarks, and $SU$ means that the matrices are unitary with unit determinant 
({\em i.e.,} ``special unitary'').  This link matrix describes how the color 
of a quark changes as it moves from one site to the next.  For example, as a 
quark is transported along a link of the lattice it can change its color 
from, say, red to green; hence, a red quark at one end of the link can 
exchange colors with a green quark at the other end.  The action functional 
for the purely gluonic part of QCD is
\begin{equation}
S_G = \beta \sum_P\, ( 1 - \frac{1}{3} Re Tr U_P ),
\end{equation}
where
\begin{equation}
U_P = U_{i,\mu}\, U_{i+\mu,\nu}\, U_{i+\nu,\mu}^\dagger\, 
U_{i,\nu}^\dagger
\end{equation}
is the product of link matrices around an elementary square or plaquette on 
the lattice --- see Figure~\ref{fig1.qcd}.  Essentially, all of the time in 
QCD simulations of gluons is spent multiplying these $SU(3)$ matrices 
together.  The main component of this is the $a \times b + c$ kernel (which 
most supercomputers can do very efficiently).  As the action involves 
interactions around plaquettes, in order to satisfy detailed balance, we can 
update only half the links in any one dimension simultaneously, as shown in
Figure~\ref{fig2.qcd} (in two dimensions for simplicity).  The partition 
function for full lattice QCD including quarks is
\begin{equation}
Z = \int D\psi D\bar\psi DU \exp(-S_G - \bar\psi (\not\!\!D +m) \psi),
\end{equation}
where $\not\!\!D +m$ is a large sparse matrix the size of the lattice squared.
Unfortunately, since the quark or fermion variables $\psi$ are anticommuting 
Grassmann numbers, there is no simple representation for them on the computer.
Instead, they must be integrated out, leaving a highly non-local fermion 
determinant:
\begin{equation}
Z = \int DU \det(\not\!\!D +m) \exp(-S_G).
\end{equation}
This is the basic integral one wants to evaluate numerically.

\begin{figure}
\vspace*{1.0in}
\caption{A Lattice Plaquette}
\label{fig1.qcd}
\end{figure}

\begin{figure}
\vspace*{1.0in}
\caption{Updating the Lattice}
\label{fig2.qcd}
\end{figure}

The biggest stumbling block preventing large QCD simulations, including 
quarks, is the presence of the determinant $\det(\not\!\!D +m)$ in the 
partition function.  There have been many proposals for dealing with the 
determinant.  The first algorithms tried to compute the change in the 
determinant when a single link variable was updated \cite{Weingarten:81a}.  
This turned out to be prohibitively expensive.  So instead, the approximate 
method of pseudo-fermions \cite{Fucito:81a} was used.  
Today, however, the preferred approach is the so-called ``Hybrid Monte 
Carlo'' algorithm \cite{Duane:87a} which is exact.  
The basic idea is to invent some dynamics for the variables in the system in 
order to evolve the whole system forward in (simulation) time, and then do a 
Metropolis accept/reject for the entire evolution on the basis of the total 
energy change.  The great advantage is that the whole system is updated in 
one fell swoop.  The disadvantage is that if the dynamics are not correct, 
the acceptance will be very small.  Fortunately, (and this is one of very few 
fortuitous happenings where fermions are concerned) good dynamics can be 
found\@: the Hybrid algorithm \cite{Duane:85a}.  This is a neat combination 
of the deterministic microcanonical method \cite{Callaway:83a}, 
\cite{Polonyi:83a} and the stochastic Langevin method \cite{Parisi:81a}, 
\cite{Batrouni:85a}, which yields a quickly evolving, ergodic algorithm for 
both gauge fields and fermions.  The computational kernel of this algorithm 
is the repeated solution of systems of equations of the form,
\begin{equation}
(\not\!\!D +m) \phi = \eta,
\end{equation}
where $\phi$ and $\eta$ are vectors that live on the sites of the lattice.
To solve these equations, one typically uses a conjugate gradient or one of 
its cousins, since the fermion matrix $(\not\!\!D +m)$ is sparse.  For more 
details, see \cite{Gupta:88a}.  Such iterative matrix algorithms have as 
their basic component the $a \times b + c$ kernel, so again computers which 
do this efficiently will run QCD well.

\subsection{Concurrent QCD Machines}

Where is the computer power needed for QCD going to come from?  Today, the 
biggest resources of computer time for research are the conventional 
supercomputers at the NSF and DOE centers.  The centers are continually 
expanding their support for lattice gauge theory, but it may not be long 
before they are overtaken by several dedicated efforts involving concurrent 
computers.  It is a revealing fact that the development of most 
high-performance parallel computers---the Caltech Cosmic Cube, the Columbia 
Machine, IBM's GF11, APE in Rome, the Fermilab Machine---was actually 
motivated by the desire to simulate lattice QCD.  Geoffrey~Fox and
Chuck~Seitz, at Caltech, built the first hypercube computer, the Cosmic Cube 
or Mark~I, in 1983.  It had 64 nodes, each of which was an 
Intel 8086/87 microprocessor with $\rm 128\,KB$ of memory, giving a total of 
about $\rm 2\,Mflops$ (measured for QCD).  This was quickly upgraded to the
Mark~II hypercube with faster chips, twice the memory per 
node, and twice the number of nodes in 1984.  Now QCD is running on the 
latest Caltech hypercube, the 128-node Mark~IIIfp\index{Mark~IIIfp} (built by 
JPL), at $\rm 600\,Mflops$ sustained \cite{Ding:90b}.  Each node of the 
Mark~IIIfp hypercube contains two Motorola 68020 microprocessors, one for 
communication and the other for calculation, with the latter supplemented by 
one 68881 coprocessor and a 32-bit Weitek floating point processor.  
Norman~Christ and Anthony~Terrano, at Columbia, built their first parallel 
computer for doing lattice QCD calculations in 1984.  It had 16 nodes, each 
of which was an Intel 80286/87 microprocessor, plus a TRW 22-bit floating 
point processor with $\rm 1\,MB$ of memory, giving a total peak performance 
of $\rm 256\,Mflops$.  This was improved in 1987 using Weitek rather than TRW 
chips so that 64 nodes give $\rm 1\,Gflops$ peak.  Recently, at Columbia, 
they have finished building their third machine\@: a 256-node $\rm 
16\,Gflops$ lattice QCD computer.  Don~Weingarton initiated the GF11
project in 1984 at IBM.  The GF11 is a SIMD machine comprising 576 Weitek 
floating point processors, each performing at $\rm 20\,Mflops$ to give the 
total $\rm 11\,Gflops$ peak implied by the name.  First preliminary results 
from this project were report at Lattice `91.  The APE (Array Processor with 
Emulator) computer is basically a collection of 3081/E processors (which were 
developed by CERN and SLAC for use in high energy experimental physics) with 
Weitek floating point processors attached.  However, these floating point 
processors are attached in a special way---each node has four multipliers and 
four adders, in order to optimize the $a \times b + c$ calculations, which 
form the major component of all lattice QCD programs.  This means that each 
node has a peak performance of $\rm 64\,Mflops$.  The first small 
machine---Apetto---was completed in 1986 and had four nodes yielding a peak 
performance of $\rm 256\,Mflops$.  Currently, they have a second generation 
of this machine with $\rm 1\,Gflops$ peak from 16 nodes.  By 1992, the APE 
collaboration hopes to have completed the $\rm 100\,Gflops$ 4096-node 
``Apecento'', based on specialized VLSI chips that are software compatible 
with the original APE\@.  Not to be outdone, Fermilab is also using its 
high-energy experimental physics emulators in constructing a lattice QCD 
machine for 1991 with 256 of them arranged as a $2^5$ hypercube of crates, 
with eight nodes communicating through a crossbar in each crate.  All 
together, they expect to get $\rm 5\,Gflops$ peak performance.  These 
performance figures are summarized in Table~\ref{tab2.qcd}.  (The ``real'' 
performances are the actual performances obtained on QCD codes.)

\begin{table}
\begin{center}
\begin{tabular}{|l|l|r|r|}
\hline
\multicolumn{1}{|c|}{Computer}&\multicolumn{1}{c|}{Year}&
\multicolumn{1}{c|}{Peak}&\multicolumn{1}{c|}{Real}
\\
\hline
\medStrut
Caltech I    & 1983 & 3 & 2" \\
Caltech II   & 1984 & 9 & 6" \\
Caltech III  & 1989 & 2000 & 600" \\
Columbia I   & 1984 & 256 & 20" \\
Columbia II  & 1987 & 1000 & 200" \\
Columbia III & 1990 & 16000 & 6000" \\
IBM GF11     & 1990 & 11000 & 10000* \\
APE I        & 1986 & 256 & 20" \\
APE II       & 1988 & 1000 & 200" \\
APE III      & 1992 & 100000 & 20000* \\
Fermilab     & 1991 & 5000 & 1200* \\
\hline
\end{tabular}
\end{center}
*All real times are measured except these predicted ones.
\caption{Peak and Real Performances in Mflops of ``Homebrew'' QCD
Machines}
\label{tab2.qcd}
\end{table}

Major calculations have also been performed on commercial SIMD machines, 
first on the ICL Distributed Array Processor (DAP) at Edinburgh University 
during the period 1982--1987, and now on the TMC Connection Machine (CM-2); 
and on commercial distributed memory MIMD machines like the nCUBE hypercube 
at Caltech.  In fact, currently, the Connection Machine is the most powerful 
commercial QCD machine available\@:  running full QCD at a sustained rate of 
approximately $\rm 5\,Gflops$ and $\rm >5\,Gflops$ after optimization, on a 
$\rm 64\,K$ CM-2 \cite{Baillie:89e}, \cite{Brickner:91b}.

It is interesting to note that when the various groups began building their
``homebrew'' QCD machines, it was clear they would out-perform all
commercial (traditional) supercomputers; however, now that commercial 
parallel supercomputers have come of age \cite{Fox:89n} the situation is not 
so obvious.  To emphasize this, we describe QCD calculations on both the 
home-grown Caltech hypercube and on the commercially available Connection 
Machine.

\subsection{QCD on the Caltech Hypercubes}

To make good use of MIMD distributed memory machines like hypercubes, one 
should employ domain decomposition.  That is, the domain of the problem 
should be divided into subdomains of equal size, one for each processor in 
the hypercube; and communication routines should be written to take care of 
data transfer across the processor boundaries.  Thus, for a lattice 
calculation, the $N$ sites are distributed among the $P = 2^d$ processors 
using a decomposition, which ensures that processors assigned to adjacent 
subdomains are directly linked by a communication channel in the hypercube 
topology.  Each processor then independently works through its subdomain of 
$n = N/P$ sites updating each one in turn, only communicating with 
neighboring processors when doing boundary sites.  This communication 
enforces ``loose synchronization,'' which stops any one processor from racing 
ahead of the others.  Load balancing is achieved with equal-size domains.  If 
the nodes contain at least two sites of the lattice, all the nodes can 
update in parallel, satisfying detail balance, since loose synchronicity 
guarantees that all nodes will be doing black, then red sites alternately.
The characteristic timescale of the communication, $t_{comm}$, corresponds to 
roughly the time taken to transfer a single $SU(3)$ matrix from one node to 
its neighbor.  Similarly, we can characterize the calculational part of the 
algorithm by a timescale, $t_{calc}$, which is roughly the time taken to 
multiply together two $SU(3)$ matrices.  For all hypercubes built {\em 
without} floating point accelerator chips $t_{comm} << t_{calc}$ and, hence, 
QCD simulations are extremely ``efficient,'' where efficiency is defined by 
the relation
\begin{equation}
\epsilon  =  \frac{T_1}{k\, T_k},
\end{equation}
where $T_k$ is the time taken for $k$ processors to perform the given 
calculation.  Typically, such calculations have efficiencies in the range 
$\epsilon \geq .90$, which means they are ideally suited to this type of 
computation since doubling the number of processors nearly halves the 
total computational time required for solution.  However, as we shall see 
(for the Mark~IIIfp hypercube, for example), the picture changes dramatically 
when fast floating point chips are used; then $t_{comm} \simeq t_{calc}$ and 
one must take some care in coding to obtain maximum performance.

Rather than describe every calculation done on the Caltech hypercubes,
we shall concentrate on one calculation that has been done several times
as the machine evolved---the heavy quark potential\index{quark potential} 
calculation (``heavy'' because the quenched approximation is used).

QCD provides an explanation of why quarks are confined inside hadrons in that
lattice calculations reveal that the inter-quark potential rises linearly as
the separation between the quarks increases.  Thus, the attractive force
between quarks grows linearly with the separation.  This force, called the 
``string tension'', is carried by the gluons, which form a ``string'' between 
the quarks.  On the other hand, at short distances, quarks and gluons are 
``asymptotically free'' and behave like electrons and photons, interacting 
via a Coulomb-like force.  Thus, the quark potential $V$ is written as
\begin{equation}
V(R) = -\frac{\alpha}{R} + \sigma R,
\end{equation}
where $R$ is the separation of the quarks, $\alpha$ is the coefficient of the 
Coulombic potential and $\sigma$ is the string tension.  Experimentally, in 
fitting charmonium data to this Coulomb plus linear potential, Eichten
{\em et al.} \cite{Eichten:80a} estimated that $\alpha=0.52$ and 
$\sigma$=0.18GeV$^2$.  Thus, a goal of the lattice calculations is to 
reproduce these numbers.  Enroute to this goal, it is necessary to show that 
the numbers from the lattice are ``scaling'', that is, if one calculates a 
physical observable on lattices with different spacings then one gets the 
same answer.  This means that the artifacts due to the finiteness of the 
lattice spacing have disappeared and continuum physics can be extracted.

The first heavy quark potential calculation using a Caltech hypercube was 
performed on the 64-node Mark~I in 1984 on a $12^3\times 16$ lattice with 
$\beta$ ranging from $5.8$ to $7.6$ \cite{Otto:84a}.  $\alpha$ was found to 
be $0.25$ and the string tension (converting to the dimensionless ratio) 
$\sqrt{\sigma}/\Lambda = 106$.  The numbers are quite a bit off from the 
charmonium data but the string tension did appear to be scaling, albeit in 
the narrow window $6.0 \le \beta \le 6.4$.

The next time around, in 1986, the 128-node Mark~II hypercube was used on a 
$20^4$ lattice with $\beta = 6.1, 6.3, 6.5, 6.7$ \cite{Flower:86b}.  The 
dimensionless string tension decreased somewhat to 83, but clear violations 
of scaling were observed\@: the lattice was still too coarse to see continuum 
physics.

Therefore, the latest (1989) calculation using the Caltech/JPL 32-node 
Mark~IIIfp hypercube concentrated on one $\beta$ value, $6.0$, and 
investigated different lattice sizes\@:  $24^3\times 10$, $24^3\times 12$, 
$12^3\times 24$, $16^3\times 24$ \cite{Ding:90b}.  Scaling was not 
investigated, however, the values of $\alpha=0.58$ and $\sqrt{\sigma}/\Lambda 
= 77$, {\em i.e.,\/} $\rm \sigma$=0.15\,GeV$^2$ compare favorably with the 
charmonium data.  This work is based on about 1300~CPU hours on the 32-node
Mark~IIIfp hypercube, which has a performance of roughly twice a CRAY~X-MP 
processor.  The whole 128-node machine performs at $\rm 0.6\,Gflops$.  As 
each node runs at $\rm 6\,Mflops$, this corresponds to a speedup of 100, and 
hence, an efficiency of 78\%.  These figures are for the most highly 
optimized code.  The original version of the code written in C ran on the 
Motorola chips at $\rm 0.085\,Mflops$ and on the Weitek chips at $\rm 
1.4\,Mflops$.  The communication time, which is roughly the same for both, is 
less than a 2\% overhead for the former but nearly 30\% for the latter.  When 
the computationally intensive parts of the calculation are written in 
assembly code for the Weitek, this overhead becomes almost 50\%.  This $\rm 
0.9\,msec$ of communication, shown in lines two and three in 
Table~\ref{tab3.qcd}, is dominated by the hardware/software message startup 
overhead (latency), because for the Mark~IIIfp the node to node communication 
time, $t_{comm}$, is given by
\begin{equation}
t_{comm} \simeq ( 150+2*W) \quad \mu\, {\rm sec},
\end{equation}
where $W$ is the number of words transmitted.  To speed up the communication, 
we update all even (or odd) links (eight in our case) in each node, allowing
us to transfer eight matrix products at a time, instead of just sending one 
in each message.  This reduces the $\rm 0.9\,msec$ by a factor  of
\begin{equation}
\frac{8*(150+18*2)}{ 150+8*18*2} = 3.4
\end{equation}
to $\rm 0.26\,msec$.  On all hypercubes with fast floating point chips---and 
on most hypercubes without these chips for less computationally intensive 
codes---such vectorization of communication is often important.  In 
Figure~\ref{fig3.qcd}, the speedups for many different total lattice sizes 
are shown.  For the largest lattice size, the speedup is 100 on the 128-node. 
The speedup is almost linear in number of nodes.  As the total lattice volume 
increases, the speedup increases, because the ratio of 
calculation/communication increases.  For more information on this 
performance analysis, see \cite{Ding:90c}.

\begin{table}
\begin{center}
\begin{tabular}{|l|c|c|c|c|}
\hline
\multicolumn{1}{|c|}{Programming Level}&\multicolumn{1}{c|}{Calc.}&
\multicolumn{1}{c|}{Comm.}&\multicolumn{1}{c|}{Total}&
\multicolumn{1}{c|}{Mflops} \\
&\multicolumn{1}{c|}{Time}&\multicolumn{1}{c|}{Time}&
\multicolumn{1}{c|}{Time}& \\
\hline
\medStrut
1.\ \ Motorola MC68020/68881 in C & 52\phantom{.}"" & 0.86 & 53\phantom{.}"" &
0.085 \\
2.\ \ Weitek XL all in C             & "2.25 & 0.90 & "3.15 & 1.4"" \\
3.\ \ Weitek XL parts in Assembly    & "0.94 & 0.90 & "1.84 & 2.4"" \\
4.\ \ Weitek XL Assembly, vec. comm. & "0.94 & 0.26 & "1.20 & 3.8"" \\
5.\ \ Weitek XL Assembly, no comm.   & "0.94 & 0.0" & "0.94 & 4.8"" \\
\hline
\end{tabular}
\caption{Link Update Time (msec) on Mark~IIIfp Node for Various Levels
of Programming}
\label{tab3.qcd}
\end{center}
\end{table}

\begin{figure}
\vspace*{2.0in}
\caption{Speedups for QCD on the Mark~IIIfp}
\label{fig3.qcd}
\end{figure}

\subsection{QCD on the Connection Machine}

The Connection Machine Model CM-2 is also very 
well-suited for large-scale simulations of QCD.  The CM-2 is a 
distributed-memory, single-instruction multiple-data (SIMD) massively 
parallel processor comprising up to 65536 ($\rm 64\,K$) processors 
\cite{Hillis:85a} \cite{Hillis:87a}.  Each processor consists of an 
arithmetic-logic unit (ALU), eight or $\rm 32\,Kbytes$ of random-access 
memory (RAM) and a router interface to perform communications among the 
processors.  There are 16 processors and a router per custom VLSI chip, with 
the chips interconnected as a 12-dimensional hypercube.  Communications
among processors within a chip work essentially like a crossbar interconnect.
The router can do general communications, but we require only local ones for 
QCD so we use the fast nearest-neighbor communication software called NEWS\@.
The processors deal with one bit at a time.  Therefore, the ALU can compute 
any two Boolean functions as output from three inputs, and all data paths are
one-bit wide.  In the current version of the Connection Machine (the CM-2),
groups of 32 processors (two chips) share a 32-bit (or 64-bit) Weitek 
floating point chip, and a transposer chip, which changes 32 bits stored 
bit-serially within 32 processors into 32 32-bit words for the Weitek, and 
vice versa.

The high-level languages on the CM, such as *Lisp and CM-Fortran, compile into
an assembly language called Parallel Instruction Set (Paris).  Paris regards 
the $\rm 64\,K$ bit-serial processors as the fundamental units in the 
machine.  However, floating-point computations are not very efficient in the 
Paris model.  This is because in Paris, 32-bit floating point numbers are 
stored ``field-wise''; that is, successive bits of the word are stored at 
successive memory locations of each processor's memory.  However, 32 
processors share one Weitek chip, which deals with words stored 
``slice-wise''---stored across the processors, one bit in each.  Therefore, 
to do a floating-point operation, Paris loads in the field-wise operands, 
transposes them slice-wise for the Weitek (using the transposer chip), does 
the operation, and transposes the slice-wise result back to field-wise for 
memory storage.  Moreover, every operation in Paris is an atomic process; 
that is, two operands are brought from memory and one result is stored back 
to memory, so no use is made of the Weitek registers for intermediate 
results.  Hence, to improve the performance of the Weiteks, a new assembly 
language called CM Instruction Set (CMIS) has been written, which models the 
local architectural features much better.  In fact, CMIS ignores the 
bit-serial processors and thinks of the machine in terms of the Weitek chips. 
Thus, data can be stored slice-wise, eliminating all the transposing back and 
forth.  CMIS allows effective use of the Weitek registers, creating a memory 
hierarchy, which combined with the internal buses of the Weiteks, offers 
increased bandwidth for data motion.

When the arithmetic part of the program is rewritten in CMIS (just like
on the Mark~IIIfp when this was rewritten in assembly code) the
communications become a bottleneck.  Therefore, we need also to speed
up the communication part of the code.  On the CM-2, this is done using
the ``bi-directional multi-wire NEWS'' system.  As explained above, the
CM chips (each containing 16 processors) are interconnected in a
12-dimensional hypercube.  However, since there are two CM chips for
each Weitek floating-point chip, the floating-point hardware is
effectively wired together as an 11-dimensional hypercube, with two
wires in each direction.  This makes it feasible to do simultaneous
communications in both directions of all four space-time directions in
QCD---bi-directional multi-wire NEWS---thereby reducing the
communication time by a factor of eight.  Moreover, the data
rearrangement necessary to make use of this multi-wire NEWS further
speeds up the CMIS part of the code by a factor of two.

Currently, the Connection Machine is the most powerful commercial QCD machine 
available\@: the ``Los Alamos collaboration'' is running full QCD at a 
sustained rate of $\rm 5.3\,Gflops$ on a $\rm 64\,K$ CM-2.  As was the 
case for the Mark~IIIfp hypercube, in order to obtain this performance, one 
must resort to writing assembly code for the Weitek chips {\em and\/} for the 
communication.  Our original code, written entirely in the CM-2 version
of *Lisp, achieved around $\rm 1\,Gflop$ \cite{Baillie:89e}.  As shown in 
Table~\ref{tab4.qcd}, this code spends 34\% doing communication.  When we 
rewrote the most computationally intensive part in the assembly language 
CMIS, this rose to 54\%.  Then when we also made use of ``multi-wire
NEWS'' (to reduce the communication time by a factor of eight) it fell
to 30\%.

\begin{table}
\begin{center}
\begin{tabular}{|l|c|c|c|r|}
\hline
\multicolumn{1}{|c|}{Programming Level}&\multicolumn{1}{c|}{Calc. time}&
\multicolumn{1}{c|}{Comm. time}&\multicolumn{1}{c|}{Total Time}&
\multicolumn{1}{c|}{Mflops}
\\
\hline
\medStrut
All in *Lisp & 8.7 & 4.5 & 13.2 & 900 \\
Inner Loop in CMIS & 3.3 & 3.9 & "7.2 & 1600 \\
Multi-wire CMIS & 1.7 & 0.5 & "2.2 & 5300
\\
\hline
\end{tabular}
\caption{Fermion Update Time (sec) on $\rm 64\,K$ Connection Machine for
Various Levels of Programming}
\label{tab4.qcd}
\end{center}
\end{table}

\subsection{Status and Prospects}

The status of lattice QCD may be summed up as\@: underway.  Already there 
have been some nice results in the context of the quenched approximation, but 
the lattices are still too coarse and too small to give definitive results.
Results for full QCD are going to take orders of magnitude more computer 
time, but we now have an algorithm---Hybrid Monte Carlo---which puts real 
simulations within reach.

When will the computer power be sufficient?  In Figure~\ref{fig4.qcd}, we 
plot the horsepower of various QCD machines as a function of the year they 
started to produce physics results.  The performance plotted in this case is 
the real sustained rate on actual QCD codes.  The surprising fact is that the 
rate of increase is very close to exponential, yielding a factor of ten every 
two years!  On the same plot, we show our estimate of the computer power 
needed to re-do correct quenched calculations on a $128^4$ lattice.  This 
estimate is also a function of time, due to algorithm improvements.

\begin{figure}
\vspace*{2.0in}
\caption{Megaflops for QCD Calculations}
\label{fig4.qcd}
\end{figure}

Extrapolating both lines, we see the outlook for lattice QCD is rather bright.
Reasonable results for the phenomenologically interesting physical 
observables should be available within the quenched approximation in the 
mid-90's.  With the same computer power, we will be able to re-do today's 
quenched calculations using dynamic fermions (but still on today's size of 
lattice).  This will tell us how reliable the quenched approximation is.  
Finally, results for the full theory with dynamical fermions on a $128^4$ 
lattice should follow early in the next century (!), when computers are two 
or three orders of magnitude more powerful.

\end{document}

