\newcommand{\comment}[1]{}

\documentstyle[11pt,fullpage,twoside,./cmu-titlepage,times,./epsf,./fancyheadings]{article}

\setcounter{topnumber}{4}
\renewcommand{\topfraction}{1}
\setcounter{bottomnumber}{4}
\renewcommand{\bottomfraction}{1}
\setcounter{totalnumber}{6}
\renewcommand{\textfraction}{0}
\renewcommand{\floatpagefraction}{0.8}

% \renewcommand{\baselinestretch}{0.93}

% \setlength{\topsep}{0pt}
% \setlength{\partopsep}{0pt}
% \setlength{\itemsep}{0pt}

%\def\headrule{{\if@fancyplain\headrulewidth\plainheadrulewidth\fi
%  \hrule\@height\headrulewidth\@width\headwidth}}

\raggedbottom
\pagestyle{plain}

\title{
Architectural Implications of a Family of Irregular Applications
}

\author{David O'Hallaron, Jonathan Richard Shewchuk, and Thomas Gross}

%\date{}

\trnumber{CMU-CS-97-189}

\support{ 
Effort sponsored in part by the Advanced Research Projects Agency and
Rome Laboratory, Air Force Materiel Command, USAF, under agreement
number F30602-96-1-0287, in part by the National Science Foundation
under Grant CMS-9318163, and in part by a grant from the Intel
Corporation. The U.S. Government is authorized to reproduce and
distribute reprints for Governmental purposes notwithstanding any
copyright annotation thereon.  The views and conclusions contained
herein are those of the authors and should not be interpreted as
necessarily representing the official policies or endorsements, either
expressed or implied, of the Advanced Research Projects Agency, Rome
Laboratory, or the U.S. Government. Authors' email addresses: 
\{droh@cs.cmu.edu
}

\keywords{computer architecture, 
finite element, parallel computing, sparse matrix}

\abstract{\small 
Irregular applications based on sparse matrices are at the core of
many important scientific computations. Since the importance of such
applications is likely to increase in the future, high-performance
parallel and distributed systems must provide adequate support for
such applications.  We characterize a family of irregular scientific
applications and derive the demands they will place on the
communication systems of future parallel systems. Running time of
these applications is dominated by thousands of sparse matrix vector
product (SMVP) operations. Using some simple performance models of the
SMVP, we investigate requirements for bisection bandwidth, sustained
bandwidth on each PE, burst bandwidth during block transfers, and
block latencies for PEs under different assumption about sustained
computational throughput.  The general conclusions are that bisection
bandwidth is not important for these applications, that PEs will need
300 MBytes/sec and 600 MBytes/sec of sustained and burst bandwidth
respectively to run irregular codes at high efficiency. Further, since
the block transfers between PEs tend to be small even for large
applications, block latency costs cannot be amortized by transferring
large messages. Other techniques for reducing observed latency must be
pursued.
}

\begin{document}
\bibliographystyle{acm}

\maketitle

% Print page numbers for table of contents using small roman numbers.
\renewcommand{\thepage}{\roman{page}}
\tableofcontents
\newpage
%\listoffigures
%\newpage

\setlength{\parskip}{\medskipamount}

\section*{About this Report}

This report is an extended version of a paper appearing in the Fourth
International Symposium on High Performance Computer Architecture, Las
Vegas, NV, Feb, 1998.

The Spark98 suite~\cite{spark98}, a collection of 10 portable
sequential and parallel SMVP kernels written in C and based on the
sf10 and sf5 meshes described in this paper, is available at
www.cs.cmu.edu/$\sim$quake/spark98.html.  The complete set of
partitioned San Fernando meshes is available from
www.cs.cmu.edu/$\sim$quake/meshsuite.html.

\cleardoublepage

% Revert to arabic numbers, and also reset counter to page 1.
\renewcommand{\thepage}{\arabic{page}}
\setcounter{page}{1}

%\topmargin -14.4pt
\pagestyle{fancy}
\addtolength{\headheight}{2.4pt}
\renewcommand{\sectionmark}[1]{\markboth{}{#1}}
\renewcommand{\subsectionmark}[1]{}
\lhead[\rm\thepage]{}
\chead[{\sl O'Hallaron, Shewchuk, and Gross}]
{{\sl \rightmark}}
\rhead[]{\rm\thepage}
\cfoot[]{}
\thispagestyle{empty}

\renewcommand{\thefootnote}{\arabic{footnote}}

\section{Introduction}
\label{sec:intro}
The peak performance of the commodity processors in modern parallel
systems is doubling every couple of years.  Clearly, communication
performance will have to improve as processor performance increases,
but the question is how much and in what respects.  It is crucial to
understand communication requirements because some parts of the
communication system will not be commodity, and are therefore
expensive.

In general it is important to understand the communication
requirements of real applications~\cite{iscanbody}, and these
requirements are especially difficult to characterize for the important
class of irregular scientific applications that manipulate sparse
matrices.  Such applications typically model processes in nature
like the flow of air over a wing, or the distribution of stresses and
strains in a bridge.  Characterizations of 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 models the performance of a family of large irregular
applications and uses the models to get some insight into the
following question: As microprocessors in parallel systems continue to
get faster, how fast must the communication system be in order to keep
up? Of course, other classes of applications will have different
communication requirements, and similar studies should be performed
for these as well. Nonetheless, the results in this paper provide very
clear answers for an application class that is a significant user of
high performance computer cycles.

%distinguish between quake -- the system and sf* -- specific simulations 
Specifically, the paper models a family of applications (called the
{\em Quake applications}) that simulate the motion of the ground
during strong earthquakes. The specific input set models the San
Fernando Valley of Southern California. The Quake applications,
described in Section~\ref{sec:quake}, were developed over a period of
three years by civil engineers at Carnegie Mellon as part of an
ongoing project to model earthquakes in the entire Los Angeles
basin~\cite{superquake96}.

A nice property of the Quake applications is that 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 Quake applications
that relates computation and communication, we derive requirements for
(1) bisection bandwidth, (2) sustained per-PE communication bandwidth,
and (3) block transfer latency and bandwidth under various assumptions
about efficiency and local MFLOPS rates. The models are derived in
Section~\ref{sec:model} and validated in Section~\ref{sec:validation}
on the Cray T3E. 

In general, we find that bisection bandwidth is not a critical issue;
overall communication performance is determined by the sustained
bandwidth at each PE. For simulations like the Quake applications with
flop/communication word ratios as low as 50:1 and running on PEs
that can sustain a local computation throughput of 200 MFLOPS on
sparse matrix computations, each PE will need to sustain a
communication rate of about 300 MBytes/sec and burst rates 
of 600 MBytes/sec on each PE to run with good efficiency.

Another key conclusion is that the maximum size of the individual
transfer units in irregular applications tends to be small.  Even if
transfer units and bandwidths are made as large as possible, latency
remains a significant component of the overall communication
performance. These results make clear the importance of decreasing
observed communication latency, either through thinner interfaces or
latency hiding techniques. These conclusions are described in more
detail in Section~\ref{sec:requirements}.

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

The Quake applications were developed to study the seismic
response of the San Fernando Valley during strong
earthquakes~\cite{superquake96}.  The goal is not to predict the
occurrence of earthquakes, but rather to understand the following kinds
of questions when an earthquake does occur: Which areas will be worst
stricken?  Which seismic frequencies will be amplified most by the
soil?

There are four Quake applications, denoted sf10, sf5, sf2, and sf1.
The ``sf'' is an abbreviation for San Fernando, and the number
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 Quake application consists of (1) a three-dimensional
unstructured finite element model of the ground beneath San Fernando,
and (2) a parallel
finite element simulation that propagates seismic waves through the
model for 60 seconds of simulated time.

\subsection{Finite element model}

Each Quake application employs a three-dimensional unstructured mesh
composed of thousands or millions of tetrahedra (i.e., pyramids with
triangular bases). The model for sf10 is shown in
Figure~\ref{fig:sf10mesh-all}.  The model corresponds to a volume of
earth roughly 50 km x 50 km x 10 km.  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 smaller elements.

Figure~\ref{fig:general} lists the topological 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 for each of three dimensions.  As a general
rule, each node in the mesh, a simulation uses about 1 KBytes of
memory at runtime to accommodate the storage of several vectors and
sparse matrices. 
\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}

\subsection{Parallel finite element simulation}

During each of the 6000 time steps, a Quake finite element 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$ (here $n$ is the number of nodes),
and ${\bf K}$ is a sparse $3n \times 3n$ {\em stiffness matrix}.
The vectors ${\bf x}$ and ${\bf y}$ have length $3n$ because each
represents three degrees of freedom---$x$, $y$, and $z$
displacements---for each node of the mesh.
${\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 that are 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{sfprops96}.

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$ processors, Archimedes
partitions the mesh into $p$ disjoint sets of elements. Each set is
called a {\em subdomain} and is assigned to some processor. 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{GiMiTe95,barnard92,chaco95}.

\subsection{Parallel SMVP}

The running time of the Quake applications is dominated by the
execution of SMVP operations, which consume over 80\% of the total
running time in the sequential case. Furthermore, the SMVP operations
are the only operations (besides I/O) that require interprocessor
communication. So even though the Quake applications are nontrivial,
their performance is characterized by the performance of
the SMVP operations. Therefore we can abstract away most of the
complexity and focus on the performance of a simple well-defined
parallel SMVP kernel.

To compute the global SMVP operation ${\bf y = Kx}$ on a set of processors,
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 processors
induced by the partition of elements among processors.  If a node $i$
resides in several processors (because $i$ is a vertex of several
elements mapped to different processors), the values $x_i$ and $y_i$
are replicated on those processors.  The matrix ${\bf K}$ is
distributed so that $K_{ij}$ resides on any processor 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 processors.  Xs represent nonzero $3 \times 3$
         submatrices.  Note that some nodes are shared by both processors,
         as are some stiffness matrix entries (corresponding to shared edges).}
\label{fig:cutmatrix} 
\end{figure}

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

\section{SMVP performance model}
\label{sec:model}

Each global SMVP consists of a simple sequence, executed
simultaneously on each processor: a local SMVP (the computation phase)
followed by an operation in which each processor 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.\footnote{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 application do not do this, so the models
in this section will follow suit. By not modeling any
overlap, we obtain conservative bandwidth and latency estimates,
% add a few reasons why we don't do it ....
and avoid possibly de-optimizing the program by complicating the runtime system
.}
Further, we assume (1) that during the communication phase, the summing
operations are overlapped with the exchanges, and (2) that the time to perform
the exchanges determines the running
time of the communication phase.  These assumptions are reasonable, 
as the summing
operations can be handled by the processor and the exchanges by a
network interface.  

Given these assumptions, the running time for each global SMVP is
\begin{displaymath}
T_{smvp} = T_{comp} + T_{comm},
\end{displaymath}
where $T_{comp}$ is the time required for all processors to finish their
computation phases and $T_{comm}$ is the time required for all processors
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
\multicolumn{2}{|c|}{Parameters for Equation~(\ref{eqn:macromodel})}\\
\hline
$F$ &flops per PE per SMVP\\
$C_{max}$ &maximum communication words per PE per SMVP\\
$E$ & target efficiency\\
$T_f$ & amortized time per flop (inverse of sustained flops)\\
\hline
\hline
\multicolumn{2}{|c|}{Parameters for Equation~(\ref{eqn:micromodel})}\\
\hline
$B_{max}$ &maximum communication blocks per PE per SMVP\\ 
$C_{max}$ &maximum communication words per PE per SMVP\\
$T_l$ &time per communication block (block latency)\\
$T_w$ &time per additional block word (inverse of burst bandwidth)\\
\hline
\hline
\multicolumn{2}{|c|}{Computed quantities}\\
\hline
$T_{smvp}$ &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\\
$T_c$ & amortized time per communication word
        (inverse of sustained bandwidth)\\
$\beta$ & upper bound on relative overestimate of communication time\\
\hline
\end{tabular}
\end{center}
\caption{Summary of symbols used in this paper.}
\label{fig:symbols}
\end{figure}

\subsection{Model of computation phase}

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 the SMVP operation because the
(minimum) number of flops each processor must perform can be counted
exactly and is independent of the compiler and its optimization
level. If a PE's 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 $FT_f$. 

Modern graph partitioners do an excellent job of distributing
computation evenly across PEs~\cite{sfprops96}, so we will assume that
each PE performs $F = 2m$ flops and thus the running time for the
computation phase 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. For example, measurements
of the local SMVPs from the Quake applications show a steady $T_f =
30$ ns for the Cray T3D (150 MHz Alpha 21064, cc -O3) and $T_f = 14$
ns for the Cray T3E (300 MHz Alpha 21164, cc -O3).

\subsection{High level model of communication phase}

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.)  
\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 running time of the communication phase is given by
\begin{displaymath}
T_{comm} = C_{max} T_c
\end{displaymath}
where $C_{max}$ is the maximum number of words communicated by any PE
during the communication phase and $T_c$ is the average time per
communication word.  The rate $T_c^{-1}$ is the {\em sustained
bandwidth} across the links of the network interface during the
communication phase. It includes all data transferred in and out of a
PE, as well as all software and hardware overheads.

We define the {\em efficiency} $E$ to be the proportion of SMVP time devoted 
to computation; that is, $E = T_{comp}/T_{smvp}$.
Given this definition,
$T_{comm}/T_{comp} = (C_{max} T_c)/(FT_f) = (1-E)/E$, and hence
\begin{equation}
T_c = \frac{F}{C_{max}} \left( \frac{1-E}{E} \right) T_f.
\label{eqn:macromodel}
\end{equation}
Equation~(\ref{eqn:macromodel}) describes at a high level the relation
between sustained computation bandwidth ($T_f^{-1}$) and sustained
communication bandwidth ($T_c^{-1}$).  The model is interesting
because it cleanly separates the various factors that relate
computation performance to communication performance: The {\em
computation/communication ratio} $F/C_{max}$ is a property of the
application and the partitioner; $T_f$ is a property of the processor
architecture and compiler; $E$ is a target efficiency imposed by the
user. This separation of factors is similar in spirit to the familiar CPI 
model for uniprocessor performance~\cite{hp95}.

In Section~\ref{sec:requirements} we use
Equation~(\ref{eqn:macromodel}) to investigate the following question:
Given a target efficiency $E$ for an SMVP running on PEs with local
computational rates of $T_f^{-1}$, what is the required communication
performance of the system? The answers to this question provide
insight, for this class of applications, into the sustained
performance we need from our communication systems as
processor speeds continue to increase.

\subsection{Low level model of communication phase}

Equation~(\ref{eqn:macromodel}) provides a high level model of sustained
performance in terms of streams of words being transferred between
memory and the network interface. However, in a lower level and more
realistic view of communication, the unit of work during 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 transfer
unit, 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
message in a message passing system~\cite{mpisuper}, a bulk
asynchronous data transfer between two PEs' memories~\cite{t3esync},
or a page in a software DSM system~\cite{treadmarks96,treadmarkcompiler}.

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 ({\em the inverse
burst bandwidth}), and $T_w^{-1}$ is the {\em burst bandwidth}. Notice
that the block latency does not include any delay through an
interconnection network. We are only modeling the overhead of
transferring data between the network interface and local memory
which reflects previous findings that most of the cost of
communication on modern systems is incurred at the individual
PEs~\cite{iscamemory}. In essence, we assume that the interconnection
network has infinite capacity and zero latency. In
Section~\ref{sec:validation} we argue empirically
that this is a reasonable
assumption for the SMVP running on tightly coupled systems.

Modern graph partitioners do an excellent job of balancing computation
($F$) and an excellent job of minimizing the global volume of
communication. Unfortunately they do less well in balancing the total
number of blocks $B_i$ and the total volume in words sent and received
by each PE $i$~\cite{sfprops96}.  As a simplifying assumption, we
assume that the PE that transfers the most words (denoted $C_{max}$)
is also the PE that transfers the most blocks (denoted
$B_{max}$). Clearly, such a PE would have the longest communication
phase, and thus
\begin{displaymath}
T_{comm} \leq B_{max}T_l + C_{max}T_w.
\end{displaymath}
Finally, since $T_c = T_{comm}/C_{max}$, we have
\begin{equation}
T_c \leq \frac{B_{max}}{C_{max}}T_l + T_w.
\label{eqn:micromodel}
\end{equation}
$T_l$ and $T_w$ are properties of the communication system, including
the architecture, hardware implementation, and software libraries.
For message-passing systems, $B_{max}$ is purely a property of the
application and the partitioner, but on shared-memory systems it is
also a property of the architecture. Thus,
Equation~(\ref{eqn:micromodel}) characterizes sustained communication
bandwidth during the communication phase of the SMVP in terms of basic
properties of the communication system and the application. In
Section~\ref{sec:requirements}, we will use this model to explore
bandwidth and latency tradeoffs in communication systems as the base
microprocessors continue to improve.

The models in Equations (\ref{eqn:macromodel}) and
(\ref{eqn:micromodel}) are similar in some ways and different in
others to the familiar LogP model~\cite{logp96}. LogP is a general
performance model for bulk-synchronous parallel (BSP)
computations~\cite{bulksync90}, where a program is viewed as a series
of {\em supersteps} separated by barrier synchronizations. During a
superstep, each PE performs local computation and transfers a limited
number of messages.  The models we have developed are for a restricted
class of BSP programs where each superstep (SMVP) consists of a
distinct computation phase followed by a distinct communication
phase. In general our models view performance at a lower level of
abstraction than LogP. The parameters $T_f$, $T_w$, $F$, $B_{max}$,
and $C_{max}$ have no counterparts in LogP. On the other hand, our
$T_l$ parameter is similar to the overhead parameter $o$ in LogP.

\subsection{Error bound}

In (\ref{eqn:micromodel}) we assume that the PE with the most words to
communicate also has the most blocks to communicate.  In practice this
assumption may not be true. The question is, by how much do we
overestimate $T_c$ (and thus $T_{comm}$) as a result?  For a given
SMVP instance, this question is easy to answer if we know the machine
parameters $T_l$ (block latency) and $T_w$ (inverse burst bandwidth).
However, we would like to establish a worst-case bound that applies to
all machines, regardless of their parameters.

The length of the communication phase is determined by the PE
that takes the longest; that is,
\[
T_{comm} = \max_{i\ \rm{is\ a\ PE}} (B_iT_l + C_iT_w).
\]
Our simplifying assumption is to calculate $T_{comm}$ as if there is some
PE $i$ such that $B_i = B_{max}$ and $C_i = C_{max}$.
Let $\beta$ be the largest factor by which this assumption might
cause our model to overestimate $T_{comm}$.
Figure~\ref{fig:beta} shows that $\beta$ is not much larger than
one for any of the Quake
applications.  Hence, (\ref{eqn:micromodel})  appears to be a
reasonable model.
The rest of this section describes how the numbers in
Figure~\ref{fig:beta} are obtained,
and may safely be skipped by the reader not interested in technical detail.

\begin{figure}[tbp]
\begin{center}
\begin{tabular}{|c||r|r|r|r|}
\hline 
\multicolumn{1}{|c||}{subdomains} & 
\multicolumn{1}{c|}{sf10} & 
\multicolumn{1}{c|}{sf5} &
\multicolumn{1}{c|}{sf2} &
\multicolumn{1}{c|}{sf1} \\ 
\hline
\hline
4	&1.00	&1.00	&1.00	&1.00\\
8	&1.00	&1.00	&1.00	&1.00\\
16	&1.09	&1.10	&1.07	&1.00\\
32	&1.01	&1.01	&1.15	&1.00\\
64	&1.03	&1.08	&1.11	&1.05\\
128	&1.03	&1.04	&1.04	&1.11\\
\hline
\end{tabular}
\end{center}
\caption{Measured error bounds $\beta$ for the Quake applications.}
\label{fig:beta}
\end{figure}

Because the sustained bandwidth we require is $T_c = T_{comm} / C_{max}$, our
assumption will cause us to overestimate $T_c$ by no more than
\begin{eqnarray*}
\beta & = & \frac{B_{max}T_l + C_{max}T_w}
                 {\max_i (B_iT_l + C_iT_w)} \\
      & = & \frac{B_{max} + \alpha C_{max}}
                 {\max_i (B_i + \alpha C_i)} \\
      & = & \min_i \frac{B_{max} + \alpha C_{max}}
                        {B_i + \alpha C_i},
\end{eqnarray*}
where $\alpha = T_w / T_l$ is a parameter representing the ratio between
the overhead per word and the overhead per block.

To ensure that our model is reasonably accurate for
any machine parameters, we wish to find the value of $\alpha$ that
yields the most pessimistic value of $\beta$.  For argument's sake,
suppose there are only three PEs, having $(B, C)$ parameters
of $(B_j, C_j)$, $(B_{max}, 0)$, and $(0, C_{max})$, respectively (the
latter two are the most conservative assumption possible).
Because of the last two PEs, the denominator
$\max_i (B_i + \alpha C_i)$ cannot be less than $B_{max}$,
nor can it be less than $\alpha C_{max}$; hence, $\beta$ cannot
be greater than two.

Now, consider the first PE $j$, and the function
\[
\beta_j(\alpha) = \frac{B_{max} + \alpha C_{max}}
                       {B_j + \alpha C_j}.
\]
It is apparent upon
inspection that the derivative $\rm{d}\beta_j / \rm{d}\alpha$ is
either everywhere positive, everywhere negative, or uniformly zero
over $\alpha > 0$; hence, there are no critical points, and
$\beta_j$ is maximized at an extreme value of $\alpha$.

However, the notion of ``extreme values of $\alpha$'' is not
straightforward, and does not necessarily correspond to the
values $\alpha = 0$ and $\alpha = \infty$.  Obviously,
the PE that maximizes the value of $B_i + \alpha C_i$
depends on $\alpha$.  As $\alpha$ is swept from zero to
infinity, the PE (called the {\em active PE} in
the terminology of optimization) that maximizes this expression may
change several times; each $\alpha$ value where such a change
occurs determines a candidate for $\beta$.

Returning to our three-PE example, we note that when
$\alpha = 0$, the active PE is $(B_{max}, 0)$ (because
it maximizes $B_j + \alpha C_j$), yielding $\beta = 1$, and
when $\alpha \rightarrow \infty$, the active PE becomes
$(0, C_{max})$, also yielding $\beta = 1$.  Hence, the value
of $\alpha$ that maximizes $\beta$ will occur between these
extremes.  As $\alpha$ increases from zero, the active PE
changes from $(B_{max}, 0)$ to $(B_i, C_i)$ when
$B_{max} = B_i + \alpha C_i$, and changes from
$(B_i, C_i)$ to $(0, C_{max})$ when $B_i + \alpha C_i = \alpha C_{max}$.
In each of these two cases, we solve for $\alpha$ and substitute
the result into $\beta$, yielding
\[
\beta\left(\alpha = \frac{B_{max} - B_i}{C_i}\right)
  = 1 + \frac{C_{max}(B_{max} - B_i)}{C_iB_{max}}
\]
and
\[
\beta\left(\alpha = \frac{B_i}{C_{max} - C_i}\right)
  = 1 + \frac{B_{max}(C_{max} - C_i)}{B_iC_{max}}.
\]
The larger of these two values is an upper bound on $\beta$ (over
all values of $\alpha$).

If there are more than three PEs, such a bound may be computed
for each PE.  Each such bound is valid, so we use the smallest.
Hence, our model does
not overestimate $T_c$ by more than
\begin{equation}
\beta = 1 + \min_{i\ \rm{is\ a\ PE}}
            \left[ \max \left\{
            \frac{C_{max}(B_{max} - B_i)}{C_iB_{max}},
            \frac{B_{max}(C_{max} - C_i)}{B_iC_{max}} \right\} \right].
\label{eqn:beta}
\end{equation}
In general, the bound in (\ref{eqn:beta}) is not tight.  The tightest
bound (when $\alpha$ is unknown) may be derived by computing the
convex hull of the set of points $(B_i, C_i)$ (where $i$ ranges over
all PEs), and checking each edge of the convex hull as a possible
maximizer of $\beta$.  (Each edge of the convex hull represents a
change in the active PE.)

The bound in (\ref{eqn:beta}) is equal to one when $C_{max}$ and
$B_{max}$ pertain to the same PE, and larger otherwise (but never
larger than two).

\section{Model validation}
\label{sec:validation}

The parameters in Equation~(\ref{eqn:macromodel}) are straightforward
to determine: $F$ and $C_{max}$ are static application properties,
$T_f$ is a constant property of the compiler and the processor, and
the efficiency $E$ is imposed by the user.  Similarly for
Equation~(\ref{eqn:micromodel}), the parameters $B_{max}$ and
$C_{max}$ are static application properties. However, the block
latency $T_l$ and inverse burst bandwidth $T_w$ are system properties
that require some care to estimate.  In this section we  summarize
the properties of the Quake SMVP instances, describe a method for
estimating $T_l$ and $T_w$, and then plug these various parameters
into our models to predict running time on a Cray T3E.

\subsection{Properties of the Quake applications}

Figure~\ref{fig:properties} summarizes the properties of the various
Quake SMVP instances.  (See O'Hallaron and Shewchuk~\cite{sfprops96}
for a complete
characterization of the applications.)  Note that the values of $B_{max}$
and $C_{max}$ in this table are always even, because each message from PE
$i$ to PE $j$ is matched by a message from $j$ to $i$ of equal length.
(The values of $C_{max}$ are also divisible by three, because there are
three degrees of freedom per node.)  The values of $B_{max}$ indicate the
degree of subdomain adjacency; for instance, if $B_{max} = 46$, then
each subdomain shares nodes with at most $23$ other subdomains, and
some PE must communicate with $23$ other PEs.  These values assume
that blocks may be arbitrarily large, and hence each PE sends at most
one block to each other PE.  On a  fine-grained shared memory machine, where
the unit of interchange between PEs may be as small as a cache line,
$B_{max}$ may be much larger.
\begin{figure}[tbp]
\begin{center}
\begin{tabular}{|c|l||r|r|r|r|}
\hline 
\multicolumn{2}{|c||}{subdomains} & 
\multicolumn{1}{c|}{sf10} & 
\multicolumn{1}{c|}{sf5} &
\multicolumn{1}{c|}{sf2} &
\multicolumn{1}{c|}{sf1} \\ 
\hline
\hline
	&$F$        		&453,924	&1,899,396	&24,640,110	&162,372,024\\
4	&$C_{max}$  		&2,352		&7,746		&55,338		&186,162\\
	&$B_{max}$		&6		&6		&6		&6	\\
	&$F/C_{max}$		&193		&245		&445		&872	\\
\hline
	&$F$        		&235,566	&970,740	&12,414,006	&81,602,442\\
8	&$C_{max}$  		&2,550		&7,080		&35,148		&151,764\\
	&$B_{max}$		&12		&12		&10		&14	\\
	&$F/C_{max}$		&92		&137		&353		&538\\
\hline
	&$F$        		&122,742	&496,872	&6,278,076	&41,116,374\\
16	&$C_{max}$  		&2,208		&5,292		&28,482		&119,280\\
	&$B_{max}$		&18		&20		&16		&18	\\
	&$F/C_{max}$		&56		&94		&220		&345\\
\hline
	&$F$        		&64,980		&257,004	&3,191,436	&20,740,734\\	
32	&$C_{max}$  		&2,172		&4,476		&24,018		&87,228	\\
	&$B_{max}$		&30		&30		&26		&26	\\
	&$F/C_{max}$		&30		&57		&133		&238\\
\hline
	&$F$        		&34,956		&134,424	&1,632,708	&10,511,586\\
64	&$C_{max}$  		&1,764		&4,296		&20,520		&73,062	\\
	&$B_{max}$		&38		&40		&36		&38	\\
	&$F/C_{max}$		&20		&31		&80		&144	\\
\hline
	&$F$        		&18,954		&70,956		&838,224	&5,332,806\\
128	&$C_{max}$  		&1,740		&3,360		&16,260		&51,048	\\
	&$B_{max}$		&62		&52		&50		&46	\\
	&$F/C_{max}$		&11		&21		&52		&104	\\
\hline
\end{tabular}
\end{center}
\caption{Quake SMVP properties.
($F$: flops per PE, $C_{max}$: maximum communication words per PE,
$B_{max}$: maximum communication blocks per PE, $F/C_{max}$:
computation/communication ratio.)}
\label{fig:properties}
\end{figure}

There are some interesting points to make about the numbers in
Figure~\ref{fig:properties}. First, conventional wisdom holds that
sparse codes like the SMVP are communication intensive.  However, this
is not always the case. As we see for sf2, which is a reasonably large
problem, the computation/communication ratios $F/C_{max}$ vary from
large (450 for sf2/4) to moderate (50 for sf2/128).\footnote{We use the
notation sf$x$/$y$ to denote an instance of application sf$x$ partitioned into $y$
subdomains.}  This common misconception about sparse codes is probably
due to the fact that researchers have not had the opportunity to run
large enough problems.

Second, while the computation/communication ratios are reasonably high
for large problems, as the problem sizes grow by a factor of ten, we
see that the computation/communication ratios grow only by a factor of
two.  This is not surprising; consider that a good partition of an
$n$-node 3D mesh will produce ${\cal O}(n^{2/3})$ shared nodes (for
the same reason that an ${\cal O}(n)$-node cube has a surface area of
${\cal O}(n^{2/3})$ nodes).  Hence, the computation/communication
ratio is ${\cal O}(n^{1/3})$, and a factor-of-ten increase in $n$ will
yield roughly a factor-of-two increase in that ratio.  The point is
that while large SMVPs indeed have reasonable computation/communication
ratios, these ratios do not increase quickly with increasing problem
size, as they do for cubic problems like dense matrix multiply. Thus,
we cannot rely on simply increasing the problem size to guarantee
good efficiency.

Third, even when blocks are as large as possible (i.e., each PE sends
at most one block to any other PE), the block sizes (whose mean is
approximately $C_{max}/B_{max}$) are surprisingly small. For example,
the largest mesh (sf1) running on 128 PEs has an average message size
of only 500 words~\cite{sfprops96}. Thus, we cannot rely on large
blocks to amortize high block latencies.

Finally, for sf1/128 each PE communicates with up to 20\% of the other
PEs. Thus, the Quake SMVP occupies an interesting middle ground
between difficult applications like the 2D FFT that require an all-to-all
communication, and simple applications like regular grid problems wherein
PEs communicate with a maximum of four neighbors.

\subsection{Estimating $T_l$ and $T_w$}

Here we describe the method we use to estimate the parameters $T_l$ and $T_w$.
Starting with the parallel SMVP from the sf2 Quake application, we
scale the size of each message by some factor $x$ and then measure and
plot the elapsed time $y(x)$ of the communication phase as a function
of $x$. Performing this step for a range of scaling factors yields a
curve. The parameter $T_l$ is determined from the $y$-intercept of
this curve, and $T_w$ is determined from its slope.\footnote{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
distinguish between bandwidth-related and latency-related overheads.}

Given some curve $y(x)$ 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(x)$ of
transferring an additional word. In particular,
\begin{equation}
T_l = \frac{y(0)}{B_{max}}
\label{eqn:tl}
\end{equation}
and 
\begin{displaymath}
T_w(x) = \frac{y(x) - y(0)}{xC_{max}}
\end{displaymath}
For a particular $x$, $T_w(x)$ is directly proportional to the
derivative of $y(x)$ evaluated at $x$. For simplicity, we will
estimate this derivative for each SMVP instance at $x=1$ using a
simple backward difference approximation~\cite{pdeintro85}:
\begin{equation}
T_w = \frac{y(1) - y(0)}{C_{max}}
\label{eqn:tw}
\end{equation}
Figure~\ref{fig:comm.scaling} graphs the communication time for
instances of the sf2 running on 8, 16, 32, and 64 PEs of a Cray T3E,
as a function of different message size scaling factors.  The
communication times scale nearly linearly with the size of the messages,
which indicates that network congestion is not an issue. If it were,
then we would see concave curves rather than nearly straight lines.
The approximation of Equation~(\ref{eqn:tw}) to the
derivative at $x=1$ is reasonable only because of this linear scaling. 

Another interesting aspect is that the volume of data transferred by
each PE does not strongly influence the running time until $x \geq 2$.
The reason for this is that even a large application like sf2
transfers fairly short blocks, even when those blocks are as large as
possible. Thus the block latency is especially influential for these
applications and cannot be amortized by long messages.
\begin{figure}[tbp]
\centerline{
\epsfxsize=3in
\epsfbox{t3e.comm.scale.epsf}
}
\caption{Scaling of the sf2 SMVP on the Cray T3E}
\label{fig:comm.scaling}
\end{figure}

Applying Equations~(\ref{eqn:tl}) and (\ref{eqn:tw}) to each of the
curves in Figure~\ref{fig:comm.scaling} yields the estimates for $T_l$
and $T_w$ shown in Figure~\ref{fig:tltw}.  
We use the averages of these values as our 
final estimates for block latency and inverse burst bandwidth
for the Cray T3E: $T_l = 22 ~\mu$s and $T_w = 55$ ns.
\begin{figure}[tbp]
\centerline{
\epsfxsize=3in
\epsfbox{t3e.tltw.epsf}
}
\caption{Measured values of $T_l$ and $T_w$ on the Cray T3E.
Note that $T_l$ is measured in different units (a thousand times larger)
than $T_w$.}
\label{fig:tltw}
\end{figure}

\subsection{Predicted vs actual performance}

Plugging the application properties from Figure~\ref{fig:properties}
and the estimates of $T_l = 22 ~\mu$s and $T_w = 55$ ns into
Equation~(\ref{eqn:micromodel}) yields the predicted running times 
shown in Figure~\ref{fig:comm.validate}. While
not exact, they provide reasonable predictions of actual running time.
\begin{figure}[tbp]
\centerline{
\epsfxsize=3in
\epsfbox{t3e.comm.validate.epsf}
}
\caption{Predicted vs actual time of the sf2 SMVP communication phase 
on the Cray T3E ($T_l = 22 ~\mu$s, $T_w = 55$ ns).}
\label{fig:comm.validate}
\end{figure}

While we have shown in this section that is possible to obtain
reasonable estimates of $T_l$ and $T_w$, these parameters are
functions of complex hardware and software systems, and thus are
really not constants, as we see clearly in Figure~\ref{fig:tltw}.
Characterizing these parameters can still be considered an open
problem for some application scenarios, 
although for our purposes, the constant estimates for $T_l$
and $T_w$  seem to be reasonable.

\section{Communication requirements for the SMVP}
\label{sec:requirements}

In this section, we use our performance models to address the
following question: As the base commodity microprocessors in parallel
systems get faster, what kind of performance will be needed from their
communication systems?

We will investigate this question for two hypothetical machines: A
machine that runs the local SMVP at a rate of 100 MFLOPS ($T_f = 10
$ ns), and another that runs at a rate of 200 MFLOPS ($T_f = 5 $ ns).
(All specifications are computed assuming 64-bit floating-point
words.)  Computational rates of 100 to 200 MFLOPS may seem low, but
they reflect the reality that the actual performance of irregular
applications is far less than peak rated performance, largely because
of irregular reference patterns and because the data structures are
too large to fit in cache. For example, a single PE of a
Cray T3E (300 MHz Alpha 21164, cc -O3) runs the local SMVP from the
Quake applications at approximately 70 MFLOPS ($T_f = 14$ ns), which
is only 10\% of the peak rated performance of 600 MFLOPS.

The sf2 SMVP will serve as a running example. It is a good example for
a number of reasons. First, it is clearly large enough to qualify as a
real problem by any standard (the sparse coefficient matrix is $1.14$M
$\times 1.14$M in size). Second, even though it is quite large, it has
a good range of computation/communication ratios ($F/C_{max}$), from a
high of 450:1 when partitioned into four subdomains, down to 50:1 for
128 subdomains. 

\subsection{Absolute performance}

Figure~\ref{fig:sf2tmv} shows the running time of the sf2 SMVP for the
two hypothetical machines.  From Figure~\ref{fig:sf2tmv} we see that
if efficiency is 50\%, then each sf2 SMVP runs for $0.5$ seconds on
four 100-MFLOP PEs, and an entire sf2 simulation runs for about 2.5
hours. (Recall that each Quake application simulates 60 seconds of
shaking using 6,000 timesteps, with three SMVPs per timestep, for a
total of 18,000 SMVPs per simulation.)  If efficiency is 90\%, then
each sf2 SMVP runs for about 10 ms on 128 100-MFLOP PEs, and the
entire simulation runs for about 3 minutes. The point is that
different combinations of efficiency and numbers of PEs can result in
orders of magnitude differences in the total running times of 
simulations.

\begin{figure*}[tbp]
\begin{center}
\begin{tabular}{c c}
\epsfxsize=3in
\epsfbox{sf2.tmv.tf10.epsf}
&
\epsfxsize=3in
\epsfbox{sf2.tmv.tf5.epsf}\\
&\\
(a) $T_f = 10$ ns (100 MFLOPS). & 
(b) $T_f = 5$ ns (200 MFLOPS).\\
\end{tabular}
\end{center}
\caption{Running time ($T_{smvp}$) of the sf2 SMVP.}
\label{fig:sf2tmv}
\end{figure*}


\subsection{Bisection bandwidth}

Bisection bandwidth is a popular measure of communication system
performance, but it proves to be unimportant for Quake SMVPs.  To
compute the bisection bandwidth requirements for a Quake SMVP running
on $p$ PEs, we create a symmetric $p \times p$ matrix $m$ such that
$m_{ij}$ is the number of 64-bit words transferred from PE $i$ to PE
$j$. If we assume that PEs $0, \ldots, p/2-1$ are on one side of the
bisection and PEs $p/2, \ldots, p-1$ are on the other side, then
\begin{displaymath}
V = 2 \sum_{i=0}^{p/2-1} \sum_{j=p/2}^{p-1} m_{ij}
\end{displaymath}
words cross the bisection during the communication phase, and the
sustained bisection bandwidth is $V/(C_{max}T_c)$.
Figure~\ref{fig:sf2bi} shows the required bisection bandwidths given
different assumptions of PE performance and overall efficiency.  The
worst case of 700 MBytes/sec ($E = 0.9$ and $T_f^{-1} = 200$ MFLOPS)
is quite modest, on the order of the bandwidth on a couple of links in
a modern system. Thus, it does not appear that bisection bandwidth
will be an issue for SMVPs as local PE performance increases.

\begin{figure*}[tbp]
\begin{center}
\begin{tabular}{c c}
\epsfxsize=3in
\epsfbox{sf2.bi.dp.tf10.epsf}
&
\epsfxsize=3in
\epsfbox{sf2.bi.dp.tf5.epsf}\\
(a) $T_f = 10$ ns (100 MFLOPS). &
(b) $T_f = 5$ ns (200 MFLOPS).\\
\end{tabular}
\end{center}
\caption{Sustained bisection bandwidths required for sf2.}
\label{fig:sf2bi}
\end{figure*}


\subsection{Sustained PE bandwidth}

Figure~\ref{fig:sf2tc} plots the sustained bandwidth for each PE
($T_c^{-1}$) that is required for the sf2 SMVP under different
assumptions of PE performance and overall efficiency.  The required
bandwidths are computed using model (\ref{eqn:macromodel}) and the
properties from Figure~\ref{fig:properties}. Notice that the required
bisection bandwidths are only a factor of two higher than the required
sustained PE bandwidths. Thus, bisection bandwidth is not an issue;
sustained PE bandwidth is the key communication performance measure.

\begin{figure*}[tbp]
\begin{center}
\begin{tabular}{c c}
\epsfxsize=3in
\epsfbox{sf2.tc.dp.tf10.epsf}
&
\epsfxsize=3in
\epsfbox{sf2.tc.dp.tf5.epsf}\\
(a) $T_f = 10$ ns (100 MFLOPS). & (b) $T_f = 5$ ns (200 MFLOPS).\\
\end{tabular}
\end{center}
\caption{Sustained PE bandwidth ($T_c^{-1}$) required for sf2.}
\label{fig:sf2tc}
\end{figure*}

On a system with 100-MFLOP PEs, maintaining a sustained rate of
120 MBytes/sec per PE during the communication phase is sufficient to
run all instances of the sf2 SMVP at 90\% efficiency.  On 
systems with 200-MFLOP PEs, a sustained PE bandwidth of about 300
MBytes/sec will be required to run all of the sf2 SMVPs at 90\%
efficiency. 

A sustained PE bandwidth of 300 MBytes/sec may seem like an easy
performance target until one remembers that it includes all overheads,
including software overheads, local strided read and write copies, and
other latencies; we are not concerned here with peak link bandwidth
ratings.  For example, even though the optimal throughput of strided
copies on the Cray T3D is 30--40 MBytes/sec~\cite{iscamemory}, current
implementations of sf2 achieve at best a measured and sustained
bandwidth of 10 MBytes/sec using the C interface to the
vendor-supplied MPI library~\cite{superquake96}.

Even more daunting is the goal of running the Quake applications at
50\% to 80\% efficiency on high-speed networks of workstations.  This
goal is attractive because it would make it possible to obtain enough
memory to run much larger simulations, but demands sustained per-PE
bandwidths of about 100 MBytes/sec, which is roughly an OC-24 rate.

\subsection{Bandwidth and latency tradeoffs}

There is a interesting tradeoff between the burst bandwidth and
latency that are necessary to achieve a target sustained PE bandwidth.
If block latency is made smaller, then burst bandwidth can be larger,
and vice-versa.  Equations~(\ref{eqn:macromodel}) and
(\ref{eqn:micromodel}) quantify these tradeoffs, which are shown in
Figure~\ref{fig:sf2bwlat} for the sf2 SMVP running on 128 200-MFLOP
PEs. Any point along a diagonal line represents a combination of burst
bandwidth ($T_w^{-1}$) and latency ($T_l$) that meets the sustained PE
bandwidth requirement ($T_c^{-1}$) from Figure~\ref{fig:sf2tc}(b).
The curves are generated by first using
Equation~(\ref{eqn:macromodel}) to compute the required inverse
sustained bandwidth $T_c$, and then using
Equation~\ref{eqn:micromodel} to determine the appropriate
combinations of $T_l$ and $T_w$.
\begin{figure*}[tbp]
\begin{center}
\begin{tabular}{c c}
\epsfxsize=3in
\epsfbox{sf2.128.lat.large.epsf}
&
\epsfxsize=3in
\epsfbox{sf2.128.lat.small.epsf}\\
&\\
(a) Arbitrarily large block size. & (b) Four-word block size.\\
\end{tabular}
\end{center}
\caption{burst bandwidth and latency tradeoffs for sf2/128 (200 MFLOP PE).}
\label{fig:sf2bwlat}
\end{figure*}

Figure~\ref{fig:sf2bwlat}(a) shows the burst bandwidth and latency
tradeoffs under the assumption that blocks are as large as possible,
and thus each PE sends at most one block to any other PE.  This is the
situation in message passing systems or DSMs that aggregate
blocks~\cite{treadmarkcompiler}.  The interesting point about this
graph is that latency matters for the SMVP. Even if burst bandwidth is
driven to infinity, observed block latency must not exceed $3 ~\mu$s if
the code is to run at 90\% efficiency.

We can also use model~(\ref{eqn:micromodel}) to get a sense of the
latency and bandwidth tradeoffs if blocks are small, fixed-sized
objects like cache lines.  To model the tradeoff for fixed-sized
blocks consisting of four 64-bit words, we set $B_{max} = C_{max}/4$
and then plot block latency as a function of burst bandwidth as
before. This is shown in Figure~\ref{fig:sf2bwlat}(b). Here, block
latency is even more crucial, as  expected. For example, if
burst bandwidth is infinite, then observed block latency must not
exceed 100 ns if the SVMP is to run at 90\% efficiency.

For a given application, a reasonable network design goal is to
achieve values for $T_l$ and $T_w^{-1}$ such that half of the time for
the communication phase is due to block latency and the other half to
burst bandwidth.  We refer to such a bandwidth as the {\em
half-bandwidth} and its corresponding latency as the {\em
half-bandwidth latency}.  These figures are useful targets for
architects of communication systems because overengineering a
network's bandwidth cannot relieve the latency constraint by more than
a factor of two, and vice-versa.  Figure~\ref{fig:sf2halfbw} plots the
half-bandwidths and half-bandwidth latencies for the entire space of
sf2 SMVPs, for both the maximal block size and the four-word block
size.
\begin{figure}[tbp]
\centerline{
\epsfxsize=4in
\epsfbox{sf2.halfbw.epsf}
}
\caption{Half-bandwidths and latencies for the sf2 SMVP.
         Note that the horizontal axis indicates decreasing bandwidth,
         for easy comparison with Figure~\protect\ref{fig:sf2bwlat}.}
\label{fig:sf2halfbw}
\end{figure}

The bottom portion of the graph shows the requirements when blocks are
fixed-sized four-word transfer units, as might occur in a
shared-memory architecture.  For the easiest case, running at 50\%
efficiency on 100-MFLOP PEs, we need a burst bandwidth of about 3
MBytes/sec with a block latency of about 10 $\mu$s.  In the most difficult
case, running at 90\% efficiency on 128 200-MFLOP PEs, we would need a
burst bandwidth of about 600 MBytes/sec and a block latency of about
70 ns.

The upper portion of the graph shows the requirements when blocks are
made as large as possible (i.e., each PE sends at most one block to
any other PE), as on a message-passing multiprocessor or network.
For the simplest case, running at 50\% efficiency on four
100-MFLOP PEs, we would need a burst bandwidth of about 3 MBytes/sec
and a half-bandwidth latency of about 8 ms.  However, the most
demanding case, where we are running at 90\% efficiency on 128
200-MFLOP PEs, requires a burst bandwidth of 600 MBytes/sec and a
block latency of about $2 ~\mu$s.

\section{Concluding remarks}
\label{sec:conclusions}

Irregular applications have been poorly understood in the past, in
large part because researchers and system builders have not had access
to large realistic example applications. This paper represents a step
toward understanding the implications of large-scale irregular codes
on the architecture of high-end systems.

We have chosen to study unstructured finite element simulations of
strong earthquake-induced ground motion in the San Fernando Valley of
Southern California because they constitute real and significant
engineering applications currently in daily use, and just as
important, their performance is characterized by the performance of a
simple SMVP kernel.  By examining the behavior of the SMVP, we gain
an understanding of the requirements for the communication systems of
future high-end computers and networks as the performance of commodity
PEs continues to increase. Of course, the results we have presented
are not sufficient; other types of applications must be characterized
as well. Still, sparse irregular codes like the Quake SMVPs are
crucial applications, and the results provide an important bit of
understanding about the communication performance of future systems.
Here are the main points:
\begin{enumerate}
\item The running time of the Quake unstructured finite element 
applications is dominated by SMVP operations, which also happen to be
the only operations in the code (aside from I/O) that require
communication.

\item Because of its irregular access pattern and large size, 
local floating-point performance for SMVP lags well behind peak
ratings. Achieved per-PE performance of only 10\% of peak is typical.

\item Blocks are small, even for big SMVP instances whose blocks
are as large as possible. Thus block latency is crucial and
cannot be amortized by large block transfers. Scaling up the size of 
an application does not help -- large messages 
will not save us for this class of computations!

\item Bisection bandwidth is not a critical issue for the SMVP.
Rather, overall performance is determined by the performance of the
communication system at each PE.

\item Systems with sustained (not peak) floating-point performance
of 200 MFLOPS per PE will need about 300 MBytes/sec of sustained communication
bandwidth per PE to run the SMVP codes efficiently. This sustained
bandwidth must include all overheads.

\item Even if blocks are maximally aggregated and burst bandwidths
are infinite, block latencies are a significant component of
communication performance. For a 200 MFLOP processor, if we strive for
a balanced system where half the communication time is due to latency,
then a block
latency of $2 ~\mu$s and a burst bandwidth of about 600 MBytes/sec are
required to run at 90\% efficiency.

\item Similar burst bandwidths are required if the transfer units are 
small fixed-sized blocks (like cache lines). However, the block
latencies must be more than an order of magnitude smaller, on the order
of 7 ns.  

\end{enumerate}

The bandwidth requirements per PE for future systems are aggressive,
especially since they must include all hardware and software
overheads. Yet they seem achievable. More worrisome are the block
latency requirements ($2 ~\mu$s for large blocks, 7 ns for small
blocks), which appear much more difficult to achieve.  The numbers in
this paper provide a strong quantitative argument for incorporating
latency reduction techniques in future systems, either with more
efficient interfaces, latency hiding techniques, or both.

\section*{Acknowledgements}
Many thanks to our Quake project partners from the CMU Dept of Civil
Engineering:  Jacobo Bielak, Omar Ghattas, Loukas Kallivokas, and
Hesheng Bao, who developed the San Fernando simulations and generously
shared their extensive knowledge. Joe Brandenburg at Intel sparked our
interest in characterizing large finite element simulations.
Performance numbers were obtained on the Cray T3E at Pittsburgh
Supercomputing Center.


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

\end{document}

