\section{Algorithm Descriptions}
\label{sec:algorithms}
The primary algorithmic components of \BS\ are
\begin{itemize}
\item
the computation of a matrix ordering allowing for the 
scalable inversion of the triangular systems arising
from incomplete matrix factorizations;
\item
the automatic extraction of clique and identical node (i-node)
information from unstructured systems that allows for the use of 
higher-level BLAS; and
\item
the computation of incomplete sparse matrix factorizations for 
use as preconditioners in the iterative solution of linear systems.
\end{itemize}
The matrix-ordering algorithms are used to ensure that, for many 
problems, the performance per processor remains roughly constant
as the problem size and number of processors are increased.
Computational efficiency is obtained through the use of 
higher-level BLAS, efficient matrix storage schemes, and 
precomputed communication data structures.
For symmetric matrices, incomplete Cholesky factorizations may
be computed; and for general matrices, incomplete LU factorizations
may be obtained.
These preconditioners are designed so that their application 
is both efficient and as scalable as possible.

In addition, \BS\ contains implementations of several well-known
iterative methods.
A preconditioned conjugate gradient method can be used for
symmetric, positive definite systems.
For symmetric indefinite systems, the preconditioned SYMMLQ 
algorithm may be used.\footnote{The SYMMLQ 
algorithm requires a positive 
definite preconditioner, and this requirement can be a serious
limitation if the matrix being solved is very indefinite.
By ``very indefinite,'' we mean that the matrix has many negative
and many positive eigenvalues.}
For nonsymmetric systems, the GMRES method may be used.
In addition, the PETSc software package \cite{petsc} contains an interface to
the \BS\ preconditioners; this interface may be used 
if other iterative methods are required.
For basic information on iterative methods,
see \cite{template}.  
\sindex{conjugate gradients}
\sindex{GMRES}
\sindex{SYMMLQ}
\sindex{indefinite systems}

You can select from a number of different preconditioners.
The first is a simple diagonal scaling of the matrix, which
can be used by itself or in conjunction with another preconditioner.
The other preconditioning options are incomplete Cholesky or LU
factorization, SSOR ($\omega = 1)$, and block Jacobi 
(where the blocks are the cliques of the graph associated
with the sparse matrix).
Perhaps the most generally applicable selection is the incomplete 
factorization with diagonal scaling.\footnote{Two 
possible exceptions to this recommendation are (1) if the matrix
has no or very small cliques and identical nodes
(in which case the factorization may be very slow) and
(2) if the space for the incomplete factorization is not available.}
The incomplete factorization is the algorithm for which \BS\ was 
designed; this approach has proved useful in a wide variety of 
practical problems.
\sindex{diagonal scaling}
\sindex{incomplete Cholesky}
\sindex{incomplete LU}
\sindex{identical nodes (i-nodes)}
\sindex{cliques}
\sindex{SSOR}
\sindex{block Jacobi}

\BS\ does not partition a matrix across the processors.
Instead, \BS\ assumes that the given partitioning is a good one.
As such, its performance may be limited by the quality of the
partitioning.
We assume that the right-hand side and the solution vector are 
\sindex{matrix partitioning}
partitioned in the same manner as the rows of the sparse matrix.
See \cite{chaco2} for more information on partitioning heuristics.

We achieve parallelism in the conjugate gradient, SYMMLQ, and
GMRES implementations by partitioning the vectors used in these
algorithms in the same manner that the rows of the matrix are 
partitioned across the processors.
For these algorithms it is (for the most part) a simple 
\sindex{parallel inner products}
matter of executing inner products and {\tt daxpy}
operations in parallel. 

\subsection{Processor Performance}
To achieve good performance on each node, we reorder
the matrix and use a layered data structure
to allow the use of the higher-level dense BLAS.  
This reordering is particularly important on machines that 
use high-performance RISC chips on which good performance 
can be achieved only by using such operations.
The reordering of the matrices is based on the identification 
of two structures that commonly arise in the graph associated
with the matrix nonzeros:  identical nodes and cliques.

Identical nodes (i-nodes) typically exist when multiple degrees
of freedom are associated with each vertex in the graph.
Cliques are found in many graphs associated with sparse matrices,
but cliques typically are found in graphs where multiple degrees
of freedom are associated with each vertex and the local connectivity
of the graph is large.
For example, if one uses a second-order, three-dimensional
finite element in a typical structural engineering problem
(with three degrees of freedom per vertex), clique
sizes of up to 81 can be found.  In general, the larger the cliques
or identical nodes, the better the performance.
This technique has been used with great success in direct matrix
factorization methods.
\sindex{BLAS}
\sindex{processor performance}
\sindex{cliques}
\sindex{identical nodes (i-nodes)}

We illustrate these structures in Figure~\ref{fig:inode} where
we depict a subsection of a graph
that would arise from a two-dimensional, bilinear, multicomponent
finite-element model with three degrees of freedom per 
discretization point.
We illustrate the three degrees of freedom by the three
dots at each node point; the linear, quadrilateral elements imply 
that the twelve degrees of freedom sharing the four node points
of each face are completely connected.
In the figure we show edges only between the nodes, these
edges represent the complete interconnection of all the
vertices on each element or face.

