BandedMatrices.jl
BandedMatrices.jl copied to clipboard
Add support (reverse) Cuthill-McKee algorithm
Discussed in discourse. Wikipedia: Cuthill-McKee algorithm
This would be great to have, but I don't have the time to do this. Would you be interested in putting together a PR?
I'm curious whether the 5 point Finite difference stencil in https://github.com/JuliaMatrices/BandedMatrices.jl/blob/master/examples/finitedifference_2d.jl has better bandwidth after applying that algorithm.
My dirty implementation is like this:
# References
# https://en.wikipedia.org/wiki/Cuthill%E2%80%93McKee_algorithm
# http://ciprian-zavoianu.blogspot.com/2009/01/project-bandwidth-reduction.html
using DataStructures
function cuthillMcKee(x::BitArray{2})
(n, m) = size(x)
n == m || throw(DomainError("does not square matrix"))
for i = 1:n, j = 1:n
if x[i, j] != x[j, i]
throw(DomainError("does not symmetric matrix"))
end
if i == j && x[i, j] != true
throw(DomainError("diagonal elements are not true"))
end
end
permutation = Int[]
x_order = sum(x, dims=2)[:, 1]
while length(permutation) < n
rest = setdiff(1:n, permutation)
queue = Deque{Int}()
push!(queue, rest[argmin(x_order[rest])])
while !isempty(queue)
i = popfirst!(queue)
if i in permutation
continue
end
push!(permutation, i)
neighbor = setdiff(findall(x[i, :]), permutation)
neighbor_order = x_order[neighbor]
p = sortperm(neighbor_order)
neighbor_added = neighbor[p]
if isempty(neighbor_added)
continue
end
push!(queue, neighbor[p]...)
end
end
permutation
end
function reverseCuthillMcKee(x::BitArray{2})
reverse(cuthillMcKee(x))
end
x = [
1 0 0 0 1 0 0 0
0 1 1 0 0 1 0 1
0 1 1 0 1 0 0 0
0 0 0 1 0 0 1 0
1 0 1 0 1 0 0 0
0 1 0 0 0 1 0 1
0 0 0 1 0 0 1 0
0 1 0 0 0 1 0 1
]
p = reverseCuthillMcKee(Bool.(x))
x_ordered = x[p, p]
# 8×8 Array{Int64,2}:
# 1 1 0 0 0 0 0 0
# 1 1 0 0 0 0 0 0
# 0 0 1 1 1 0 0 0
# 0 0 1 1 1 0 0 0
# 0 0 1 1 1 1 0 0
# 0 0 0 0 1 1 1 0
# 0 0 0 0 0 1 1 1
# 0 0 0 0 0 0 1 1
JuliaFEM has another implementation. https://github.com/JuliaFEM/NodeNumbering.jl/blob/master/src/RCM.jl
I had it lying around so I published a small package to this effect: CuthillMcKee.jl. Perhaps you may find it useful.