Is there a function that permutes a sparse matrix (or its incidence graph) into block triangular form (BTF)?
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)]
1 Like
@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:")
1 Like