BandedMatrices.jl icon indicating copy to clipboard operation
BandedMatrices.jl copied to clipboard

Add support (reverse) Cuthill-McKee algorithm

Open Paalon opened this issue 6 years ago • 3 comments

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.

dlfivefifty avatar Oct 17 '18 09:10 dlfivefifty

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

Paalon avatar Oct 22 '18 19:10 Paalon

I had it lying around so I published a small package to this effect: CuthillMcKee.jl. Perhaps you may find it useful.

rleegates avatar Apr 13 '19 22:04 rleegates