%
% Converted from RTF format using rtf2TeX
% Comments and bugs to Robert Lupton (rhl@astro.princeton.edu)
%
\font\rm = cmr10  % New York
\headline={\tenrm \kern 36pt \number\pageno\hfil }
\footline={\tenrm }
\input TeX_defs.tex
%
\hsize=7in
\vsize=9in
\hoffset=-0.25in
\parindent=0pt
%
\font\rm = cmr10  % Times
\leftskip=0pt plus 1fill \rightskip=0pt plus 1fill
Parallel Finite-Difference Migration


{\pointsize{10} Stewart A. Levin}


Mobil Exploration and Producing Technical Center

P.O. Box 650232, Dallas, Texas 75265-0232


\pointsize{0} {\it ABSTRACT\/}


\leftskip=0pt \rightskip=0pt
     Finite-difference migration in the natural {\it x---t \/} 
coordinates, e.g. the 15-degree implicit wave-equation method, has 
been in common use for over a decade of geophysical data processing.  
As such, it predated the invention of large vector computers, such as 
the Cray-1S introduced in 1977, and so was organized and optimized to 
achieve the sensible objectives of reducing the number of computations 
and I/O operations as much as possible.  Outboard array processors, 
employed to accelerate computations, changed the relative weights of 
com\-putation and I/O but still fit comfortably into existing program 
designs.  Since the mid 80's the author$^{ 1,2}$ has revisited the 
implementation of 15-degree implicit finite-difference migration on 
vector and parallel computers now in common use and found that 
impressive performance gains can be achieved by completely 
reorganizing the computations.  This report summarizes the highlights 
of these efforts and gives benchmark results for a variety of vector 
and parallel platforms. 

\leftskip=0pt \rightskip=-72pt

\leftskip=0pt plus 1fill \rightskip=0pt plus 1fill
{\it 1. INTRODUCTION\/}

\leftskip=0pt \rightskip=-72pt

\leftskip=0pt \rightskip=0pt
     15-degree finite-difference migration consists primarily of 
repeatedly solving the recursive system


\leftskip=0pt \rightskip=-229.3pt

\leftskip=0pt \rightskip=0pt
$$(1+(\beta - \alpha) \delta_{xx} ) [ P( t + \Delta t , \tau + \Delta \tau , z ) + 
P(t,\tau,z)] = ( 1 + ( \beta + \alpha ) \delta_{xx} ) [
P (t + \Delta t , \tau , z ) + P ( t , \tau + \Delta \tau , z ) ]
\eqno(1)$$

\leftskip=0pt \rightskip=-414pt

\leftskip=0pt \rightskip=-297pt

\leftskip=0pt \rightskip=-229.3pt

\leftskip=0pt \rightskip=0pt

where $\alpha$  is given by the formula



\leftskip=0pt \rightskip=0pt
$$\alpha = { v^2 \Delta t  \Delta \tau } \over { 8 ( \Delta x )^2 }
\eqno(2)$$

\leftskip=0pt \rightskip=-287.9pt

\leftskip=0pt \rightskip=0pt

\leftskip=0pt \rightskip=0pt
where {\it v\/}  is wave propagation velocity and $
\beta$  is typically chosen to be around one-sixth.  
Initial condition $P(t, \tau = 0 , z)$ and boundary condition
$P( t = t_{\max} ,  \tau, z ) \equiv 0 $ are specified and the output image is the 
wavefield $P( t = \tau , \tau , z)$.


{\pointsize{10}      There are two basic tools one can use to 
implement 15-degree implicit finite-difference on a vector or parallel 
platform.  The first is a concurrent organization, that derives from 
now-standard wavefront methods.  This exposes a level of parallelism 
on the order of 100 hidden in the conventional computations.  The 
second tool is some variant of cyclic reduction to introduce 
additional parallelism into bi- or tridiagonal solutions required at 
each stage of the conventional or wavefront differencing kernel.  The 
latter is not necessarily an advantage to use as, asymptotically,  it 
uses twice as many arithmetic operations as the best serial algorithm. 
  In each of the following sections, I describe how I employed these 
techniques to implement 15-degree implicit finite-difference on 
several different platforms.}

