Block-triangular-form permutation?

question

#1

Is there a function that permutes a sparse matrix (or its incidence graph) into block triangular form (BTF)?


#2

If you have a DLL of libbtf from SuiteSparse, this should work:

import Base.SparseArrays.CHOLMOD: SuiteSparse_long

function btf{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti})
    @assert A.m == A.n
    @assert Ti == SuiteSparse_long
    P = zeros(Ti,A.n)
    Q = zeros(Ti,A.n)
    R = zeros(Ti,A.n+1)
    Work = zeros(Ti,5*A.n)
    work = 1.0;
    maxwork = -1.0;
    # convert to C indices
    colptr = A.colptr-1
    rowval = A.rowval-1
    nmatch = one(Ti)

    ccall((:btf_l_order,"libbtf"),
           SuiteSparse_long, (SuiteSparse_long, Ptr{SuiteSparse_long},
                              Ptr{SuiteSparse_long}, Float64, Ptr{Float64},
                              Ptr{SuiteSparse_long}, Ptr{SuiteSparse_long},
                              Ptr{SuiteSparse_long}, Ptr{SuiteSparse_long},
                              Ptr{SuiteSparse_long}),
           A.n, colptr, rowval, maxwork, &work, P, Q, R, &nmatch, Work)
    # convert to Julia indices
    for i in eachindex(P)
        P[i] += 1
        Q[i] += 1
    end
    for i in eachindex(R)
        R[i] += 1
    end
    P,Q,R
end

p,q,r = btf(A)
A_btf = A[p,abs(q)]

#3

