\documentstyle[11pt]{article}
\def\lea{\mathrel{\raise .4ex\hbox{\rlap{$<$}\lower 1.2ex\hbox{$\sim$}}}}
\def\gea{\mathrel{\raise .4ex\hbox{\rlap{$>$}\lower 1.2ex\hbox{$\sim$}}}}
\newcommand{\medStrut}{{\rule{0cm}{4.0ex}}}
\newdimen\digitwidth
\setbox1=\hbox{0}
\digitwidth=\wd1
\catcode`"=\active
\def"{\kern\digitwidth}

\title{Parallel Computing Comes of Age: Supercomputer Level Parallel
Computations at Caltech}

\author{Geoffrey C. Fox \\
	Caltech Concurrent Computation Program \\
	Mail Code 206-49 \\
	Pasadena, California  91125}
	
\begin{document}
\maketitle

\begin{abstract}
Parallel supercomputers are now in regular use at Caltech for several
major scientific calculations.  We use this experience to abstract a set
of lessons for applications, decomposition, performance, hardware and
software.  We consider hypercubes, transputer arrays, networks of work
stations, and the SIMD Connection Machine CM-2 and AMT DAP.  These are
contrasted, where possible, with CRAY and other high performance
conventional computers.  Applications covered are lattice gauge theory,
string dynamics, statistical and condensed matter physics, plasma
physics, astronomical data analysis, quantum chemistry, graphics ray
tracing, grain dynamics, astrophysical particle dynamics, computer chess,
and Kalman filters.
\end{abstract}

\section{Introduction}

The heart of this article is a review of twelve applications which have
used parallel computers recently and achieved supercomputer level
performance.  we have operationally defined this being equivalent to
calculations involving a few hundred or more CRAY-XMP hours.  This
review is contained in Sec.~3 to Sec.~4 and is followed by a set of
lessons that we have learned from these and other experiences using
concurrent computers at Caltech.  Hopefully, our experience will be
useful to designers of hardware and software for parallel computers, as
well as potential users of such systems.

In previous reviews \cite{Fox:85c}, \cite{Fox:87c}, \cite{Fox:88b}, and
\cite{Fox:88b}, we have described the broad range of applications and
algorithms that successfully parallelize.  However, one of the most
exciting developments of the last year is the emergence of several
parallel machines that have comparable or greater performance than a
CRAY.  The hardware used in the applications of Secs.~3 to 14 is detailed
in Table~1.  There are other interesting machines, but table~1
represents those available at Caltech.  The accompanying performance
evaluation of \cite{Messina:89a} has implemented some of our
applications on a broad range of machines.  We will refer to this paper
when appropriate.

\begin{table}
\caption{Parallel Supercomputers Used in Secs.~3--14}
\vspace{7pt}
\begin{center}
\begin{tabular}{||p{1.00in}| p{1.75in}| p{2.50in}||}
\hline
\multicolumn{1}{||c|}{Machine}&\multicolumn{1}{c|}{Configuration}
&\multicolumn{1}{c||}{Key Characteristics}\\
\hline
\medStrut
NCUBE hypercube & 1024 nodes --- SANDIA & Scalar nodes with about 0.1 
megaflop per node\\
& "576 nodes --- Caltech &\\
\hline
Transputer Array & 32 nodes --- Caltech & Scalar nodes, each with about 4
times performance of NCUBE node\\
(MEIKO$\ldots$) & Large system coming into use at Edinburgh &\\
\hline
Mark IIIfp hypercube & 128 nodes --- Caltech (mainly used in 32 node
chunks so far)& Each node is a (short vector) pipelined FPU.  1--2
megaflops with rather disappointing compiler.  5--8 megaflops/node in
assembly language (easier than microcode used in previous WEITEK chip
sets).\\
\hline
Connection & 64K nodes --- Los Alamos & 64K single bit node system\\
CM-2 & 16K nodes --- ANL/Caltech & peaks out at about 1 gigaflop.\\
& 64K 1 bit processors is really 2048 32 bit processor for floating
point work & Probably will improve.\\
\hline
AMT DAP 510 & 1024 nodes --- ANL & Mesh of single bit processors ---
each faster than those of CM-2, but fewer of them\\
\hline
\end{tabular}
\end{center}
\end{table}

In Sec.~2, we describe data parallelism or decomposition as a universal
source of parallelism.  We introduce a classification of problems by
their temporal or computational structure \cite{Fox:88b}, which will be
helpful in succeeding sections in understanding which types of problems
perform well on which architectures.

\section{A Classification of Problems}
 
\subsection{Domain Decomposition}

All the applications considered here and essentially all other
successful uses of highly parallel machines have obtained their
parallelism from ``data parallelism'' \cite{Hillis:86a} or equivalently
``domain decomposition'' \cite{Fox:88a}.  This implies that a
computation can be considered as an algorithm applied to a data set.
Concurrency is obtained by acting simultaneously on different parts of
the data set.  This data set can be a collection of grid points, galaxies,
possible chess moves, neurons, database records, and so on.  Whatever
the architecture of the computer, parallelism comes from evolving these
entities simultaneously.  The three major commercial architectures
obtain their parallelism from this source, but implement it differently.
\begin{itemize}
\item Multicomputer (MIMD, distributed memory):\ Data and processing
distributed,
\item Shared Memory MIMD:\ Processing only distributed, and the user
does not need to make explicit the underlying data parallelism,
\item Distributed Memory SIMD:\ Data and synchronous processing
distributed.
\end{itemize}
Radical architectures such as dataflow may exploit other sources of
parallelism, but this has yet to be proved practical.

We should say that some functional parallelism has also been exploited
as in battle simulation of Sec.~14.  However, this does not lead to the
essentially unlimited parallelism typical of large data sets.  It is a
useful way of augmenting domain decomposition.

How easily, and under what circumstances, can one simultaneously evolve
different parts of the data set and so use domain decomposition to achieve
parallelism?  We can extend the previous idea of a data domain into a
physical, or space-time, picture of applications \cite{Fox:88b}, 
\cite{Fox:88kk}, \cite{Fox:88tt}, and \cite{Fox:88uu}.  This allows a 
classification of the ``architecture'' of problems and a quantification of the 
matching of problems and computer architectures.

\subsection{Space-Time Picture}

We will term the underlying data set, as discussed above, as the
``space'' associated with the problem.  Any implementation (simulation)
of this problem on a computer involves work or computation performed on
the elements of this ``space''.  this computation can be labeled with a
(discrete in practice) variable which we call ``time''.  Let us consider
some examples.

\begin{description}
\item[(a)] In the simulation of a three dimensional physical system, the
data domain or ``space'' is a (discretized version of) conventional
physical space.  ``Time'' is again just discretized time.  We have, of
course, designed our nomenclature so that it reduces to the conventional
one for such physical simulations.

\item[(b)] In say the calculation of the distribution of electromagnetic
waves in a radar antenna problem, we would step through time solving
Maxwell's wave equation.  This can be contrasted with iterative
solutions to elliptic equations which could involve a similar spatial
operator (say, $\nabla^{2}$).  We can consider successive
over-relaxation or Gauss Seidel iterative methods for solving the sparse
matrix problems.  In these cases, the ``space'' is again a conventional
three dimensional space but now iteration count plays the role of
``time''.

\item[(c)] Finally, consider full $M\times M$ matrix problems ---
inversion, multiplication, eigenvalues, etc..  Now the ``space'' is more
abstract being the (two dimensional array formed by the) set of $M^{2}$
matrix elements.  A typical algorithm is LU decomposition where one
successively eliminates rows and columns from the matrix; ``time'' is
naturally defined as the label ($i=1\ldots M$) of the eliminated
variable.
\end{description}

In Figure~1, we illustrate the resultant space-time-lattice which we
associate with a given problem.  A computation can be considered as a
{\em complex system\/} defined by a lattice in space-time
\cite{Fox:86a}.  The algorithm defines an interconnect between the
members of the data set and this leads to a {\em computational graph\/}
which we define by slicing the complex system at a fixed time.  In
\cite{Fox:87c}, \cite{Fox:88f}, and \cite{Fox:88kk}, we show that this
picture leads a string or world line formalism to describe such
systems.  Here, we discuss the time or temporal structure of the problem
and its relationships to the architecture of the underlying parallel
computer \cite{Fox:88b}.

\subsection{Temporal Classification of Problems}

We divide the temporal structure of a problem into three general classes
which we describe in order of increasing generality, and illustrate in
Figure~2.

{\em Synchronous\/} problems have a computational graph that is
independent of time even on a microscopic time scale.  One gets such a
structure when nature is viewed microscopically in terms of a set of
similar fundamental quantities evolved with identical algorithms.  A
good example is the QCD simulations described in Sec.~3.

{\em Loosely Synchronous\/} problems only show a temporally constant
computation graph when this is averaged over a macroscopic time
interval.  In between macroscopic synchronization points, the different
parts of the space evolve independently with no microscopic
synchronization.  This is seen clearly in spatially irregular physical
simulations where we step in time from $t_{0}$ to $t_{0}+\delta,
t_{0}+2\delta\ldots$.  These times provide the macroscopic
synchronization but if the problem is spatially irregular, quite
different computations could be involved for each data point in the time
evolution.  This temporal structure if natural when one models nature
macroscopically in terms of nonfundamental (and correspondingly
different in space and time) objects.  An example would be a
heterogeneous circuit simulation or the clustering algorithm of Sec.~7.

{\em Asynchronous\/} problems show no regularity in time and each
spatial member evolves in time with no natural mutual synchronization.
An example shown in Figure~2(c) is an event driven simulation.

We will find that nearly all successful parallel implementations
correspond to {\em synchronous\/} or {\em loosely synchronous\/}
problems.  Microscopic or macroscopic algorithmic synchronization
naturally leads ot paralle speedup.

We will use the term {\em properly asynchronous\/} to denote {\em
asynchronous\/} problems that are not {\em loosely synchronous.
Properly loosely synchronous\/} will denote cases which are {\em loosely
synchronous\/} but not {\em synchronous\/}.  We will use the term {\em
loosely synchronous\/} to include both {\em synchronous\/} or {\em
properly loosely synchronous\/} cases.

We should emphasize that we are discussing the property of problems and
not the software implementation. For instance, most hypercubes have, as
standard, {\em asynchronous\/} communication systems even though these
are typically used for {\em loosely synchronous\/} problems.  However,
OCCAM on the Transputer and CrOS~III on the Caltech hypercubes are {\em
synchronous\/} problems, but it's greater generality naturally supports
debugging and other user environment features that lead to the whole
implementation being {\em asynchronous\/}.

This classification is further explored in \cite{Fox:88b} while Figure~3
summarizes our current lessons on the correspondence between computer
and problem architectures.  At a detailed level, this involves the
relationship between the spatial structure of the problem and the
topology of the parallel computer.  In particular, we introduce a
further class of {\em embarrassingly parallel\/} or {\em spatially
disconnected} problems.  In these, the basic entities in the
computational graph given by domain decomposition are computationally
independent.  These give excellent parallel performance for {\em
asynchronous\/} problems which when nontrivially spatially connected, as
in the example of Figure~2(c), provide major challenges.  The article on
time-warp [Weiland, this issue of CPandE], describes an important
approach to this difficult problem class.  The relevance of this
classification is illustrated in Table~2 which summarizes 84 parallel
applications or algorithms surveyed in \cite{Fox:88b}.  90\% of these
applications fell into the classes {\em synchronous, loosely
synchronous\/}, or {\em embarrassingly parallel\/}, which naturally give
parallelism scales to many nodes.

\begin{table}
\caption{Summary of 84 Separate Applications on Parallel Computers}
\vspace{7pt}
\begin{center}
\begin{tabular}{||r|p{2.5in}||}
\hline
\multicolumn{1}{||c|}{No.}&\multicolumn{1}{c||}{Application Field}\medStrut
\\
\hline
\medStrut
9  & Biology\\
4  & Chemistry \& Chemical Engineering\\
14 & Engineering\\
10 & Geology \& Earth/Space Science\\
13 & Physics\\
5  & Astronomy \& Astrophysics\\
11 & Computer Science\\
18 & Numerical Algorithm\\
\hline
\end{tabular}
\end{center}

\begin{center}
\begin{tabular}{||l|c|c|l||}
\hline
\multicolumn{1}{||l|}{Application Classification}
& \multicolumn{1}{c|}{No.}
& \multicolumn{1}{c|}{Fraction}& \multicolumn{1}{l||}{\ }\medStrut
\\
\hline
Synchronous (S)         & 34 & 0.40 & Total S+LS\\
\cline{1-3}
Loosely Synchronous (LS)  & 30 & 0.36 & 0.76\\
(not synchronous)       &    &      &     \\
\hline
Embarrassingly Parallel (EP) & "6 & 0.07 & Total S+LS+EP\\
-- runs on SIMD         &    &      & 0.90\\
-- needs MIMD           & "6 & 0.07 & Clear Scaling\\
\hline
Truly Asynchronous      & "8 & 0.10 & Unclear Scaling\\
\hline
\end{tabular}
\end{center}
\end{table}

\section{Lattice Monte Carlo Simulations}

\subsection{General Remarks}

Lattice theories represent one of the most computationally intense class
of problems \cite{Baillie:89b}.  They arise from numerical approaches to 
statistical physics or, in the path integral approach, to quantum field 
theories.  One has an effective partition function
\begin{equation}
Z = \int d\varphi_{i} \exp (- S [\varphi_{i}] )
\end{equation}
where even in today's modest problems, the integral running over the
fields $\varphi_{i}(x)$ can have over one million dimensions (degrees of
freedom).  Measurements or observables are then found from
\begin{equation}
\langle O\rangle\ = \int d\varphi_{i} O(\varphi_{i}) \exp (- S [\varphi_{i}] 
)/Z
\end{equation}
In the Monte Carlo method, one replaces the integral (2) by the sum
\begin{equation}
\langle O\rangle = \lim_{N \rightarrow \infty} \sum_{i} O(\varphi_{i})/N_{c}
\end{equation}
over $N_{c}$ configurations.  These configurations are generated
successively by making a series of small changes --- usually for single
sites \cite{Metropolis:53a}.  This ensures that one keeps configurations
$\varphi_{i}$ distributed according to the function
$\exp[-S(\varphi_{i}])$ which emphasizes the minute region of phase
space which is not exponentially suppressed.  For the currently hardest
calculations --- dynamical QCD with a $16^{4}$ lattice --- it takes
50--100 hours of time on CM-2 running at one gigaflop to produce a
statistically distinct (uncorrelated) configuration.  Very many such
configurations are needed in the averages eq.(3).  The slow evolution
produced by the successive application of a single site update, has
encouraged development of cluster update methods, but these are only
known at the present for the simpler theories.  This is a critical issue;
not only could clustering speed up the computation, it could alter the
necessary architecture as we will illustrate later.

In general, parallelism is straightforward for these problems --- one uses 
domain decomposition of the underlying  regular space.  There have been 
several algorithmic improvements recently --- which can usually be viewed as 
better importance sampling in the Monte Carlo integral.  The need to
continually update algorithms favors relatively general purpose machines with 
flexible high level software. This suggests to me, in the long run, that 
commercial parallel machines will be more successful than the many special 
purpose computers constructed within the high energy physics community. 

\subsection{Discrete Spin:  Ising and Potts Models (C. Baillie, P. 
Coddington)} 

These theories have actions $S$ given by
\begin{equation}
S = \sum_{\langle ij\rangle} \sigma_{i} \sigma_{j}
\end{equation}
over nearest neighbors $ij$ in a space of dimension $d$.  The spins
$\sigma_{i}$ are one bit (Ising) or several bits each (Potts).  These
problems are ideal for bit serial machines like the AMT DAP or
Connection Machine CM-2 --- although one cannot use the WEITEK floating
point units on the latter for these bit oriented problems.

We are currently using the AMT DAP 510 to investigate a $256^{3}$
three-dimensional Ising model.  It is conventional to rate machines by
the speed at which they generate new configurations --- measured in spin
updates/second.  The DAP~510 compares well with the fastest
supercomputers in at $0.6\times 10^{9}$ spin updates/sec compared
with the Hitachi S-820/80, which is only 40\% faster \cite{Ito:88a}.
However, this performance rating is misleading as our physics
calculation is dominated by the calculation by eq.(3) of observables.
In fact, we have only crudely coded the update stage on the DAP~510 so
that is a factor of 100 slower than the value quoted above.  Even this
slow update only consumes 1\% of the total time.

\subsection{Continuous Spin:  Two Dimensional {\bf $XY$} and {\bf $O(3)$} 
Models (R. Gupta (Los Alamos), J. Apostolakis, C. Baillie)}

These systems are given by actions with the form of eq.(4) with spins 
$\sigma_{i}$ that are $N\times N$ matrices of the group $O(N)$.  We finished, 
last summer, a major $O(2)$ (or $XY$) calculation using the novel 
``over-relaxation'' single site update algorithm \cite{Gupta:88a}.  This used 
the 128 node FPS T~Series hypercube at Los Alamos and realized two megaflops 
per node.  These calculations are being continued for the $O(3)$ and similar 
groups by my graduate student, John Apostolakis.  The high statistics of the 
$XY$ study allowed the refutation of the conclusions of a recent paper
\cite{Seiler:88a} and a confirmation of the theoretical predictions of 
Kosterlitz and Thouless \cite{Kosterlitz:73a}.  The T~Series  hypercube is a 
poorly designed machine and is only suitable for a small class of regular 
calculations; it has poor scalar compared to vector performance and slow 
communication.  We used a $256\times 256$ lattice allowing the other 
dimension to be viewed as a large vector (multiples of 128) achieving good 
performance from the inflexible WEITEK chip based node.  Some of our 
calculations have been marred by hardware glitches, requiring that one, for 
instance, avoid the hardware vector divide and long vector ($> 256$) 
instructions. 

Niedermayer and Wolff \cite{Niedermayer:88a}, \cite{Wolff:89a} have introduced effective 
clustering methods which currently we believe are unsuitable for the FPS 
architecture.  We expect to use the MEIKO transputer array to continue these 
calculations.  Note that the SIMD CM-2 would perform well on the algorithm we 
used for the T~Series, but we currently believe that a true MIMD architecture 
may be needed for the clustering calculation.  The regular vectors needed by 
the T~Series, make this old hypercube essentially SIMD in character; more 
precisely in the language of Sec.~2, it requires synchronous problems for 
good performance.  The currently known clustering algorithms are 
{\em properly loosely synchronous\/}.  We are currently experimenting with 
SIMD implementations which, although inefficient, may still have sufficiently 
good performance \cite{Baillie:89c}. 

\subsection{Pure Gauge QCD (S. Otto, J. Flower, H. Ding, C. Baillie)}

In this case, the action $S$ takes the form
\begin{equation}
S = \beta \sum_{\stackrel{\rm plaquettes}{p}} (1 - ReTrUp)
\end{equation}
where $Up$ is the product of SU(3) link matrices around the elementary 
plaquettes --- these are the eight $1\times 1$ loops including a given link 
which joins two sites.  This calculation can be both vectorized and 
parallelized even for the so called random block lattice \cite{Chiu:86a}, 
\cite{Chiu:88e}.  It can achieve good performance on essentially all 
architectures with MIMD, SIMD, or vector characteristics.  The large number of 
floating calculations in eq.(5) (explicitly 3347 for each link) dwarfs 
overheads such as communication. 

Even a modest calculation on a $16^{4}$ lattice requires 262,144 degrees of 
freedom, 19 megabytes of memory and of order $10^{15}$ floating point
operations.  This problem was originally tackled \cite{Otto:84a} on the first 
64 node Cosmic Cube Hypercube \cite{Seitz:85a} with a $12^{3}\times 16$ 
lattice and repeated on the 128 node Mark~II hypercube with a larger $20^{4}$ 
lattice.  These machines had peak performance below 5~mflop.  We have 
implemented this application on several more powerful machines with the
approximate performance given below. 

\begin{tabbing}
Connection Machine CM-2 (64K nodes)\ \ \= 900 mflops \kill
CRAY XMP (1 processor)\> 60 mflops\\
NCUBE (1024 nodes) \> 80 mflops\\
JPL Mark~IIIfp hypercube (128 nodes)\> 500 mflops\\
Connection Machine CM-2 (64K nodes)\> 900 mflops\\
\end{tabbing}

Table~3 gives more details of the performance of the Mark~IIIfp
hypercube on this application.

\begin{table}
\caption{Table 3}
\end{table}

In Figure~4, we show physics results found by Ding on the Mark~IIIfp 
hypercube on  $24^{3}\times 10$ lattice \cite{Ding:89a}.  This implementation 
involved 8,000 lines of C and WEITEK XL assembly code --- the commercial 
compiler for this pipelined chip set is poor.  The CM-2 results correspond to 
3,000 lines of *LISP written by Brickner \cite{Brickner:89a}.  One needs 
eight virtual processors (corresponding to eight separate calculations) to 
get good performance. 

The flexibility of parallel machines is illustrated by Chiu's calculations 
\cite{Chiu:88e}.  These use a random block lattice where, in the largest
simulation, a different $40^{4}$ lattice is calculated on each node of the 
1024 node NCUBE hypercube at SANDIA; an {\em embarrassingly parallel\/} 
application in the language of Sec.~2. 

\subsection{Dynamical Fermion QCD (C. Baillie, R. Brickner (Los Alamos), R. 
Gupta (Los Alamos), G. Kilcup (Brown), A. Patel (CERN), S. Sharpe
(Washington))}

This ``ultimate'' QCD calculation includes the quark (fermion) degrees
of freedom $\varphi$ at the lattice sites with action
\begin{equation}
S = - \varphi^{\dagger} (D + m)^{-1} \varphi
\end{equation}
Technically, one uses a Wilson fermion  formulation and this is dominated by 
inverting the $65536\times 65536$ (on a $16^{4}$ lattice) sparse fermion 
matrix $D + m$.  This is currently done with a conjugate gradient or minimal 
residue method.  The CRAY-2 uniprocessor code runs at 100 megaflops and the 
full size Connection Machine CM-2 has about one gigaflop performance
\cite{Baillie:89e}.  We expect this to improve as the CM-2 support software 
for floating point calculations is developed.  We anticipate this problem to 
need $10^{16}$ flops or $\sim 3,000$ CM-2 hours for even initial studies. 

Currently, there are no effective clustering algorithms known for either the 
pure gauge or fermion QCD problem.  The regular ``dumb'' algorithm used so 
far does not require a sophisticated architecture.  This could change when an
effective clustering method is introduced. 

DOE has awarded two collaborations large ($\sim 6,000$ hours) blocks of time 
for ``grand challenge'' QCD calculations on the CRAY-2 and ETA-10 for QCD 
calculations.  NSF has certainly already devoted much more time than this to 
this problem at the supercomputer centers.  Conventional and parallel 
supercomputers are compared in Table~4. With several months running on the 
CM-2 or Mark~IIIfp, parallel supercomputers are very competitive with the 
ETA-10 and CRAY-2.  The scalar node MIMD architecture of the NCUBE is not 
well optimized for this highly regular vectorizable problem.  A {\em
synchronous\/} algorithm is the classification of Sec.~2.

\newdimen\decwidth
\setbox2=\hbox{.}
\decwidth=\wd2
\catcode`@=\active
\def@{\kern\decwidth}