\leftskip=0pt \rightskip=-72pt
\pointsize{0} 

\leftskip=0pt plus 1fill \rightskip=0pt plus 1fill
{\it 2. FPS 120 ARRAY PROCESSOR\/}

\leftskip=0pt \rightskip=-72pt

\leftskip=0pt \rightskip=0pt
     The FPS 120B was one of the most successful of the attached array 
processors.  Its peak speed was 12 Mflop/sec using independent, 
pipelined multiplier (3 stage) and adder  (2 stage) units.  For almost 
a decade it was the geophysical workhorse at the Stanford Exploration 
Project to which the author belonged.  During this period, the 
conventional 15-degree implicit kernel was microcoded for speed and 
achieved a performance of about 2 Mflops/sec on the standard memory 
model  (1~memory access every 2 cycles).  The limiting factor was 
memory access  and neither the wavefront method nor cyclic reduction 
were helpful.  For a fast-memory model , such as the one formerly at 
Mobil, the memory access rate was doubled and the microcoded serial 
algorithm was not optimal.  Here the performance automatically doubled 
in all parts of the code except the key recurrence


\leftskip=0pt \rightskip=0pt
$$ y_i = x_i + a_i y _{i-1} \eqno(3)$$

\leftskip=0pt \rightskip=-292.5pt

\leftskip=0pt \rightskip=0pt


where the 5 cycle mult-add can't be pipelined due to the recurrence.  
Either the wavefront method or cyclic reduction can be applied to 
speed this up; it is instructive to compare their utility.  In the 
wavefront approach, two independent recurrences are interleaved in the 
arithmetic pipelines.  This  nearly doubles the performance, implying 
2 results every 6 cycles instead of 1 every 5.  Cyclic reduction also 
allows us to interleave two independent recurrences derived from 
equation (3):

\leftskip=0pt \rightskip=0.5pt

\leftskip=0pt \rightskip=0pt
$$y_{2i} = ( x_{2i} + a_{2i} x_{2i-1}) + b_{2i} y_{2i-2}
\eqno(4a)$$

\leftskip=0pt \rightskip=0.5pt

\leftskip=0pt \rightskip=-319pt

\leftskip=0pt \rightskip=-300.7pt

