
\section{Parallel Sparse Matrix Algorithms}

    Sparse matrix computations are some of the most challenging problems 
for parallel machines, but they are key to many power system applications,
including load flow and transient stability analysis.
Significant progress has already been made with concurrent sparse matrix
solvers.  Various organizations have efforts in progress to develop general 
sparse matrix solvers.  There has already been substantial work in
this area, especially by Michael Heath at Oak Ridge National Laboratory and
Ian Duff in Europe and groups at Sandia.
One of our graduate students is presently working on parallel sparse matrix
techniques especially suited for power systems applications. 
Another graduate student in our team has worked with a group at Lawrence 
Livermore National Laboratory on a related gas pipeline simulation.  
Several of these efforts are currently being monitored.  
During the study phase, parallel sparse matrix solvers will be developed that
interface with the differential-algebraic equation (DAE) solvers also being 
investigated.  

The in-house research at Syracuse University will develop specialized sparse 
matrix solvers that exploit hierarchical structure in the underlying problem.  
Sparse matrices result when solving the differential-algebraic equations that 
represent the generator
differential equations and the power distribution network when examining system
dynamics that result from a fault in a transient stability analysis.  
During the first study phase, in addition to monitoring work by Heath on 
general sparse matrix solvers, three competing techniques
will be examined to solve sparse power system matrices:
\begin{enumerate}
\item direct methods
\item iterative methods, and
\item composite direct/iterative methods.
\end{enumerate}
Critical to any of these approaches are techniques to reorder the linear sparse
matrices to exploit the underlying hierarchy.  For direct techniques, 
reordering enhances parallelism, by identifying independent portions of 
matrices that can be computed with minimal interprocessor communications. 
Also for direct techniques, optimal reordering will minimize the number of
calculations by limiting the added calculations due to fillin.  
For iterative techniques, such as waveform relaxation, reordering according to
the hierarchical nature of the underlying problem is required for rapid 
convergence.
Reordering can be used to improve load balancing for either direct or iterative
techniques.  

Diakoptics theory and the actual matrix sparsity structure (the power 
distribution network structure) will be used to 
develop algorithms to reorder the sparse matrices resulting from power system
applications into block-bordered diagonal form. 
Several layers of hierarchy may exist in
state or regional power networks.  
The differential-algebraic equations for a power system transient stability
problem can be represented in block bordered diagonal form.
There are distinct portions of the matrix that
represent the generator differential equations and the power distribution 
network.  By reordering the power distribution network, additional independent 
matrix partitions are available.  
After the matrix is reordered according to the underlying
structure, additional reordering may be required to more effectively balance 
the number of calculations.  Reordering to balance the load for individual 
processors often requires the trade-off where the number of equations in the 
border is expanded in the process of
equalizing the number of floating point operations in the independent blocks.  
Load balancing could result in slightly more calculations, however, the effect 
of additional calculations should be mitigated by a better distribution of all 
calculations over the multiple processors. 

	For pure direct methods, there are four sequential steps:
\begin{enumerate}
\item parallel factorization of independent blocks
\item parallel factorization of coupling equations 
\item parallel forward reduction, and 
\item parallel backward substitution.
\end{enumerate}
By reordering the matrix into block-bordered diagonal form, many of the 
operations that would be sequential in a general solver (for example, 
parallel forward reduction and parallel backward
substitution) are now highly parallel.  This will simplify the
development of scalable, transportable algorithms.  Iterative waveform 
relaxation techniques would solve the tightly coupled equations in the blocks 
first, then include the coupling equations.
Iterative techniques with rapid convergence could be
used concurrently in the independent blocks, while concurrent solutions of
the entire set of equations would require slower convergent algorithms.  
Composite direct/iterative
techniques could solve the equations in the blocks by direct methods, then 
iteratively update the answers to consider the interconnections.  
The direct solution of independent blocks can be viewed as a preconditioning 
step for the overall iterative solution.  The benefit of this technique is 
that fillin in the borders of the block-bordered diagonal
matrix is eliminated, possibly reducing the number of calculations.  


    Parallel sparse matrix solvers are critically dependent on a proper
match between particular problem and machine architecture and "general 
approaches" are often less than optimal.
It is essential to have matrices of representative sizes to successfully 
research  this real-life problem. 
We will work with Niagara Mohawk to obtain
the actual structure of the 
transmission network for three interconnected power grids:
\begin{enumerate}
\item Niagria Mohawk's network
\item the New York State Power Pool network, and
\item the Eastern US regional power network.
\end{enumerate}
It would not be surprising if our research indicates that different 
techniques are applicable for the different size networks. 
The problem size will also influence the applicable hardware.  

	We expect that transient stability analysis for 300 generators and a
2000 bus network system would perform
well on the 32 node MIMD computers that NPAC presently has. 
However, a "small," say 10 generator and 50 bus network, system would not
have enough parallelism to easily make use of existing parallel computers.
Meanwhile, the large 10 teraflop ($10^{13}$ floating point operations
per second) parallel supercomputers available in the year 2000 will probably
contain 10,000 nodes.  This will be too much parallelism for a simulation of
Niagara Mohawk's system but could be well suited for the transmission
system covering a large geographic area.  The coordination of grids over a
large region may be an exciting possibility opened up by the parallel
supercomputers of the future.  

	An additional source of parallelism in the power system
transient stability problem, is the examination of multiple, independent
fault events simultaneously.
Such a technique could significantly improve parallel computer speed-up and
efficiency by matching Niagria Mohawk's
requirements to available computer architectures.  The solution of the class of
concurrent problems where there is precedence in the order of the calculations,
such as the solution of linear equations, generally has
monotonically decreasing efficiency as the number of processors is increased.
However, transient stability analyses for independent faults can be processed 
totally independently.  To better describe this phenomenon, an example will be
presented.  

Presently, scalable general sparse matrix solvers developed by Heath yield
compounded speedups of 1.5 for every time the number of processors is doubled.
For a transient stability analysis on a 64 ($2^{6}$) processor MIMD parallel 
computer, eight independent faults can be analyzed using eight processors per
fault in the same time that only three independent faults can be analyzed if 
all 64 processors are dedicated to analyzing a single fault at a time. 
The use of all 64 processors, however, could yield the solution to a single 
transient stability fault in only 30\% of the time required for eight
processors, If time permits, it is cost effective to spread the multiple 
independent cases to separate groups of processors, and
exploit the requirement to examine multiple faults to enhance efficiency.
It is our belief that the specialized sparse matrix algorithms 
being developed for this effort can significantly improve upon the speedups 
reported in this example for a general sparse matrix solver.  

