\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{Parallel Sparse Matrix Algorithms}
\label{sec:Sparse} 
% \label{sparse}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Since power networks are sparsely connected (each bus is connected to
only a few others) simulations of power systems generally
require the solution of sparse linear equations,
which is done using sparse matrix techniques.
Sparse matrix solvers have complex and irregular data structures,
and are therefore much more difficult to parallelize than dense matrix
solvers, which are very regular, and thus easily distributed among 
multiple processors. The communication of data between processors in a
parallel sparse matrix solver is also much more complicated, 
and the amount of computation required is much less, than for a parallel
dense solver. Solving a dense $N \! \times \! N$ matrix takes order $N^3$
calculations, whereas for sparse matrices this is generally less than 
order $N^2$. Power systems matrices are extremely sparse, with
computational complexities as low as $N^{1.3}$.

Obtaining substantial speed-ups for a parallel sparse
matrix solver has thus proven to be an extremely difficult task, and
can usually only be achieved by using problem-specific methods that 
identify and effectively utilize any structure present in the sparse matrices.
General methods, although they exist, are usually less than optimal.
Power networks tend to be hierarchical in structure, and this can be
utilized to order the matrices into a structured 
{\it block-diagonal-bordered} form, where non-zero data elements occur
only in the borders and in blocks along the diagonal of the matrices
(see Fig.~{BBDM}).
Parallel block-diagonal-bordered sparse linear
solvers offer the potential for regularity often absent from other
parallel sparse solvers.  
We have performed a theoretical examination of the parallelism available
in block-diagonal-bordered matrices and developed software to
order power systems matrices into this form.

In this portion of the report, we review the factorization and
ordering of sparse matrices.  We then discuss the
available parallelism in block-diagonal-bordered form matrices and 
describe a three-step preprocessing phase we have developed to order
sparse power systems matrices into block-diagonal-bordered form and 
distribute the work-load evenly among processors ({\it load balancing}).
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.

For power system applications, 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 parallel sparse matrix solver.
This fact will be evident in the preliminary empirical data we have
collected on the Thinking Machines CM-5 parallel computer on which we
have developed our programs. 


\subsection{LU and Choleski Factorization \label{Choleski}}

Sparse matrix equations are solved using either {\it direct} or 
{\it iterative} methods. We have concentrated on implementing direct methods,
since these are currently used in power systems programs such as ETMSP.
However we are also studying the possible use of iterative methods,
which may be more effective for parallel processing.

The standard direct methods used to solve matrix equations are 
Choleski factorization for symmetric matrices (which occur in many power 
systems problems, such as calculating power flows in a network),
and LU factorization for non-symmetric matrices (which occur in transient
stability analysis).

We are considering the direct solution of the linear system
\begin{equation}
Ax=b,
\end{equation}
where $A$ is an $N \! \times \! N$ sparse matrix. 
This is done using the standard method of LU factorization \cite{Duff}.
The sparse matrix $A$ can be numerically factored into a product of
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.
The above equation is solved by setting $Ux=y$, and
substituting $y$ for $Ux$.  The numerical
solution for the resulting equation $Ly=b$ can be 
found by forward reduction, and the numerical solution for $x$ is
calculated by backward substitution in the equation $Ux=y$.
Forward reduction and backward substitution utilize the fact that
triangular linear systems can be readily solved numerically by solving
for the first value and substituting that value into subsequent equations,
in a recursive manner \cite{Duff}.

If the matrix $A$ is a symmetric positive definite sparse matrix,
then a special form of LU factorization known as 
Choleski factorization can be used,
that exploits the symmetry and inherent numerically stable
characteristics of this matrix form \cite{Duff}.  
Symmetric positive definite matrices occur in power system load flow analysis.
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, 
\end{equation}
%
where $L^{T}$ is the transpose of $L$, i.e. $L^{T}(i,j) = L(j,i)$, and 
is therefore an upper triangular matrix.
The equations are then solved by forward reduction and backward substitution
in the same way as for LU factorization.
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}.


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

The most significant aspect of parallel sparse matrix 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 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
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.  
In that case, the diagonal blocks can be factored completely independently
\cite{DPK-par-dir}.  
Provided that the matrix can be clustered into block-diagonal-bordered form,
then the parallel sparse matrix algorithm can reap additional
benefits, such as the elimination of task graphs for
distributed-memory multi-processor implementations. 
% DEFINE TASK GRAPH??
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 matrix solvers
\cite{shar-mem}.

