Skip to content

Commit 25a2f14

Browse files
samuelsonricsamuelsonric
authored
Add support for reverse Cuthill-McKee algorithm. (#470)
* Add function `symrcm`. * Add documentation and unit tests.
1 parent c51d3f1 commit 25a2f14

File tree

6 files changed

+165
-1
lines changed

6 files changed

+165
-1
lines changed

‎Project.toml‎

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,14 +9,17 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
99
PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
1010

1111
[weakdeps]
12+
CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8"
1213
SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
1314

1415
[extensions]
1516
BandedMatricesSparseArraysExt = "SparseArrays"
17+
CliqueTreesExt = "CliqueTrees"
1618

1719
[compat]
1820
Aqua = "0.8"
1921
ArrayLayouts = "1.6.0"
22+
CliqueTrees = "0.5"
2023
Documenter = "1"
2124
FillArrays = "1.3"
2225
GenericLinearAlgebra = "0.3"
@@ -31,6 +34,7 @@ julia = "1.10"
3134

3235
[extras]
3336
Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
37+
CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8"
3438
Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
3539
GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
3640
InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
@@ -40,4 +44,4 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
4044
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
4145

4246
[targets]
43-
test = ["Aqua", "Documenter", "GenericLinearAlgebra", "InfiniteArrays", "Random", "SparseArrays", "Test", "Quaternions"]
47+
test = ["Aqua", "CliqueTrees", "Documenter", "GenericLinearAlgebra", "InfiniteArrays", "Random", "SparseArrays", "Test", "Quaternions"]

‎docs/src/index.md‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -162,6 +162,44 @@ To loop over the nonzero elements of a BandedMatrix, you can use `colrange(A, c)
162162

163163
Use `Symmetric(::BandedMatrix)` to work with symmetric banded matrices.
164164

165+
## Bandwidth minimization
166+
167+
The [reverse Cuthill-McKee algorithm](https://en.wikipedia.org/wiki/Cuthill–McKee_algorithm) permutes a symmetric matrix into a band matrix with small bandwidth.
168+
This algorithm is implemented by the function `symrcm`.
169+
170+
```julia
171+
julia> using BandedMatrices, CliqueTrees
172+
173+
julia> matrix = [
174+
1 0 1 0 0 0 0 1 1
175+
0 1 1 0 0 1 1 1 0
176+
1 1 1 0 0 0 0 0 0
177+
0 0 0 1 0 0 0 1 1
178+
0 0 0 0 1 0 1 1 0
179+
0 1 0 0 0 1 1 0 0
180+
0 1 0 0 1 1 1 0 0
181+
1 1 0 1 1 0 0 1 0
182+
1 0 0 1 0 0 0 0 1
183+
];
184+
185+
julia> bandedmatrix, perm, invp = symrcm(matrix);
186+
187+
julia> bandedmatrix
188+
9×9 LinearAlgebra.Symmetric{Int64, BandedMatrix{Int64, Matrix{Int64}, Base.OneTo{Int64}}}:
189+
1 1 1 0 ⋅ ⋅ ⋅ ⋅ ⋅
190+
1 1 1 1 0 ⋅ ⋅ ⋅ ⋅
191+
1 1 1 0 1 1 ⋅ ⋅ ⋅
192+
0 1 0 1 0 1 0 ⋅ ⋅
193+
⋅ 0 1 0 1 0 1 0 ⋅
194+
⋅ ⋅ 1 1 0 1 1 1 0
195+
⋅ ⋅ ⋅ 0 1 1 1 0 1
196+
⋅ ⋅ ⋅ ⋅ 0 1 0 1 1
197+
⋅ ⋅ ⋅ ⋅ ⋅ 0 1 1 1
198+
199+
julia> bandedmatrix == matrix[perm, perm]
200+
true
201+
```
202+
165203
## Banded matrix interface
166204

167205
Banded matrices go beyond the type `BandedMatrix`: one can also create

‎ext/CliqueTreesExt.jl‎

Lines changed: 35 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,35 @@
1+
module CliqueTreesExt
2+
3+
using BandedMatrices
4+
using Base.Sort: Algorithm
5+
using CliqueTrees: CliqueTrees, permutation, RCMGL
6+
using CliqueTrees.SparseArrays
7+
using LinearAlgebra
8+
9+
function BandedMatrices.symrcm(matrix::AbstractMatrix, alg::Algorithm)
10+
BandedMatrices.symrcm(sparse(matrix), alg)
11+
end
12+
13+
function BandedMatrices.symrcm(matrix::SparseMatrixCSC{T}, alg::Algorithm) where T
14+
order, index = permutation(matrix; alg=RCMGL(alg))
15+
lower = tril!(permute(matrix, order, order))
16+
17+
bandwidth = maximum(axes(lower, 2)) do j
18+
p = last(nzrange(lower, j))
19+
return rowvals(lower)[p] - j
20+
end
21+
22+
banded = BandedMatrix{T}(undef, size(lower), (bandwidth, 0))
23+
fill!(banded.data, zero(T))
24+
25+
@inbounds for j in axes(lower, 2)
26+
for p in nzrange(lower, j)
27+
i = rowvals(lower)[p]
28+
banded.data[i - j + 1, j] = nonzeros(lower)[p]
29+
end
30+
end
31+
32+
return Symmetric(banded, :L), order, index
33+
end
34+
35+
end

‎src/BandedMatrices.jl‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,7 @@ export BandedMatrix,
5252
bandwidth,
5353
BandError,
5454
band,
55+
symrcm,
5556
Band,
5657
BandRange,
5758
bandwidths,
@@ -87,6 +88,7 @@ include("symbanded/BandedCholesky.jl")
8788
include("symbanded/SplitCholesky.jl")
8889
include("symbanded/bandedeigen.jl")
8990
include("symbanded/tridiagonalize.jl")
91+
include("symbanded/symrcm.jl")
9092

9193
include("tribanded.jl")
9294

‎src/symbanded/symrcm.jl‎

Lines changed: 42 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,42 @@
1+
"""
2+
symrcm(matrix[, alg::Algorithm])
3+
4+
Use the Reverse Cuthill-Mckee algorithm to reduce the bandwith of a matrix.
5+
6+
# Examples
7+
```julia
8+
julia> using BandedMatrices, CliqueTrees
9+
10+
julia> matrix = [
11+
1 0 1 0 0 0 0 1 1
12+
0 1 1 0 0 1 1 1 0
13+
1 1 1 0 0 0 0 0 0
14+
0 0 0 1 0 0 0 1 1
15+
0 0 0 0 1 0 1 1 0
16+
0 1 0 0 0 1 1 0 0
17+
0 1 0 0 1 1 1 0 0
18+
1 1 0 1 1 0 0 1 0
19+
1 0 0 1 0 0 0 0 1
20+
];
21+
22+
julia> bandedmatrix, perm, invp = symrcm(matrix);
23+
24+
julia> bandedmatrix
25+
9×9 LinearAlgebra.Symmetric{Int64, BandedMatrix{Int64, Matrix{Int64}, Base.OneTo{Int64}}}:
26+
1 1 1 0 ⋅ ⋅ ⋅ ⋅ ⋅
27+
1 1 1 1 0 ⋅ ⋅ ⋅ ⋅
28+
1 1 1 0 1 1 ⋅ ⋅ ⋅
29+
0 1 0 1 0 1 0 ⋅ ⋅
30+
⋅ 0 1 0 1 0 1 0 ⋅
31+
⋅ ⋅ 1 1 0 1 1 1 0
32+
⋅ ⋅ ⋅ 0 1 1 1 0 1
33+
⋅ ⋅ ⋅ ⋅ 0 1 0 1 1
34+
⋅ ⋅ ⋅ ⋅ ⋅ 0 1 1 1
35+
36+
julia> bandedmatrix == matrix[perm, perm]
37+
true
38+
```
39+
"""
40+
function symrcm(matrix::AbstractMatrix)
41+
symrcm(matrix, Base.Sort.DEFAULT_UNSTABLE)
42+
end

‎test/test_symbanded.jl‎

Lines changed: 43 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@ module TestSymBanded
33
using ArrayLayouts
44
using BandedMatrices
55
import BandedMatrices: MemoryLayout, SymmetricLayout, HermitianLayout, BandedColumns, isbanded
6+
using CliqueTrees: CliqueTrees
67
using FillArrays
78
using GenericLinearAlgebra
89
using LinearAlgebra
@@ -421,4 +422,46 @@ end
421422
end
422423
end
423424

425+
@testset "symrcm" begin
426+
matrix = Int16[
427+
1 0 1 0 0 0 0 1 1
428+
0 1 1 0 0 1 1 1 0
429+
1 1 1 0 0 0 0 0 0
430+
0 0 0 1 0 0 0 1 1
431+
0 0 0 0 1 0 1 1 0
432+
0 1 0 0 0 1 1 0 0
433+
0 1 0 0 1 1 1 0 0
434+
1 1 0 1 1 0 0 1 0
435+
1 0 0 1 0 0 0 0 1
436+
]
437+
438+
bandedmatrix, perm, invp = symrcm(matrix)
439+
@test invp == invperm(perm)
440+
@test bandedmatrix == matrix[perm, perm]
441+
@test eltype(bandedmatrix) == eltype(matrix)
442+
443+
# LFAT5
444+
matrix = Float32[
445+
0.608806 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -0.304403 0.0
446+
0.0 1.57088 0.0 -94.2528 0.78544 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
447+
0.0 0.0 15080.4 0.0 0.0 -94.2528 0.0 0.0 94.2528 0.0 -7540.22 0.0 0.0 0.0
448+
0.0 -94.2528 0.0 15080.4 0.0 94.2528 0.0 0.0 0.0 0.0 -7540.22 0.0 0.0 0.0
449+
0.0 0.78544 0.0 0.0 3.14176 0.78544 0.0 0.0 0.0 0.0 -94.2528 0.0 0.0 0.0
450+
0.0 0.0 -94.2528 94.2528 0.78544 3.14176 0.0 0.0 0.0 0.78544 0.0 0.0 0.0 0.0
451+
0.0 0.0 0.0 0.0 0.0 0.0 1.25664f7 0.0 0.0 0.0 0.0 -6.2832f6 0.0 -6.2832f6
452+
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.608806 0.0 0.0 0.0 0.0 -0.304403 0.0
453+
0.0 0.0 94.2528 0.0 0.0 0.0 0.0 0.0 1.57088 0.78544 0.0 0.0 0.0 0.0
454+
0.0 0.0 0.0 0.0 0.0 0.78544 0.0 0.0 0.78544 3.14176 94.2528 0.0 0.0 0.0
455+
0.0 0.0 -7540.22 -7540.22 -94.2528 0.0 0.0 0.0 0.0 94.2528 15080.4 0.0 0.0 0.0
456+
0.0 0.0 0.0 0.0 0.0 0.0 -6.2832f6 0.0 0.0 0.0 0.0 1.25664f7 0.0 0.0
457+
-0.304403 0.0 0.0 0.0 0.0 0.0 0.0 -0.304403 0.0 0.0 0.0 0.0 0.608806 0.0
458+
0.0 0.0 0.0 0.0 0.0 0.0 -6.2832f6 0.0 0.0 0.0 0.0 0.0 0.0 1.25664f7
459+
]
460+
461+
bandedmatrix, perm, invp = symrcm(matrix)
462+
@test invp == invperm(perm)
463+
@test bandedmatrix == matrix[perm, perm]
464+
@test eltype(bandedmatrix) == eltype(matrix)
465+
end
466+
424467
end # module

0 commit comments

Comments
 (0)