\input epsf
\input{psadobe}
\makeatletter
\long\def\unmarkedfootnote#1{{\long\def\@makefntext##1{##1}\footnotetext{#1}
}}\makeatother
\newcommand{\comment}[1]{}
\documentstyle[times,11pt]{cmu-art}

\begin{document}
\bibliographystyle{acm}

\title{
\begin{flushleft}
\vspace{-0.9in}
\scriptsize {\em
Scientific Programming, v6, n1, Spring, 1997, pp. 59--72
}
\end{flushleft}
\vspace{0.4in}
\Large\bf 
Performance Issues in HPF Implementations of Sensor--Based Applications
}

\author{David R. O'Hallaron, Jon Webb, Jaspal Subhlok\\
\\
School of Computer Science\\
Carnegie Mellon University\\
Pittsburgh, PA 15213\\
}
\date{}
\thispagestyle{empty}
\maketitle

\unmarkedfootnote{
This research was sponsored in part by the Advanced Research Projects
Agency/CSTO under two contracts: one monitored by SPAWAR (contract
number N00039-93-C-0152), the other monitored by Hanscom Air Force
Base (contract number F19628-93-C-0171), in part by the Air Force
Office of Scientific Research under contract F49620--92--J--0131, in
part by the National Science Foundation under Grant ASC-9318163, and
in part by grants from the Intel Corporation.  Authors' email
addresses: droh@cs.cmu.edu, webb+@cmu.edu, jass@cs.cmu.edu.  }

\begin{abstract}
Applications that get their inputs from sensors are an important and
often overlooked application domain for High Performance Fortran
(HPF).  Such {\em sensor--based applications} typically perform
regular operations on dense arrays, and often have latency and
throughput requirements that can only be achieved with parallel
machines. This paper describes a study of sensor--based applications,
including the fast Fourier transform, synthetic aperture radar
imaging, narrowband tracking radar processing, multibaseline stereo
imaging, and medical magnetic resonance imaging. The applications are
written in a dialect of HPF developed at Carnegie Mellon, and are
compiled by the Fx compiler for the Intel Paragon.  The main results
of the study are (1) that it is possible to realize good performance
for realistic sensor--based applications written in HPF, and (2) the
performance of the applications is determined by the performance of
three core operations: independent loops (i.e., loops with no
dependences between iterations), reductions, and index permutations.
The paper discusses the implications for HPF implementations and
introduces some simple tests that implementors and users can use to
measure the efficiency of the loops, reductions, and index
permutations generated by an HPF compiler.
\end{abstract}

\section{Introduction}
\label{sec:intro}

There is an important class of computer applications that manipulate
inputs from the physical environment.  The inputs are continuously
collected by one or more sensors and then passed on to the computer,
where they are manipulated and interpreted.  The sensors are devices
like cameras, antennas, and microphones. The manipulation of the
sensor inputs is variously referred to as signal processing or image
processing, depending on the dimensionality of the inputs. We refer to
the entire class of applications as {\em sensor--based applications}
to emphasize this common quality of processing inputs from the natural
world.

Sensor--based applications have traditionally been found in military
domains like radar and sonar, and there is an increasing interest in
commercial domains, such as medical imaging, surveillance, and
real--world modeling.  For example, a real--world modeling application
could use a stereo algorithm to acquire depth information from
multiple cameras and then use the information to build realistic 3D
models of the environment.  The models could then be used for
applications like virtual 3D conferencing, building walk--throughs, or
experiencing a sporting event from the point of view of one of the
players.

Sensor--based applications are an interesting and often overlooked
application domain for a parallel language like HPF~\cite{hpf93}. The
computations, which typically consist of regular operations on dense
arrays, are naturally expressed in HPF. Furthermore, there are often
stringent latency and bandwidth requirements that demand parallel
processing.  For example, a stereo program that extracts depth
information from multiple cameras can process only a few frames per
second on a powerful RISC workstation, which is well below the
standard video rate of 30 frames per second.  If the results of a
sensor--based computation are used to control some process, then there
is also be a minimal latency that can be tolerated. For example, an
online medical imaging application that gathers and processes multiple
images might automatically adjust the scanner to compensate for
movement by the patient. The importance of minimizing latency, rather
than just maximizing throughput, is one of the key properties that
distinguishes sensor--based applications from batch--oriented
scientific computations~\cite{webb93bb}.

This paper describes the results of an empirical study of the
performance of a set of HPF sensor--based applications on a commercial
parallel computer. The results were obtained using a prototype
compiler, developed at Carnegie Mellon, for a dialect of HPF running
on an Intel Paragon. The paper makes several main points: First,
contrary to the fears of many in the HPF community, performance for
the HPF applications we studied is good.  Second, a few core
computational patterns (parallel DO loops, reductions, and index
permutations) dominate sensor--based applications.  HPF implementors
can realize great benefits by focusing on these patterns.  Third,
there are some simple tests that HPF programmers and developers can
use to evaluate the efficiency of the parallel DO loops, reductions,
and index permutations that are crucial to the effective execution
of sensor--based applications.  Fourth, since the data sets in
sensor--based applications are often fixed by properties of the
sensors, scalability is an important issue.  Finally, the same
patterns that appear in sensor--based computations also appear in
scientific applications.  In particular, we will briefly examine an Fx
regional air quality modeling code and an Fx earthquake ground motion
modeling based on the method of boundary elements.

In Section~\ref{sec:fx} we give a brief overview of the prototype HPF
compiler (the Fx compiler) that was used in the study.
Section~\ref{sec:sensorapps} describes the applications that we
implemented in Fx and their performance on the Intel Paragon.
Sections~\ref{sec:loops}, ~\ref{sec:reductions}, and
~\ref{sec:permutations} describe some key issues in generating
efficient code for HPF DO loops, reductions, and permutations, and
introduce some simple tests for measuring the efficiency of these
operations. Section~\ref{sec:scale} discusses the issue of scalability
in sensor--based applications. Finally, Section~\ref{sec:scientific}
shows how the same DO loops, reductions, and permutations that are
crucial to sensor--based also appear in scientific computations.

