
\documentstyle[12pt]{article}
\def\registeredsymbol{{\ooalign{\hfil\raise.05ex\hbox{\hskip.02ex$\rm
                                                      \scriptscriptstyle R$}\hfil\crcr
                                \hbox{$\scriptstyle\mathchar"20D$}}}}
\def\registered{\ifmmode^\registeredsymbol\else$^\registeredsymbol$\fi}
\setlength{\parskip}{7.5pt}
\newlength{\boxwidth}
\setlength{\boxwidth}{\textwidth}
\addtolength{\boxwidth}{-20 pt}
\newcommand{\boxer}[1]{\begin{center}
                        \fbox{\fbox{
                                \begin{minipage}{\boxwidth}  #1
                                \end{minipage}}}
                        \end{center}}

\newcommand{\RI}{{\uppercase\expandafter{\romannumeral1}} }
\newcommand{\RII}{{\uppercase\expandafter{\romannumeral2}} }
\newcommand{\RIII}{{\uppercase\expandafter{\romannumeral3}} }
\newcommand{\RIV}{{\uppercase\expandafter{\romannumeral4}} }
\newcommand{\RV}{{\uppercase\expandafter{\romannumeral5}} }
\begin{document}

\title{Hypercube Algorithms for Direct N-body Solvers
for Different Granularities}

\author{
Jean-Philippe Brunet\thanks{
Thinking Machines Corporation,
245 First Street, Cambridge, MA 02142, brunet@think.com, jill@think.com}
\\
Alan Edelman\thanks
{Department of Mathematics, University
of California, Berkeley, CA 94720, edelman@math.berkeley.edu}
\\
Jill P. Mesirov$^*$
}

\date {(New improved version -- Friday September 27)}

\maketitle

\begin{abstract}

We compare and contrast two algorithms for the N-body
problem, particularly where $N$ is in the range  
for which direct methods outperform indirect methods.
With fewer bodies than processors, the so-called ``replicated
orrery" on a three dimensional grid
has been used successfully on the Connection Machine 
CM2 architecture.  With more bodies, the ``rotated 
and translated Gray codes" is an ideal direct algorithm
for machines such as the CM2 in that it takes optimal advantage
of the communications bandwidth of the machine.

A classical Latin square can be used to abstractly denote
any direct N-body calculation.  Computational windows
superimposed on this square illustrate the granularity of
the computation.
This point of view, naturally illustrates 
the continuity of a range of algorithms in the following
order: ``massively parallel," ``replicated orrery",
``orrery", ``rotated/translated Gray codes", and ``serial".
 

\end{abstract}

\section{{\protect\large {\protect\bf Introduction}}}


In a wide variety of applications areas, including astronomy,
molecular 
dynamics
\cite{heller,raine,smith,windemuth}, 
and fluid dynamics \cite{sethian}, a serious computational
bottleneck occurs from the need to determine the mutual interaction of
$N$ bodies.  Directly computing the $N^2$ pairwise interactions
given by
$$
\mbox{Force(body $i$)} = \sum_j \mbox{Force(body $i$ due to body $j$)}
$$
on a
conventional computer can take an unacceptably long amount of time for
values of $N$ which are exactly in the range of interest for the
application.

Indirect methods are one approach towards  
reducing  computational complexity.
Examples of this are the particle/mesh \cite{hockney}, local correction
\cite{anderson},  hierarchical \cite{apel,barnes}, and  multipole
\cite{greengard} methods.
The structure of such algorithms has been investigated by 
\cite{Katzenelson}.

Another
approach is to keep to the naive direct method, but to design algorithms
that take advantage of  advanced computer architectures to
obtain results  in a reasonable amount of time
through the use of parallel computation and the minimization
of communication.
We have chosen to explore direct methods in this work.  Although this
clearly has asymptotically higher complexity, these algorithms remain
usueful for two important reasons 1) they are faster for a large
range of $N$ on certain architectures and
2) they are simpler to code and maintain
thus miniminizing the important commodity of programming effort.
Though the indirect algorithms must ultimately be advantagous
for large values of $N$, for many problems $N$ is in the range where
direct body algorithms are superior.
See \cite{sethian} for further discussion of the advantages of a
direct algorithm over an indirect algorithm.



The algorithms for direct $N$-body solvers that  we will discuss in this
paper are designed for uses on a SIMD  parallel architecture with a hypercube
network topology. In fact, we feel that hypercube architectures are
ideally suited to the direct method.  For non-direct methods on the
hypercube, see \cite{Johnsson,Fox,B-H}.

We focus on two different algorithms, but we take a point of view
that we have not yet seen in the literature:
we show that the two algorithms fit along a continuous spectrum.
Perhaps we can make an analogy with the notion that X-rays and
radio waves properly fit along one axis.  However, instead of
using frequency, the important parameter is {\em  granularity} measured
here by $P/N$.


For the case where the number
of processors, $P$, satisfies $N \leq P \leq N^2$ we describe an
optimal mapping of a family of {\em replicated orreries} onto a
hypercube.  For the case where $P < N$, we present a multiwire
all-to-all broadcast algorithm which utilizes ``rotated and translated
Gray codes'' to make optimal use of the communication bandwidth of the
hypercube.  Thus, in this approach, the communication time of the
algorithm is negligible for large values of $N$.

Both of these algorithms have been implemented on the Connection
Machine CM-2.  There are two programming models for the CM-2 and each
algorithm lends itself more naturally to one of them.
The older model (which is becoming more and more obsolete)
is known as the ``fieldwise model", while the more recent model
is the ``slicewise model." 
From a conceptual point of view, the value of $P$ is 32 times
greater for the fieldwise model as it is for the slicewise model.
Both
approaches, and thus both programming models, were used for the
$N$-body solver in the implementation of a random vortex method for
numerical simulation of fluid flow \cite{sethian}; we will analyze the
computational complexity of each approach. 

Section 2 discusses N-body problems in an abstract mathematical context. 
In Sections 3 and 4 we present two concrete algorithms with complexity
analyses:  the replicated orrery
and a multiwire
algorithms.  The CM-2
implementation of both algorithms is briefly  described in Section 5.  A sample
application is described in Section 6 and the performance of the CM-2
implementation of both approaches is discussed.


\section{{\protect\large {\protect\bf Latin squares and N-body problems }}}

Before we analyze in detail two particular approaches to solving
direct N-body problems that have been successful, we would like
to abstract out the important mathematical features of an N-body
computation on a parallel machine.  In particular, the two
important questions are {\bf where} (which processor?) and {\bf when}
(which timestep?) does each interaction occur. 

Implicit in the knowledge of the ``where and when'' of the computation is
the communication that occurs.  Since we know which data needs to be
where at a particular moment, the only degree of freedom might be
how soon the data arrives at a particular processor before
a computation starts.

As a mathematical abstraction, we consider it enlightening
to display the ``where and when" of an N-body computation
on a {\bf Latin square}.   The basic idea is quite simple, 
however almost any N-body problem on any parallel machine can
be defined as a Latin square. 

A {\bf Latin square} is an $N$ by $N$ table of entries
$1,2,\ldots,N$ with the property that each number $i$
appears exactly once in each column and once in each row. 
Perhaps the earliest know reference is a puzzle
requiring that the 16 lettered cards (A,K,Q,J) of an ordinary 52 card deck be
arranged so that no letter appears twice in the same column or row
nor along a diagonal \cite{LATIN}.
Classically Latin squares have had many uses including the design
of experiments,  displays of multiplication tables of finite groups,
and other recreational puzzles (See~\cite{BBC,LATIN}).
Their first systematic study was performed by Euler.\cite{LATIN}

\begin{figure}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\multicolumn{1}{c}{1\hspace{-.7in} Body: \hspace{.15in} } &
\multicolumn{1}{c}{2} &
\multicolumn{1}{c}{3} &
\multicolumn{1}{c}{4} \\
\hline 

2 & 3 & 4 & 1 \\ \hline
3 & 4 & 1 & 2 \\ \hline
4 & 1 & 2 & 3 \\ \hline
1 & 2 & 3 & 4 \\ \hline
\end{tabular}
\end{center}
\setlength{\baselineskip 11pt}
\caption{
\footnotesize
A Latin square of Order 4. To indicate an N-body calculation,
we let  the entry in row $i$ and column $j$
denote an interaction between that body and body $j$.
A serial algorithm could proceed from left to right and from top to
bottom along this square, one box at a time.
}
\label{AA}
\end{figure}


As an example, Figure \ref{AA} displays a Latin square of order 4.
In the context of an N-body computation, the entry $a$ in row
$i$ and column $j$ can symbolize a force law interaction of body
$a$ with body $j$.  Thus in Figure \ref{AA} above, 
the first column depicts the interactions of body 1 with
bodies 2,3,4, and 1 respectively.  We remark that in practice,
self-interactions are usually omitted.  Also in practice,
symmetric interactions may be avoided, giving a speedup of
a factor at most equal to 2.


\begin{figure}
\begin{center}
\begin{tabular}{|l|c|l|r|r|}
\multicolumn{3}{c}{ }&\multicolumn{2}{c}{Example}  \\ 
\multicolumn{3}{c}{ }&\multicolumn{1}{r}{P}& \multicolumn{1}{r}{N}   \\ \hline 
Case \RI & $P=N^2$ & Massively Parallel & 16 & 4 \\
Case \RII& $N<P<N^2$ & Replicated Orrery  & 8 & 4 \\
Case \RIII & $P=N$ & Orrery             & 4 & 4 \\
Case \RIV & $1<P<N$ & Rot/Trans Gray Codes  & 2  & 4 \\
Case \RV & $P=1$ & Serial & 1 & 4 \\ \hline 
\end{tabular}
\end{center}
\label{BB}
\caption{Possible granularities (Example for $N=4$)}
\end{figure}



\newcommand{\text}[2]{\put (#1){\makebox(0,0){#2}}}
\begin{figure}

\setlength{\unitlength}{0.7cm}
\begin{picture}(0,5)(-5,-5)
\multiput(-8,0)(5,0){5}{
\multiput(0,0)(0,-1){5}{ \multiput(0,0)(.2,0){20}{\line(1,0){.1}}}
\multiput(0,0)(1,0){5}{ \multiput(0,0)(0,-.2){20}{\line(0,-1){.1}}}
\text{0.5,-0.5}{$p_1$}
}

\multiput(-8,0)(5,0){4}{
\text{1.5,-0.5}{$p_2$}}

\multiput(-8,0)(5,0){3}{
\text{2.5,-0.5}{$p_3$}
\text{3.5,-0.5}{$p_4$}}

\multiput(-8,0)(5,0){2}{
\text{0.5,-1.5}{$p_5$}
\text{1.5,-1.5}{$p_6$}
\text{2.5,-1.5}{$p_7$}
\text{3.5,-1.5}{$p_8$}}

\multiput(-8,0)(5,0){1}{
\text{0.5,-2.5}{$p_9$}
\text{1.5,-2.5}{$p_{10}$}
\text{2.5,-2.5}{$p_{11}$}
\text{3.5,-2.5}{$p_{12}$}
\text{0.5,-3.5}{$p_{13}$}
\text{1.5,-3.5}{$p_{14}$}
\text{2.5,-3.5}{$p_{15}$}
\text{3.5,-3.5}{$p_{16}$}}

\put(-8,-4){\framebox(4,4){}}
\put(-3,-2){\framebox(4,2){}}
\put( 2,-1){\framebox(4,1){}}
\put( 7,-1){\framebox(2,1){}}
\put(12,-1){\framebox(1,1){}}

\text{-6,-4.5}{\RI:$P=16$}
\text{-1,-4.5}{\RII:$P=8$}
\text{ 4,-4.5}{\RIII:$P=4$}
\text{ 9,-4.5}{\RIV:$P=2$}
\text{14,-4.5}{\RV:$P=1$}
\end{picture}
\setlength{\baselineskip 11pt}
\caption{
\footnotesize
Computational windows superimposed on the  Latin square ($N=4$).
The entries of the square are omitted, while
the processor numbers for the first timestep of the computation
are included.
Later timesteps are obtained by sweeping the window accross the
square, so as to completely cover the square.
}
\label{CC}
\end{figure}

In the sections to follow, we study the five possibilities
given in Figure \ref{BB}.  Spefically the replicated orrery is described
in detail in Section 3, while the rotated and translated Gray codes is
described in Section 4. 
Before we give the details of the algorithms, however,
we can already describe some necessary features based on the granularity
alone.

For each of the five cases in Figure \ref{BB},
Figure \ref{CC} exhibits the simultaneity of one parallel
computation
when $N=4$.
For clarity of exposition, we omit the entries of the square, but
instead superimose 
a computational ``window" upon the square.  The size of this
window is exactly the number of processors.  
In each window, there is one intereaction for each processor.
The computation proceeds by sweeping the window accross the square.
Of course, if  the number of bodies is held constant, then the
size of the window grows with the number of processors, and the
number of timesteps (equals the number of windows needed to cover
the square) decreased.

Notice in Figure \ref{CC} that in Case \RI there is only one window --
all N-body computations occur in parallel.  On the other hand in
Case \RV, the algorithm is serial, one interaction occurs
in any given time-step.
The processors are labeled $p_i$ in the windows, and the computation
proceeds by shifting the windows and the corresponding processors
over the Latin square.


\begin{figure}
\begin{tabular}{|cccc|cccc|}
\multicolumn{8}{c}{body} \\
\multicolumn{1}{c}{1} &
\multicolumn{1}{c}{2} &
\multicolumn{1}{c}{3} &
\multicolumn{1}{c}{4} &
\multicolumn{1}{c}{5} &
\multicolumn{1}{c}{6} &
\multicolumn{1}{c}{7} &
\multicolumn{1}{c}{8} \\ 
\multicolumn{8}{c}{interacts with} \\
\hline
2 & 3 & 4 & 1 & 8 & 5 & 6 & 7 \\ \hline
8 & 5 & 6 & 7 & 2 & 3 & 4 & 1 \\ \hline\hline
3 & 4 & 1 & 2 & 7 & 8 & 5 & 6 \\ \hline
7 & 8 & 5 & 6 & 3 & 4 & 1 & 2 \\ \hline\hline
4 & 1 & 2 & 3 & 6 & 7 & 8 & 5 \\ \hline
6 & 7 & 8 & 5 & 4 & 1 & 2 & 3 \\ \hline\hline
1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\ \hline
5 & 6 & 7 & 8 & 1 & 2 & 3 & 4 \\ \hline
\end{tabular}
\begin{picture}(10,10)(-10, 5)
\put(30,-30){\framebox(80,80){}}
\text{26,56}{1,5}
\text{36,46}{$p_1$}
\text{114,56}{2,6}
\text{104,46}{$p_2$}
\text{114,-36}{3,7}
\text{104,-23}{$p_3$}
\text{26,-36}{4,8}
\text{36,-23}{$p_4$}
\text{70,-55}{2 Dimensional Hypercube}
\put(130,0){\framebox(80,20){}}
\text{140,10}{$p_1$}
\put(150,0){\line(0,1){20}}
\text{160,10}{$p_2$}
\put(170,0){\line(0,1){20}}
\text{180,10}{$p_3$}
\put(190,0){\line(0,1){20}}
\text{200,10}{$p_4$}
\text{170,-15}{Comp Window}
\end{picture}
\setlength{\baselineskip 11pt}
\caption{
\footnotesize
Latin square and computational windows for an 8-body problem
with 4 processors.  Each 1 by 4 rectangle denotes the four interactions
of one time step, one in each processor.
The entries of the Latin square are obtained by circulating
a dynamic copy of the first body counterclockwise, while circulating
the second body clockwise on the square.  A static copy of each
body remains in the processor.
Communication is indicated by the vertical space between pairs of lines.
After the first circulation, body 1 can interact with bodies
2 and 8, while body 5 can also interact with bodies 2 and 8 without
further communication.  Thus a total of four windows can pass over the
square before another communication need occur.
}
\label{DD}
\end{figure}

We now see that once the relationship between $P$ and $N$ is given,
all that we need to specify a direct N-body computation
is 1) a choice of Latin square and 2) a mapping from the processors
to the window.  We explore this in further detail  
in later sections, but we can give an illustrative example,
as we do
in Figure \ref{DD}.
This figure
illustrates an example when $N=8$ and $P=4$. 
The square on the left is readily checked to be a Latin Square.
The construction of the square is quite simple.  We assume that
there are two bodies per node as illustrated on the 2 dimensional
hypercube in the figure.
One copy of the bodies remains static, while one circulates, the
first body circulates counterclockwise, while the second body circulates
clockwise.





\section{{\protect\large {\protect\bf $N \leq P \leq N^2$:
Orreries and Replicated Orreries}}}




One standard approach towards solving the N-body 
has $N$ processors connected as a ring.  Processor $i$ accumulates the
forces on body $i$.  At each time step, $N$ of the pairwise
interactions are computed, then the appropriate data are passed around
the ring so that in $N$ steps all $N^2$ interactions have been
computed and accumulated into the correct processor, see Figure 1.
This particular parallel implementation of the N-body algorithm has
been referred to as the ``digital orrery'' \cite{sussman}
as the ``systolic loop"\footnote{The use of the term ``systolic loop"
has also been used to refer to the indirect algorithm in \cite{B-H}.}
\cite{windemuth}, and
as ``systolic loop double" \cite{raine}.

As a Latin Square, this is the addition table of the
Abelian cyclic group of order $N$:
the square is circulant, and the $(i,j)$th
entry is essentially $i-j$ modulo $n$.
This case is illustrated in Figure \ref{CC} as Case III,
which shows that one row of interactions is computed in each timestep.


With more than $N$ processors one can replicate the orrery,
perform a different piece of the computation on each copy, and
accumulate the summands for the forces at the end \cite{barnes1}.  For
example, with $2N$ processors the computation takes place on two rings
of processors (see Figure 2) each of which passes its data only $N/2$
times.  More generally, with $MN$ processors we use $M$
rings of $N$ processors each.  At each computation step $MN$
interactions are computed; $N/M$ steps with cyclic message passing
of the data are required to compute the $N^2$ interactions.  Thus the
arithmetic complexity of the algorithm is $O(N/M)$.

An efficient processor topology comes from mapping the replicated
orrery onto a  3-dimensional cube of processors as observed in
\cite{MT}.  We assume the cube is embedded in a larger dimensional
hypercube in such a way that we can distribute a value to, or sum the
contents of, $M$ processors along an axis of the
cube with $O(\log M)$ communications \cite{JH2}.
 
To understand how that data mapping works, assume that the hypercube
is configured as an $N/M$ by $M$ by $M$ cube (see Figure 3).  The data
is initially resident on the front face of the cube and is sent to the
top of the cube.  The data on the front face is spread back, in $log
M$ time, along the $y$-dimension of the cube, and the data on top is
spread down along the $z$-dimension of the cube.  The data is now
distributed through the cube so that each plane parallel to the
$xz$-axis contains a copy of the orrery.  The computation proceeds
by computing a pairwise interaction and doing  1-dimensional cyclic
message passing along the $x$-axis.  This process is repeated $N/M$
times until all $N^2$ pairwise interactions have been completed.  Now
the corresponding results for each individual body are accumulated, in
$log M$ time, along the $y$-axis back to the front face of the cube.
In most cases we are restricted to $M$ and $N$ being powers of two.
This ensures the correct wrap-around for cyclic message passing, and
efficient logarithmic time spreading and accumulating.  Observe also
that as the number of replications, $M$, increases from 1 to $N$, and
thus the number of processors increases from $N$ to $N^2$, the cube
grows out from the $x$-axis and then projects onto the $yz$-plane.
Thus, the dimensions of the cube can be changed to accomodate 
a family of parallel $N$-body solvers where the number of processors
ranges from $N$ to $N^2$.  In particular, with $MN$ processors,
the solver runs in $O(N/M) + O(log M)$ time.
Note that this gives the correct
complexity bounds for the extremes, i.e., $O(N)$ for the case $M=1$,
and $O(1) + O(log N)$ when $M=N$.

Since $P$, the number of processors, is usually fixed given some
particular piece of hardware, there is another useful way to view the
complexity of the algorithm.  Note that $P = MN$, thus, the complexity
bound $O(N/M) + log(M)$ becomes $O(N^2/P) + log(P/N)$.

We remark that the algorithms described above are essentially
grid algorithms.  Once the data is set up, all communication
is performed through cyclic message passing.  Though this
is the easiest conceptually, there are other communications
ideas that make far better use of the hypercube bandwidth.
In the next section, we describe a multiwire algorithm that
makes optimal use of the bandwidth when $N>P$ or to be more
precise, when $N$ is a multiple of $dP$, where $d$ is the dimension
of the hypercube.  There are other communications primitives 
that make better use
of the hypercube bandwidth in the orrery, 
such as the ``multinode broadcast under the MLA assumption''
described in \cite{Bertsekas-etc} and an equivalent algorithm
related to us by Ho \cite{ho}.
We have chosen not to persue this here.

Figure \ref{EE} depicts the Latin Square for a replicated
orrery case for $N=8$ and $P=16$.  Notice how the cube
unfolds onto the computational window.

\begin{figure}
\vspace{0.5in}
Figure will go here
\vspace{0.5in}
\caption{Latin Square for replicated orrery, $N=8$,$P=16$}
\label{EE}
\end{figure}




\section{{\protect\large {\protect\bf $P < N$: Multiwire Algorithm}}}

\subsection{{\protect\large {\protect\bf The Approach}}}

In the case where $P < N$ the data mapping is much simpler.  We divide
the $N$ bodies evenly over the $P=2^d$ processors of a $d$-dimensional
hypercube. For simplicity, assume the number of bodies in each node,
$N/P$, is an integer multiple of $d$.  To compute the $N^2$
interactions, the data for each body must be transmitted to all of the
other processors.  We want to do this using as much of the
communications bandwidth of the hypercube as possible.  If the
processors have enough memory, then this transmission can be performed
initially, and afterwards, all the arithmetic can be performed
locally.  This type of communication in a hypercube has been called
{\bf all to all broadcasting}
\cite{Johnsson-Ho} and {\bf multinode broadcast} \cite {Bertsekas-etc}.

We do not follow this strategy exactly.   We can
interleave pieces of the computation with the communications and have no
need to store the data for all the bodies at each node.  Thus, we can use
the idea of ``rotated and translated Gray codes" described in
\cite{Johnsson-Ho}.  Although this method has ``limited potential for
pipelining'', it is just right for the $N$-body problem.
A Gray code is rotated to construct as many conflict-free paths as
possible leaving a node, visiting all others, and returning to its
starting point.  By translating the starting point of each of these
paths to every node of the hypercube, the data in each node circulates
to every other one in as short a time as possible.

More precisely the rotated and translated Gray codes produce $d2^d$
time-wise edge disjoint Hamiltonian paths through the hypercube, $d$
of them\\
starting at each node.  A {\bf Hamiltonian path} in a graph
visits all the nodes in the graph only once.  {\bf Time-wise edge
disjoint} means that, although the paths may share edges of the
hypercube, no two paths traverse the same edge in the same direction
on the same communication step.  Consequently, assuming
bidirectionality, there is no contention for communication channels.
The data is circulated on these paths, but not retained, avoiding the
memory waste.  Furthermore, from a programming point of view, this
approach is very simple and one can take advantage of the symmetry of
the hypercube.  The combination of the simplicity and the optimal use
of communication bandwidth makes this communication pattern one of the
fastest which a hypercube architecture can perform.





\subsection{Complexity Analysis}
The complexity analysis of the algorithm is straightforward.  The
computation time, $T_{arith}$, scales like ${N^2}/{2^d}$, where $N$ is
the total number of bodies and $d$ is the dimension of the hypercube.
The communication time, $T_{comm}$, scales like $\lceil \frac{N}{2^dd}
\rceil(2^d-1)$; where $\lceil x \rceil$ denotes the ceiling of $x$,
i.e., the smallest integer greater than or equal to $x$. 
Since there are
$N/{2^d}$ bodies per node and $d$ wires leaving each node, the first
factor in $T_{comm}$ represents the number of communication operations
needed to empty a node of data.  The second factor is the
length of the path through the hypercube that the data must traverse.

Given a value of $N$ such that $\frac{N}{2^d} \gg {d}$, $T_{comm}$ is 
at first a linearly decreasing function of $d$. Then, as the number 
of dimensions increases, the number of  
required communication operations will eventually decrease to one and the 
second factor in the expression of $T_{comm}$ starts to dominate.
Because at least one such communication operation must be performed
even if there are less than $d$ data items at a node we get,
for a fixed value of $N$,
\begin{center}
\[\lim_{d \rightarrow \infty} T_{comm}(d) \simeq 2^d. \] 
\end{center}
The function $T_{comm}$ therefore exhibits a minimum which roughly
occurs when the number of bodies per node is equal to the size of the 
bandwidth $d$. As the number 
of bodies per node decreases
 one is making less and less effective use of the communications
bandwidth of the hypercube and $T_{comm}$ becomes a more significant
portion of the total execution time.  In fact, when $\frac{N}{2^d}
 \leq 1$,
i.e., the number of bodies is much smaller than the number of
processors, then the replicated orrery approach becomes more cost
effective.  Note that in the boundary case, $N=P$, the total
computational complexity reduces to $O(N)$ which agrees with the
orrery algorithm.

\subsection{Rotated and Translated Gray Codes}

A $d$-dimensional hypercube is a graph with $2^d$ nodes labeled by the
$d$-bit binary representation of the integers $0$ to $2^d - 1$.  Each bit
in the $d$-bit representation is associated with a different dimension of
the cube.  There is an edge between two nodes $i$ and $j$ in the hypercube
if and only if their binary representations differ in only one bit.  We can
think of the edges as traversing different dimensions of the hypercube.

A {\bf Gray code} is a circuit of all binary $d$-tuples, such that only one
coordinate position changes at each step.  Thus, a Gray code represents a
Hamiltonian path on the hypercube.  There are many ways to construct Gray
codes, but the most famous is the binary reflected Gray code
\cite{graycoderef}.  To simplify the terminology, we will refer to this code
as {\bf the} Gray code.  We define the {\bf transition sequence}
\cite{graycoderef} corresponding to a Gray code as the list of positions
which change at each step of the code or the list of dimensions traversed
on the hypercube at each step of the circuit.  In Figure 4, we show the
binary reflected Gray code for $d=3$, the corresponding transition
sequence, and the associated Hamiltonian path on the cube.


Given the Gray code beginning at $0$, we can ``translate" it to node $i$ by
taking the {\it exclusive or} of $i$ with the $d$-tuples in the original
code.  The transition sequence of the translated Gray code is identical to
the original.  Figure 5 shows the translated code starting at node $5 =
101$.  Furthermore, given the Gray code beginning at $0$, we can ``rotate"
it by performing circular shifts of the bits in the original code.  In the
tables below, we list the three rotated Gray codes and the corresponding
transition sequences when $d=3$.  Notice that in general, the transition
sequence of the rotated code $i$ is obtained by adding $i$ modulo $d$ to
the transition sequence of the binary reflected Gray code.

\vspace{0.2in}
\begin{center}
\begin{tabular}{|l|cccccccc|}  \hline
\multicolumn{9}{|c|}{Rotated Gray Codes} \\ \hline
Code 0  & 000 & 001 & 011 & 010 & 110 & 111 & 101 & 100 \\
Code 1  & 000 & 010 & 110 & 100 & 101 & 111 & 011 & 001 \\
Code 2  & 000 & 100 & 101 & 001 & 011 & 111 & 110 & 010 \\ \hline
\end{tabular}
\end{center}
\begin{center}
\begin{tabular}{|l|ccccccc|}\hline
\multicolumn{8}{|c|}{Transition Sequences} \\ \hline
Code 0 & 0 & 1 & 0 & 2 & 0 & 1 & 0  \\
 Code 1 & 1 & 2 & 1 & 0 & 1 & 2 & 1  \\
Code 2 & 2 & 0 & 2 & 1 & 2 & 0 & 2  \\ \hline
\end{tabular}
\end{center}
\vspace{0.2in}


By taking all of the $2^d$ possible translations of the $d$ rotated
Gray codes beginning at $0$, we can generate a set, $H$, of $d2^d$
Hamiltonian paths on the hypercube.  This set $H$ is more elegantly
described as the paths generated by applying a group of $d2^d$
symmetries (rotations and translations) of the hypercube to the binary
reflected Gray code\footnote{The use of the binary reflected Gray code
is merely for convenience of implementation; any other Hamiltonian
path will suffice.} 
starting from node $0$.  These symmetries are the
subgroup of all the hypercube symmetries that rotate the dimension
numbers modulo $d$.  It is straightforward to verify that two distinct
symmetries of a hypercube transform a directed edge into two distinct
directed edges.  It follows that the collection of paths $H$ have the
nice property of being time-wise edge disjoint: i.e., if $j \in {1,2,
\ldots,2^d-1}$, then there is no overlap among the $j$th edges of the paths
in $H$. This implies that if we have $d$ packets of
information in each node of the hypercube, we can circulate the data, one
packet to a Hamiltonian path, in $2^d-1$ communication steps.  Furthermore,
we can perform the arithmetic on the data in between each communication
step eliminating the need for data storage.

\section{{\protect\large {\protect\bf CM-2 implementations}}}


\subsection{Implementation of the replicated orrery}

Using the fieldwise programming model for the CM-2, we typically have
$N \leq P \leq {N^2}$, thus the replicated orrery $N$-body algorithm
is most appropriate.  We describe the details of the implementation
here.

In configuring the CM-2 as $k$-dimensional grids, we are required to
restrict the lengths in each dimension to a power of $2$. 
Thus, the data mapping of the replicated orrery onto a 3-dimensional grid
proceeds as follows. 
Pick $n$ such that $2^n$ is the first power of $2$ greater than or equal 
to $N$. 
Then we will actually use an orrery of size $2^n$ with only $N$ valid entries. 
We are assuming that $N \leq  P \leq {N^2}$, so we can write $P = 2^ {n+m}$, 
where $m \leq n$.
We will want to replicate the orrery $M = 2 ^ m$ times. 
We configure a $2^ {n+m}$ processor CM-2 as an $2^ {n-m}$ 
by $2^m$ by  $2^m$
cube ($N/M$ by $M$ by $M$). 
This configuration allows us to use versions of the parallel prefix 
{\bf scan}
operations to spread or accumulate data along any of the 3 coordinate directions in
logarithmic time, and ensures the correct wrap-around for grid communications.
The cube mapping, complete with sends, spreads, and accumulates, is drawn 
in Figure 3. 

There will be two copies of the data for each body, one dynamic and
one static. 
The static data is initially resident on the front face of the cube. 
It is sent to the top of the cube by a low density 
{\bf send}.
The data on the front face is spread back along the $y$-dimension of the 
cube, and the data on top is spread down along the $z$-dimension of the cube. 
The data is now distributed through the cube so that each plane parallel to 
the $xz$-axis contains a copy of the orrery. 
The computation now proceeds by computing a pairwise interaction according to 
the relevant force law and doing a 1-dimensional wrap-around 
{\bf NEWS}
communication along the $x$-axis.
This process is repeated $2^ {n-m}$ times until all $N ^ 2$ pairwise
interactions have been completed. 
Now the contributing forces for each individual body are accumulated by doing 
a 
{\bf plus-scan}
along the $y$-axis back to the front face of the cube.  

Because of the power of 2 restrictions noted above, the time required
to do the computation for any value of $N$ between $2 ^ n$ and
$2 ^ {n+1}-1$ is the same. 
As $N$  increases above the next power of 2,
we can reconfigure the grid to the appropriate size.  


\subsection{Implementation of the multiwire algorithm}

Using  the slicewise model of the CM-2,  we are generally in the case
where $P < N$, and so we use the multiwire algorithm described above.


If we have $d$ words at every node, and if we wish to circulate the data
using the Hamiltonian paths described in the  Section 3, one might
expect that we would either need to attach some labeling information to be
transmitted with the data or else use a look up
table that specifies the dimension  over which data needs to be sent.  Neither
of these operations is without cost.  In our implementation, as we
shall now describe, we avoid both of these types of methods.  In fact, the
information in the transition sequence for the binary reflected Gray code
is all that is needed.

Using the slicewise model, the CM-2 architecture allows for the loading of
as many as $2d$ words to be sent over the $d$ dimensions simultaneously,
two words in each dimension.  Typically, one loads one set of $d$ words for
the $d$ dimensions followed by another set of $d$ words.  For simplicity,
we will concentrate on one set with the understanding that the other set
follows exactly the same pattern.  Another important feature of the CM-2
architecture is that it is particularly efficient to load consecutive
words, or more generally, words with a constant stride, to be sent over
consecutive dimensions.

To illustrate how the CM-2 works, let us first assume that we only have
one word in each processor and each of these words will follow the
Hamiltonian path specified by the transition sequence of the Gray code:
${g}=\{0,1,0,2,0,1,\ldots\}$. Further assume we have only one memory
location (location 0) in each node where the word is stored.  Then the
algorithm is simply:

\boxer{
\begin{eqnarray*}
\lefteqn {\mbox{For }  k=1,2,\ldots,2^d-1} \\
& & \mbox{In each node, the data at location 0 is} \\
& & \mbox{1. Loaded to be sent over dimension } g_k \\
& & \mbox{2. Transmitted} \\
& & \mbox{3. Stored in location 0} \\
\end{eqnarray*}
}

Now, if we have $d$ words in each node at memory locations $0$ through
$d-1$, we can quite readily make sure that the word at location $j$ in each
node follows the $j$th transition sequence.  The method is to store this
word at location $j$ as it passes through each node like a traveler who
prefers the same hotel in each city.  Then, at step $k$, we send the data
in each node at location $j$ over dimension $g_k + j \pmod{d}$. To take
advantage of the Connection Machine's ability to load consecutive words, we
in fact use the algorithm in the box below.

\begin{table}
\boxer{
\begin{eqnarray*}
\lefteqn {\mbox{For }  k=1,2,\ldots,2^d-1} \\
& & \mbox{In each node} \\
& & \mbox{1. Data at locations $0$ through $d-g_k-1$ is loaded to} \\
& & \mbox{\ \ \ be sent over dimensions $g_k$ through $d-1$ respectively.} \\
& & \mbox{2. Data at location $d-g_k$ through $d-1$ is loaded to} \\
& & \mbox{\ \ \ be sent over dimensions $0$ through $g_k-1$ respectively.} \\
& & \mbox{3. Data is transmitted.} \\
& & \mbox{4. Data received over dimensions $g_k$ through $d-1$ is}\\
& & \mbox{\ \ \ stored at locations $0$ through $d-g_k-1$ respectively.}\\
& & \mbox{5. Data received over dimensions $0$ through $g_k-1$ is}\\
& & \mbox{\ \ \ stored at locations $d-g_k$ through $d-1$ respectively.}\\
\end{eqnarray*}
}
\end{table}

Notice that steps 2 and 5 are vacuous if $g_k=0$, which is half of the
time.  The actual code to perform this operation is not much longer than
the boxed pseudo-code.

In the real algorithm, we assume there is a multiple of $2d$ bodies in each
node, and each body may contain more than one word of information.  In the
example described in the next section, there are two words of coordinate
information and also a coefficient, giving three words for each body.
Since we can transmit $2d$ words simultaneously, we loop over all the data
until it is all transmitted.  After the data is transmitted, we perform a
computation such as the one described in the next section, and then we
repeat for $k=1, 2,\ldots, 2^d-1$.

The ideas we presented in this section have now been coded up
as a Connection Machine Scientific Software Library Primitive.
See \cite{BJ} for details.


\section{A sample application}


\subsection{Description}

As an application example, we computed the velocities of a collection of
interacting point vortices in two dimensions. This $N$-body problem is at
the heart of the vortex method to solve the Navier-Stokes equation for
incompressible viscous flows \cite{sethian}. A Lagrangian perspective of
the vorticity-stream function formulation of the Navier-Stokes equation
leads to the introduction of discrete vortex elements which interact
according to a Biot-Savart type force law.  To avoid potential numerical
instability caused by very close encounters of particles, finite size
vortex elements should be used.  However, in these timing runs, we have
considered the interaction of point vortices.

In two dimensions, the $(u,v)$ velocity components of a point
vortex at position $(x,y)$ are given by the following formulae:

\begin{center}
\[u(i) = \frac{1}{2\pi} \sum_{j=1}^{N} \frac{c(j)\Delta y(i,j)}{\Delta
r(i,j)^2}, \;\;\;\; i=1,...,N\]
\[v(i) = \frac{-1}{2\pi} \sum_{j=1}^{N} \frac{c(j)\Delta x(i,j)}{\Delta
r(i,j)^2}, \;\;\;\;i=1,...,N\]
\end{center}
\begin{eqnarray*}
 \mbox{with}  & \Delta x(i,j) = x(j) - x(i) & \\
&   \Delta y(i,j) = y(j) - y(i) & \\
 & \Delta r(i,j)^2 = \Delta x(i,j)^2 + \Delta y(i,j)^2  & \\
\mbox{and} & c(j) = \mbox{ circulation of point vortex $j$}. &
\end{eqnarray*}

\subsection{Timings and analysis}

%!! Mention that these two methods are not comparable.

We performed several calculations for both the replicated orrery and
the multiwire algorithm with a varying number of particles on various
machine sizes. Timing results are shown in the two tables below.
Timings for the replicated orrery algorithm were run on 4K and 8K
processor configurations.  The multiwire timings were obtained on CM-2
configurations with $2^{d+5}$ processors for $d = 4,5,7,8,9$, i.e.,
with $2^d$ floating point nodes.  Typically, 14,000 vortices interact
in 2 seconds on a 512 floating-point node machine (16K processors)
using the multiwire algorithm, and 17 seconds on an 8K processor
machine using the replicated orrery.  For an equivalent number of
bodies on a fixed configuration, the multiwire algorithm outperforms
the replicated orrery by a factor of 3 to 5.  Some loss in performance
of the replicated orrery algorithm is due to the necessity of
transposing fieldwise data back and forth for the floating point
units. Counting the divide operation as one operation, the execution
rate of the multiwire algorithm is about 5.2 Gflops on a fully
configured CM-2 (2048 nodes and 64K processors).
\footnote{On the 32-bit floating point unit of the CM-2, a divide
operation requires the explicit implementation of two Newton-Raphson
iterations of the initial value which is provided by an internal
look-up table.  This amounts to a total of six atomic operations which
we count as one}.


\vspace{0.2in}
\begin{center}
\begin{tabular}{|c|c|c|c|}  \hline
\multicolumn{4}{|c|}{Replicated Orrery} \\ \hline
Processors & Bodies & $T_{total}$ & $T_{comm}$\\ \hline
4K & 1024 & .15	& 40\% \\
   & 2048 & .6	& 30\% \\
   & 4096 & 2.5	& 30\% \\
   & 8192 & 8.5	& 24\% \\
   & 16384 & 30 & 20\% \\ \hline

8K & 1024 & .08	& 44\% \\
   & 2048 & .3	& 30\% \\
   & 4096 & 1.3	& 28\% \\
   & 8192 & 5.1	& 29\% \\
   & 16384 & 16.7 & 24\% \\ \hline
\end{tabular}
\end{center}


\begin{center}
\begin{tabular}{|c|c|c|c|}\hline
\multicolumn{4}{|c|}{Multiwire Algorithm} \\ \hline
Processors  & Bodies & $T_{total}$ & $T_{comm}$  \\  \hline
512 ($d=4$) & 1792 & 1.0 & 4\%  \\
	    & 3584 & 4.0 & 2\%  \\ \hline
1K ($d=5$) & 1792 & .52 & 7\%  \\
	   & 3584 & 2.0 & 4\%  \\ 
	   & 7168 & 7.8 & 2\%  \\ 
	   & 14336 & 30.8 & 1\%  \\ \hline
4K ($d=7$) & 1792 & .16 & 15\%  \\
	   & 3584 & .54 & 9\%  \\ 
	   & 7168 & 2.0 & 5\%  \\ 
	   & 14336 & 7.8 & 3\%  \\ \hline
8K ($d=8$) & 1792 & .11 & 28\%  \\
	   & 3584 & .3 & 10\%  \\ 
	   & 7168 & 1.07 & 6\%  \\ 
	   & 14336 & 4.04 & 3\%  \\ \hline
16K ($d=9$) & 3584 & .2 & 50\%  \\ 
	   & 7168 & .6 & 16\%  \\ 
	   & 14336 & 2.1 & 9\%  \\ \hline
\end{tabular}
\end{center}
\vspace{0.2in}


Keeping the number of bodies constant, the time to compute the
velocity components appears inversely proportional to the number of
floating-point nodes of the hypercube.  This is the kind of behavior
we expect given the complexity analysis of both algorithms.  Note
that, for a fixed number of processors, the proportion of time spent
in the communications for the replicated orrery decreases very slowly
as the number of bodies increases.  In the multiwire algorithm the
proportion of communication time is much smaller and the decrease is
more rapid.  We get a more accurate picture of what is going on in the
multiwire algorithm from Figure 6 which plots the actual speedup of
the implementation against the optimal speedup of $2^d$ (the number of
processors).  We calculated the speedup by taking the ratio of
$T_{arith}$ times the number of processors and the total execution
time of the calculation.  Note that again the curves behave as
predicted by the complexity analysis.  For a fixed hypercube
dimension, as the number of bodies increases the speedup approaches
optimality.  In particular, over the range of test cases we ran,
$T_{comm}$ varies from a few percent to fifty percent of the total
execution time as the number of bodies approaches the number of
processors.

\section{Acknowledgements}
We would like to thank Washington Taylor for his help in the first
implementation of the replicated orrery and Adam Greenberg for running
a set of timings of the replicated orrery for us.

\def\\{\medskip\noindent\hangindent .3in{}}
\setlength{\baselineskip 11pt}
\frenchspacing
\begin{thebibliography}{99}

\bibitem{anderson} Anderson, C.R.,  1986, ``A method of local
corrections for computing the velocity field due to a distribution of
vortex blobs'', {\em SIAM J. Numer. Anal.} {\bf 22}, 413-440.

\bibitem{apel} Appel, A.W., 1985, ``An efficient program for many-body
simulation'', {\em SIAM J. Sci. Stat. Comput.} {\bf 6}, 85-103.


\bibitem{sussman} Applegate, J.F., M.R. Douglas, Y. Gursel, P. Hunter,
C. Seitz, and G.J. Sussman, 1985,
``A digital orrery",
{\em IEEE Trans. Comput.} {\bf 84}, 822.

\bibitem{barnes1} Barnes, J., ``Multiple-timestep N-body algorithms
for parallel computers'', 1987, Institute for Advanced Study,
Princeton, (preprint--never published).

\bibitem{barnes} Barnes, J., and P. Hut, ``A Hierarchical
$O(NlogN)$ force calculation algorithm'', 1986, {\em Nature} {\bf
324}, 446-449.


\bibitem{B-H}
Barnes, J. and Hillis, W.D. Programming a highly parallel computer. {\sl
Nature} {\bf 326}, (1987), 27--30.

\bibitem{Bertsekas-etc}
Bertsekas, D.P., Ozveren, C., Stamoulis, G.D., Tseng, P., and Tsitsiklis,
1991,
J.N. Optimal communication algorithms for hypercubes. 
{\sl Journal of Parallel and Distributed Computing}
{\bf 11}, 263--275. 

\bibitem{BBC}
Beck, A., Bleicher, M.N., Crowe,  D.W.,
{\sl Excursions into Mathematics}, Worth Publisher, Inc., N.Y., 1976.


\bibitem{parpref-ref} 
Blelloch, G.E. ``Scans as primitive parallel operations",
{\sl Proc. Int. Conf. on Parallel Processing}, 1987.

\bibitem{BEM}
Brunet, J.P., Edelman, A., Mesirov, J.P.,
{\sl Proceedings Supercomputing '90}.
IEEE Computer Society Press, Los Alamitos, CA.
748--752.

\bibitem{BJ}
Brunet, J.P., Johnsson, L.,
``All-to-all broadcast and applications on the Connection
Machine", TMC-189/BA.

\bibitem{LATIN}
D\'{e}nes, J., Keedwell, A.D.,
{\sl Latin Squares and their Applications},
Academic Press, N.Y. 1974.



\bibitem{Fox} Fox, G.C., Hipes, P., and Salmon, J. Practical parallel
supercomputing: examples from chemistry and physics. {\sl Proceedings
Supercomputing '89}. ACM Press, Reno, NV, 1989,  58--70.

\bibitem{greengard} Greengard, L., and V. Rokhlin, 1986, ``A fast
algorithm for particle simulations'', 
{\sl J. Comp. Phys} {\bf 73}, 325--348.


\bibitem{heller}
Heller, H., Grubm\"{u}ller, H., Schulten, K.,
``Molecular dynamics simulation on a parallel computer"
{\sl Molec. Simul.} {\bf 5}, 133--165.


\bibitem{ho} Ho, C.T., personal communication, 1991.

\bibitem{hockney} Hockney, R.W., and J.W. Eastwood, ``Computer
simulation using particles'', McGraw-Hill, New York, 1981.

\bibitem{JH2}  Johnsson, S.L., ``Communication efficient basic linear algebra
computations on hypercube architectures'', 1987, {\sl Journal of
Parallel and Distributed Computing} {\bf 4}, pp. 133-172.

\bibitem{Johnsson-Ho}
Johnsson, S.L. and Ho, C.T. Optimum broadcasting and personalized
communication in hypercubes. {\sl IEEE Transactions on Computers} {\bf
38}, (1989), 1249--1268.

\bibitem{Katzenelson}
Katzenelson, J., 1989, ``Computational structure of the N-body problem'',
{\sl SIAM J. Sci. Stat. Comput.} {\bf 10}, 787--815. 

\bibitem{MT} Mesirov, J.P., and W. Taylor, 1988, ``A Family of
N-body Solvers'', unpublished.

\bibitem{raine},
Raine, A.R.C., Fincham, D., Smith, W.,
``Systolic loop methods for molecular dynamics simulation
using multiple transputers",
{\sl Computer Physics Comm.} {\bf 55}, 13--30.

\bibitem{graycoderef}
Reingold, E.M., Nievergelt, J., and Deo, N. {\sl Combinatorial Algorithms},
Prentice-Hall, Englewoods Cliffs, N.J., 1977.

\bibitem{sethian} Sethian, J.,  J. Brunet, A.
Greenberg, and  J.P. Mesirov, 1989, ``Two-dimensional, viscous,
incompressible flow on a massively parallel processor'', (to appear
{\em Journal of Comp. Physics}).

\bibitem{smith}
Smith, W.,
``Molecular dynamics on hypercube parallel computers",
{\sl Computer Physics Comm.} {\bf 62}, 229-248.

\bibitem{windemuth} Windemuth, A., Schulten, K.,
``Molecular dynamics simulation on the Connection Machine",
{\sl Molec. Simu.}, {\bf 5}, 353-361.


\bibitem{Johnsson}
Zhao, F. and Johnsson, S.L. ``The parallel multipole method on the Connection
Machine''. 
{\sl SIAM J. Sci. Stat. Comput.} {\bf ??}, 1991.


\end{thebibliography}

\end{document}

