NASA Contractor Report 189724

ICASE Report No. 92-56

/33

2/sICASE

A GENERAL PURPOSE SUBROUTINE FOR FAST

FOURIER TRANSFORM ON A DISTRIBUTED

MEMORY PARALLEL MACHINE

A. Dubey

M. Zubair

C. E. Grosch

(NASA-CR-189726) A GENERAL PURPOSE

SUBROUTINE FOR FAST FOURIER

TRANSFORM ON A DISTRIBUTED MEMORY

PARALLEL MACHINE (ICASE) 15 p

G3/61

N93-16681

Unclas

0133988

NASA Contract Nos. NAS1-18605 and NAS1-19480

November 1992

Institute for Computer Applications in Science and Engineering

NASA Langley Research Center ......

Hampton, Virginia 23681-0001

Operated by the Universities Space Research Association

N/t A

National Aeronautics and

Space Administration .....................

Langley Research Center

Hamplon, Virginia 23665-5225

A GENERAL PURPOSE SUBROUTINE FOR

FAST FOURIER TRANSFORM

ON A DISTRIBUTED MEMORY PARALLEL MACHINE

A. Dubey, M. Zubair and C.E. Grosch

Department of Computer Science

Old Dominion University

Norfolk, VA

ABSTRACT

One issue which is central in developing a general purpose FFT subroutine on a distributed memory parallel machine is the data distribution. It is possible that different users

would like to use the FFT routine with different data distributions. Thus there is a need to

design FFT schemes on distributed memory parallel machines which can support a variety

of data distributions. In this paper we present an FFT implementation on a distributed

memory parallel machine which works for a number of data distributions commonly encountered in scientific applications. We have also addressed the problem of rearranging the data

after computing the FFT. We have evaluated the performance of our implementation on a

distributed memory parallel machine Intel iPSC/860.

1Research for the last two authors was supported in part by the National Aeronautics and Space Administration under NASA Contract Nos. NASI-18605 and NAS1-19480 while the authors were in residence

at the Institute for Computer Applications in Science and Engineering (ICASE), NASA Langley Research

