\section{Parallel Implementation of
\newline Matrix-Vector Product}

Here, we consider how a matrix-vector product $y = Ax$ might be implemented
for a vertex-partitioned graph and an element-partitioned graph. We assume
that $A$, $x$, and $y$ are distributed over the processors according to the
partition of $G$. Recall that $G$ represents the sparsity structure of $A$.
We will show that element-partitioning has distinct advantages.

After we have partitioned $G$ into $p$ blocks, let $G_c$ be the subgraph of
$G$ corresponding to the block that is mapped to processor $c$ by the 
placement algorithm (see Section \ref{placement}).
Let $G_a$ and $G_b$ be two such graphs, whose corresponding supervertices in
$G'$ are adjacent.

\subsection{Vertex-Partitioned Problems}

Suppose $G$ has been vertex-partitioned.  Then, for each vertex $u$ in $G_a$,
processor $a$ stores the value $x_u$ and the nonzero entries of row $u$ of
$A$. In order to calculate $y_u$, processor $a$ must have the entry of $x$ for
each vertex adjacent to $u$ in $G$.  In other words, processor $a$ must
allocate storage for every vertex adjacent to $G_a$ in $G$, as well as for
every vertex in $G_a$.  Thus, processors $a$ and $b$ share two ``layers'' of
vertices at the interface between $G_a$ and $G_b$.

The multiplication is performed in two phases.  The first phase is a
communication phase in which every processor receives needed entries of $x$
from its neighbors in graph $G'$, so that each processor has up-to-date
entries for every vertex it uses.  The second phase is a computation phase in
which each processor $p$ calculates $y_u$ for each vertex $u$ in $G_p$.

\subsection{Element-Partitioned Problems}

Now, suppose instead that $G$ has been element-partitioned.  Then, for each
vertex $u$ in $G_a$, processor $a$ stores the entry $x_u$. For each edge $(u,
v)$ in $G_a$, processor $a$ stores the value $A_{uv}$ if $a$ is the only
processor storing both $u$ and $v$. Otherwise, the value $A_{uv}$ is split
between several processors.

To understand how the matrix $A$ is stored, it is necessary to have some
understanding of how the sparse matrices used in finite element methods are
generated.  Typically, a separate calculation is made for each element.
Any element whose vertices include both $u$ and $v$ will make an additive
contribution to $A_{uv}$.  The process of summing the contributions of
individual elements to determine the value of $A$ is called {\it assembly}.

Figure \ref{elem} is an example of element-partitioning.
Suppose that the vertices $s$ and $t$ are each in both $G_a$
and $G_b$ (remember that in an element partition, vertices on the boundary
are contained in both subgraphs).  Suppose that edge $(s, t)$ is part of
two elements:  $e$ in $G_a$ and element $e'$ in $G_b$.  Then, the
contribution of $e$ to $A_{st}$ is stored on processor $a$, and the
contribution of $e'$ to $A_{st}$ is stored on processor $b$.  Let
$^pA_{st}$ denote the contribution to $A_{st}$ stored on processor $p$.  Then,
$A_{st} = \sum_{p}\,^pA_{st}$.  We say that $^pA$ is {\it partially assembled},
because it has been assembled everywhere except at processor boundaries.
Clearly, we could {\it fully assemble} $A$ by communicating the shared values
between processors and summing; but as we shall see, a fully assembled $A$
would be less useful.

\begin{figure}[htb]
\epsfxsize=3in
\centerline{\epsfbox{/afs/cs/usr/eschwabe/sc-paper/mist.eps}}
\caption{\label{elem} An element-partitioned mesh.}
\end{figure}

Again, the multiplication is performed in two phases; however, the
computation phase is first.  Each processor $p$ calculates a local
matrix-vector product of the form $^py =\,^pAx$.

The communication phase follows.  For each vertex $u$ shared by two
or more processors, these processors exchange their respective entries for
$^py_u$, and take the sum.  This gives us correct values because
$y = Ax = (\sum_{p}\,^pA)x = \sum_{p}\,^pAx = \sum_{p}\,^py$.
We see that vectors, like matrices, have both
partially assembled and fully assembled forms; and the product of a
partially assembled matrix and a fully assembled vector is a partially
assembled vector.  (There is no similar property for fully assembled
matrices.)

\subsection{Comparison of Partitioning
\newline Approaches}

When using the element-partitioning approach, processors $a$ and $b$
share only one ``layer'' of vertices at the
interface between $G_a$ and $G_b$, so less storage space is required for
a vector than when using the vertex-partitioning approach.  The sparse matrix
$A$, which dominates memory usage, requires less storage in the
element-partitioning approach if $A$ is symmetric, and requires less
storage in the vertex-partitioning approach if $A$ is asymmetric (for
reasons we omit here).  The element-partitioning approach makes initial
matrix generation much easier, because $A$ is built element-by-element,
so no communication is required.

So far, these are small differences.  There are two far greater advantages
to the element-partitioning approach.  First, in nonlinear problems where the
matrix $A$ changes over time, $A$ can be recomputed without communication.

Second, the element-partitioning approach
allows assembly of a vector to be delayed, and thus may reduce
the amount of communication which must be done.  For instance, suppose
one wishes to calculate $z = Ax + By$, where $A$ and $B$
are sparse matrices with the same sparsity structure.  (This arises in
the solution of second-order differential equations.)  The products
$Ax$ and $By$ each produce a partially assembled vector.
These two partially assembled vectors can be summed locally to form a
third partially assembled vector; this vector, if necessary, can be
assembled to form the result $z$.  Note that only one assembly, and
thus only one communication step, is needed for two matrix-vector
products.  This is called the {\em delayed assembly optimization}.

Because of the communication reduction and speed advantages of
matrix generation and the delayed assembly optimization, we strongly
prefer the element-partitioning approach.

\subsection{Code Generation Issues}

Our goal is a code generator for finite element methods,
that takes as its input the mesh $M$ and
a machine-independent description of a
finite element algorithm, and automatically partitions
the problem and generates multiprocessor code to execute the algorithm.

A code generator
can handle the delayed assembly optimization in two ways.  One is to
determine at compile time when vector
assembly should take place.  This requires
data dependence analysis, and cases will arise where it is not
applicable because the behavior of an algorithm is not always predictable.
The other solution is to defer all assembly decisions until runtime.  When
an operation requires a fully assembled vector as input, it checks the
vector provided and calls an assembly subroutine if necessary.  However,
this creates another problem:
inner products and elementwise multiplication require that at least one
of their operands is fully assembled.  If two partially assembled vectors
are provided, one must be chosen arbitrarily for assembly, with no
knowledge of which would be the optimal choice.  We have not yet surveyed
enough different finite element problems to resolve this issue.

We have not yet implemented a code generator.  We have, however, implemented
code for the iWarp multiprocessor to solve a second-order
time-varying equation used for earthquake simulation, with encouraging
results.

