\documentstyle[12pt]{article}
\title{Implicit $15^\circ$ finite-difference migration on the Connection Machine}
\author{Stewart A. Levin}
\date{\sl draft \rm \today}
\begin{document}
\maketitle
\begin{abstract}
Claerbout's implicit $15^\circ$ finite-difference migration is a logical
candidate for interactive migration programs.  Its unconditional
stability readily allows one to trade accuracy for speed.  Arbitrary
gridded velocity functions can be accomodated.  In addition it is
vectorizable, and hence parallelizable as previously shown by the
author.  This report shows how to implement the algorithm on the
massively parallel Connection Machine and presents timing results.
\end{abstract}
\section{Introduction}

\hspace*{\parindent}
Claerbout's implicit $15^\circ$ finite-difference migration (T15mig)
is a second-order, centered
Crank-Nicholson approximation to the paraxial wave equation
\begin{equation}
{ {\partial^2 P } \over {\partial t \partial \tau } }
= - { v^2 \over 8 }
{ {\partial^2 P } \over {\partial x^2} }
\end{equation}
which is used to migrated dips out to about $25$ or $30^\circ$.
Here $P$ is the zero-offset wavefield, $t$ is time, $x$ is distance
along the surface, and $\tau$ is vertical traveltime.

To discretize (1) for computation, denote the temporal sampling
interval by $\Delta t$, the trace spacing by $\Delta x$, and 
select a vertical traveltime interval $\Delta \tau$ that is some
integer multiple of $\Delta t$.  Define
\begin{equation}
\matrix{P_1 &\equiv& P(x, t, \tau - \Delta \tau)\hfill\cr
    P_2 &\equiv& P(x, t + \Delta t , \tau - \Delta \tau )\hfill \cr
    P_3 &\equiv& P(x, t + \Delta t , \tau )\hfill \cr
    P_4 &\equiv& P(x, t , \tau )\hfill \cr}
\end{equation}
as the four vectors entering in the 2 by 2 finite difference stencil.
Then the stencil may be written
\begin{equation}
{\bf T_1 } [ P_4 + P_2 ] = { \bf T_2 } [ P_1 + P_3 ]
\end{equation}
where
\begin{equation}
\matrix{\alpha\hfill & = & (v^2 \Delta t \Delta \tau )/(32 \Delta x^2) \cr
{\bf T }\hfill & = & \hbox{\sl tridiag \rm (-1,2,-1)} \cr
{\bf T_1 }\hfill & = & {\bf  I + T } ( \alpha - \beta ) \cr
{\bf T_2 }\hfill & = & {\bf I - T } ( \alpha - \beta ) \cr}
\end{equation}
and $\beta$ is a tunable parameter used to improve accuracy and
typically given a value around 1/6'th.

For migration, the unknown in the stencil is $P_4$.  Referring to
Figure 1, the computations are started by initializing the first
column ($\tau = 0 $) to be the stacked section and the bottom of each column
($\tau > t_{\hbox{\sl max}}$) to zeros.
The migrated image appears on the right-hand
``diagonal'' edge of the grid.
Computational efficiency is achieved by vectorization along
$45^\circ$ diagonals as depicted in the figure.
Vector lengths vary from top to bottom, with an average in
real life of about 100.
Using this parallelism, the author achieved full vectorization
of the computation on the Cray and Convex and a 55-fold speedup
on a 64 node NCUBE parallel processor array.

\section{The Connection Machine}

