From ec1ac49c5c6069c506bce384333e6baa77f5eecf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lars=20G=C3=B6ttgens?= Date: Tue, 26 Sep 2023 15:35:05 +0200 Subject: [PATCH 1/2] Add sparse block diagonal matrices --- docs/src/sparse/intro.md | 1 + src/Sparse/Matrix.jl | 47 +++++++++++++++++++++++++++++++++++++++- 2 files changed, 47 insertions(+), 1 deletion(-) diff --git a/docs/src/sparse/intro.md b/docs/src/sparse/intro.md index c5922ba2e0..6015c9cb1e 100644 --- a/docs/src/sparse/intro.md +++ b/docs/src/sparse/intro.md @@ -127,6 +127,7 @@ There are special constructors: identity_matrix(::Type{SMat}, ::Ring, ::Int) zero_matrix(::Type{SMat}, ::Ring, ::Int) zero_matrix(::Type{SMat}, ::Ring, ::Int, ::Int) +block_diagonal_matrix(xs::Vector{<:SMat{T}}) where {T} ``` Slices: ```@docs diff --git a/src/Sparse/Matrix.jl b/src/Sparse/Matrix.jl index 54912301da..912bdce5b3 100644 --- a/src/Sparse/Matrix.jl +++ b/src/Sparse/Matrix.jl @@ -1235,9 +1235,54 @@ Return a sparse $n$ times $m$ zero matrix over $R$. """ zero_matrix(::Type{SMat}, R::Ring, n::Int, m::Int) = sparse_matrix(R, n, m) +similar(A::SMat) = sparse_matrix(base_ring(A), nrows(A), ncols(A)) + +similar(A::SMat, m::Int, n::Int) = sparse_matrix(base_ring(A), m, n) + +################################################################################ +# +# (Block) diagonal matrices +# +################################################################################ + +function diagonal_matrix(V::Vector{<:SMat{T}}) where {T} + return block_diagonal_matrix(V) +end + +function diagonal_matrix(x::SMat{T}, xs::SMat{T}...) where {T} + return block_diagonal_matrix([x, xs...]) +end + +@doc raw""" + block_diagonal_matrix(xs::Vector{SMat}) + +Return the block diagonal matrix with the matrices in `xs` on the diagonal. +Requires all blocks to have the same base ring. +""" +function block_diagonal_matrix(xs::Vector{<:SMat{T}}) where {T} + @req length(xs) > 0 "At least one matrix is required" + R = base_ring(xs[1]) + @req all(x -> base_ring(x) == R, xs) "All matrices must have the same base ring" + rows = sum(nrows(N) for N in xs) + cols = sum(ncols(N) for N in xs) + M = similar(xs[1], rows, cols) + i_offset = 0 + j_offset = 0 + for x in xs + for (i, x_row) in enumerate(x) + M_row = sparse_row(R, x_row.pos .+ j_offset, x_row.values) + setindex!(M, M_row, i_offset + i) + end + i_offset += nrows(x) + j_offset += ncols(x) + end + + return M +end + ################################################################################ # -# Julias concatenatin syntax +# Julias concatenation syntax # ################################################################################ From d44095a27455e2babe1a4b4a790c4db9cc6aef65 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Lars=20G=C3=B6ttgens?= Date: Tue, 26 Sep 2023 15:51:23 +0200 Subject: [PATCH 2/2] Add tests --- src/Sparse/Matrix.jl | 2 +- src/Sparse/Row.jl | 2 +- src/Sparse/ZZRow.jl | 2 +- test/Sparse/Matrix.jl | 11 +++++++++++ 4 files changed, 14 insertions(+), 3 deletions(-) diff --git a/src/Sparse/Matrix.jl b/src/Sparse/Matrix.jl index 912bdce5b3..c510ebaa63 100644 --- a/src/Sparse/Matrix.jl +++ b/src/Sparse/Matrix.jl @@ -1241,7 +1241,7 @@ similar(A::SMat, m::Int, n::Int) = sparse_matrix(base_ring(A), m, n) ################################################################################ # -# (Block) diagonal matrices +# Block diagonal matrices # ################################################################################ diff --git a/src/Sparse/Row.jl b/src/Sparse/Row.jl index 79a10c479f..8b12dad349 100644 --- a/src/Sparse/Row.jl +++ b/src/Sparse/Row.jl @@ -89,7 +89,7 @@ end Constructs the sparse row $(a_i)_i$ over $R$ with $a_{i_j} = x_j$, where $J = (i_j)_j$ and $V = (x_j)_j$. """ -function sparse_row(R::Ring, pos::Vector{Int}, val::Vector{T}; sort::Bool = true) where T +function sparse_row(R::Ring, pos::Vector{Int}, val::AbstractVector{T}; sort::Bool = true) where T if sort p = sortperm(pos) pos = pos[p] diff --git a/src/Sparse/ZZRow.jl b/src/Sparse/ZZRow.jl index 64074b63e1..10ff5cd026 100644 --- a/src/Sparse/ZZRow.jl +++ b/src/Sparse/ZZRow.jl @@ -249,7 +249,7 @@ function sparse_row(R::ZZRing, A::Vector{Tuple{Int64, ZZRingElem}}; sort::Bool = return SRow(R, l, a) end -function sparse_row(R::ZZRing, pos::Vector{Int64}, val::Vector{ZZRingElem}; sort::Bool = true) +function sparse_row(R::ZZRing, pos::Vector{Int64}, val::AbstractVector{ZZRingElem}; sort::Bool = true) if sort p = sortperm(pos) pos = pos[p] diff --git a/test/Sparse/Matrix.jl b/test/Sparse/Matrix.jl index 47dd6c0b6b..a21589269c 100644 --- a/test/Sparse/Matrix.jl +++ b/test/Sparse/Matrix.jl @@ -264,6 +264,17 @@ using Hecke.SparseArrays @test D == sparse_matrix(FlintZZ, [0 0 0; 0 0 0; 0 0 0; 0 0 0]); @test D == sparse_matrix(FlintZZ, 4, 3) + # Block diag matrix + + D1 = sparse_matrix(FlintZZ, [1 5; 0 1]) + D2 = sparse_matrix(FlintZZ, [2 3 8; 4 0 0]) + D = @inferred block_diagonal_matrix([D1, D2]) + @test D == sparse_matrix(FlintZZ, [1 5 0 0 0; 0 1 0 0 0; 0 0 2 3 8; 0 0 4 0 0]) + D = @inferred diagonal_matrix([D1, D2]) + @test D == sparse_matrix(FlintZZ, [1 5 0 0 0; 0 1 0 0 0; 0 0 2 3 8; 0 0 4 0 0]) + D = @inferred diagonal_matrix(D1, D2) + @test D == sparse_matrix(FlintZZ, [1 5 0 0 0; 0 1 0 0 0; 0 0 2 3 8; 0 0 4 0 0]) + # Concatenation syntax D = sparse_matrix(FlintZZ, [1 5 3; 0 -10 0; 0 1 0])