Add sparse block diagonal matrices

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

src/Sparse/Matrix.jl

@@ -1235,9 +1235,54 @@
 """
 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
 #
 ################################################################################
 

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]

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]

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])