Implementing Cuthill-McKee

#1

I am trying to implement Cuthill McKee on a number of large, sparse matrices. In particular, I would like to get the mapping from original row order to post-CM (or post-reverseCM) row order.

The only package that I can find that implements Cuthill McKee does this using nodes. In particular, the example linked above wants a Dict for the sparsity pattern (and another dict representing the node-types).

First, is there a cleaner way that anyone knows to do this? And, if not, is there a clean way to output the sparsity pattern as a dict?

As of now, I suspect I will just be writing a lot of for-loops, but I am hoping that there is a more elegant solution. Thank you!

0 Likes

#2

@dlfivefifty have you done this in ApproxFun yet?

0 Likes

#3

I implemented it a few years ago. Functio rcm at https://github.com/pjabardo/Makhno.jl/blob/master/src/numbering.jl

I remember it worked. I don’t know if it is efficient. Of course it was taylored for a specific data structure. I will try to convert it to use CSC sparse matrix this week end.

0 Likes

#4

No I haven’t. Would be nice to have in BandedMatrices.jl

0 Likes

#5

Thanks everyone. I found an acceptable intermediate solution for my particular problem that avoids the need for CM/RCM.

@Paulo_Jabardo, if you do ever convert your code to work with CSC Sparse Matrices, please do let me know.

@dlfivefifty, if I end up writing my own reasonably efficient code for RCM, I’ll try to contribute it to BandedMatrices.jl.

3 Likes

#6

If someone is still interested in this, I had it lying around so I made a quick package. Probably quite rough, so anyone is welcome to make it prettier. See: https://github.com/rleegates/CuthillMcKee.jl

3 Likes

#7

Thanks for sharing Robert.

I do have some nagging doubts about the source algorithm (Zavoianu).
I looked at the Zavoianu version some time ago, and it differed in some ways
from other variants of the RCM I’ve seen in the literature.

I implemented a variant of the version from Carey’s Computational Grids book in FinEtools.

image

One of my modifications enables renumbering of grids with disconnected components.

Disclaimer: I haven’t optimized the speed for larger matrices yet.
For smaller matrices it runs okay.

0 Likes

#8

Oh, I wasn’t aware of this, thanks for the heads up! Do you know what the suspected error is? The output sparsity looks okay to me. This RCM algorithm was part of my FE toolbox for quite some time and I never had a problem, but perhaps there is a minor issue that I may have missed? I do suspect the matrices were SPD and the graph connected in that part of the code. I don’t have a reason to believe unconnected components present a problem though, I do detect and count them, see the symrcm function docs. The permutation “looks fine” in either case. I updated the README.md with some verification and some images.

1 Like

#9

Zavoianu’s algorithm calls for a constant maintenance of the queue Q. The child C should add its neighbors immediately, in a sorted order, to Q. The queue would then grow very long, and would also need to be sorted all the time. I don’t do any of that.

0 Likes

#10

What a coincidence, LightGraphs is getting that algorithm too; https://github.com/JuliaGraphs/LightGraphs.jl/pull/1128

3 Likes

#11

That is pretty cool. That would have saved me some time.
By the way, where is it, I can’t find it. Did you rename it?

0 Likes

#12

Should be in that pull request, that I linked (it is not merged yet)

1 Like

#13

Perhaps it would be good to agree on a representation of adjacency and degree and then make a collection of bandwidth reduction algorithms in some central place? It could be called BandwidthReduction.jl or something similar? Then other packages, e.g. FinEtools.jl, LightGraphs.jl, BandedMatrices.jl could simply call upon these algorithms.

I must say that I don’t quite agree with the assertion made in the LightGraphs pull-request discussion that graph permutation algorithms are essentially a preprocessing task and that their performance is therefore a secondary consideration. This may be the case for a graph which remains static throughout a computation. My use case, for example, is the following: I often work with hierarchical grids, e.g. hierarchical simplicial finite elements, sparse grids with wavelet bases, etc. which are continuously refined throughout the computation based on some form of error measure. Depending on the method of integration over these grids, it is often possible to compute integrals (the weights attached to edges in the graph) involving only the additional basis functions (the new vertices in the graph) while leaving the “old” components of the operator unchanged. For this purpose I use some form of a dynamic graph structure for insertion and deletion. In a perfect world, in addition to having a performant renumbering algorithm, I desire some form of “online” updating algorithm for the numbering of the basis functions in the sparse operator representation instead of just appending them to the operator (in the lower right block) and then permuting.

5 Likes