\makeatletter
\long\def\unmarkedfootnote#1{{\long\def\@makefntext##1{##1}\footnotetext{#1}
}}\makeatother
\newcommand{\comment}[1]{}
\documentstyle[times,11pt,epsf,fullpage]{article}

\setcounter{topnumber}{4}
\renewcommand{\topfraction}{1}
\setcounter{bottomnumber}{4}
\renewcommand{\bottomfraction}{1}
\setcounter{totalnumber}{6}
\renewcommand{\textfraction}{0.1}
\renewcommand{\floatpagefraction}{0.9}

\begin{document}
\bibliographystyle{acm}

\title{
\begin{flushleft}
\vspace{-0.9in}
\scriptsize {\em
Working paper, December, 1996.
}
\end{flushleft}
\vspace{0.4in}
\Large\bf 
Communication Performance of Irregular Applications and\\
the Implications for Graph Partitioners
}

\author{David R. O'Hallaron and Jonathan Richard Shewchuk\\
\\
School of Computer Science\\
Carnegie Mellon University\\
Pittsburgh, PA 15213\\
\{droh, jrs\}@cs.cmu.edu
}
\date{}

\thispagestyle{empty}
\maketitle

\unmarkedfootnote{
This research was sponsored in part by the Advanced Research Projects
Agency/CSTO monitored by SPAWAR under contract N00039-93-C-0152, 
in part by the National Science Foundation under
Grant CMS-9318163, and in part by a grant from the Intel Corporation.
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.
}

\begin{abstract}
We derive a simple performance model for a family of sparse irregular
applications and use the model to predict running times on the Cray
T3D.  Based on insights gained from applying the model to the
applications, we identify some limitations of modern graph
partitioners and suggest some improvements. In particular, we suggest
an approach whereby a partitioner could use our model to attempt to
minimize running time on a specific system.
\end{abstract}

\section{Introduction}
\label{sec:intro}
In general it is important to understand the communication
requirements of real applications~\cite{iscanbody}, and these
requirements are especially difficult to understand for the important
class of irregular scientific applications that manipulate sparse
matrices.  Irregular applications include those that model processes in nature
like the flow of air over a wing, or the distribution of stresses and
strains in a bridge.  Characterizing realistic irregular applications
is difficult, primarily because they simulate complex
physical phenomena operating over complex physical domains, and thus
tend to be large, complex, and poorly understood.  

This paper characterizes the communication performance of a family of
large irregular applications and uses this characterization to
identify some limitations in modern graph partitioning
approaches.  Specifically, we are concerned with a family of {\em
Quake applications} that model the motion of the ground during strong
earthquakes in the San Fernando Valley of Southern California. The
Quake applications, described in Section~\ref{sec:quake}, were
developed over a period of two years by civil engineers at Carnegie
Mellon as part of an ongoing project to model the entire Los Angeles
basin~\cite{superquake96}.

A nice property of the Quake applications is that their running time is
dominated by a sparse matrix-vector product (SMVP) operation that is
repeated thousands of times. Further, the SMVP is the only operation
in the application that requires the transfer of data between
processors.  Thus, we can focus on understanding the behavior of the SMVP
operation and abstract away much of the complexity of the application.

Starting from a simple performance model of the SMVP operation, we
describe an approach for deriving the application-specific and
system-specific parameters in the model, and then use them to predict
the running time on the Cray T3D of one of the Quake applications.
Our performance model for the SMVP is described in
Section~\ref{sec:model}, and the model is validated in
Section~\ref{sec:validation}.

Then in Section~\ref{sec:implications} we use the model to identify a
limitation in current graph partitioners and to suggest a different
approach. In general, modern graph partitioners attempt to minimize
the global volume of data transferred, subject to the constraint that
each PE has an equal amount of work. This is a convenient metric for
algorithm designers, but it is not clear how closely related it is to
actual running time of parallel programs, because of two
shortcomings: (1) It
ignores the number of neighbors that each PE must communicate with,
and (2) it does not attempt to balance communication across
PEs.  Shortcoming (1) is significant because, according to our model, real
applications can be either dominated by block latency, dominated by
burst bandwidth during block transfers, or somewhat balanced between
the two.  The choice depends on the properties of the underlying system.
In cases where the communication phase is latency-dominated, the
partitioner should attempt to minimize the maximum
number of neighbors of each subgraph
(and thus the number of distinct messages that
each PE must send).  Shortcoming (2) is
significant because
the communication time can vary significantly from PE to PE
when using modern graph partitioners, and the running time of the
slowest PE determines the overall running time.