Figure~\ref{BBDM} illustrates both the 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 U}$. The gray areas are those containing non-zero 
elements. This figure depicts four diagonal blocks, 
and the processors (P1, P2, P3, and P4) assigned to each portion of the
matrix data are listed with the data block.
The last block (at bottom right) is distributed over all the processors,
and solved using standard parallel dense matrix techniques.

\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{Ordering and Load Balancing \label{3step}}

A power network can be looked at as a mathematical {\it graph}, 
with the buses being the {\it nodes} of the graph 
and the lines connecting these buses being the {\it edges} of the graph.
Symmetric sparse matrices can be represented by graphs with the non-zero
elements corresponding to edges in the graph \cite{Duff,Heath}.
In the power systems case, this corresponds to the fact that every pair 
of connected buses has an associated equation for the power flowing 
between them, and this coupling is represented by an element of the matrix.

In sequential sparse matrix solvers, an {\it ordering} of the matrix
is done before the factorization.
Ordering a symmetric sparse matrix is actually little more
than changing the labels associated with nodes in the graph,
however, this simple task can drastically reduce the amount of
calculations involved when factoring a sparse matrix by affecting the
order of calculation and consequently the amount of {\it fillin}  
(the elements of the matrix that are initially zero, but become non-zero 
during factorization). 
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 \cite{Duff,Heath}. 

In parallel sparse matrix solvers, ordering can also be used to 
put the matrices into a more regular and structured form so that they
can more easily be distributed over processors, in such a way as to 
balance the computational load between processors.
We have used minimum degree ordering in conjunction with a
node-tearing-based ordering technique to order the sparse matrices
into block-diagonal-bordered form \cite{DPK-par-dir}.

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. 

Ordering requires first running a {\it clustering} routine that 
examines the connectivity of the data and splits it into elements that
can be placed in the diagonal blocks, and elements that are placed in the
last block. This involves various tunable parameters, such as the size
of the last block and the maximum size of the diagonal blocks.
For this report, optimizing over these parameters was performed by hand --
various values of the parameters were examined and the
block-diagonal-bordered form sparse matrix with the best load balance
and least numbers of operations was 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 we have used in performance studies can vary
from $O(N^{1.1})$ to $O(N^{2.0})$.  The computational complexity for
factoring these entire sparse matrices varies 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{enumerate}
\item factor independent blocks,
\item update the last block using data from the borders,
\item factor the last block.
\end{enumerate}
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 dense
matrix algorithm to factor the last diagonal block.


\subsection{Performance 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 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
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 speed-up and sources of parallel processing
overhead.  The block-diagonal-bordered Choleski solver has been run on
this test matrix using one to sixteen processors of the
Thinking Machines CM-5 parallel computer at NPAC. 
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 one of the power systems matrices we have studied
(the BCSPWR09 matrix presented in sections \ref{results-cluster}
and \ref{results-solve}),
although it has a substantially higher computational complexity,
$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}
\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 speed-up of each distinct portion of the
algorithm.  Data was collected in such a manner as not to impact the
overall measures of performance.  Performance of the multiprocessor
algorithms are illustrated using graphs plotting speed-up
versus the number of processors.

A graph of speed-up 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 speed-up 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 speed-up graph shows that the block-diagonal-bordered
Choleski solver has significant potential for efficient operations,
because speed-ups of over 11 were obtained with 16 processors.
This equates to a processor utilization efficiency of 70\%.  When
examining the runtime timing data for each processor, near perfect
speed-up (speed-up equal to the number of processors) is not achieved
because of overhead that occurs during those portions 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 Choleski factorization of the last block
\item the update of data in the $b$-vector 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.
\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 speed-up calculated from the empirical performance data is
presented in Fig.~\ref{TEST-perf-2}.  In this graph, the 
factorization of the diagonal blocks shows perfect speed-up, which is
not unexpected because there is no communications in this step.
However, the measured speed-up 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.
\label{TEST-perf-2}}
\end{center}
\end{figure}

