\input{caption}
\makeatletter
\long\def\unmarkedfootnote#1{{\long\def\@makefntext##1{##1}\footnotetext{#1}}}
\makeatother

\newcommand{\ignore}[1]{}

\documentstyle[popl,picture]{article}
\begin{document}
\bibliographystyle{acm}
\pagestyle{myheadings}
\markboth
{\mbox{Whelchel, Withers, O'Hallaron, and Lieu - DRAFT - \today\ }}
{\mbox{Whelchel, Withers, O'Hallaron, and Lieu - DRAFT - \today\ }}

\title{Computational Aspects of the Pipelined Phase-Rotation FFT}

\author{
\begin{tabular}[c]{ c@{\extracolsep{2em}} c} 
Lang P. Withers, Jr., John E. Whelchel & David R. O'Hallaron, Peter J. Lieu \\ 
\\
{E-Systems, Inc., Melpar Div.} & {School of Computer Science}\\
{7700 Arlington Boulevard}     & {Carnegie Mellon University}\\     
{Falls Church, VA 22046}       & {Pittsburgh, PA 15213}
\end{tabular}
}

\date{} 
\maketitle{}
\thispagestyle{empty}
\unmarkedfootnote{
Supported in part by the Advanced Research Projects Agency,
Information Science and Technology Office, under the title "Research
on Parallel Computing," ARPA Order No. 7330.  Work furnished in
connection with this research is provided under prime contract
MDA972--90--C--0035 issued by ARPA/CMO to Carnegie Mellon University,
and in part by the Air Force Office of Scientific Research under
Contract F49620--92--J--0131.

%\vspace {1in}  %ACM copyright notice
}

\begin{abstract}
The phase-rotation FFT has a simple, constant-geometry,
parallel-pipeline architecture.  This paper considers computational
aspects of implementing the FFT on a parallel computer system.
Shuffle address and twiddle recipes, described directly in terms of
the pipeline, are provided, and a recent full-bandwidth implementation
on the iWarp parallel computer system is described.  A new,
parallel-pipeline equivalent of the index-reversing shuffle is added
to complete the original phase-rotation FFT design.
\end{abstract}

\section {Introduction}
The Fast Fourier Transform (FFT) is an important algorithm with many
applications in signal processing and scientific computing. The
Whelchel phase-rotation FFT \cite{whelchel90} is a pipelined version
of the Pease constant-geometry FFT \cite{pease68}, which itself is
derived from the original Cooley-Tukey FFT \cite{cooley65}.  

The phase-rotation FFT of radix $r$ is designed for a pipeline of $r$
parallel data channels.  At each time step, in each stage, the
pipeline carries the next $r$ data points, one from each channel, into
a DFT kernel.  Unlike earlier pipelined FFT's
\cite{mcclellan78,corinthios71}, the phase-rotation FFT has the key
property that {\em no data is switched across channels, except within
the DFT kernel and at the input and output}.

The phase-rotation FFT extends easily to higher radices, reducing
memory and latency while preserving high throughput and parallel
shuffling simplicity of lower radix versions.  The phase-rotation FFT
has also been extended to a vector-radix, multidimensional
parallel-pipeline FFT with the same qualities of the one-dimensional
algorithm, and without transposes \cite{withers91}.

This paper reports some results from a project to implement the 1D
phase-rotation FFT on a parallel computer system.  There are three
main results:
\begin{itemize}

\item The bit reversing shuffle step in the original version
of the phase-rotation FFT \cite{whelchel90} is a potential pipeline
bottleneck. We describe a new version that corrects this problem 
by using a parallel-pipeline bit-reversing step.

\item While the structure of the phase-rotation FFT is extremely
simple, we soon learned that generating the appropriate twiddles and
shuffle indices from the original tensor algebraic formulation
\cite{whelchel90} is quite difficult, even for the designers of the
algorithm! To try to help the implementer, we have reformulated the
phase-rotation FFT. Following \cite{vanloan92}, we present a new set
of recipes, in a {\sc Matlab}-like format, for generating the twiddles
and shuffle indices directly in terms of the parallel pipeline.

\item Finally, we describe an implementation of the phase-rotation
FFT on the iWarp, a private-memory multicomputer system developed by
Intel and Carnegie Mellon \cite{iwarp,iwarpcomm}. In particular, we
demonstrate that the 1D radix-2 phase-rotation FFT can run at the full
40 Mbytes/sec rate of the iWarp physical links.
\end{itemize}

Section~\ref{sec:basic} introduces the phase-rotation concept. 
Section~\ref{sec:phasedef} formally defines the improved FFT 
algorithm. Section~\ref{sec:recipes} gives the recipes for generating
the twiddles and shuffle indices in terms of the parallel pipeline.
Finally, Section~\ref{sec:implementation} describes the full-bandwidth
implementation on iWarp.

\section{The basic idea} 
\label{sec:basic}
This section introduces the concept of the phase-rotation FFT.
Starting with the Pease constant-geometry FFT, we informally derive
the pipelined phase-rotation FFT, identifying the key insights along
the way.

\subsection{Constant-geometry FFT}
Figure~\ref{fig:phase}(a) shows the flowgraph for a radix $r$ $N$-point
decimation-in-frequency (DIF) constant-geometry FFT, with $r=2$ and
$N=r^n=8$. There are $n$ stages. Each stage consists of $N/r$ {\em
kernels}, each of which is an operator that inputs two complex numbers
and outputs two complex numbers.
\begin{figure}
\centering
\makebox{
\scaledpict 50mm by 29mm (phase scaled 1000)
}
\caption{Derivation of the phase-rotation FFT. {\rm (a) Initial constant-geometry FFT. (b) Pipelined constant geometry FFF. (c) Pipelined FFT based on cyclic rotations. (d) Phase-rotation FFT.}
} 
\label{fig:phase}
\end{figure}

Each stage in the constant-geometry performs an identical perfect
stride-by-$s$ shuffle of its data vector, where $N=rs$.  If the data
vector is regarded as an $s \times r$ array, then the perfect shuffle
is a transpose operation. For example, the following transpose is a
stride-by-4 perfect shuffle, for $N=8$ points and radix $r=2$:
\begin{center}
\begin{tabular} {l c c c}                      
$\left[ 
\begin{array} {c c} 
 0 & 4 \\
 1 & 5 \\
 2 & 6 \\
 3 & 7 
 \end{array} \right]$
&
 $\stackrel{T}{\longrightarrow}$
&
 $\left[ \begin{array} {c c c c} 
 0 & 1 & 2 & 3 \\
 4 & 5 & 6 & 7
\end{array} \right]$\\ 
\end{tabular} 
\end{center}
The data items, labeled by their indices in the original column vector,
are regarded as equivalent to a $4 \times 2$ array composed by a
stride-by-4 unstacking of the 8-point column vector.  After the
transpose, the $2 \times 4$ array is equivalent to a new 8-point
column vector composed by a stride-by-2 stacking. As we shall see,
this transpose creates difficulties when we try to pipeline the
constant-geometry FFT. And it is precisely these difficulties that the
phase-rotation FFT addresses.

\subsection{Pipelining the FFT}
The constant-geometry FFT can be pipelined by computing each stage on
a single processor. For example, Figure~\ref{fig:phase}(b) shows the
pipeline for a single stage. The pipeline consists of a sequence of
operators connected by {\em pipeline segments}. Each pipeline segment
consists of $r$ parallel {\em channels}. Each channel consists of a
{\em stream} of $N/r$ data points, which are labeled in this example
by their indices from the original column vectors in
Figure~\ref{fig:phase}(a). For each pipeline segment, the $r$ data points
in the same position in each stream are known as an $r$-{\em frame},
or simply, a {\em frame}. For example, in Figure~\ref{fig:phase}(b),
the first frame in the pipeline segment between $S$ and $F$ is (0,4),
the second frame is (1,5), and so on.

Each time step, the twiddle operators ($W$) collectively read a frame,
perform an element-wise complex multiplication, and write the
resulting frame.  Notice that each stream is operated on
independently.  Similarly, the kernel operator ($F$) reads a frame,
computes the radix-$r$ kernel, and writes the resulting frame.
However, in this case, the streams are not independent; each data item
in the output frame is a function of every data point in the input
frame.

The twiddle and kernel operators pipeline nicely because during each
time step they read and write a single frame.  However, the pipelined
shuffle operator ($S$) is less well behaved.  It must read multiple
input frames before it can begin to produce the first output frame.
For example, for the radix-2 FFT in Figure~\ref{fig:phase}(b), the shuffle
operator must read half the frames before it can begin producing any
output frames. For the radix-4 FFT, the situation is worse, the shuffle
operator must read {\em all} of the frames before it can produce a
single output frame.

The pipelined shuffle operator creates a number of problems for the
implementer.  In general the shuffle operator must read and store the
entire pipeline segment before it can begin to produce any outputs.
Since double buffering must be used to guarantee full bandwidth, each
shuffle operator must have enough RAM to store $2N$ data points.
Further, this memory must be $r$-ported. As the FFT radix increases,
building such memories become increasingly expensive. Notice that the
kernel operator suffers from a similar problem, but to a much lesser
extent because it operates only on single frames.

\subsection{The phase-rotation concept}
The phase-rotation FFT attacks the problems caused by the perfect
shuffle in the pipelined constant-geometry FFT. The idea is based on a
simple factorization of the perfect shuffle operator. Once in this
form, the DFT cyclic shift theorem can be applied to replace cyclic
rotations with phasor multipliers, which are then merged into the
twiddles. The result is a pipelined FFT where the only communication
between streams occurs at each kernel and during the bit reversal at
the end. Further, the communications that do exist only involve data
points contained in the same frame. The rest of this section provides
the intuition for these central ideas.

When represented in matrix form, the shuffle operator $S$ in
Figure~\ref{fig:phase}(b) can be factored into the double matrix
product of the form
\begin{equation}
{\bf S} = {\bf C}_{fast} {\overline{\overline{\bf S}}} {\bf C}_{slow}
\label{eqn:decomp}
\end{equation}
where ${\overline{\overline{\bf S}}}$ is a {\em parallel-pipeline
shuffle} that shuffles independently along each input stream, and
${\bf C}_{slow}$ and ${\bf C}_{fast}$ are frame-wise (i.e., involving
only data items in a single frame) cyclic rotations.  We will formally
define ${\overline{\overline{\bf S}}}$, ${\bf C}_{slow}$, and ${\bf
C}_{fast}$ in Section~\ref{sec:phasedef}. For now, we will rely on the
example in Figure~\ref{fig:decomp} to get the idea across.
\begin{figure}
\centering
\makebox{
\begin{tabular} {l c c c}                      
& 
 $\left[ \begin{array} {c c c c} 
 0 & 2 & 4 & 6 \\
 1 & 3 & 5 & 7
 \end{array} \right]$ 
&
 $\stackrel{{\bf S}}{\longrightarrow}$ 
&
 $\left[ \begin{array} {c c c c} 
 0 & 1 & 2 & 3 \\
 4 & 5 & 6 & 7
 \end{array} \right]$ 
\\ 
&
\\
& 
$\downarrow \; {\bf C}_{slow}$ & & $\uparrow \; {\bf C}_{fast}$ 
\\ 
& 
\\
&
 $\left[ \begin{array} {c c c c} 
 0 & 2 & 5 & 7 \\
 1 & 3 & 4 & 6             %% needs hand-drawn circles, squares 
 \end{array} \right]$ 
&
 $\stackrel{ \overline{\overline{\bf S}} }{\longrightarrow}$ 
&
 $\left[ \begin{array} {c c c c} 
 0 & 5 & 2 & 7 \\
 4 & 1 & 6 & 3             %% needs hand-drawn circles, squares 
\end{array} \right]$ 
\end{tabular} 
}
\caption{Interpretation of ${\bf S} = {\bf C}_{fast} {\overline{\overline{\bf S}}} {\bf C}_{slow}$}
\label{fig:decomp}
\end{figure}
In this example, a pipeline segment is represented as a matrix. Each
row in the matrix is a stream and each column is a frame.  Notice that
the operators in (\ref{eqn:decomp}) are applied from right to left,
and not left to right. 

The first step in Figure~\ref{fig:decomp} is a cyclic rotation ${\bf
C}_{slow}$, which rotates each frame.  This rotation is frame-wise, in
the sense that only data points contained in the same frame are
rotated across the streams.  Next, a parallel-pipeline shuffle
permutes the data in each stream.  Notice that no data points need to
be transferred between streams in this step. The last step is a cyclic
rotation ${\bf C}_{fast}$, which also performs a frame-wise rotation.
If we apply (\ref{eqn:decomp}) to the pipelined FFT in
Figure~\ref{fig:phase}(b), we get a pipelined FFT based on cyclic
rotations, as shown in Figure~\ref{fig:phase}(c).

The basic rotation in Figure~\ref{fig:phase}(c) that is applied at
slowly varying, then fast-varying rates, is represented by the $r
\times r$ cyclic (circular) shift permutation matrix ${\bf C}_{r}$,
made by permuting the rows of the identity matrix down by one row,
and moving the bottom row up to the top.  
\ignore{
For example,
\begin{eqnarray*}
{\bf C}_{4} = \left(
\begin{array}{ c c c c }
0&0&0&1\\
1&0&0&0\\
0&1&0&0\\
0&0&1&0
\end{array} \right)\;\;.
\end{eqnarray*}
}%end ignore
The key insight of the phase-rotation FFT is that the cyclic shift
theorem for the DFT can be applied to the cyclic shift operators in
Figure~\ref{fig:phase}(c).  In matrix form, the cyclic shift theorem for
a DFT is the relation
\begin{equation}
{\bf F}_{r}{\bf C}_{r} = {\bf D}_{r}{\bf F}_{r},
\label{eqn:shifttheorem}
\end{equation} 
where 
${\bf D}_{r}= diag( 1, \omega, \omega^{2}, ..., \omega^{r-1} )$
is a set of twiddles, and the DFT matrix of size $r$ is
\begin{displaymath}
{\bf F}_{r} = \frac {1}{\sqrt{r}} ( \omega_{r}^{jk})_{j,k=0}^{r-1} ,
\end{displaymath}
where $\omega = e^{-\frac{2 \pi i}{r}}$.  For the pipelined FFT,
(\ref{eqn:shifttheorem}) says that phasor multipliers after a DFT
kernel give the same effect as physical data rotations before the DFT
kernel.  Likewise, physical rotations after the kernel are equivalent
to phasor multipliers before it. This meaning of
(\ref{eqn:shifttheorem}) is shown graphically in
Figure~\ref{fig:theorem} for a pipelined radix-2 kernel.
\begin{figure}
\centering
\makebox{
\scaledpict 50mm by 29mm (theorem scaled 1000)
}
\caption{Interpretation of ${\bf F}_{r}{\bf C}_{r} = {\bf D}_{r}{\bf F}_{r}$}
\label{fig:theorem}
\end{figure}

The key point is that the cyclic rotations in Figure~\ref{fig:phase}(c)
can be replaced by constant phasor multipliers.  These phasors can
then be absorbed by the twiddle factors on either side of the kernel,
leaving only a parallel-pipeline shuffle.  The result is the pipelined
phase-rotation FFT, which is shown in Figure~\ref{fig:phase}(d).  This
completes our informal derivation of the phase-rotation FFT.

The structure of the phase-rotation FFT in Figure~\ref{fig:phase}(d) is
quite similar to the original pipelined FFT in
Figure~\ref{fig:phase}(b). The twiddle operators ($W^\prime$) are
identical to the original twiddle operators ($W$), except now the
twiddles incorporate the original twiddles, phasors for the ${\bf
C}_{fast}$ operator from the previous stage, and phasors for the ${\bf
C}_{slow}$ operator from the next stage. The kernel operators are
identical as well. The important difference is that the troublesome
perfect shuffle operator has now been replaced by a parallel-pipeline
shuffle that requires no communication across the streams.

There are several other important things to note about the
phase-rotation FFT in Figure~\ref{fig:phase}(d).  First, there are no
additional multiplications or additions, compared to the original
pipelined FFT. Second, the only communication across streams occurs at
the kernel, and this communication is constrained, in the sense that
only data points within a single frame need to be switched across
channels.

\section{Improved phase-rotation FFT}
\label{sec:phasedef}
In this section we give a formal definition of an improved version of
the original phase-rotation FFT \cite{whelchel90}. The new version
replaces the bit reversing permutation in the original phase rotation
FFT with a parallel-pipeline shuffle, followed by a frame-wise cyclic
rotation.  The advantage of this new approach is that during the bit
reversal step at the end, all communication between streams is limited
to data points within a single frame.

For radix $r$ and $N=r^n$ points ($n>1$), the 1-dimensional
phase-rotation FFT is a matrix factorization of the $N$-point DFT
matrix ${\bf F}_{N}$.  Starting with the Pease constant-geometry
factorization, we replace its perfect shuffles ${\bf S}$ by ${\bf S} =
{\bf C}_{fast} \overline{\overline{{\bf S}}} {\bf C}_{slow}$.
Similarly, at the left end we replace the radix-$r$
index-digit-reversing permutation ${\bf Q} = {\bf Q}_{N,r}$ of $N$
data points by ${\bf Q} = {\bf C}_{slow}^T \overline{\overline{\bf
Q}}{\bf C}_{slow}$, where $\overline{\overline{\bf Q}}$ is another
parallel-pipeline shuffle that will be defined formally in
Section~\ref{sec:recipes}. The phase-rotation FFT is then defined by:
\begin{eqnarray}
{\bf F}_{N} &=& {\bf Q} \cdot 
\prod_{j=1}^{n} \left( {\bf F}{\bf S}{\bf T}_{j} \right)
 = \cdots 
\left( \mbox{ \small
\begin{tabular}[c]{c} 
\small vigorous\\
\small algebraic\\
\small shuffling
\end{tabular}  } 
\right) \cdots 
\nonumber\\ 
%&=& 
%\left( \prod_{j=1}^{n-1} {\bf F}
%{\bf C}_{fast} \overline{\overline{{\bf S}}} 
% {\bf C}_{slow} {\bf T}_{j} \right) {\bf F}{\bf S}{\bf Q} 
% = \cdots \nonumber\\
&=& {\bf C}_{slow}^T \cdot \overline{\overline{\bf Q}}{\bf D}_{fast}^{\prime} 
\left[  
\prod_{j=1}^{n} 
\left( {\bf F} \overline{\overline{{\bf S}}}{\bf D}_{j}^{\prime} \right) 
\right]
\cdot {\bf C}_{slow} . 
\label{eqn:phase} 
\end{eqnarray}
%
As before, we let $s = N/r$ and $r^{\prime} = N/r^2$.  ${\bf F}$ is a
direct (tensor, Kronecker) product ${\bf I}_{s} \otimes {\bf F}_{r} =
\mbox{diag}({\bf F}_{r}, {\bf F}_{r}, ...,{\bf F}_{r})$.  We interpret
this as a kernel DFT ${\bf F}_{r}$ operating on $s$ successive
frames of $r$ points placed in the pipeline.  For $j=1:n$, the other
parts of (\ref{eqn:phase}) are defined by
\begin{eqnarray}
{\bf C}_{slow} &=&  \bigoplus_{k=0}^{r-1} 
   \left( {\bf I}_{r^{\prime}} \otimes {\bf C}_{r}^{k} \right ) \nonumber\\ 
{\bf C}_{fast} &=&  {\bf I}_{r^{\prime}} \otimes 
   \left(   \bigoplus_{k=0}^{r-1} ({\bf C}_{r}^{T})^{k} \right ) \nonumber\\ 
\omega_j &=& exp \left({-\frac{2 \pi i}{r^{j}}} \right) \nonumber\\
{\bf D}_{r} &=& 
 diag( 1, \omega_{1}, \omega_{1}^{2}, ..., \omega_{1}^{r-1} ) 
 \nonumber\\ 
{\bf D}_{r^{j+1}} &=& 
 diag( 1, \omega_{j+1}, \omega_{j+1}^{2}, ..., \omega_{j+1}^{r^{j}-1} ) 
 \nonumber\\ 
{\bf D}_{slow}^{-1} &=&  \bigoplus_{k=0}^{r-1} 
   \left( {\bf I}_{r^{\prime}} \otimes {\bf D}_{r}^{-k} \right) 
   \nonumber\\ 
{\bf D}_{slow}^{\prime} &=&
              {\bf C}_{fast}^T {\bf D}_{slow}^{-1} {\bf C}_{fast} \nonumber\\
{\bf D}_{slow}^{\prime \prime} &=& 
     \overline{\overline{{\bf S}}}^T 
     {\bf D}_{slow}^{\prime} 
     \overline{\overline{{\bf S}}} \nonumber\\
{\bf D}_{fast}^{-1}  &=&  {\bf I}_{r^{\prime}} \otimes 
   \left(   \bigoplus_{k=0}^{r-1} {\bf D}_{r}^{-k} \right) 
   \nonumber\\ 
{\bf D}_{fast}^{\prime} &=&
              {\bf C}_{slow}^T {\bf D}_{fast}^{-1}  {\bf C}_{slow} \nonumber\\
{\tilde{\bf T}}_{j} &=& {\bf I}_{ \frac{N}{r^{(j+1)}} } \otimes 
   \left(  \bigoplus_{k=0}^{r-1}{\bf D}_{r^{j+1}}^{k} \right) \nonumber\\
%
{\bf T}_{j} &=& {\bf S}^{j\;T}{\tilde{\bf T}}_{j}{\bf S}^{j} \nonumber\\
{\bf T}_{j}^{\prime} &=& 
                 {\bf C}_{slow} {\bf T}_{j} {\bf C}_{slow}^T \nonumber\\
{\bf D}_{1}^{\prime} &=&      
     \left( \overline{\overline{{\bf S}}}^T 
            {\bf D}_{slow}^{-1} 
            \overline{\overline{{\bf S}}} \right) \cdot 
{\bf T}_{1}^{\prime}{\bf D}_{fast}^{-1}  \nonumber\\
{\bf D}_{j}^{\prime} &=& 
   {\bf D}_{slow}^{\prime \prime} {\bf T}_{j}^{\prime} {\bf D}_{fast}^{-1} \;\;,
\;\;\;\; j=2:n-1 \nonumber\\ 
{\bf D}_{n}^{\prime} &=& 
{\bf D}_{slow}^{\prime \prime} {\bf T}_{n}^{\prime} = 
{\bf D}_{slow}^{\prime \prime} \;\;. 
\label{eqn:phasederivation}
\end{eqnarray}
The direct sums are of the form 
\begin{displaymath}
\bigoplus_{k=0}^{r-1} {\bf A}_{k} =
\mbox{diag}( {\bf A}_{0}, {\bf A}_{1},..., {\bf A}_{r-1} ). 
\end{displaymath}
See \cite{withers91} for more on the basic definitions and relations
used to derive (\ref{eqn:phase}), as well as the generalization to
higher dimension FFT's.

Note that the stages in (\ref{eqn:phase}) are counted in reverse time
order by the index $j$. This is in keeping with the fact that
(\ref{eqn:phase}) is a decimation-in-frequency (DIF) version of the FFT.
The transpose of (\ref{eqn:phase}), with the product
$\prod_{j=n}^{1}$, is the decimation-in-time (DIT) version of the
phase-rotation FFT.

A ${\bf C}_{slow}$ shuffle and its inverse remain at the input and
output ends of the pipeline, respectively.  As we have seen, ${\bf
C}_{slow}$ is a completely frame-wise rotation.  It rotates
(commutates) the data within each successive frame (column $r$-vector)
of the $r \times s$ pipeline segment for a stage.  There is also an
implicit frame-wise broadcast within each FFT kernel engine, when an
$r$-point DFT is somehow computed.  So in the phase-rotation FFT, data
motion is all parallel, except for frame-wise motions at I/O and at
every FFT kernel.  The simplicity of the phase-rotation FFT is that no
data point ever moves both down and across the pipeline in one
time-step.

\section{Pipeline recipes}
\label{sec:recipes}
While the structure of the pipelined phase-rotation FFT is extremely
simple, experience has taught us that generating the appropriate
twiddles and shuffle indices from the matrix formulations of
(\ref{eqn:phase}) and (\ref{eqn:phasederivation}) is difficult and
confusing.  To address this problem, we have developed a collection of
recipes, in a {\sc matlab}-like format, for generating the
phase-rotation twiddles and shuffle indices off-line.  The recipes are
defined for any 1D phase-rotation FFT with a power of two radix.

As we saw in (\ref{eqn:phase}), the pipelined phase-rotation FFT
performs a typical ``twiddle, shuffle, kernel'' cycle at each stage.
Only the twiddles vary from stage to stage, and there is a
digit-reversing shuffle at the end.  To implement this FFT using
parallel $r \times s$ pipeline segments (one per stage), we insert the
$N$-vector of input data ${\bf x}$ into the pipeline as an $r \times
s$ array $X$: the first $r$ points of ${\bf x}$ go into the first
frame (column) $X$, the second $r$ points go into the second frame,
and so on.  We must also have a shuffle address and a twiddle factor
ready for each point in the pipeline.  In other words, we would like
to fill one $r \times s$ copy $A$ of the pipeline segment with
addresses, and another copy $D$ with twiddles.

With this information in hand, the behavior of the processors that
execute each stage, ${\bf F} \overline{\overline{{\bf S}}}{\bf
D}_{j}^{\prime}$, in the pipeline is completely determined.  Using the
current frame of addresses, each processor will fetch the current
$r$-frame of data $X(0$:$r-1, A(0$:$r-1,t))$, and the current
$r$-frame of twiddles $D(0$:$r-1, A(0$:$r-1,t))$ (pointwise in
parallel), multiply these two frames pointwise, and then do an
$r$-point DFT ${\bf F}_r$ of the twiddled data frame.  

The twiddle and shuffle in this section are recipes are
``in place'' in the sense that they work inside the $r \times s$
pipeline segments that will contain the desired addresses and
twiddles. They are not ``in place'' in the usual sense, as we will
freely use an input and an output copy of a pipeline segment.  This
approach avoids constructing and operating with large $N \times N$
matrices (each containing only $N$ non-zero elements).  Each
parallel-pipeline function recipe is given a name similar to that of
the $N \times N$ matrix factor in the FFT (\ref{eqn:phase}) that it
effectively implements.  

\subsection{Shuffle recipes}
As a convention, pipeline addresses (pipeline array row and column
indices) run 0:$r-1$ and 0:$s-1$, respectively.  To do
parallel-pipeline shuffles, we only need the horizontal (column)
addresses, since the data inside each pipe will only jump within that
stream (row).  The cross-stream shuffles, Cslow and Cfast, are
implemented using $\pi_{r}$, a cyclic rotation of a frame (a vertical
slice of the parallel pipeline) that has the effect of ${\bf y}_r =
{\bf C}_r{\bf x}_r$.  $\pi_{r}$ takes a column $r$-vector ${\bf x}_r =
(x_0, x_1, x_2, ..., x_{r-1})^T \mapsto {\bf y}_r = (x_{r-1}, x_0,
x_1, ..., x_{r-2})^T$.
\begin{tabbing}
ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill
\>    {\bf function} $Y$ = Cslow$( X )$ \\
\>    $col = 0$\\
\>    for $k = 1:r$\\
\>    \>  for $j = 1:r^{\prime}$\\
\>    \>    \>  $Y(:,col) = \pi_{r}^{k} (X(:,col))$\\
\>    \>    \>  $col = col + 1$\\
\>    \>  end \\
\>    end
\end{tabbing}       
\begin{tabbing}
ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill
\>    {\bf function} $Y$ = Cfast$( X )$ \\
\>    $col = 0$\\
\>    for $j = 1:r$\\
\>    \>  for $k = 1:r^{\prime}$\\
\>    \>    \>  $Y(:,col) = \pi_{r}^{k} (X(:,col))$\\
\>    \>    \>  $col = col + 1$\\
\>    \>  end \\
\>    end
\end{tabbing}       
The inverses of Cslow and Cfast are formed by simply reversing
$\pi_{r}$.  Next, we define some perfect shuffles.
\begin{tabbing}
ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill
\>    {\bf function} $Y$ = S$( X ) \;\;$ !stride by $s$\\
\>    $col = 0$\\
\>    for $row = 0:r-1$\\
\>    \>  for $k1 = 0:r:s-r$\\
\>    \>    \>  $k2 = k1+r-1$\\
\>    \>    \>  $Y(row,k1:k2) = X(:,col)$\\
\>    \>    \>  $col = col + 1$\\
\>    \>  end \\
\>    end
\end{tabbing}       
\begin{tabbing}
ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill
\>    {\bf function} $Y$ = S$^{-1}( X )\;\;$ !stride by $r$\\
\>    $col = 0$\\
\>    for $row = 0:r-1$\\
\>    \>  for $k1 = 0:r:s-r$\\
\>    \>    \>  $k2 = k1+r-1$\\
\>    \>    \>  $Y(:,col) = X(row,k1:k2)$\\
\>    \>    \>  $col = col + 1$\\
\>    \>  end \\
\>    end 
\end{tabbing}       
To implement the parallel-pipeline shuffles, 
$\overline{\overline{{\mbox{S}}}}$, 
$\overline{\overline{{\mbox{S}}}}^{-1}$, and 
$\overline{\overline{{\mbox{Q}}}}$, 
we will use the parallel-pipeline addresses $A$, which are 
computed by the following function:
\begin{tabbing}
ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill
\>    {\bf function} $A$ = 
           $\overline{\overline{\mbox{S}}}$\_addresses$(r,s)$ \\
\>    ${\bf a} = ( 0, r^{\prime}, ..., (r-1)r^{\prime})^T$  \> \> \>\\
\>    $col = 0$\\
\>    for $j = 1:r^{\prime}$\\
\>    \>  for $k = 1:r$\\
\>    \>    \>  $A(:,col) = {\bf a}$\\
\>    \>    \>  $col = col + 1$\\
\>    \>    \>  ${\bf a} = \pi_{r}({\bf a})$ \> \\
\>    \>  end \\
\>    \>  ${\bf a} = {\bf a} + {\bf 1}_{r}$ \\ %% \;\;\;(mod\;\;\; r)$ \> \> \\
\>    end 
\end{tabbing}       
Looking closely, one can see Cfast$^{-1}$ at work producing the
addresses $A$ in the last function.  The addresses $A$ can also be
generated by loading a pipeline segment with simple $r \times r$
address blocks $B_{rr}$, and then applying Cfast$^{-1}$ to the
pipeline segment.  The first block to load is $B_{rr} =$
diag$(0$:$r^{\prime}$:$s-1)*{\bf 1}_{rr}$, where ${\bf 1}_{rr}$ is the
$r \times r$ matrix whose entries are all 1's.  The next block is
always $B_{rr} = B_{rr} + {\bf 1}_{rr}$, until the pipeline segment
contains $r^{\prime}$ blocks and is full.
\begin{tabbing}
ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill
\>    {\bf function} $Y = \overline{\overline{{\mbox{S}}}}( X )$ \\
\>    $A$ = $\overline{\overline{{\mbox{S}}}}$\_addresses$(r,s)$ \\
\>    for $row = 0:r-1$\\
\>    \>  $Y(row,:) = X(row,A(row,:))$ \\
\>    end 
\end{tabbing}       
\begin{tabbing}
ttttabbb \=tab\=tabbb \=tabbb1234512345123451234512\= \kill
\>    {\bf function} $Y = \overline{\overline{\mbox{S}}}^{-1}( X )$ \\
\>    $A$ = $\overline{\overline{{\mbox{S}}}}$\_addresses$(r,s)$ \\
\>    $[AA,I] =$ sort($A$)  \\
\>    for $row = 0:r-1$\\
\>    \>  $Y(row,:) = X(row,I(row,:))$ \\
\>    end 
\end{tabbing}       
In the above functions, sort($A$) sorts each row of an array $A$ in
ascending order. It returns the row-sorted array $AA$ and the
corresponding array of addresses $I$ where the successive row elements
were found in $A$.  After we have sorted the addresses $A$ for
$\overline{\overline{{\mbox{S}}}}$, $I$ has the addresses for
$\overline{\overline{\mbox{S}}}^{-1}$.

The pipeline addresses for $\overline{\overline{\mbox{Q}}}$ are
obtained by block-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}
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 the issues that arise when the
phase-rotation FFT is implemented on a 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 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
Carnegie Mellon and Intel Corporation \cite{iwarp,iwarpcomm}.  iWarp
systems are 2-dimensional tori of iWarp nodes, ranging in size from 4
to 1024 node. 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.

The iWarp 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, centered around a $32
\times 128$-bit register file, that runs at a maximum rate of 20
MFLOPS. 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 the 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.  iWarp nodes communicate with other nodes using
unidirectional point-to-point structures called {\em pathways}. Each
pathway is a sequence of queues that can be defined dynamically at
runtime.  Figure~\ref{fig:pathways} shows a pair of such pathways.
\begin{figure}
\centering
\makebox{
\psfig{figure=pathways.eps}
}
\caption{iWarp communication structures.}
\label{fig:pathways}
\end{figure}

Data traveling along a pathway passes from queue to queue {\em
automatically}, without disturbing the computations on intermediate
nodes. For example, in Figure~\ref{fig:pathways}, data traveling over
the pair of pathways does not disturb the computation on node 1.  The
latency from queue to queue is small, ranging from 100-300
nanoseconds.

Multiple pathways can share the same link. For example, in
Figure~\ref{fig:pathways}, two pathways share the link from node 1 to
node 2. In this case, the pathways share the link bandwidth in a
round-robin fashion, one word at a time.  If only one pathway is
sending data over a link link, then it gets the entire link bandwidth.
If multiple pathways are sending data over a link, then the link is
utilized at the full 40 Mbytes/sec, and each pathway is guaranteed a
proportional fraction of the bandwidth.

User programs can directly access the queues, one word at a time, by
reading and writing special registers in the register file called {\em
gates}.  To an iWarp instruction, a gate is just another register in
the register file. The important point is that a program can read or
write a word in a queue with the latency of a register access. A
single instruction can read and write up to 4 words from queues, with
a maximum aggregate bandwidth of 160 Mbytes/sec. Gates are accessed
from user-level C programs using \verb.send_word(). and
\verb.recv_word(). primitives.

\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 objective is a
pipelined radix-2 FFT that runs at the maximum 40 Mbytes/sec rate of
an iWarp communication link. We will also consider a mapping that uses
fewer nodes per stage, at the price of lower bandwidth. In the final
paper, we plan to treat the radix-4 case as well.

The simplest 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).
\begin{figure}
\centering
\makebox{
\psfig{figure=mappings.eps}
}
\caption{Mapping strategies.}
\label{fig:mappings}
\end{figure}

Each iWarp node in PHASE5 executes a small program that implements its
flowgraph operator. Each twiddle node ($W^\prime$) repeatedly reads a
complex number from its input pathway, multiplies by the appropriate
twiddle (precomputed off-line using the recipes in
Section~\ref{sec:recipes}), 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:recipes}) 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.