The general problem is that modern graph partitioners do not
attempt to minimize actual running time.  In
Section~\ref{sec:implications} we suggest a system-specific approach
based on our simple performance model. With this approach, the
partitioner uses the model in Section~\ref{sec:model} to attempt
to minimize overall running time, using the method outlined in
Section~\ref{sec:validation} to determine the system-specific
bandwidth and latency parameters.  We are not outlining a new
partitioning algorithm, but we are suggesting a system-dependent
objective function that future partitioners should seek to optimize.

\section{The Quake applications}
\label{sec:quake}

The paper is motivated by our experience with a family of four {\em
Quake applications}, denoted sf10, sf5, sf2, and sf1. Each is a
finite-element simulation of the motion of the ground in the San
Fernando Valley during a strong earthquake~\cite{superquake96}.  The
digit in each name indicates the highest frequency wave (in seconds)
that the simulation is able to resolve.  For example, sf10 resolves
waves with 10 second periods, sf5 resolves waves with 5 second
periods, and so on.  Each application simulates 60 seconds of shaking
as shock waves travel through a model of the San Fernando Valley.
Each model employs a three-dimensional unstructured finite element
mesh composed of thousands or millions of tetrahedra (i.e., pyramids
with triangular bases). The model for sf10 is illustrated in
Figure~\ref{fig:sf10mesh-all}.  The model corresponds to a volume of
earth roughly 50 km x 50 km x 10 km, and the smallest distance between
nodes is roughly 50 m.  Beverly Hills is in the lower right-hand
corner. The town of San Fernando is in the midst of the darkly shaded
region near the upper left corner.

\begin{figure}[htbp]
\centerline{
\epsfxsize=4in
\epsfbox{sf10.ele.all.bw.epsf}
}
\caption{Finite element mesh for the sf10 model of the San Fernando Valley.}
\label{fig:sf10mesh-all}
\end{figure}

Each tetrahedron in Figure~\ref{fig:sf10mesh-all} is called an {\em
element}, and the vertices of the tetrahedra are called {\em
nodes}.  Some finite element simulations use structured meshes
constructed from regular grids; however, the Quake applications
require {\em unstructured} meshes, which can accommodate
the wildly varying density
of the soils in the valley. Wavelengths are shorter in softer soils,
thus softer soils need a higher density of elements.

Figure~\ref{fig:general} lists the properties of the San Fernando
meshes. When the wave period is halved, its frequency doubles,
and the number of nodes increases
by a factor of nearly eight---a factor of two in each of three dimensions.
As a general rule, for each node in the mesh a simulation uses
about 2 KBytes of memory at runtime to accommodate the storage
of several vectors and sparse matrices.  Thus sf10
requires about 15 MBytes and sf1 requires about 5 GBytes.

\begin{figure}[tbp]
\begin{center}
\begin{tabular}{|l||r|r|r|r|}
\hline
\multicolumn{1}{|c||}{property} & 
\multicolumn{1}{|c|}{sf10} & 
\multicolumn{1}{|c|}{sf5} & 
\multicolumn{1}{|c|}{sf2} & 
\multicolumn{1}{|c|}{sf1}\\
\hline
\hline
nodes        & 7,294     & 30,169   & 378,747    & 2,461,694\\
edges        & 44,922 	 & 190,377  & 2,509,064  & 16,684,112\\
faces        & 72,698	 & 311,514  & 4,198,057  & 28,202,581\\
elements     & 35,025	 & 151,239  & 2,067,739  & 13,980,162\\
\hline
\end{tabular}
\end{center}
\caption{Properties of the Quake applications.}
\label{fig:general}
\end{figure}