\begin{table}
\caption{Approximate Dynamical Fermion QCD Performance for the ``Grand
Challenge''}
\vspace{7pt}
\begin{center}
\begin{tabular}{||p{2.25in}|c|p{1.75in}||}
\hline
\multicolumn{1}{||c|}{Machine}&\multicolumn{1}{c|}{Performance}&
\multicolumn{1}{c||}{Time Allocated}
\\
&\multicolumn{1}{c|}{(megaflops)} &\multicolumn{1}{c||}{in 
1989 and comments}
\\
\hline
\medStrut
ETA-10 (1 processor)   & ""@350         & 1 year\\
CRAY-XMP (1 processor) & ""@100         & 1 year\\
Mark IIIfp (128 nodes) & ""@750         & \\
NCUBE (1024 nodes)     & ""@100         & Not competitive\\
CM-2 (64K nodes)       & "\llap{$>$}1,000 & \\
Several SIMD Coarse    & "1,000\rlap{$\rightarrow$} & limited by\\
Grain Special Purpose Computers & 10,000 & software\\
\hline
\end{tabular}
\end{center}
\end{table}

\section{Section 4}

\section{High {\bf $T_{c}\/$} Superconductivity (T. Barnes, C.
Kotchan (Toronto))}

A Toronto group has been using the Caltech NCUBE to study the Quantum 
Anisotropic Heisenberg Model \cite{Barnes:88a}, \cite{Barnes:88b}, and 
\cite{Barnes:89a}.  This problem is related to the systems of Sec.~3.2, 
except one is studying the dynamical and not the statistical properties of a 
spin system.  The Hamiltonian $H$ is given by a nearest neighbor sum $i,j$ 
over a two dimensional grid: 
\begin{equation} H = \sum_{\langle i,j\rangle} 
\left[ s^{z}_{i} s^{z}_{j} + g ( s^{x}_{i} s^{x}_{j} + s^{y}_{i} s^{y}_{j} ) 
\right] 
\end{equation} 
(The coupling ratio, $g=0$, is the Ising model, $g=\infty$ is the $XY$ model) 

