update broadcasting

Project.toml

@@ -19,9 +19,9 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-BandedMatrices = "0.11"
+BandedMatrices = "0.12"
 BlockArrays = "0.10"
 FillArrays = "0.7"
-LazyArrays = "0.11"
+LazyArrays = "0.12"
 MatrixFactorizations = "0.1,0.2"
 julia = "1"

src/BandedBlockBandedMatrix.jl

@@ -468,7 +468,10 @@ end
 ######################################
 # BandedMatrix interface  for Blocks #
 ######################################
-@inline bandwidths(V::BandedBlockBandedBlock) = subblockbandwidths(parent(V))
+@inline function bandwidths(V::BandedBlockBandedBlock) 
+    inblockbands(V) && return subblockbandwidths(parent(V))
+    (-720,-720)
+end
 
 
 
@@ -477,7 +480,8 @@ blocks(V::BandedBlockBandedBlock)::Tuple{Int,Int} = Int(first(parentindices(V)).
                                                     Int(last(parentindices(V)).block)
 
 
-function bandeddata(V::BandedBlockBandedBlock)
+function bandeddata(V::BandedBlockBandedBlock{T}) where T
+    inblockbands(V) || return Array{T}(undef, 0, size(V,2))
     A = parent(V)
     u = A.u
     K_sl, J_sl = parentindices(V)

src/BlockBandedMatrices.jl

@@ -26,7 +26,7 @@ import LazyArrays: AbstractStridedLayout, ColumnMajor, @lazymul, MatMulMatAdd, M
                     triangularlayout, UpperTriangularLayout, TriangularLayout, MatLdivVec,
                     triangulardata, subarraylayout, @lazyldiv, @lazylmul,
                     AbstractColumnMajor, DenseColumnMajor, ColumnMajor,
-                    DiagonalLayout, apply!, materialize!, MulAdd
+                    DiagonalLayout, apply!, materialize!, MulAdd, mulapplystyle, MulAddStyle
 
 import BlockArrays: BlockSizes, nblocks, blocksize, blockcheckbounds, global2blockindex,
                         Block, BlockSlice, getblock, unblock, setblock!, globalrange,

src/broadcast.jl

@@ -357,3 +357,9 @@ function similar(bc::Broadcasted{BandedBlockBandedStyle, <:Any, typeof(+),
 
     BandedBlockBandedMatrix{T}(undef, _combine_blocksizes(A,B), (max(Al,Bl), max(Au,Bu)), (max(Aλ,Bλ), max(Aμ,Bμ)))
 end
+
+
+####
+# Special case for Diagonal multiplicartion
+####
+

src/linalg.jl

@@ -165,7 +165,8 @@ end
 similar(M::MulAdd{<:DiagonalLayout,<:AbstractBlockBandedLayout}, ::Type{T}) where T = similar(M.B,T)
 similar(M::MulAdd{<:AbstractBlockBandedLayout,<:DiagonalLayout}, ::Type{T}) where T = similar(M.A,T)
 
-
+mulapplystyle(::DiagonalLayout, ::AbstractBlockBandedLayout) = MulAddStyle()
+mulapplystyle(::AbstractBlockBandedLayout, ::DiagonalLayout) = MulAddStyle()
 
 
 

src/triblockbanded.jl

@@ -98,68 +98,71 @@ end
 end
 
 
-
-@inline function materialize!(M::MatLdivVec{<:TriangularLayout{'U',UNIT,<:AbstractBlockBandedLayout},
-                                   <:AbstractStridedLayout}) where UNIT
-    U,dest = M.A,M.B
-    T = eltype(dest)
-
-    A = triangulardata(U)
-    @assert hasmatchingblocks(A)
-
-    @boundscheck size(A,1) == size(dest,1) || throw(BoundsError())
-
-    # impose block structure
-    b = PseudoBlockArray(dest, BlockSizes((cumulsizes(blocksizes(A),1),)))
-
-    Bs = blocksizes(A)
-    N = nblocks(Bs,1)
-
-    for K = N:-1:1
-        b_2 = view(b, Block(K))
-        Ũ = _triangular_matrix(Val('U'), Val(UNIT), view(A, Block(K,K)))
-        apply!(\, Ũ, b_2)
-
-        if K ≥ 2
-            KR = blockcolstart(A, K):Block(K-1)
-            V_12 = view(A, KR, Block(K))
-            b̃_1 = view(b, KR)
-            b̃_1 .=  (-one(T)).*Mul(V_12, b_2) .+ b̃_1
+for UNIT in ('U', 'N')
+    @eval begin
+        @inline function materialize!(M::MatLdivVec{<:TriangularLayout{'U',$UNIT,<:AbstractBlockBandedLayout},
+                                        <:AbstractStridedLayout})
+            U,dest = M.A,M.B
+            T = eltype(dest)
+
+            A = triangulardata(U)
+            @assert hasmatchingblocks(A)
+
+            @boundscheck size(A,1) == size(dest,1) || throw(BoundsError())
+
+            # impose block structure
+            b = PseudoBlockArray(dest, BlockSizes((cumulsizes(blocksizes(A),1),)))
+
+            Bs = blocksizes(A)
+            N = nblocks(Bs,1)
+
+            for K = N:-1:1
+                b_2 = view(b, Block(K))
+                Ũ = _triangular_matrix(Val('U'), Val($UNIT), view(A, Block(K,K)))
+                apply!(\, Ũ, b_2)
+
+                if K ≥ 2
+                    KR = blockcolstart(A, K):Block(K-1)
+                    V_12 = view(A, KR, Block(K))
+                    b̃_1 = view(b, KR)
+                    materialize!(MulAdd(-one(T), V_12, b_2, one(T), b̃_1))
+                end
+            end
+
+            dest
         end
-    end
 
-    dest
-end
+        @inline function materialize!(M::MatLdivVec{<:TriangularLayout{'L',$UNIT,<:AbstractBlockBandedLayout},
+                                        <:AbstractStridedLayout})
+            L,dest = M.A, M.B
+            T = eltype(dest)
+            A = triangulardata(L)
+            @assert hasmatchingblocks(A)
 
-@inline function materialize!(M::MatLdivVec{<:TriangularLayout{'L',UNIT,<:AbstractBlockBandedLayout},
-                                   <:AbstractStridedLayout}) where UNIT
-    L,dest = M.A, M.B
-    T = eltype(dest)
-    A = triangulardata(L)
-    @assert hasmatchingblocks(A)
+            @boundscheck size(A,1) == size(dest,1) || throw(BoundsError())
 
-    @boundscheck size(A,1) == size(dest,1) || throw(BoundsError())
+            # impose block structure
+            b = PseudoBlockArray(dest, BlockSizes((cumulsizes(blocksizes(A),1),)))
 
-    # impose block structure
-    b = PseudoBlockArray(dest, BlockSizes((cumulsizes(blocksizes(A),1),)))
+            Bs = blocksizes(A)
+            N = nblocks(Bs,1)
 
-    Bs = blocksizes(A)
-    N = nblocks(Bs,1)
+            for K = 1:N
+                b_2 = view(b, Block(K))
+                L̃ = _triangular_matrix(Val('L'), Val($UNIT), view(A, Block(K,K)))
+                b_2 .= Ldiv(L̃, b_2)
 
-    for K = 1:N
-       b_2 = view(b, Block(K))
-       L̃ = _triangular_matrix(Val('L'), Val(UNIT), view(A, Block(K,K)))
-       b_2 .= Ldiv(L̃, b_2)
-
-       if K < N
-           KR = Block(K+1):blockcolstop(A, K)
-           V_12 = view(A, KR, Block(K))
-           b̃_1 = view(b, KR)
-           b̃_1 .=  (-one(T)).*Mul(V_12, b_2) .+ b̃_1
-       end
-    end
+                if K < N
+                    KR = Block(K+1):blockcolstop(A, K)
+                    V_12 = view(A, KR, Block(K))
+                    b̃_1 = view(b, KR)
+                    materialize!(MulAdd(-one(T), V_12, b_2, one(T), b̃_1))
+                end
+            end
 
-    dest
+            dest
+        end
+    end
 end
 
 

test/test_bandedblockbanded.jl

@@ -1,5 +1,5 @@
 using BlockArrays, BandedMatrices, BlockBandedMatrices, FillArrays, SparseArrays, Test, LazyArrays , LinearAlgebra
-import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolrange, blockrowrange, colrange, rowrange, isbandedblockbanded
+import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolrange, blockrowrange, colrange, rowrange, isbandedblockbanded, bandeddata
 
 @testset "BandedBlockBandedMatrix" begin
     @testset "constructors" begin
@@ -190,7 +190,11 @@ import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolrange, blockrowran
         @test A[1,2] == view(A,Block(1,2))[1,1] == 11
         @test A[1,3] == view(A,Block(1,2))[1,2] == view(A,Block(1,2))[2] == 19
 
-        @test view(A, Block(3),Block(1)) ≈ [0,0,0]
+        @test bandwidths(view(A, Block(3),Block(1))) == (-720,-720)
+        @test isempty(bandeddata(view(A, Block(3),Block(1))))
+
+        @test A[Block(3,1)] == view(A, Block(3),Block(1)) == zeros(3,1)
+        @test A[Block(3,1)] ≈ view(A, Block(3),Block(1)) ≈ zeros(3,1)
         @test_throws BandError view(A, Block(3),Block(1))[1,1] = 4
         @test_throws BoundsError view(A, Block(5,1))
     end