\section{Fx overview}
\label{sec:fx}
The Fx project was started in Fall 1991 with the goal of learning how
to generate efficient code for programs written in the emerging HPF
standard.\footnote{Although the first meeting of the HPF Forum was not
until January 1992, preproposals from Rice, Vienna, and ICASE were
already circulating during Summer 1991, so the general form of the HPF
programming model was already clear by Fall 1991.} The input language
is a dialect of subset HPF and consists of F77 with HPF data layout
statements, array assignment statements with support for general
CYCLIC(k) distributions in an arbitrary number of array
dimensions~\cite{stichnoth93a,stichnoth94}, an index permutation
intrinsic, and a parallel DO loop that is integrated with arbitrary
user--defined associative reduction operators~\cite{yang93}. Fx also
provides a mechanism for mixing task and data parallelism in the same
program~\cite{gross94a,subhlok93,subhlok94a}.  The initial target was
the Intel iWarp. Fx was later ported to the IBM SP/2, the Intel
Paragon, and workstation clusters.

Much of the early work on Fx was driven by the 2D fast Fourier
transform (FFT) and algorithms for global image processing. The system
was generating inefficient versions of the FFT by Spring 1992 and
efficient versions of the FFT by Fall 1992. By Spring 1993 task and
data parallelism were integrated into the compiler, and by Fall 1993
arbitrary user--defined reductions were integrated into the parallel
DO loop. The first significant Fx applications were written in Summer
1994: a multibaseline stereo program (from scratch) and a spotlight
synthetic aperture radar program (from F77 code supplied by Sandia
Labs). By Spring 1995 there were a number of additional Fx
applications, including a stereo front-end for a real--world modeling
application developed by the second author (from scratch), a magnetic
resonance image reconstruction program (from scratch), an earthquake
ground motion modeling program (from F77 code developed at Southern
California Earthquake Center), and a regional air quality modeling
program (from F77 code derived from the CIT airshed model).

\section{Sensor--based applications}
\label{sec:sensorapps}

Sensor--based applications typically have the high--level two--stage
pipelined structure shown in Figure~\ref{fig:sbcomp}.
\begin{figure}[htbp]
\centerline{\epsfbox{figs/sbcomp.eps}}
\caption{The structure of sensor--based applications.}
\label{fig:sbcomp}
\end{figure}
The front--end accepts a stream of inputs from one or more sensors and
manipulates these inputs into some desired form. The back--end
interprets the results of the front--end and either displays them or
initiates some action. For example, in a radar tracking application,
the front--end might transform input phase histories from an
antenna array into an image in the spatial domain, and the back--end
would manipulate this image to name, identify and track objects of
interest.

The front--end processing typically consists of numerous, regular,
data parallel operations on dense arrays, requires high MFLOPS rates,
and the operations performed are usually data--independent.
Computations such as the fast Fourier transform (FFT), convolution,
scaling, thresholding, data reduction, and histogramming are common
operations. The back--end processing is typically more dynamic,
irregular, and data--dependent, with real--time scheduling of
processes. In this paper, we are concerned with the front--end
processing, where HPF on a parallel system is most appropriate.  For
the remainder of the paper, when we refer to sensor--based
applications we are referring to the front--end. 

One of the nice qualities of sensor--based computing is that many
applications have similar computational patterns. The similarities
allow us to focus on a few small application kernels, with the
assurance that anything that we learn about compiling these small
programs will accrue benefits in larger, more realistic programs.  Two
examples that capture many of the key computational patterns, and were
of tremendous help in the development of the Fx compiler, are the 2D
FFT (FFT2) and the image histogram (HIST). The high--level parallel
structure of these computations is shown in
Figure~\ref{fig:exemplars}.

\begin{figure}[htbp]
\centerline{\epsfbox{figs/fft2.eps}}
\centerline{(a) FFT2 -- 2D fast Fourier transform}
\vspace*{2.0\baselineskip}
\centerline{\epsfbox{figs/hist.eps}}
\centerline{(b) HIST -- Image histogram}
\vspace*{\baselineskip}
\caption{Sensor--based computation exemplars}
\label{fig:exemplars}
\end{figure}

Figure~\ref{fig:exemplars} depicts the course--grained parallelism
that is available in FFT2 and HIST. The vertical lines depict
independent operations on array columns and the horizontal lines
depict independent operations on rows. The FFT2 program is collection
of independent local 1D FFTs on the columns of an array, followed by a
collection of independent local 1D FFTs on the rows.  If these columns
or rows are distributed across the processors, then each column or row
operation can run independently.  In HPF, the FFT2 example can be
written as:
\begin{verbatim}
    COMPLEX a(N,N),b(N,N)
    !HPF$ DISTRIBUTE (*,BLOCK):: a,b

    !HPF$ INDEPENDENT
    DO k=1,N
     call fft(a(:,k))
    ENDDO
    b = TRANSPOSE(a)
    !HPF$ INDEPENDENT
    DO k=1,N
     call fft(b(:,k))
    ENDDO
    a = TRANSPOSE(b)
\end{verbatim}
The Fortran 90 TRANSPOSE intrinsic is important because it allows each
1D FFT to operate locally and in place on a contiguous vector. The
same effect could also be achieved with the Fortran 90 RESHAPE
intrinsic, using the ORDER argument to swap the array
indices~\cite{fortran90}:
\begin{verbatim}
    b = RESHAPE(a, (/N,N/), (/0/), (/2,1/))	
\end{verbatim}
An alternative is to distribute one array by columns and the other
array by rows, and then redistribute from columns to rows using an
array assignment statement. This approach results in rows that are
stored locally but not contiguously. Since it would require extra
copies to get the rows in contiguous form, this approach will not be
considered here.

The HIST example consists of a collection of independent local
histograms on the columns of an integer array, followed by a
plus--reduction operation that adds the local histogram vectors to
form the final result. This might be written in HPF as:
\begin{verbatim}
    INTEGER a(N,N),h(M,N),r(M)
    !HPF$ DISTRIBUTE (*,BLOCK):: a

    h = 0
    !HPF$ INDEPENDENT
    DO j=1,N
     DO i=1,N
      h(i,j) = h(a(i,j),j) + 1
     ENDDO
    ENDDO
    r = SUM(h,2)
\end{verbatim}
The FFT2 and HIST examples capture the core computational patterns in
sensor--based computing: parallel DO loops, reductions, and index
permutations. FFT2 is a pair of parallel DO loops, each of which is
followed by an index permutation (the TRANSPOSE intrinsic).  HIST is a
parallel DO loop followed by a reduction. These two patterns occur
repeatedly in the sensor--based applications we have studied.