This can be studied as a three dimensional lattice theory, but Barnes has 
developed an equivalent random walk approach for solving Schr\"{o}dinger's 
equation with an imaginary time $T$:
\begin{equation} 
-\frac{\partial\psi}{\partial T} = H\psi (T) 
\end{equation} 
Typically, an $8\times 8$ lattice is evolved separately on each node by the 
NCUBE with approximately $10^{6}$ independent evolutions needed. This 
algorithm is reminiscent of the (far more complex) neutron transport 
calculations studied at DOE laboratories.  This {\em embarrassingly
parallel\/} algorithm was easily implemented in C on the NCUBE.  The 256 node 
NCUBE hypercube achieved three times the performance of the original FORTRAN 
implementation for the CRAY XMP, which clearly vectorized poorly.  Current 
calculations correspond to several hundred hours of CRAY~XMP time.

\section{Exchange Energies in {\bf He$^{3}$} at a Temperature of
{\bf $0.1$~mK$^\circ$} (S. Callahan, M. Cross)} 

Callahan's Caltech condensed matter Ph.D. involved a Monte Carlo method to 
calculate exchange energies in solid He$^{3}$ \cite{Callahan:88a}, 
\cite{Callahan:88b}.  Rather modest systems were used with 54 to 128 
particles arranged in a three-dimensional spatial mesh which is further 
extended in time.  Use of the 512 node NCUBE required parallelism in several 
aspects of the problem.  The forces are not nearest neighbor and 
decomposition of their calculations over space leads to a factor of four in 
parallelism.  Decomposing time (direction of path) leads to another factor of 
16.  This 64 fold data parallelism is combined with $2\rightarrow 8$ 
independent runs ({\em i.e.,\/} decompose space of random configurations in 
the integral of eq.(2)). 

In an unfair comparison, the 64 node NCUBE has an efficiency of only 64\%, 
but outperforms the CRAY XMP by over a factor of two.  However, this used an 
early implementation of the C language for the CRAY which only realized a few 
megaflops on this problem \cite{Salmon:87a}.  Note more positively that our 
implementation used Salmon's CUBIX environment \cite{Fox:88a} allowing the 
identical code of about 2500 lines of C to run on either the NCUBE hypercube 
or CRAY! 

