error in broadcasting with fewer bands (#305)

Author: jishnub

src/generic/Band.jl

@@ -73,7 +73,7 @@ BandError(A::AbstractMatrix) = BandError(A, max(size(A)...)-1)
 
 function showerror(io::IO, e::BandError)
     A, i = e.A, e.i
-    print(io, "attempt to access $(typeof(A)) with bandwidths " *
+    print(io, "BandError: attempt to access $(typeof(A)) with bandwidths " *
               "($(bandwidth(A, 1)), $(bandwidth(A, 2))) at band $i")
 end
 

src/generic/broadcast.jl

@@ -78,7 +78,7 @@ function checkzerobands(dest, f, A::AbstractMatrix)
     l, u = bandwidths(A)
 
     if (l,u) ≠ (d_l,d_u)
-        for j = 1:n
+        for j = rowsupport(A)
             for k = max(1,j-u) : min(j-d_u-1,m)
                 iszero(f(A[k,j])) || throw(BandError(dest,j-k))
             end
@@ -93,21 +93,12 @@ function _banded_broadcast!(dest::AbstractMatrix, f, src::AbstractMatrix{T}, _1,
     z = f(zero(T))
     iszero(z) || checkbroadcastband(dest, size(src), bandwidths(broadcasted(f, src)))
     m,n = size(dest)
+    m == n == 0 && return dest
 
     d_l, d_u = bandwidths(dest)
     s_l, s_u = bandwidths(src)
-    if d_l < min(s_l,m-1)
-        # check zeros
-        for j = 1:n, k = max(1,j+d_l+1):min(j+s_l,j+d_l,m)
-            iszero(f(inbands_getindex(src, k, j))) || throw(BandError(dest))
-        end
-    end
-    if d_u < min(s_u,n-1)
-        # check zeros
-        for j = 1:n, k = max(1,j-d_u,j-s_u):min(j-d_u-1,m)
-            iszero(f(inbands_getindex(src, k, j))) || throw(BandError(dest))
-        end
-    end
+
+    checkzerobands(dest, f, src)
 
     _banded_broadcast_anylayout!(dest, src, f, z, (s_l, s_u), (d_l, d_u), m)
 
@@ -252,13 +243,15 @@ end
     min_su_du = min(s_u, d_u)
     min_sl_dl = min(s_l, d_l)
     for j = rowsupport(dest)
-        for k = max(1,j-d_u):min(j-s_u-1,m)
+        # if s_u < d_u, set extra bands in dest above the diagonal to zero
+        for k = max(1,j-d_u):min(j+min(d_l, -(s_u+1)),m)
             inbands_setindex!(dest, z, k, j)
         end
 
         _banded_broadcast_loop_overlap!(dest, src, f, (min_sl_dl, min_su_du), m, j)
 
-        for k = max(1,j+s_l+1):min(j+d_l,m)
+        # if s_l < d_l, set extra bands in dest below the diagonal to zero
+        for k = max(1,j+max(s_l+1, -d_u)):min(j+d_l,m)
             inbands_setindex!(dest, z, k, j)
         end
     end
@@ -276,12 +269,17 @@ function _banded_broadcast!(dest::AbstractMatrix, f, (src,x)::Tuple{AbstractMatr
     z = f(zero(T), x)
     iszero(z) || checkbroadcastband(dest, size(src), bandwidths(broadcasted(f, src,x)))
     m,n = size(dest)
+    m == n == 0 && return dest
+
+    f_x = Base.Fix2(f, x)
+
+    # if dest has fewer bands (either above or below the diagonal) than dest,
+    # then f.(x, dest) must be zero in these bands
+    checkzerobands(dest, f_x, src)
 
     d_l, d_u = bandwidths(dest)
     s_l, s_u = bandwidths(src)
-    (d_l ≥ min(s_l,m-1) && d_u ≥ min(s_u,n-1)) || throw(BandError(dest))
 
-    f_x = Base.Fix2(f, x)
     _banded_broadcast_anylayout!(dest, src, f_x, z, (s_l, s_u), (d_l, d_u), m)
 
     dest
@@ -291,12 +289,19 @@ function _banded_broadcast!(dest::AbstractMatrix, f, (x,src)::Tuple{Number,Abstr
     z = f(x, zero(T))
     iszero(z) || checkbroadcastband(dest, size(src), bandwidths(broadcasted(f, x,src)))
     m,n = size(dest)
+    m == n == 0 && return dest
+
+    f_x = Base.Fix1(f, x)
+
+    # if dest has fewer bands (either above or below the diagonal) than dest,
+    # then f.(x, src) must be zero in these bands
+    checkzerobands(dest, f_x, src)
 
     d_l, d_u = bandwidths(dest)
     s_l, s_u = bandwidths(src)
-    (d_l ≥ min(s_l,m-1) && d_u ≥ min(s_u,n-1)) || throw(BandError(dest))
 
-    f_x = Base.Fix1(f, x)
+    checkzerobands(dest, f, src)
+
     _banded_broadcast_anylayout!(dest, src, f_x, z, (s_l, s_u), (d_l, d_u), m)
 
     dest

test/test_broadcasting.jl

@@ -49,6 +49,17 @@ import BandedMatrices: BandedStyle, BandedRows
         end
     end
 
+    function test_empty(f!)
+        # the function f! must not set the RHS to zero for the error check to work
+        @testset "empty dest" begin
+            D = brand(4, 4, -2, 1) # empty
+            B = brand(size(D)...,1,2) # non-empty and non-zero bands
+            @test_throws BandError f!(D, B)
+            f!(D, zero(B))
+            @test all(iszero, D)
+        end
+    end
+
     @testset "identity" begin
         n = 100
         A = brand(n,n,2,2)
@@ -87,6 +98,16 @@ import BandedMatrices: BandedStyle, BandedRows
         B = brand(n,n,1,-1)
         A .= B
         @test A == B
+
+        @testset "empty dest" begin
+            test_empty((D,B) -> D .= B)
+        end
+
+        @testset "adjtrans" begin
+            @testset "empty dest" begin
+                test_empty((D,B) -> D' .= B')
+            end
+        end
     end
 
     @testset "lmul!/rmul!" begin
@@ -187,6 +208,11 @@ import BandedMatrices: BandedStyle, BandedRows
         B .= 2.0 .\ A
         @test B == 2.0 \ A == 2.0 \ Matrix(A)
 
+        @testset "empty dest" begin
+            test_empty((D,B) -> D .= 2 .* B)
+            test_empty((D,B) -> D .= B .* 2)
+        end
+
         @testset "trans-adj" begin
             A = brand(5,5,1,1)
             Ã = copy(A)
@@ -207,6 +233,15 @@ import BandedMatrices: BandedStyle, BandedRows
             @test A == 2Ã
             rmul!(A', 1/2)
             @test A == Ã
+
+            @testset "empty dest" begin
+                test_empty((D,B) -> D' .= 2 .* B')
+                test_empty((D,B) -> D' .= B' .* 2)
+                test_empty((D,B) -> D' .= 2 .* B' .* 2)
+                test_empty((D,B) -> D' .= 2 .* B)
+                test_empty((D,B) -> D' .= B .* 2)
+                test_empty((D,B) -> D' .= 2 .* B .* 2)
+            end
         end
     end