During each of the 6000 time steps, a Quake simulation executes three
sparse matrix-vector product (SMVP) operations of the
form ${\bf y = Kx}$, where ${\bf x}$ and ${\bf y}$ are vectors of
length $3n$ (representing three degrees of freedom---$x$, $y$, and $z$
displacements---for each node of the mesh), and ${\bf K}$ is a sparse
$3n \times 3n$ {\em stiffness matrix}.
${\bf K}$ can be likened to an adjacency matrix of the nodes of the mesh;
${\bf K}$ contains a $3 \times 3$ submatrix for each pair of nodes
connected by an edge of the mesh (including self-edges).
The stiffness matrix is extremely sparse; each node is connected to an
average of $13$ neighbors (in addition to itself), so
each row of ${\bf K}$ contains an average of $14 \times 3 = 42$
nonzero numbers~\cite{meshsuite96}.

The simulations are parallelized using a domain-specific tool chain
for finite element problems called Archimedes~\cite{superquake96}.  To
generate a simulation that will run on $p$ PEs, Archimedes
partitions the mesh into $p$ disjoint sets of elements. Each set is
called a {\em subdomain} and is assigned to some PE. The
partitioner is based on a recursive geometric bisection
algorithm~\cite{miller92} that divides the elements equally among the
subdomains while attempting to minimize the total number of nodes that
are shared by subdomains, and hence the total communication volume.
The geometric partitioning algorithm has provable upper bounds on the
separator sizes and in practice usually generates partitions that are
as good as those produced by other modern partitioning
algorithms~\cite{barnard92,GiMiTe95,chaco95}.

To compute ${\bf y = Kx}$ on a set of PEs, we must consider the
data distribution by which vectors and matrices are stored.  The
vectors ${\bf x}$ and ${\bf y}$ are stored in a distributed fashion
according to the mapping of nodes to PEs induced by the partition
of elements among PEs.  If a node $i$ resides in several
PEs (because $i$ is a vertex of several elements mapped to
different PEs), the values $x_i$ and $y_i$ are replicated on those
PEs.  The matrix ${\bf K}$ is distributed so that $K_{ij}$
resides on any PE on which nodes $i$ and $j$ both reside.
Figure~\ref{fig:cutmatrix} demonstrates this method of distributing data.

\begin{figure}[tbp]
\centerline{\epsfbox{cutmatrix.ps}}
\caption{A finite element mesh and corresponding stiffness matrix ${\bf K}$,
         distributed among two PEs.  Xs represent nonzero $3 \times 3$
         submatrices.  Note that some nodes are shared by both PEs,
         as are some stiffness matrix entries (corresponding to shared edges).}
\label{fig:cutmatrix} 
\end{figure}

Given this method of distributing data, the multiplication ${\bf y =
Kx}$ is performed in two steps.  First, each PE computes a
local matrix-vector product over the subdomain that resides on that
PE. Second, PEs that share nodes communicate and combine
their nodal ${\bf y}$ values into correct global values for each node.
In Figure~\ref{fig:cutmatrix}, PEs 1 and 2 must communicate to
resolve the values of the shared nodes 4, 5, and 6.


\section{Performance model}
\label{sec:model}

As we have seen, each global SMVP consists of a simple sequence,
executed simultaneously on each PE: a local SMVP (the
computation phase) followed by an operation in which each PE
exchanges and sums data with each of its neighbors (the communication
phase).  The computation and communication phases are synchronous,
i.e., each phase begins with a global barrier synchronization.  (It is
in principle possible, with much difficulty,
to improve performance by overlapping the
computation and communication phases, but the real implementations of
the Quake SVMP do not do this, so the models in this section
will follow suit.)  Further, we assume that during the communication
phase, the summing operations are overlapped with the exchanges, and
that the running time of the communication phase is determined by the
time to perform the exchanges.  This is a reasonable assumption, as
the summing operations can be handled by the PE and the exchanges by a
network interface.  Given these assumptions, the running time for each
global SMVP is
\begin{displaymath}
T_{mv} = T_{comp} + T_{comm}
\end{displaymath}
where $T_{comp}$ is the time required for all PEs to finish their
computation phases and $T_{comm}$ is the time required for all PEs
to finish their communication phases. (For reference,
Figure~\ref{fig:symbols} summarizes all symbols introduced in this
section.)