Figure~\ref{fig:sensorapps} shows a collection of sensor--based
applications. All but ABI (Figure~\ref{fig:sensorapps}(b)) have been
implemented in Fx, and could be ported to HPF with small changes.  The
STEREO program, developed at Carnegie Mellon, extracts depth
information using the images from multiple video
cameras~\cite{kang95}. The RADAR program was adapted from a C program
developed by MIT Lincoln Labs to measure the effectiveness of various
multicomputers for their radar applications~\cite{shaw92}.  The SAR
program was adapted from a Fortran 77 program developed by Sandia
National Laboratories\cite{plimpton91}. The MR program was developed
from an algorithm by Doug Noll at Pitt Medical Center~\cite{noll92}.

\begin{figure}[htbp]
\centerline{\epsfbox{figs/mstereo.eps}}
\centerline{(a) STEREO -- Multibaseline stereo.}
\vspace*{2.0\baselineskip}
\centerline{\epsfbox{figs/abi.eps}}
\centerline{(b) ABI -- Sonar adaptive beam interpolation.}
\vspace*{2.0\baselineskip}
\centerline{\epsfbox{figs/radar.eps}}
\centerline{(c) RADAR -- Narrowband tracking radar.}
\vspace*{2.0\baselineskip}
\centerline{\epsfbox{figs/sar.eps}}
\centerline{(d) SAR -- Synthetic aperture radar image reconstruction.}
\vspace*{2.0\baselineskip}
\centerline{\epsfbox{figs/mri.eps}}
\centerline{(e) MR -- Magnetic resonance image reconstruction.}
\vspace*{2.0\baselineskip}

\caption{Other Fx/Paragon sensor--based applications}
\label{fig:sensorapps}
\end{figure}

A striking aspect of Figure~\ref{fig:sensorapps} is the number of
parallel DO loops that operate independently along one dimension or
another of the array.  Each application contains at least one of these
loops.  The pointwise scaling operation in RADAR is also another form
of parallel DO loop, which is usually expressed as an array
assignment.  Another common pattern from the FFT2 example is: (1)
operate along one dimension, then (2) operate along another dimension.
This pattern, which occurs in FFT2, ABI, RADAR, SAR, and MR, is
typically implemented with a TRANSPOSE intrinsic between (1) and (2).
Reductions are found in HIST, STEREO, ABI, and RADAR. The point is
that FFT2 and HIST capture the basic computational structure of a wide
range of sensor--based computations.

There has been some concern about the performance that can be expected
from HPF programs~\cite{knies94}. However, in our experience, the
performance of HPF--like programs compiled by the Fx compiler for the
Paragon is good, even for moderately sized problems.
Figure~\ref{fig:perf} shows the absolute performance of
representatively sized FFT, HIST, and SAR programs.  FFT1 is a
parallel 1D FFT program, FFT2 is the FFT exemplar, and FFT3 is a
parallel 3D FFT program; each is computed in a way similar to FFT.
The programs in Figure~\ref{fig:perf} scale reasonably well (although
not linearly) and running time is not dominated by communication
overhead. In the case of the SAR program, communication accounts for
less 10\% of the running time.
\begin{figure}[htbp]
\begin{minipage}[t]{3in}
\centerline{\epsfxsize=2.75in \epsfbox{figs/fftc1perf.eps}}
\centerline{(a) $256K \times 1$ 1D FFT.}
\vspace*{2.0\baselineskip}
\end{minipage} \
\begin{minipage}[t]{3in}
\centerline{\epsfxsize=2.75in \epsfbox{figs/fftc2perf.eps}}
\centerline{(b) $512 \times 512$ 2D FFT.}
\vspace*{2.0\baselineskip}
\end{minipage}

\begin{minipage}[t]{3in}
\centerline{\epsfxsize=2.75in \epsfbox{figs/fftc3perf.eps}}
\centerline{(c) $64 \times 64 \times 64$ 3D FFT.}
\vspace*{2.0\baselineskip}
\end{minipage}\
\begin{minipage}[t]{3in}
\centerline{\epsfxsize=2.75in \epsfbox{figs/histperf.eps}}
\centerline{(d) $1K \times 1K$ image histogram.}
\vspace*{2.0\baselineskip}
\end{minipage}

\begin{minipage}[t]{3in}
\centerline{\epsfxsize=2.75in \epsfbox{figs/sarperf.eps}}
\centerline{(e) $512 \times 512$ synthetic aperture radar}
\end{minipage}

\caption{Fx/Paragon performance.}
\label{fig:perf}
\end{figure}

While certainly not exhaustive, Figure~\ref{fig:perf} offers some hope
that good performance can be expected from sensor--based applications
written in HPF.  In the remaining sections, we will discuss the issues
involved in achieving good performance.

\section{Parallel loops}
\label{sec:loops}

DO loops in typical sensor--based applications can be efficiently
parallelized using a variation on the simple FORTRAN D copy--in
copy--out model~\cite{fortrand92}. The computation in the main body of
the program is modeled as a single thread operating on a global data
space. Each iteration in a loop is modeled as a separate thread
operating on its own local data space. When control reaches the loop,
the contents of the global data space are (conceptually) copied in to
each of the local data spaces.  Each loop iteration then works
independently on its local copy of the global data space. When all of
the loop iterations have terminated, the contents of the local data
spaces are (conceptually) copied out of the local data spaces back
into the global data space. If multiple iterations write to the same
address in the local address space, then the values are merged with a
user--defined binary associative reduction operator before copying
back to the global address space.  This loop model, called the {\em
PDO model}, is described in detail in~\cite{yang93}.

A parallel DO loop based on the PDO model can be characterized in
terms of the addresses that it references. Let $R_k$ denote the set of
addresses read by iteration $k$ and $W_k$ the set of addresses written
by iteration $k$.  If each $R_k$ is disjoint, then the loop has {\em
disjoint reads}, otherwise the loop has {\em overlapped reads}.
Similarly, if each $W_k$ is disjoint, then the loop has {\em disjoint
writes}, otherwise the loop has {\em overlapped writes}.

The most common form of loop in sensor--based applications has
disjoint reads and writes.  Every application in
Figures~\ref{fig:exemplars} and \ref{fig:sensorapps} has at least one
loop with disjoint reads and writes, and FFT2
(Figure~\ref{fig:exemplars}(a)), SAR (Figure~\ref{fig:sensorapps}(d)),
and MR (Figure~\ref{fig:sensorapps}(e)) consist exclusively of these
kinds of loops. For example, the FFT2 exemplar consists of HPF DO
loops of the form:

\begin{verbatim}
    !HPF$ INDEPENDENT
    DO k=1,N
     CALL fft(a(:,k))
    ENDDO
\end{verbatim}

