Implement vcat(::AbstractBandedMatrix...)

ext/BandedMatricesSparseArraysExt.jl

@@ -2,37 +2,53 @@
 
 using BandedMatrices
 using BandedMatrices: _banded_rowval, _banded_colval, _banded_nzval
-using SparseArrays
+using SparseArrays, FillArrays
 import SparseArrays: sparse
 
 function sparse(B::BandedMatrix)
 	sparse(_banded_rowval(B), _banded_colval(B), _banded_nzval(B), size(B)...)
 end
 
 function BandedMatrices.bandwidths(A::SparseMatrixCSC)
-    l,u = -size(A,1),-size(A,2)
-
-    m,n = size(A)
+    l = u = -max(size(A,1),size(A,2))
+    n = size(A)[2]
     rows = rowvals(A)
     vals = nonzeros(A)
+
+    if isempty(vals)
+        return bandwidths(Zeros(1))
+    end
+
     for j = 1:n
         for ind in nzrange(A, j)
             i = rows[ind]
             # We skip non-structural zeros when computing the
             # bandwidths.
             iszero(vals[ind]) && continue
-            ij = abs(i-j)
-            if i ≥ j
-                l = max(l, ij)
-                u = max(u, -ij)
-            elseif i < j
-                l = max(l, -ij)
-                u = max(u, ij)
-            end
+            u = max(u, j-i)
+            l = max(l, i-j)
         end
     end
 
     l,u
 end
 
+#Treat as n x 1 matrix
+function BandedMatrices.bandwidths(A::SparseVector)
+    l = u = -size(A,1)
+    rows = rowvals(A)
+
+    if isempty(rows)
+        return bandwidths(Zeros(1))
+    end
+
+    for i in rows
+        iszero(i) && continue
+        u = max(u, 1-i)
+        l = max(l, i-1)
+    end
+
+    l,u
+end
+
 end

src/generic/AbstractBandedMatrix.jl

@@ -401,3 +401,57 @@
         throw(ArgumentError("dimension must be ≥ 1, got $dims"))
     end
 end
+
+###
+# vcat
+###
+
+function LinearAlgebra.vcat(x::AbstractBandedMatrix...)
+    #avoid unnecessary steps for singleton
+    if length(x) == 1
+        return x[1]
+    end
+
+    #instantiate the returned banded matrix with zeros and required bandwidths/dimensions
+    m = size(x[1], 2)
+    l,u = -m, typemin(Int64)
+    n = 0
+    isempty = true
+
+    #Check for dimension error and calculate bandwidths
+    for A in x
+        if size(A, 2) != m
+            sizes = Tuple(size(b, 2) for b in x)
+            throw(DimensionMismatch("number of columns of each matrix must match (got $sizes)"))
+        end
+
+        l_A, u_A = bandwidths(A)
+        if l_A + u_A >= 0
+            isempty = false
+            u = max(u, min(m - 1, u_A) - n)
+            l = max(l, min(size(A, 1) - 1, l_A) + n)
+        end
+
+        n += size(A, 1)
+    end
+
+    type = promote_type(eltype.(x)...)
+    if isempty
+        return BandedMatrix{type}(undef, (n, m), bandwidths(Zeros(1)))
+    end
+    ret = BandedMatrix(Zeros{type}(n, m), (l, u))
+
+    #Populate the banded matrix
+    row_offset = 0
+    for A in x
+        n_A = size(A, 1)
+
+        for i = 1:n_A, j = rowrange(A, i)
+            ret[row_offset + i, j] = A[i, j]
+        end
+
+        row_offset += n_A
+    end
+
+    ret
+end

src/interfaceimpl.jl

@@ -56,6 +56,28 @@
 sublayout(::AbstractTridiagonalLayout, ::Type{<:Tuple{AbstractUnitRange{Int},AbstractUnitRange{Int}}}) =
     BandedLayout()
 
+#Implement bandwidths for OneElement structure
+function bandwidths(o::OneElement)
+    pos = FillArrays.nzind(o)
+    if length(pos) == 1
+        n = length(o)
+        if pos[1] > n
+            bandwidths(Zeros(o))
+        else
+            (pos[1] - 1, -pos[1] + 1)
+        end
+    elseif length(pos) == 2
+        n,m = size(o)
+        if pos[1] > n || pos[2] > m
+            bandwidths(Zeros(o))
+        else
+            (pos[1]-pos[2],pos[2]-pos[1])
+        end
+    end
+end
+
+LinearAlgebra.vcat(x::Union{OneElement, AbstractBandedMatrix}...) = vcat(BandedMatrix.(x)...)
+
 ###
 # rot180
 ###