\hspace*{\parindent}
The Connection Machine (CM-2) at DRL is a massively parallel computer with
16K processors, allocable in multiples of 8K,
each having a small amount of local memory.
These processors are usually configured as some Cartesian
grid with fast data links between neighboring nodes of the grid.
General interprocessor communication is also supported,
though at a greatly reduced performance.
Control is handled through an attached host, a VAX 6310 here at DRL, and
the processors are conceptually a ganged, fully synchronized (SIMD)
array.  This means that each instruction issued by the host is
executed by every active processor in the grid.
The major characteristics of the CM-2 that affect implementation of T15mig are:
\begin{description}
\item[parallelism] The algorithm needs to achieve a parallelism 
of at least 8K in order to use the CM efficiently.  For a serial
program (parallelism of 1), the CM-2 is about as fast as an IBM AT.
\item[locality] The algorithm should use nearest-neighbor grid
communications (NEWS) when possible.  General communication can
be easily 100 times slower.
\item[alignment] Computational arrays should be aligned
whenever possible so that corresponding
nodes are in the same processor thus avoiding communication altogether.
\item[power-of-two] All axes of Cartesian grids are required to be powers
of two.  A significant fraction of the processors may be wasted due to padding.
Computations are also further slowed if processors have to be disabled
to accomodate the padding.
\end{description}

\section{T15mig}

\hspace*{\parindent}
It is clear that the Cray and NCUBE implementations of T15mig will
not map ideally to the CM-2.  Indeed, these fall short by a factor of
100 in the amount of parallelism needed to take advantage of the
CM-2.  The locality of data and its alignment to the
computational grid was a small issue on the Cray and NCUBE.  Additionally,
padding requirements are significantly more restrictive on the CM-2.


\subsection{Data and Computational Grids}

\hspace*{\parindent}
To implement T15mig on the CM-2, we need a computational grid to hold the
working arrays $P_j$ along with the various tridiagonal coefficients,
an input data grid to hold the unmigrated time section and an output
data grid to hold the migrated time section.
As in previous implementations
for the Cray and the NCUBE, it is possible to save space and equivalence
the input and output data grids.

In T15mig the computational grid is
smaller than the input/output grid, typically by a factor
of $jtau$ $\equiv \Delta \tau / \Delta t$ $ \approx 10$.
It is best to fold the input grid
into several sheets as illustrated in Figure 2
so that each sheet is the same size as the computational
grid.
This permits us to align the input and computational grids
and improve performance.  Of course one {\it could } align
the grids by leaving the input section as one sheet
and padding the computational grid to the size of the input.
However this layout leads us to perform about {\it jtau\/} times as many
calculations as the multi-sheeted layout.

With this alignment, data movement into the computational grid is
now local, i.e. nearest neighbor moves and in-processor copies.  Since
data must feed from the input grid through the computational grid
and back to the output grid, we need to decide which grids to actually move.
That is, should the input grid shift forward in order to feed the
computational grid or should the computational grid shift backwards
in order to consume the input grid?  This depends upon how much data
has to be shifted.  For T15mig, we will find that one of
the $P_j$'s will need to be shifted in the former and two in the latter.
If the computational
grid shifts backwards,
we need to shift the attached tridiagonal coefficients
too.
If the input data shifts forward, we need to determine how many sheets
have to be shifted at each computational step.  We defer answers to
these questions until later.

In terms of padding, the computational grid still must be padded to
powers of two.  Accepting this, one can turn this to advantage by
reducing {\it jtau\/} to preserve the same padded dimensions but achieve
greater accuracy.

\subsection{Tridiagonal Equations and Parallelism}

\hspace*{\parindent}
We've now got a handle on how locality, alignment and power-of-two
should be incorporated in a CM-2 version of T15mig.  This leaves
parallelism to be considered.  Here we must diverge from the Cray
or NCUBE algorithms which are built around the parallelism of about 100
as indicated in Figure 1.  We need to achieve another 100-fold improvement
in parallelism in order to use the CM-2 effectively.  The bottleneck
is in the tridiagonal solve of $\bf T_1$ on the left hand side of (3).
The answer is to employ a more parallel method of solving
the tridiagonal systems.
Actually, just as in previous T15mig implementations,
it is advantageous to prefactor as much as possible of
the tridiagonal matrices $\bf T_1$ in (3).  The standard
product of
lower bidiagonal, diagonal and and upper bidiagonal matrices,
as typified by the well-known $LDU$ factorization, is one
reasonable approach.
This formulation reduces the
number of arithmetic operations in the individual solves by about
a factor of two.  However ingenuity may be needed to compute
these factors on the CM -- it might be better to use the front-end
host.