Each invocation of the \verb.fft(). subroutine performs an inplace 1D
FFT on the $kth$ column of an array, reading and writing only elements
in the $kth$ column.  Although \verb.fft(). is a complicated
subroutine with a complex pattern of array references (and might even
be an assembly language library routine), the pattern of array
references between loop iteration is extremely simple: the $kth$
iteration references the $kth$ column.  This dichotomy of complicated
intra--iteration reference patterns and simple inter--iteration
reference patterns is a recurring theme in sensor--based computing,
with important implications for HPF implementations.

Another important form of DO loop has overlapped reads and disjoint
writes. Loops of this form are typically used to perform convolution
operations such as the error computation in STEREO
(Figure~\ref{fig:sensorapps}(a)).  A similar pattern occurs in
relaxation algorithms from scientific computing. For example, a simple
1D convolution is of the form:
\begin{verbatim}
    REAL a(N),b(N),h(3)
    !HPF$ INDEPENDENT
    DO k=2,N-1
      b(k) = a(k-1)*h(1) + a(k)*h(2) + a(k+1)*h(3)
    ENDDO
\end{verbatim}
Finally, loops with overlapped writes are typically used by
sensor--based applications to implement reductions. We discuss this
important class of loops in Section~\ref{sec:reductions}. The
remainder of this section discusses only loops with disjoint writes.

\subsection{Implications for HPF implementations}

Generating efficient code inside parallel loops is the key to
achieving good performance with sensor--based applications written in
HPF. And since parallel loops with disjoint writes are so common,
occurring in every application we have studied, generating efficient
code for these loops is especially important.

Although loops with disjoint writes are often dismissed as
``embarrassingly parallel'', it is nontrivial in general to generate
efficient parallel code for them. There are a number of reasons, all
complicated by the fact that loop bodies of real applications
typically contain a lot of code, with complex intra--iteration
reference patterns, calls to external library routines, and even
inlined assembly language inserts.

First, the HPF compiler must somehow determine that the write sets are
disjoint and that addresses that are written by one iteration are not
read by another iteration.  The HPF INDEPENDENT directive is very
helpful here.  The INDEPENDENT directive informs the compiler that no
address is written by one iteration of a DO loop and then read or
written by another iteration.

Second, the HPF compiler must ensure that the read and write sets are
aligned with the loop iterations before the iterations are executed.
By {\em aligned}, we mean that the read and write sets of each loop
are available in local memory.  Aligning the data sets before
executing a loop is key to achieving good performance because this
allows the programmer to use arbitrary sequential code in the loop
body, including calls to efficient sequential math libraries written
in assembly language.  For example, in a Paragon HPF implementation,
the \verb.fft(). routine called by the FFT2 loop might be an assembly
language routine hand--crafted for the i860 microprocessor.

Finally, the HPF compiler must compute local loop bounds and translate
global array indices in the loop body to local indices. If not handled
properly, these computations can be a significant source of runtime
overhead.

The Fx compiler relies on a new PDO keyword to assert that write sets
are disjoint and that addresses that are written by one iteration are
not read by another. (The HPF INDEPENDENT directive conveys the same
information and is more compatible with standard F90 compilers.)
Also, since the bodies of loops with disjoint writes can be
arbitrarily complex, consisting of calls to externally compiled
library routines or hand--coded assembly routines, it is not always
possible to rely on compile--time analysis to align read and write
sets with loop iterations.  Thus, the Fx compiler relies on hints from
the user that describe the read and write sets for each loop
iteration. This is an important point because it allows Fx to generate
efficient code for loops with disjoint writes, regardless of the
complexity of the loop body.

\subsection{Loop efficiency test}
\label{sec:loopeff}
There is a simple test, called the {\em loop efficiency test}, that
implementers and users can use to measure the overhead introduced by
HPF compilers in the loop bodies of parallel loops.  Consider the
following canonical parallel DO loop:
\begin{verbatim}
    REAL a(N,N)
    !HPF$ DISTRIBUTE (*,BLOCK):: a

    !HPF$ INDEPENDENT
    DO k=1,N
      a(:,k) = k
    ENDDO
\end{verbatim}
This loop requires no communication at runtime, but is somewhat subtle
to translate because the {\em lhs} instance of $k$ must be converted
from a global index to a local index, but the {\em rhs} instance must
remain a global index.  If we compile and run the loop on $P$ nodes,
where $P$ divides $N$ evenly, then the total running time is bounded
from below by the running time of the following sequential DO loop:
\begin{verbatim}
    REAL sa(N,N/P)

    DO k=1,N/P
      sa(:,k) = k
    ENDDO
\end{verbatim}
If $L_{p}(N,P)$ is the running time of the canonical parallel DO loop
with $N$ iterations and $P$ nodes, and $L_{s}(N,P)$ the running time
of the corresponding sequential DO loop with $N/P$ iterations, then
$E_{loop}(N,P) = L_{s}(N,P)/ L_{p}(N,P)$ is the {\em loop efficiency}
of the parallel DO loop generated by the HPF compiler.  Loop
efficiency is a useful measure because it provides a way to isolate
the runtime loop overheads that are introduced by the compiler,
without having to instrument the generated code. It is important to
realize that loop efficiency is not the same as the speedup over the
single--node version of the parallel code. Instead, we are comparing
to the performance of the tightest {\em sequential} version of the
loop.

\begin{figure}[htbp]
\epsfxsize=4in
\centerline{\epsfbox{figs/loopeff2.eps}}
\caption{Fx/Paragon loop efficiency ($N \times N$ canonical loop).}
\label{fig:loopeff}
\end{figure}

If we record loop efficiency for different values of $N$ and $P$, we
get an interesting family of curves.  Figure~\ref{fig:loopeff} shows
the results for the Fx version of the canonical loop on the Intel
Paragon.

The family of curves in Figure~\ref{fig:loopeff} provides some
intestesting insight into the quality of the parallel loops generated
by the compiler. Loop efficiency is bounded from above by the curve
for 1 node and bounded from below by the curve for 32 nodes, so in
general it is only necessary to plot two curves.  In
Figure~\ref{fig:loopeff} the loop efficiency for 1 node is almost
unity. For a single node the compiler introduces almost no overhead,
which tells us that the loop is nearly as tight as the corresponding
sequential loop. The curve for 32 nodes converges to near unity, which
tells us that the overheads are being amortized across loop
iterations. Further, the 32--node curve converges quite rapidly, with
$E_{loop} > 0.5$ at $N=256$ and $E_{loop} > 0.8$ at $N=512$. Thus the
Fx compiler is introducing minimal overheads that are quickly
amortized.  This conclusion is confirmed by inspecting the F77 code
generated by Fx for $N=1024$ and $P=8$:
\begin{verbatim}
    IF (fxcellid.LT.8) THEN
     fxloopstart0 = (MAX(((fxcellid * 128) + 1),1))
     fxlmidx0 = IFXLM(fxadesc,1,fxloopstart0)
     DO k = fxloopstart0, MIN(((fxcellid * 128) + 128),1024), 1
      DO fxindex1 = 1, 1024, 1
       a(fxindex1,fxlmidx0) = (k)
      ENDDO
      fxlmidx0 = fxlmidx0 + 1
     ENDDO
    ENDIF