Callahan's thesis involved a total of about 250 hours computation on the
512 node NCUBE --- another supercomputer level calculation.

\section{Astrophysical Particle Dynamics (P. Quinn, J. Salmon, M.
Warren)}

N-body calculations have been revolutionized by a clustering technique
introduced by Appel \cite{Appel:85a} and developed significantly by
Barnes, Hut and Greengard \cite{Barnes:86a}, \cite{Greengard:87a}.  The
basic idea is simple; consider a cluster of $M$ stars for which we need
to calculate the interaction with a single star (far) outside the
cluster.  This straightforwardly requires $O(M)$ steps, but can gain a
factor of $M$ by ignoring the details of the cluster and just computing
with its center of mass.  As implemented by Barnes and Hut, one can
apply this idea recursively generating a tree (quad tree in two
dimensions) as illustrated in Figure~6 with, at most, one particle in
the cluster at the lowest level of tree \cite{Warren:88b},
\cite{Warren:88c}.  The naive calculation takes a time for each
simulation (time) step for a system of $N_{p}$ particles
\begin{equation}
T^{\,\mbox{naive}}_{\,\mbox{step}} = \frac{1}{2} N_{p}^{2} 
t_{2\mbox{particle}}
\end{equation}
while explicit implementation, shows that the clustering method takes
time
\begin{equation}
T^{\/\mbox{cluster}}_{\/\mbox{step}} = (20-50) N_{p}\log_{2}N 
t_{2\mbox{particle}}
\end{equation}
The cluster method has superior performance for $N_{p}\gea 1000$ particles.  
The current limit of O(10,000) particles for the $O(N_{p}^{2})$ algorithm is 
increased by an order of magnitude for the Barnes-Hut method.  The 
possibility of large $N_{p}$ of $O(10^{6})$ particles opens up several 
important astrophysical calculations including 
\begin{itemize} 
\item Study of the growth of fluctuations in the early universe 
\item Dynamics of globular clusters where one finds in nature of
$O(10^{5-6})$ stars and a difficult calculation as very short range 
interactions (binary stars) are critical 
\item Galactic structure and the collision of galaxies \cite{Hut:88a}. 
\end{itemize} 
We have just finished a calculation of the last type, illustrated in 
Figure~7, which used about 200 hours on the 512 node NCUBE \cite{Salmon:89a}.

This computation has several interesting features.  The cluster tree is
rebuilt each time step; a stage which is negligible in the sequential
version, but which appears to take of order 30\% of the concurrent
execution time in the $N=180,000$ particle simulation of \cite{Salmon:87a}.  
Initially, we found load imbalance, but this was solved by dynamically 
redistributing particles at each time step.  The information for this was 
found from the ``workload'' at the previously calculated time step.  As 
illustrated in Figure~6, orthogonal recursive bisection is used to 
distribute the particles.  Communication is required to fetch those parts
of the tree that are stored outside the node and will be needed for
updating particles within the node.  This ensures maximum re-use of the
communicated data and low communication overhead --- about 10\%.  This
approach does, however, use 75\% of the available NCUBE memory and limit the 
simulation size.  In spite of this, we are able to consider, on the current 
NCUBE, large problems that are difficult to implement on the limited memory 
CRAY-XMP and future such hypercubes with several megabytes of memory per node 
will allow much larger values of $N$.

The program was implemented with 3,400 lines of C code for the NCUBE, and some
timing information is given in Table~5 for a single time step.  This has
a comparison with FORTRAN CRAY code.  Figure~8 shows the measured
speedup on the NCUBE as a function of the number $N_{p}$ of particles.

\begin{table}
\caption{Performance on Astrophysical Particle Dynamics}
\vspace{7pt}
\begin{center}
\begin{tabular}{||l|c|c||}
\hline
& \multicolumn{1}{|c|}{$10^{4}$ Particles} &
\multicolumn{1}{c||}{$10^{5}$ Particles}\medStrut
\\
\hline
\medStrut
CRAY-XMP & 10 secs & 130 secs \\
Some Optimization & & \\
\hline
256 node NCUBE & & \\
\ \ --- time & 21 secs & 118 secs \\
\ \ --- efficiency & 24\% & 56\% \\
\hline
\end{tabular}
\end{center}
\end{table}

We see that the 256 node NCUBE outperforms the CRAY-XMP as reported in
\cite{Hernquist:87a} for $10^{5}$ particles, even though the hypercube 
efficiency is quite low.  The CRAY efficiency is even lower!

We can compare this for results from \cite{Messina:89a} on the $O(N_{p}^{2})$ 
particle dynamics algorithm for a vortex approach to fluid dynamics.  This 
vectorizes well on CRAY and runs with $>95$\% efficiency on the NCUBE; this 
time the CRAY-XMP is four times the performance of the 256 node NCUBE.

Considering SIMD architectures, we are not certain how to implement the 
clustering algorithm to get good performance on machines like the CM-2 
\cite{Hillis:87b}.  This is related to the difficulties we saw in Sec.~3.3 
with Monte Carlo clustering and in Sec.~4 with adaptive grids on SIMD 
machines.  Note that the ``regular'' clustering/multiscale algorithms such as 
multigrid or the FFT run quite well on SIMD machines; the difficulties in 
implementations for synchronous machines occurs for geometrically irregular 
cluster algorithms. This application is of the {\em properly loosely
synchronous\/} type.

\section{Astronomical Data Analysis (Anderson, Gorham, Kulkarni,
Prince)}

This group has pioneered the use of the NCUBE for astronomical data analysis 
\cite{Fox:88v}.  Our Caltech NCUBE system has a small (4 disks) parallel disk 
farm connected to the main hypercube; these disks are controlled by a SUN-4 
which also has additional peripherals, including the necessary tape drives. 

\subsection{Radio Astronomy}

In the most exciting work, radio data from the Arecibo radio telescope was
taken on December 26, 1988 in particularly advantageous circumstances in
that the holiday spirit reduced the ambient interference --- especially
from a nearby naval base.  Data is taken with a 0.5 millisecond time
interval and Fourier transformed (a large $2^{24}$ one-dimensional FFT)
to look for peaks corresponding to radio pulses from the rapidly rotating 
neutron star.  Two new pulsars have been discovered, using the NCUBE.
Both are located in the globular cluster M~15, making a total of three
known pulsars in this globular cluster.  The discovery of these pulsars
has prompted a reanalysis of current ideas concerning the origin of
neutron stars in globular clusters.

The computation involves both the FFT which is efficiently implemented
on the hypercube, and an I/O intensive stage taking a total time
comparable to the FFT.  This first I/O dominated stage is overlapped
with a calculation which corrects for frequency dispersion in the
interstellar medium.  The measured I/O performance of the system is
modest at 40~Kbytes/sec/drive including all overheads. Our ESMD disk
drives on the NCUBE are rated at a factor of 25 higher performance and
this suggests the need for better I/O software on the NCUBE.

The striking peak corresponding to the neutron star rotation period of 30.5 
milliseconds is shown in Figure~9 for the second pulsar discovered by the 
NCUBE.  This pulsar is part of a binary system and required a further compute 
intensive ``acceleration correction'' to remove the orbital effects of the 
binary system.  This additional computation has negligible I/O to disk but 
substantial internode communication. Processing 90 minutes of data taken at 
Arecibo takes about two hours for the dispersion correction and FFT stages 
(which discovered the first pulsar), while about 40 hours of 512 node NCUBE 
time were needed for the pulsar shown in figure~9.  The major I/O and large 
memory requirements of this calculation make it hard to compare the NCUBE 
with a CRAY or IBM~3090 performance. 

\subsection{Optical Astronomy}

Traditionally the resolution of ground based optical telescopes is
limited by atmospheric turbulence to about one second of arc (1/3600th
of a degree).  However, it is possible to eliminate the effects of turbulence 
using interferometry techniques familiar from radio astronomy, and
achieve resolution of 30 milliarcseconds.  Using the 200'' Mount Palomar 
telescope, one divides the total aperture into approximately one thousand 
15~cm disks.  The correlations between these disks are summarized in $10^{6}$ 
Fourier coefficients --- the bispectral function. These are averaged over 
many samples or frames lasting from 10--100~ms over which the turbulence is 
essentially constant.  This technique has been implemented on the NCUBE where 
a 20 minute Palomar exposure on an asteroid was analyzed for 10 hours on a 
256 node subcube and obtained the best image resolution ever seen for an 
asteroid which are illustrated in Figure~10; the improvement in resolution is 
about a factor of 30 over traditional methods in each linear dimension.
A similar analysis resolves binary stars \cite{Fox:89i}.

This novel method is still in its infancy; it illustrates graphically
how powerful computers can open up new approaches to scientific
problems.

\section{Quantum Chemistry Reaction Dynamics (Wu, Hipes, Kuppermann,
Cuccaro)}

\subsection{Overview}

