From 1fac8024c568a7e9ee6651e73c1d77fe77bb5ab1 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Tue, 4 Mar 2025 12:13:34 -0500 Subject: [PATCH 1/2] Add function `symrcm`. --- Project.toml | 4 ++++ ext/CliqueTreesExt.jl | 35 +++++++++++++++++++++++++++++++++++ src/BandedMatrices.jl | 1 + src/symbanded/symrcm.jl | 39 +++++++++++++++++++++++++++++++++++++++ 4 files changed, 79 insertions(+) create mode 100644 ext/CliqueTreesExt.jl create mode 100644 src/symbanded/symrcm.jl diff --git a/Project.toml b/Project.toml index 571b0076..0b9fba90 100644 --- a/Project.toml +++ b/Project.toml @@ -9,14 +9,17 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a" [weakdeps] +CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" [extensions] BandedMatricesSparseArraysExt = "SparseArrays" +CliqueTreesExt = "CliqueTrees" [compat] Aqua = "0.8" ArrayLayouts = "1.6.0" +CliqueTrees = "0.5" Documenter = "1" FillArrays = "1.3" GenericLinearAlgebra = "0.3" @@ -31,6 +34,7 @@ julia = "1.10" [extras] Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a" InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c" diff --git a/ext/CliqueTreesExt.jl b/ext/CliqueTreesExt.jl new file mode 100644 index 00000000..8907e260 --- /dev/null +++ b/ext/CliqueTreesExt.jl @@ -0,0 +1,35 @@ +module CliqueTreesExt + +using BandedMatrices +using Base.Sort: Algorithm +using CliqueTrees: CliqueTrees, permutation, RCMMD +using CliqueTrees.SparseArrays +using LinearAlgebra + +function BandedMatrices.symrcm(matrix::AbstractMatrix, alg::Algorithm) + BandedMatrices.symrcm(sparse(matrix), alg) +end + +function BandedMatrices.symrcm(matrix::SparseMatrixCSC{T}, alg::Algorithm) where T + order, index = permutation(matrix; alg=RCMMD(alg)) + lower = tril!(permute(matrix, order, order)) + + bandwidth = maximum(axes(lower, 2)) do j + p = last(nzrange(lower, j)) + return rowvals(lower)[p] - j + end + + banded = BandedMatrix{T}(undef, size(lower), (bandwidth, 0)) + fill!(banded.data, zero(T)) + + @inbounds for j in axes(lower, 2) + for p in nzrange(lower, j) + i = rowvals(lower)[p] + banded.data[i - j + 1, j] = nonzeros(lower)[p] + end + end + + return Symmetric(banded, :L), order, index +end + +end diff --git a/src/BandedMatrices.jl b/src/BandedMatrices.jl index 36959ba6..71ffa269 100644 --- a/src/BandedMatrices.jl +++ b/src/BandedMatrices.jl @@ -87,6 +87,7 @@ include("symbanded/BandedCholesky.jl") include("symbanded/SplitCholesky.jl") include("symbanded/bandedeigen.jl") include("symbanded/tridiagonalize.jl") +include("symbanded/symrcm.jl") include("tribanded.jl") diff --git a/src/symbanded/symrcm.jl b/src/symbanded/symrcm.jl new file mode 100644 index 00000000..6bd15530 --- /dev/null +++ b/src/symbanded/symrcm.jl @@ -0,0 +1,39 @@ +""" + symrcm(matrix[, alg::Algorithm]) + +Use the Reverse Cuthill-Mckee algorithm to reduce the bandwith of a matrix. + +```julia +julia> using BandedMatrices, CliqueTrees + +julia> matrix = [ + 0 1 1 0 0 0 0 0 + 1 0 1 0 0 1 0 0 + 1 1 0 1 1 0 0 0 + 0 0 1 0 1 0 0 0 + 0 0 1 1 0 0 1 1 + 0 1 0 0 0 0 1 0 + 0 0 0 0 1 1 0 1 + 0 0 0 0 1 0 1 0 + ]; + +julia> bandedmatrix, perm, iperm = BandedMatrices.symrcm(matrix); + +julia> bandedmatrix +8×8 LinearAlgebra.Symmetric{Int64, BandedMatrix{Int64, Matrix{Int64}, Base.OneTo{Int64}}}: + 0 1 1 0 ⋅ ⋅ ⋅ ⋅ + 1 0 1 0 1 ⋅ ⋅ ⋅ + 1 1 0 1 0 1 ⋅ ⋅ + 0 0 1 0 0 1 0 ⋅ + ⋅ 1 0 0 0 0 1 0 + ⋅ ⋅ 1 1 0 0 1 1 + ⋅ ⋅ ⋅ 0 1 1 0 1 + ⋅ ⋅ ⋅ ⋅ 0 1 1 0 + +julia> bandedmatrix == matrix[perm, perm] +true +``` +""" +function symrcm(matrix::AbstractMatrix) + symrcm(matrix, Base.Sort.DEFAULT_UNSTABLE) +end From 56696ba6fcb4fb846e4dcd8c24ccb0c043a128cb Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Wed, 5 Mar 2025 13:08:34 -0500 Subject: [PATCH 2/2] Add documentation and unit tests. --- Project.toml | 2 +- docs/src/index.md | 38 ++++++++++++++++++++++++++++++++++++ ext/CliqueTreesExt.jl | 4 ++-- src/BandedMatrices.jl | 1 + src/symbanded/symrcm.jl | 39 ++++++++++++++++++++----------------- test/test_symbanded.jl | 43 +++++++++++++++++++++++++++++++++++++++++ 6 files changed, 106 insertions(+), 21 deletions(-) diff --git a/Project.toml b/Project.toml index 0b9fba90..9dbd8be1 100644 --- a/Project.toml +++ b/Project.toml @@ -44,4 +44,4 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["Aqua", "Documenter", "GenericLinearAlgebra", "InfiniteArrays", "Random", "SparseArrays", "Test", "Quaternions"] +test = ["Aqua", "CliqueTrees", "Documenter", "GenericLinearAlgebra", "InfiniteArrays", "Random", "SparseArrays", "Test", "Quaternions"] diff --git a/docs/src/index.md b/docs/src/index.md index bdff4fa9..a041cecb 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -162,6 +162,44 @@ To loop over the nonzero elements of a BandedMatrix, you can use `colrange(A, c) Use `Symmetric(::BandedMatrix)` to work with symmetric banded matrices. +## Bandwidth minimization + +The [reverse Cuthill-McKee algorithm](https://en.wikipedia.org/wiki/Cuthill–McKee_algorithm) permutes a symmetric matrix into a band matrix with small bandwidth. +This algorithm is implemented by the function `symrcm`. + +```julia +julia> using BandedMatrices, CliqueTrees + +julia> matrix = [ + 1 0 1 0 0 0 0 1 1 + 0 1 1 0 0 1 1 1 0 + 1 1 1 0 0 0 0 0 0 + 0 0 0 1 0 0 0 1 1 + 0 0 0 0 1 0 1 1 0 + 0 1 0 0 0 1 1 0 0 + 0 1 0 0 1 1 1 0 0 + 1 1 0 1 1 0 0 1 0 + 1 0 0 1 0 0 0 0 1 + ]; + +julia> bandedmatrix, perm, invp = symrcm(matrix); + +julia> bandedmatrix +9×9 LinearAlgebra.Symmetric{Int64, BandedMatrix{Int64, Matrix{Int64}, Base.OneTo{Int64}}}: + 1 1 1 0 ⋅ ⋅ ⋅ ⋅ ⋅ + 1 1 1 1 0 ⋅ ⋅ ⋅ ⋅ + 1 1 1 0 1 1 ⋅ ⋅ ⋅ + 0 1 0 1 0 1 0 ⋅ ⋅ + ⋅ 0 1 0 1 0 1 0 ⋅ + ⋅ ⋅ 1 1 0 1 1 1 0 + ⋅ ⋅ ⋅ 0 1 1 1 0 1 + ⋅ ⋅ ⋅ ⋅ 0 1 0 1 1 + ⋅ ⋅ ⋅ ⋅ ⋅ 0 1 1 1 + +julia> bandedmatrix == matrix[perm, perm] +true +``` + ## Banded matrix interface Banded matrices go beyond the type `BandedMatrix`: one can also create diff --git a/ext/CliqueTreesExt.jl b/ext/CliqueTreesExt.jl index 8907e260..21eabdf3 100644 --- a/ext/CliqueTreesExt.jl +++ b/ext/CliqueTreesExt.jl @@ -2,7 +2,7 @@ module CliqueTreesExt using BandedMatrices using Base.Sort: Algorithm -using CliqueTrees: CliqueTrees, permutation, RCMMD +using CliqueTrees: CliqueTrees, permutation, RCMGL using CliqueTrees.SparseArrays using LinearAlgebra @@ -11,7 +11,7 @@ function BandedMatrices.symrcm(matrix::AbstractMatrix, alg::Algorithm) end function BandedMatrices.symrcm(matrix::SparseMatrixCSC{T}, alg::Algorithm) where T - order, index = permutation(matrix; alg=RCMMD(alg)) + order, index = permutation(matrix; alg=RCMGL(alg)) lower = tril!(permute(matrix, order, order)) bandwidth = maximum(axes(lower, 2)) do j diff --git a/src/BandedMatrices.jl b/src/BandedMatrices.jl index 71ffa269..bfcc7bf7 100644 --- a/src/BandedMatrices.jl +++ b/src/BandedMatrices.jl @@ -52,6 +52,7 @@ export BandedMatrix, bandwidth, BandError, band, + symrcm, Band, BandRange, bandwidths, diff --git a/src/symbanded/symrcm.jl b/src/symbanded/symrcm.jl index 6bd15530..058e7176 100644 --- a/src/symbanded/symrcm.jl +++ b/src/symbanded/symrcm.jl @@ -3,32 +3,35 @@ Use the Reverse Cuthill-Mckee algorithm to reduce the bandwith of a matrix. +# Examples ```julia julia> using BandedMatrices, CliqueTrees julia> matrix = [ - 0 1 1 0 0 0 0 0 - 1 0 1 0 0 1 0 0 - 1 1 0 1 1 0 0 0 - 0 0 1 0 1 0 0 0 - 0 0 1 1 0 0 1 1 - 0 1 0 0 0 0 1 0 - 0 0 0 0 1 1 0 1 - 0 0 0 0 1 0 1 0 + 1 0 1 0 0 0 0 1 1 + 0 1 1 0 0 1 1 1 0 + 1 1 1 0 0 0 0 0 0 + 0 0 0 1 0 0 0 1 1 + 0 0 0 0 1 0 1 1 0 + 0 1 0 0 0 1 1 0 0 + 0 1 0 0 1 1 1 0 0 + 1 1 0 1 1 0 0 1 0 + 1 0 0 1 0 0 0 0 1 ]; -julia> bandedmatrix, perm, iperm = BandedMatrices.symrcm(matrix); +julia> bandedmatrix, perm, invp = symrcm(matrix); julia> bandedmatrix -8×8 LinearAlgebra.Symmetric{Int64, BandedMatrix{Int64, Matrix{Int64}, Base.OneTo{Int64}}}: - 0 1 1 0 ⋅ ⋅ ⋅ ⋅ - 1 0 1 0 1 ⋅ ⋅ ⋅ - 1 1 0 1 0 1 ⋅ ⋅ - 0 0 1 0 0 1 0 ⋅ - ⋅ 1 0 0 0 0 1 0 - ⋅ ⋅ 1 1 0 0 1 1 - ⋅ ⋅ ⋅ 0 1 1 0 1 - ⋅ ⋅ ⋅ ⋅ 0 1 1 0 +9×9 LinearAlgebra.Symmetric{Int64, BandedMatrix{Int64, Matrix{Int64}, Base.OneTo{Int64}}}: + 1 1 1 0 ⋅ ⋅ ⋅ ⋅ ⋅ + 1 1 1 1 0 ⋅ ⋅ ⋅ ⋅ + 1 1 1 0 1 1 ⋅ ⋅ ⋅ + 0 1 0 1 0 1 0 ⋅ ⋅ + ⋅ 0 1 0 1 0 1 0 ⋅ + ⋅ ⋅ 1 1 0 1 1 1 0 + ⋅ ⋅ ⋅ 0 1 1 1 0 1 + ⋅ ⋅ ⋅ ⋅ 0 1 0 1 1 + ⋅ ⋅ ⋅ ⋅ ⋅ 0 1 1 1 julia> bandedmatrix == matrix[perm, perm] true diff --git a/test/test_symbanded.jl b/test/test_symbanded.jl index 2cf802a3..a8362025 100644 --- a/test/test_symbanded.jl +++ b/test/test_symbanded.jl @@ -3,6 +3,7 @@ module TestSymBanded using ArrayLayouts using BandedMatrices import BandedMatrices: MemoryLayout, SymmetricLayout, HermitianLayout, BandedColumns, isbanded +using CliqueTrees: CliqueTrees using FillArrays using GenericLinearAlgebra using LinearAlgebra @@ -421,4 +422,46 @@ end end end +@testset "symrcm" begin + matrix = Int16[ + 1 0 1 0 0 0 0 1 1 + 0 1 1 0 0 1 1 1 0 + 1 1 1 0 0 0 0 0 0 + 0 0 0 1 0 0 0 1 1 + 0 0 0 0 1 0 1 1 0 + 0 1 0 0 0 1 1 0 0 + 0 1 0 0 1 1 1 0 0 + 1 1 0 1 1 0 0 1 0 + 1 0 0 1 0 0 0 0 1 + ] + + bandedmatrix, perm, invp = symrcm(matrix) + @test invp == invperm(perm) + @test bandedmatrix == matrix[perm, perm] + @test eltype(bandedmatrix) == eltype(matrix) + + # LFAT5 + matrix = Float32[ + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + 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 + -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 + 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 + ] + + bandedmatrix, perm, invp = symrcm(matrix) + @test invp == invperm(perm) + @test bandedmatrix == matrix[perm, perm] + @test eltype(bandedmatrix) == eltype(matrix) +end + end # module