@Ralph_Smith Thanks, works indeed. Discovered your reply just now. Meanwhile I wrote an interface to Trilinos BTF (Fortran) version (Pothen & Fan, 1990) (source: https://github.com/trilinos/Trilinos/tree/master/packages/epetraext/src/btf/pothen) that also handles non-square matrices. Please find the Julia code below. Thanks again!

"""
  pprint(::AbstractArray)

Array pretty printing.
"""
function pprint(A::AbstractArray, pre = "")
  print(pre, typeof(A), size(A), "[")
  for i in 1:size(A)[1]
    println("");print("  ")
    for a in A[i,:]
      print(a,"\t")
    end
  end
  println("]")
end

# BTF

libdir = "..."

const libBTF = libdir * "libBTF"


"""
  BTF

Block Triangular Form data:

  form ∈ ["L","U"] -- triangular form, L(ower) or U(pper)

  rowperm   -- the new to old permutation array for the rows
  colperm   -- the old to new permutation array for the columns

  nhrows, nhcols, hrzcmp
         -- number of rows, columns and connected components
            in the horizontal (underdetermined) block

  nsrows, sqcmpn
         -- number of rows (and columns) and strong components
            in the square (exactly determined) block

  nvrows, nvcols, vrtcmp
         -- number of rows, columns and connected components
            in the vertical (overdetermined) block

  rcmstr -- index of first row in a diagonal block
            (component starting row)
            where
                (rcmstr(1), ..., rcmstr(hrzcmp)) give the
                     indices for the components in the
                      horizontal block
                (rcmstr(hrzcmp+1), ..., rcmstr(hrzcmp+sqcmpn))
                     give the indices for the components in the
                     square block
                (rcmstr(hrzcmp+sqcmpn+1), ...,
                 rcmstr(hrzcmp+sqcmpn+vrtcmp)) give the indices
                     for the components in the vertical block
                 rcmstr(hrzcmp+sqcmpn+vrtcmp+1) is equal to
                     nrows+1 for convenience

  ccmstr -- index of first column in a diagonal block
            (component starting column)
            where
                (ccmstr(1), ..., ccmstr(hrzcmp)) give the
                     indices for the components in the
                      horizontal block
                (ccmstr(hrzcmp+1), ..., ccmstr(hrzcmp+sqcmpn))
                     give the indices for the components in the
                     square block, making this block itself
                     in block lower triangular form
                (ccmstr(hrzcmp+sqcmpn+1), ...,
                 ccmstr(hrzcmp+sqcmpn+vrtcmp)) give the indices
                     for the components in the vertical block
                 ccmstr(hrzcmp+sqcmpn+vrtcmp+1) is equal to
                     ncols+1 for convenience
"""
type BTF

  form::String

  rowperm::Array{Int32}
  colperm::Array{Int32}

  nhrows::Int32
  nhcols::Int32
  hrzcmp::Int32

  nsrows::Int32
  sqcmpn::Int32

  nvrows::Int32
  nvcols::Int32
  vrtcmp::Int32

  rcmstr::Array{Int32}
  ccmstr::Array{Int32}

  function BTF(form,
               rowperm, colperm,
               nhrows, nhcols, hrzcmp,
               nsrows, sqcmpn,
               nvrows, nvcols, vrtcmp,
               rcmstr, ccmstr)

    assert(form in ["L","U"])
    _ = new()
    _.form = form
    _.rowperm, _.colperm = rowperm, colperm
    _.nhrows, _.nhcols, _.hrzcmp = nhrows, nhcols, hrzcmp
    _.nsrows, _.sqcmpn = nsrows, sqcmpn
    _.nvrows, _.nvcols, _.vrtcmp = nvrows, nvcols, vrtcmp
    _.rcmstr, _.ccmstr = rcmstr, ccmstr
    return _
  end

end

"""
    permBTF{Tv}(A::SparseMatrixCSC{Tv,Int32}
                ;
                form = "L", msglvl = 0, output = 6)

Finds the lower block triangular form of the square submatrix
in the general block triangular form. The square submatrix
consists entirely of matched rows and columns. Therefore,
with each row matched to its matching column, the submatrix
has a nonzero diagonal, as required by duff's algorithm.

From a graph-theoretic standard, this is the same as considering
the subgraph induced by sr and sc, if non-matching edges
are directed from rows to columns, and matching edges are shrunk
into single vertices, the resulting directed graph has strongly
connected components. This is the Dulmage-Mendelsohn decomposition.

Uses mmc13e (modified version from harwell mc13e by alex pothen and
chin-ju fan; bcs modifications, john lewis, sept. 1990)

mmc13e uses Tarjan's algorithm to find the strongly connected
components by depth-first search. All the pairs have been visited
will be labeled in the order they are visited, and associated a
lowlink for each vertex, stored in the stack - lowlnk.
"""
function permBTF{Tv}(A::SparseMatrixCSC{Tv,Int32}
                     ;
                     form = "L", msglvl = 0, output = 6)

    msglvl::Int32 = msglvl
    output::Int32 = output

    At = A'

    # Arrange for Lower or Upper form
    #
    if     form == "L"
      rowstr,colidx,colstr,rowidx = A.colptr,A.rowval,At.colptr,At.rowval
    elseif form == "U"
      colstr,rowidx,rowstr,colidx = A.colptr,A.rowval,At.colptr,At.rowval
    else
      assert(form in ["L","U"])
    end

    #

    nrows::Int32, ncols::Int32 = size(A)

    w = zeros(Int32,10*max(nrows,ncols))

    rnto, cnto = zeros(Int32,nrows), zeros(Int32,ncols)

    # Receiving values as Fortran scalar arguments must go through scalar
    # pointers (i.e. memory must be reserved in Julia).

    # horizontal block:  rows, columns, connected components
    nhrows = Int32[0]  # each rhs is a memloc;
    nhcols = Int32[0]  # ...don't do a = b = Int[0]: share same memloc.
    hrzcmp = Int32[0]
    # square block:  rows=columns, connected components
    nsrows = Int32[0]
    sqcmpn = Int32[0]
    # vertical block:  rows, columns, connected components
    nvrows = Int32[0]
    nvcols = Int32[0]
    vrtcmp = Int32[0]

    rcmstr, ccmstr = zeros(Int32,nrows+1), zeros(Int32,ncols+1)

    ccall((:genbtf_, libBTF),
            Void,
            (
            Ref{Int32}, Ref{Int32},
            Ref{Int32}, Ref{Int32}, Ref{Int32}, Ref{Int32},
            Ptr{Int32},
            Ptr{Int32}, Ptr{Int32},
            Ptr{Int32}, Ptr{Int32}, Ptr{Int32},
            Ptr{Int32}, Ptr{Int32},
            Ptr{Int32}, Ptr{Int32}, Ptr{Int32},
            Ptr{Int32}, Ptr{Int32},
            Ptr{Int32}, Ptr{Int32}
            ),
            nrows, ncols,
            #colstr,   rowidx,    rowstr,   colidx,
            #At.colptr, At.rowval, A.colptr, A.rowval,
            rowstr,   colidx,   colstr,    rowidx,
            #A.colptr, A.rowval, At.colptr, At.rowval,
            w     ,
            rnto  , cnto,
            nhrows, nhcols, hrzcmp,
            nsrows, sqcmpn,
            nvrows, nvcols, vrtcmp,
            rcmstr, ccmstr,
            &msglvl, &output)

    p,q = (form=="L")? (rnto, cnto) : (cnto, rnto)

    btf = BTF(form,
              p, q,
              nhrows[1], nhcols[1], hrzcmp[1],
              nsrows[1], sqcmpn[1],
              nvrows[1], nvcols[1], vrtcmp[1],
              rcmstr, ccmstr)

    return btf
end

"""
  permBTF{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}
                 ;
                 form = "L", msglvl = 0, output = 6)

Finds the lower block triangular form of the square submatrix
in the general block triangular form. Also see:

  permBTF{Tv}(A::SparseMatrixCSC{Tv,Int32}
              ;
              form="L", msglvl = 0, output = 6)
"""
function permBTF{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}
                        ;
                        form = "L", msglvl = 0, output = 6)
    return permBTF(convert(SparseMatrixCSC{Tv,Int32}, A),
                   form = form, msglvl = msglvl, output = output)
end

# test

using CxSupport.ArrayUtils: pprint

# Values +1 for illustration.
A = sparse([ 1, 1, 1, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 7, 8 ],
           [ 1, 5, 6, 2, 5, 6, 3, 5, 6, 7, 4, 6, 8, 1, 3, 4, 5, 2, 4, 5, 1, 2 ],
           [ 3, 1, 1, 3, 1, 1, 3, 1, 1, 1, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 ]
           )

pprint(A,"A:")

btf = permBTF(A)
println(btf)
B = A[btf.rowperm,btf.colperm]
pprint(B,"B:")

btf = permBTF(A,form="U")
println(btf)
B = A[btf.rowperm,btf.colperm]
pprint(B,"B:")