The pipeline addresses for $\overline{\overline{\mbox{Q}}}$ are
obtained by 
%% block-
blockwise      % % changed 3/14pm 
perfect shuffles (along the length of the pipeline)
of the addresses for $\overline{\overline{\mbox{S}}}$:
\begin{tabbing}
ttttabbb \=tab\=tab\=tab\=tab\=tab\=12345123451234512\= \kill
\>    {\bf function} $Y = \overline{\overline{\mbox{Q}}}( X, n )$ \\ 
\>    $A$ = $\overline{\overline{{\mbox{S}}}}$\_addresses$(r,s)$ \\
\\
\>    if $n > 2$ \\ 
\>    \>  for $ns = (n-2):-1:1 \;$ \\
%%                                 !perfect shuffle index exponent 
\>    \>    \>  $stride= r^ns$\\ 
\>    \>    \>  $block = r^{n-2-ns}\;\;\;$ ! block length \\ 
\>    \>    \>  $col2=0$\\ 
\>    \>    \>  for $k_1=1:stride$\\ 
\>    \>    \>  \> $col1 = (k_1-1)*block$\\ 
\>    \>    \>  \> for $k=1:r$\\  
\>    \>    \>  \> \> for $j=1:block$\\
\>    \>    \>  \> \> \> $B(:,col2)=A(:,col1)$\\
\>    \>    \>  \> \> \> $col1 = col1 + 1$\\
\>    \>    \>  \> \> \> $col2 = col2 + 1$\\
\>    \>    \>  \> \> end\\
\>    \>    \>  \> \> $col1 = col1 + (stride-1)*block$\\ 
\>    \>    \>  \> end\\
\>    \>    \>  end\\
\>    \>    \>  $A=B$\\
\>    \>  end\\
\>    end \\ 
\\
\>    for $row = 0:r-1$\\
\>    \>  $Y(row,:) = X(row,A(row,:))$ \\
\>    end  
\end{tabbing}       

\subsection{Twiddle recipes}
\label{sec:twiddlerecipes}
Every twiddle matrix ${\bf D}$ is diagonal, so it operates on a data
vector as a point-to-point vector multiply. Given some permutation
matrix ${\bf P}$, a new twiddle matrix ${\bf PDP}^T$ is equivalent to
a rediagonalizing of the vector shuffle of the diagonal of ${\bf D}$,
that is, ${\bf PDP}^T$ = diag$({\bf P}*$diag$({\bf D}))$.  (This is a
{\sc Matlab} notation: diag() puts the diagonal of a matrix in a
vector, and puts a vector in the diagonal of a matrix.)  Since we want
to perform shuffles within pipeline arrays, we reshape the twiddle
$N$-vector diag$({\bf D})$ as an $r \times s$ pipeline array $D$, just
as we originally reshaped the data vector.  Then we shuffle the
pipelined twiddles, to effect the equivalent of the vector shuffle
${\bf P}*$diag$({\bf D})$.  So we interpret the ${\bf PDP}^T$ operator
as an in-pipeline shuffle of the pipelined twiddles $D$, which are
then in position to operate on the pipelined data $X$ directly by
point-to-point multiplication, $Y = D.*X$. (As mentioned, the data will
actually be twiddled frame-by-frame in the pipelined implementation.)

We will interpret the twiddles expressed in (\ref{eqn:phasederivation}) this
way.  Each twiddle function below returns an $r \times s$ array $D$ of
twiddle factors (the actual twiddling of the data is not included):
\begin{tabbing}
ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill
\>    {\bf function} $D_{slow}$ = Dslow\_twiddles$( r, s )$ \\ 
\>    $\omega_{j} = exp(-2\pi i /r)$\\ 
\>    $t = 0$\\ 
\>    for $j = 0:(r-1)$\\
\>    \>  for $k = 0:(r^{\prime}-1)$\\
\>    \>    \>  $D_{slow}(:,t) = 
      (1, \omega_r^k, \omega_r^{2k}, ..., \omega_r^{(r-1)k})^T$ \\ 
\>    \>    \>  $t = t + 1$\\
\>    \>  end \\
\>    end 
\end{tabbing}       
\begin{tabbing}
ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill
\>    {\bf function} $D_{fast}$ = Dfast\_twiddles$( r, s )$ \\ 
\>    $\omega_{j} = exp(-2\pi i /r)$\\ 
\>    $t = 0$\\ 
\>    for $k = 0:(r^{\prime}-1)$\\
\>    \>  for $j = 0:(r-1)$\\
\>    \>    \>  $D_{fast}(:,t) = 
      (1, \omega_r^j, \omega_r^{2j}, ..., \omega_r^{(r-1)j})^T$ \\ 