\leftskip=0pt \rightskip=0pt
$$y_{2i+1} = (x_{2i+1} + a_{2i+1} x_{2i} + b_{2i+1} y_{2i-1}
\eqno(4b)$$

\leftskip=0pt \rightskip=-328pt

\leftskip=0pt \rightskip=-309.7pt
where


\leftskip=0pt \rightskip=0pt
$$b_i = a_i a_{i-1}\eqno(5)$$

\leftskip=0pt \rightskip=-279pt

\leftskip=0pt \rightskip=0pt

can be precomputed at negligible cost.  This formulation leads to 2 
results every 8 cycles, not as good as using the wavefront approach.  
A  superior arithmetic improvement can be obtained by dropping the 
second recurrence (4b) and using the original formula (3) to 
back-substitute for $y_{2i-1}$:


\leftskip=0pt \rightskip=0pt
$$y_{2i-1} = x_{2i-1} + a_{2i-1} y_{2i-2} 
\eqno(4b^\prime)$$

\leftskip=0pt \rightskip=-301pt
(4b)

\leftskip=0pt \rightskip=0pt

\leftskip=0pt \rightskip=0pt
which permits 2 results with 7 memory operations  nearly the same rate 
as the wavefront method.  There is no benefit here in taking the 
recurrence any deeper as that only increases the memory load per 
result.  There is also no benefit to increasing the number of 
differencing stars solved simultaneously in the wavefront 
organization.  The final performance figures were 1/6th of peak for 
standard memory and 1/3rd of peak for fast memory.


\leftskip=0pt plus 1fill \rightskip=0pt plus 1fill
{\it 3. CONVEX\/}


\leftskip=0pt \rightskip=0pt
     The Convex C-1 was the first vector computer on which I 
implemented the wavefront method.  The original program from Stanford 
University used an optimized, scalar code for solving the key 
recurrence (3).  Loops were unrolled in order to coax the compiler 
into reusing the output of the recurrence before storing it to memory. 
  In the vectorized versions, both the wavefront and the cyclic 
reduction techniques were tried.  Cyclic reduction improved 
performance by about 3:1.  The wavefront approach, which I implemented 
by solving hundreds of differencing stars at a time, proved even 
faster, cutting run times by a factor of 5 over the serial version 
when the key recurrence was assembly coded for extra speed.  This 
worked out to be 1/3rd of the peak speed of the machine, consistent 
with the memory-to-arithmetic load of the algorithm on this 
single-pipe memory architecture.


{\pointsize{10}      (After first publishing the code for the fast 
Convex version in a 1986 Stanford Exploration Project sponsor's report, 
I was told by Bee Bednar of Amerada Hess that he ran it on an IBM 
3090.  After reconfiguring it to increase the available memory, the 
standalone test sustained over 90 percent of the peak speed of the 
machine!  It took a few months, however, to rewrite it to save memory 
by working in strips of 64 equations at a time, the machine's vector 
length, and using external disk storage.)}


\leftskip=0pt plus 1fill \rightskip=0pt plus 1fill
\pointsize{0} {\it 4. CRAY X-MP\/}

\leftskip=0pt \rightskip=0pt

     After implementing a vector version on the Convex, the next 
challenge was to tackle a vector parallel version on a multi-headed 
Cray XMP.   The first step was to implement the pure vector version on 
a single-headed X-MP.  This proved quite easy and the pure Fortran 
version worked quite well.  Indeed, the Fortran worked so well that it 
consistently beat every attempt at handcoded assembly language I 
tried.  My test 1024x256 problem ran in 16 seconds, compared to 180 
seconds on the Convex C-1.  Of the factor of approximately 10 
performance ratio a factor of 5 is due to cycle time and the remaining 
factor of 2 arises from the extra memory bandwidth of the Cray which 
allows up to two memory reads and one memory write on each cycle.  The 
ability to read twice as many memory operands per cycle lifted the 
Cray performance to about 2/3rds of the Cray peak vector speed on this 
application.  In principle, this application is no longer memory 
limited on this machine and ultimate performance approaches the peak 
of 11/14ths or 78 percent purely on the basis on flop count of 
equation (1).  Memory and arithmetic pipeline overhead do, however, 
reduce the measured performance.


    The goal of vector-parallel performance can be achieved by 
augmenting the vector algorithm in two ways.  The first is to 
microtask or autotask the inner loops.  This approach  expands the 
wavefront parallelism by simply increasing the implied machine vector 
length to 64N where N is the number of processors assigned to the 
task.  This indeed reduces the runtime of the test problem to just 
under 5 seconds, pushing migration into the interactive realm.  (For 
comparison, the FPS 120B took about 20 minutes to perform the same 
calculations.)  To reduce run times even more, I modified the test 
problem to perform approximately 8 times fewer calculations, albeit 
with a corresponding 8-fold reduction in parallelism.  In this case, 
the wavefront parallelism was 128, or twice the machine vector length, 
rather than the original 1024.  Therefore to apply more than 2 
processors at a time to the problem, it is necessary to use cyclic 
reduction to increase parallelism further.  Indeed, I obtained time on 
a 4 processor X-MP at the San Diego supercomputer center for this 
purpose.  The idea was to use two stages of the cyclic reduction 
(4a,b) to produce four independent recursions to be solved 
independently.  This was not accomplished for the simple reason that 
the operating system at that time did not have full multitasking 
support.  In hindsight, though, the idea was flawed  each level of 
parallel reduction doubles the number of mult-adds needed to solve the 
key recurrence (3).  Thus doubling the number of processors will 
produce results no faster.  Instead a serial reduction, i.e. of the 
form ($4b'$) rather than ($4b$), is required to avoid doubling the 
arithmetic load.  In this case the mult-add load is increased by 50 
percent while the number of processors is doubled.  Theoretically, one 
can hope to reduce the recurrence cost by 1/3rd.  Unfortunately, there 
is no good way of transferring data from the vector registers of one 
processor to that of another except by storing it into memory.  This 
limits the recurrence speed to at best 8 results every 9 cycles, worse 
than not using cyclic reduction at all.


\leftskip=0pt plus 1fill \rightskip=0pt plus 1fill
{\it 5. AMETEK 64-NODE HYPERCUBE\/}


\leftskip=0pt \rightskip=0pt
     The first purely parallel machine on which I implemented a 
parallel version of 15-degree migration was the Ametek 64-node 
hypercube multiprocessor.  The programming paradigm was C (or Fortran) 
on each node with an explicit message-passing library for 
interprocessor communication.  A host computer controlled cube 
initialization and external com\-munication.  The program was written 
over a period of two weeks based on library reference manuals for the 
machine.  It was then brought down to CalTech on a tape and loaded 
onto a newly-acquired machine owned by Rob Clayton.  We loaded the 
operating system at lunch time and by 6 PM the test program was 
running.  Run times for the 1024 by 256 problem were about 20 minutes, 
i.e. almost the same speed as the FPS 120B runs.  This worked out to 
an efficiency of 1/15th the rated machine peak speed.  Some of this 
(about a factor of two) was  elegant but inefficient programming on my 
part.  The rest was poor compiler optimization of floating point 
arithmetic.  This was my first, but not my last, observation of low 
efficiency on parallel platforms.  One should keep this in mind when 
comparing price-performance numbers of parallel and vector computer 
architectures.


\leftskip=0pt plus 1fill \rightskip=0pt plus 1fill
{\it 6. SAXPY 1-M\/}

\leftskip=0pt \rightskip=0pt

     The Saxpy 1-M was a 1 Gflop/sec design consisting of 32 ganged 32 
Mflop/sec pipelined arithmetic processors.  These processors were 
interconnected by both a global bus and individual systolic 
connections to their right-hand neighbor and were supported by a large 
gigabyte-sized memory and a MicroVAX control processor.  
Unfortunately, the program was never run on an actual machine before 
the company folded, but a thorough analysis and simulator results 
permitted a fairly good estimate of its performance$^{ 3}$. 


     In this version of the algorithm,  the systolic connections of 
the processors fit quite elegantly into the computation.  Each 
processor handled a separate differencing star and then all performed 
a right shift as if in a bucket-brigade as computations advanced to 
the next set of time levels.  In addition, cyclic reduction was 
employed to improve the vector performance of the in-processor 
recurrence.  Here was a case in which  cyclic reduction did contribute 
significantly in keeping the 9 stage pipeline nearly full in exchange 
for an overall doubling of arithmetic in the bidiagonal solve.  The 
final performance figure was estimated to be 375 Mflop/sec, i.e. 
3/8ths of the machine peak.  (The actual Mflop/sec was about 500, but 
this included the doubled arithmetic load of cyclic reduction). 
Allowing for miscellaneous overheads, this drops to somewhere around 
1/3rd of peak performance, quite comparable to the Convex pure-vector 
efficiency. 


\leftskip=0pt plus 1fill \rightskip=0pt plus 1fill
{\it 7. NCUBE 64-NODE HYPERCUBE\/}

\leftskip=0pt \rightskip=0pt

     After the Ametek experience, I had an opportunity to repair 
deficiencies in my message-passing hypercube migration code and rerun 
it  on an Ncube-1.  This model was a 64 node machine with a 7 MHz 
clock and an aggregate peak floating point rate of about 22 Mflop/sec. 
 Unlike the Ametek, it took almost 3 months to get around hardware and 
software problems and achieve a working code.  The most interesting 
lesson was that while a linear geometry is easily obtained by breaking 
a ring, this act can break your program.   The ``gotcha'' is that the 
processor at the end of the ring doesn't have to worry about sending or 
receiving from one of its neighbors and can therefore flood its other 
neighbor's message buffers during repeated nearest-neighbor shifts.


     The initial results for the test problem indicated a sustained 
speed of 8 percent of the peak speed of the machine.  Parallel speedup 
was quite good, however, about 55 out of a possible 64.  Detailed 
instrumentation brought to light a curious anomaly: the program 
appeared to slow down a fraction of the way through the calculations.  
This was found to be a quirk of the Ncube processor design which 
detected and short-circuited addition with or multiplication by zero.  
When this was accounted for, the sustained speed dropped to 4 percent 
of peak performance.   As with the Ametek, this performance loss was 
almost entirely due to lack of floating point optimization of the node 
compiler I used.


     One additional variation I tested on the Ncube was to devise a 
multicircuit through the hypercube network that employed all the 
communication paths in the array rather than only the third that the 
normal Gray-code path used.  Theoretically this makes three times the 
bandwidth  available to my ring-structured algorithm and can be used 
in a number of ways to improve performance.  In my case the goal was 
to get all processors doing useful computations (as opposed to working 
on dummy zero vectors) as quickly as possible.  At the same time, I 
wanted to keep the communication-to-compute ratio low.  This 
multicircuited mapping promised both.


     To implement this mapping I took advantage of the Ncube's ability 
to send and receive messages on more than one communication channel at 
a time.  By breaking each outgoing message into sufficiently large 
packets, 256 bytes or greater,  I was able to successfully hide the 
additional message startup overhead and reduce communication time by 
half.  Unfortunately, this was not a significant savings here because 
the slow arithmetic already accounted for 95 percent of the run time.  
Based upon hardware timings, I estimate that handcoded assembler could 
speed up the arithmetic by nearly a factor of ten.  In that event, the 
arithmetic drops to 75 percent of the run time and halving the 
communication time would be significant.


\leftskip=0pt plus 1fill \rightskip=0pt plus 1fill
{\it 8. THINKING MACHINES CM-2\/} 

\leftskip=0pt \rightskip=0pt

     The next opportunity for parallel finite-difference programming 
was on an 16K CM-2 at Mobil Research and Development.    The big 
challenge here was to boost parallelism in the code by another factor 
or 100 over the Cray or Ncube versions in order to use the machine 
relatively efficiently. To accomplish this, a combination of both 
wavefront parallelism and cyclic reduction was employed.  In this 
implementation it was found that the LDU tridiagonal factorization 
used in all the versions of the code was not advantageous.  Instead, 
the conventional {\it sericr\/}$^{ 4}$ version of cyclic reduction was 
applied directly to the tridiagonal equations once, and the 
internally-generated coefficients stored for later reapplication.  
This prefactoring still saves about 50 percent of the arithmetic, just 
as LDU does for the conventional marching algorithm.  David Race of 
Thinking Machines Corporation managed to nicely parallelize the 
coefficient computation as well.


     The resulting implicit finite-difference program was implemented 
with (low level) PARIS calls.  The algorithm is expressible in 
CM-Fortran, but at that time, the compiler was not yet capable of 
correctly handling a few key constructs.   The standard 1024x256 test 
problem ran in approximately 50 seconds, about 2 $^{ 1}$/$_{ 2}$ times 
slower than the Cray.  This run time increased only a small percentage 
when using only an 8K half of the machine.  This reflects the limited 
parallelism available for that problem size even with cyclic reduction 
techniques.  By this measure, the CM-2 delivered approximately 3 
percent of peak efficiency on this application, comparable to the 
measured performance on other parallel machines.


     The most interesting discovery made in the course of the CM-2 
implementation, was that cyclic reduction in this setting can be 
significantly optimized.  The reason for this is that the original 
tridiagonal system is diagonally dominant.  This causes the 
progressive coefficients to become vanishingly small as the reduction 
proceeds.  By monitoring these coefficients, one can terminate the 
reduction early and save significantly. Not only is the arithmetic 
load reduced, but also the parallelism is kept high.  By 
experimentation with the test case, I found that the cutoff point is 
3-4 reduction stages.


\leftskip=0pt plus 1fill \rightskip=0pt plus 1fill
{\it 9. SILICON GRAPHICS POWERSTATION 8\/}

\leftskip=0pt \rightskip=0pt

     A recent trend in the computer arena is called the attack of the 
killer micros.  These are typically workstations based on RISC 
technology that are capable of performing calculations at rates 
previously only achieved on large mainframes.  Not long ago, Silicon 
Graphics produced a line of workstations designed to beef up the 
floating point performance of their high end graphics workstations.  
This PowerStation series used multiple processors accessing a shared 
memory to tackle this goal.  I used 15-degree implicit 
finite-difference to benchmark their 8-processor model.


     In the initial work, I reduced the test problem size to 256 by 
512 and used the vectorized version of the algorithm from the Convex.  
For a single processor, this code took 5 $^{ 1}$/$_{ 2}$ minutes to 
run.  When automatic parallelization was turned on, all 8 processors 
were fully active, but wall clock time decreased to only 4 minutes.  
The ineffective parallelization was determined to be due to data 
layout.   The data access pattern was guaranteeing that every 
processor was loading all pages of the data into its cache.  This 
8-fold increase in memory traffic greatly degraded performance.


     To combat the cache access problem, all arrays in the original 
code were transposed.  Do loops were correspondingly interchanged to 
match.  With manual directives to enforce parallelization on the outer 
loops, run times decreased to 2 minutes. 


     Finally, I ran my original test problem with 4 times the number 
of computations.  Arithmetically, therefore, one would expect to see a 
run time 4 times longer, i.e. 8 minutes using all 8 processors.  The 
actual run time was 19 $^{ 1}$/$_{ 2}$ minutes, a factor of 2 $^{ 
1}$/$_{ 2}$ slower than anticipated.  Here again cache memory was the 
culprit.  With the larger data dimensions, each processor's piece of 
the job was simply too big to fit into the approximately 300 kilobytes 
of per-processor cache.  Suddenly, all 8 processors were competing for 
data in main memory, hence the degradation.   For this situation, it 
is possible to rearrange calculations to fit within cache limitations 
by making multiple passes over the data.  This trick was learned from 
my nCube experience where each processor only had 128K of memory.  
With this reorganization it is possible to do the test case migration 
in slightly over 8 minutes on the 8 processor Silicon Graphics.  
Relative to the Convex C210, this is about a factor of 5 slower.

\leftskip=0pt plus 1fill \rightskip=0pt plus 1fill
 

{\it 10. SUMMARY\/}

\leftskip=0pt \rightskip=0pt

     Implicit 15-degree finite-difference has proven to be a very nice 
testbed for experimenting with and benchmarking vector and parallel 
systems.  Each system needed a different version of the algorithm to 
achieve good performance.  A number of useful lessons did arise from 
these experiments.  Wavefront parallelism worked well on all platforms 
because it revealed the parallelism inherent in the original, 
arithmetically optimized algorithm.  Cyclic reduction was helpful in 
increasing parallelism when the wavefront parallelism was not enough 
to keep all processors busy.   No actual implementation of the 
algorithm on a distributed-memory parallel machine, however, ever 
achieved even 10 percent efficiency.  Finally, distributed memory 
organization proved helpful in designing good code for a cache-based machine.


\leftskip=0pt plus 1fill \rightskip=0pt plus 1fill
{\it 11. REFERENCES\/}


\parindent=-9pt 
\leftskip=9pt \rightskip=0pt
1. S.A. Levin, Parallel Finite-Difference Migration, {\it 56th Annual 
Mtg. Extended Abstracts\/}, S7.7, p. 440-7, Society of Exploration 
Geophysicists, Tulsa, 1986.

2. S.A. Levin, A novel hypercube mapping with application to parallel 
seismic migration, Fifth SIAM Conference on Parallel Processing for 
Scientific Computing, Chicago, 1991.

{\pointsize{10} 3. S.A. Levin and T. Parks, Parallel Finite-Difference 
Migration, Second SIAM Conference on Parallel Processing, Norfolk, 1985.}

\pointsize{0} 4. R.W. Hockney and C.R. Jesshope, {\it Parallel 
Computers\/}, p. 286-9, Adam Hilger Ltd., 1981.


\bye
