\newcommand{\comment}[1]{}
\input epsf
\makeatletter
\long\def\unmarkedfootnote#1{{\long\def\@makefntext##1{##1}\footnotetext{#1}}}
\makeatother
\documentstyle[prog,cmu-titlepage]{cmu-art}
\include{psfig}
%\include{caption}

\renewcommand{\floatpagefraction}{1}
\renewcommand{\topfraction}{1}
\renewcommand{\bottomfraction}{1}
\renewcommand{\textfraction}{0}
%\renewcommand{\baselinestretch}{1}

\newtheorem{definition}{Definition}
\newtheorem{theorem}{Theorem}
\newtheorem{lemma}{Lemma}
\newcommand{\ignore}[1]{}

%\textfloatsep = 0.2in
%\intextsep = 0.2in
%\newcommand{\Tau}{{\large\tau\normalsize\mit}}
\newcommand{\Tau}{\tau}
\newcommand{\Mu}{\mu}
\newcommand{\Chi}{{\cal X}}
\begin{document}
\bibliographystyle{acm}
\title{The CMU Task Parallel Program Suite}
\author{Peter Dinda \and Thomas Gross \and David O'Hallaron \and Edward Segall \and James Stichnoth Jaspal Subhlok \and Jon Webb \and Bwolen Yang}
\trnumber{CMU-CS-94-131}

\date{March 1994}

%\vspace{0.0 in}
\thispagestyle{empty}

%\vspace*{-0.2 in}
\abstract{
The idea of exploiting both task and data parallelism in programs is
appealing. However, identifying realistic, yet manageable example
programs that can benefit from such a mix of task and data parallelism
is a major problem for researchers. We address this problem by
describing a suite of five applications from the domains of
scientific, signal, and image processing that are of reasonable size,
are representative of real codes, and can benefit from exploiting task
and data parallelism. The suite includes fast Fourier transforms,
narrowband tracking radar, multibaseline stereo imaging, and airshed
simulation.  Complete source code for each example program is
available from the authors.  }

\unmarkedfootnote{
This research was sponsored by the Advanced Research Projects
Agency/CSTO monitored by SPAWAR under contract N00039-93-C-0152, and
also by the Air Force Office of Scientific Research under contract
F49620--92--J--0131.  The views and conclusions contained in this
document are those of the authors and should not be interpreted as
representing the official policies, either express or implied, of the
U.S. government.  

Contact author: David O'Hallaron, droh@cs.cmu.edu.
}

%\normalsize
\keywords{Parallel programming, task parallelism, functional parallelism,
parallelizing compilers, program mapping
}
\maketitle

\section{Introduction}
There is growing interest in the idea of exploiting both task and data
parallelism
~\cite{agrawal93,chandy93,chapman94,cheung93,crovella93,foster92,subhlok93}.
There are a number of practical reasons for this interest.  For many
applications, especially in the domains of image and signal
processing, data set sizes are limited by physical constraints and
cannot be easily increased~\cite{subhlok93}. In such cases the amount
of available data parallelism is limited. For example, in the
multibaseline stereo application described later in this report, the
size of an image is determined by the circuitry of the video cameras
and the throughput of the camera interface.  Increasing the image size
means buying new cameras and building a faster interface, which may
not be feasible. Since the data parallelism is limited, additional
parallelism must come from tasking.

Another reason for the increased interest in task parallelism is that
simulations are becoming increasingly sophisticated as they attempt to
capture interactions among different physical phenomena.  The
phenomena might represent different scientific disciplines, and
different parts of the simulation might even be written by different
groups.  For example, the airshed model described later in this report
characterizes the formation of air pollution as the interaction
between the wind blowing and reactions among various chemical species.
It is natural to model such interactions using tasking, where one task
models the wind blowing, and the other task models the chemical
reactions.  Further, if the codes are written by different groups,
task parallelism may be the only feasible way to integrate the codes.

Applications that can benefit from a mix of task and data parallelism
tend to be somewhat complicated because they are typically composed of
a collection of nontrivial functions, each of which is a sequence of
data parallel operations.  Identifying and building representative
example programs that are a manageable size is a major stumbling block
for computer science researchers who are not application domain
experts.  We address this problem by describing a set of realistic
example programs from the domains of scientific, signal, and image
processing that can benefit from a mix of task and data parallelism:
\begin{enumerate}
\item 1D fast Fourier transform. 
\item 2D fast Fourier transform. 
\item Narrowband tracking radar. 
\item Multibaseline stereo. 
\item Airshed simulation.  
\end{enumerate}
We identified these applications in the course of developing an
integrated model of task and data parallelism for the Fx parallelizing
Fortran compiler~\cite{gross94,subhlok93,subhlok93b,subhlok94,yang93}
and have found them to be extremely helpful.

Complete Fortran 77 sources of the programs are available from the
authors. Each program is fewer than 500 lines of code.  The Fortran 77
sources are useful for a number of reasons.  First, the sources
provide an unambiguous specification of each application, including
input and output data sets. Second, there are many models and dialects
for task parallelism.  Fortran 77 represents a lowest common
denominator of sorts, available to everyone, for describing the
applications. Finally, the source code clearly identifies the obvious
sources of course-grained parallelism in the form of calls to
subroutines.  This enabled us, in all but one case, to port the
Fortran 77 programs to the Fx system with only minor modifications
(the only exception being the airshed simulation, which requires a
dynamic model of task parallelism that Fx does not currently support).
We invite other researchers who are interested in task parallelism to
use these Fortran 77 sources as the basis for writing the applications
in their favorite task parallel dialect.

Section~\ref{sec:overview} briefly outlines a simple space of the
different ways that task and data parallelism can be used in the same
program. Sections~\ref{sec:fft}--\ref{sec:airshed} describe the
example programs and discuss briefly how they can mapped onto a
parallel system using a mix of task and data parallelism.
Section~\ref{sec:distribution} describes where to find the online
Fortran 77 codes and provides some more detail on their structure.

\section{Combinations of task and data parallelism}
\label{sec:overview}
Each program in our task parallel suite has the following form:
\tt
\begin{tabbing}
ttttttt\=ttt\=ttt\= ttt\=ttt\= \kill
\>     do i = 1,m \\
\>\>     call T1()\\
\>\>     call T2()\\
\>\>     call T3()\\
\>     enddo\\
\end{tabbing}
\rm
There is an outer loop that iterates over $m$ input data sets.  The
body of the loop consists of calls to three {\em task-subroutines}.
Each task-subroutine typically consists of a sequence of data-parallel
statements (e.g., array assignment statements or DOALL-like parallel
loops) that can run on multiple processors. Depending on the
application, there may or may not be dependences across iterations of
the outer loop and among the task-subroutines.  We can graphically
represent the program using a task graph, as shown in
Figure~\ref{fig:example}. Each node corresponds to a task-subroutine,
and each arc corresponds to a data dependence between a pair of
task-subroutines.
\begin{figure}[ht]
\centerline{\epsfbox{example.eps}}
\caption{Example task graph}
\label{fig:example}
\end{figure}

We can view different task and data parallel mappings for an
application as points in a two-dimensional space~\cite{subhlok94}.
The first axis corresponds to different {\em clusterings} of
task-subroutines, where each task-subroutine in cluster runs on the
same set of processors, and each cluster runs on a unique set of
processors.  In some cases, assigning two task subroutines to the same
cluster reduces communication overhead, at the cost of reduced
parallelism. In these cases, there is a tension between reducing
communication cost and increasing parallelism.
\label{sec:fft}
\begin{figure}[ht]
\centerline{\epsfbox{forms.eps}}
\caption{Combinations of task and data parallel mappings}
\label{fig:forms}
\end{figure}

If the input data sets for a cluster are independent, then that
cluster can be {\em replicated}.  For example, one replicated instance
of a cluster could process even numbered data sets and the other
replicated instance could process the odd numbered data sets.  The
degree of replication is captured by the second axis. Replication is
beneficial for task-subroutines that do not scale well.

Some examples of different combinations of task and data parallelism
are shown in Figure~\ref{fig:forms}.  Figure~\ref{fig:forms}(a) shows
the usual data parallel mapping where all task-subroutines are
clustered onto all processors, with no replication.
Figure~\ref{fig:forms}(b) shows a task parallel mapping where each
task-subroutine gets its own cluster, again with no replication. It
may be desirable in some cases to replicate the data-parallel
clusters, as in Figure~\ref{fig:forms}(c).  Finally, a mix of task and
data parallelism, using both clustering and replication, is shown in
Figure~\ref{fig:forms}(d).

\section{Fast Fourier transform}
The fast Fourier transform (FFT) converts an input data set from the
temporal/spatial domain to the frequency domain, and vice versa.
While it is an important algorithm in its own right, with numerous
applications in scientific, signal, and image processing, it is an
especially interesting example program for studying task parallelism
because it exemplifies a common computational pattern, i.e.,
manipulate the data row-wise, then manipulate the data column-wise.
This pattern appears in diverse applications such as medical
imaging~\cite{noll92}, synthetic aperture radar
imaging~\cite{plimpton91}, ADI solvers, sonar beamforming, and the
radar and airshed programs described later in this document.

More precisely, a discrete Fourier transform (DFT) is a complex
matrix-vector product
\begin{displaymath}
y = F_n x 
\end{displaymath}
where $x$ and $y$ are complex vectors, and $F_n = (f_{pq})$ is an $n
\times n$ matrix such that $f_{pq} = \omega_n^{pq}$, where 
\begin{displaymath}
\omega_n = \cos(2\pi/n) - i \sin(2\pi/n) = e ^ {-2 \pi i /n } 
\end{displaymath}
and $i = \sqrt{-1}$.  A fast Fourier transform (FFT) is an efficient
procedure for computing the DFT that exploits structure in $F_n$ to
compute the DFT matrix-vector product in $O(n \log n)$ time.  Higher
dimensional FFTs are also defined. In general, an $m$-dimensional FFT
operates on an $m$-dimensional matrix with $n$ elements in each
dimension in $O(n^m \log n)$ time.  See \cite{vanloan92} for an
excellent description of the numerous FFT algorithms that have been
developed over the past 40 years.

\subsection{1D fast Fourier transform}
If $n=n_1 n_2$, then a 1D~FFT can be computed using a collection of
smaller independent 1D~FFTs~\cite{bailey90,gentleman66,vanloan92}.
Starting with a view of the input vector $x$ as an $n_1 \times n_2$
array $A$, stored in column major order, perform $n_1$ independent
$n_2$-point FFTs on the rows of $A$, scale each element $a_{pq}$ of
the resulting array by a factor of $\omega_n^{pq}$, and then perform
$n_2$ independent $n_1$-point FFTs on the columns of $A$. The final
result vector $y$ is obtained by concatenating the rows of the
resulting array.  Figure~\ref{fig:fft1} shows the task graph for the
1D~FFT.
\begin{figure}[ht]
%\epsfxsize=4in
\centerline{\epsfbox{fft1.eps}}
\caption{1D FFT task graph for one input vector}
\label{fig:fft1}
\end{figure}

Input and output are sequences of vectors reshaped as 2D arrays.
Nodes labeled {\tt trans} perform a transpose operation, nodes labeled
{\tt col FFTs} perform a set of 1D FFTs on the columns of its input
array, and the node labeled {\tt scale} multiplies each element of its
input array by a constant.  To exploit locality in the memory
subsystem, the program implements each set of row-wise operations as a
transpose followed by a set of column-wise operations.  This specific
order is an artifact of the fact that Fortran 77 stores matrices in
column-major order. If the example program were written in C, each
column-wise operation would be implemented as a transpose followed by
a set of row-wise operations.

Mapping the 1D~FFT onto a parallel system is easy in some ways and
challenging in other ways.  The problem is easy in the sense that the
column-wise FFTs and the scaling operation are perfectly parallel and
easily represented by conventional data parallel constructs.  The
column-wise FFTs are naturally expressed with DOALL-like parallel loop
construct such as the HPF independent DO statement~\cite{hpf93}. The
scaling operation can be expressed with a Fortran 90 array assignment
statement.  However, the problem is challenging because the transpose
operation requires an efficient redistribution of data, usually in
the form of a complete exchange where each processor must send data to
every other processor. Since the input arrays are independent, both
replication and clustering of the task graph are possible.

\subsection{2D fast Fourier transform}
Computing the 2D~FFT is similar to computing the 1D~FFT.  Given an
$n_1 \times n_2$ input array $A$, perform $n_2$ independent
$n_1$-point FFTs on the columns of $A$, followed by $n_1$ independent
$n_2$-point 1D~FFTs on the rows of A.  Figure~\ref{fig:fft2} shows
the task graph for the 2D~FFT.  As with the 1D~FFT, row-wise FFTs are
replaced by a transpose followed by a set of column-wise FFTs. Notice
that the 2D~FFT is simpler than the 1D~FFT.  No scaling is required
and there is one fewer transpose operation. Again, the input and output
are sequences of arrays, and since these arrays are independent, both
replication and clustering of the task graph are possible.

\begin{figure}[htb]
%\epsfxsize=4in
\centerline{\epsfbox{fft2.eps}}
\caption{2D FFT task graph for one input array}
\label{fig:fft2}
\end{figure}

\section{Narrowband tracking radar}
\label{sec:radar}
The narrowband tracking radar benchmark was developed by researchers
at MIT Lincoln Laboratories to measure the effectiveness of various
multicomputers for their radar applications~\cite{shaw92}.  It is an
interesting program for studying task parallelism because of its hard
real-time requirements, and because the size of the input data set is
limited by physical properties of the radar sensor.  The task graph
for the radar application is shown in Figure~\ref{fig:radar}.
\begin{figure}[tbh]
%\epsfxsize=4in
\centerline{\epsfbox{radar.eps}}
\caption{Radar task graph for one input array}
\label{fig:radar}
\end{figure}

The program inputs data from a
single sensor along $c=4$ independent {\em channels}.  Every 5
milliseconds, for each channel, the program receives $d=512$ complex
vectors of length $r=10$, one after the other in the form of an $r
\times d$ complex array $A$ (assuming the column major ordering of
Fortran).  At a high-level, each input array $A$ is processed in the
following way: (1) {\em Corner turn} the $r \times d$ input array to
form a $d \times r$ array.  (2) Perform $r$ independent $d$-point
FFTs.  (3) Convert the resulting complex $d \times r$ array to a real
$w \times r$ subarray, $w=40$, by replacing each element $a + ib$ in
the $w \times r$ subarray with its scaled magnitude $\sqrt{a^2 +
b^2}/d$.  (4) Threshold each element $a_{jk}$ of the subarray using a
cutoff that is a function of $a_{jk}$ and the sum of the subarray
elements.  Elements that are above the threshold are set to unity;
elements below the threshold are set to zero.

The corner turn operation is equivalent to a transpose, so it can
potentially induce a complete exchange where each processor
communicates with every other processor. As with the 1D~FFT, the
column FFTs, scaling, and thresholding operations can be naturally
expressed using conventional data parallel constructs. Further, the
reduction operation requires an efficient reduction mechanism.
However, the most interesting computational
property of the radar benchmark is the fact that the size parameters
$r$, $d$, $c$, and $w$ are determined by mother nature and the
properties of current sensor technology.  The luxury of simply
increasing the data set size simply does not exist in this case.  The
amount of available low-level data parallelism is limited, so
additional parallelism must come from higher-level task parallelism.
Like the FFT examples, input data sets are independent, so both 
replication and clustering of the task graph are possible.

\section{Multibaseline stereo}
\label{sec:stereo}
The multibaseline stereo uses an algorithm developed at Carnegie
Mellon that gives greater accuracy in depth through the use of more
than two cameras~\cite{okutomi93}. It is an interesting program for
studying task parallelism because it contains significant amounts of
both inter-task and intra-task communication\cite{webb93b}, and
because, like the radar example, the size of the input data sets
cannot be easily increased. Our implementation is adapted from a
previous data-parallel implementation written in a specialized image
processing language~\cite{webb93a}.

Figure~\ref{fig:stereo} shows the task graph for the stereo
program.  Input consists of three $m \times n$ images
acquired from three horizontally aligned, equally spaced cameras. One
image is the {\em reference image}, the other two are {\em match
images}. For each of 16 disparities, $d = 0,\ldots,15$, the first
match image is shifted by $d$ pixels, the second image is shifted by
$2d$ pixels. A {\em difference image} is formed by computing the sum
of squared differences between the corresponding pixels of the
reference image and the shifted match images.  Next, an {\em error
image} is formed by replacing each pixel in the difference image with
the sum of the pixels in a surrounding $13 \times 13$ window. A {\em
disparity image} is then formed by finding, for each pixel, the
disparity that minimizes error.  Finally, the depth of each pixel is
displayed as a simple function of its disparity.

\begin{figure}
%\epsfxsize=4in
\centerline{\epsfbox{stereo.eps}}
\caption{Multibaseline stereo task graph for one input data set}
\label{fig:stereo}
\end{figure}

The stereo program requires efficient mechanisms for broadcasting and
reducing large data sets. The computation of the difference images
requires simple pointwise operations on the three input images and can
thus be naturally expressed with Fortran 90 array statements.  The
computation of the error images is somewhat more interesting, being
similar to a convolution operation.  The convolution can be modeled as
a DOALL where the loop iterations operate on overlapping
regions of the image, which means that processors must communicate
before the loop iterations can begin executing. As with the FFT and
radar programs, the data sets are independent, so both replication and
clustering of the task graph are possible.

\section{Airshed simulation}
\label{sec:airshed}
The airshed simulation is significantly more complex than the previous
examples.  The multiscale airshed model captures the formation,
reaction, and transport of atmospheric pollutants and related chemical
species~\cite{mcrae82,mcrae92}. It is an interesting application
because it requires a dynamic task parallel model, and because
different parts of the application exhibit widely varying amounts of
DOALL parallelism.

The airshed application simulates the behavior of the airshed model
when it is applied to $s$ chemical species, distributed over domains
containing $p$ grid points in each of $l$ atmospheric layers.  Typical
values are $s = 35$ species, $500 \leq p \leq 5000$ grid points, and
$l = 5$ atmospheric layers.  Because of the multiscale grid, the
entire northeastern United States can be modeled with problems in this
size range. A total of about 200 chemical reactions are modeled.

The program computes in two principle phases: (1) horizontal transport
(using a finite element method with repeated application of a direct
solver), followed by (2) chemistry/vertical transport (using an
iterative, predictor-corrector method).  Figure~\ref{fig:airshed}
depicts the task graph for one hour of simulated time.  Input is an $l
\times s \times p$ concentration array. Initial conditions are input
from disk ({\tt inputhour}), and in a preprocessing phase for the
horizontal transport phases to follow, the finite element stiffness matrix for
each layer is assembled and factored ({\tt pretrans}). The atmospheric
conditions captured by the stiffness matrix are assumed to be constant
during the simulated hour, so this step is performed just once per
hour.  This is followed by a sequence of steps --- the number of steps
is one of the initial conditions --- where each step consists of a
horizontal transport phase, followed by a chemistry/vertical transport
phase, followed by another horizontal transport phase.  Each
horizontal transport phase performs $ls$ backsolves, one for each
layer and species. All may be computed independently; however, for
each layer $l$, all backsolves use the same factored matrix $A_l$.
The chemistry/vertical transport phase performs an independent
computation for each of the $p$ grid points. Output for the hour is an
updated concentration array, which is then input to the next hour.
\begin{figure}
%\epsfxsize=4in
\centerline{\epsfbox{airshed.eps}}
\caption{Task graph for one hour of the airshed simulation}
\label{fig:airshed}
\end{figure}

A number of interesting issues arise when we map the airshed to a
parallel system. In the other example programs we have discussed, the
number of tasks is known at compile time. However, in the airshed
program, the number of transport/chemistry/transport steps for each
hour is not known until runtime, which implies a dynamic model of task
parallelism. Also, since the output concentration array of one hour is
the input to the next hour, replication of the task graph is not
feasible, as it was with the previous example programs.

Another issue is that the preprocessing phase, the transport phase,
and the chemistry phase have very different levels of obvious DOALL
parallelism because the sizes of the different dimensions of the
concentration array differ by orders of magnitude. For example, the
preprocessing phase independently computes stiffness matrices for each
layer; unfortunately there are only 5 layers, so the obvious DOALL
approach will use 5 processors. To get better utilization, we must 
parallelize the computation for each layer, or we must try to employ
task parallelism to pipeline the computation for each layer, or both.

The issues involved in mapping the transport phase are particularly
interesting. Since there are $l=5$ layers and $s=35$ species, the
transport phase could be easily implemented with doubly nested DOALLs
that consist of 175 independent loop iterations. For moderate sized
parallel systems, with say 64 processors, this approach might work
well.  However, for larger systems, with say 512 processors, this
approach uses only a fraction of the processors. As with the
preprocessing phase, we can get better utilization by either
parallelizing the sparse finite element computation (a difficult task)
or trying to use task parallelism to pipeline the computation.

The final issue stems from the fact that the preprocessing, horizontal
transport, and chemistry/vertical transport phases each operate on
different dimensions of the concentration array. To exploit locality
in the memory hierarchy, an implementation will most likely insert the
appropriate transpose operation before each phase. On a parallel
system, this can induce a complete exchange where each processor
communicates with every other processor. Again, as with the FFT and
radar examples, we see the need for an efficient complete exchange
mechanism.

\section{Distribution}
\label{sec:distribution}
Complete Fortran 77 sources for the application described in
this report are available via anonymous FTP from {\tt warp.cs.cmu.edu}
in file {\tt fx-codes/tpsuite/tpsuite.tar}. World Wide Web clients
like Mosaic and Cello can use the following URL:
\begin{quote}
{\tt ftp://warp.cs.cmu.edu/usr/anon/fx-codes/tpsuite/tpsuite.tar}
\end{quote}
Each source program is a vanilla Fortran 77 code that operates on a
sequence of inputs and produces a sequence of outputs. There are no
data files; all input data sets are produced automatically by the
program.  Each program has fewer than 500 lines of code and consists
of a single source and include file. The include file contains size
constants that can be changed if the researcher wants to measure
scalability. Each program (except for radar) checks its output
automatically. The radar code prints a few lines of output which can
be easily verified by the user; directions are provided in the source
code.  Sample outputs of the programs running on a DEC 3000 Alpha
workstation are also provided.  Each program (except the airshed)
computes physically meaningful results. To keep the program at a
manageable size, we have provided a version of the airshed simulation that
uses a synthetic workload for the innermost loops.

\section{Concluding remarks}
We have described a suite of realistic programs that can benefit from
a mix of task and data parallelism. Researchers in the areas of
mapping, parallelizing compilers, and parallel programming
environments are invited to use these programs to test and validate
their ideas for exploiting task and data parallelism in applications.


\section*{Acknowledgements}
Many thanks to Takeo Kanade and his group at Carnegie Mellon for
developing the original stereo code, to Dennis Ghiglia and his imaging
group at Sandia National Laboratories who, in generously providing us
with source code for their spotlight SAR application, showed us the
importance of the 2D FFT as a research tool, and to G. Shaw and fellow
researchers at MIT Lincoln Labs for providing C source code for the
radar benchmark in their technical report.

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


\end{document}