\>    \>    \>  $t = t + 1$\\
\>    \>  end \\
\>    end 
\end{tabbing}       
The inverses of $D_{slow}$ and $D_{fast}$ are just their complex
conjugates, and are generated simply by replacing $\omega_{j}$ by
$\omega_{j}^{-1}$.
For stages $j = 1$:$n$ (counted down from $n$), we generate pipelined 
twiddles $\tilde{T}_j$ by 
\begin{tabbing}
ttttabbb \=tab\=tab\=tab\=1234512345123451234512\= \kill
\>    {\bf function} $\tilde{T}_j$ = 
                     $\tilde{\mbox{T}}$\_twiddles$( r, s, j )$ \\ 
\>    $\omega_{j} = exp(2\pi i /r^{j+1})$\\ 
\>    $\omega_{j}^{\prime} = \omega_{r^{j+1}}^r$\\ 
\>    for $k = 0:(r-1) \;\;$ ! direct sum loop\\ 
\>    \>  $t_1 = k \cdot r^{j-1}$\\ 
\>    \>  for $p = 0:(r-1)$\\ 
\>    \>    \>  $\tilde{T}_j(p,t_1) = \omega_{j}^{kp} $\\ 
\>    \>  end \\ 
\>    \>  $t_1 = t_1 + 1$\\ 
\>    \>  $t_2 = t_1 + r^{j-1}$\\ 
\>    \>  for $t = t_1:t_2\;\;$ ! fill next column from last \\ 
\>    \>    \>  $\tilde{T}_j(:,t) = \omega_{j}^{\prime k} \cdot 
                 \tilde{T}_j(:,t-1)$ \\ 
\>    \>  end \\ 
\>    end\\ 
\\
\>   if $j < n$ \\
\>    \>    $t_2 = r^j$\\
\>    \>    for $k = 0:(N/r^{j+1}) \;\;$ \\
\>    \>    \>  $t_1 = t_2$\\
\>    \>    \>  $t_2 = k \cdot r^j$\\
\>    \>    \>  $t = 0$\\
\>    \>    \>  for $t_0 = t_1:t_2$\\ 
\>    \>    \>    \>$\tilde{T}_j(:,t_0) = \tilde{T}_j(:,t)\;\;$ ! copy
columns\\ 
\>    \>    \>    \>$t = t + 1$\\ 
\>    \>    \>  end \\ 
\>    \>    end\\ 
\>    end 
\end{tabbing}       
The rest of the twiddle arrays can now be defined in terms of the
shuffles:
\begin{tabbing}
ttttabbb \=tab\=tabbbb\=tabbb \=tabbb1234512345123451234512\= \kill
\>    $D_{slow}^{\prime}$ \> \> = S$^{-1}( D_{slow}^{-1} )$ \\
\>    $D_{slow}^{\prime \prime}$ \> \> = Cslow$( D_{slow}^{\prime} )$ \\
\\
\>    $D_{fast}^{\prime}$ \> \> =  Cslow$( D_{fast}^{-1} )$ \\
\\
\>    $T_j$ \> \> = S$^{-1}( \tilde{T}_j )$ \\
\>    $T_j^{\prime}$ \> \> = Cslow$( T_j )$ \\
\\
\>    $D_1^{\prime}$ \> = 
           $\overline{\overline{\mbox{S}}}(D_{slow}^{-1}) 
                                  .* T_1^{\prime}.*D_{fast}^{-1}$\\
\>    if $1<j<n$\\
\>    \> $D_j^{\prime}$ \> = 
         $D_{slow}^{\prime \prime}.* T_j^{\prime}.*D_{fast}^{-1}$\\ 
\>    end \\ 
\>    $D_n^{\prime}$ \> = $D_{slow}^{\prime \prime}$ 
\end{tabbing}       

\section{Implementation issues} 
\label{sec:implementation}
In this section we describe issues that arise when the phase-rotation
FFT is implemented on a real parallel system.  In particular, we
describe implementation approaches for the radix-2 FFT on the iWarp
system. The main result is a scalable implementation of the pipelined
phase-rotation FFT that runs at the full 40 Mbytes/second rate of the
iWarp physical links.

\subsection{iWarp}
The iWarp is a private-memory multicomputer developed jointly by Intel
and Carnegie Mellon~\cite{iwarp}.  iWarp systems are 2-dimensional
tori of {\em nodes}, ranging in size from 4 to 1024 nodes. Each node
consists of an iWarp {\em component}, up to 16 Mbytes of off-chip
local memory, and a set of 8 unidirectional communication {\em links}
that physically connect the node to four neighboring nodes.  Each
component is a VLSI chip that contains a {\em processing agent} and a
{\em communication agent}. The processing agent is a general-purpose
load-store microprocessor that runs at a maximum rate of 20 MFLOPs.
Thus, a {\em clock}, or cycle time, is 50 ns.  The local memory is
accessed at a rate of 160 Mbytes/sec.  Each link runs at 40
Mbytes/sec, for a maximum aggregate bandwidth of 320 Mbytes/sec per
node.