\begin{figure}[tbp]
\begin{center}
\small
\begin{tabular}{|l||l|}
\hline
Symbol & Description\\
\hline
\hline
$T_{mv}$   &running time for sparse matrix-vector product (SMVP)\\
$T_{comp}$ &running time for computation phase of SMVP\\
$T_{comm}$ &running time for communication phase of SMVP\\
\hline
\hline
$F$ &flops per processing element (PE)\\
$C$ &number of communication words per PE\\
$B$ &number of blocks per PE\\ 
$\beta$ & communication balance\\
\hline
\hline
$T_f$ &time per flop (inverse of sustained FLOPS)\\
$T_c$ &time per communication word (inverse of sustained bandwidth)\\
$T_l$ &time per block (block latency)\\
$T_w$ &time per additional block word (inverse of block bandwidth)\\
\hline
\end{tabular}
\end{center}
\caption{Symbols.}
\label{fig:symbols}
\end{figure}

The unit of work for the computation phase is a floating-point
operation ({\em flop}), which is either a scalar add or multiply.
Although the flop is not an ideal way to measure work in every application,
it is reasonable for our SMVP operation because the (minimum) number of
flops each PE must perform can be counted exactly and is
independent of the compiler and its optimization level. If a PEs local
matrix has $m$ nonzeros, then the local SMVP requires at
least $F = 2m$ flops; one add and one multiply for each nonzero. 
If each flop requires an average time of $T_f$, then the running time
for a local SMVP is 
\begin{displaymath}
T_{comp} = FT_f. 
\end{displaymath}
Since $T_f$ includes all hardware and software overheads (e.g., loads,
stores, various miss penalties, pipeline stalls, etc.) it is difficult
to predict.  However, $T_f$ can be precisely measured for a particular
application running on a particular system. 

The hardware part of the communication system on a processing element
(PE) is modeled as a network interface (NI) and an input link and
output link, as shown in Figure~\ref{fig:pe}.  (We use {\em processing
element} instead of the common term {\em node} to avoid confusion with
the nodes of finite element meshes.)  The unit of work for the
communication phase is the transfer of a {\em block} between the
network and the local memory system of the PE.  The time for this
transfer is called the {\em transfer time}. A block is a group of
words that move together from one PE to another. Blocks may be
fixed-sized or variable-sized. For example, a block might be a cache
line in a CC-NUMA system~\cite{dash92}, a bulk asynchronous data
transfer between two PEs' memories~\cite{t3esync}, or a page in
a software DSM system~\cite{treadmarkcompiler}.
\begin{figure}[tbp]
\centerline{
%\epsfxsize=3in
\epsfbox{pe.epsf}}
\caption{Processing element (PE) model (P: processor, M: memory system,
         NI: network interface).}
\label{fig:pe}
\end{figure}

The transfer time for a block $i$ of $l_i$ words is $T_l + l_iT_w$,
where $T_l$ is the constant {\em block latency}, $T_w$ is the constant
marginal cost of transferring each additional word, and $T_w^{-1}$ is
the {\em burst bandwidth}. Notice that the block latency does not
include any delay through a network. We are only modeling the overhead
of transferring data between the network and local memory, which reflects
previous findings that most of the cost of communication on modern
systems is incurred at the individual PEs~\cite{ISCA-MEMORY}.