\end{verbatim}
The parallel loop overhead consists of a few statements before the
loop that compute the local loop bounds, and a function call that
computes the initial local index value. The only overhead in the loop
body is a statement that increments the local index value.  A similar
approach to index conversion is first described
in~\cite{chatterjee93}.  It is hard to imagine a tighter loop.

In summary, the loops in sensor--based applications can often be
implemented as independent parallel loops, without any communication
between loop iterations. A primary goal of an HPF implementation
should be to generate parallel loop bodies that are as efficient as
their sequential counterparts. In particular, implementors should
focus on minimizing the overhead of parallel loops with disjoint reads
and writes.  The loop efficiency test we have introduced in this
sections provides a simple way to characterize these overheads.

\section{Reductions}
\label{sec:reductions}

In sensor--based applications, loops with overlapped writes are used
primarily to implement reductions. For example,
\begin{verbatim}
    DO k=1,N
     v = v + a(:,k)
    ENDDO
\end{verbatim}
A common pattern in sensor--based applications is to operate
independently on the columns of an array, and then reduce the columns
into a single column by adding them together. The HIST, STEREO, ABI,
and RADAR program all perform this type of simple reduction.  However
there are other sensor--based applications that require a mechanism
for the programmer to define generalized reduction operations. For
example, a connected components algorithm can be written as a parallel
loop over the rows of the image, where each iteration computes a
segment table for its row. This is followed by a generalized reduction
step that merges the segment tables.

\subsection{Implications for HPF implementations}

For most sensor--based applications the HPF SUM intrinsic is
sufficient. However, the first version of HPF provides no support for
operations like connected components that require generalized
reductions.  Achieving good performance in these cases will require
sophisticated compiler analysis to recognize the
reductions~\cite{ghuloum95}.  In Fx, we avoid this analysis by
incorporating a mechanism for defining arbitrary binary associative
reductions into the parallel loop construct~\cite{yang93}.

\subsection{Reduction efficiency test}
\label{sec:redeff}
A user or implementer can measure the quality of the parallel
reduction loops generated by an HPF compiler using a test similar to
the loop efficiency test in Section~\ref{sec:loopeff}. Consider the
following loop that adds the columns of an $N \times N$ array.
\begin{verbatim}
    REAL a(N,N),v(N)
    !HPF$ DISTRIBUTE (*,BLOCK):: a

    v = 0.0
    DO k=1,N
     v = v + a(:,k)
    ENDDO
\end{verbatim}
The performance of this loop on $P$ nodes is bounded from below
by the performance of the following sequential HPF DO loop:
\begin{verbatim}
    REAL sa(N,N/P),v(N)

    v = 0.0
    DO k=1,N/P
     v = v + sa(:,k)
    ENDDO
\end{verbatim}
If $R_{p}(N,P)$ is running time of the parallel reduction of an $N
\times N$ array on $P$ nodes and $R_{s}(N,P)$ the running time of the
corresponding sequential reduction of an $N \times N/P$ array, then
$E_{reduce}(N,P) = R_{s}(N,P)/R_{p}(N,P)$ is the {\em reduction
efficiency} of the parallel reduction generated by the
compiler.\footnote{An alternative formulation of the reduction
efficiency test is to use the HPF SUM intrinsic for the parallel
reduction.  In this case, the sequential reduction must use the same
local computational kernel as the SUM intrinsic.} Reduction efficiency
exposes the runtime overheads that are incurred by performing the
reduction in parallel.  Unlike loop efficiency, which is completely
determined by the compiler, reduction efficiency is a function of
overheads due to the compiler, as well as overheads due to the
underlying communication system. Thus, in general it may be impossible
to achieve a reduction efficiency of unity.

\begin{figure}[htbp]
\epsfxsize=4in
\centerline{\epsfbox{figs/redeff.eps}}
\caption{Fx/Paragon reduction efficiency ($N \times N \rightarrow N$
plus-reduction).}
\label{fig:redeff}
\end{figure}

Figure~\ref{fig:redeff} shows the results for a simple plus--reduction
using the Fx compiler on Paragon. As with the loop efficiency graph,
reduction efficiency is bounded from above by the curve for 1 node and
from below by the curve for 32 nodes. Surprisingly, the reduction
efficiency for a single node decreases as the problem size increases.
This suggests a problem in the Fx implementation of the reduction.
Ideally, the efficiency on a single node should be close to unity.
Another point of concern is the slow convergence of the curves for 16
and 32 nodes. Since the local computation step of each parallel
reduction grows as roughly $N^2/P$ and the communication step grows as
roughly $N\log P$, we might expect these curves to converge faster
than they do. Yet even for a relatively large $N$, the reduction
efficiency is below 50\%.  While the reduction efficiency does not
pinpoint the source of the overhead, it does point out an opportunity
for improvement in the Fx implementation.

\section{Index permutations}
\label{sec:permutations}

As we saw in Section~\ref{sec:sensorapps}, the following computational
pattern occurs in many sensor--based applications: (1) operate
independently along one dimension of an array, then (2) operate
independently along another dimension. The FFT2, ABI, RADAR, SAR, and
MR programs all exhibit this pattern. For example, FFT2 performs a
local FFT on each column of an array, then performs a local FFT on
each row.  In order to exploit locality, this pattern is usually
implemented with an {\em index permutation} (also referred to as a
transpose or corner turn) between steps (1) and (2):
\begin{verbatim}
    DO k=1,N
     call fft(a(:,k))
    ENDDO
    b = TRANSPOSE(a)
    DO k=1,N
     call fft(b(:,k))
    ENDDO
    a = TRANSPOSE(b)
\end{verbatim}
In general, an index permutation is an assignment from a
$d$--dimensional array $a$ to a $d$--dimensional array $b$ such that
after the assignment,
\begin{displaymath}
b(i_1, i_2, \ldots , i_d) = a(\pi (i_1, i_2, \ldots , i_d)), 
~~~~~1 \leq i_1, i_2 ,\ldots, i_d \leq N
\end{displaymath}
where each $i_k$ is a unique array index and $\pi$ is some permutation
operator. For example, after the TRANSPOSE step in the FFT2 program,
$b(i,j) = a(j,i)$, for $1 \leq i,j \leq N$. 

