On Algorithms For Permuting Large Entries to the Diagonal of a Sparse Matrix.
We consider bipartite matching algorithms for computing permutations of a sparse matrix so that the diagonal of the permuted matrix has entries of large absolute value. We discuss various strategies for this and consider their implementation as computer codes. We also consider scaling techniques to further increase the relative values of the diagonal entries. Numerical experiments show the effect of the reorderings and the scaling on the solution of sparse equations by a direct method and by preconditioned iterative techniques.
Introduction
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. Permuting large nonzero entries onto the diagonal of a sparse matrix 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 of this kind can be useful whether direct or iterative methods are used
for solution (see Olschowka and Neumaier (1996) and Duff and Koster (1997)).
The work in this report is a continuation of the work reported by Duff and Koster
(1997) who presented an algorithm that maximizes the smallest entry on the diagonal
and relies on repeated applications of the depth first search algorithm MC21 (Duff 1981)
in the Harwell Subroutine Library (HSL 1996). In this report, we will be concerned with
other bipartite matching algorithms for permuting the rows and columns of the matrix so
that the diagonal of the permuted matrix is large. The algorithm that is central to this
report computes a matching that corresponds to a permutation of a sparse matrix such
that the product (or sum) of the diagonal entries is maximized. This algorithm is already
mentioned and used in Duff and Koster (1997), but is not fully described. In this report,
we describe the algorithm in more detail. We also consider a modified version of this
algorithm to compute a permutation of the matrix that maximizes the smallest diagonal
entry. We compare the performance of this algorithm with that of Duff and Koster (1997).
We also investigate the influence of scaling of the matrix. Scaling can be used before or
after computation of the matching to make the diagonal entries even larger relative to
the off-diagonals. In particular, we look at a sparse variant of a bipartite matching and
scaling algorithm of Olschowka and Neumaier (1996) that first maximizes the product of
the diagonal entries and then scales the matrix so that these entries are one and all other
entries are no greater than one.
The rest of this report is organized as follows. In Section 2, we describe some concepts
of bipartite matching that we need for the description of the algorithms. In Section 3,
we review the basic properties of algorithm MC21. MC21 is a relatively simple algorithm
that computes a matching that corresponds to a permutation of the matrix that puts as
many entries as possible onto the diagonal without considering their numerical values.
The algorithm that maximizes the product of the diagonal entries is described in Section
4. In Section 5, we consider the modified version of this algorithm that maximizes the
smallest diagonal entry of the permuted matrix. In Section 6, we consider the scaling of the
reordered matrix. Computational experience for the algorithms applied to some practical
problems and the effect of the reorderings and scaling on direct and iterative methods of
solution are presented in Sections 7 to 7.2. The effect on preconditioning is also discussed.
Finally, we consider some of the implications of this current work in Section 8.
matching
be a general n \Theta n sparse matrix. With matrix A, we associate a bipartite
graph E) that consists of two disjoint node sets V r and V c and an edge set
E, where (u; v) 2 E implies that . The sets V r and V c have cardinality n and
correspond to the rows and columns of A respectively. Edge (i; only if a ij 6= 0.
We define the sets ROW
c . These sets correspond to the positions of the entries in row i and column j of the
sparse matrix respectively. We use both to denote the absolute value and to signify
the number of entries in a set, sequence, or matrix. The meaning should always be clear
from the context.
A subset M ' E is called a matching (or assignment) if no two edges of M are incident
to the same node. A matching containing the largest number of edges possible is called a
maximum cardinality matching (or simply maximum matching). A maximum matching
is a perfect matching if every node is incident to a matching edge. Obviously, not every
bipartite graph allows a perfect matching. However, if the matrix A is nonsingular, then
there exists a perfect matching for GA . A perfect matching M has cardinality n and
defines an n \Theta n permutation matrix
so that both PA and AP are matrices with the matching entries on the (zero-free) diagonal.
Bipartite matching problems can be viewed as a special case of network flow problems (see,
for example, Ford Jr. and Fulkerson (1962)).
The more efficient algorithms for finding maximum matchings in bipartite graphs make
use of augmenting paths. Let M be a matching in GA . A node v is matched if it is incident
to an edge in M . A path P in GA is defined as an ordered set of edges in which successive
edges are incident to the same node. A path P is called an M-alternating path if the
edges of P are alternately in M and not in M . An M-alternating path P is called an M -
augmenting path if it connects an unmatched row node with an unmatched column node.
In the bipartite graph in Figure 2.1, there exists an M-augmenting path from column node
8 to row node 8. The matching M (of cardinality 7) is represented by the thick edges. The
black entries in the accompanying matrix correspond to the matching and the connected
matrix entries to the M-augmenting path. If it is clear from the context which matching
M is associated with the M-alternating and M-augmenting paths, then we will simply
refer to them as alternating and augmenting paths.
Let M and P be subsets of E. We define
If M is a matching and P is an M-augmenting path, then M \Phi P is again a matching,
and jM \Phi is an M-alternating cyclic path, i.e., an alternating path
whose first and last edge are incident to the same node, then M \Phi P is also a matching
and jM \Phi

Figure

2.1: Augmenting path
r1368425791368In the sequel, a matching M will often be represented by a pointer array
Augmenting paths in a bipartite graph G can be found by constructing alternating
trees. An alternating tree subgraph of G rooted at a row or column
node and each path in T is an M-alternating path. An alternating tree rooted at a
column node j 0 can be grown in the following way. We start with the initial alternating
tree (;; fj 0 g; ;) and consider all the column nodes j 2 T c in turn. Initially . For
each node j, we check the row nodes i 2 COL(j) for which an alternating path from i
to j 0 does not yet exist. If node i is already matched, we add row node i, column node
to T . If i is not matched, we extend T by row node i and
edge (and the path in T from node i to the root forms an augmenting path). A key
observation for the construction of a maximum or perfect matching is that a matching M
is maximum if and only if there is no augmenting path relative to M .
Alternating trees can be implemented using a pointer array c such that,
given an edge (i; is either the root node of the tree, or the edges
are consecutive edges in an alternating path towards the root.
Augmenting paths in an alternating tree (provided they exist) can thus easily be obtained
from p and m.
Alternating trees are not unique. In general, one can construct several alternating
trees starting from the same root node that have equal node sets, but different edge
sets. Different alternating trees in general will contain different augmenting paths. The
matching algorithms that we describe in the next sections impose different criteria on the
order in which the paths in the alternating trees are grown in order to obtain augmenting
paths and maximum matchings with special properties.
Matching
The asymptotically fastest currently known algorithm for finding a maximum matching is
by Hopcroft and Karp (1973). It has a worst-case complexity of O(
the number of entries in the sparse matrix. An efficient implementation of this algorithm
can be found in Duff and Wiberg (1988). The algorithm MC21 implemented by Duff
(1981) has a theoretically worst-case behaviour of O(n- ), but in practice it behaves more
like O(n Because this latter algorithm is simpler, we concentrate on this in the
following although we note that it is relatively straightforward to use the algorithm of
Hopcroft and Karp (1973) in a similar way to how we will use MC21 in later sections.
MC21 is a depth-first search algorithm with look-ahead. It starts off with an empty
matching M , and hence all column nodes are unmatched initially. See Figure 3.1. For
each unmatched column node j 0 in turn, an alternating tree is grown until an augmenting
path with respect to the current matching M is found (provided one exists). A set B is
used to mark all the matched row nodes that have been visited so far. Initially,
First, the row nodes in COL(j 0 ) are searched (look-ahead) for an unmatched node i 0 .
If one is found, the singleton path is an M-augmenting path. If there is
no such unmatched node, then an unmarked matched node i
is marked, the nodes i 0 and
, and the edges (i are added to
the alternating tree (by setting
The search then continues with column node
. For node j 1 , the row nodes in COL(j 1 ) are first checked for an unmatched node.
If one exists, say then the path forms an augmenting
path. If there is no such unmatched node, a remaining unmarked node i 1 is picked from
is set to j 1 ,
, and the search moves to node j 2 .
This continues in a similar (depth-first search) fashion until either an augmenting path
(with nodes j 0 and i k unmatched) or
until for some k ? 0, COL(j k ) does not contain an unmarked node. In the latter case,
MC21 backtracks by resuming the search at the previously visited column node j k\Gamma1 for
some remaining unmarked node i 0
Backtracking for
if MC21 resumes the search at column node j 0 and COL(j 0 ) does not contain an unmarked
node, then an M-augmenting path starting at node j 0 does not exist. In this case, MC21
continues with the construction of a new alternating tree starting at the next unmatched
column node. (The final maximum matching will have cardinality at most n \Gamma 1 and hence
will not be perfect.)

Figure

3.1: Outline of MC21.
do
repeat
if there exists i 2 COL(j) and i is unmatched then
else
if there exists
else
until iap 6= null or
if iap 6= null then augment along path from node iap to node j 0
end for
Weighted matching
In this section, we describe an algorithm that computes a matching for permuting a sparse
matrix A such that the product of the diagonal entries of the permuted matrix is maximum
in absolute value. That is, the algorithm determines a matching that corresponds to a
permutation oe that maximizes n
Y
This maximization multiplicative problem can be translated into a minimization
additive problem by defining matrix
log a
where a is the maximum absolute value in column j of matrix A. Maximizing
(4.1) is equal to minimizing
log
log a oe i
log ja ioe i
(log a oe i
log ja ioe i
Minimizing (4.2) is equivalent to finding a minimum weight perfect matching in an edge
weighted bipartite graph. This is known in literature as the bipartite weighted matching
problem or (linear sum) assignment problem in linear programming and combinatorial
optimization. Numerous algorithms have been proposed for computing minimum weight
perfect matchings, see for example Burkard and Derigs (1980), Carpaneto and Toth (1980),
Carraresi and Sodini (1986), Derigs and Metz (1986), Jonker and Volgenant (1987), and
Kuhn (1955). A practical example of an assignment problem is the allotment of tasks to
in the cost matrix C represents the cost or benefit of assigning person i
to task j.
be a real-valued n \Theta n E) be the
corresponding bipartite graph each of whose edges (i; . The weight
of a matching M in GC , denoted by c(M ), is defined by the sum of its edge weights, i.e.,
(i;j)2M
A perfect matching M is said to be a minimum weight perfect matching if it has smallest
possible weight, i.e., c(M) - c(M 0 ), for all possible maximum matchings M 0 .
The key concept for finding a minimum weight perfect matching is the so-called shortest
augmenting path. An M-augmenting path P starting at an unmatched column node j is
called shortest if c(M \Phi P possible M-augmenting paths P 0
starting at node j. We define
as the length of alternating path P . A matching M is called extreme if and only if it does
not allow any alternating cyclic path with negative length.
The following two relations hold. First, a perfect matching has minimum weight if it
is extreme. Second, if matching M is extreme and P is a shortest M-augmenting path,
then M \Phi P is extreme also. The proof for this goes roughly as follows. Suppose M \Phi P is
not extreme. Then there exists an alternating cyclic path Q such that c((M \Phi P
is extreme, there must exist a
subset forms an M-augmenting path and is shorter than P . Hence, P is
not a shortest M-augmenting path. This contradicts the supposition.
These two relations form the basis for many algorithms for solving the bipartite
weighted matching problem: start from any (possibly empty) extreme matching M and
successively augment M along shortest augmenting paths until M is maximum (or perfect).
In the literature, the problem of finding a minimum weight perfect matching is often
stated as the following linear programming problem. Find matrix
minimizing X
subject to X
If there is a solution to this linear program, there is one for which x ij 2 f0; 1g and there
exists a permutation matrix X such that 1g is a minimum weight perfect
matching (Edmonds and Karp 1972, Kuhn 1955). Furthermore, M has minimum weight
if and only if there exist dual variables u i and v j with
Using the reduced weight matrix
the reduced weight c(M) of matching M equals
the reduced length l(P ) of any M-alternating path P equals
(i;j)2P nM
and if M \Phi P is a matching, the reduced weight of M \Phi P equals
Thus, finding a shortest augmenting path in graph GC is equivalent to finding an
augmenting path in graph G C , with minimum reduced length.
edge contains no alternating paths P with negative length,
leading subpath P 0 of P .
Shortest augmenting paths in a weighted bipartite graph E) can be
obtained by means of a shortest alternating path tree. A shortest alternating path tree T
is an alternating tree each of whose paths is a shortest path in G. For any node i
we define d i as the length of the shortest path in T from node i to the root node (d
if no such path exists). T is a shortest alternating path tree if and only if d
for every edge (i; nodes i, j,
An outline of an algorithm for constructing a shortest alternating path tree rooted at
column node j 0 is given in Figure 4.1. Because the reduced weights c ij are non-negative,
and graph G C contains no alternating paths with negative length, we can use a sparse
variant of Dijkstra's algorithm (Dijkstra 1959). The set of row nodes is partitioned into
three sets B, Q, and W . B is the set of (marked) nodes whose shortest alternating paths
and distances to node j 0 are known. Q is the set of nodes for which an alternating path
to the root is known that is not necessarily the shortest possible. W is the set of nodes for
which an alternating path does not exist or is not known yet. (Note that since W is defined
implicitly as V r n (B [Q), it is not actually used in Figure 4.1.) The algorithm starts with
shortest alternating tree and extends the tree until an augmenting
path is found that is guaranteed to be a shortest augmenting path with respect to the
current matching M . Initially, the length of the shortest augmenting path lsap in the tree
is set to infinity, and the length of the shortest alternating path lsp from the root to any
node in Q is set to zero. On each pass through the main loop, another column node j is
chosen that is closest to the root j 0 . Initially
Each row node i 2 COL(j) whose shortest alternating path to the root is not known
yet (i 62 B), is considered. If P j 0 !j!i , the shortest alternating path from the root node
0 to node j (with length lsp) extended by edge (i; j) from node j to node i (with length
longer than the tentative shortest augmenting path in the tree (with length lsap),
then there is no need to modify the tree. If P j 0 !j!i has length smaller than lsap, and i
is unmatched, then a new shorter augmenting path has been found and lsap is updated.
If i is matched and P j 0 !j!i is also shorter than the current shortest alternating path to
(with length d i ), then a shorter alternating path to node i has been found and the tree
is updated, d i is updated, and if node i has not been visited previously, i is moved to Q.
Next, if Q is not empty, a node i 2 Q is determined that is closest to the root. Since all
weights c ij in the bipartite graph are non-negative, there cannot be any other alternating
path to node i that is shorter than the current one. Node i is marked (by adding it to
B), and the search continues with column node j This continues until there are
no more column nodes to be searched or until no new augmenting path can be
found whose length is smaller than the current shortest one (line lsap - lsp).
The original Dijkstra algorithm (intended for dense graphs) has O(n 2 ) complexity.
For sparse problems, the complexity can be reduced to O(- log n) by implementing the
set Q as a k-heap in which the nodes i are sorted by increasing distance d i from the root
(see for example Tarjan (1983) and Gallo and Pallottino (1988)). The running time of
the algorithm is dominated by the operations on the heap Q of which there are O(n)
delete operations, O(n) insert operations, and O(-) modification operations (these are
necessary each time a distance d i is updated). Each insert and modification operation
runs in O(log k n) time, a delete operation runs in O(k log k n) time. Consequently, the
algorithm for finding a shortest augmenting path in a sparse bipartite graph has run time
n) and the total run time for the sparse bipartite weighted algorithm is
n). If we choose 2, the algorithm uses binary heaps and we obtain a
time bound of O(n(- +n) log 2 n). If we choose we obtain a bound
of O(n- log -=n n).
The implementation of the heap Q is similar to the implementation proposed in Derigs
and Metz (1986). Q is a pair (Q is an array that contains all the row
nodes for which the distance to the root is shortest (lsp), and
By separating the nodes in Q that are closest to the root, we may reduce the number
of operations on the heap, especially in those situations where the cost matrix C has
only few different numerical values and many alternating paths have the same length.
Deleting a node from Q for which d i is smallest (see Figure 4.1), now consists of choosing
an (arbitrary) element from Q 1 . If Q 1 is empty, then we first move all the nodes in Q 2
that are closest to the root to Q 1 .
After the augmentation, the reduced weights c ij have to be updated to ensure that
alternating paths in the new G have non-negative length. This is done by modifying the

Figure

4.1: Construction of a shortest augmenting path.
while true do
dnew
if dnew ! lsap then
unmatched then
lsap := dnew; isap :=
else
choose
if lsap - lsp then exit while-loop;
if lsap 6= 1 then augment along path from node isap to node j
dual vectors u and v. If is the shortest alternating path tree that was
constructed until the shortest augmenting path was found, then u i and v j are updated as
The updated dual variables u and v satisfy (4.3) and the new reduced weights c ij are
non-negative.
The running time of the weighted matching algorithm can be decreased considerably
by means of a cheap heuristic that determines a large initial extreme matching M . We
use the strategy proposed by Carpaneto and Toth (1980). We calculate
Inspecting the sets COL(j) for each column node j in turn, we determine a large initial
matching M of edges for which for each remaining unmatched
column node j, every node i 2 COL(j) is considered for which and that is
matched to a column node other than j, say j . If an unmatched row node
can be found for which c i
in M is replaced
by (i; having repeated this for all unmatched columns, the search for
shortest augmenting paths starts with respect to the current matching.
Finally, we note that the above weighted matching algorithm can also be used for
maximizing the sum of the diagonal entries of matrix A (instead of maximizing the product
of the diagonal entries). To do this, we again minimize (4.2), but we redefine matrix C as
0; otherwise:
Maximizing the sum of the diagonal entries is equal to minimizing (4.2), since
a oe i
(a oe i
5 Bottleneck matching
We describe a modification of the weighted bipartite matching algorithm from the previous
section for permuting rows and columns of a sparse matrix A such that the smallest ratio
between the absolute value of a diagonal entry and the maximum absolute value in its
column is maximized. That is, the modification computes a permutation oe that maximizes
min
1-i-n
a oe i
where a j is the maximum absolute value in column j of the matrix A. Similarly to the
previous section, we transform this into a minimization problem. We define the matrix
a j
Then maximizing (5.1) is equal to minimizing
1-i-n
a oe i
a oe i
1-i-n
Given a matching M in the bipartite graph E), the bottleneck value of
M is defined as
(i;j)2M
The problem is to find a perfect (or maximum) bottleneck matching M for which c(M) is
minimal, i.e. c(M) - c(M 0 ), for all possible maximum matchings M 0 . A matching M is
called extreme if it does not allow any alternating cyclic path P for which c(M \PhiP
The bottleneck algorithm starts off with any extreme matching M . The initial
bottleneck value b is set to c(M ). Each pass through the main loop, an alternating tree
is constructed until an augmenting path P is found for which either c(M \Phi P
or c(M \Phi P small as possible. The initializations and the main loop
for constructing such an augmenting path are those of Figure 4.1. Figure 5.1 shows the
inner-loop of the weighted matching algorithm of Figure 4.1 modified to the case of the
bottleneck objective function. The main differences are that the sum operation on the path
lengths in Figure 4.1 is replaced by the "max" operation and, as soon as an augmenting
path P is found whose length lsap is less than or equal to the current bottleneck value
b, the main loop is exited, P is used to augment M , and b is set to max(b; lsap). The
bottleneck algorithm does not modify the edge weights c ij .
Similarly to the implementation discussed in Section 4, the set Q is implemented as a
now the array Q 1 contains all the nodes whose distance to the root is less
than or equal to the tentative bottleneck value b. Q 2 contains the nodes whose distance
to the root is larger than the bottleneck value but not infinity. Q 2 is again implemented
as a heap.

Figure

5.1: Modified inner loop of Figure 4.1 for the construction of a bottleneck
augmenting path.
dnew :=
if dnew ! lsap then
unmatched then
lsap := dnew; isap :=
if lsap - b then exit while-loop;
else
A large initial extreme matching can be found in the following way. We define
as the smallest entry in row i and column j, respectively. A lower bound b 0 for the
bottleneck value is
An extreme matching M can be obtained from the edges (i; j) for which c ij - b 0 ; we scan
all nodes in turn and for each node i 2 COL(j) that is unmatched and for which
is added to M . Then, for each remaining unmatched column node j,
every node i 2 COL(j) is considered for which c ij - b and that is matched to a column
node other than j, say j . If an unmatched row node i 1 2 COL(j 1 ) can be
found for which c i
is replaced by (i; having
done this for all unmatched columns, the search for shortest augmenting paths starts with
respect to the current matching.
Other initialization procedures can be found in the literature. For example, a slightly
more complicated initialization strategy is used by Finke and Smith (1978) in the context
of solving transportation problems. For every use
as the number of admissible edges incident to row node i and column node j respectively.
The idea behind using g i and h j is that once an admissible edge (i; j) is added to M , all
the other admissible edges that are incident to nodes i and j are no longer candidates
to be added to M . Therefore, the method tries to pick admissible edges such that the
number of admissible edges that become unusable is minimal. First, a row node i with
minimal g i is determined. From the set ROW (i) an admissible entry (i; (provided one
exists) is chosen for which h j is minimal and (i; j) is added to M . After deleting the edges
and the edges (k; j), k 2 COL(j), the method repeats the same for
another row node i 0 with minimal g i 0 . This continues until all admissible edges are deleted
from the graph.
Finally, we note that instead of maximizing (5.1) we also could have maximized the
smallest absolute value on the diagonal. That is, we maximize
min
1-i-n
and define the matrix C as
Note that this problem is rather sensitive to the scaling of the matrix A. Suppose for
example that the matrix A has a column containing only one nonzero entry whose absolute
value v is the smallest absolute value present in A. Then, after applying the bottleneck
algorithm, the bottleneck value b will be equal to this small value. The smallest entry on
the diagonal of the permuted matrix is maximized, but the algorithm did not have any
influence on the values of the other diagonal values. Scaling the matrix prior to applying
the bottleneck algorithm avoids this.
In Duff and Koster (1997), a different approach is taken to obtain a bottleneck
matching. Let A ffl denote the matrix that is obtained by setting to zero in A all entries
a ij for which ja ij denote the matching obtained by removing
from matching M all the entries (i; j) for which ja ij Throughout the
algorithm, fflmax and fflmin are such that a maximum matching 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 matching for the matrix A ffl is computed using a variant of MC21. If
this matching has size jM j, then fflmin is set to ffl, 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, M 0 is the computed bottleneck matching and ffl the
corresponding bottleneck value.
6 Scaling
Olschowka and Neumaier (1996) use the dual solution produced by the weighted matching
algorithm to scale the matrix. Let u and v be such that they satisfy relation (4.3). If we
define
then we have
Equality holds when that is (i; In words, D 1 AD 2 is a matrix whose
diagonal entries are one in absolute value and whose off-diagonal entries are all less than
or equal to one. Olschowka and Neumaier (1996) call such a matrix an I-matrix and use
this in the context of dense Gaussian elimination to reduce the amount of pivoting that is
needed for numerical stability. The more dominant the diagonal of a matrix, the higher
the chance that diagonal entries are stable enough to serve as pivots for elimination.
For iterative methods, the transformation of a matrix to an I-matrix is also of interest.
For example, from Gershgorin's theorem we know that the union of all discs
contains all eigenvalues of the n \Theta n matrix A. Disc K i has center at a ii and radius that is
equal to the sum of the absolute off-diagonal values in row i. Since the diagonal entries of
an I-matrix are all one, all the n disks have center at 1. The estimate of the eigenvalues
will be sharper as A deviates less from a diagonal matrix. That is, the smaller the radii of
the discs, the better we know where the eigenvalues are situated. If we are able to reduce
the radii of the discs of an I-matrix, i.e. reduce the off-diagonal values, then we tend to
cluster the eigenvalues more around one. In the ideal case, all the discs of an I-matrix
have a radius smaller than one, in which case the matrix is strictly row-wise diagonally
dominant. This guarantees that many types of iterative methods will converge (in exact
even simple ones like the Jacobi and Gauss-Seidel method. However, if at
least one disc remains with radius larger than or close to one, zero eigenvalues or small
eigenvalues are possible.
A straightforward (but expensive) attempt to decrease large off-diagonal entries of a
matrix is by row and column equalization (Olschowka and Neumaier 1996). Let A be
an I-matrix. We define matrix simplicity we assume
that A contains no zero entries.) Equalization consists of repeatedly equalizing the largest
absolute value in row i and the largest absolute values in column i:
for k := do
for to n do
For
and thus, if we define d 1
the algorithm minimizes the
largest off-diagonal absolute value in matrix D 1 AD 2 . The diagonal entries do not change.
Note that the above scaling strategy does not guarantee that all off-diagonal entries
of an I-matrix will be smaller than one in absolute value, for example if the I-matrix A
contains two off-diagonal entries a kl and a lk , k 6= l, whose absolute values are both one.
7 Experimental results
In this section, we discuss several cases where the reorderings algorithms from the previous
section 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.
The set of matrices that we used for our experiments are unsymmetric matrices taken
from the Harwell-Boeing Sparse Matrix Test Collection (Duff, Grimes and Lewis 1992)
and from the sparse matrix collection at the University of Florida (Davis 1997).
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 computer used for the experiments is a SUN UltraSparc with 256 Mbytes of main
memory. The algorithms are implemented in Fortran 77.
We use the following acronyms. MC21 is the matching algorithm from the Harwell
Subroutine Library for computing a matching such that the corresponding permuted
matrix has a zero free-diagonal (see Section 3). BT is the bottleneck bipartite matching
algorithm from Section 5 for permuting a matrix such that the smallest ratio between
the absolute value of a diagonal entry and the maximum absolute value in its column is
maximized. BT' is the bottleneck bipartite matching algorithm from Duff and Koster
(1997). MPD is the weighted matching algorithm from Section 4 and computes a
permutation such that the product of the diagonal entries of the permuted matrix
is maximum in absolute value. MPS is equal to the MPD algorithm, but after the
permutation, the matrix is scaled to an I-matrix (see Section 6).

Table

7.1 shows for some large sparse matrices the order, number of entries, and
the time for the algorithms to compute a matching. The times for MPS are not listed,
because they are almost identical to those for MPD. In general, MC21 needs the least time
to compute a matching, except for the ONETONE and TWOTONE matrices. For these
matrices, the search heuristic that is used in MC21 (a depth-first search with look-ahead)
does not perform well. This is probably caused by the ordering of the columns (variables)
and the entries inside the columns of the matrix. A random permutation of the matrix
prior to applying MC21 might lead to other results. There is not a clear winner between
the bottleneck algorithms BT and BT', although we note that BT' requires the entries
inside the columns to be sorted by value. This sorting can be expensive for relatively
dense matrices. MPD is in general the most expensive algorithm. This can be explained
by the more selective way in which this algorithm constructs augmenting paths.
7.1 Experiments with a direct solution method
For direct methods, putting large entries on the diagonal suggests that pivoting down the
diagonal might be more stable. Indeed, stability can still not be guaranteed, but 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.
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

Table

7.1: Times (in seconds) for matching algorithms. Order of matrix is n and
number of entries - .
Matrix
GOODWIN 7320 324784 0.27 2.26 4.17 1.82
(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 7.2 where we can see that even using MC21 can be very
beneficial although the other algorithms can show significant further gains.
In

Table

7.3, we show the effect of this on the number of entries in the factors. Clearly
this mirrors the results in Table 7.2.
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. Here, we define the structural symmetry for a matrix A as the number of entries a ij for
which a ji is also an entry, divided by the total number of entries. The structural symmetry
after the permutations is shown in Table 7.4. The matching orderings in some cases
increase the symmetry of the resulting reordered matrix, which 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 matching algorithms
may actually cause problems since they could reorder a symmetrically structured matrix
with a zero-free diagonal, whereas MC21 will leave it unchanged.

Table

7.2: Number of delayed pivots in factorization from MA41. An "-" indicates
that MA41 needed more than 200 MBytes of memory.
Matrix Matching algorithm used
None MC21 BT MPD MPS
GOODWIN 536 1622 427 53 41

Table

7.3: Number of entries (10 3 ) in the factors from MA41.
Matrix Matching algorithm used
None MC21 BT MPD MPS
ONETONE2 14,083 2,876 2,298 2,170 2,168
GOODWIN 1,263 2,673 2,058 1,282 1,281

Table

7.4: Structural symmetry after permutation.
Matrix Matching algorithm used
None MC21 BT MPD/MPS
GEMAT11
Finally, Table 7.5 shows the effect on the solution times of MA41. We sometimes
observe a dramatic reduction in time for the solution when preceded by a permutation.

Table

7.5: Solution time required by MA41.
Matrix Matching algorithm used
None MC21 BT MPD MPS
GOODWIN 3.64 14.63 7.98 3.56 3.56
Our implementations of the algorithms described in this paper have been used
successfully by Li and Demmel (1998) to stabilize sparse Gaussian elimination in a
distributed-memory environment without the need for dynamic pivoting. Their method
decomposes the matrix into an N \Theta N block matrix A[1 : by using the notion of
unsymmetric supernodes (Demmel, Eisenstat, Gilbert, Li and Liu 1995). The blocks are
mapped cyclically (in both row and column dimensions) onto the nodes (processors) of a
two-dimensional rectangular processor grid. The mapping is such that at step k of the
numerical factorization, a column of processors factorizes the block column A[k : N; k], a
row of processes participates in the triangular solves to obtain the block row U
and all processors participate in the corresponding multiple-rank update of the remaining
The numerical factorization phase in this method does not use (dynamic) partial
pivoting on the block columns. This allows for the a priori computation of the nonzero
structure of the factors, the distributed data structures, the communication pattern, and a
good load balancing scheme, which makes the factorization more scalable on distributed-memory
machines than factorizations in which the computational and communication
tasks only become apparent during the elimination process. To ensure a solution that
is numerically stable, the matrix is permuted and scaled before the factorization to
make the diagonal entries large compared to the off-diagonal entries, any tiny pivots
encountered during the factorization are perturbed, and a few steps of iterative refinement
are performed during the triangular solution phase if the solution is not accurate enough.
Numerical experiments demonstrate that the method (using the implementation of the
MPS algorithm) is as stable as partial pivoting for a wide range of problems.
7.2 Experiments with iterative solution methods
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.
7.2.1 Preconditioning by incomplete factorizations
In incomplete factorization preconditioners, pivots are often taken from the diagonal and
fill-in is discarded if it falls outside a prescribed sparsity pattern. (See Saad (1996) for
an overview.) Incomplete factorizations are used so that the resulting factors are more
economical to store, to compute, and to solve with.
One of the reasons why incomplete factorizations can behave poorly is that pivots
can be arbitrarily small (Benzi, Szyld and van Duin 1997, Chow and Saad 1997). Pivots
may even be zero in which case the incomplete factorization fails. Small pivots allow
the numerical values of the entries in the incomplete factors to become very large, which
leads to unstable and therefore inaccurate factorizations. In such cases, the norm of the
residual
U will be large. (Here, -
L and -
U denote the computed incomplete
A way to improve the stability of the incomplete factorization, is to preorder the
matrix to put large entries onto the diagonal. Obviously, a successful factorization still
cannot be guaranteed, because nonzero diagonal entries may become very small (or even
zero) during the factorization, but the reordering may mean that zero or small pivots
are less likely to occur. Table 7.6 shows some results for the reorderings applied prior to
incomplete factorizations of the form ILU(0), ILU(1), and ILUT and the iterative methods
GMRES(20), BiCGSTAB, and QMR. In some cases, the method will only converge after
the permutation, in others it greatly improves the convergence.
However, we emphasize that permuting large entries to the diagonal of matrix will
not always improve the accuracy and stability of incomplete factorization. An inaccurate
factorization can also occur in the absence of small pivots, when many (especially large)
fill-ins are dropped from the incomplete factors. In this respect, it may be beneficial to
apply a symmetric permutation after the matching reordering to reduce fill-in. Another
kind of instability in incomplete factorizations, which can occur with and without small
pivots, is severe ill-conditioning of the triangular factors. (In this situation, jjRjj F need
not be very large, but jjI \Gamma A( -
will be.) This is also a common situation when
the coefficient matrix is far from diagonally dominant.

Table

7.6: Number of iterations required by some preconditioned iterative methods after
permutation.
Matrix and method Matching algorithm
QMR 72 21 12 12
MAHINDAS
WEST0497
We also performed a set of experiments in which we first permuted the columns of the
matrix A by using a reordering computed by one of the matching algorithms, followed by
a symmetric permutation of A generated by the reverse Cuthill-McKee ordering (Cuthill
and McKee 1969) applied to A . The motivation behind this is that the number
of entries that is dropped from the factors can be reduced by applying a reordering of
the matrix that reduces fill-in. In the experimental results, we noticed that the additional
permutation sometimes has a positive as well as a negative effect on the performance of the
iterative solvers. Table 7.7 shows some results for the three iterative methods from Table
7.6 preconditioned by ILUT on the WEST matrices from the Harwell-Boeing collection.

Table

7.7: Number of iterations required by some ILUT-preconditioned iterative methods
after the matching reordering with and without reverse Cuthill-McKee.
Matrix and method Matching algorithm Matching algorithm
without RCM with RCM
7.2.2 Experiments with a block iterative solution method
The Jacobi method is not a particularly current or powerful method so we focussed our
experiments on 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 the sparse direct method MA27 (HSL 1996).
We show the effect of this in Table 7.8 on the problem MAHINDAS from Table 7.6.
The matching algorithm was followed by a reverse Cuthill-McKee algorithm to obtain a
block tridiagonal form. The matrix was partitioned into 2, 4, 8, and 16 blocks of rows and
the accelerations used were block CG algorithms with block sizes 1, 4, and 8. The block
rows are chosen of equal (or nearly equal) size.

Table

7.8: Number of iterations of block Cimmino algorithm for the matrix MAHINDAS.
Acceleration Matching algorithm
# block rows
None MC21 BT MPD MPS
In general, we noticed in our experiments that the block Cimmino method often was
more sensitive to the scaling (in MPS) and less to the reorderings. The convergence
properties of the block Cimmino method are independent of row scaling. However, the
sparse direct solver MA27 (HSL 1996) used for solving the augmented systems, performs
numerical pivoting during the factorizations of the augmented matrices. Row scaling
might well change the choice of the pivot order and affect the fill-in in the factors and the
accuracy of the solution. Column scaling should affect convergence of the method since it
can be considered as a diagonal preconditioner. For more details see (Ruiz 1992).
8 Conclusions and future work
We have considered, in Sections 3-4, techniques for permuting a sparse matrix so that
the diagonal of the permuted matrix has entries of large absolute value. We discussed
various criteria for this and considered their implementation as computer codes. We also
considered in Section 6 possible scaling strategies to further improve the weight of the
diagonal with respect to the off-diagonal values.
In Section 7, we then indicated several cases where such a permutation (and scaling)
can be useful. These include the solution of sparse equations by a direct method and
by an iterative technique. We also considered its use in generating a preconditioner for
an iterative method. The numerical experiments show that for a multifrontal solver and
preconditioned iterative methods, the effect of these reorderings can be dramatic. The
effect on the block Cimmino iterative method seems to be less dramatic. For this method,
the discussed scaling tends to have a more important effect.
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. One reason for this is that increasing the size of the diagonal only is not always
sufficient to improve the performance of the method. For example, for the incomplete
preconditioners that we used for the numerical experiments in Section 7, it is not only the
size of the diagonal but also the amount and size of the discarded fill-in plays an important
role. We have thus started experimenting with combining the strategies mentioned in
Sections 3-4 and, particularly for generating a preconditioner and the block Cimmino
approach, with combining our unsymmetric ordering with symmetric orderings.
Another interesting extension to the discussed reorderings is a block approach to
increase the size of diagonal blocks instead of only the diagonal entries and use for example
a block Jacobi preconditioner on the permuted matrix. This is of particular interest for
the block Cimmino method. One could also build other criteria into the weighting for
obtaining a bipartite matching, for example, to incorporate a Markowitz cost so that
sparsity would also be preserved by the choice of the resulting diagonal as a pivot. Such
combination would make the resulting ordering suitable for a wider class of sparse direct
solvers.
Finally, we notice in our experiments with MA41 that one effect of the matching
algorithm was to increase the structural symmetry of unsymmetric matrices. We are
exploring further the use of ordering techniques that more directly attempt to increase
structural symmetry.

Acknowledgments

We are grateful to Michele Benzi of Los Alamos National Laboratory and Miroslav Tuma
of the Czech Academy of Sciences for their assistance on the preconditioned iterative
methods and Daniel Ruiz of ENSEEIHT for his help on block iterative methods.



--R



Orderings for incomplete factorization preconditioning of nonsymmetric problems
Assignment and Matching Problems: Solution Methods with FORTRAN-Programs


Experimental study of ILU preconditioners for indefinite matrices
Reducing the bandwidth of sparse symmetric matrices
University of Florida sparse matrix collection
A supernodal approach to sparse partial pivoting



The design and use of algorithms for permuting large entries to the diagonal of sparse matrices


Users' guide for the Harwell-Boeing sparse matrix collection (Release 1)








Making sparse Gaussian elimination scalable by static pivoting

Solution of large sparse unsymmetric linear systems with a block iterative method in a multiprocessor environment
Iterative methods for sparse linear systems
Data structures and network algorithms

--TR

--CTR
Laura Grigori , Xiaoye S. Li, A new scheduling algorithm for parallel sparse LU factorization with static pivoting, Proceedings of the 2002 ACM/IEEE conference on Supercomputing, p.1-18, November 16, 2002, Baltimore, Maryland
Abdou Guermouche , Jean-Yves L'Excellent , Gil Utard, Impact of reordering on the memory of a multifrontal solver, Parallel Computing, v.29 n.9, p.1191-1218, September
Abdou Guermouche , Jean-Yves L'excellent, Constructing memory-minimizing schedules for multifrontal methods, ACM Transactions on Mathematical Software (TOMS), v.32 n.1, p.17-32, March 2006
Chi Shen , Jun Zhang , Kai Wang, Distributed block independent set algorithms and parallel multilevel ILU preconditioners, Journal of Parallel and Distributed Computing, v.65 n.3, p.331-346, March 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
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
Xiren Wang , Wenjian Yu , Zeyi Wang , Xianlong Hong, An improved direct boundary element method for substrate coupling resistance extraction, Proceedings of the 15th ACM Great Lakes symposium on VLSI, April 17-19, 2005, Chicago, Illinois, USA
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
Timothy A. Davis, A column pre-ordering strategy for the unsymmetric-pattern multifrontal method, ACM Transactions on Mathematical Software (TOMS), v.30 n.2, p.165-195, June 2004
Michele Benzi, Preconditioning techniques for large linear systems: a survey, Journal of Computational Physics, v.182 n.2, p.418-477, November 2002
bipartite weighted matching;shortest path algorithms;sparse matrices;iterative methods;preconditioning;direct methods