Avoid recursion in _banded_muladd! for vectors (#293)

* Avoid recursion in _banded_muladd for vector

* Handle zero-sized arrays in _banded_gbmv

* Add comma in comment

* Add tests

* condense _banded_gbmv! branches

Author: jishnub

src/generic/matmul.jl

@@ -26,11 +26,18 @@ banded_gbmv!(tA, α, A, x, β, y) =
 
 
 @inline function _banded_gbmv!(tA, α, A, x, β, y)
-    if x ≡ y
-        banded_gbmv!(tA, α, A, copy(x), β, y)
+    #= Some BLAS implementations throw warnings
+    with zero-sized arrays, so we handle
+    these cases separately.
+    =#
+    length(y) == 0 && return y
+    if length(x) == 0
+        _fill_lmul!(β, y)
     else
-        banded_gbmv!(tA, α, A, x, β, y)
+        xc = x ≡ y ? copy(x) : x
+        banded_gbmv!(tA, α, A, xc, β, y)
     end
+    return y
 end
 
 
@@ -43,11 +50,15 @@ function _banded_muladd!(α::T, A, x::AbstractVector, β, y) where T
     l, u = bandwidths(A)
     if -l > u # no bands
         _fill_lmul!(β, y)
-    elseif l < 0
-        _banded_muladd!(α, view(A, :, 1-l:n), view(x, 1-l:n), β, y)
-    elseif u < 0
+    elseif l < 0 # with u >= -l > 0, that is, all bands lie above the diagonal
+        # E.g. (l,u) = (-1,2)
+        # set lview = 0 and uview = u + l >= 0
+        _banded_gbmv!('N', α, view(A, :, 1-l:n), view(x, 1-l:n), β, y)
+    elseif u < 0 # with -l <= u < 0, that is, all bands lie below the diagnoal.
+        # E.g. (l,u) = (2,-1)
+        # set lview = l + u >= 0 and uview = 0
         y[1:-u] .= zero(T)
-        _banded_muladd!(α, view(A, 1-u:m, :), x, β, view(y, 1-u:m))
+        _banded_gbmv!('N', α, view(A, 1-u:m, :), x, β, view(y, 1-u:m))
         y
     else
         _banded_gbmv!('N', α, A, x, β, y)

test/test_banded.jl

@@ -111,6 +111,15 @@ Base.similar(::MyMatrix, ::Type{T}, m::Int, n::Int) where T = MyMatrix{T}(undef,
             @test A*v ≈ Matrix(A)*v
             @test A'*w ≈ Matrix(A)'*w
         end
+
+        @testset "empty" begin
+            let B=BandedMatrix((0=>ones(0),), (10,0)), v = ones(size(B,2))
+                @test B * v == zeros(size(B,1))
+            end
+            let B=BandedMatrix((0=>ones(0),), (0,10)), v = ones(size(B,2))
+                @test B * v == zeros(size(B,1))
+            end
+        end
     end
 
     @testset "Banded * Dense" begin