Our discussion will center on parallel algorithms that implements the first
order recursion
\begin{equation}
y_i = x_i + a_i y_{i-1}
\end{equation}
that arises in the $LDU$ factored version of the tridiagonal matrix
$\bf T_1$.  We will also cover exceptions where alternate prefactoring
is superior.
There are a number of choices available in the literature.

\subsubsection{Cyclic Reduction}

\hspace*{\parindent}
The first choice is cyclic reduction due to \cite{b&g:pois}.
This is derived by
substituting for $y_{i-1}$ in (5) to yield a new recursion
\begin{equation}
\matrix{
y_i &=& (x_i + a_i x_{i-1} ) + (a_i a_{i-1} ) y_{i-2}\hfill\cr
    &=& x_i' + a_i' y_{i-2}\hfill\cr }
\end{equation}
which links every other index $i$.  Proceeding recursively with
off-end values defined as $y_i = 0 $ and $a_i = 1$, we get
relations between $y_i$ and $y_{i-4}$, $y_{i-8}$, $y_{i-16}$
and so forth.  Eventually every $y_i$ reduces to a relation
involving a known zero value of $y$ and the recursion terminates.

The number of stages in this reduction is at most
$\lceil {\log}_2 nx \rceil$
and may be substantially fewer if only approximate answers are required.
In particular, if the $a_i$'s have smaller than unit magnitude, the
coefficients of $y_{i-2^k}$ diminish geometrically with $k$.  In T15mig
it behooves us to monitor these off-diagonal coefficients.  Experience
shows that it is permissible to ignore coefficients less than $10^{-3}$
and that in practice this generally terminates the recursion after
3 or 4 stages, saving significantly over a full ${\log}_2 nx$ reduction
and making the cyclic reduction arithmetically competitive with ordinary
elimination
while allowing $nx$ times as much parallelism as the conventional
algorithm.

In T15mig it is not useful to create the conventional $LDU$ factor
prior to applying cyclic reduction.  Instead we apply cyclic reduction
directly to the original tridiagonal equations once and store the coefficients
for later reapplication.  This is equivalent to, and faster than,
precomputing the
coefficients in the bidiagonal cyclic reduction of $L$ and $U$.

This method has been coded in the initial version of T15mig.  For
512 traces of 1024 time points, and a $\Delta \tau$ of 8 samples,
factoring takes a second or so, and migration takes 12 seconds.
This is 3 times slower than the Cray X-MP for this size data.
Of the 12 seconds, 9 seconds is spent solving tridiagonal equations,
8 seconds of which is devoted to power-of-two communication.


\subsubsection{Recursive Doubling}

\hspace*{\parindent}
A close relative of cyclic reduction, recursive doubling (\cite{hss:recdub})
interleaves
reduction and backsubstitution.  The first step is to use equation (5)
to derive $y_1$ from $y_0$.  The reduction (6) is applied to compute
$y_2$ from $y_0$ and $y_3$ from $y_1$ (with a parallelism of 2.)  Another
reduction step gives $y_4, \ldots, y_7$ from $y_0, \ldots , y_3$.
Eventually, after
$\lceil {\log}_2 nx \rceil$
steps, we have computed the complete
solution.  Unlike cyclic reduction, however, there is no intermediate
stopping point at which diagonal dominance can be invoked to terminate
the recursion early.  For this reason, cyclic reduction is preferable
to recursive doubling for T15mig.

\subsubsection{Product Expansion}