The key feature of iWarp is its communication system, which is
summarized in Figure~\ref{fig:pathways}.  Each communication agent
contains a set of 20 hardware FIFO {\em queues}. Each queue can hold
up to 8 32-bit words, and can be accessed by user programs at the cost
of a register access. iWarp nodes communicate with other nodes using
unidirectional point-to-point structures called {\em pathways}. Each
pathway is a sequence of queues. Pathways can be created and destroyed
dynamically at runtime.  Data traveling along a pathway passes from
queue to queue {\em automatically}, without disturbing the
computations on intermediate nodes.  Multiple pathways can share the
same link by multiplexing in a round-robin fashion, one word at a
time. If each pathway has data to send, they are guaranteed a
proportional fraction of the link bandwidth. If only one pathway has
data to send, then it gets all of the link bandwidth.
\begin{figure}
%\epsfxsize=3.25in
\centerline{\epsfbox{pathways.eps}}
\caption{iWarp communication concepts.}
\label{fig:pathways}
\end{figure}

\subsection{Mapping strategies on iWarp}
The problem is to develop a mapping of the flowgraph in
Figure~\ref{fig:phase}(d) to an iWarp array.  The simplest mapping
strategy is to assign each flowgraph node to a unique processor node
of a linear array, route the flowgraph arcs through this array, and
then embed the resulting linear array in the iWarp torus. This
approach, called the PHASE5 mapping because it uses 5 iWarp nodes for
each FFT stage, is shown in Figure~\ref{fig:mappings}(a).

Each iWarp node in PHASE5 executes a small {\em node program} that
implements its flowgraph operator. Each twiddle node ($D^\prime$)
repeatedly reads a complex number from its input pathway 
multiplies it by the appropriate twiddle (precomputed off-line
using the recipes in Section~\ref{sec:twiddlerecipes}), and sends the
result to its output pathway.  Each shuffle
operator ($\overline{\overline{S}}$) repeatedly reads a complex data
item from its input pathway, stores it in memory, and uses the
appropriate shuffle index (again precomputed off-line using the
recipes in Section~\ref{sec:shufflerecipes}) to send an appropriate
double-buffered data point to the output pathway. The kernel node
($F$) repeatedly reads two complex numbers from its input pathways,
performs the radix-2 DFT kernel operation, and outputs two complex
numbers to its output pathways.

\begin{figure}
\epsfxsize=3.25in
\centerline{\epsfbox{mappings.eps}}
\caption{Strategies for mapping one stage of the FFT onto a linear array.
\rm{(a) PHASE5 mapping. (b) PHASE3 mapping.}}
\label{fig:mappings}
\end{figure}

Another approach, the PHASE3 mapping, combines the twiddle and shuffle
operators on a single node, as shown in Figure~\ref{fig:mappings}(b),
so that each stage requires 3 nodes instead of 5 nodes.  As we shall
see, the communication and computation throughputs of the two mappings
are identical.  The advantage of the PHASE3 mapping is that it is more
node-efficient, requiring fewer nodes per stage than the PHASE5
mapping.  The advantage of the PHASE5 mapping is its simplicity. Each
node is assigned exactly one operator from the flowgraph.

Figure~\ref{fig:phasenet} shows a working implementation of a
16K-point radix-2 phase-rotation FFT on a 64-node iWarp array at
Carnegie Mellon.  The implementation is based on the PHASE3 mapping
>from Figure~\ref{fig:mappings}(b).  The large squares are iWarp nodes,
labeled with the corresponding operator and stage number, where $D$ is
a twiddle/shuffle pair and $F$ is a kernel.  The small squares are
queues.  The arrows are iWarp pathways. As an artifact of our display
program intermediate queues are not drawn. Each of the 14 FFT stages
uses 3 nodes, with an additional 3 nodes for the parallel-pipeline
digit-reversing step at the end.

\subsection{Performance}
Each iteration of each node program in the PHASE3 and PHASE5 mappings
runs in at most 8 clocks. At the peak rate of 40 Mbytes/sec, each link
can produce and consume a 32-bit floating-point number every 2 clocks.
Further, each data point in the pipeline is a complex number
consisting of a pair of 32-bit floating-point words. As a result, each
pathway consumes exactly half of the available link bandwidth.  Since
each link is shared by two pathways, and since the iWarp communication
agent gives each pathway an equal share of the link bandwidth, without
disturbing the computations on intermediate nodes, each link is fully
utilized.  The result is a radix-2 FFT that runs at the full 40
Mbytes/sec rate of an iWarp link, regardless of the number of points
in the FFT! Since each sample consists of 8 bytes, the FFT runs at a
constant rate of $5~ \mbox{Msamples/sec}.$ Given a sufficient number
of nodes, the iWarp phase-rotation FFT's will produce arbitrarily
large FFT's at this rate.  Perhaps even more important, the
performance is the same on smaller FFT's.