If a PE transfers $B$ blocks with a total volume of $C = \sum l_i$
communication words, then the running time of the communication phase
on that PE is at least
\begin{displaymath}
T_{comm} = BT_l + CT_w.
\end{displaymath}
Assume for the moment that $F$, $C$, and $B$ are the same for each
PE. Then the running time for the SMVP is
\begin{displaymath}
T_{mv} = T_{comp} + T_{comm} = FT_f + BT_l + CT_w
\end{displaymath}
and each word is transferred in an average time of
\begin{displaymath}
T_c = \frac{T_{comm}}{C} = \frac{B}{C}T_l + T_w.
\end{displaymath}
The rate $T_c^{-1}$ is the {\em sustained bandwidth} during the
communication phase. The sustained bandwidth includes all block
transfers in and out of a PE, as well as all software and hardware
overheads. We will use sustained bandwidth as the performance metric
for the communication phase.

Modern graph partitioners do an excellent job of balancing computation
($F$), and an excellent job of minimizing the global volume of
communication, but they do less well in balancing the amount of
communication ($C$ and $B$) on each PE~\cite{meshsuite96}.  So let us
relax our assumption that $C$ and $B$ are the same on each PE.  As $C$
and $B$ vary among PEs, we assume that the length of the
communication phase is determined by the PE for which
$T_{comm}$ is largest.  As a simplifying assumption, we assume that
the PE with the largest communication volume in words (denoted
$C_{max}$) is the PE that transfers the greatest number of
blocks (denoted $B_{max}$); clearly, this {\em bottleneck PE} has the
largest $T_{comm}$ and its sustained bandwidth is given by
\begin{equation}
T_c = \frac{B_{max}}{C_{max}}T_l + T_w 
\label{eqn:model1}
\end{equation}

How reasonable is our assumption that the PE that transfers the most
words also transfers the most blocks?  Intuitively, we would expect the two
measures to be correlated, and our assumption is conservative---it never
underestimates the sustained bandwidth we require---although it may
overestimate it.  Regardless of the tradeoff between burst bandwidth
and latency, one can show (with some difficulty) that
we will never overestimate the sustained bandwidth by more than
\begin{equation}
\beta = 
1 + \min_{i\ \rm{is\ a\ PE}} \left[
\max \left\{
\frac{(B_{max} - B_i) C_{max}}{B_{max} C_i},
\frac{(C_{max} - C_i) B_{max}}{B_i C_{max}}
\right\} \right].
\label{eqn:beta}
\nonumber
\end{equation}
The metric $\beta$ is equal to one when $C_{max}$ and $B_{max}$
pertain to the same PE, and larger otherwise (but never larger than
two).\footnote{The final version of the paper will include a proof as
an appendix.}  In our sf meshes, partitioned among various numbers
of processors, we have never seen $\beta$ exceed $1.15$, and it is
often unity; hence, we are confident that our model approximates
reality well.

\section{Model validation on Cray T3D}
\label{sec:validation}
In this section we outline a simple method for estimating the
system-specific parameters $T_l$ and $T_w$, and then we use the method
to derive the parameters for a Cray T3D using MPI to handle
communication.  Using the parallel SMVP in a real
application (the sf2 Quake code), we scale the size of each message by
some factor $x$ and then measure and plot the elapsed time $y$ of the
communication phase as a function of $x$.  (Note that changing the
message sizes in this manner does not necessarily yield a working
application.  However, by measuring the change in communication time
as a function of changing message length, we can discriminate between
bandwidth-related and latency-related overheads.)  The parameter $T_l$ is
determined from the $y$-intercept of the resulting line, and $T_w$ is
determined from its slope. These parameters values can then be substituted
into Equation (\ref{eqn:model1}) to predict the sustained bandwidth
(and thus the total running time) of the communication phase of an SMVP
operation.

Figure~\ref{fig:mvperf} lists the absolute running times for the
computation and communication phases of various sf2 SMVP
operations. Notice that for the T3D, the communication phase
represents less than 10\% of the total running time. However, as we
shall see, this percentage becomes significantly larger if we replace
our multicomputer with a network of workstations, and hence incur a
much larger block latency $T_l$.