\hspace*{\parindent}
A second choice for parallel solution of the linear recurrence (5)
is to rewrite it in the matrix form
\begin{equation}
({\bf I} - {\bf A }) y = x
\end{equation}
where ${\bf A}$ is the lower subdiagonal.  Employing the product expansion
\begin{equation}
{\bf I } = ( {\bf I} - {\bf A } ) \prod_{k=0}^{\lfloor {\log}_2 nx \rfloor}
( {\bf I} + {\bf A}^{2^k} )
\end{equation}
leads to a solution
\begin{equation}
y = \prod_{k=0}^{\lfloor {\log}_2 nx \rfloor}
( {\bf I} + {\bf A}^{2^k} ) x
\end{equation}
for the recurrence.

As with cyclic reduction, the product can be terminated early if the
entries in the shifted diagonal matrix $ { \bf A}^ {2^k} $ are
negligible.  Again, this typically happens around $k = 3$ in T15mig.

\subsubsection{Series Expansion}

\hspace*{\parindent}
A third choice for parallel solution of the linear recurrence (5) is
to use the series expansion
\begin{equation}
y = \sum_{k=0}^{nx-1} {\bf A}^k x
\end{equation}
obtained, for example, by multiplying (9) out.
This is really an explicit construction
of progressive subdiagonals
of the inverse to $ { \bf I }
- { \bf A }$.
For T15mig we can terminate
the series early at around the $ 2^3 $rd term.  This yields a banded,
approximate inverse and shows that each $y_i$ really depends only upon
a few preceding values $ x_i , \ldots , x_{i-M}$ for some small value of $M$.

Compared to the previous two choices, this method has significantly more
arithmetic, using $M$ terms compared to ${\log}_2 M$ terms.  Its uses 
different, generally faster,
communication primitives on the CM-2 than do the other two, so it can
still be of benefit.

\subsubsection{Partitioned Reduction}

\hspace*{\parindent}
Partitioned methods break the original system into blocks, perform
ordinary elimination on the blocks followed by a fixup step to
handle the fill-in that arises at block boundaries.  There are
several variations --- the ones with the most parallelism eliminate
the blocks in parallel, solve for one element in each fill-in block (by,
for example, cyclic reduction), and eliminate the fill-in blocks in
parallel.  This is how the eventual CMSSL
(Connection Machine Scientific Subroutine Library)
tridiagonal routine will operate (\cite{slj:ens,slj:cmsl}).
The size of the blocks
is chosen so that there are at least as many blocks as there are
processors in the CM-2.  This lets us exploit the full parallelism
of the machine.  The number of arithmetic operations in the partition
methods is about twice that of straight elimination and so comparable
to cyclic reduction operation counts.  Indeed, we may think of this
as cyclic reduction by powers greater than 2.

\subsubsection{Scan Reduction}

\hspace*{\parindent}
David Race of TMC has pointed out that the CM-2 has primitive operations
for recursive sums or products along successive segments of a vector.
Algebraically, the first order recursion (5) can be manipulated into
the form of a running sum by defining
\begin{equation}
\matrix{
{x_i}' = {{x_i} \bigm/
{ \prod_{k=0}^{i} a_k }
} \cr
{y_i}' = {{y_i} \bigm/
{ \prod_{k=0}^{i} a_k }
} \cr}
\end{equation}
to yield the simple running sum
\begin{equation}
y_i' = x_i' + y_{i-1}'
\end{equation}
which can be implemented with a single add-scan.  In addition, the running
products in the transformation (11) can also be computed with one
multiply-scan.

Internally, the scan operations are implemented on a binary tree and
hence perform in logarithmic time, comparable to cyclic reduction
or product expansion.  Its advantage is that it has been subject
to several rounds of optimization at TMC and so runs significantly faster
than, say, $\log nx$ individual communication-plus-add operations.
The disadvantage of scan implementation for T15mig is that the product
of the $a_i$'s approaches zero geometrically and so will quickly lead to
overflow, or at least severe round-off and scaling problems.

\subsubsection{Hybrid Methods}