Although most of the sensor--based applications that we have studied
transpose 2D arrays, there are important cases where index
permutations of higher--dimensional matrices are necessary. In
particular, 2D arrays of complex variables are often implemented as 3D
arrays of real variables, and for $d > 1$, a $d$--dimensional FFT must
permute two indices of a $d$--dimensional complex array between each
local FFT step.

\subsection{Implications for HPF implementations}

Efficient index permutation is crucial to achieving good performance
in sensor--based applications. In general, an index permutation
induces a complete exchange (i.e., an all--to--all personalized
communication), where each node sends data to every other node.  The
Fortran 90 TRANSPOSE and RESHAPE intrinsics~\cite{fortran90} adopted
by HPF provide an opportunity to optimize this important operation,
but unfortunately TRANSPOSE is only defined for 2D arrays.  So for the
general case, an HPF implementation must either provide a highly tuned
RESHAPE intrinsic, or provide a generalized TRANSPOSE extrinsic (this
is being considered for HPF-2), or be able to generate efficient code
for index permutations that are implemented with a combination of
array assignments and DO loops:
\begin{verbatim}
    REAL a(N,N),b(N,N)
    !HPF$ DISTRIBUTE (*,BLOCK):: a
    !HPF$ DISTRIBUTE (BLOCK,*):: b

    b = a
    !HPF$ INDEPENDENT
    DO k=1,N
     a(:,k) = b(k,:)
    ENDDO
\end{verbatim}
The Fx compiler provides an index permutation intrinsic with the same
basic functionality as a FORTRAN 90 RESHAPE intrinsic with an ORDER
argument.  The advantage of using an intrinsic is that an intrinsic
can leverage off of the existing code in the compiler for generating
array assignment statements. Writing an extrinsic with the same
functionality means duplicating the compiler's array assignment code
in the run--time library.  Furthermore, capturing the index
permutation in an intrinsic allows the compiler to exploit significant
optimizations on systems with toroidal
interconnects~\cite{hinrichs94}.

\subsection{Permutation efficiency test}

Just as with loops and reductions, there is a simple test for
measuring the efficiency of HPF index permutations. The execution time
of a parallel index permutation of an $N \times N$ array (using either
the TRANSPOSE intrinsic, a permutation extrinsic, or an assignment
statement and a DO loop) is bounded from below by the time to
sequentially permute an $N \times N/P$ array on a single node:
\begin{verbatim}
    REAL sa(N/P,N),sb(N,N/P)

    DO k=1,N/P
     sa(k,:) = sb(:,k)
    ENDDO
\end{verbatim}
If $T_{p}(N,P)$ is the running time of the parallel index permutation
of an $N \times N$ array and $T_{s}(N,P)$ is the running time of the
corresponding sequential permutation of an $N \times N/P$ array, then
$E_{permute}(N,P) = T_{s}(N,P)/T_{p}(N,P)$ is the {\em permutation
efficiency} of the parallel index permutation generated by the HPF
compiler.\footnote{As with the reduction efficiency test, if an
intrinsic or extrinsic is used for the parallel permutation, then care
must be taken to use the same local copy mechanism in the sequential
version.} Permutation efficiency is a rough measure of the percentage
of effective local memory bandwidth that is realized by the parallel
permutation. Like reduction efficiency, permutation efficiency is
influenced by overheads due to the compiler, as well as overheads due
to the underlying communication system.

\begin{figure}[htbp]
\epsfxsize=4in
\centerline{\epsfbox{figs/permeff.eps}}
\caption{Fx/Paragon permutation efficiency ($N \times N$ 2D transpose).}
\label{fig:permeff}
\end{figure}

Figure~\ref{fig:permeff} shows the results for a 2D transpose of an $N
\times N$ array using the Fx compiler on Paragon. The graph provides
a couple of interesting insights. There is substantial overhead even
for the single--node version of the parallel permutation, which achieves
only 85\% of the effective local memory bandwidth for large problems.
The inefficiency is due to the fact that the Fx--generated code
unnecessarily checks for communications that never occur (because the
transpose is on a single node).

The multiple--node versions of the parallel permutation use only 30\%
of the effective local memory bandwidth for large problem sizes. This
suggests that the parallel permutation on the Paragon is
communication--bound, and that further improvements will require a new
message passing layer.  The inefficiency is largely due to overhead
from the underlying communication system and it is tempting for us to
wash our hands of responsibility for its performance.  However, in our
experience, significant performance benefits can be realized in
compiler--generated code by tailoring the runtime communication
libraries~\cite{stichnoth94,ISCA-MEMORY}.  HPF implementors need to be
aware of the communication overheads for a particular target machine.
The reduction and permutation efficiency tests are a simple and
useful way to expose the performance impact of these overheads.

\section{Scalability}
\label{sec:scale}

Sensor--based applications are typically implemented as 
collections of functions
that process continuous streams of data sets. The sizes of the data
sets are determined by external factors such as the type of sensor,
the number of sensors, and the frequencies of interest. For example,
the image size of the STEREO application is fixed at $240 \times 256$
by the camera system and cannot be modified by the programmer. The
magnetic resonance scanner used by the MR application processes $512
\times 512$ images (oversampled from $256 \times 256$ input). The
radar subsystem used by the RADAR application produces $512 \times
10$ data sets.

The fixed size of the data sets is an important property of
sensor--based applications that distinguishes them from typical
scientific applications.  Since the data set sizes are fixed, the
amount of available computational work on each node decreases as the
number of nodes increases. If a data parallel function performs a
nontrivial amount of internal communication, then the efficiency of
the function will tend to decrease as the number of nodes increases.
This behavior is shown in Figure~\ref{sec:scale} for a $512 \times
512$ local FFT loop, a $512 \times 512$ image histogram, and a $512
\times 512$ transpose. The local FFT function contains no
communication, and thus scales perfectly with the number of nodes.
However, the histogram and transpose functions contain internal
communication and their efficiency decreases significantly as the
number of nodes increases.

\begin{figure}[htbp]
\epsfxsize=4in
\centerline{\epsfbox{figs/scale.eps}}
\caption{Scalability of various Fx/Paragon functions}
\label{fig:scale}
\end{figure}