\begin{figure} [htb]
   \centerline{
		\hbox{
		\psfig{figure=grid1.ps,height=2.0in} \hspace{0.25in}
		\psfig{figure=grid2.ps,height=2.0in}
	  }
 }
   \centerline{
		\hbox{
		\psfig{figure=grid3.ps,height=2.0in}
	  }
 }
\caption{\label{fig:inode}A subgraph generated by
a two-dimensional, bilinear finite-element model with three
degrees of freedom per discretization point.
The geometric partition shown by the dotted lines yields 
an assignment of the vertices in the enclosed subregion
to one processor.
The upper right figure shows one possible decomposition
of the graph into local cliques, and the lower figure
shows the corresponding clique (or supernode) graph.
}
\end{figure}

The dashed lines in the figure represent a geometric partitioning
of the grid; we assume that the vertices in the central region
are all assigned to one processor.
%We make several observations about the local structure of this subgraph. 
We note that the adjacency structures of the vertices
at the same geometric node
(i.e., the nonzero structure of the associated variables)
are identical.\sindex{identical nodes (i-nodes)}
\BS\ takes advantage of these so-called i-nodes by maintaining only one
copy of the indexing for each set of identical rows and by storing 
nonzeros for these rows in a dense matrix form that can be used
with the BLAS.

A second data structure layer is determined by the {\em clique}
structure of the graph.
Recall that a clique is a completely connected subgraph.
In the upper right of Figure~\ref{fig:inode}, 
the grouping of vertices by the dotted lines partitions
the graph into cliques.
\BS\ associates a ``supernode'' with each clique and
orders the unknowns associated with each clique consecutively.
Based on this ordering, the submatrix corresponding to the
unknowns in the same clique corresponds to a dense matrix.
In the sparse matrix, these dense matrices manifest themselves
as dense blocks on the matrix diagonal.
The \BS\ data structures
explicitly store these dense blocks to allow for greater
efficiency in using the BLAS.

\subsection{Multicolor Ordering}
Following the graph reductions that are accomplished after
the identification of i-nodes and cliques in the sparse matrix
structure, \BS\ must construct a matrix ordering based on the 
reduced graph, illustrated at the bottom of Figure~\ref{fig:inode}.
To achieve scalable parallel performance in the incomplete
factorization and SSOR preconditioners, \BS\ colors the reduced 
graph using a parallel coloring heuristic \cite{jp:pcolor}.
The graph reduction and coloring phases of the computation
are efficient and typically require a small amount of 
time relative to the incomplete matrix factorization routine.

The advantage of using a graph coloring is that the number of
colors required is essentially a {\em local} property of the
graph; thus, the number of colors is roughly independent
of the number processors used for a fixed discretization scheme.
This fact allows for the scalable performance of the \BS\
package \cite{jp:scalable}.
The trade-off for this scalability is that these multicolor
orderings may not be the optimal orderings to choose 
for minimizing the number of iterations required for the 
iterative solver; the convergence of the iterative solver
with respect to different orderings is highly problem dependent.
However, the combination of coloring a general symmetric sparse 
matrix and the incomplete Cholesky algorithm has proved to be 
successful for solving large problems on scalable parallel 
architectures \cite{jp:bell_prize}.
In \cite{jp:copper} we address the issue of convergence of 
this combination of algorithms for the model problem.
\sindex{scalable performance}
\sindex{graph coloring}

\subsection{Communication Efficiency}
The computational kernel of the \BS\ package is the solution of 
triangular matrix systems and triangular matrix multiplication.
Triangular matrices of the same structure are solved and multiplied
in parallel many times in the use of one of the iterative solvers.
The key observation is that the required interprocessor communication 
in triangular system solve or multiply is the same every time; the 
only difference is in the floating-point values that are sent and
received.

To take advantage of this repeated communication pattern, the \BS\
routines in {\tt BMcomp\_msg()}\findex{BMcomp_msg()}
are used to build reusable data structures.
During each communication phase each processor needs to send and receive 
data from some subset of the processors (a small subset, if the 
partitioning is good). 
Each message to be sent is an array of floating-point values and an 
array of integer values (matrix indices).
A key observation is that these integer values do not change; only the 
floating-point values change.

Prior to any triangular matrix solution and multiplication, each processor
declares (via function calls in {\tt BMcomp\_msg()})\findex{BMcomp_msg()}
the messages that it will be sending and receiving and specifies the integer 
index values associated with these messages.
The routines in {\tt BMcomp\_msg()} build data structures by assembling the
declarations from each processor.
In addition, the integer index values associated with the messages
are sent to the appropriate processor and stored there; these values will
not be sent in the future.  
Since the size of the required messages is now known, message buffers
can be preallocated.

We currently perform no special communication pattern optimizations in 
{\tt BMcomp\_msg()} although these routines have been set up 
for this possibility.
The optimizations currently performed are
\begin{enumerate}
\item
presending of integer values;
\item
preallocation of message buffers, allowing ``forced'' messages;
{\tt MPI\_Irsend()}, to be used; and
\item
no requests are made between processors for non-local data
(processors use local information to know what messages to expect 
and what to send).

\end{enumerate}
For an example of the message setup process, see the routine
{\tt BSsetup\_forward()}.\findex{BSsetup_forward()}
For an example of how the message data structures are used, see the
routine {\tt BSforward()}.\findex{BSforward()}

\newpage