\begin{figure}[htb]
\begin{center}
\small
\begin{tabular}{|l||r|r|r|}
\hline 
\multicolumn{1}{|c|}{} & 
\multicolumn{1}{c|}{16 PEs} & 
\multicolumn{1}{c|}{32 PEs} &
\multicolumn{1}{c|}{64 PEs}\\
\hline
\hline
comp (ms)	&195.8	&99.3	&51.4\\
comm (ms)	&4.9	&4.5 	&4.1\\
total (ms)	&17.9	&109.3	&59.2\\
\hline
\end{tabular}
\end{center}
\caption{Running time of sf2 matrix-vector product on Cray T3D (ms)}
\label{fig:mvperf}
\end{figure}

Figure~\ref{fig:comm.scaling} graphs the measured running times on T3D
of the sf2 SMVP communication phase, as a function of different
message size scaling factors. The numeric values are given in
Figure~\ref{fig:comm.scaling.table}.
\begin{figure}[htbp]
\centerline{
%\epsfxsize=4in
\epsfbox{comm.scaling.epsf}
}
\caption{Running time for the communication phase of the sf2 SMVP on 
the Cray T3D}
\label{fig:comm.scaling}
\end{figure}

\begin{figure}[htbp]
\begin{center}
\begin{tabular}{|c||r|r|r|}
\hline 
\multicolumn{1}{|c||}{}&\multicolumn{3}{c|}{Communication times (us)}\\
\cline{2-4}
\multicolumn{1}{|c||}{scaling factor} & 
\multicolumn{1}{c|}{16 PEs}& 
\multicolumn{1}{c|}{32 PEs}&
\multicolumn{1}{c|}{64 PEs}\\
\hline
\hline
$0$		&462	&577	&848	\\
$1/4$		&1606	&1592	&1699	\\
$1/2$		&2789	&2517	&2255	\\
$1/\sqrt{2}$	&3510	&3302	&3011	\\
$1$		&4939	&4524	&4092	\\
$\sqrt{2}$	&6263	&5418	&5454	\\
$2$		&8499	&7776	&7098	\\
$3$		&13251	&11513	&10535	\\	
$4$		&16152	&14231	&13003	\\
\hline
\end{tabular}
\end{center}
\caption{Raw data for Figure~\protect{\ref{fig:comm.scaling}}}
\label{fig:comm.scaling.table}
\end{figure}

Each curve in Figure~\ref{fig:comm.scaling} corresponds to the running
time of the sf2 application as a function of the scaling factor.  Let
$x$ be a scaling factor and let $y(x)$ be the running time as a
function of $x$. There are several interesting aspects of
Figure~\ref{fig:comm.scaling}.  First, we see that the communication
times for different numbers of PEs scale nearly linearly, which
indicates that network congestion is not an issue. If it were, then we
would see concave curves rather than nearly straight lines.

Second, given some curve $y$ and some scaling factor $x$,
$y(0)$ is the running time associated with the block latency $T_l$ and
$y(x)-y(0)$ is the running time associated with the marginal cost $T_w$
of transferring an additional word. In particular,
\begin{equation}
T_l = \frac{y(0)}{B_{max}}
\label{eqn:tl}
\end{equation}
and 
\begin{equation}
T_w = \frac{y(x) - y(0)}{xC_{max}}
\label{eqn:tw}
\end{equation}
For a particular $x$, $T_w$ is directly proportional to the slope of
$y$ evaluated at $x$. 