If efficient use of processing nodes is a goal (as it is in embedded
systems where additional nodes increase the cost, size, power, and
weight of the system) then we want to use a smaller number of nodes
for functions like the histogram and transpose. But if we have a large
parallel system with many nodes, how can we effectively use the
remaining nodes?  One approach that has been proposed is to use a mix
of task and data
parallelism~\cite{subhlok93,CFKK94,chapman94,gross94a}.

Task parallelism can significantly improve the performance of
applications with functions that do not scale well. For example, using
a mix of task and data parallelism doubled the throughput (compared to
the most efficient data parallel code) of the $240 \times 256$ Fx
STEREO program so that it was able to run in
real--time~\cite{subhlok94a}. Since HPF does not currently support
task parallelism (although it is being considered by the HPF--2
committee) there is the risk that HPF sensor--based codes with smaller
data sets will not run efficiently.  This puts additional pressure on
HPF developers to maximize the loop, reduction, and permutation
efficiencies identified earlier.

\section{Relation to scientific codes}
\label{sec:scientific}

Although this paper does not specifically address scientific HPF
codes, other groups have used Fx to implement two nontrivial physical
simulations: QUAKE, a 3D earthquake ground motion simulation (based on
the method of boundary elements) and AIR, a regional air quality
modeling program~\cite{kumar95}. Both codes are legacy F77 codes of
about $10,000$ lines that were ported to the Paragon version of Fx.
QUAKE is especially interesting because it was ported in a few weeks
by a seismologist from the Southern California Earthquake Center who
had never written a parallel program.

Figure~\ref{fig:scicodes} shows the computational structure of the
QUAKE and AIR programs.  QUAKE is a single perfectly parallel DO
loop. AIR is a sequence of DO loops that operate on different
dimensions of an array, with 3D transposes interspersed between the
loops. The interesting thing is that at a high level, the structure of
these two moderately large scientific codes is almost identical to the
FFT!  So again, we see that parallel loops and index permutations are
important operations for HPF developers to optimize.

\begin{figure}[htbp]
\centerline{\epsfbox{figs/bemquake.eps}}
\centerline{(a) QUAKE -- Boundary element earthquake ground motion modeling.}
\vspace*{2.0\baselineskip}
\centerline{\epsfbox{figs/airshed.eps}}
\centerline{(b) AIR -- Finite element air quality modeling.}
\vspace*{\baselineskip}
\caption{Complex scientific codes with simple data parallel structure}
\label{fig:scicodes}
\end{figure}

The QUAKE and AIR programs reinforce an important point that we
touched on in Section~\ref{sec:loops}: complicated programs with
complicated inner loops can nonetheless have a simple data parallel
structure that is straightforward to parallelize. The AIR program
takes this to extremes: each iteration of the parallel DO loop in each
of the horizontal transport steps solves an independent {\em sparse
and irregular finite element problem}. We normally assume that HPF is
not a good platform for sparse irregular codes, but AIR is an example
of a sparse code that is quite well suited for HPF.

\section{Conclusions}
\label{sec:conclusions}

We have identified sensor--based applications as an important class of
applications that are generally well suited to HPF. Sensor--based
applications operate on dense arrays of data and they reference this
data in regular ways that are naturally expressed using HPF loops and
array assignment statements. Typical sensor--based applications
include fast Fourier transforms, radar and sonar processing, computer
vision, and medical imaging.

In the process of studying realistic sensor--based applications
compiled by a subset HPF compiler developed at Carnegie Mellon, we
observed that performance is determined by the efficiency of three key
operations: independent DO loops, reductions, and index permutations.
Of course, good performance on these operations is important for many
types of applications.  The interesting point is that for
sensor--based applications, good performance for a small set of
operations is generally sufficient for achieving good overall
performance.

Implementing an efficient compiler for all of HPF is difficult and
expensive. We believe that focusing attention and effort on the
performance of independent DO loops, reductions, and index
permutations will allow for the rapid development of practically
efficient compilers.  Based on our experience with the Fx system, HPF
developers who focus on these operations will reap large rewards in
performance.

Given the importance of loops, reductions, and index permutations for
sensor--based applications, we have introduced some simple and useful
tests that users and implementers can use to measure the efficiency of
the code generated by an HPF compiler for these key operations.  The
tests are written entirely at the user level, but can provide good
insight into the quality of code generated by an HPF compiler. For
example, the Fx compiler, as well as several commercial HPF compilers,
exhibit a measured loop efficiency of close to 100\% for the canonical
DO loop. In other words, the loop bodies generated by these compilers
are nearly as good as their sequential counterparts.  On the other
hand, several commercial HPF compilers exhibit a maximum loop
efficiency of only 50\% for the same DO loop, which suggests a problem
in the way these compilers generate loops.  A future paper will
address the performance of commericial HPF compilers in more detail.

Finally, although HPF is a good match in general with sensor--based
applications, scalability can be an issue because of the fixed sizes
of the data sets. If the data sets are small, then a mix of task and
data parallelism is often required to get good performance.

\section*{Acknowledgements}
Keith Bromley at the Naval Oceans Systems Center encouraged us to
search for similarities in signal and image processing applications.
Dennis Ghiglia at Sandia Labs generously provided us with F77 SAR
code, which Peter Lieu ported to Fx. Jim Wheeler at GE taught us about
underwater sonar applications.  Doug Noll at Pitt Medical Center
developed the MR algorithm, and Claudson Bornstein, Bwolen Yang, and
Peter Lieu implemented it in Fx.  Yoshi Hisada from the USC Southern
California Earthquake Center took a chance and implemented his 3D
boundary element ground motion algorithm in Fx.  Ed Segall, Chang-Hsin
Chang, and Peter Lieu ported the air quality modeling application to
Fx. Thomas Gross, Jim Stichnoth, Bwolen Yang, and Peter Dinda made
major contributions to the Fx compiler. Thanks also to the anonymous
referees.

\bibliography{/afs/cs/project/iwarp/member/droh/bib/refs}
\end{document}

%.bbl file
\begin{thebibliography}{10}

\bibitem{fortran90}
{\sc {American National Standards Institute, Inc.}}
\newblock {\em ANSI Systems Programming Language Fortran(Fortran 90),Draft S8,
  Version 114 (X3.9-1990)}, Jan. 1990.

\bibitem{CFKK94}
{\sc Chandy, M., Foster, I., Kennedy, K., Koelbel, C., and Tseng, C.}
\newblock Integrated support for task and data parallelism.
\newblock {\em International Journal of Supercomputer Applications 8}, 2
  (1994), 80--98.

