\documentstyle[11pt,titlepage]{article}
%\nonstopmode
\setcounter{secnumdepth}{5}
\setcounter{tocdepth}{5}
\renewcommand{\baselinestretch}{1.30}
\textwidth=6.0in
\textheight=9.0in
\topmargin=-0.5in
\oddsidemargin=0.25in

\input{psfig}

\begin{document}

\section{Executive Summary \label{ex}}

To date, we have performed a theoretical examination of the available
parallelism in block-diagonal-bordered matrices \cite{DPK-par-dir} and
developed diakoptic-based node-tearing software to order real power
system matrices into block-diagonal form
\cite{DPK-scale,DPK-par-dir}.  We have developed a parallel
block-diagonal-bordered Choleski solver and implemented it on the
Thinking Machines CM-5 using a message-passing paradigm.  The results
of our research to date show that the performance of our parallel
Choleski algorithm is limited by the limited amount of calculations
when solving load-flow matrices \cite{DPK-par-dir}.  An overview of
this work is included in this report.  The sparse matrices obtained
when performing transient stability analysis have additional available
calculations that result from solving the generator equations.  These
equations are naturally in block-diagonal-bordered form and when the
calculations in these mutually independent diagonal blocks are
combined with the independent calculations available after ordering
the network equations, we believe that parallel LU factorization
algorithms can be developed that yield reasonably efficient
performance for eight or more processors.

Tasks remaining in our analysis of block-diagonal-bordered sparse
matrix solvers include porting the existing sparse Choleski solver to
a portable communications library, so the algorithm can be tested on
various distributed-memory multiprocessors and distributed computing
environments.  We also plan to develop and implement two additional
block-diagonal-bordered sparse matrix solvers, a parallel direct
solver based on parallel block-diagonal-bordered LU factorization and
an iterative parallel block-diagonal-bordered Gauss-Seidel solver.
The LU factorization-based solver will be designed for transient
stability analysis, where portions of the matrix are diagonal blocks
of generator equations and the remaining partitions of the matrix are
the generator equations.  Our past work on ordering the network
equations will permit additional mutually independent diagonal blocks
to improve performance.  The performance of the two direct techniques
will be compared to the performance of the block-diagonal-bordered
Gauss-Seidel-based solver.


\section{Parallel Sparse Matrix Algorithms \label{sparse}}

We have examined the applicability of parallel direct
block-diagonal-bordered sparse solvers for real power system load-flow
applications that require the solution of symmetric positive definite
sparse matrices.  Parallel block-diagonal-bordered sparse linear
solvers offer the potential for regularity often absent from other
parallel sparse solvers.  Load flow analysis entails the solution of
non-linear systems of simultaneous equations, which are performed by
repeatedly solving sparse linear equations.  For power system
applications, however, the limited size of the matrices and load
imbalance due to limited parallelism in the matrix structure
significantly limits the number of processors that can be used
efficiently for a single parallel Choleski solver.  This fact will be
evident in the empirical data collected on the CM-5.  Our research
into specialized matrix ordering techniques has shown that it is
possible to order actual power system matrices readily into
block-diagonal-bordered form, but load imbalance becomes excessive
beyond four processors, limiting potential parallelism for a single
parallel Choleski solver within an application.  Nevertheless, other
dimensions exist in electrical power system applications that can be
exploited to efficiently make use of large-scale multi-processors.  

In this portion of the report, we review Choleski factorization and
the ordering of positive definite matrices.  We then discuss the
available parallelism in block-diagonal-bordered form matrices and we
discuss the three step preprocessing phase required to place spare
power systems matrices in block-diagonal-bordered form.  We then 
discuss the performance of our parallel block-diagonal-bordered
Choleski solver and the performance of clustering software required to
order matrices into block-diagonal-bordered form.


\subsection{Choleski Factorization \label{Choleski}}

We are considering the direct solution of the linear system
\begin{equation}
Ax=b,\label{Axeqb}
\end{equation}
where $A$ is an $N \times N$ sparse matrix. The sparse
matrix $A$ can be numerically factored into two separate triangular
matrices, one sparse matrix being lower triangular, $L$, and the other
sparse matrix being upper triangular, $U$:
\begin{equation}
Ax=LUx=b,
\end{equation}
A lower triangular matrix, $L$, has all zeros above the diagonal and
an upper triangular matrix, $U$, has all zeros below the diagonal \cite{Duff}.

If the matrix $A$ is an $N \times N$ symmetric positive definite
sparse matrix, then a special form of LU factorization can be used
that exploits the symmetry and inherently numerical stable
characteristics of this matrix form \cite{Duff}.  Symmetric positive definite
matrices occur in power system load flow analysis.  Our analysis to
examine the available parallelism in block-diagonal-bordered LU
factorization can be readily extended to block-diagonal-bordered
Choleski factorization with minimal modifications \cite{DPK-par-dir}.