While the details are beyond the scope of this paper, each node
program in the PHASE5 mapping performs each of its iterations in at
most 8 clocks. At its 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 2 32-bit floating point words. {\em As a result, each
pathway in the PHASE5 mapping requires 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! The computational throughput is 8 MFLOPS per node,
again independent of the size of the FFT.

A potential disadvantage of the PHASE5 mapping is that it is
node-inefficient, requiring 5 iWarp nodes for each FFT stage, plus an
additional 5 nodes for the bit reversal at the end. 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. Each node program in the
PHASE3 mapping executes in at most 11 clocks per iteration, for a
communication throughput of 30 Mbytes/sec, and a computational
throughput of 10 MFLOPS per node. So on the iWarp, the PHASE3 mapping
trades off higher computational throughput for lower communication
throughput.

\subsection{Example FFT on iWarp}
Figure~\ref{fig:example} shows a working implementation of a 2K-point
radix-2 phase-rotation FFT running on a 64-node iWarp array at
Carnegie Mellon. The large squares are iWarp nodes, labeled with the
corresponding operator and stage number. The bit reversal at the end
is treated as another stage. The small squares are queues.  The arrows
are iWarp pathways. The implementation uses a PHASE5 mapping.  Each of
the 11 FFT stages uses 5 nodes, with an additional 5 nodes for the
bit-reversing step at the end, for a total of 60 nodes.

The measured communication throughput of this FFT on the iWarp
hardware is 40 Mbytes/second (5 Msamples/sec). The average
computational throughput for each node is 8 MFLOPS, for an aggregate
computational throughput of $8 \cdot 60 = 480$ MFLOPS.
\begin{figure}
\centering
\makebox{
\psfig{figure=example.eps}
}
\caption{2K phase-rotation FFT on iWarp}
\label{fig:example}
\end{figure}

\section{Concluding remarks}

\section*{Acknowledgements}
Thanks to thank Tom Warfel and LeeAnn Tzeng for their help with
the iWarp implementation, and Doug Noll and Doug Smith for discussions
that led to the node-efficient mapping.

\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}