Kuppermann's group has been developing a fundamental approach to the 
understanding of chemical reactions \cite{Cuccaro:88a}, \cite{Kuppermann:86a}.  
A goal is the description of reactions like $F + H_{2} \rightarrow FH + H$ 
which are the basis of an important chemical laser.  This is a difficult 
computation involving a factor of a thousand more computation than the 
prototypical initial example. 
\begin{equation} 
X + H_{2} (v,j,m) \rightarrow \\
\left\{ \begin{array}{ll} X + H_{2} (v', j', m') & \mbox{nonreactive} \\ H + 
XH (v'', j'', m'') & \mbox{reactive} \end{array} \right. 
\end{equation} 
where 

\begin{tabbing}
$v=0,1,2,\ldots$; \=  labels vibrational energy content\\*
$j=0,1,2,\ldots$; \>  labels rotational energy content\\*
$m=0, \pm 1, \pm 2,\ldots, \pm j$; \= labels spatial orientation\\*
\= of molecule\\*
\end{tabbing}
Operationally, this problem involves solving Schr\"{o}dinger's equation
$H\psi=E\psi$ for the wave function $\psi\/$ in the novel hyperspherical
basis set.  Here $H$, the Hamiltonian, is a second order linear
differential operation in six variables and $E$ is the total collision
energy.  This  computation breaks down into two phases which we describe 
separately below \cite{Hipes:88a}, \cite{Hipes:88b}.

\subsection{Phase I:\ \ Calculate Basis Functions and Matrices}

\begin{enumerate}
\item Construct a primitive basis set composed of product of analytic
and numerical functions --- each processor solves an independent
tridiagonal eigenvalue problem using bisection.
\item Evaluate 2D integrals using the primitive basis functions which
requires spline interpolation onto the quadrature grid.  Each processor 
calculates a subset of the matrix of such integrals.
\item Assemble the integrals into a real symmetric matrix.  Parallel
reduction by Householder transformations.  Redundant tridiagonal QR
algorithm in each processor to get eigenvalues/vectors.
\item Calculate more matrices of integrals using primitives and coefficients 
from step 3.  These matrix elements are distributed among the
processors.
\item I/O to store matrices
\item Repeat 51 times 
\end{enumerate}

Steps 1, 2, and 4 are trivially parallel and involve independent
computation.  Step 3 needs significant parallel algorithms adapted from work 
by Patterson at JPL \cite{Patterson:88a}.

This has been fully implemented for $H+H_{2}$ reactive scattering on the
32 node Mark~IIIfp hypercube with 3,000 lines of coding.  The total
runtime time was 9.8 hours of which 2.2 hours were I/O (stage 5).  This will
be improved soon with the high performance CIO (Concurrent I/O) hardware
on the Mark~IIIfp.  The same calculation on the SCS-40 (which has about 25\%
the performance of the CRAY XMP) took 71 hours.

This initial phase is followed by:

\subsection{Phase II:\ \ Integration of Coupled Linear Systems of Ordinary 
Differential Equations}

This phase uses Johnson's algorithm \cite{Johnson:73a} --- a fourth order 
special purpose integrator for chemistry simulations.  The resultant 
algorithm is dominated by matrix inversion (and not LU decomposition!) with 
some matrix multiplication.  I/O is also needed to initialize with the
results of Phase~I which determined the matrix elements.  This involved
2,000 lines of code and 74 coupled ordinary differential equations
integrated for 250 steps.  31 energies were calculated simultaneously to
reduce I/O overheads and the resultant calculation took two hours on the
Mark~IIIfp with an additional I/O overhead of 10\%.  A typical
transition probability curve is compared with the corresponding SCS-40 
calculations in Figure~11 and are benchmarked in \cite{Messina:89a}.
The vector machines do well on this calculation and these systematics
are abstracted in Table~6.

\begin{table}
\caption{Idealized Performance of High Performance Computers.  The listed 
numbers are approximate megaflops.}
\vspace{7pt}
\begin{center}
\begin{tabular}{||l|c|c|c|c||}
\hline
   &   &   &   \multicolumn{1}{c|}{Mark IIIfp} & \multicolumn{1}{c||}{NCUBE}\\
   & \multicolumn{1}{c|}{CM-2} & \multicolumn{1}{c|}{CRAY XMP} 
   & \multicolumn{1}{c|}{Hypercube} & \multicolumn{1}{c||}{Hypercube}\\
   & \multicolumn{1}{c|}{64K}  & \multicolumn{1}{c|}{1 processor} 
   & \multicolumn{1}{c|}{128 nodes} & \multicolumn{1}{c||}{1024 nodes}\\
   & \multicolumn{1}{c|}{``huge'' vector} 
   & \multicolumn{1}{c|}{``long'' vector} 
   & \multicolumn{1}{c|}{``short'' vector} 
   & \multicolumn{1}{c||}{scalar}\\
\hline
Super Regular & 3000 & 200 & 750 & 100\\
(e.g., large full matrix) &&     &     &\\
\hline
Typical Regular & 1000 & "75 & 500 & 100\\
(e.g., QCD)     &      &     &     &\\
\hline
Irregular & Fails & "10 & 100 & "50\\
(e.g., chess, clustering) & & & & \\
\hline
\end{tabular}
\end{center}
\end{table}

The NCUBE looks good compared to CRAY's on problems like those in Secs.~4, 5, 
6, 7, 8, 12, and 13 which are irregular and vectorize poorly.  Matrix algebra 
and, to a lesser extent, synchronous algorithms like QCD (Sec~3) show the 
vector machines at their best \cite{Messina:89a}.  This is not surprising --- 
they were designed for this problem class! Table~6 attempts to show a 
progression from SIMD  to long vector supercomputers to hypercubes with short 
vector nodes to scalar node MIMD machines.  This corresponds to increasingly 
general purpose machines.  As commented in Sec.~3.5, machines like the NCUBE 
are not competitive in peak performance on regular problems, but they offer 
good performance on range of problems and their speed degrades slowly as one 
increases the irregularity of the problem. 

\section{Plasma Physics (Liewer, Zimmerman (JPL), \mbox{Decyk}, Dawson 
(UCLA))}

Plasma Physics computations represent interesting challenges for distributed 
memory machines because the PIC or particle in the cell algorithm used 
involves two distinct decompositions.  Our example problem involves 
calculation of the orbits of plasma electrons in their own electromagnetic 
field as considered by a JPL-UCLA group on the JPL Mark~IIIfp hypercube 
\cite{Liewer:89a}, \cite{Liewer:88e}, \cite{Liewer:88b}.  In the first stage 
of the calculation, one finds the field using an FFT; this involves a 
decomposition with an equal number of mesh points on each node of the 
hypercube.  Then, one transforms to a separate decomposition shown in 
Figure~12 with equal number of particles in each node; this latter is the 
particle update or ``push'' part of the computation where the particle 
positions are evolved in the field.  Each stage can be efficiently 
implemented on the hypercube, but transforming between the two distinct 
decompositions must be done at each time step.  A general strategy for this 
has been discussed by Walker \cite{Walker:89a} using the 
{\em crystal\_accumulator} algorithm \cite{Fox:88a}.  This does not attempt 
to localize the calculation for each particle to a single node, but rather 
distributes it with calculations done ``on the fly'' as information is routed 
through the hypercube.  This method was originally developed by Furmanski for 
neural network simulations \cite{Furmanski:87a}, \cite{Fox:88g}.  It can only 
be implemented well on machines like the NCUBE and transputer arrays; the 
Mark~IIIfp and Ametek S2010, where calculation and communication subsystems 
are separated, do not support the {\em crystal\_accumulator\/} well.  If 
important, this algorithm requires communication subsystems that support a 
combination of messages in a similar fashion to fetch and add, as proposed, 
for combining networks in shared memory parallel computers 
\cite{Gottlieb:86a}. 

Currently, we only have measurements for the more straightforward
strategy where information is routed to the destination node and then
combined; rather than being combined en route.  This is almost certainly
the best algorithm for the Mark~III hypercube which has separate 
communication and calculation subsystems.  The combining overhead would be 
severe for Walker's approach on this and similar machines, where the
interface between communication and calculation on the node introduces a 
significant latency.  In Table~7, we compare the performance of (a) the push 
stage and (b) the total code for a variety of machines.  The 64~node 
Mark~IIIfp hypercube is about twice the performance of the CRAY XMP on the 
push stage, but only comparable for the total calculation.  This indicates 
that either the FFT or movement between decompositions is inefficient on the 
hypercube.  Note that both the CRAY-XMP and Mark~IIIfp are running far from 
their peak vector performance; each is about a factor of eight below peak.  
This implementation involves 4,000 lines of FORTRAN.

\begin{table}
\caption{Performance of One-Dimensional BEPS1 Plasma Physics Code
$\lbrack$Liewer:89a$\rbrack$}
\vspace{7pt}
\begin{center}
(a)\ \ {\bf Comparison of Push Times per Particle on Various
Computers}\\
(Particle Update Only)
\vspace{0.25in}
\begin{tabular}{||l|c|l||}
\hline
\multicolumn{1}{||c|}{Computer}
& \multicolumn{1}{c|}{Push Time} & \multicolumn{1}{c||}{\ }\\
& \multicolumn{1}{c|}{$\mu\/$secs} &\\
\hline
\medStrut
Mark IIIfp (64 processor) & ""0.8 & ($\sim$50 megaflops)\\
CRAY XMP/48 (1 processor) &     & \\
\ \ Vectorized & ""1.5 & ($\sim$25 megaflops)\\
\ \ Scalar     & ""4.1 &\\
CRAY 2 (1 processor)   &  & \\
\ \ Vectorized & ""2.1 &\\
\ \ Scalar     & "10.1 &\\
IBM 3090 VF &  & \\
\ \ Vectorized & ""2.9 &\\
\ \ Scalar     & ""6.0 &\\
Mark III (64 processor) & ""3.9 &\\
Alliant FX/8 & "12.6 &\\
VAX 11/750, F.P.A. & 200.9 &\\
Convex C-1 (vector) & "19.5 &\\
\hline
\end{tabular}
\end{center}

\bigskip
\begin{center}
(b)\ \ {\bf Total Run Time Comparison}\\
(Particle and Grid Updates)
\vspace{0.25in}
(720,896 particles, 1024 grid points, 1000 time steps)
\vspace{0.25in}
\begin{tabular}{lc}
Mark IIIfp (64 nodes) & 1062 secs\\
CRAY 2 (1 processor)  & 1714 secs\\
\end{tabular}
\end{center}
\end{table}

\section{Grain Dynamics by Lattice Gas (Gutt)}

A very interesting use of the NCUBE was recently reported in Gary Gutt's 
Ph.D. thesis at Caltech \cite{Gutt:89a}.  Earlier, Werner's Ph.D. research 
had studied grain dynamics using the first hypercubes \cite{Werner:87a}, 
\cite{Werner:88a}.  This work put in detailed Newtonian dynamics to study the 
motion of sand and other granular material.  This is an interesting 
alternative to conventional continuum approximations to material dynamics.  
Gutt proposes an intermediate model for such systems using cellular automata 
or lattice gas techniques that have already been applied to fluids 
\cite{Frisch:86a}.  Gutt's automata are quite dense (of order one for every 
two lattice sites) and one must store the relative displacements of each 
automata from the lattice site positions.  Thus, this automata method does 
not use binary arithmetic, but rather 32 bit arithmetic.  To improve 
performance, integer and not floating point arithmetic is used.  Gutt's 
thesis used about 200 hours of NCUBE 512 node time with the largest 
simulation involving $0.5\times 10^{6}$ grains on a $8064\times 128$ 
lattice.  This involves Poiseuille flow down a pipe driven by gravity.  We do 
not have a CRAY implementation of this code, but it is possible that 
irregularities in lattice site occupancy would make vectorization difficult.  
The parallelization on the NCUBE is straight forward and efficient.

\section{Computer Chess (Felten, Otto, Morison)}

Computer chess involves constructing a tree of possible moves and
dynamically pruning it with the $\alpha-\beta$ technique as illustrated
in Figure~13 \cite{Felten:88g}, \cite{Felten:88h}, \cite{Felten:88i}.
On a parallel machine, the tree (``data domain'') is decomposed over the
nodes.  We found that a real time graphics display (using the NCUBE parallel
graphics subsystem) was critical in achieving a factor of five better
performance.  This allowed us to change the algorithm for processor
assignment and improve the load balance.  The measured speed up is shown
in Figure~14 with a speed up of 101(170) seen on 256(512) nodes for trees
of depth appropriate for a middle game.  Note that the speed up
increases as the problem gets bigger, {\em i.e.,\/} as one spends a longer
time on each move.  This illustrates that for problems with real time
constraints, increasing processor performances increases parallelization
efficiency by increasing the size of the problem that can be solved in a given
fixed time.

This was a very difficult code to design and develop as the algorithm is
{\em asynchronous\/} in the classification of Sec.~2.  Probably, the 
technical implementation and in particular the debugging were harder than the 
algorithmic issues.  The result is 8,000 lines of C code using the commercial 
NCUBE operating system VERTEX with a ``special shared memory'' enhancement to 
allow concurrent access and update of a distributed database --- the so-called
transposition table of currently evaluated positions.  In contrast, the
earlier Caltech hypercube scientific calculations in this review, use our 
internal {\em loosely synchronous\/} communications system CrOS which is 
faster than VERTEX \cite{Fox:88a}.

Let us consider the future of computer chess.  We estimate that the 256
node NCUBE has a U. S. chess rating of 2,100 at present, or in a more
familiar unit of megamoves searched per second, we have the results in 
Table~8.

\begin{table}
\caption{(a) Chess Positions Searched per Second}
\vspace{7pt}
\begin{center}
\begin{tabular}{||p{2.25in}|c||}
\hline
\multicolumn{1}{||c|}{Machine} & \multicolumn{1}{c||}{Performance}\\
\multicolumn{1}{||c|}{\ } & \multicolumn{1}{c||}{($10^{6}$
moves/sec)}\\
\hline
\medStrut
{\bf Special Purpose} & \\
\ \ BELLE (1980) & 0.075\\
\ \ HITECH (1985) & 0.25"\\
\ \ Deep Thought (1988) & 2.5""\\
\ \ New AT and T Machine (1989?) & 5@"""\\
{\bf General Purpose} & \\
\ \ 256 node NCUBE & 0.025\\
\ \ 4 Processor CRAY XMP & 0.15"\\
\hline
\end{tabular}
\end{center}

\begin{center}
(b) Ratings of the Best Chess Players\\
\vspace{0.25in}
\begin{tabular}{||l|c||}
\hline
\multicolumn{1}{||c|}{Player} & \multicolumn{1}{c||}{Rating}\medStrut
\\
\hline
\medStrut
HITECH & 2450\\
Deep Thought & 2600\\
Kasparov & 2850\\
1000$\times\/$Deep Thought & 2850\rlap{?}\\
\hline
\end{tabular}
\end{center}
\end{table}

The NCUBE does quite well compared to a CRAY (a 512 node NCUBE is
approximately one head of a CRAY XMP), but neither is competitive with 
special purpose machines.  However, we can consider using the same 
parallelization technique developed for the NCUBE, but for an array of 
special purpose chess chips --- not the general purpose microprocessor used 
on the NCUBE.  We conjecture that a system of some 8,000 special purpose 
chips like those in Deep Thought, would achieve a speed up of 1,000 and
be very competitive with Kasparov as shown in Table~8.  Each node processor 
of our world chess champion would have for chess about 100 times the power of 
a NCUBE node.  The communication overhead for the NCUBE case is about 10\%; 
the ``Deep Thought'' chips would need an internode bandwidth about ten times 
that of NCUBE to keep a manageable (50\%) communication overhead.  Such a 
system seems quite practical, but outside my group's resources.  Rather, we 
are concentrating on a different approach; can one improve chess programs by 
the use of neural network based position evaluators?

\section{Ray-tracing in Computer Graphics (J. Salmon, J. Goldsmith
(JPL), and S. Huang)}

The 512 processor NCUBE was used to render a one minute long animation 
entitled ``Self-Portrait'', which shows a hypercube grow from zero through 
four dimensions and then rotate in four dimensions. \cite{Huang:89a},
\cite{Goldsmith:88a}.  The sequence consists of approximately 600 frames, 
each of which was rendered at $2048\times1944$ pixel resolution and filtered 
to $512\times486$ for transfer to videotape.  The most complex scenes contain
169 polygons, 16 spheres, 32 cylinders and three light sources.  Over the 
course of the sequence, approximately 8~billion rays were shaded and 
approximately 200~billion ray intersection tests were performed.  The
entire calculation required about 300 hours on the 512 processor NCUBE
system, and a typical frame is shown in Figure~15.

Ray-tracing is generally implemented as an {\em embarrassingly parallel\/}
algorithm.  In order to render an image which consists of a large number of
pixels, each pixel is considered separately and independently.  The 
calculation for each pixel is, itself, quite complex, involving anywhere from 
tens to thousands of floating point operations, depending on the complexity 
of the scene.  Since computer graphics images generally have $512\times512$ 
or more pixels, a  decomposition in which individual pixels are assigned to 
processors can make efficient use of hardware with up to $O(10^{5})$ 
processors.  A minor subtlety arises because even within a single image, the
complexity of individual pixels can vary by several orders of magnitude.  
Pixels at the edge of the scene, whose rays may miss the model entirely, are 
very easy to compute, while those near the ``interesting'' part of the image 
may involve hundreds of intersection and shading calculations.  If pixels 
were allocated to processors according to the ``tiled decomposition'' of 
Figure~16(a), the large fluctuations in the time to compute different pixels 
would translate into severe load imbalance.  The solution to this dilemma is 
to use the ``scattered decomposition'' shown schematically in Figure~16(b).  
Now each processor computes representative pixels from all over the image, 
and severe load imbalance is extremely unlikely.  In fact, we can estimate 
that the load imbalance is proportional to $n^{-1/2} \left( 
\frac{\sigma_{pix}}{\mu_{pix}} \right)$, where $\mu_{pix}$ and $\sigma_{pix}$ 
are the mean and standard deviations of the distribution of pixel computation 
times, and $n$ is the number of pixels per processor \cite{Salmon:88b}.

A much more significant difficulty with the {\em embarrassingly parallel\/}
approach to ray-tracing is that it requires that the entire database be
available to each and every processor.  This is no problem, in principle, for 
shared memory machines in which the database is stored once in shared memory, 
but it is a major restriction for distributed memory machines, as it requires 
that a separate copy of the database reside in each processor's local 
memory.  Even for shared memory architectures, the multiple access of nodes 
to a single database may introduce contention and performance degradation.  
Thus, embarrassingly parallel distributed memory ray-tracers are only capable 
of rendering models which fit, in their entirety, within the memory of a 
single processor.

To address this problem, we are forced to distribute the database as well as 
the computation.  It is necessary, at this point, to have some understanding 
of how the ray-tracing algorithm proceeds.  For each pixel, a ``ray'' is 
generated which passes from the ``view-point'' through the pixel on the 
screen, and into the model space.  Once this is known, a shading calculation 
is performed, which usually involves projecting a number of additional rays 
from the point of intersection toward the light-sources in the model and in 
the directions given by Snell's and Hero's Laws.

By far, the most time-consuming part of the computation is the determination 
of which primitive object is struck by a ray.  A complex scene can have up to 
a million primitive objects, so a naive linear search is clearly out of the 
question.  Several related techniques are common for reducing the complexity
of the intersection to $O(log N_{obj})$.  The Caltech group elected to 
construct a tree of ``bounding volumes'', in which the leaves are primitive 
objects, and each internal node is a rectangular box, large enough to contain 
all of its descendants.  The intersection calculation begins by testing the 
ray with the root bounding box.  If the root bounding box intersects the ray, 
then each of its children is tested recursively.  Whenever a bounding box 
fails to intersect the tree, all its children may immediately be pruned as 
candidates for intersection as well.

The tree of bounding volumes is a natural data structure to distribute to the 
processors.  A copy of nodes near the root of the tree is duplicated in each 
processor.  Thus, each processor can begin an intersection calculation with 
data available locally.  We call this common data structure the ``forest'' 
because its terminal elements are pointers to sub-trees which are stored only 
once, in one of the processors of the parallel ensemble.  When a processor 
needs to know the results of an intersection calculation between a ray and a 
sub-tree, a remote procedure call is made, the remote processor which 
maintains the appropriate part of the database performs the calculation and 
sends back the result in a message.  

The communication traffic in this program is highly asynchronous and 
unpredictable.  This is one of the few algorithms we have studied which would 
not perform well under a loosely synchronous message passing operating
system.  In fact, it is not just the communication, but also in the
language of Sec.~2, the computation as a whole that is {\em asynchronous\/}.  
A processor must divide its time between its own work (traversing the forest 
for pixels which it is responsible for computing) and that requested by 
others.  In order to manage the highly asynchronous computation, the MOOSE 
operating system \cite{Salmon:88b}, \cite{Koller:88a}, was developed to
allow numerous lightweight processes to run on each processor.  Each sub-tree
that is assigned to a processor is represented by a separate lightweight
process, with its own queue of requests.  A process is active when its input 
queue is non-empty.  Otherwise, the process is idle, and does not consume CPU 
resources.  Note that computer chess described in the previous section
was also {\em asynchronous\/}, but involved larger grain chunks which
were irregular but static.  This application did not need MOOSE, but
rather the conventional asynchronous message passing system VERTEX
supplied with the NCUBE.

\section{Kalman Filters (T. Gottschalk)}

A JPL team headed by D. Curkendall \cite{Meier:89a} has developed a 
sophisticated battle management simulation.  This includes threat generation 
(launch missiles), tracking, engagement planning (launch anti-missiles), and
graphics.  The total of 200,000 lines of code is one of the largest
single parallel computer projects --- much larger than the sum of codes
described so far in the previous eleven subsections!  The current, so called, 
SIM88 project completely simulates up to 250 objects launched from six 
sites.  There is an interesting hybrid approach to the simulation with each 
component (tracking satellite, planning platform) functionally decomposed 
with a very coarse grained object oriented model.  Traditional data 
parallelism is used with functions (objects) assigned to a subcube of the 
hypercube.  Correspondingly, a hybrid software model CENTAUR for the 
Mark~IIIfp hypercube supported general but slow inter object communication 
and fast loosely synchronous CrOS communication within objects 
\cite{Lee:89a}.  Formally, the simulation is an {\em asynchronous\/} event 
driven simulation at the object level, but the predictable and coarse grain 
nature of the object to object communication allowed efficient {\em loosely 
synchronous\/} implementation within a simple conservative framework 
\cite{Chandy:81a}.

Let us focus on the tracking component which used a parallel implementation 
of a traditional multi-object Kalman filter tracker \cite{Gottschalk:87f}, 
\cite{Gottschalk:88a}.  Much of the calculation can be done independently in 
each node when one distributes the tracks.  There are some significant 
overheads when tracks share measurements and load balance is an issue 
addressed by dynamically redistributing the data at each measurement
cycle.  This is a loosely synchronous algorithm implemented as 4,000
(20,000) lines of C in the SIM~87(88) versions.

\cite{Messina:89a} has a complete discussion of the performance of the
SIM~87 code on a variety of machines showing good performance for shared
memory machines.  These avoid the data shuffling overhead for
overlapping tracks.  This algorithm was designed for the so called boost
phase when there is modest parallelism coming from the decomposition of
a total of a few hundred objects.  The 480 target problem shown in
Figure~17 \cite{Baillie:88f}, does not perform well on the larger NCUBE
systems due to the limited parallelism.

We note that the communication overhead is perhaps overestimated in Figure~17 
as it depends critically on the sophistication of the track model.  The newer 
SIM~88 implementation has a fully parallel tracking component that takes 
three times longer than SIM~87.  The communication overhead is reduced by 
this factor for SIM~88.

We expect the situation to change in the harder midcourse phase where up
to $10^{5}$ real or decoy objects can be anticipated.  This will need
new algorithms --- perhaps neural networks \cite{Fox:89h} --- and the
parallelism issues will be different, but we expect it to be easier and
not harder to exploit massively parallel machines.  In the boost phase, it
appears that existing parallel machines with a few megaflops of
performance will allow real time tracking; mid course will require true
supercomputers.
 
\section{Lessons}

Here we collect together some lessons we have drawn from ``using real
parallel computers to solve real problems with real software''.

This information is obtained from both the supercomputer applications of
Sec.~3 to 14 and for instance the earlier surveys of \cite{Fox:85c}, 
\cite{Fox:88b}, \cite{Fox:88ll}, and \cite{Fox:89i}.

\subsection{Application Lessons}

\noindent{\bf 15.1.1}\ \ One can achieve high performance on essentially all 
scientific computations which are 
\begin{itemize}
\item Large (necessary condition)
\item {\em Loosely synchronous\/} --- MIMD (sufficient condition)
\item or \hspace{1.0cm}{\em synchronous\/} --- SIMD
\end{itemize}
More research is needed to clarify Figure~3 and Table~2 and see how far one 
can violate synchronization and still get good results for irregular problems 
on SIMD machines.  Note that even for loosely synchronous problems, 
synchronous communication is usually sufficient; communication is typically 
only necessary at the macroscopic synchronization points of the algorithm
\cite{Fox:88a}, \cite{Fox:88b}, \cite{Fox:88ll}.

\medskip

\noindent{\bf 15.1.2}\ \ Domain Decomposition or data parallelism is a 
universal source of parallelism that scales to large numbers of
nodes \cite{Fox:88a}, \cite{Hillis:86a}.

\medskip

\noindent{\bf 15.1.3}\ \ These results are true on a broad range of computer 
architectures (SIMD, MIMD, shared, distributed memory, hypercubes, 
(transputer) meshes $\ldots$ ).

\medskip

\noindent{\bf 15.1.4}\ \ University successes on parallel computers have come 
with 1,000--10,000 line codes written from scratch for a particular machine.

\medskip

\noindent{\bf 15.1.5}\ \ It is not clear how to extrapolate these successes 
to up to $10^{6}$ line commercial codes where you have perhaps less knowledge 
as to the inner workings of program.  Certainly, one must establish and use 
standard methods in parallel software to justify expense of a new parallel 
implementation.  These standards must apply across a range of architectures.

\medskip

\noindent{\bf 15.1.6}\ \ In many cases ({\em cf.\/} Sec.~5), it is easier to 
decompose for a parallel machine than to vectorize for a conventional 
supercomputer.  This is especially true for small (university) codes.

We can superficially abstract from this that universities should
purchase parallel computers and industry vector supercomputers!

\medskip

\noindent{\bf 15.1.7}\ \ The importance of parallel machines for artificial 
intelligence (AI) is unclear to me.  If the AI is implemented with neural 
networks, then the relevance and use of parallel computers is clear.  More 
traditional AI systems parallelize less easily although chess (Sec.~2.11) is 
a good example which both parallelizes and needs high performance.  How many
other such AI applications are there?

\subsection{Performance Lessons}

\noindent{\bf 15.2.1}\ \ One can get high performance on essentially all 
scientific computations.  As shown in \cite{Fox:84c}, \cite{Fox:85c},
\cite{Fox:88a}, \cite{Gustafson:88a}, and \cite{Fox:89b}
\begin{itemize}
\item Performance scales linearly in number of nodes at constant grain
size (problem size proportional to machine size)
\item Fixed problem size does not scale; this can be viewed as Amdahl's
Law
\end{itemize}

\noindent{\bf 15.2.2}\ \ Some initial disappointments can be traced to 
imbalance in early commercial machines, such as the iPSC/1 and FPS T~Series 
hypercubes.

\medskip

\noindent{\bf 15.2.3}\ \ We saw in Secs.~5--7, and Sec.~13 that machines 
like an NCUBE or transputer arrays, with many simple nodes, look particularly 
attractive compared to the CRAY~XMP class machines on irregular problems 
where one finds 
\begin{itemize}
\item it is more natural and easier to decompose than to vectorize
\item the NCUBE efficiency is ``low''; maybe 50\%
\item but the CRAY efficiency even lower; maybe 5\%.
\end{itemize}

Note that the average CRAY-XMP performance in computer center operation is 
about 25 megaflops with a 12\% efficiency.  The NCAR CRAY realizes a
sustained 50 megaflops, which is perhaps the peak average performance.  
We have too high a standard for the efficiency of parallel machines!

\medskip

\noindent{\bf 15.2.4}\ \ Scalar node MIMD machines are natural general 
purpose machines with reasonable performance over a range of problems.  As 
seen in Table~6, hypercubes have high efficiency on regular (e.g., full 
matrix and QCD) problems, but so does the CRAY XMP class machine.  Hypercubes 
with vector nodes or SIMD machines are attractive for regular problems.

\medskip

\noindent{\bf 15.2.5}\ \ Different programming methodologies and lack of 
standards handicap performance studies.  As FTN-8X is yet to be implemented 
uniformly, it is hard even to port between CRAY XMP and the ETA-10.  This, for
instance, is handicapping the ``grand challenge'' mentioned in Sec.~3.5.  Of 
course, porting between parallel and sequential machines is hard and 
\cite{Messina:89a} essentially re-implemented from scratch many algorithms 
for their performance evaluation.  The different software methodologies for 
shared and distributed memory machines cannot be avoided as some sense it is 
the more convenient environment that motivates shared memory machines.  It 
would be unreasonable to require shared memory machines to always use message
passing!

\subsection{Decomposition Lessons}

Here we refer to issues concerning the decomposition or dividing up of 
problems to minimize communication and equalize load on processors.

\medskip

\noindent{\bf 15.3.1}\ \ Three years ago, I thought decomposition or load 
balancing was a key problem, but as described in Sec~4 and
\cite{Williams:89b}, it is surprisingly easy!
\begin{itemize}
\item Usually, the application scientist can specify it from the natural
geometric structure of the problem
\item Several heuristic methods provide automatic decomposition.  These
include recursive bisection, simulated annealing and neural networks.
\end{itemize}

\noindent{\bf 15.3.2}\ \ Current hardware trends have emphasized transparent 
message routing where the user need not be aware of machine topology.  This 
is clearly convenient, but Sec~15.3.1 indicates it is not strictly necessary 
for a broad class of problems.  We can note that:
\begin{itemize}
\item Most problems can be mapped with software to any reasonable bandwidth
topology with modest routing overheads.
\item In particular, for current Caltech codes (say those described in
Sec.~3--14) there is an average overhead of less than $5\%$ due to routing.  
This should be compared to overheads of perhaps 50\% due to poor compilers 
(e.g., for NCUBE and WEITEK), 50\% as system is overall rather flaky and 
25\% due to node to neighboring hypercube node communication.
\item We must emphasize that these software solutions have yet to be
``packaged'' nicely for general use.  Not every programmer is
comfortable with simulated annealing.  Thus, automatic routing hardware
is certainly convenient in the ``real world''.
\end{itemize}

\noindent{\bf 15.3.3}\ \ In current systems, message start up time, which
includes hardware and software effects is a much 
more serious overhead than either node to node through routing, or channel 
transmission between neighboring nodes.

\subsection{Hardware Lessons}

\noindent{\bf 15.4.1}\ \ So far, high performance computations on moderately 
or massively parallel machines ($\gea 8$ nodes) has been confined to 
distributed memory machines.  The comparison between distributed and shared 
memory architecture is hard because of the lack of comparable machines and
experience.

\medskip

\noindent{\bf 15.4.2}\ \ Five years ago, there were many university projects 
building novel machines, but in the future commercial systems will
dominate the parallel field as they now do with conventional supercomputers.

\medskip

\noindent{\bf 15.4.3}\ \ The U. S. entrepreneurial environment will guarantee 
a wide range of architectures even without new university projects.  Portable 
software will be at a premium to exploit these machines.

\medskip

\noindent{\bf 15.4.4}\ \ Many or perhaps all the current commercial parallel 
systems are disappointing in some ways.  For instance, our NCUBE is now in 
full use as a ``production supercomputer'' but this took two years and at 
least \$150K costs in software development at Caltech.  If we remember Bill
Joy's law that sequential computers improve a factor of two per year in
cost-performance, we see that a two-year delay translates into a factor
of four loss in cost effectiveness compared to the conventional
competition.  Novel computers are bound to need extra time to develop
viable software than their sequential competition; this certainly
handicaps their ability to compete.

\medskip

\noindent{\bf 15.4.5}\ \ As implied by Sec.~15.4.4, systems integration is 
not yet well addressed in the parallel machines.  This includes issues such 
as,
\begin{itemize}
\item General multiuser operating systems especially for distributed
memory machines; debugging.
\item Adherence to standards.
\item Input/Output for disks and graphics.  The architecture of the I/O
system gets surprisingly little attention in the literature.
\item High performance appropriate hosts (not PC's or workstations) for
a parallel supercomputer.
\end{itemize}

\medskip

\noindent{\bf 15.4.6}\ \ For MIMD machines, we have already discussed some 
issues in Secs.~15.1.1, 15.2.3, and 15.2.4.
\begin{itemize}
\item Currently, machines like the NCUBE or transputer with scalar (floating 
point) nodes seem more successful than machines with vector nodes.
The poor old user finds it hard to vectorize and decompose problems!
One such optimization is enough.
\item Given that such scalar node systems are particularly attractive
for irregular problems, maybe one should consider adding specialized
support for the data structures like linked lists needed for irregular
problems.
\item Several megabytes (but not arbitrarily large) memory per node is
needed to hold program, decomposed data, databases and reused communication 
data (see Sec.~7).
\item One can expect the differences between shared and distributed
memory architectures to lessen as both are based on low latency
networks.  Either local memory (for machines like the hypercube) or
caches on shared memory machines will require data locality for good
performance.
\end{itemize}

\noindent{\bf 15.4.7}\ \ SIMD machines can support at least 50\% of 
university scientific applications as shown in Table~2.  They currently give 
the peak performance for regular problems (see Tables~4, 6).  Perhaps the
commercial applications are more irregular and will show a lower
fraction appropriate for the SIMD architecture.

\noindent{\bf 15.4.8}\ \ In Figure~18, we show a possible structure of an 
integrated high performance novel architecture computer environment.  Simple 
scalar node MIMD machines support general problems with either vector or SIMD
architectures as an ``accelerator'' for regular problems.

\subsection{Software Lessons}

\noindent{\bf 15.5.1}\ \ Whereas the role of universities in developing 
hardware systems may be limited in the future, we expect universities to have 
a critical role in the software for parallel machines where we cannot hope 
for the commercial systems to be adequate.

\medskip

\noindent{\bf 15.5.2}\ \ A key question is:  ``What is the appropriate 
productive standard programming environment for parallel machines?''  This 
could be based on \cite{Fox:88u}, \cite{Fox:88w}:
\begin{itemize}
\item New languages and approaches (e.g., graphical techniques)
\item Compiler generated parallelism
\item Application specific high level environments
\item Explicit user decomposition
\end{itemize}

\noindent{\bf 15.5.3}\ \ Note that essentially all successful reasonable 
performance use of parallel machines have used explicit user decomposition 
which is low level and machine dependent.  We expect we must find more 
portable attractive methods if parallel computers are to take over from the
conventional architectures.

\medskip

\noindent{\bf 15.5.4}\ \ The asynchronous applications of Secs.~12--14,
need more sophisticated operating environments than the loosely
synchronous or embarrassingly parallel applications of the earlier
sections.  Each of the last three applications required different
software systems (shared memory VERTEX, MOOSE, and CENTAUR respectively)
for good performance.  The simple, but high performance, synchronous
CrOS communication system \cite{Fox:88a} is adequate for Secs.~3 to 11.

\medskip

\noindent{\bf 15.5.5}\ \ Approaches like LINDA (``shared message space'') or 
the new language, OCCAM, appear not to address enough of the issues to
be the solution of the question posed in Sec.~15.5.2.

\medskip

\noindent{\bf 15.5.6}\ \ When we started work on the hypercube in 1981, I had 
great confidence that they would be successful as I could make a simple 
performance model and prove ``mathematically'' that for dedicated users, 
hypercubes would work \cite{Fox:83a}.  I don't have any way of making a 
similar prediction for the software environment!

\begin{center}
{\bf Acknowledgements}\\
\end{center}

I would like to thank Paul Messina and the many others whose research is
reviewed in this paper.

\newpage
\bibliography{/home/sirius/kea/bib/doc.bib,/home/sirius/kea/bib/outside.bib,/home/sirius/kea/bib/spcp.bib}
\bibliographystyle{c3p}
\end{document}