For a brief review, the symmetric positive definite sparse matrix $A$
can be numerically factored into a single lower triangular matrix $L$:
\begin{equation}
Ax=L L^{T}x=b, \label{lin-sys}
\end{equation}
A lower triangular matrix, $L$, has all zeros below the diagonal.
Equation~\ref{lin-sys} is solved by setting $L^{T}x=y$, and
substituting $y$ for $L^{T}x$.  The numerical solution for $Ly=b$ is
found by forward reduction, and the numerical solution for $x$ is
calculated by backward substitution in the equation $L^{T}x=y$.
Triangular linear systems can be readily solved numerically by solving
for the first value in the triangular linear system and substituting
that value into subsequent equations \cite{Duff}.

Symmetric sparse matrices can be represented by graphs with elements
in equations corresponding to undirected edges in the graph
\cite{Duff,Heath}.  Ordering a symmetric sparse matrix is actually little more
than changing the labels associated with nodes in an undirected graph,
however, this simple task can drastically effect the amount of
calculations involved when factoring a sparse matrix by affecting the
order of calculation and consequently the amount of fillin.  For
symmetric positive definite matrices, there is much latitude in the
order to perform the calculations, because there is no requirement for
pivoting for numerical stability and the only effect of modifying the
order of calculations might result from changes in round-off errors
\cite{Heath}.  There are several notable techniques to minimize fillin, 
with one of the commonly used techniques being minimum-degree
ordering.  We use minimum degree ordering in conjunction with
node-tearing-based ordering technique to generate
block-diagonal-bordered form sparse matrices \cite{DPK-par-dir}.

Modifying the ordering of a sparse matrix is simple to perform using a
permutation matrix $P$ of all zeros and ones that simply generates
elementary row and column exchanges.  Applying the permutation matrix
$P$ to the original linear system in equation~\ref{Axeqb} yields the
linear system
\begin{equation}
(PAP^{T})(Px) = (Pb),\label{PAPxeqb}
\end{equation}
that is solved by factoring $PAP^{T}$ into the Choleski factor
$\bar{L}$ in $\bar{L}{\bar{L}}^{T}$ and then performing forward
reduction, backward substitution, and undoing the permutation on the
$x$ vector.  These steps would require the solutions of:
\begin{equation}
\bar{L} = Pb,
{\bar{L}}^{T}z = y,
x=P^{T}z.\label{fr-bs-undo}
\end{equation}
As long as the symmetric matrix $A$ is ordered with the permutation
matrix $P$ to $PAP^{T}$, the resultant matrix after ordering remains
symmetric positive definite \cite{Duff}.

\subsection{Available Parallelism in Block-Diagonal-Bordered Form Matrices \label{Shur}}

The most significant aspect of parallel sparse Choleski factorization
is that the sparsity structure can be exploited to offer more
parallelism than is available with dense matrix solvers.  Parallelism
in dense matrix factorization is achieved by distributing the data in
a manner that the calculations in one of the {\bf for} loops can be
performed in parallel.  Sparse factorization algorithms have
inadequate calculations in any row or column for efficient
parallelism; however, sparse matrices offer additional parallelism as
a result of the nature of the data and the precedence rules governing
the order of calculations.  Instead of just parallelizing a single
{\bf for} loop as in parallel dense matrix factorization, entire
independent portions of a sparse matrix can be factored in parallel
--- especially when the sparse matrix has been ordered into
block-diagonal-bordered form.  Our analysis of available parallelism
is based on partitioning the matrix and calculating the Shur
complement \cite{DPK-par-dir}.  This work has some similarities with
elimination graphs and super-nodes described in \cite{Heath}; however,
this work adds a two-dimensional flavor to these concepts.  Provided
that the matrix can be clustered into block-diagonal-bordered form,
then the parallel sparse Choleski algorithm can reap additional
benefits, such as the elimination of task graphs for
distributed-memory multi-processor implementations.  Minimizing or
eliminating task graphs is significant because the task graph can
contain as much information as the representation of the sparse matrix
for more conventional parallel sparse Choleski solvers
\cite{shar-mem}.

Figure~\ref{BBDM} illustrates both the Choleski factorization steps
and the reduction/substitution steps for a block-diagonal-bordered
sparse matrix.  In this figure, the strictly lower diagonal portion of
the matrix is ${\bf L}$, while the strictly upper diagonal portion of
the matrix is ${\bf L}^{T}$.  This figure depicts four diagonal
blocks, and processor assignments (P1, P2, P3, and P4) are listed with
the data block.