(:enter, Hampton, VA 23681-0001.

1 Introduction

Fourier transform plays an important role in mathematical and numerical analysis. Some

of the applications which use the Fourier transform are: digital filtering, auto and cross

correlation, solution of partial differential equations etc. The fast Fourier transform algorithm (FFT) computes the transform of an n-component sequence of complex numbers in

O(n log n) time. The implementation of fast Fourier transform on parallel machines has been

well studied, for example [6, 3, 4, 7, 5, 8, 11]. However, there has been very little effort in

investigating implementations which can provide a basis for a general purpose subroutines

on distributed memory machines.

One issue which is central in developing a general purpose subroutine on a distributed

memory parallel machine is the data distribution [2]. It is possible that different users wish

to use the FFT routine with different data distributions. Typically, users determine their

data distribution based on the over all application requirements, which could be different for

different users. Thus there is a need to design FFT schemes on distributed memory parallel

machines which can support a variety of data distributions.

There are two possible approaches to this problem. The first one is to design an FFT

subroutine for a specific data distribution which gives optimal performance, along with a set

of basic communication subroutines to convert a user data distribution to the specific data

distribution. This approach has the problem of rearranging the user data initially which is

quite costly on distributed memory parallel machines. The second approach is to design an

FFT routine which works well for arbitrary data distributions. We have designed a scheme

which takes the second approach, but only supports a set of common data distributions as

opposed to arbitrary data distributions. A common set of data distributions, referred to

as block scattered distributions, has been identified by Walker and Dongarra [14] as very

useful for distributed memory parallel machines. Block scattered distributions encompass

the two most common data distributions; the linear data distribution and the scattered data

distribution. For a one dimensional data set, a block scattered distribution is specified by

the block size. The data is divided into a set of equal sized blocks. A block j is mapped to

node (j rood p), where p is the number of nodes. For example, two data distributions for a

one dimensional array of i6 data values on a 4 node machine with two different block sizes

are shown in Figure l(a) and l(b).

In this paper we give an FFT implementation on a distributed menaory parallel machine

for block scattered distributions of data with different block sizes. We have also addressed the

problem of rearranging the data after computing FFT on the same machine. The motivation

for rearrangement comes from problems such as solution of partial differential equations using
spectraltechniqueswhich requirethe final data distribution to beidentical to the initial one.

Finally, weevaluatedthe performanceof our implementation on the Intel iPSC/860.

The rest of the paper is organizedasfollows. In the secondsection wediscussthe issues

involvedin designingaparallel FFT algorithm. In the third sectionweproposea parallel FFT

algorithm for variable block size. The fourth section describesthe rearrangementproblem.

In the fifth sectionthere is a discussionof the experimentalresultsand in the sixth section

we give conclusions.

2 Fourier Transform

The DFT, X(k), of an N-point sequence x(r) is defined as,

N-1

X(E) = __, x(r)e -j2=_k/N, 0 <_ k < N- 1,

r=0

where j = v/Z- 1. We assume, for convenience, that N is a power of 2.

(1)

FFT Algorithm

There are two distinct classes of FFT algorithms; namely decimation in time, and decimation

in frequency. We have used a decimation in frequency FFT algorithm and it is described

here briefly. (For details one can refer to [1]). To begin with, the N-point sequence x(r) is

divided into two halves, Xl(r) and x2(r) so that the the transformed sequence can be written

as

(N/2-1)

X(2k) = [x,(,') + x2(r)] N (2)

r=O

(N/2-1)

X(2k+l) = _ [Xl(F)--.T,2(7")]aJrNaJN 2kr, k=O, 1..X/2-1. (3)

r=O

where ornega_u = e -j2_rr/N.

These equations represent two N/2-point DFT's of sequences [xa(r) + x2(r)] and [x,(r) -

x2(r)]J u. The process is then repeatedly applied to the two subsequences. An example

of an eight point decimation in frequency FFT algorithm is shown in Figure 2(a). The

fundamental unit of computation that we use in our algorithm is a butterfly (see Figure

2(b)).

2.1 Parallel Implementation

A typical implementation of FFT on a distributed memory machine results in a sequence

of butterflies at each node interspersed with internode communication. Depending upon the
initial distribution, data for someof tile butterflies is availablelocally and for others, off-node

data are required. There are two approaches for computing butterflies which need off-node

data. The first approach splits a butterfly between two nodes, and in the second approach

a complete butterfly is computed on a node. The parallelism in the later case is achieved

by distributing different butterflies on different nodes. For example, consider a simple case

of computing two butterflies on a two node machine as shown in Figure 3(a). Notice that

both butterflies need off-node data. The two approaches are illustrated in Figure 3(b) and

Figure 3(c) respectively. It is obvious from these figures that the first approach has certain

disadvantages. These are:

High communication volume. The first approach requires twice the inter-node communication volume as compared to the second approach.

Unbalanced computational load . The first approach results in additional computation

on some of the nodes. For example, the multiplication by w in the computation of both the

butterflies of Figure 3(b) is done on node P1, unlike the second approach (See Figure 3(c).)

Extra storage. The first approach requires twice the storage of the second approach.

For these reasons our parallel implementation of the FFT algorithm is based on the second

approach.

3 Algorithm

As stated earlier, the main feature of our FFT scheme is that it works for block scattered

data distribution with variable block sizes. That is, the same algorithm can be used for

different data distributions without any initial rearrangement of the data. The algorithm

consists of three phases: the first and the third phase compute butterflies for which the data

is locally available, and the second phase computes butterflies for which off-node data is

required. As a result, internode communication is required only during the second phase.

Depending upon the block size the work distribution for the first and third phases will differ.

In the extreme cases one of these two phases will not be executed. For a block size of one

we need to execute only the first two phases, while for the block consisting of all the data

on a node only the last two phases are executed. For all other block sizes all three phases

of the algorithm are executed. Given the number of processors, the amount of work in the

second phase remains constant for all block sizes. When all the data on a node forms a single

block, it must be treated as a special case in phase 2 of the algorithm. This is so because

to compute butterflies needing off-node data the block must be divided into two sub-blocks,

which does not happen for other block sizes.

An FFT algorithm with data size N has log2(N ) distinct stages of computation. Each
of these stagescompute N/2 butterflies. In our FFT scheme, as in most other parallel

FFT schemes, all the nodes participate in computing a stage by operating on different data

points. We describe the algorithm by considering an N point FFT on a p-node machhle,

with n = N/p as tile number of data points mapped per node and b as the block size. The

distribution of work in the three phases is as follows

(i) The first log2(n/b ) stages are computed in the first phase. Tile butterflies in these

stages require data available locally from different blocks.

(ii) The next log2(p) stages are computed in the second phase. The butterflies need offnode data, hence internode communication takes place.

(iii) Tile last Io92(b ) stages constitute phase three. Tile butterflies in phase three are again

computed with local data, but from within a block.

The algorithm has a computational kernel dftstep which is common to all the three

phases. The kernel is common to all nodes and computes all the butterflies of a stage

mapped onto a node. It assumes that the co's values have been precomputed and arranged

so that they are available in the right order as needed. A FORTRAN call to the kernel can

be made as follows:

call dftstep(a, w, offset, groups, dist, wincr)

where the arguments are:

a Array of input sequence.

w Array of trigonometric coefficients _o.

offset Tile distance between the two elements of a butterfly.

groups The number of similar sets of butterflies.

dist The distance between two groups.

winer The stride for w.

The kernel dftstep essentially computes n/2 butterflies. How these butterflies are formed

is determined by the three arguments offset, groups and dist. For example consider the

data in p0 in Figure l(a). If a stage requires that the butterflies be formed by combining

x0 with x8 and x, with x9 then a call to dftstep will have offset = 2, groups = 1 and any

value in dist. On the other hand, if a stage requires combination of x0 with x, and x8 with

x9 then the values of these parameters will be offset = 1, 9roups = 2 and dist = 2.

The pseudo code given below describes the FFT algorithm using 'dftstep'. { This code

is executed on each node }

4
begin{phase 1}

offset= n/2 {distance between two points of a butterfly }

groups = 1 {only one subgroup in the first stage }

wincr = 1 { stride for w }

for i= 1 to log2(n/b ) do

dftstep(a,w,offset,groups,offset*2,wincr)

offset = offset/2

groups = groups*2

end for

end {phase 1}

begin {phase 2 }

offset = n/2

groups = offset/b

for i = 1 to log2(p) do

{negh is node with i _h bit differing}

{exchange half the data with negh }

negh = mynode • 2 d-i+1

exchange(a(k),negh) {if bit i = 0 then k = 1, else k = n/2+l }

dftstep(a,w,offset,groups,b,wincr)

if(b = n) then {special case}

negh = rnynode _) 2 d

exchange(a(k),negh)

end if

end for

end {phase 2 }

begin {phase 3 }

for i = 1 to log2(b ) do

groups = groups*2

dftstep (a,w,offset,groups,offset*2,wincr)

offset = offset/2

wincr = 2*wincr

end for

end {phase 3 }

The call to procedure exchange initiates a send first and then posts a receive for incoming

data. This protocol is followed to ensure concurrent communication [10]. The working of
the three phasesis shownin Figure 4 with tile initial data distribution of Figure l(a). In

this problema total of 4 stagesarerequired. The first phaseis computedin stageoneof the

algorithm. The butterflies in this stageareformed by the correspondingelementsof the two

blocks. The secondphaseis computed in stages2 and :3which require exchangesof data.

The fourth stageforms the third phase,which is computedby combiningthe data within a

block.

4 Rearrangement

In general FFT algorithms generate tile resultant sequence in an order different from the

original one. As a result, the data is in the the wrong node. Also all the data to be sent to

one node may not be contiguous. In circuit switch or message passing environments, sending

one data at a time to a node is extremely expensive due to the overheads [9]. Hence it is

desirable that all the data destined for one particular node be collected in one place and

sent together. In addition, at the destination node the data from different nodes may need

to be interspersed. In view of these problems we require the following steps in the process

of rearrangement at each node.

(i) For each data point on a node determine the destination. The destination includes the

node number, the block number and the displacement within the block.

(ii) Collect all the data destined for one node in one set to avoid multiple sends.

(iii) Send the collected data to the appropriate nodes and receive data being sent by other

nodes.

(iv) Place each data item in its destination block with the correct displacement.

The data rearrangement problem in the worst case is equivalent to the complete exchange

problem [13, 12] (see Figure 5(a)). That is each node needs to send data to all the other

nodes in the machine. However, that is not always the case. For example, the FFT shown

in Figure 5(b) requires data exchange between only a subset of nodes. Each node sends data

to and receives data from only two out of four nodes. This situation usually arises when the

number of blocks mapped per node is of the same order as the number of of data points per

block.

5 Experimental results

We evaluated the performance of our implementation on Intel iPSC/860 for different block

sizes and for different data sizes. The Intel iPSC/860 is a distributed memory machine

6
which can have up to 128nodes. Internode communicationis done through a hypercube

interconnectionnetwork. We carried out all of our experimentson 64 nodesof a 128nodes

machine.Tile results for different block sizesaresummarizedin Figure 6. Tile two curvesin

Figure 6 representperformanceper node in Mflops for the completecodeand for tire FFT

sectionalone. Recall that the completecodeconsistsof the FFT computation followed by

the data rearrangement. Theseresults indicate that while tile performanceof FFT section

varies by small amounts, tile overall performanceshowsmorenoticeable differences.Both

curves tend to peak in tile middle. The variation in tile overall performanceis due to tile

data rearrangement.Dependingon the block size, tile data rearrangement may or may not

be a complete exchange problem. Also, observe from Figure 7 that the data rearrangement

time forms a large enough fraction of the total execution time to have influence on the overall

performance.

The reason for the slight variation in the FFT section performance lies in tile relative

distribution of work between the three phases of the algorithm. In Figure 8 we have plotted

tile effect of block size on relative distribution of work in the three phases of the FFT section.

As expected, the fraction of time taken by the first phase is maximum for tile smallest block

size and steadily decreases as the block size increases. The third phase exhibits a reverse

trend. TILe region where both these phases have approximately equal work is also tile region

which shows higher performance in Figure 6. Also notice from Figure 8 that the fraction of

time used in the second phase remains almost constant for all block sizes. The the second

phase takes rnore time than the other two since it also involves internode coImnunication.

The effect of data size on the communication and computation time of the fit section of

the code is shown in Figure 9. To plot this we picked tile best performance for every data

size. With very small data sizes ahnost the entire time is taken up by the communication.

As the data sizes increase, computation starts taking larger fractions of the execution time.

The computation fraction tends to saturate when the data sizes become sufficiently large.

To explain the reason for this behavior we consider the communication characteristics of the

machine

Ill Intel iPSC/860 tile cost of communication is determined by the expression

t¢om_ = 164 + 0.398c_ + 29.9fl

where c_ is the number of bytes in the message and fl is tile distance between two ,Lodes [12].

The first term in the equation is the setup overhead. As c_ increases for fixed fl, the fraction

of total time used in setup decreases. Notice from Figure 9 that the saturation occurs for

c_ = 211. The value of fl is 1. For this data size the contribution from tile overhead term in

the expression is about 5%.
6 Conclusions

In this paper we gave an FFT implementation on a distributed memory parallel machine

which works for a number of data distributions commonly encountered in scientific applications.

We evaluated the performance of our implementation on the Intel iPSC/860. The results

of our experiments indicate that the variation in block sizes has more effect on the performance of the rearrangement section than on the FFT section. For certain block sizes the data

rearrangement cost is significantly lower than for others. In the FFT section, the variation

in performance is not significant enough to make the choice of block size a critical issue. On

a 64 node machine we obtained a peak performance of 203 Mflops, that is 3.17 Mflops per

node. (This figure does not include the initialization costs which are incurred only once for

a given size FFT) If we include the data rearrangement, the performance decreases to 2.78

Mflops per node.

We believe that the data rearrangement can be made more efficient. At present we have

adapted one of the complete exchange schemes. It is not clear whether it is the best option

for the data rearrangement problems arising in FFT or other similar algorithms. Also, the

performance can be further improved (expected to be about 25%) by using an optimized (at

assembly level) version of dftstep.
4

3.5

3

2.5

MF|ops

per 2

node

1.,5

1

0.,5

0

0

I I I I I I I

I ..l I J I [__ !

2 4 6 8 10 12 14

log2(b)

Figure 6. The variation in performance as blocksize is increased.

100

8O

60

Time %

40

2O

I I t I I I I I

!

fit fraction

rearrange fraction

0 , 1 _ , l _ i

/

0 2 4 6 8 10 12

log2(b)

14

Figure 7. The fraction of time taken by FFT and Rearrangement.
100

8O

60

Time%

40

2O

0

0

I l I I I

phase 3

phase 2

1 I I I I

2 4 6 8 10 12 14

log2(b)

Figure 8. The distribution of work in the three phases with different blocksizes.

100

80

60

Time %

40

20

l 4 t I

communication _'sat ulation

0 I " I I I

4 6 8 10 12 14

log2(n)

point

16

Figure 9. The relative contributions of computation and communication in the FFT section

as the datasize is varied.

10
References

[1]

[2]

[3]

[4]

L.R. Rabiner and B. Gold, Theory and Application of Digital Signal Processing,

Prentice-HM1, 1975.

C.L. Lawson, R.J. Hanson, and D.R. Kincaid, and F.T. Krogh, "Basic linear algebra subprogram for FORTRAN usage," ACM Trans. Math. Sotw. 5, 3, 1979,

pp308-323.

R.W. Hockney and C.R. Jesshope, Parallel Computers, Adam Higler, 1981.

L.H. Jamieson,P.T. Mueller, and H.J. Siegel, "FFT Algorithm for SIMD Parallel

Processing Systems," Journal of Parallel and Distributed Computing, vol.3, 1986,

pp.48-71.

P.N. Swarztrauber, "Multiprocessor FFTs," Parallel Computing 5(1987), pp.197-

210.

[6]

[7]

[8]

[9]

[10]

[11]

[12]

D.H. Bailey, "A High-Performance Fast Fourier Transform Algorithm for the

CRAY-2", The Journal of Supercomputing, vol.1, 1987, pp.43-60.

R.A. Kamin III and G.B. Adams III, "Fast Fourier Transform Algorithm Design

and Tradeoffs on the CM-2," in Proceedings of the Conference on Scientific Applications on The Connection Machine, September 1988, pp.I34-160.

R.M. Chamberlain, "Gray Codes, Fast Fourier Transforms and Hypercubes," Parallel Computing 6 (1988), pp.225-233.

Thomas Schmiermund and Steven R. Seidel, "A Communication Model for the

Intel iPSC/2," Computer Science Technical Report CS-TR 90-02, Michigan Tech,

Univ, April 1990.

Steven R. Seidel, Ming-Horng Lee and Shivi Fotedar "Concurrent Bidirectional

Communication on the Intel iPSC/860 and iPSC/2," Computer Science Technical

Report CS-TR 90-06, Michigan Tech, Univ, November 1990.

D.H. Bailey, and P.O. Frederickson, "Performance Results for Two of the NAS

Parallel Benchmarks," in Supercomputing '91. (pp.166-173).

Shahid H. Bokhari, "Complete Exchange on the iPSC/860," ICASE Report No.

91-4, NASA Contractor Report No. 187498.

11
[13] Shahid H. Bokhari, "Multiphase Complete Exchange on a Circuit Switched Hypercube," (ICASE Report No. 91-5), Proceedings of the 1991 International Conference

on Parallel Processing, pp. 525-, 1991.

[14] D. W. Walker and a. Dongarra, "Design issues for a scalable library algebra subroutines," Draft Preprint, Oak Ridge National Laboratory October 1991.

12
,,m

REPORT DOCUMENTATION PAGE Fo,_ Approved OMB No 0704-0188

PuOhC feDor_l_ Duraen for lhqs _'_t=ec_on of _nfcrrr a,_l_r :_ PSt_'_a_ ,'3 _erage " Dour oet _esoonse JrcIL.oln_ the t_me _Or re_,ew_ng _nstrc.d_ons, $earch_r,g e_tst.ng data source_

, i . e ee , _.c; , _ _ -.,e , ; !,_ c." e , ,nf '_ t _ e r r bt" _ _ e t

Davis HigD_va,,, Suite !2C4. ,_ri_rgtcn,#,_ 22202-_']CP _rd t5 ','__!_-e _( ",¢_aoeme_i a_O Buc,_e_ P=mpe,c_ot.,_e_u.'_,lo_ Project C,_7C-'-0_8), Wash ngton. DC 20503

1. AGENCY USE ONLY (Leave blank) 2. REPORTDATE 13. REPORTTYPE AND DATES COVERED

November 1992 Contractor

4. TITLE AND SUBTITLE

A GENERAL PURPOSE SUBROUTINE FOR FAST FOURIER TRANSFORM

ON A DISTRIBUTED MEMORY PARALLEL MACHINE

6.AUTHOR{S)

A. Dubey, M. Zubair, and C.E. Orosch

7. PERFORMINGORGANIZATION NAME(S) AND ADDRESS(ES)

Institute for Computer Applications in Sclence

and Engineering

Mall Stop 132C, NASA Langley Research Center

Hampton, VA 23681-0001

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Adudnlstratlon

Langley Research Center

Hampton, VA 23681-0001

11.SUPPLEMENTARY NOTES

Langley Technical Monitor:

Final Report

Michael F. Card

i 12a. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Report

5. F'UNDING NUMBERS

C NASI-18605

C NASI-19480

WU 505-90-52-01

8. PERFORMING ORGANIZATION

REPORTNUMBER

ICASE Report No. 92-56

10. SPONSORING/MONITORING

AGENCY REPORT NUMBER

NASA CR-189724

ICASE Report No. 92-56

Submitted to International Conf.'

on Mathematlcal Modelling &

Scientific Comp., Dec. 7-11,

1992, Banza_ore

12b.DISTRIBUTION CODE

Subject Category 61

13. ABSTRACT(Maximum200wor_)

One issue which is central in developing a general purpose FFT subroutine on a distributed memory parallel machine is the data distribution. It is possible that

different users would llke to use the FFT routine with different data distributions.

Thus there is a need to design FFT schemes on distributed memory parallel machines

which can support a variety of data distributions. In this paper we present an FFT

Implementatlon on a distributed memory parallel machine which works for a number of

data distributions commonly encountered in scientific applications. We have also

addressed the problem of rearranging the data after computing the FFT. We have

evaluated the performance of our implementation on a distributed memory parallel

machine Intel iPSC/860.

!14. SUBJECTTERMS

distributed memory parallel machine, Fourier transform

17. SECURITYCLASSIFICATION

OF REPORT

Unclassified

NSN 7540-01-280-5500

18 SECURITYCLASSIFICATION

OF THIS PAGE

Unclasslfied

15. NUMBER OF PAGES

14

16. PRICECODE

A03

19. SECURITYCLASSIFICATION 20. LIMITATION OF ABSTRACT

OF ABSTRACT

Star_dard Form 298 (Rev 2-89)

Prescr_t:x_ b_ ANti 5_d Z]9-1B

298. t02

NASA-Langley, 1992