Both $T_l$ and $T_w$ can be computed quite easily once we know the
application-specific parameters $C_{max}$ and $B_{max}$.  These
parameters for the sf2 SMVP, partitioned for execution on 16, 32, and
64 PEs, are summarized in Figure~\ref{fig:comm.parameters} for a range
of scaling factors between zero and four.  The parameters are
determined directly from properties of the partitioned meshes that
induce the SMVP operations~\cite{meshsuite96}.
\begin{figure}[htbp]
\begin{center}
\begin{tabular}{|c||r|r||r|r||r|r|}
\hline 
\multicolumn{1}{|c||}{} & 
\multicolumn{2}{c||}{16 PEs}& 
\multicolumn{2}{c||}{32 PEs}&
\multicolumn{2}{c|}{64 PEs}\\
\cline{2-7}
\multicolumn{1}{|c||}{Scaling factor} & 
\multicolumn{1}{c|}{$B_{max}$}&\multicolumn{1}{c||}{$C_{max}$}&
\multicolumn{1}{c|}{$B_{max}$}&\multicolumn{1}{c||}{$C_{max}$}&
\multicolumn{1}{c|}{$B_{max}$}&\multicolumn{1}{c|}{$C_{max}$}\\
\hline
\hline
$0$		&16	&0	&26	&0	&36	&0	\\
$1/4$		&16	&7120	&26	&6004	&36	&5130	\\
$1/2$		&16	&14241	&26	&12009	&36	&10260	\\
$1/\sqrt{2}$	&16	&19937	&26	&16812	&36	&14363	\\
$1$		&16	&28482	&26	&24018	&36	&20520	\\
$\sqrt{2}$	&16	&39874	&26	&33625	&36	&28727	\\
$2$		&16	&56964	&26	&48036	&36	&41040	\\
$3$		&16	&85446	&26	&72054	&36	&61560	\\
$4$		&16	&113928	&26	&96072	&36	&82080	\\
\hline
\end{tabular}
\end{center}
\caption{Application parameters for the sf2 SMVP}
\label{fig:comm.parameters}
\end{figure}

We can use Equations~(\ref{eqn:tl}) and (\ref{eqn:tw}), and the
application-specific parameters in Figure~\ref{fig:comm.parameters} to
arrive at reasonable estimates for $T_l$ and $T_w$ for the T3D. The
average $T_l$ thus derived for 16, 32, and 64 PEs is roughly 25 ns,
and we use that as an estimate of the block latency on the
T3D. Similarly, the average $T_w$ for 16, 32, and 64 PEs is roughly
160 ns, and we use that as the estimate for the marginal cost of
transferring a block word. While not exact, they give excellent
predictions of actual running time for the communication phase of the
sf2 SMVP operation, as shown in Figure~\ref{fig:comm.validate}.
\begin{figure}[htbp]
\centerline{
%\epsfxsize=4in
\epsfbox{comm.validate.epsf}
}
\caption{Predicted vs actual communication performance for the sf2 SMVP
($T_l = 25 us$, $T_w = 160 ns$)}
\label{fig:comm.validate}
\end{figure}

\section{Implications for Graph Partitioners}
\label{sec:implications}

Modern graph partitioners typically try to minimize global
communication volume (i.e., the sum of the communication volumes, $C$,
on each PE) subject to the constraint that each PE gets an equal
amount of computational work. While this is a convenient metric for
algorithm designers, it is not clear how closely this metric is
related to the actual running time of parallel programs. This is
mainly due to the facts that (1) it ignores the number of neighbors
that each PE must communicate with, and (2) it does not attempt to
someow balance communication across PEs.

The qualitative significance of (1) and (2) for real applications can
be be seen quite clearly by using the model in (\ref{eqn:model1}) to
separate the relative impact of block latency ($T_l$) and burst
bandwidth ($T_w^{-1}$) on the running time.  Consider
Figure~\ref{fig:comm.imbalance}, which shows the time due
to block latency ($BT_l$) versus the time due to burst bandwidth
($CT_w$) for the sf2 SMVP operation partitioned for execution on 32
PEs.
\begin{figure}[htbp]
\centerline{
%\epsfxsize=4in
\epsfbox{comm.imbalance.epsf}
}
\caption{Impact of block latency and burst bandwidth on communication
time for sf2 running on 32 PEs ($T_w = 160 ns$)} 
\label{fig:comm.imbalance}
\end{figure}
Each point in Figure~\ref{fig:comm.imbalance} represents the latency
and bandwidth components ($BT_l, CT_w$) of the total communication
time $BT_l + C_Tw$ for a PE involved in the SMVP operation. There are
three clusters of points, with each cluster corresponding to a
different assumption about the block latency $T_l$.  (1) $T_l = 25
us$: This represents the low-latency case, and is the measured block
latency for the Cray T3D using MPI.  (2) $T_l = 200 us$: This
represents an average-latency case, such as we might find on a
high-overhead message passing system like some PVM implementations.
(3) $T_l = 1000 us$: This represents the high-latency case, such as we
might find on a network of workstations where a message passing
operation requires a supervisor call. Notice for the low-latency case
that the absolute communication time is a small fraction of the total
SMVP time. However, for the medium and high latency cases the
communication time is a significant fraction of the overall SMVP time,
and thus it is important to reduce this time.

