BlockTriangularForm

maxtrans(nrow, ncol, Ap::Vector{Ti}, Ai::Vector{Ti}, maxwork) where {Ti}

Finds a permutation of the columns of a matrix so that it has a zero-free diagonal. The input is an m-by-n sparse matrix in compressed column form. The array Ap of size n+1 gives the starting and ending positions of the columns in the array Ai. Ap[1] must be one. The array Ai contains the row indices of the nonzeros of the matrix A, and is of size Ap[n + 1]. The row indices of column j are located in Ai[Ap[j] ... Ap[j+1]-1]. Row indices must be in the range 0 to m-1. Duplicate entries may be present in any given column. The input matrix is not checked for validity (row indices out of the range 0 to m-1 will lead to an undeterminate result - possibly a core dump, for example). Row indices in any given column need not be in sorted order. However, if they are sorted and the matrix already has a zero-free diagonal, then the identity permutation is returned.

The output of maxtrans is an array Match of size n. If row i is matched with column j, then A(i,j) is nonzero, and then Match[i] = j. If the matrix is structurally nonsingular, all entries in the Match array are unique, and Match can be viewed as a column permutation if A is square. That is, column k of the original matrix becomes column Match[k] of the permuted matrix. This can be expressed as (for non-structurally singular matrices):

nmatch, Match = maxtrans(A)
B = A[:, Match]

If row i is not matched to any column, then Match[i] is == -1. The maxtrans routine returns the number of nonzeros on diagonal of the permuted matrix nmatch and Match.

Notes

This algorithm is based on the paper "On Algorithms for obtaining a maximum transversal" by Iain Duff, ACM Trans. Mathematical Software, vol 7, no. 1, pp. 315-330, and "Algorithm 575: Permutations for a zero-free diagonal", same issue, pp. 387-390. Algorithm 575 is MC21A in the Harwell Subroutine Library. This code is not merely a translation of the Fortran code into C. It is a completely new implementation of the basic underlying method (depth first search over a subgraph with nodes corresponding to columns matched so far, and cheap matching). This code was written with minimal observation of the MC21A/B code itself. See comments below for a comparison between the maxtrans and MC21A/B codes.

This routine operates on a column-form matrix and produces a column permutation. MC21A uses a row-form matrix and produces a row permutation. The difference is merely one of convention in the comments and interpretation of the inputs and outputs. If you want a row permutation, simply pass a compressed-row sparse matrix to this routine and you will get a row permutation (just like MC21A). Similarly, you can pass a column-oriented matrix to MC21A and it will happily return a column permutation.

order(n, Ap::Vector{Ti}, Ai::Vector{Ti}, maxwork = -1) where {Ti}

Permutes a square matrix into upper block triangular form. It does this by first finding a maximum matching (or perhaps a limited matching if the work is limited), via the maxtrans function. If a complete matching is not found, order completes the permutation, but flags the columns of PAQ to denote which columns are not matched. If the matrix is structurally rank deficient, some of the entries on the diagonal of the permuted matrix will be zero. order then calls strongcomp to find the strongly-connected components.

On output, P and Q are the row and column permutations, where i = P[k] if row i of A is the kth row of P*A*Q, and j = unflip(Q[k]) if column j of A is the kth column of P*A*Q. If Q[k] < 1, then the (k,k)th entry in P*A*Q is structurally zero.

The vector R gives the block boundaries, where block b is in rows/columns R[b]:R[b+1]-1 of the permuted matrix, and where b ranges from 1 to the number of strongly connected components found.

strongcomp!(n, Ap::Vector{Ti}, Ai::Vector{Ti}, Q = nothing) where {Ti}

Finds the strongly connected components of a graph, returning a symmetric permutation. The matrix A must be square, and is provided on input in compressed-column form (see maxtrans). The diagonal of the input matrix A (or A*Q if Q is provided on input) is ignored.

If Q is not nothing on input, then the strongly connected components of A*Q are found. Q may be flagged on input, where Q[k] < 0 denotes a flagged column k. The permutation is j = unflip(Q[k]). On output, Q is modified (the flags are preserved) so that P*A*Q is in block upper triangular form.

If Q is nothing, then the permutation P is returned so that P*A*P' is in upper block triangular form.

The vector R gives the block boundaries, where block b is in rows/columns R[b]:R[b+1]-1 of the permuted matrix, and where b ranges from 1 to the number of strongly connected components found.