\begin{figure}
\epsfxsize=3.25in
\centerline{\epsfbox{phasenet.eps}}
\caption{16K-point pipelined phase-rotation FFT running at 40 Mbytes/sec
(350 MFLOPS) on iWarp}
\label{fig:phasenet}
\end{figure}

Another way to characterize performance is by computational
throughput, expressed in millions of floating-point operations per
second (MFLOPS). However, there is a subtlety involved in using MFLOPS
as a performance measure. 
The iWarp phase-rotation FFT performs 16
floating-point operations %%% (FLOPs) 
per iteration per stage (2 adds and 4 multiplies 
%% by each twiddle operator, 
   by each of the two twiddle operators,    % % changed this para. 3/15/94
... 
and 4 adds by the kernel operator).  
   These 16 floating-point operations per iteration 
   reduce to 10, when one of the twiddles is always 1 
   and can be omitted, as in the radix-2 Pease and Cooley-Tukey FFTs. 
   This reduced figure results in 
the standard formula for computing FFT MFLOPS, 
$5N \log_{2} N$ floating-point operations 
per N-point FFT~\cite{carlson92}. 
   While one of the two twiddles is always $\pm 1$ in every 
   iteration of the radix-2 phase-rotation FFT, we have not yet 
   discovered an efficient, load-balanced pipeline mapping 
   that takes advantage of this fact. 
Therefore, 
   to compare the phase rotation FFT fairly with other FFTs, 
   we count its 16 FLOPs per iteration as equivalent to only 10 FLOPs. 
%%----nicer writing below but saving space...-----------------------
%% the phase rotation FFT trades off ease of pipelining for
%% number of floating point operations.  
%% to make a fair comparison with
%% other FFT algorithms, 
%% we must compute the phase rotation performance 
%% using the standard of 10 
%% floating-point operations per iteration per stage, 
%% even though the phase-rotation FFT is actually performing 16
%% floating-point operations per iteration per stage.

Since each node program executes its computation in at most 8 clocks,
and since each clock is $50$ ns, each stage of the iWarp
phase-rotation FFT runs at a rate of $25~\mbox{MFLOPS}$ for an
aggregate performance over all $\log N$ stages of $25 \log N~
\mbox{MFLOPS}$.  For example, the 16K-point FFT in
Figure~\ref{fig:phasenet} achieves a measured performance of $25
\times 14 = 350$ MFLOPS (single precision) on the iWarp systems at CMU. 
By comparison, a highly optimized 16K-point FFT has been measured at
237 MFLOPS (double precision) on a single-processor Cray
Y-MP~\cite[p.114]{carlson92}.  The numbers are not directly comparable
because of the different floating-point precisions, but they do
suggest that the absolute performance of the phase-rotation FFT is
quite good.

\section{Concluding remarks}
We have described an improved version of the Whelchel pipelined
phase-rotation FFT, developed recipes for generating the appropriate
twiddles and shuffle indices off-line and directly in terms of the
parallel pipeline, outlined mapping approaches for the radix-2 case on
the iWarp parallel computer, and presented measured performance
results of an implementation on iWarp.

The improvement on the original phase-rotation FFT is significant in
that it eliminates a potential pipeline bottleneck during the digit
reversing step at the end.  The twiddle and shuffle recipes should be
helpful to the programmer who wants to implement the pipelined phase
rotation FFT. The iWarp implementation validates a simple and
realistic approach for building scalable pipelined FFT's on a
programmable parallel system. Further, the implementation demonstrates
that, given a balanced parallel computer architecture with word-level
access to the communication links, it is possible to build FFT's that
run at the full link bandwidth of the links, even when the FFT's are
relatively small.

Other parallel systems are being considered as targets for the
multidimensional phase-rotation FFT. For example, the Maspar MP2
provides indirect addressing and routing capabilities, which would
facilitate fetching data into and out of the kernel FFTs and
performing data communication with large kernels.  As another example,
the Cray T3D multicomputer, like iWarp, provides direct, low-latency,
word-level access to the communication system, which would support the
fine grained parallelism found in the phase-rotation FFT.

\section*{Acknowledgements}
We would like to thank Tom Warfel, LeeAnn Tzeng, Doug Noll, and Doug
Smith their help and suggestions.

\small
\bibliography{/afs/cs/project/iwarp/member/droh/bib/compiler,/afs/cs/project
/iwarp/member/droh/bib/defs,/afs/cs/project/iwarp/member/droh/bib/me,/afs/cs
/project/iwarp/member/droh/bib/refs}

\end{document}

