
\section{Parallel Direct Sparse Solver Research}

We have performed extensive research into parallel direct methods for
block-diagonal-bordered sparse matrices --- LU factorization and
Choleski factorization algorithms developed with special consideration
for irregular sparse matrices originating in the electrical power
systems community.  We have focused on developing parallel linear
solvers optimized for sparse matrices from the power systems
community --- in particular, we have examined linear solvers for
matrices resulting from power distribution system networks.  These
matrices are some of the most sparse matrices encountered in real-life
applications, and these matrices also are irregular.  Recently,
\cite{Kumar-scale,Kumar-int-pt} have reported scalable Choleski
solvers, but they are for matrices that have more row/columns, that
have more nonzero elements per row/column, and that are more regular
than power systems matrices.  When scalability of sparse linear
solvers is examined using real, irregular sparse matrices, the
available parallelism in the sparse matrix and load-imbalance overhead
can be as much the reason for poor parallel efficiency as the parallel
algorithm or implementation \cite{DPK-scale,Roth-SC94}.

\subsection{Block-Diagonal-Bordered Power System Matrices}

Power system distribution networks are generally hierarchical with
limited numbers of high-voltage lines transmitting electricity to
connected local networks that eventually distribute power to
customers.  In order to ensure reliability, highly interconnected
local networks are fed electricity from multiple high-voltage sources.
Electrical power grids have graph representations which in turn can be
expressed as matrices --- electrical buses are graph nodes and matrix
diagonal elements, while electrical transmission lines are graph edges
which can be represented as non-zero off-diagonal matrix elements.  We
show that it is possible to identify the hierarchical structure within
a power system matrix using only the knowledge of the interconnection
pattern by {\em tearing} the matrix into partitions and coupling
equations that yield a block-diagonal-bordered matrix.
Node-tearing-based partitioning identifies the basic network structure
that provides parallelism for the majority of calculations within the
direct solution of a linear system.

In this paper we examine the applicability of parallel direct
block-diagonal-bordered sparse solvers for real power system
applications that require either the solution of symmetric positive
definite sparse matrices or location symmetric sparse matrices that
result from solving problems relating to power systems networks.
Variations of this technique could be used to solve other power system
sparse linear systems such as those that result from solving
linearized differential-algebraic equations that result from transient
stability analysis or small-signal stability assessments.  The
implementations we describe in this paper work directly with the
equations resulting from the power systems network, the smallest class
of power systems matrix.

The implementations we developed can be used to solve symmetric
positive definite load flow analysis Jacobian matrices or position
symmetric network Jacobian matrices from transient stability analysis.
In spite of only examining direct linear solver implementations that
solve relatively small network-related matrices, we have been able to
obtain good parallel speedups --- as much as factors of 10 for
factorization of linear systems of complex variables.  We expect that
even better performances would be possible for parallel
implementations designed to solve a single system of linear equations
that represent a combination of the generator dynamical equations and
network equations in a linearized form of the differential-algebraic
equations from transient stability analysis or small-signal analysis.
For these problems, there is additional parallel calculations with no
additional parallel communications overhead.

\subsection{Block-Diagonal-Bordered Direct Linear Solvers}

Direct block-diagonal-bordered sparse linear solvers exhibit distinct
advantages when compared to general direct parallel sparse algorithms
for irregular matrices.  Task assignments for numerical factorization
on distributed-memory multi-processors depend only on the assignment
of data to blocks, and data communications are significantly reduced
with uniform and structured communications.  

Factorization algorithms for block-diagonal-bordered form matrices
require a specialized ordering step coupled to an explicit load
balancing step in order to generate this matrix form and to uniformly
distribute the computational workload for an irregular matrix
throughout a distributed-memory multi-processor.  This ordering
relates to more general sparse Choleski algorithms that use
elimination trees --- however, this algorithm develops an elimination
tree with only two-levels of supernodes.  Matrix orderings are
performed using a diakoptic technique based on node-tearing-nodal
analysis, with load balancing to optimize performance when factoring
the diagonal blocks and borders, the lower layer of supernodes in the
elimination tree.

Block-diagonal-bordered sparse matrix algorithms require modifications
to the normal preprocessing phase described in numerous papers on
parallel Choleski factorization \cite{gang4-3,Heath,BlockPar-3,icase2}.  
Each of these papers use the paradigm to {\em order} the sparse matrix
and then perform {\em symbolic factorization} in order to determine
the locations of all fillin values so that static data structures can
be utilized for maximum efficiency when performing numerical
factorization.  We modify this commonly used sparse matrix
preprocessing phase to include an explicit {\em load balancing} step
coupled to the {\em ordering} step so that the workload is uniformly
distributed throughout a distributed-memory multi-processor and
parallel algorithms make efficient use of the computational resources.

