The Design and Use of Algorithms for Permuting Large Entries to the Diagonal of Sparse Matrices.
We consider techniques for permuting a sparse matrix so that the diagonal of the permuted matrix has entries of large absolute value.  We discuss various criteria for this and consider their implementation as computer codes. We then indicate several cases where such a permutation can be useful.  These include the solution of sparse equations by a direct method and by an iterative technique.  We also consider its use in generating a preconditioner for an iterative method. We see that the effect of these reorderings can be dramatic although the best a priori strategy is by no means clear.
Introduction
We study algorithms for the permutation of a square unsymmetric sparse matrix A
of order n so that the diagonal of the permuted matrix has large entries. This can
be useful in several ways. If we wish to solve the system
where A is a nonsingular square matrix of order n and x and b are vectors of length
n, then a preordering to place large entries on the diagonal can be useful whether
direct or iterative methods are used for solution.
For direct methods, putting large entries on the diagonal suggests that pivoting
down the diagonal might be more stable. There is, of course, nothing rigorous in this
and indeed stability is not guaranteed. However, if we have a solution scheme like
the multifrontal method of Duff and Reid (1983), where a symbolic phase chooses
the initial pivotal sequence and the subsequent factorization phase then modifies
this sequence for stability, it can mean that the modification required is less than if
the permutation were not applied.
For iterative methods, simple techniques like Jacobi or Gauss-Seidel converge
more quickly if the diagonal entry is large relative to the off-diagonals in its row or
column and techniques like block iterative methods can benefit if the entries in the
diagonal blocks are large. Additionally, for preconditioning techniques, for example
for diagonal preconditioning or incomplete LU preconditioning, it is intuitively
evident that large diagonals should be beneficial.
We consider more precisely what we mean by such permutations in Section 2, and
we discuss algorithms for performing them and implementation issues in Section 3.
We consider the effect of these permutations when using direct methods of solution
in Section 4 and their use with iterative methods in Sections 5 and 6, discussing
the effect on preconditioning in the latter section. Finally, we consider some of the
implications of this current work in Section 7.
Throughout, the symbols jxj should be interpreted in context. If x is a scalar,
the modulus is intended; if x is a set, then the cardinality, or number of entries in
the set, is understood.
Permuting a matrix to have large diagonals
2.1 Transversals and maximum transversals
We say that an n \Theta n matrix A has a large diagonal if the absolute value of each
diagonal entry is large relative to the absolute values of the off-diagonal entries in
its row and column. We will be concerned with permuting the rows and columns
of the matrix so the resulting diagonal of the permuted matrix has this property.
That is, for the permuted matrix, we would like the ratio
to be large for all j, 1 - j - n. Of course, it is not even possible to ensure that
this ratio is greater than 1.0 for all j as the simple example
shows. It
is thus necessary to first scale the matrix before computing the permutation. An
appropriate scaling would be to scale the columns so that the largest entry in each
column is 1.0. The algorithm that we describe in Section 2.2 would then have the
effect of maximizing (2.1).
For an arbitrary nonsingular n \Theta n matrix, it is a necessary and sufficient
condition that for a set of n entries to be permuted to the diagonal, no two can
be in the same row and no two can be in the same column. Such a set of entries
is termed a maximum transversal, a concept that will be central to this paper and
which we now define more rigorously.
We let T denote a set of (at most n) ordered index pairs (i; j), 1 -
which each row index i and each column index j appears at most once. T is called
a transversal for matrix A, if a ij 6= 0 for each (i; . T is called a maximum
transversal if it has largest possible cardinality. jT j is equal to n if the matrix is
nonsingular. If indeed jT defines an n \Theta n permutation matrix P with
so that P T A is the matrix with the transversal entries on the diagonal.
In sparse system solution, a major use of transversal algorithms is in the
first stage of permuting matrices to block triangular form. The matrix is first
permuted by an unsymmetric permutation to make its diagonal zero-free after which
a symmetric permutation is used to obtain the block triangular form. An important
feature of this approach is that the block triangular form does not depend on which
transversal is found in the first stage (Duff 1977). A maximum transversal is also
required in the generalization of the block triangular ordering developed by (Pothen
and Fan 1990).
2.2 Bottleneck transversals
We will consider two strategies for obtaining a maximum transversal with large
transversal entries. The primary strategy that we consider in this paper is to
maximize the smallest value on the diagonal of the permuted matrix. That is, we
compute a maximum transversal T such that for any other maximum transversal T 1
we have
min
(i;j)2T
Transversal T is called a bottleneck transversal 1 , and the smallest value ja ij j, (i;
T , is called the bottleneck value of A. Equivalently, if jT smallest value
on the diagonal of P T A is maximized, over all permutations P, and equals the
bottleneck value of A.
An outline of an algorithm that computes a bottleneck transversal T 0 for a matrix
A is given below. We assume that we already have an algorithm for obtaining a
maximum transversal and denote by MT routine that returns a maximum
transversal for a matrix A, starting with the initial "guess" transversal T . We let
A ffl denote the matrix that is obtained by setting to zero in A all entries ja ij j for
which denote the transversal obtained by removing
from transversal T all the elements
Algorithm BT
Initialization:
Set fflmin to zero and fflmax to infinity.
while (there exist do
begin
choose
(We discuss how this is chosen later)
then
fflmin := ffl;
else
endif
Complete transversal for permutation;
(Needed if matrix structurally singular)
M is a maximum transversal for A, and hence jM j is the required cardinality
of the bottleneck transversal T 0 that is to be computed. If A is nonsingular, then
Throughout the algorithm, fflmax and fflmin are such that a maximum
transversal of size jM j does not exist for A fflmax but does exist for A fflmin . At each
step, ffl is chosen in the interval (fflmin; fflmax), and a maximum transversal for the
matrix A ffl is computed. If this transversal has size jM j, then fflmin is set to ffl,
1 The term bottleneck has been used for many years in assignment problems, for example
Glicksberg and Gross 1953)
otherwise fflmax is set to ffl. Hence, the size of the interval decreases at each step and
ffl will converge to the bottleneck value. After termination of the algorithm, T 0 is
the computed bottleneck transversal and ffl the corresponding bottleneck value. The
value for ffl is unique. The bottleneck transversal T 0 is not usually unique.
Algorithm BT makes use of algorithms for finding a maximum transversal. The
currently known algorithm with best asymptotic bound for finding a maximum
transversal is by Hopcroft and Karp (1973). It has a worst-case complexity of
O(
n- ), where - is the number of entries in the matrix. An efficient implementation
of this algorithm can be found in Duff and Wiberg (1988). The depth-first search
algorithm implemented by Duff (1981) in the Harwell Subroutine Library code MC21
has a theoretically worst-case behaviour of O(n- ), but in practice it behaves more
like O(n Because this latter algorithm is far simpler, we concentrate on this
in the following although we note that it is relatively straightforward to modify and
use the algorithm of Hopcroft and Karp (1973) in a similar way.
A limitation of algorithm BT is that it only maximizes the smallest value on
the diagonal of the permuted matrix. Although this means that the other diagonal
values are no smaller, they may not be maximal. Consider, for example, the 3 \Theta 3
ffiC A (2.2)
with ffi close to zero. Algorithm BT applied to this matrix returns
either the transversal f(1; 1); (2; 2); (3; 3)g or f(2; 1); (1; 2); (3; 3)g. Clearly, the latter
transversal is preferable. The modifications that we propose help to do this by
choosing large entries when possible for the early transversal entries.
It is beneficial to first permute the matrix to block triangular form and then to
use BT on only the blocks on the diagonal. This can be done since all entries in any
maximum transversal must lie in these blocks. Furthermore, not only does this mean
that BT operates on smaller matrices, but we also usually obtain a transversal of
better quality inasmuch as not only is the minimum diagonal entry maximized but
this is true for each block on the diagonal. Thus for matrix (2.2), the combination
of an ordering to block triangular form followed by BT would yield the preferred
transversal f(2; 1); (1; 2); (3; 3)g.
There are other possibilities for improving the diagonal values of the permuted
matrix which are not the smallest. One is to apply a row scaling subsequent to an
initial column scaling of the matrix A. This will increase the numerical values of all
the nonzero entries in those rows for which the maximum absolute numerical value
is less than one. A row scaling applied to the matrix (2.2) changes the coefficient
a 33 from ffi to 1:0, and now algorithm BT will compute f(2; 1); (1; 2); (3; 3)g as the
bottleneck transversal of the matrix (2.2). Unfortunately, such a row scaling does
not always help, as can be seen by the matrixB @
1:0 ffiC A
with the maximum transversals
all legitimate bottleneck transversals. Indeed the BT algorithm is very dependent on
scaling. For example, the matrix
has bottleneck transversal f(2; 1); (1; 2)g
whereas, if it is row scaled to
, the bottleneck transversal is f(1; 1); (2; 2)g.
Another possibility for improving the size of the diagonal values is to apply
algorithm BT repeatedly. Without loss of generality, suppose that, after application
of BT, entry a nn has the smallest diagonal value. Algorithm BT can then be
applied to the (n \Gamma 1) \Theta (n \Gamma 1) leading principal submatrix of A, and this could be
repeated until (after k steps) the (n \Gamma leading principal submatrix of
A only contains ones (on assumption original matrix was row and column scaled).
Obviously, this can be quite expensive, since algorithm BT is applied O(n) times
although we have a good starting point for the BT algorithm at each stage. We call
this algorithm the successive bottleneck transversal algorithm. Because of this and
the fact that we have found that it usually gives little improvement over BT, we do
not consider it further in this paper.
2.3 Maximum Product transversals
An algorithm yielding the same transversal independent of scaling is to maximize
the product of the moduli of entries on the diagonal, that is to find a permutation
oe so that
Y
a ioe i j (2.3)
is maximized. This is the strategy used for pivoting in full Gaussian elimination by
Olschowka and Neumaier (1996) and corresponds to obtaining a weighted bipartite
matching. Olschowka and Neumaier (1996) combine a permutation and scaling
strategy. The permutation, as in (2.3), maximizes the product of the diagonal entries
of the permuted matrix. (Clearly the product is zero if and only if the matrix is
structurally singular.) The scaling transforms the matrix into a so-called I-matrix,
whose diagonal entries are all one and whose off-diagonal entries are all less than or
equal to one.
Maximizing the product of the diagonal entries of A is equivalent to minimizing
the sum of the diagonal entries of a matrix that is defined as follows (we
here assume that denotes an n \Theta n nonnegative nonsingular matrix):
log a
where a is the maximum absolute value in column j of matrix A.
Minimizing the sum of the diagonal entries can be stated in terms of an
assignment problem and can be solved in O(n 3 ) time for full n \Theta n matrices or in
O(n- log n) time for sparse matrices with - entries. A bipartite weighted matching
algorithm is used to solve this problem. Applying this algorithm to C produces
vectors u, v and a transversal T , all of length n, such that
If we define
then, the scaled matrix is an I-matrix. We do not do this scaling
in our experiments but, unlike Olschowka and Neumaier, we use a sparse bipartite
weighted matching whereas they only considered full matrices.
The worst case complexity of this algorithm is O(n- log n). This is similar to BT,
although in practice it sometimes requires more work than BT. We have programmed
this algorithm, without the final scaling. We have called it algorithm MPD (for
Maximum Product on Diagonal) and compare it with BT and MC21 in the later
sections of this paper. Note that on the matrixB @
the MPD algorithm obtains the transversal f(1; 1); (2; 2); (3; 3)g whereas,
for example for Gaussian elimination down the diagonal, the transversal
would be better. Additionally, the fact that scaling does
influence the choice of bottleneck transversal could be deemed a useful characteristic.
3 Implementation of the BT algorithm
We now consider implementation details of algorithm BT from the previous section.
We will also illustrate its performance on some matrices from the Harwell-Boeing
Collection (Duff, Grimes and Lewis 1989) and the collection of Davis (1997). A code
implementing the BT algorithm will be included in a future release of the Harwell
Subroutine Library (HSL 1996).
When we are updating the transversal at stage (?) of algorithm BT, we can
easily accelerate the algorithm described in Section 2 by computing the value of the
minimum entry of the transversal, viz.
min
(i;j)2T
and then setting fflmin to this value rather than to ffl. The other issue, crucial
for efficiency, is the choice of ffl at the beginning of each step. If, at each step,
we choose ffl close to the value of fflmin then it is highly likely that we will find
a maximum transversal, but the total number of steps required to obtain the
bottleneck transversal can be very large. In the worst case, we could require
steps when the number of nonzero entries in A ffl reduces by only one at each iteration.
The algorithm converges faster if the size of the interval (fflmin; fflmax) reduces
significantly at each step. It would therefore appear sensible to choose ffl at each
step so that the interval is split into two almost equal subintervals, that is ffl -
(fflmin+fflmax)=2. However, if most of the nonzero values in A that have a magnitude
between fflmin and fflmax, are clustered near one of these endpoints, the possibility
exists that only a few nonzero values are discarded and the algorithm again will
proceed slowly. To avoid this, ffl should be chosen as the median of the nonzero
values between fflmin and fflmax.
We now consider how a transversal algorithm like MC21 can be modified to
implement algorithm BT efficiently. Before doing this, it is useful to describe briefly
how MC21 works. Each column of the matrix is searched in turn (called an original
column) and either an entry in a row with no transversal entry presently in the row
is found and this is made a transversal entry (a cheap assignment) or there is no
such entry and so the search moves to a previous column whose transversal entry is
in one of the rows with an entry in the original column. This new column is then
checked for a cheap assignment. If one exists, then this cheap assignment and the
entry in the original column in the row of the old transversal entry, replace that as
transversal entries thereby extending the length of the transversal by 1. If there is
no cheap assignment, then the search continues to other columns in a depth first
search fashion until a chain or augmenting path of the form
is found where there are no transversal entries in row i and every odd member of the
path is a transversal entry. The assignment is made in column j and the transversal
extended by 1 by replacing all transversal entries in the augmenting path with the
even members of this path.
Transversal selection algorithms like MC21 do not take into account the numerical
values of the nonzero entries. However, it is clear that the algorithm BT will
converge faster if T is chosen so that the value of its minimum entry is large. We
do this by noting that, when constructing an augmenting path, there are often
several candidates for a cheap assignment or for extending the path. MC21 makes an
arbitrary choice and we have modified it so the candidate with largest absolute value
is chosen. Note that this is a local strategy and does not guarantee that augmenting
paths with the highest values will be found.
The second modification aims at exploiting information obtained from previous
steps of algorithm BT. Algorithm BT repeatedly computes a maximum transversal
ffl ). The implementation of MC21 in the Harwell Subroutine Library
computes T from scratch, so we have modified it so that it can start with a partial
transversal. This can easily be achieved by holding the set of columns which contain
entries of the partial transversal and performing the depth search search through
that set of columns.
Of course, there are many ways to implement the choice of ffl. One alternative
is to maintain an array PTR (of length -) of pointers, such that the entries in the
first part of PTR point to those entries in A that form matrix A ffl max , the first two
parts of PTR point to the entries that form A ffl min
, and the elements in the third
part of PTR point to all the remaining (smaller) entries of A. A new value for ffl
can then be chosen directly (O(1) time) by picking the numerical value of an entry
that is pointed to by an element of the second part of PTR. After the assignment
in algorithm BT to either ffl min or ffl max , the second part of PTR has to be permuted
so that PTR again can be divided into three parts. An alternative is to do a global
(using a fast sorting algorithm) on all the entries of A, such that the elements
of PTR, point to the entries in order of decreasing absolute value. Then again PTR
can be divided into three parts as described in the previous alternative. By choosing
(in O(1) time) ffl equal to the numerical value of the entry pointed to by the median
element of the second part of PTR, ffl will divide the interval (ffl min
of close-to-equal size. Both alternatives have the advantage of being able to choose
the new ffl quickly, but require O(-) extra memory and (repeated) permutations of
the pointers.
We prefer an approach that is less expensive in memory and that matches our
transversal algorithm better. Since MC21 always searches the columns in order, we
facilitate the construction of the matrices A ffl , by first sorting the entries in each
column of the matrix A by decreasing absolute value. For a sparse matrix with
a well bounded number of entries in each column, this can be done in O(n) time.
The matrix A ffl is then implicitly defined by an array LEN of length n with LEN[j]
pointing to the first entry in column j of matrix A whose value is smaller than ffl,
which is the position immediately after the end of column j of matrix A ffl . Since the
entries of a column of A ffl are contiguous, the repeated modification of ffl by algorithm
BT, which redefines matrix A ffl , corresponds to simply changing the pointers in the
array LEN.
The actual choice of ffl at phase (?) in algorithm BT is done by selecting in
matrix A ffl min an entry that has an absolute value X such that ffl min
The columns of A ffl min
are searched until such an entry is found and ffl is set to its
absolute value. This search costs O(n) time since, for each column, we have direct
access to the entries with absolute values between ffl min and ffl max through the pointer
array LEN.
As mentioned before, by choosing ffl carefully, we can speed up algorithm BT
considerably. Therefore, instead of choosing an arbitrary entry from the matrix to
define ffl, we can choose a number (k say) of entries lying between ffl min and ffl max at
random, sort them by absolute value, and then set ffl to the absolute value of the
median element. 2 In our implementation we used
The set of matrices that we used for our experiments are unsymmetric matrices
taken from the sparse matrix collections Duff, Grimes and Lewis (1992) and Davis
(1997).

Table

3.1 shows the order, number of entries, and the time to compute a
bottleneck transversal for each matrix. All matrices are initially row and column
scaled. By this we mean that the matrix is scaled so that the maximum entry in
each row and in each column is one.
The machine used for the experiments in this and the following sections is a 166
MHz SUN ULTRA-2. The algorithms are implemented in Fortran 77.
Matrix n - Time in secs
GOODWIN 7320 324784 0.27 2.26 1.82
ONETONE2 36057 227628 2.63 0.53 0.42

Table

3.1: Times for transversal algorithms. Order of matrix is n and number of
entries - .
2 This is a technique commonly used to speed up sorting algorithms like quicksort.
4 The solution of equations by direct methods
MCSPARSE, a parallel direct unsymmetric linear system solver developed by
Gallivan, Marsolf and Wijshoff (1996), uses a reordering to identify a priori large and
medium grain parallelism and to reorder the matrix to bordered block triangular
form. Their ordering uses an initial nonsymmetric ordering that enhances the
numerical properties of the factorization, and subsequent symmetric orderings
are used to obtain a bordered block triangular matrix (Wijshoff 1989). The
nonsymmetric ordering is effectively a modified version of MC21. During each search
phase, for both a cheap assignment and an augmenting path, an entry a ij is selected
only if its absolute value is within a bound ff, 0 - ff - 1, of the largest entry in
column j. Instead of taking the first entry that is found by the search that satisfies
the threshold, the algorithm scans all of the column for the entry with the largest
absolute value.
The algorithm starts off with an initial bound ff = 0:1. If a maximum transversal
cannot be found, then the values in each column are examined to determine the
maximum value of the bound that would have allowed an assignment to take place
for that column. The new bound is then set to the minimum of the bound estimates
from all the failed columns and the algorithm is restarted. If a bound less than a
preset limit is tried and a transversal is still not found, then the bound is ignored
and the code finds any transversal. In our terminology (assuming an initial column
scaling of the matrix) this means that a maximum transversal of size n is computed
for the matrix A ff .
In the multifrontal approach of Duff and Reid (1983), later developed by Amestoy
and Duff (1989), an analysis is performed on the structure of A A T to obtain
an ordering that reduces fill-in under the assumption that all diagonal entries will
be numerically suitable for pivoting. The numerical factorization is guided by an
assembly tree. At each node of the tree, some steps of Gaussian elimination are
performed on a dense submatrix whose Schur complement is then passed to the
parent node in the tree where it is assembled (or summed) with Schur complements
from the other children and original entries of the matrix. If, however, numerical
considerations prevent us from choosing a pivot then the algorithm can proceed, but
now the Schur complement that is passed to the parent is larger and usually more
work and storage will be needed to effect the factorization.
The logic of first permuting the matrix so that there are large entries on the
diagonal, before computing the ordering to reduce fill-in, is to try and reduce the
number of pivots that are delayed in this way thereby reducing storage and work for
the factorization. We show the effect of this in Table 4.1 where we can see that even
using MC21 can be very beneficial although the BT algorithm can show significant
further gains and sometimes the use of MPD can cause further significant reduction
in the number of delayed pivots. We should add that the numerical accuracy of
the solution is sometimes slightly improved by these permutations and, in all cases,
good solutions were found.
Matrix Transversal algorithm used
None MC21 BT MPD
GOODWIN 536 1622 358 53

Table

4.1: Number of delayed pivots in factorization from MA41. An "-" indicates
that MA41 requires a real working space larger than 25 million words (of 8 bytes).
In

Table

4.2, we show the effect of this on the number of entries in the factors.
this mirrors the results in Table 4.1 and shows the benefits of the transversal
selection algorithms. This effect is seen in Table 4.3 where we can sometimes observe
a dramatic reduction in time for solution when preceded by a permutation.
Matrix Transversal algorithm used
None MC21 BT MPD
ONETONE2 14082683 2875603 2167523 2169903
GOODWIN 1263104 2673318 1791112 1282004

Table

4.2: Number of entries in the factors from MA41.
In addition to being able to select the pivots chosen by the analysis phase, the
multifrontal code MA41 will do better on matrices whose structure is symmetric
or nearly so. The transversal orderings in some cases increase the symmetry of the
resulting reordered matrix. This is particularly apparent when we have a very sparse
system with many zeros on the diagonal. In that case, the reduction in number of off-diagonal
entries in the reordered matrix has an influence on the symmetry. Notice
that, in this respect, the more sophisticated transversal algorithms may actually
cause problems since they could reorder a symmetrically structured matrix with a
zero-free diagonal, whereas MC21 will leave it unchanged.
Matrix Transversal algorithm used
None MC21 BT MPD
GOODWIN 3.64 14.63 6.00 3.56

Table

4.3: Time (in seconds on Sun ULTRA-2) for MA41 for solution of system.
5 The solution of equations by iterative methods
A large family of iterative methods, the so-called stationary methods, has the
iteration scheme
is a splitting of A, and M is chosen such that a system of the
is easy to solve. If M is invertible, (5.1) can be written as
We have
where ae is the spectral radius, so that, if jjM convergence of the
iterates x (k) to the solution A \Gamma1 b is guaranteed for arbitrary x (0) . In general, the
smaller jjM the faster the convergence. Thus an algorithm that makes
entries in M large and those in N small should be beneficial.
The most simple method of this type is the Jacobi method, corresponding to the
splitting denotes the diagonal, L the strictly
lower triangular part, and U the strictly upper triangular part of the matrix A.
However, this is not a particularly current or powerful method so we conduct our
experiments using the block Cimmino implementation of Arioli, Duff, Noailles and
Ruiz (1992), which is equivalent to using a block Jacobi algorithm on the normal
equations. In this implementation, the subproblems corresponding to blocks of rows
from the matrix are solved by a direct method similar to that considered in the
previous section. For similar reasons, it can be beneficial to increase the magnitude
of the diagonal entries through unsymmetric permutations.
We show the effect of this in Table 5.1, where we see that the number of iterations
for the solution of the problem MAHINDAS 7682). The convergence
tolerance was set to 10 \Gamma12 . The transversal selection algorithm was followed by a
reverse Cuthill McKee algorithm to obtain a block tridiagonal form. The matrix
was partitioned in 2, 4, 8, and 16 block rows and the acceleration used was a block
CG algorithm with block sizes of 1, 4, and 8.
Acceleration
# block rows None MC21 BT MPD

Table

5.1: Number of iterations of block Cimmino algorithm on MAHINDAS.
In every case, the use of a transversal algorithm accelerates the convergence of
the method, sometimes by a significant amount. However, the use of the algorithms
to increase the size of the diagonal entries does not usually help convergence further.
The convergence of block Cimmino depends on angles between subspaces which is
not so strongly influenced by the diagonal entries.
6 Preconditioning
In this section, we consider the effect of using a permutation induced by our
transversal algorithms prior to solving a system using a preconditioned iterative
method. We consider preconditionings corresponding to incomplete factorizations
of the form ILU(0), ILU(1), and ILUT and study the convergence of the iterative
methods GMRES(20), BiCGSTAB, and QMR. We refer the reader to a standard
text like that of Saad (1996) for a description and discussion of these methods. Since
the diagonal of the permuted matrix is "more dominant" than the diagonal of the
original matrix, we would hope that such permutations would enhance convergence.
We show the results of some of our runs in Table 6.1. The maximum number
of iterations was set to 1000 and the convergence tolerance to 10 \Gamma9 . It is quite
clear that the reorderings can have a significant effect on the convergence of the
preconditioned iterative method. In some cases, the method will only converge after
the permutation, in others it greatly improves the convergence. It would appear
from the results in Table 6.1 and other experiments that we have performed, that
the more sophisticated MPD transversal algorithm generally results in the greatest
reduction in the number of iterations, although the best method will depend on the
overall solution time including the transversal selection algorithm.
7 Conclusions and future work
We have described algorithms for obtaining transversals with large entries and have
indicated how they can be implemented showing that resulting programmes can be
written for efficient performance.
While it is clear that reordering matrices so that the permuted matrix has a
large diagonal can have a very significant effect on solving sparse systems by a wide
range of techniques, it is somewhat less clear that there is a universal strategy that is
best in all cases. We have thus started experimenting with combining the strategies
mentioned in this paper and, particularly for the block Cimmino approach, with
combining our unsymmetric ordering with a symmetric ordering. One example that
we plan to study is a combination with the symmetric TPABLO ordering (Benzi,
Choi and Szyld 1997).
It is possible to extend our techniques to orderings that try to increase the size
of not just the diagonal but also the immediate sub and super diagonals and then
use the resulting tridiagonal part of the matrix as a preconditioner.
One can also build other criteria into the weighting for obtaining a bipartite
matching, for example, to incorporate a Markowitz count so that sparsity would
also be preserved by the choice of the resulting diagonal as a pivot.
Finally, we noticed in our experiments with MA41 that one effect of transversal
selection was to increase the structural symmetry of unsymmetric matrices. We are
thus exploring further the use of ordering techniques that more directly attempt to
increase structural symmetry.

Acknowledgments

We are grateful to Patrick Amestoy of ENSEEIHT, Michele Benzi of CERFACS, and
Daniel Ruiz of ENSEEIHT for their assistance with the experiments on the direct
methods, the preconditioned iterative methods, and the block iterative methods
respectively. We would also like to thank Alex Pothen for some early discussions on
bottleneck transversals, and John Reid and Jennifer Scott for comments on a draft
of this paper.
Matrix and method Transversal algorithm
BiCGSTAB 123 21 11
QMR 101 26 17
QMR 72 19 12
MAHINDAS
WEST0497

Table

6.1: Number of iterations required by some preconditioned iterative methods.



--R



Threshold ordering for preconditioning nonsymmetric problems


The design and use of a frontal scheme for solving sparse unsymmetric equations



Users' guide for the Harwell-Boeing sparse matrix collection (Release I)
A production line assignment problem





Iterative methods for sparse linear systems
Symmetric orderings for unsymmetric sparse matrices
--TR

--CTR
Iain S. Duff , Jennifer A. Scott, Stabilized bordered block diagonal forms for parallel sparse solvers, Parallel Computing, v.31 n.3+4, p.275-289, March/April 2005
Kai Shen, Parallel sparse LU factorization on second-class message passing platforms, Proceedings of the 19th annual international conference on Supercomputing, June 20-22, 2005, Cambridge, Massachusetts
Olaf Schenk , Klaus Grtner, Two-level dynamic scheduling in PARDISO: improved scalability on shared memory multiprocessing systems, Parallel Computing, v.28 n.2, p.187-197, February 2002
Olaf Schenk , Andreas Wchter , Michael Hagemann, Matching-based preprocessing algorithms to the solution of saddle-point problems in large-scale nonconvex interior-point optimization, Computational Optimization and Applications, v.36 n.2-3, p.321-341, April     2007
Olaf Schenk , Klaus Grtner, Solving unsymmetric sparse systems of linear equations with PARDISO, Future Generation Computer Systems, v.20 n.3, p.475-487, April 2004
Patrick R. Amestoy , Iain S. Duff , Jean-Yves L'excellent , Xiaoye S. Li, Analysis and comparison of two general sparse solvers for distributed memory computers, ACM Transactions on Mathematical Software (TOMS), v.27 n.4, p.388-421, December 2001
Kai Shen, Parallel sparse LU factorization on different message passing platforms, Journal of Parallel and Distributed Computing, v.66 n.11, p.1387-1403, November 2006
Xiaoye S. Li, An overview of SuperLU: Algorithms, implementation, and user interface, ACM Transactions on Mathematical Software (TOMS), v.31 n.3, p.302-325, September 2005
Anshul Gupta, Recent advances in direct methods for solving unsymmetric sparse systems of linear equations, ACM Transactions on Mathematical Software (TOMS), v.28 n.3, p.301-324, September 2002
Anwar Hussein , Ke Chen, Fast computational methods for locating fold points for the power flow equations, Journal of Computational and Applied Mathematics, v.164-165 n.1, p.419-430, 1 March 2004
Jack Dongarra , Victor Eijkhout , Piotr uszczek, Recursive approach in sparse matrix LU factorization, Scientific Programming, v.9 n.1, p.51-60, January 2001
Belur V. Dasarathy, Editorial: Identity fusion in unsupervised environments, Information Fusion, v.7 n.2, p.157-160, June, 2006
Xiaoye S. Li , James W. Demmel, SuperLU_DIST: A scalable distributed-memory sparse direct solver for unsymmetric linear systems, ACM Transactions on Mathematical Software (TOMS), v.29 n.2, p.110-140, June
Nicholas I. M. Gould , Jennifer A. Scott , Yifan Hu, A numerical evaluation of sparse direct solvers for the solution of large sparse symmetric linear systems of equations, ACM Transactions on Mathematical Software (TOMS), v.33 n.2, p.10-es, June 2007
Michele Benzi, Preconditioning techniques for large linear systems: a survey, Journal of Computational Physics, v.182 n.2, p.418-477, November 2002
sparse matrices;iterative methods;maximum transversal;preconditioning;direct methods