Bringing Nonzero Elements to the Diagonal

The pymc21 Module

A Python interface to the HSL subroutine MC21AD.

References

[Duff81a]

I. S. Duff, On Algorithms for Obtaining a Maximum Transversal, ACM Trans. Math. Software, 7, pp. 315-330, 1981.

[Duff81b]

I. S. Duff, Algorithm 575: Permutations for a zero-free diagonal, ACM Trans. Math. Software, 7, pp. 387-390, 1981.

Functions available

pyorder.pymc21.pymc21.nonzerodiag(nrow, colind, rowptr)

Given the sparsity pattern of a square sparse matrix in compressed row (csr) format, attempt to find a row permutation so the row-permuted matrix has a nonzero diagonal, if this is possible. This function assumes that the matrix indexing is 0-based. Note that in general, this function does not preserve symmetry. The method used is a depth-first search with lookahead described in [Duff81a] and [Duff81b].

Parameters:

nrow:The number of rows of the input matrix. colind:An integer array (or list) of length nnz giving the column indices of the nonzero elements in each row. rowptr:An integer array (or list) of length nrow+1 giving the indices of the first element of each row in colind.

Returns:

perm:An integer array of length nrow giving the variable permutation. If irow and jcol are two integer arrays describing the pattern of the input matrix in triple format, perm[irow] and jcol describe the permuted matrix. nzdiag:The number of nonzeros on the diagonal of the permuted matrix.

Examples

Basic Usage

This first example calls the Fortran subroutine directly with the matrix data in compressed sparse column format.

Warning

Keep in mind that when calling the Fortran subroutines directly, all indices in the matrix data must be 1-based, i.e., row indices range from 1 through nrow and column indices range from 1 through ncol.

import numpy as np
from pyorder.pymc21 import mc21module
n = 4
icn = np.array([1,4,3,4,1,4,2,4], dtype=np.int32)
ip = np.array([1,3,5,7], dtype=np.int32)
lenr = np.array([2,2,2,2], dtype=np.int32)
iperm,numnz = mc21module.mc21ad(icn,ip,lenr)
print 'iperm = ', iperm, ' (1-based)'
print 'numnz = ', numnz

Python Interface

In this second example, Python usage is intended. Matrix indices must therefore be 0-based. The permutation vector returned by nonzerodiag() is now also 0-based.

"""
Illustrate usage of the pymc21 module, using an input matrix in Harwell-Boeing
or Rutherford-Boeing format. Supply a file name as input argument on the
command line and uncomment below as appropriate.
"""

import sys
import numpy as np
from pyorder.pymc21.pymc21 import nonzerodiag
from pyorder.tools.hrb import HarwellBoeingMatrix, RutherfordBoeingData

if len(sys.argv) < 2:
    sys.stderr.write('Supply input matrix as argument\n')
    sys.exit(1)

fname = sys.argv[1]
#M = HarwellBoeingMatrix(fname, patternOnly=True, readRhs=False)
M = RutherfordBoeingData(fname, patternOnly=True, readRhs=False)

if M.nrow != M.ncol:
    sys.stderr.write('Input matrix must be square\n')
    sys.exit(1)

perm, nzdiag = nonzerodiag(M.nrow, M.ind, M.ip)
print 'Number of nonzeros on diagonal after reordering: ', nzdiag

try:
    from pyorder.tools.spy import FastSpy
    import pylab
    (irow, jcol) = M.find()
    left = pylab.subplot(121)
    FastSpy(M.nrow, M.ncol, irow, jcol, sym=M.issym, ax=left.get_axes())

    right = pylab.subplot(122)
    FastSpy(M.nrow, M.ncol, perm[irow], jcol, sym=M.issym, ax=right.get_axes())
    pylab.show()
except:
    pass

Reorganizing a square matrix so it has a nonzero diagonal is useful for instance when factorizing a submatrix of a basis in the Simplex method for linear programming.

The effect of nonzerodiag() is illustrated on the following example, shl_0 from the University of Florida Sparse Matrix Collection. I used the matrix in Rutherford-Boeing format.