As an example of Fourier domain algorithms I discuss a prestack phase-shift migration (Gazdag, 1978). In a phase-shift migration the prestack data are first transformed into the frequency-wavenumber domain by a multidimensional FFT. Then, the transformed data are downward continued by multiplication with a complex exponential. Finally, the downward continued data are imaged at zero time and zero offset. The imaging at zero time and zero offset can be simply computed by summing the downward continued data for all frequencies, and all offset wavenumbers, without computing the inverse Fourier transform for all times and offsets.

The downward continuation is naturally parallel because the different
wavenumber components of the data can be downward continued independently;
they interact only when
they are summed together for imaging the result at zero offset.
This summation is easily coded in Fortran 90 by the intrinsic function
** SUM**; the compiler takes care of issuing the correct instructions
for the necessary interprocessor communications.
Global communications, as the sum over one or more axes of an array,
are usually very fast because they can use the whole bandwidth
of the communication network linking the processors.
Consequently, the time spent in communications is only
a small percentage of the total run time (about 5-8 % in this particular
case).

Not all the wavenumber components need to be downward continued,
because some of them correspond to evanescent waves.
The zeroing of the evanescent waves cannot be expressed with the usual
serial ** IF** statement, because all the wavenumbers are
processed in parallel.
Therefore the parallel conditional statement

The kernel of phase-shift prestack migration is expressed in Fortran 90 by the following code fragment:

C C loop over migrated time C DO ITAU = 1,NTAU-1 C C compute Vel*Vel*Ks*Ks and Vel*Vel*Kg*Kg C VELSQ = VELTAU(ITAU)*VELTAU(ITAU) KSVELSQ = VELSQ*KS2SQ KGVELSQ = VELSQ*KG2SQTAUSLICE = 0.

C C loop over frequencies C DO IW = 1,NW CMOMEGA = OMEGA(IW)

C C downward propagate the non-evanescent waves C WHERE((CMOMEGA .GT. KSVELSQ).AND.(CMOMEGA .GT. KGVELSQ))

KTAU = -DTAUD2*(SQRT(CMOMEGA-KSVELSQ)+SQRT(CMOMEGA-KGVELSQ)) DATA(IW,:,:) = CMPLX(COS(KTAU),SIN(KTAU))*DATA(IW,:,:)

C C zero the evanescent waves C ELSEWHERE DATA(IW,:,:) = 0. END WHERE

C C image at zero time C TAUSLICE = TAUSLICE+DATA(IW,:,:)

END DO

C C image at zero offset C MIGRES(ITAU+1,:) = SUM(TAUSLICE,DIM=1)

END DO

Notice that in a phase-shift algorithm the different
frequencies components could have also been processed in parallel,
but, by serially looping along the frequency axis, I avoid
the redundant evaluation, and storage, of some working
arrays (** KSVELSQ** and

Prestack Stolt migration (Clayton and Stolt, 1981) can be implemented similarly to phase shift; also for Stolt migration the wavenumber components can be processed in parallel. However, in an efficient implementation of Stolt migration, the frequency axis must be local to each processor for avoiding interprocessor communications during the Stolt transformation of the frequency axis. If the frequency axis is kept local, the Stolt transformation can be simply implemented with the usual interpolation algorithms used on serial computers.

12/18/1997