Unfortunately, the calculations in the diagonal blocks for the
triangular solution are of a lesser computational complexity than
the calculations in the factorization step.  Consequently, poorer
parallel performance (lower speed-up) would be expected in forward
reduction and backward substitution than can be obtained in
block-diagonal-bordered factorization.
Parallel performance of those sections of the triangular solution algorithms
that require pipelined communications will also have less speed-up potential
than the triangular solution of data in the diagonal blocks.


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, all-to-all communications in both 
      factorization and forward reduction when
      sums of matrix-matrix products or vector-matrix
      products from data in the lower border are used to update the last
      diagonal block or the part of the $b$-vector associated with the
      last diagonal block.
\end{enumerate}
In both instances, the pipelined communications appear to 
{\em data-starve}, that is, there is not enough calculation to 
hide the cost of communications.
This is not surprising in the forward reduction of the last block, but
this also occurs in the factorization of the last block.
Experiments with adjusting the block size offered little improvement.
Further research into other asynchronous communications
techniques such as active messages \cite{Active-1} may be able to
significantly improve performance in this area.

We have discussed the effects on performance of communications overhead,
however other sources of parallel processing overhead are possible \cite{Fox}.
They include indexing overhead and load imbalance overhead. 
This analysis has been performed on data that had no load imbalance.
The sources of load imbalance overhead are a function of the interconnections
in the data. In the preprocessing phase, we attempt to minimize load 
imbalance in our parallel algorithm, however we will show in section
\ref{results-solve} that even with attempts at load balancing, 
this source of overhead cannot be entirely removed.

Indexing overhead is the extra computation 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 solver, thus the costs incurred in this stage are not addressed here. 


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

The performance of our block-diagonal-bordered Choleski solver has been
analyzed using power systems matrices taken from Boeing-Harwell set of
sample sparse matrix data \cite{BH-guide}, 
as well as data for the Niagara Mohawk power network obtained from 
the NMPC Operations group.
We used two matrices from the Boeing-Harwell series, which represent
real electrical power system networks from the Western and Eastern 
US respectively:  
\begin{itemize}
\item BCSPWR09 -- 1723 nodes and 4117 edges in the graph
\item BCSPWR10 -- 5300 nodes and 13571 edges in the graph
\end{itemize}

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 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 diagonal block size ($max_{DB}$),
defines how the node-tearing algorithm clusters the matrix into blocks.  
An example of applying the node-tearing algorithm to the
BCSPWR09 matrix is presented in figure~\ref{NT-3} for a maximum diagonal
block size of 128 nodes.  Original non-zero 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 3838, with a substantial
portion of the fillin occurring in the lower right-hand corner, or the
last diagonal-block.  The Boeing-Harwell matrix BCSPWR09 is very similar
to data obtained from the Niagara Mohawk Operations group.
That matrix has 1766 nodes and 4272 edges in the graph. 
Examples of applying the node-tearing clustering algorithm to
Niagara Mohawk data is presented in figure~\ref{nimo}.