\begin{figure}
\begin{center}
\begin{picture}(450,234)
\psfig{figure=fig-BBDM.ps,width=432pt,height=234pt}
\end{picture}
\caption{Block Bordered Diagonal Form Sparse Matrix Solution Steps\label{BBDM}}
\end{center}
\end{figure}

\subsection{The Three-Step Preprocessing Phase\label{3step}}

For parallel sparse block-diagonal-bordered
matrix algorithms to be efficient when factoring irregular sparse
matrices, the following three step preprocessing phase must be
performed:
\begin{itemize}
\item {\em order} the matrix into block-diagonal-bordered form while minimizing the number of calculations,
\item {\em pseudo-factor} to identify both fillin and the number of calculations for all diagonal blocks and corresponding portions of the borders, and
\item {\em load balance} to uniformly distribute the calculations among processors.
\end{itemize}
The first step determines the block-diagonal-bordered form and the
ordering of nodes within diagonal blocks to minimize calculations; the
second step determines the locations of fillin values for static data
structures and also determines the number of calculations in
independent blocks for the load balancing step; and the third step
determines a mapping of data to processors for efficient
implementation of the algorithm for the user specified data.  These
three steps may be incorporated into an optimization framework that
uses the three-step preprocessing phase to produce matrix orderings
with optimal overall performance for a particular version of the
block-diagonal-bordered sparse matrix factorization algorithm.  For
this report, the optimization was performed by hand --- various values
of input parameters for the clustering routine were examined and the
block-diagonal-bordered form sparse matrix with the best load balance
and least numbers of operations were chosen for collecting performance
benchmarks.

The metric for load balancing or the metric for an optimization
routine to determine the most efficient overall ordering technique
must be based on the actual workload required by the individual
processors.  This number may differ substantially from the number of
equations assigned to processors because the number of calculations in
an independent sub-block is a function of the number and location of
non-zero matrix elements in that block --- not the number of equations
in a block.  For dense matrices, the computational complexity of
factorization is $O(N^{3})$, however, the computational complexity for
factoring individual diagonal matrix blocks in the real power systems
matrices used in later parallel algorithm performance studies can vary
from $O(N^{1.1})$ to $O(N^{2.0})$.  The computational complexity for
factoring those entire sparse matrices can vary from $O(N^{1.30})$ to
$O(N^{1.55})$.  Determining the actual workload requires a detailed
simulation of all processing for the factorization and triangular
solution phases, which we refer to as pseudo-factorization.

The ordering portion in the preprocessing phase must identify diagonal
matrix blocks while also attempting to minimize the amount of fillin
during factorization.  We chose to use node tearing
\cite{Duff,NodeTear}, a specialized form of diakoptics \cite{Diak} to 
order matrices into block-diagonal-bordered form.  This technique
attempts to extract the natural structure in the matrix or graph, and
generally produces many irregularly sized blocks, while minimizing the
number of coupling equations or the size of the lower border and last
diagonal block.  Load balancing techniques must be used after the node
tearing matrix ordering step to uniformly distribute the processing
load onto a multi-processor.  

To accomplish the two-fold requirement for both identifying the
location of fillin and determining the amount of calculations in each
independent block, we utilize a pseudo-factorization step as part of
the preprocessing phase.  Pseudo-factorization is merely a replication
of the numerical factorization process without actually performing the
calculations.  Counters are used to tally the numbers of calculations
to factor the independent data blocks and the numbers of calculations
to update the last block using data from the borders.  There is no way
to avoid the computational expense of this preprocessing step, because
the computational workload in factorization is not correlated with the
number of equations in an independent block.  Efficient parallel
sparse matrix solvers require that any disparities in processor
workloads be minimized in order to minimize load imbalance overhead,
and consequently, to maximize processor utilization.

The load balancing step of the preprocessing phase can be performed
with a simple pigeon-hole type algorithm that uses a metric based on
the numbers of calculations determined in the pseudo-factorization
step.  There are three distinct steps in the proposed
block-diagonal-bordered matrix solver:
\begin{itemize}
\item factor independent blocks,
\item update the last block using data from the borders,
\item factor the last block.
\end{itemize}
Load balancing, as implemented here, emphasizes the uniform
distribution of the processing workload in the first two steps.  Our
research into this area has emphasized uniformly distributing the
workload to separate processors based on the number of calculations
when factoring both the independent blocks and calculating the updates
of the last block from data in the borders \cite{DPK-scale}.The
parallel calculations in the last diagonal block are load balanced
separately, which is simple, because we are using a parallel blocked
kij-saxpy Choleski algorithm to factor the last diagonal block.

\subsection{Empirical Results --- Sparse Choleski Solver with Test Data\label{results-dg}}