test/runtests.jl

@@ -27,3 +27,4 @@ include("test_tribanded.jl")
 include("test_interface.jl")
 include("test_miscs.jl")
 include("test_sum.jl")
+include("test_cat.jl")

test/test_cat.jl

@@ -0,0 +1,37 @@
+module TestCat
+
+using BandedMatrices, LinearAlgebra, Test, Random, FillArrays, SparseArrays
+
+@testset "vcat" begin
+    @testset "banded matrices" begin
+        a = BandedMatrix(0 => 1:2)
+        @test vcat(a) == a
+
+        b = BandedMatrix(0 => 1:3,-1 => 1:2, -2 => 1:1)
+        @test_throws DimensionMismatch vcat(a,b)
+
+        c = BandedMatrix(0 => [1.0, 2.0, 3.0], 1 => [1.0, 2.0], 2 => [1.0])
+        @test eltype(vcat(b, c)) == Float64
+        @test vcat(b, c) == vcat(Matrix(b), Matrix(c))
+
+        for i in ((1,2), (-3,4), (0,-1))
+            a = BandedMatrix(ones(Float64, rand(1:10), 5), i)
+            b = BandedMatrix(ones(Int64, rand(1:10), 5), i)
+            c = BandedMatrix(ones(Int32, rand(1:10), 5), i)
+            d = vcat(a, b, c)
+            sd = vcat(sparse(a), sparse(b), sparse(c))
+            @test eltype(d) == Float64
+            @test d == sd
+            @test bandwidths(d) == bandwidths(sd)
+        end
+    end
+
+    @testset "one element" begin
+        n = rand(3:20)
+        x,y = OneElement(1, (1,1), (1,n)), OneElement(1, (1,n), (1,n))
+        b = BandedMatrix((0 => ones(n-2), 1 => -2ones(n - 2), 2 => ones(n - 2)), (n-2, n))
+        @test vcat(x,b,y) == Tridiagonal([ones(n - 2); 0], [1 ; -2ones(n - 2); 1], [0; ones(n - 2)])
+    end
+end
+
+end

test/test_interface.jl

@@ -1,10 +1,10 @@
 module TestInterface
 
-using BandedMatrices, LinearAlgebra, ArrayLayouts, FillArrays, Test
+using BandedMatrices, LinearAlgebra, ArrayLayouts, FillArrays, Test, Random
 import BandedMatrices: isbanded, AbstractBandedLayout, BandedStyle,
                         BandedColumns, bandeddata
 import ArrayLayouts: OnesLayout, UnknownLayout
-using InfiniteArrays
+using InfiniteArrays, SparseArrays
 
 struct PseudoBandedMatrix{T} <: AbstractMatrix{T}
     data::Array{T}
@@ -310,6 +310,18 @@ end
     @test layout_getindex(T,1:10,1:10) isa BandedMatrix
 end
 
+@testset "OneElement" begin
+    o = OneElement(1, 3, 5)
+    @test bandwidths(o) == (2,-2)
+    n,m = rand(1:10,2)
+    o = OneElement(1, (rand(1:n),rand(1:m)), (n, m))
+    @test bandwidths(o) == bandwidths(sparse(o))
+    o = OneElement(1, (n+1,m+1), (n, m))
+    @test bandwidths(o) == bandwidths(Zeros(o))
+    o = OneElement(1, 6, 5)
+    @test bandwidths(o) == bandwidths(Zeros(o))
+end
+
 @testset "rot180" begin
     A = brand(5,5,1,2)
     R = rot180(A)

test/test_miscs.jl

@@ -50,8 +50,17 @@ import BandedMatrices: _BandedMatrix, DefaultBandedMatrix
             @test bA isa BandedMatrix
             @test bA == A
             @test bandwidths(bA) == min.((l,u),9)
+            v = sparsevec(brand(10, 1, l, u))
+            @test bandwidths(v) == (l, min(0, u))
         end
 
+        l, u = -1, 0
+        A = brand(10, 10, l, u)
+        sA = sparse(A)
+        @test bandwidths(sA) == bandwidths(Zeros(1))
+        v = sparsevec(brand(10, 1, l, u))
+        @test bandwidths(v) == bandwidths(Zeros(1))
+
         for diags = [(-1 => ones(Int, 5),),
                      (-2 => ones(Int, 5),),
                      (2 => ones(Int, 5),),