Parallel block-diagonal-bordered sparse linear solvers offer the
potential for regularity often absent from other parallel sparse
solvers \cite{DPK-1,DPK-par-chol,DPK-SC94,DPK-par-GS}.  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, and load-balancing is sufficiently
effective that relative speedups greater than ten have been observed
in empirical performance measurements for 32 processors on a Thinking
Machines CM-5 multi-processor.  In addition to the promising speedup
encountered for only the parallel direct linear solver, other
dimensions exist in electrical power system applications that can be
exploited to efficiently make use of multi-processors with greater
than 32 processors.

\subsection{Active Messages}

The algorithms we developed require active message remote procedure
calls in order to minimize communications overhead and obtain good
relative speedup.  Active message remote procedure calls (RPCs)
provide protocol-less access to the transport layer of the
communication network on the CM-5.  The user must assume
responsibilities for all aspects of communications --- as a result,
active message RPCs provide very low-latency communications.  Active
messages provide extremely fast interprocessor communications, and
also permit a parallel-code development paradigm that greatly
simplifies the implementation of these sparse matrix algorithms.
Empirical data has been collected on both active message-based
implementations and more traditional cooperative buffer-based message
passing commands in order to illustrate the need for low latency
communications when solving matrices that are as sparse and irregular
as power systems matrices.

Active messages are an example of a distributed-memory multiprocessor
message-passing paradigm that has been developed by the supercomputer
community to provide extremely low latency communications.  All other
CM-5 message-passing software is built upon this low latency
protocol-less remote procedure call.  The empirical data presented in
this paper clearly shows that low-latency communications are required
for direct linear solvers for power systems network matrices.  This
research has been performed on a small version (only 32 processors) of
a massively parallel processing (MPP) architecture, the Thinking
Machines CM-5.  We believe that scalable parallel processing (SPP)
architectures, like the IBM SP2, will eventually provide similar
low-latency communications for short messages because there are many
parallel algorithms that can only be implemented efficiently with this
type of interprocessor communications support.  SPP hardware
developers recognize that low-latency communications increase the
utility of their computer and, consequently, improve market potential.

\subsection{Embedded Software Applications}

Our research has examined the performance of a block-diagonal-bordered
direct solvers, with implementations of both Choleski and LU
factorization, to be incorporated within electrical power system
applications.  Because we are considering software to be embedded
within a more extensive application, we examine efficient parallel
forward reduction and backward substitution algorithms in addition to
parallel factorization algorithms.  Due to the reduced amount of
calculations in the triangular solution phases of solving a system of
factored linear equations, these algorithms are often ignored when
parallel Choleski or LU factorization algorithms are presented in the
literature.

In our research, we have found that the development of parallel
factorization algorithms must consider forward reduction and backward
substitution, because the choice of the order of calculations in
factorization can greatly influence the performance of the parallel
triangular solutions.  Data structures are dependent on the order of
calculations, in order to ensure cache coheriency, and the amount of
communications in parallel forward reduction and backward substitution
is dependent on the data layout.  We have found that the results of
additional communications overhead can eliminate any potential speedup
for parallel forward reduction with column oriented data storage.
Fortunately, the LU factorization algorithm can be implemented in a
manner to eliminate this communications overhead problem, although, it
cannot be eliminated for Choleski factorization, where either forward
reduction or backward substitution must be performed with an implicit
transpose of the factored matrix.

\subsection{Overview of NPAC Technical Report SCCS-679}

Our research into direct sparse linear solvers for power systems
applications is presented in {\em Parallel Direct Methods for
Block-Diagonal-Bordered Sparse Matrices}, NPAC technical report
SCCS-679 \cite{DPK-dir}. That paper is organized as follows.  In
section 2, we describe the electrical power system applications that
are the basis for this work.  In section 3, we briefly review direct
solution techniques for factorization and forward reduction/backward
substitution.  In this section, we review the literature concerning
general parallel LU and Choleski factorization algorithms.  This is
followed by a theoretical derivation of the available parallelism in
both the factorization and forward reduction/backward substitution
phases when solving block-diagonal-bordered form sparse matrices.
Paramount to exploiting the advantages of this parallel linear solver
is the process of ordering the irregular sparse power system matrices
into this form in a manner that balances the workload among
multi-processors.  In section 5, we describe the three-step
preprocessing phase used to generate matrix ordering for
block-diagonal-bordered matrices with uniformly distributed processing
load.  We present the node-tearing algorithm developed to order
matrices into block-diagonal-bordered form in section 6.  In section
7, we describe our block-diagonal-bordered sparse LU and Choleski
algorithms that has been implemented on the CM-5.  Analysis of the
performance of these ordering techniques are presented in section 8
for actual power system network matrices from the Boeing-Harwell
series, the Electrical Power Research Institute (EPRI), and the
Niagara Mohawk Power Corporation.  We present our conclusions
concerning parallel block-diagonal-bordered direct linear solvers for
electrical power system applications in section 9.