\hspace*{\parindent}
The methods mentioned above can be combined in various ways
to improve performance.  We have already noted how cyclic
reduction combines with partitioning.

David Race has worked out a hybrid method that combines
scan reduction, partitioning, and series expansion.  The partition
size is chosen to match the number of terms in a truncated series
expansion.  The truncated series is then used to compute the value
of $y_i$ at one end of each partition.  (This is implemented with
a vector multiply and a segmented add-scan.)  The transformation
(11) is then applied to each separate partition (but in parallel.)
Then (12) is applied to the partitions using a segmented add-scan.
Finally the transformation (11) is undone (by vector multiply) and
processing is complete.  By this construction, scaling is now
automatically limited to
the reciprocal of the tolerance used to truncate the series expansion,
e.g. $\epsilon = .001$, and roundoff is controlled.

A hybrid of the product and series expansion also can be useful to
reduce long-distance communication in the product expansion method.
Here, for example, we may employ the identity
\begin{equation}
({\bf I } + {\bf A}^{2^k}) =
(({\bf I } + {\bf A}^{2^{k-1}})^2 - 2 {\bf A}^{2^{k-1}} )
\end{equation}
to reduce long-distance communication at a cost of greater arithmetic.
Since long-distance communication is twice as slow as local
communication which is in turn twice as slow as arithmetic operations,
this can lead to some savings.  Additionally, the CM-2 provides
instructions that combine arithmetic with local communication
resulting in even greater savings.

\subsubsection{Parallel LDU}

\hspace*{\parindent}
For the various expansion methods based on bidiagonal factoring
of the original tridiagonal equations, we would like to compute
the factors with a highly parallel algorithm.  For T15mig, this
is not a critical necessity, as the Cray algorithm allows for
a parallelism of at least 100 in the process.
Since the average number of reuses of the factors is around 100,
a 10\% efficiency in the factoring is an overall penalty of perhaps
10 percent in run time.  However, the technique of cyclic reduction
does permit us to inject parallelism into the factoring step.

The key step in the factoring is to compute a recursion, following
Claerbout \cite{jfc:fgdp}, that takes the form
\begin{equation}
E_k = { {-A_k} \over {B_k + C_k E_{k-1} } }
\end{equation}
which is a special case of the more general continued-fraction
recursion
\begin{equation}
\epsilon_k = { {\alpha_k + \beta_k \epsilon_{k-1} }
\over {\gamma_k + \delta_k \epsilon_{k-1} } }
\hbox{\ \ \ .}
\end{equation}
As with equation (6) above, one can express $\epsilon_k$ in terms
of $\epsilon_{k-2}$ in the same form as above with
\begin{equation}
\matrix{
\alpha_k' & = & \alpha_k \gamma_{k-1} + \beta_k \alpha_{k-1} \cr
\beta_k' & = & \alpha_k \delta_{k-1} + \beta_k \beta_{k-1} \cr
\gamma_k' & = & \gamma_k \gamma_{k-1} + \delta_k \alpha_{k-1} \cr
\delta_k' & = & \gamma_k \delta_{k-1} + \delta_k \beta_{k-1} \cr
}
\end{equation}
as the new coefficients of the linear fraction.  Proceeding recursively,
we can eventually reduce the recursion down to a suitably small series
to be expanded ``by hand'' and back substituted to complete the
calculation.  Alternatively we can use recursive doubling to
evaluate the recursion.  Unlike the apply step in the series expansion,
we do not get a geometric decrease in the coefficient of
$\epsilon_{k-2^j}$ and so cyclic reduction does not have an
early-truncation advantage.

\subsection{Data Movement}

