fix multiplication of BlockSkylineMatrix

Author: dlfivefifty

src/linalg.jl

@@ -103,23 +103,25 @@ function add_bandwidths(A::AbstractBlockBandedMatrix,B::AbstractBlockBandedMatri
     Al,Au = colblockbandwidths(A)
     Bl,Bu = colblockbandwidths(B)
 
-    l = copy(Al)
-    u = copy(Au)
+    l = copy(Bl)
+    u = copy(Bu)
 
-    for (v,Bv) in [(l,Bl),(u,Bu)]
+    for (v,Av) in [(l,Al),(u,Au)]
         n = length(v)
         for i = 1:n
-            sel = max(i-Au[i],1):min(i+Al[i],n)
+            sel = max(i-Bu[i],1):min(i+Bl[i],length(Av))
             isempty(sel) && continue
-            v[i] += maximum(Bv[sel])
+            v[i] += maximum(Av[sel])
         end
     end
 
     l,u
 end
 
-add_bandwidths(A::BlockBandedMatrix,B::BlockBandedMatrix) =
-    colblockbandwidths(A) .+ colblockbandwidths(B)
+function add_bandwidths(A::BlockBandedMatrix,B::BlockBandedMatrix)
+    l,u = blockbandwidths(A) .+ blockbandwidths(B)
+    Fill(l,nblocks(B,2)), Fill(u,nblocks(B,2))
+end
 
 function similar(M::MatMulMat{<:AbstractBlockBandedLayout,<:AbstractBlockBandedLayout}, ::Type{T}) where T
     A,B = M.factors

test/test_linalg.jl

@@ -32,29 +32,29 @@ end
     @test unsafe_load(pointer(V)) == 46
     @test unsafe_load(pointer(V) + stride(V,2)*sizeof(Float64)) == 53
 
-x = randn(size(A,2))
-    @test A*x == (similar(x) .= Mul(A,x)) ≈ Matrix(A)*x
+    x = randn(size(A,2))
+        @test A*x == (similar(x) .= Mul(A,x)) ≈ Matrix(A)*x
 
-z = randn(size(A,2)) + im*randn(size(A,2))
-A*z == (similar(z) .= Mul(A,z)) ≈ Matrix(A)*z
+    z = randn(size(A,2)) + im*randn(size(A,2))
+    A*z == (similar(z) .= Mul(A,z)) ≈ Matrix(A)*z
 
-Matrix(A)*z
-z[1]
-view(A,Block(1,1)) * z[1] + view(A,Block(1,2))*z[2:3]
+    Matrix(A)*z
+    z[1]
+    view(A,Block(1,1)) * z[1] + view(A,Block(1,2))*z[2:3]
 
-A*z
+    A*z
 
-X = randn(size(A))
-    @test A*X == (similar(X) .= Mul(A,X)) ≈ Matrix(A)*X
-    @test X*A == (similar(X) .= Mul(X,A)) ≈  Matrix(X)*A
+    X = randn(size(A))
+        @test A*X == (similar(X) .= Mul(A,X)) ≈ Matrix(A)*X
+        @test X*A == (similar(X) .= Mul(X,A)) ≈  Matrix(X)*A
 
 
-Z = randn(size(A)) + im*randn(size(A))
+    Z = randn(size(A)) + im*randn(size(A))
 
-A*Z
+    A*Z
 
 
-v = fill(1.0,4)
+    v = fill(1.0,4)
     U = UpperTriangular(view(A, Block(N,N)))
     @test Matrix(U) == U
     w = Matrix(U) \ v
@@ -72,8 +72,6 @@ v = fill(1.0,4)
     @test v == fill(1.0,size(A,1))
     @test ldiv!(U, v) === v
     @test v ≈ w
-
-
 end
 
 @testset "BandedBlockBandedMatrix linear algebra" begin
@@ -155,3 +153,26 @@ end
     @time BLAS.axpy!(1.0, A, B)
     @test B ≈ AB
 end
+
+@testset "Rectangular block *" begin
+    A = BlockBandedMatrix{Float64}(undef, (Ones{Int}(2), Ones{Int}(2)), (0,2))
+        A.data .= randn.()
+    B = BlockBandedMatrix{Float64}(undef, (Ones{Int}(2), Ones{Int}(3)), (0,2))
+        B.data .= randn.()
+
+    @test A*B ≈ Matrix(A)*Matrix(B)
+
+
+    rows = rand(1:10, 5)
+    l = rand(0:2, 5)
+    u = rand(0:2, 5)
+
+    m = sum(rows)
+
+    A = BlockSkylineMatrix(Zeros(m,m), (rows,rows), (l,u))
+    A.data .= randn.()
+
+    B = BlockSkylineMatrix(Zeros(m,m+1), (rows,[rows;1]), ([l;1],[u;1]))
+    B.data .= randn.()
+    @test A*B ≈ Matrix(A)*Matrix(B)
+end