The first step in understanding the performance potential of this
block-diagonal-bordered Choleski solver is to examine the parallel
algorithm performance with machine generated test data that has
perfect load balance for all processors.  This data has been machine
generated with equal numbers of calculations in diagonal-blocks to
eliminate requirements for load balancing.  The pattern of the
matrices is a recursive block-diagonal-bordered form as illustrated in
figure~\ref{TEST}.  This figure depicts the sparsity structure in the
machine generated test matrix, where non-zero values are black and the
remainder of the matrix is comprised of zero-values.  A bounding box
has also been placed around the example matrix.  This matrix has been
used in tests to examine speedup and sources of parallel processing
overhead.  The block-diagonal-bordered Choleski solver has been run on
this machine generated matrix using one to sixteen processors on the
Thinking Machines CM-5.  This matrix has sixteen separate major
diagonal blocks, each with 128 diagonal elements.  The last diagonal
block also has 128 diagonal elements, with eight columns per diagonal
block in the lower border.  The size of this matrix was chosen to be
similar to the BCSPWR09 test matrix, although it has a substantially
higher computational complexity than the BCSPWR09 matrix,
$O(N^{1.85})$ versus $O(N^{1.38})$

\begin{figure}
\begin{center}
\begin{picture}(216,216)
\psfig{figure=MG-test.eps,width=216pt,height=216pt}
%\psfig{figure=TBD.eps,width=216pt,height=216pt}
\end{picture}
\caption{Machine Generated Test Matrix\label{TEST}}
\end{center}
\end{figure}

Performance data has been collected for distinct operations in the
factorization and triangular solution of block-diagonal-bordered
matrices to examine the speedup of each distinct portion of the
algorithm.  Data was collected in such a manner as not impact the
overall measures of performance.  Performance of the multiprocessor
algorithms are illustrated using graphs plotting relative speedup
versus the number of processors (one to sixteen).

{\bf Definition --- (Relative Speedup)} {\it Given a
single problem with a sequential algorithm running on one processor
and a concurrent algorithm running on p independent processors,
relative speedup is defined as }
\begin{equation}
S_{p} \equiv \frac{T_{1}}{T_{p}} ,
\end{equation}
{\it where $T_{1}$ is the time to run the sequential algorithm as a
single process and $T_{p}$ is the time to run the concurrent algorithm
on p processors \cite{Fox}.}

A graph of speedup calculated from empirical performance data for the
test matrix is provided in figure~\ref{TEST-perf-1}.  This figure has a
family of three curves that show speedup for:
\begin{enumerate}
\item Choleski factorization 
\item forward reduction and backward substitution (triangular solve)
\item a combination of factorization and a single forward reduction and backward substitution
\end{enumerate}
In general this speedup graph shows that the block-diagonal-bordered
Choleski solver has significant potential for efficient operations,
because speedups of over eleven were obtained with sixteen processors.
This equates to a processor utilization efficiency of 70\%.  When
examining the runtime timing data for each processor, near perfect
speedup (speedup equal to the number of processors) is not achieved
because of overhead that occurs during those portion of the algorithm
that use asynchronous pipelined communications.  The following sections
of the block-diagonal-bordered Choleski solver use asynchronous
pipelined communications:
\begin{itemize}
\item the update of data in the last block using data from the lower border
\item the dense blocked kij-saxpy based Choleski factorization of the last block
\item the update of data in the b-vector partition using data from the lower border
\item the forward reduction of the last block
\item the backward substitution of the last block
\end{itemize}

\begin{figure}
\begin{center}
\begin{picture}(216,144)
\psfig{figure=gplotout7.ps,width=216pt,height=144pt}
\end{picture}
\caption{Speedup for Generated Test Data --- 2, 4, 8, and 16 Processors \label{TEST-perf-1}}
\end{center}
\end{figure}