\begin{figure}
\begin{center}
\begin{picture}(216,216)
% \psfig{figure=bcs09-bdb128.eps,width=216pt,height=216pt}
\end{picture}
\caption{Node-tearing-based ordering for BCSPWR09 with $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}
\end{picture}
\caption{Node-tearing-based ordering for Niagara Mohawk Operations Data 
with $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. 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 for the last
diagonal block even more desirable 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, 
since this is not as efficient to solve in parallel,
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 orderings 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{Number of nodes in the last block for node-tearing-based
ordering of 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{Load imbalance for factoring BCSPWR09 using node-tearing-based
ordering. 
\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 the diagonal blocks assigned to
a processor reordered 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.  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}
\end{picture}
\caption{Load balancing for factoring BCSPWR09 on 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.  This figure also
illustrates load balancing for four processors based on the number of
calculations to perform the factorization. 
Statistics for the reordered Boeing-Harwell matrices with various 
maximum diagonal block sizes are presented in table~\ref{NT-12}.

\begin{figure}
\begin{center}{}
\begin{picture}(216,216)
% \psfig{figure=bcs10-bdb256-4-fac.eps,width=216pt,height=216pt}
\end{picture}
\caption{Load balancing for factoring BCSPWR10 on 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{Results for clustering and ordering of Boeing-Harwell matrices.
\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) permit
reasonable efficiencies when factoring this matrix with four
processors.  However for BCSPWR09 the solver 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 extremely difficult to obtain good speed-up for an algorithm with
computational complexity so close to unity.


\subsection{Results for 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 is illustrated using
graphs with plots of speed-up versus the number of processors.

Graphs of speed-up 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 speed-up 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 speed-up graphs illustrate that for the real power
system load-flow matrices from the Boeing-Harwell series, speed-up
appears to be very limited, regardless of the number of processors
used to solve the problem (i.e. the algorithm is not scalable).  

\begin{figure}
\begin{center}
\begin{picture}(216,144)
% \psfig{figure=gplotout4.ps,width=216pt,height=144pt}
\end{picture}
\caption{Speed-up for BCSPWR09 Data.
\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{Speed-up for BCSPWR10 Data. 
\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.  Solutions with smaller values of the maximum
diagonal block size yield larger last diagonal blocks. 
However, while better speed-ups are achieved, the actual run
times are not improved.

The performance of the block-diagonal-bordered sparse Choleski solver
shows much better speed-up for the larger BCSPWR10 matrix, due to 
the increased amount of computation required. This means that parallel
processing can be more effectively used if the size of the network is 
large.



\subsection{The ETMSP Sparse Matrix Solver}

We are planning to interface our parallel sparse matrix solver into ETMSP,
to replace the sequential solver currently used.
ETMSP uses an in-house sparse matrix factorization package based on LU 
factorization. This package provides five main functions \cite{ETMSP}:
%
   \begin{itemize}
   \item Matrix reordering
   \item Symbolic factorization
   \item LU factorization
   \item Forward reduction and backward substitution
   \item Interface and support functions
   \end{itemize}
%
Minimum degree ordering is used to order the matrix once before the
first factorization, after which reordering is done only when there 
are significant changes to the matrix structure.

Symbolic factorization can be used to discover the fillin pattern of
the matrix during the factorization, so that static data structures 
for the non-zero elements can be set up to optimize memory access
(this is similar to the pseudo-factorization for our parallel solver). 
However ETMSP does not use a symbolic factorization -- it merely
leaves memory space available for adding fillin as it is generated.
This method saves time for problems where the structure of the matrix
does not change much during the calculation, which is true for the
power systems problems.

Standard LU factorization, forward reduction and backward substitution
are used. The interface function provides interface services for ETMSP
to access the sparse package. Additional functions to store and access
the solver's state are also provided. Our work to replace the ETMSP
sequential solver with a parallel sparse matrix solver mainly requires
modifying these interface routines.



\subsection{Conclusions and Recommendations}
\label{sec:SparseConc} 

We have developed a parallel block-diagonal-bordered sparse Choleski solver 
and implemented it on the Thinking Machines CM-5 parallel computer
using a message-passing paradigm.  
Preliminary results show the speed-up that can be obtained from 
solving load flow problems for power systems is very limited,
since the amount of computation required is so low that communications
overhead seriously degrades performance. Load imbalance is also a problem.

Transient stability analysis, however, uses LU factorization, and this 
should give a larger speed-up, since it requires more calculation.
A speed-up of only 5 is all that is required of a parallel sparse matrix 
solver to get an order of magnitude speed-up in ETMSP, since the generator
equations are completely independent and should give perfect speed-up
(see section \ref{sec:ParTSA}). 
We are currently working on extending the methods descibed here 
to produce a parallel
block-diagonal-bordered LU solver for transient stability analysis, 
that we believe will achieve this level of performance.
The speed-up is also dependent on the size of the problem -- 
simulations of larger networks, for example the whole of New York state
plus neighboring areas, would obtain much larger speed-ups than would
simulating only the Niagara Mohawk network. 

Using DAE solvers for the transient stability problem produces larger
matrices with additional available calculations, due to the inclusion 
of the generator equations.
These matrices are naturally in block-diagonal-bordered form.
The increased size and computational complexity of these matrices should
allow much greater speed-up for parallel LU factorization.
We believe parallel sparse matrix algorithms can be developed that can
give an order of magnitude speed-up for transient stability analysis
using a parallel DAE solver.

Tasks remaining in this work include the development and implementation of two
additional parallel block-diagonal-bordered sparse matrix solvers --
a direct solver based on LU factorization and an iterative Gauss-Seidel 
solver. Iterative techniques may be better suited to parallel processing.
These programs will be ported to the IBM SP-1 parallel computer,
and to a standardized portable parallel language that can be run on any 
distributed memory parallel computer.


\newpage
\bibliography{sparse}
\bibliographystyle{plain}

\end{document}
\endinput