\bibitem{chapman94}
{\sc Chapman, B., Mehrotra, P., {Van Rosendale}, J., and Zima, H.}
\newblock A software architecture for multidisciplinary applications:
  Integrating task and data parallelism.
\newblock Tech. Rep. 94-18, ICASE, NASA Langley Research Center, Hampton, VA,
  Mar. 1994.

\bibitem{chatterjee93}
{\sc Chatterjee, S., Gilbert, J., Long, F., Schreiber, R., and Teng, S.}
\newblock Generating local addresses and communication sets for data-parallel
  programs.
\newblock In {\em Proceedings of the ACM SIGPLAN Symposium on Principles and
  Practice of Parallel Programming\/} (San Diego, CA, May 1993), pp.~149--158.

\bibitem{ghuloum95}
{\sc Ghuloum, A., and Fisher, A.}
\newblock Flattening and parallelizing irregular, recurrent loop nests.
\newblock In {\em Proc. of the Fifth ACM SIGPLAN Symposium on Principles and
  Practice of Parallel Programming (PPoPP)\/} (Santa Barbara, CA, July 1995).

\bibitem{gross94a}
{\sc Gross, T., O'Hallaron, D., and Subhlok, J.}
\newblock Task parallelism in a {H}igh {P}erformance {F}ortran framework.
\newblock {\em IEEE Parallel \& Distributed Technology 2}, 3 (1994), 16--26.

\bibitem{hpf93}
{\sc {High Performance Fortran Forum}}.
\newblock {H}igh {P}erformance {F}ortran language specification, version 1.0.
\newblock Tech. Rep. CRPC-TR92225, Center for Research on Parallel Computation,
  Rice University, May 1993.

\bibitem{hinrichs94}
{\sc Hinrichs, S., Kosak, C., O'Hallaron, D., Stricker, T., and Take, R.}
\newblock An architecture for optimal all-to-all personalized communication.
\newblock In {\em Proc. SPAA '94\/} (Cape May, NJ, June 1994), ACM,
  pp.~310--319.

\bibitem{fortrand92}
{\sc Hiranandani, S., Kennedy, K., and Tseng, C.~W.}
\newblock Compiling {F}ortran {D} for {MIMD} distributed-memory machines.
\newblock {\em CACM 35}, 8 (Aug 1992), 66--80.

\bibitem{kang95}
{\sc Kang, S., Webb, J., Zitnick, C., and Kanade, T.}
\newblock A multibaseline stereo system with active illumination and real--time
  image acquisition.
\newblock In {\em Proceedings of the International Conference on Computer
  Vision\/} (Cambridge, MA, 1995).

\bibitem{knies94}
{\sc Knies, A., O'Keefe, M., and MacDonald, T.}
\newblock {H}igh {P}erformance {F}ortran: {A} practical analysis.
\newblock {\em Scientific Programming 3}, 3 (Fall 1994), 187--200.

\bibitem{kumar95}
{\sc Kumar, N., Russel, A., Segall, E., and Steenkiste, P.}
\newblock Parallel and distributed application of an urban regional multiscale
  model.
\newblock Carnegie Mellon Dept. of Mech. Eng. and Dept. of Computer Science,
  1995.
\newblock working paper.

\bibitem{noll92}
{\sc Noll, D., Pauly, J., Meyer, C., Nishimura, D., and Macovski, A.}
\newblock Deblurring for non-2{D} fourier transform magnetic resonance imaging.
\newblock {\em Magnetic Resonance in Medicine 25\/} (1992), 319--333.

\bibitem{plimpton91}
{\sc Plimpton, S., Mastin, G., and Ghiglia, D.}
\newblock Synthetic aperture radar image processing on parallel supercomputers.
\newblock In {\em Proceedings of Supercomputing '91\/} (Albuquerque, NM,
  November 1991), pp.~446--452.

\bibitem{shaw92}
{\sc Shaw, G., Gabel, R., Martinez, D., Rocco, A., Pohlig, S., Gerber, A.,
  Noonan, J., and Teitelbaum, K.}
\newblock Multiprocessors for radar signal processing.
\newblock Tech. Rep. 961, MIT Lincoln Laboratory, Nov. 1992.

\bibitem{stichnoth93a}
{\sc Stichnoth, J.}
\newblock Efficient compilation of array statements for private memory
  multicomputers.
\newblock Tech. Rep. CMU-CS-93-109, School of Computer Science, Carnegie Mellon
  University, Feb. 1993.

\bibitem{stichnoth94}
{\sc Stichnoth, J., O'Hallaron, D., and Gross, T.}
\newblock Generating communication for array statements: Design,
  implementation, and evaluation.
\newblock {\em Journal of Parallel and Distributed Computing 21}, 1 (Apr.
  1994), 150--159.

\bibitem{ISCA-MEMORY}
{\sc Stricker, T., and Gross, T.}
\newblock Optimizing memory system performance for communication in parallel
  computers.
\newblock In {\em Proc. 22nd Intl. Symposium on Computer Architecture\/} (Santa
  Marguerita di Ligure, Italy, June 1995), ACM, pp.~308--319.

\bibitem{subhlok94a}
{\sc Subhlok, J., O'Hallaron, D., Gross, T., Dinda, P., and Webb, J.}
\newblock Communication and memory requirements as the basis for mapping task
  and data parallel programs.
\newblock In {\em Proc. Supercomputing '94\/} (Washington, DC, Nov. 1994),
  pp.~330--339.

\bibitem{subhlok93}
{\sc Subhlok, J., Stichnoth, J., O'Hallaron, D., and Gross, T.}
\newblock Exploiting task and data parallelism on a multicomputer.
\newblock In {\em Proc. of the ACM Symposium on Principles and Practice of
  Parallel Programming (PPoPP)\/} (San Diego, CA, May 1993), pp.~13--22.

\bibitem{webb93bb}
{\sc Webb, J.}
\newblock Latency and bandwidth consideration in parallel robotics image
  processing.
\newblock In {\em Supercomputing '93\/} (Nov. 1993), pp.~230--239.

\bibitem{yang93}
{\sc Yang, B., Webb, J., Stichnoth, J., O'Hallaron, D., and Gross, T.}
\newblock {D}o\&{M}erge: Integrating parallel loops and reductions.
\newblock In {\em Proc. Sixth Workshop on Languages and Compilers for Parallel
  Computing\/} (Portland, OR, Aug. 1993), vol.~768 of {\em Lecture Notes in
  Computer Science}, Springer Verlag, pp.~169--183.

\end{thebibliography}
