Merge branch 'master' of https://github.com/JuliaMatrices/BlockBandedMatrices.jl

Author: dlfivefifty

src/BlockBandedMatrices.jl

@@ -15,7 +15,7 @@ import BandedMatrices: isbanded, bandwidths, bandwidth, banded_getindex, colrang
                         BandedColumnMajor,
                         BandedSubBandedMatrix, bandeddata, tribandeddata,
                         _BandedMatrix, colstart, colstop, rowstart, rowstop,
-                        BandedStyle,
+                        BandedStyle, _fill_lmul!,
                         _banded_colval, _banded_rowval, _banded_nzval # for sparse
 
 import Base: getindex, setindex!, checkbounds, @propagate_inbounds, convert,

src/linalg.jl

@@ -40,7 +40,7 @@ function materialize!(M::MatMulVecAdd{<:AbstractBlockBandedLayout,<:AbstractStri
     y = PseudoBlockArray(y_in, BlockSizes((cumulsizes(blocksizes(A),1),)))
     x = PseudoBlockArray(x_in, BlockSizes((cumulsizes(blocksizes(A),2),)))
 
-    lmul!(β, y)
+    _fill_lmul!(β, y)
 
     for J = Block.(1:nblocks(A,2))
         for K = blockcolrange(A,J)
@@ -52,7 +52,7 @@ end
 
 function materialize!(M::MatMulMatAdd{<:AbstractBlockBandedLayout,<:AbstractBlockBandedLayout,<:AbstractBlockBandedLayout})
     α, A, X, β, Y = M.α, M.A, M.B, M.β, M.C
-    lmul!(β, Y)
+    _fill_lmul!(β, Y)
     for J=Block(1):Block(nblocks(X,2)),
             N=blockcolrange(X,J), K=blockcolrange(A,N)
         view(Y,K,J) .= α .* Mul( view(A,K,N), view(X,N,J)) .+ view(Y,K,J)
@@ -62,7 +62,7 @@ end
 
 function materialize!(M::MatMulMatAdd{<:AbstractBlockBandedLayout,<:AbstractColumnMajor,<:AbstractColumnMajor})
     α, A, X_in, β, Y_in = M.α, M.A, M.B, M.β, M.C
-    lmul!(β, Y_in)
+    _fill_lmul!(β, Y_in)
     X = PseudoBlockArray(X_in, BlockSizes((cumulsizes(blocksizes(A),2),[1,size(X_in,2)+1])))
     Y = PseudoBlockArray(Y_in, BlockSizes((cumulsizes(blocksizes(A),1), [1,size(Y_in,2)+1])))
     for N=Block.(1:nblocks(X,1)), K=blockcolrange(A,N)
@@ -73,7 +73,7 @@ end
 
 function materialize!(M::MatMulMatAdd{<:AbstractColumnMajor,<:AbstractBlockBandedLayout,<:AbstractColumnMajor})
     α, A_in, X, β, Y_in = M.α, M.A, M.B, M.β, M.C
-    lmul!(β, Y_in)
+    _fill_lmul!(β, Y_in)
     A = PseudoBlockArray(A_in, BlockSizes(([1,size(A_in,1)+1],cumulsizes(blocksizes(X),1))))
     Y = PseudoBlockArray(Y_in, BlockSizes(([1,size(Y_in,1)+1],cumulsizes(blocksizes(X),2))))
     for J=Block(1):Block(nblocks(X,2)), N=blockcolrange(X,J)

test/runtests.jl

@@ -1,6 +1,7 @@
 using BlockBandedMatrices, Test, LinearAlgebra
 
 include("test_blockbanded.jl")
+include("test_blockskyline.jl")
 include("test_bandedblockbanded.jl")
 include("test_broadcasting.jl")
 include("test_linalg.jl")

test/test_blockskyline.jl

@@ -0,0 +1,29 @@
+
+using LazyArrays, BlockBandedMatrices, LinearAlgebra, Random, Test
+
+Random.seed!(0)
+
+@testset "BlockSkylineMatrix" begin
+
+    @testset "@jagot lmul! bug" begin
+        rows = rand(1:10, 5)
+        l = rand(-2:2, 5)
+        u = rand(-2:2, 5)
+
+        m = sum(rows)
+
+        A = BlockSkylineMatrix(Zeros(m,m), (rows,rows), (l,u))
+        A.data .= rand(size(A.data)...)
+
+        V = zeros(m,2)
+        V[:,1] .= rand(m)
+        reference = A*V[:,1]
+
+        @view(V[:,2]) .= Mul(A, @view(V[:,1]))
+        @test V[:,2] ≈ reference
+
+        V[:,2] .= NaN
+        @view(V[:,2]) .= Mul(A, @view(V[:,1]))
+        @test V[:,2] ≈ reference
+    end
+end