\hspace*{\parindent}
We now return to the question of data movement.  From the previous
paragraph, we find that approximately 10 grids are involved if we
consume the input by shifting the computational mesh backwards.
This is about the same as the number of sheets in the input grid
in a typical run, suggesting that feeding the input data into
the computational grid would not be much faster or slower.
Fortunately, we need only
feed data into the grid from one input sheet at a time.  We also
will need to shift a corresponding sheet of the output grid in order
to store answers.  This means that we need to shift at most 2
sheets of the input/output grid at each computational step.  This is
more efficient than shifting 8 computational grids and so the computational
grid should stay fixed and the input/output grid should move.
One question remains.  How will the output data be arranged in the
equivalenced input/output grid after all this shifting?

The answer to this is given in Appendix A.  There it is shown that
the loop sequence
\begin{displaymath}
\begin{array}{l}
\hbox{Insert input into $P_1$} \\
\hbox{Extrapolate $P_4$} \\
\hbox{Copy $P_4$ to $P_3$} \\
\hbox{Copy $P_1$ to $P_2$} \\
\hbox{Shift $P_4$ to $P_1$} \\
\hbox{Extract output from $P_2$}
\end{array}
\end{displaymath}
will leave the output in proper order, making for easy display or
host retrieval.  The only concern is that the input data on the
last sheet needs to be positioned at the end, rather than the beginning
of the last sheet.  The output data then appears at the beginning
of the last sheet --- exactly where we want it.  The input positioning
can be done as the data is loaded into the CM-2 by specifying a
nonzero {\it start-vector\/} in the {\bf cm\_write\_to\_news\_array}
calling sequence.

\subsection{I/O and Display}

\hspace*{\parindent}
Data must be loaded into the CM-2 and the migrated image must be
retrieved and displayed.  On our CM-2, the data is staged to
disk on the front-end host Vax 6310 before migration is invoked.
The migration velocity model is also placed on front-end disk.
When the program is invoked, the input data is retrieved from
disk and placed into the CM-2 using the above-mentioned
{\bf cm\_write\_to\_news\_array} instruction which transfers the data
at about half a megabyte a second across the CM interface.
The velocity model is repeatedly retrieved in an inner remigration
loop, the allowing the user to set up and modify migration velocities
based upon inspection of the output of the previous migration iteration.
The migrated section is both displayed on the CM-2 internal frame
buffer and placed on output using the {\bf cm\_read\_from\_news\_array}
instruction.  The CM-2 graphics library allows runtime selection
of display in an X-window on a workstation instead of the frame buffer.

\begin{thebibliography}{99}
\bibitem{b&g:pois} B.L. Buzbee, G.H. Golub and C.W. Neilson, 1970,
{\em On direct methods for solving Poisson's equation,}
SIAM J. Num. Anal., v 7, pp. 627-656.
\bibitem{jfc:fgdp} J.F. Claerbout, 1976,
{\em Fundamentals of Geophysical Data Processing,}
McGraw-Hill, New York.
\bibitem{slj:ens} S.L. Johnsson, 1987,
{\em Solving tridiagonal systems on ensemble architectures,}
SIAM J. Sci. Stat. Comput., v 8, n 3, pp. 354-392.
\bibitem{hss:recdub} H.S. Stone, 1973,
{\em An efficient parallel algorithm for the solution of a tridiagonal
linear system of equations,}
J. Assoc. Comput. Mach., v 20, pp. 27-38.
\bibitem{slj:cmsl} A.E. Trefethen, K.K. Mathur and S.L. Johnsson, 1989,
{\em Connection Machine scientific library,}
Fourth SIAM Conf. Parallel Processing for Scientific Computing,
Dec 11-13, Final Program, p. A-46.
\bibitem{vdv:vlr} H.A. van derVorst and K. Dekker, 1989,
{\em Vectorization of linear recurrence relations,}
SIAM J. Sci. Stat. Comput., v 10, n 1, pp. 27-35.
\bibitem{hhw:ptri} H.H. Wang, 1987,
{\em A parallel method for tridiagonal equations,}
ACM Trans. Math. Software., v 8, n 3, pp 170-183.
\end{thebibliography}
\end{document}
