Fix rowsupport/colsupport over ranges for rectangular block-banded ma…

Project.toml

@@ -1,6 +1,6 @@
 name = "BlockBandedMatrices"
 uuid = "ffab5731-97b5-5995-9138-79e8c1846df0"
-version = "0.13.8"
+version = "0.13.9"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

src/AbstractBlockBandedMatrix.jl

@@ -161,6 +161,7 @@ end
 
 @inline function blockbanded_colstop(A, i::Integer)
     CS = blockcolstop(A,findblock(axes(A,2),i))
+    Int(CS) > blocksize(A, 1) && return last(axes(A,1)) # empty column past last block: keep colstop monotone
     CS in blockaxes(axes(A,1), 1) || return 0
     last(axes(A,1)[CS])
 end
@@ -178,6 +179,7 @@ end
 
 @inline function blockbanded_rowstop(A, i::Integer)
     CS = blockrowstop(A,findblock(axes(A,1),i))
+    Int(CS) > blocksize(A, 2) && return last(axes(A,2)) # empty row past last block: keep rowstop monotone
     CS in blockaxes(axes(A,2), 1) || return 0
     last(axes(A,2)[CS])
 end

test/test_bandedblockbanded.jl

@@ -578,6 +578,22 @@ import ArrayLayouts: RangeCumsum
         @test copy(transpose(A)) isa Transpose
         @test copy(transpose(A)) == transpose(copy(A))
     end
+
+    @testset "rectangular row vector mul" begin
+        for (m,n) in ((9,5), (5,9)), lu in ((2,0), (0,2), (1,1), (1,-1)), λμ in ((1,1), (0,2))
+            S = BandedBlockBandedMatrix(randn(sum(1:m), sum(1:n)), 1:m, 1:n, lu, λμ)
+            a, b = BlockVector(randn(sum(1:m)), 1:m), BlockVector(randn(sum(1:n)), 1:n)
+            @test all(k -> rowsupport(S,k) ⊆ rowsupport(S,axes(S,1)), axes(S,1))
+            @test all(j -> colsupport(S,j) ⊆ colsupport(S,axes(S,2)), axes(S,2))
+            @test a' * S ≈ Vector(a)' * Matrix(S)
+            @test Vector(a)' * S ≈ Vector(a)' * Matrix(S)
+            @test transpose(a) * S ≈ transpose(Vector(a)) * Matrix(S)
+            @test transpose(S) * a ≈ transpose(Matrix(S)) * Vector(a)
+            @test S' * a ≈ Matrix(S)' * Vector(a)
+            @test S * b ≈ Matrix(S) * Vector(b)
+            @test S * Matrix(S)' ≈ Matrix(S) * Matrix(S)'
+        end
+    end
 end
 
 if false # turned off since tests have check-bounds=yes

test/test_blockbanded.jl

@@ -252,6 +252,18 @@ import BlockBandedMatrices: MemoryLayout, ColumnMajor, BroadcastStyle,
         s = BlockVector(1:21, 1:6)
         @test BlockBandedMatrices.subblockbandwidths(s) == (5, 0)
     end
+
+    @testset "rectangular row vector mul" begin
+        for (m,n) in ((9,5), (5,9)), lu in ((2,0), (0,2), (1,-1))
+            A = BlockBandedMatrix(randn(sum(1:m), sum(1:n)), 1:m, 1:n, lu)
+            a = BlockVector(randn(sum(1:m)), 1:m)
+            @test all(k -> rowsupport(A,k) ⊆ rowsupport(A,axes(A,1)), axes(A,1))
+            @test all(j -> colsupport(A,j) ⊆ colsupport(A,axes(A,2)), axes(A,2))
+            @test a' * A ≈ Vector(a)' * Matrix(A)
+            @test Vector(a)' * A ≈ Vector(a)' * Matrix(A)
+            @test A * Matrix(A)' ≈ Matrix(A) * Matrix(A)'
+        end
+    end
 end
 
 end # module