Implementing Cuthill-McKee

dehann · July 28, 2020, 1:18pm

Hi @PetrKryslUCSD, oh that’s interesting and would be great if RCM gets good performance since it’s already available as native Julia code and well connected to LightGraphs.jl.

Just for reference, I did recently come across another claim that some adaptations of GPS is faster, although its not immediately clear how everyone is comparing (I’ll have to look more in detail to be sure):
https://dl.acm.org/doi/pdf/10.1145/355993.355998?casa_token=5Glhu4oMK6EAAAAA:DskdnaX9oHkDRyeSYqKiTot0WMu8dFMZJ5t0-59QXIqGXw1_YArWEA1G-mMn3vJj0VEKfHCsThWQ

Hi @dpo, sure that is AMD.jl wrapper code. The library libamd.so that is used by AMD.jl is licensed differently, likely LGPL – following:
https://opensuse.pkgs.org/tumbleweed/opensuse-oss-x86_64/libamd2-2.4.6-35.3.x86_64.rpm.html

If the functionality is coming through SuiteSparse (e.g. CCOLAMD) then it’s unclear if that requires a commercial license from SuiteSparse.

dpo · July 28, 2020, 3:31pm

The library libamd.so that is used by AMD.jl is licensed differently, likely LGPL – following:
https://opensuse.pkgs.org/tumbleweed/opensuse-oss-x86_64/libamd2-2.4.6-35.3.x86_64.rpm.html

If the functionality is coming through SuiteSparse (e.g. CCOLAMD) then it’s unclear if that requires a commercial license from SuiteSparse.

BSD 3 clause: https://github.com/DrTimothyAldenDavis/SuiteSparse/blob/master/AMD/Doc/License.txt

PetrKryslUCSD · July 28, 2020, 3:35pm

Do you have GPS implemented at this point?

PetrKryslUCSD · August 2, 2020, 8:33pm

I have at this point an implementation of the reverse Cuthill-McKee
that is 2 times faster than the Matlab implementation of that algorithm, viz

>> S = sprand(100000, 100000, 0.01);
>> S = S + S';
>> tic; r = symrcm(S); toc
Elapsed time is 2.874957 seconds.
>>

(the above edited: I forgot to symmetrize S)
vs

julia> using FinEtools                                                                                                              
julia> using SparseArrays                                                                                                                       
julia> let                                                                                                                                                           
           S = sprand(100000, 100000, 0.01);                                                                                                                         
               @time ag = adjgraph(S);                                                                                                                               
           @time nd = nodedegrees(ag);                                                                                                                               
           @time p = revcm(ag, nd);                                                                                                                                  
       true                                                                                                                                                          
       end;                                                                                                                                                          
  0.536772 seconds (100.00 k allocations: 779.344 MiB, 3.84% gc time)                                                                                                
  0.002363 seconds (2 allocations: 781.328 KiB)                                                                                                                      
  0.444042 seconds (13 allocations: 3.245 MiB)

dpo · August 2, 2020, 10:41pm

Any chance that you might extract that function into its own package? I think it would be generally useful, independently of FEM.

PetrKryslUCSD · August 3, 2020, 12:27am

Has this happened? I guess not, since I cannot find it in the LightGraphs code…

PetrKryslUCSD · August 3, 2020, 8:27pm

That can certainly be done. The question is what the API should be. Is this something that should have a general applicability to LightGraphs?

The Matlab equivalent functionality is available to users through the call p = symrcm(S).
In other words, work directly with a sparse matrix. To accommodate this is easy.

dpo · August 3, 2020, 8:36pm

For my use case, that would be the best.

dehann · August 6, 2020, 7:27am

Hi, sorry bit behind on this thread. A sparse matrix API would probably be the best, since systems are usually large. The method is applicable to both matrix and graph calculations. We currently convert a LightGraphs.jl object into a SparseArrays.SparseMatrixCSC{Int64,Int64} – e.g.

using IncrementalInference
fg = generateCanonicalFG_Kaess()
drawGraph(fg, show=true)
DFG.getBiadjacencyMatrix(fg)

NamedTuple{(:B, :varLabels, :facLabels),Tuple{SparseArrays.SparseMatrixCSC,Array{Symbol,1},Array{Symbol,1}}}((
  [2, 1]  =  1
  [3, 1]  =  1
  [3, 2]  =  1
  [4, 2]  =  1
  [6, 2]  =  1
  [5, 3]  =  1
  [5, 4]  =  1
  [6, 4]  =  1
  [1, 5]  =  1
  [2, 5]  =  1
  [4, 5]  =  1, [:l1, :x2, :l2, :x3, :x1], [:x1f1, :x1l1f1, :x2l1f1, :x1x2f1, :x3l2f1, :x2x3f1]))

and our current use is (not yet RCS, but the default in qr which is likely AMD)

julia> getEliminationOrder(fg)
5-element Array{Symbol,1}:
 :l2
 :l1
 :x3
 :x1
 :x2

Do you have GPS implemented at this point?

No not yet, but we were about to go down this path until @PetrKryslUCSD 's earlier comments (thanks!), and I’d much rather just use the RCM implementation (from FinETools.jl). I can certainly help out, and if it’s okay to ask for collaborator access to help build and maintain such a new repo if possible please? There’s probably a better place to host, but if not I can start a repo at JuliaRobotics (with MIT or BSD license)?

I suspect in the long term more (possibly new) algorithms are likely to be built alongside RCS in a “standardized” repo, RCS would be a great way to get that started and have as reference solution!

dehann · August 6, 2020, 8:29am

Oh, I should probably add motivation for suggesting a separate package for RCM-like functionality. If it were included in LightGraphs.jl, then it’s difficult to use for bipartite graphs (as in our case). It also limits the options for studying/using/modifying such algorithms from either matrix or graphical origin. This is the problem with the ordering methods in SuiteSparse/BLAS/qr/Cholesky, any attempt to develop new ones requires wrappers or surgical compiling etc. Since AMD.jl (libamd.so) also exists as a separate repo, I think it’s probably best to have a separate repo for the pure Julia implementation – also an NP-hard problem in its own right. This may also help attract more community contributions in the future.

PetrKryslUCSD · August 6, 2020, 4:20pm

Jiawei_Duan · February 18, 2023, 4:17am

GPS paper compares a complete RCM(without the pseudoperipheral vertex lookup process), means RCM will generate level structures for all vertices which degree is under median of all vertices’ degrees. Of course it’s blazing fast…