As a result of the regularity in this test matrix, it is possible to
examine the performance of each portion of the parallel
block-diagonal-bordered Choleski factorization algorithm.  A graph
illustrating speedup calculated from the empirical performance data is
presented in figure~\ref{TEST-perf-2}.  In this graph, the speedup for
factorization of the diagonal blocks show perfect speedup, which is
not unexpected because there is no communications in this step.
However, the measured speedup for the two sections of the Choleski
factorization algorithm that require communications, is not as
impressive, due to communications overhead.  As we use more
processors, there is more communications and less work for each
processor.  It takes longer to fill communications pipelines and there
are less calculations per processor to help mitigate the cost of
communications.  Consequently, there are performance limits to this
algorithm.  These limits are not a function of a sequential component
of the parallel calculations (Amdahl's Law) \cite{Fox}, but rather the
limits are a function of what percentage of the calculations are
perfectly parallel and what percentage of the calculations have
asymptotic performance limits due to communications overhead.

\begin{figure}
\begin{center}
\begin{picture}(216,144)
\psfig{figure=gplotout8.ps,width=216pt,height=144pt}
\end{picture}
\caption{Speedup for Choleski Factorization Algorithm Components  --- 2, 4, 8, and 16 Processors \label{TEST-perf-2}}
\end{center}
\end{figure}

Parallel block-diagonal-bordered Choleski factorization algorithm
performance can be viewed as a linear combination of the individual
performance of these three portions of the algorithm.
\begin{equation}
S_{p} = (\lambda_{DB} \times S_{DB,p} + (\lambda_{U} \times S_{U,p}) + (\lambda_{LB} \times S_{LB,p}) \label{linear-comb}
\end{equation}
\begin{tabbing}
where: \= \\
       \> $S_{p}$ is the $p$ processor speedup for block-diagonal-bordered Choleski factorization. \\
       \> $S_{DB,p}$ is the $p$ processor speedup for factoring of the diagonal blocks. \\
       \> $S_{U,p}$ is the $p$ processor speedup for updating the data in the last block using data in the border. \\
       \> $S_{LB,p}$ is the $p$ processor speedup for factoring the last block. \\
       \> $\lambda_{DB}$ is the fraction of time spent factoring of the diagonal blocks. \\
       \> $\lambda_{U}$ is the fraction of time spent updating the data in the last block using data in the border. \\
       \> $\lambda_{LB}$ is the fraction of time spent factoring the last block. \\
       \> $\lambda_{DB}+\lambda_{U}+\lambda_{LB} = 1$
\end{tabbing}
Our empirical data shows that $S_{p_{DB}} \approx p$ because this step
has no communications, so the performance of this factorization
algorithm on a problem using other block-diagonal-bordered sparse
matrices can be estimated as a function of the percentage of
calculations in the diagonal blocks.  Similar analyses can be
performed for forward reduction and backward substitution.  Parallel
performance of those sections of the triangular solution algorithms
that require pipelined communications will have less speedup potential
than the triangular solution of data in the diagonal blocks.
Unfortunately, the calculations in the diagonal blocks for the
triangular solution will be of a lesser computational complexity than
the calculations in the factorization step.  Consequently, poorer
parallel performance, lower speedup, would be expected in forward
reduction and backward substitution than can be obtained in
block-diagonal-bordered factorization.

There are two basic forms of communications used in this algorithm.
\begin{enumerate}
\item Pipelined broadcasts for those portions of the algorithm that factor
or solve the last diagonal block
\item Pipelined, asynchronous, everyone-to-everyone communications in both 
factorization and forward reduction when
sums of matrix $\times$ matrix products or vector $\times$ matrix
products from data in the lower border are used to update the last
diagonal block or the b-vector partition associated with the last diagonal block.
\end{enumerate}
In both instances, the pipelined communications appears to {\em
data-starve} and not hide communications behind the calculations.
This is not surprising in the forward reduction of the last block, but
this phenomenon also occurs in the block factorization step.
Experiments with adjusting the block size offered little improvement.
While the communications have been made regular in the other
communications steps, it appears that communications times are so
large relative to calculations in the same processing step, that
little benefit can be made of hiding communications behind
calculations.  Further research into other asynchronous communications
techniques such as active messages \cite{Active-1} may be able to
significantly improve performance in this area.

This analysis has been performed on data that had no load imbalance
overhead.  Additional sources of overhead would degrade potential
performance of the algorithm \cite{Fox}.  We have discussed the
effects of communications overhead, nevertheless, other sources of
parallel processing overhead are possible.  They include indexing
overhead and load imbalance overhead.  Indexing overhead is the extra
computational cycles required to set up loops with multiple
processors.  Care has been exercised in the development of this
algorithm to minimize indexing overhead.  While it is not possible to
eliminate the cost of repeating the process of setting up loops on
each processor, the amount of calculations required to keep track of
the values as a function of individual multi-processors can be
minimized.  Many of the sources of indexing overhead have been
accounted for in the preprocessing stage that prepares a matrix form
for processing by this algorithm.  The computational cost of the
preprocessing stage is expected to be amortized over multiple uses of
the block-diagonal-bordered Choleski, thus the costs incurred in this
stage are not addressed here.  Load imbalance overhead is also
addressed in the preprocessing stage, although the sources of load
imbalance overhead are a function of the interconnections in the real
data sets.  In the preprocessing phase, we use load-balancing to
attempt to order the data as to minimize this source of overhead in
our parallel algorithm.  We will show in section~\ref{results-solve}
that even with attempts at load balancing, this source of overhead
cannot be entirely removed.

\subsection{Empirical Results --- Clustering \label{results-cluster}}

Performance of our block-diagonal-bordered Choleski solver will be
analyzed with two separate power systems matrices:
\begin{itemize}
\item Boeing-Harwell matrix BCSPWR09 --- 1723 nodes and 4117 edges in the graph \cite{BH-guide}
\item Boeing-Harwell matrix BCSPWR10 --- 5300 nodes and 13571 edges in the graph \cite{BH-guide}
\end{itemize}
Matrices BCSPWR09 and BCSPWR10 are from the Boeing Harwell series and
represent real electrical power system networks from the Western and
Eastern US respectively.  

A detailed clustering analysis has been performed on the BCSPWR09 data
to illustrate the ability of the node-tearing clustering algorithm
\cite{DPK-par-dir}.  Our block-diagonal-bordered Choleski solver
requires that the BCSPWR09 matrix be ordered into
block-diagonal-bordered form with uniformly distributed workload at
each of the processors.  A single specified input parameter, the
maximum cluster size, defines how the node-tearing algorithm clusters
the matrix.  An example of applying the node-tearing algorithm to the
BCSPWR09 matrix is presented in figure~\ref{NT-3} for maximum diagonal
block size of 128 nodes.  Original nonzero values are represented in
this figure in black, fillin locations are represented in gray, and
all remaining zero values are white.  A bounding box has been placed
around the sparse matrix.  The number of fillin for this
block-diagonal-bordered form matrices is 3,838, with a substantial
portion of the fillin occurring in the lower right-hand corner, or the
last diagonal-block.  The Boeing-Harwell data is comparable to data
obtained from Niagara Mohawk operations.  That matrix has 1766 nodes
and 4272 edges in the graph.  Examples of applying the node-tearing
algorithm to Niagara Mohawk operations data is presented in
figure~\ref{nimo}.

\begin{figure}
\begin{center}
\begin{picture}(216,216)
\psfig{figure=bcs09-bdb128.eps,width=216pt,height=216pt}
%\psfig{figure=TBD.eps,width=216pt,height=216pt}
\end{picture}
\caption{Node-Tearing-Based Ordering --- BCSPWR09 --- $max_{DB}= 128$ \label{NT-3}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(216,216)
\psfig{figure=nimo-OPS-bdb128.eps,width=216pt,height=216pt}
%\psfig{figure=TBD.eps,width=216pt,height=216pt}
\end{picture}
\caption{Node-Tearing-Based Ordering --- Niagara Mohawk Operations Data --- $max_{DB}= 128$ \label{nimo}}
\end{center}
\end{figure}

The general performance of node-tearing-based ordering is dependent on
the maximum number of nodes in a diagonal block, and the intended
number of processors to which the mutually independent diagonal blocks
will be distributed.  The number of coupling equations and the size of
the last block is dependent only on the maximum number of nodes in a
diagonal block, which is illustrated in figure~\ref{NT-4}.  Note that
the maximum size of the diagonal blocks is inversely related to the
size of the last diagonal block.  This is intuitive, because as
diagonal matrix blocks are permitted to grow larger, multiple smaller
blocks can be incorporated into a single block.  Not only will the two
blocks be consolidated into the single block, but in addition, any
elements in the coupling equations that are unique to those network
clusters could also be moved into the larger block.  Another
interesting point with the relationship between maximum size of the
diagonal block and the size of the last block, is that the percentage
of non-zeros and fillin in the last diagonal block increases
significantly as the size of the last block decreases.  This makes the
use of parallel dense Choleski factorization techniques even more
desirable for the last diagonal block when the selection of the
maximum size of the diagonal blocks is chosen to minimize the size of
the last diagonal block.

It is desirable to minimize the size of the last block, but this can
cause load imbalance as the number of processors increases.  Load
imbalance is defined as the ratio of the difference of the maximum and
minimum numbers of floating point operations divided by the maximum
number of floating point operations per processor and illustrates the
percentage of time that the processor with the least amount of work is
inactive.  Families of curves illustrating load imbalance as a
function of the number of processors is presented in
figure~\ref{NT-5}.  The load imbalance data presented in this figure
are calculated only for the number of operations required to calculate
the Choleski factor of the matrix.  For all ordering with different
size diagonal blocks, there is perfect load-balancing for four
processors, however, as the number of processors increases, so does
load imbalance.  

\begin{figure}
\begin{center}
\begin{picture}(216,144)
\psfig{figure=gplotout2.ps,width=216pt,height=144pt}
\end{picture}
\caption{Node-Tearing --- Number of Nodes in Last Block for BCSPWR09 \label{NT-4}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(216,144)
\psfig{figure=gplotout1.ps,width=216pt,height=144pt}
\end{picture}
\caption{Node-Tearing --- Load Imbalance for Factoring BCSPWR09 \label{NT-5}}
\end{center}
\end{figure}

For the BCSPWR09 matrix, we have selected a maximum diagonal block
size of 128 as the trade-off between load balance and minimum size of
the last block.  Figure~\ref{NT-6} depicts a reordered version of the
matrix presented in figure~\ref{NT-3} with diagonal blocks assigned to
a processor ordered to make those blocks contiguous.  The ordering of
the clustered matrix to reflect contiguous diagonal blocks within each
processor is performed by simple row and column exchanges, and has no
effect on either the positive definite nature of the matrix or the
number of fillin.  This is just a modification to the permutation
matrix used in equation~\ref{PAPxeqb}.  In this figure, the measure of
merit for the pigeon-hole load balancing algorithm is the number of
factorization (multiply and addition) operations.

\begin{figure}
\begin{center}
\begin{picture}(216,216)
\psfig{figure=bcs09-bdb128-4-fac.eps,width=216pt,height=216pt}
%\psfig{figure=TBD.eps,width=216pt,height=216pt}
\end{picture}
\caption{Load Balancing on Factoring --- BCSPWR09 --- 4 Processors \label{NT-6}}
\end{center}
\end{figure}

Similar clustering analyses have been performed for the BCSPWR10 load
flow matrices.  An example of the BCSPWR10 matrix after applying the
node-tearing algorithm to order into block-diagonal-bordered form is
presented in figure~\ref{NT-9}.  Figure~\ref{NT-9} illustrates the
block-diagonal-bordered form achievable when ordering the BCSPWR10
matrix with maximum block size of 256.  The number of graph nodes or
diagonal elements in the BCSPWR10 matrix is 5300 and this matrix has
476 equations in the lower border and last block.  This figure also
illustrates load balancing for four processors based on the number of
calculations to perform the factorization. Summary statistics are
presented in table~\ref{NT-12} for the two Boeing-Harwell matrices
with various maximum diagonal block sizes.

\begin{figure}
\begin{center}{}
\begin{picture}(216,216)
\psfig{figure=bcs10-bdb256-4-fac.eps,width=216pt,height=216pt}
%\psfig{figure=TBD.eps,width=216pt,height=216pt}
\end{picture}
\caption{Load Balancing on Factoring --- BCSPWR10 --- 4 Processors \label{NT-9}}
\end{center}
\end{figure}

\begin{table}
\begin{center}
\begin{tabular}{|l|r|r|r|r|r|r|r|} \hline
\multicolumn{1}{|c|}{}          & \multicolumn{1}{c|}{}          & \multicolumn{1}{c|}{}          & \multicolumn{1}{c|}{}          & \multicolumn{1}{c|}{}           & \multicolumn{1}{c|}{\bf Last}  & \multicolumn{1}{c|}{}           & \multicolumn{1}{c|}{$\triangle$} \\
\multicolumn{1}{|c|}{\bf Graph} & \multicolumn{1}{c|}{\bf Nodes} & \multicolumn{1}{c|}{\bf Edges} & \multicolumn{1}{c|}{$max_{DB}$} & \multicolumn{1}{c|}{\bf Fillin} & \multicolumn{1}{c|}{\bf Block} & \multicolumn{1}{c|}{\bf Factor} & \multicolumn{1}{c|}{\bf Solve} \\ \hline
BCSPWR09 & 1,723 & 4,117  &  32 &  2,765 & 190 & $O(N^{1.31})$ & $O(N^{1.21})$ \\
         &       &        &  64 &  3,248 & 153 & $O(N^{1.35})$ & $O(N^{1.22})$ \\
         &       &        & 128 &  3,486 & 128 & $O(N^{1.38})$ & $O(N^{1.23})$ \\ \hline
BCSPWR10 & 5,300 & 20,596 & 128 & 26,234 & 578 & $O(N^{1.51})$ & $O(N^{1.23})$ \\
         &       &        & 256 & 29,420 & 476 & $O(N^{1.55})$ & $O(N^{1.24})$ \\ 
         &       &        & 512 & 38,880 & 511 & $O(N^{1.62})$ & $O(N^{1.24})$ \\ \hline
\end{tabular}
\caption{Summary Statistics \label{NT-12}}
\end{center}
\end{table}

The BCSPWR10 matrix has substantially more edges per node than the
BCSPWR09 matrix, which results in a significantly greater
computational complexity for the factorization of the matrix.
Meanwhile, the computational complexity of the triangular solution
$(\triangle)$ is similar to the BCSPWR09 matrix.  Greater
computational complexity means more work, and if that work can be
distributed uniformly to all processors, then there is the potential
for greater parallel implementation efficiency.  Our goal is to
develop scalable algorithms for Choleski factorization; however,
parallelism in sparse Choleski algorithms is heavily dependent on the
available parallelism in the actual sparse matrix.  In
section~\ref{results-solve}, it will be shown that the additional
computations in BCSPWR10 (when compared to BCSPWR09) permits
reasonable efficiencies when factoring this matrix with four
processors.  Meanwhile, the solver for BCSPWR09 simply does not have
adequate calculations to overcome the effects of limited amounts of
calculations.

Further examination of table~\ref{NT-12} shows that the computational
complexity of the forward reduction and backward substitution phases
are similar in magnitude regardless of the choice of matrix.  It will
be difficult to obtain good speedup with an algorithm with
computational complexity so close to unity.

\subsection{Empirical Results --- Sparse Choleski Solver \label{results-solve}}

The performance of the block-diagonal-bordered Choleski algorithm was
tested with the following maximum diagonal block size for the
respective matrices:
\begin{itemize}
\item BCSPWR09 --- $max_{DB}=128$
\item BCSPWR10 --- $max_{DB}=256$
\end{itemize}
Examples of each of these block-diagonal-bordered sparse matrices were
presented in the previous section with load balancing for four
processors.  These ordered matrices have been used because they
offered the best performance for a particular load-flow matrix.
Performance of the multiprocessor algorithms are illustrated using
graphs with plots of relative speedup versus the number of processors.
Relative speedup has been defined above.

Graphs of speedup calculated from empirical performance data are
provided in figures~\ref{SP-1} and~\ref{SP-2}.  Each figure has a
family of three curves that show speedup for:
\begin{enumerate}
\item Choleski factorization 
\item forward reduction and backward substitution
\item factorization and a single forward reduction and backward substitution
\end{enumerate}
In general these speedup graphs illustrate that for the real power
system load-flow matrices from the Boeing-Harwell series, speedup
appears to be limited to approximately 2.7 regardless of the number of
processors used to solve the problem.  Careful analysis of the
performance data shows that the following data trends.

\begin{figure}
\begin{center}
\begin{picture}(216,144)
\psfig{figure=gplotout4.ps,width=216pt,height=144pt}
\end{picture}
\caption{Speedup for BCSPWR09 Data --- 2, 4, and 8 processors \label{SP-1}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\begin{picture}(216,144)
\psfig{figure=gplotout5.ps,width=216pt,height=144pt}
\end{picture}
\caption{Speedup for BCSPWR10 Data --- 2, 4, 8, and 16 processors \label{SP-2}}
\end{center}
\end{figure}

The time to factor, forward reduce, or backward substitute the data in
the last diagonal block dominates the calculations for these matrices
because of the low order of computational complexity in the diagonal
block partitions of the matrix.  These algorithms use pipelined
communications that {\em data-starve} quickly with the small sizes of
the last diagonal blocks.  On the other hand, solutions with different
values of the maximum diagonal block size yield larger last diagonal
blocks.  However, while better speedups are achieved, the actual run
times are not improved.

The performance of the block-diagonal-bordered sparse Choleski solver
using the BCSPWR10 matrix shows some improvement when compared to the
solver performance on the BCSPWR09 matrix.  The additional operations
in the BCSPWR10 matrix permit improvements in speedup and efficiencies
of 60\% for four processors, although speedup remains nearly constant
while additional processors are used.  

\subsection{Conclusions on Sparse Matrix Research \label{conc}}

In general, solving load-flow matrices from real-power systems has
proven to be a very difficult challenge for parallel sparse Choleski
solvers.  As we have illustrated, the computational complexity of the
actual data matrices is sufficiently low, even for matrix
factorization, that the speedup performance of the parallel
block-diagonal-bordered Choleski solver is limited by computational
starvation and in some instances, load imbalance.

In this paper we present research into parallel
block-diagonal-bordered sparse Choleski factorization algorithms
developed with special considerations to irregular sparse matrices
originating in the electrical power systems community.  Available
parallelism in the block-diagonal-bordered matrix structure offers
promise for simplified implementation and also offers a simple
decomposition of the problem into clearly identifiable subproblems.
Parallel block-diagonal-bordered Choleski solvers require a three step
preprocessing phase that is reusable for static matrices.

We developed an implementation that offered efficiencies of 60\% for
Choleski factorization with four processors, although the
implementation was not efficient beyond four processors with any of
the power system load-flow matrices examined.  Further examinations
into techniques to better solve the last diagonal block could
significant improve performance.  While a dense Choleski solver was
used in this implementation, due to the available parallelism in
block-diagonal-bordered form matrices, any technique can be used to
solve this sub-matrix, including iterative techniques.

\subsection*{Acknowledgments}
We thank Alvin Leung, Kamala Anupindi, Nancy McCracken, Paul
Coddington, and Tony Skjellum for their assistance in this research.
This work has been supported in part by Niagara Mohawk Power
Corporation, the New York State Science and Technology Foundation, the
NSF under co-operative agreement No. CCR-9120008, and ARPA under
contract \#DABT63-91-K-0005.

\newpage
\bibliography{paperIII}
\bibliographystyle{plain}

\end{document}
\endinput