Since the slowest PE determines the running time, we would like to
minimize the maximum running time on any PE so that the the
communication time is evenly balanced across all the PEs. If the
communication phase were perfectly balanced, then all the points in a
cluster in Figure~\ref{fig:comm.imbalance} would fall on a 45-degree
diagonal extending from the top left corner to the bottom right corner
of the graph. In practice, however, we see that the communication
times are not balanced across the PEs and that both the latency and
bandwidth components of the individual PE running times can vary
significantly.  The point is that by paying attention only to the
global communication volume, the partitioner produces an SMVP
communication phase with unbalanced communication.

Another interesting aspect of Figure~\ref{fig:comm.imbalance} is the
degree to which different assumptions on block latency can cause the
communication phase of the same application to be either dominated by
the block latency component, dominated by the burst bandwidth
component, or balanced between the two. Applications that are
latency-dominated lie in the lower triangle, those that are
bandwidth-dominated lie in the upper triangle, and those that are
balanced between the two lie along the diagonal.  For the sf2 SMVP, we
see that the low-latency case ($T_l = 25 ns$) is bandwidth-dominated,
that the high-latency case ($T_l = 1000 ns$) is latency-dominated, and
that the medium-latency case is influenced equally by both latency and
bandwidth.  Thus, latency can be a significant component
of communication performance for real applications. And by using
global communication volume as the only optimization metric, the
partitioner is ignoring the impact of an important component in the
running time of the communication phase.

The main point is that modern partitioners are pursuing an
optimization goal (minimizing the global communication volume) that
does not necessarily minimize communication time because it ignores
the impact on communication time of both unbalanced communication and
block latencies. Global communication volume is certainly a convenient
optimization criterion because it is a single measure and it corresponds
directly to nice mathematical properties like separator sizes.
Unfortunately, it doesn't reflect the reality that block latency and
communication balance matter. 

Instead of minimizing global communication volume, we recommend an
approach, based on the results in this paper, where partitioners
directly try to minimize total communication time. We have seen from
the model in (\ref{eqn:model1}) that running time can be accurately
predicted by the metric $B_{max}T_l + C_{max}T_w$. The partitioner can
compute the application parameters $B_{max}$ and $C_{max}$ directly
from the graph as it is being partitioned, and the system parameters
can be derived using the techniques outlined in
Section~\ref{sec:validation} and then passed as run-time arguments to
the partitioner.  Again, we emphasize that we are not proposing a new
partitioning algorithm. Instead we are suggesting a system-dependent
objective function that future partitioners should seek to optimize.

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

We have developed a simple performance model for sparse matrix-vector
product (SMVP) operations such as those required by unstructured
finite element applications.  The parameters that
characterize the block latency and burst bandwidth of the underlying
system can be measured in a straightforward way and used to
accurately predict the running time of the communication phase of the
SMVP operations from a real application.

Experience with the model on real applications has provided us with
the following insight into the limitations of current graph
partitioning technology: by attempting to minimize only global
communication volume, modern graph partitioners ignore the impact of
block latency and overall communication imbalance. Partitioners
will be more effective if they can reduce the maximum degree
of adjacency of the subgraphs (especially for networks of workstations),
and balance the communication volume among processors.  We propose a
attempting to directly minimize running
time on a particular system, based on some simple application-specific
and system-specific parameters.

We have been able to quantify the extent to which a real application
can be both dominated and unbalanced by block latency; our approach attacks
these problems directly. However, we have not measured
how much such a direct approach can improve running time,
compared with the existing approach of minimizing communication volume.
We intend to pursue these questions, and encourage other researchers to
do the same.

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

\end{document}
