Fast copyto! for banded blocks (#73)

* Fast copyto! for banded blocks

* Fix bounds check

* update project

* Update Project.toml

* Update .travis.yml

Author: dlfivefifty

.editorconfig

@@ -13,9 +13,9 @@ matrix:
 notifications:
   email: false
 after_success:
-  - julia -e 'Pkg.add("Documenter")'
-  - julia -e 'cd(Pkg.dir("BlockBandedMatrices")); include(joinpath("docs", "make.jl"))'
-  - julia -e 'cd(Pkg.dir("BlockBandedMatrices")); Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())'
+  - julia -e 'using Pkg; Pkg.add("Documenter")'
+  - julia -e 'using Pkg; cd(Pkg.dir("BlockBandedMatrices")); include(joinpath("docs", "make.jl"))'
+  - julia -e 'using Pkg; Pkg.add("Coverage"); using Coverage; Codecov.submit(Codecov.process_folder())'
 
 # uncomment the following lines to override the default test script
 #script:

.travis.yml

@@ -1,6 +1,6 @@
 name = "BlockBandedMatrices"
 uuid = "ffab5731-97b5-5995-9138-79e8c1846df0"
-version = "0.8.4"
+version = "0.8.5"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -18,9 +18,9 @@ Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-ArrayLayouts = "0.2.6, 0.3"
-BandedMatrices = "0.15.6"
-BlockArrays = "0.12.4"
-FillArrays = "0.8.8"
-MatrixFactorizations = "0.4"
+ArrayLayouts = "0.3.3"
+BandedMatrices = "0.15.10"
+BlockArrays = "0.12.7"
+FillArrays = "0.8.10"
+MatrixFactorizations = "0.4.1"
 julia = "1.2"

Project.toml

@@ -128,7 +128,7 @@ end
 #  RaggedMatrix interface
 ######################################
 
-@inline function colstart(A::AbstractBlockBandedMatrix, i::Integer)
+@inline function blockbanded_colstart(A, i::Integer)
     bs = blockcolstart(A,findblock(axes(A,2),i))
     if isempty(axes(A,1)) 
         1
@@ -139,26 +139,26 @@ end
     end
 end
 
-@inline function colstop(A::AbstractBlockBandedMatrix, i::Integer)
+@inline function blockbanded_colstop(A, i::Integer)
     CS = blockcolstop(A,findblock(axes(A,2),i))
     CS == Block(0) && return 0
     last(axes(A,1)[CS])
 end
 
-@inline rowstart(A::AbstractBlockBandedMatrix, i::Integer) =
+@inline blockbanded_rowstart(A, i::Integer) =
     first(axes(A,2)[blockrowstart(A,findblock(axes(A,1),i))])
 
-@inline function rowstop(A::AbstractBlockBandedMatrix, i::Integer)
+@inline function blockbanded_rowstop(A, i::Integer)
     CS = blockrowstop(A,findblock(axes(A,1),i))
     CS == Block(0) && return 0
     last(axes(A,2)[CS])    
 end
 
-@inline blockbanded_colsupport(A, j::Integer) = colrange(A, j)
-@inline blockbanded_rowsupport(A, j::Integer) = rowrange(A, j)
+@inline blockbanded_colsupport(A, j::Integer) = blockbanded_colstart(A, j):blockbanded_colstop(A, j)
+@inline blockbanded_rowsupport(A, j::Integer) = blockbanded_rowstart(A, j):blockbanded_rowstop(A, j)
 
-@inline blockbanded_rowsupport(A, j) = isempty(j) ? (1:0) : rowstart(A,minimum(j)):rowstop(A,maximum(j))
-@inline blockbanded_colsupport(A, j) = isempty(j) ? (1:0) : colstart(A,minimum(j)):colstop(A,maximum(j))
+@inline blockbanded_rowsupport(A, j) = isempty(j) ? (1:0) : blockbanded_rowstart(A,minimum(j)):blockbanded_rowstop(A,maximum(j))
+@inline blockbanded_colsupport(A, j) = isempty(j) ? (1:0) : blockbanded_colstart(A,minimum(j)):blockbanded_colstop(A,maximum(j))
 
 @inline colsupport(::AbstractBlockBandedLayout, A, j) = blockbanded_colsupport(A, j)
 @inline rowsupport(::AbstractBlockBandedLayout, A, j) = blockbanded_rowsupport(A, j)

src/AbstractBlockBandedMatrix.jl

@@ -318,7 +318,7 @@ end
 end
 
 ## structured matrix methods ##
-function _bandedblockbanded_replace_in_print_matrix(A, i, j, s)
+function layout_replace_in_print_matrix(::AbstractBandedBlockBandedLayout, A, i, j, s)
     bi = findblockindex.(axes(A), (i,j))
     I,J = block.(bi)
     i,j = blockindex.(bi)
@@ -327,9 +327,6 @@ function _bandedblockbanded_replace_in_print_matrix(A, i, j, s)
     -l ≤ Int(J-I) ≤ u && -λ ≤ j-i ≤ μ ? s : Base.replace_with_centered_mark(s)
 end
 
-Base.replace_in_print_matrix(A::BandedBlockBandedMatrix, i::Integer, j::Integer, s::AbstractString) = 
-    _bandedblockbanded_replace_in_print_matrix(A, i, j, s)
-
 
 ############
 # Indexing #

src/BandedBlockBandedMatrix.jl

@@ -7,7 +7,7 @@ import Base: getindex, setindex!, checkbounds, @propagate_inbounds, convert,
                         unsafe_convert, fill!, length, first, last,
                         eltype, getindex, to_indices, to_index,
                         reindex, _maybetail, tail, @_propagate_inbounds_meta,
-                        ==, axes, copyto!, similar, OneTo, replace_in_print_matrix
+                        ==, axes, copyto!, similar, OneTo
 
 import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, DefaultArrayStyle, Broadcasted, broadcasted,
                         materialize, materialize!
@@ -26,19 +26,19 @@ import ArrayLayouts: BlasMatLmulVec,
                     triangulardata, sublayout, sub_materialize,
                     AbstractColumnMajor, DenseColumnMajor, ColumnMajor,
                     DiagonalLayout, MulAdd, mul, colsupport, rowsupport,
-                    _qr, _factorize, _copyto!, zero!
+                    _qr, _factorize, _copyto!, zero!, layout_replace_in_print_matrix
 
 import BlockArrays: blocksize, blockcheckbounds, BlockedUnitRange, blockisequal, DefaultBlockAxis,
                         Block, BlockSlice, getblock, unblock, setblock!, block, blockindex,
                         _blocklengths2blocklasts, BlockIndexRange, sizes_from_blocks, BlockSlice1,
                         blockcolsupport, blockrowsupport, blockcolstart, blockcolstop, blockrowstart, blockrowstop,
-                        AbstractBlockLayout
+                        AbstractBlockLayout, BlockLayout, blocks
 
 import BandedMatrices: isbanded, bandwidths, bandwidth, banded_getindex, colrange,
                         inbands_setindex!, inbands_getindex, banded_setindex!,
                         banded_generic_axpy!,
                         BlasFloat, banded_dense_axpy!, MemoryLayout,
-                        BandedLayout, BandedColumnMajor,
+                        BandedLayout, BandedColumnMajor, BandedColumns,
                         BandedSubBandedMatrix, bandeddata, tribandeddata,
                         _BandedMatrix, colstart, colstop, rowstart, rowstop,
                         BandedStyle, _fill_lmul!, bandshift,

src/BlockBandedMatrices.jl

@@ -26,8 +26,6 @@ function bb_blockstarts(axes, l::AbstractVector{Int}, u::AbstractVector{Int})
     b_start
 end
 
-bb_blockstarts(b_size, l::Integer, u::Integer) = bb_blockstarts(b_size, Fill(l, nblocks(b_size,2)), Fill(u, nblocks(b_size,2)))
-
 function bb_blockstrides(b_axes, l::AbstractVector{Int}, u::AbstractVector{Int})
     N, M = blocksize.(b_axes,1)
     L,U = maximum(l), maximum(u)
@@ -233,18 +231,12 @@ BlockBandedMatrix(A::Union{AbstractMatrix,UniformScaling},
 BlockBandedMatrix(A::AbstractMatrix, lu::NTuple{2,Int}) = BlockBandedMatrix(A, BlockBandedSizes(axes(A), lu...))
 
 function convert(::Type{BlockSkylineMatrix}, A::AbstractMatrix)
-    @assert isblockbanded(A)
     block_sizes = BlockSkylineSizes(axes(A), colblockbandwidths(A)...)
 
-    ret = BlockSkylineMatrix{eltype(A)}(undef, block_sizes)
-    for J = blockaxes(ret,2), K = blockcolsupport(ret, J)
-        view(ret, K, J) .= view(A, K, J)
-    end
-    ret
+    copyto!(BlockSkylineMatrix{eltype(A)}(undef, block_sizes), A)
 end
 
 function convert(::Type{BlockBandedMatrix}, A::AbstractMatrix)
-    @assert isblockbanded(A)
     convert(BlockSkylineMatrix, A)
 end
 
@@ -320,13 +312,6 @@ end
     return v
 end
 
-## structured matrix methods ##
-function Base.replace_in_print_matrix(A::BlockSkylineMatrix, i::Integer, j::Integer, s::AbstractString)
-    bi = findblockindex.(axes(A), (i,j))
-    I,J = Int.(block.(bi))
-    -A.block_sizes.l[J] ≤ J-I ≤ A.block_sizes.u[J] ? s : Base.replace_with_centered_mark(s)
-end
-
 ############
 # Indexing #
 ############
@@ -405,8 +390,8 @@ const BlockBandedBlock{T} = SubArray{T,2,<:BlockSkylineMatrix,<:Tuple{<:BlockSli
 
 
 # gives the columns of parent(V).data that encode the block
-blocks(V::BlockBandedBlock)::Tuple{Int,Int} = first(first(parentindices(V)).block.n),
-first(last(parentindices(V)).block.n)
+_parent_blocks(V::BlockBandedBlock)::Tuple{Int,Int} = 
+    first(first(parentindices(V)).block.n),first(last(parentindices(V)).block.n)
 
 ######################################
 # Matrix interface  for Blocks #
@@ -417,16 +402,16 @@ MemoryLayout(::Type{<:BlockBandedBlock}) = ColumnMajor()
 
 function Base.unsafe_convert(::Type{Ptr{T}}, V::BlockBandedBlock{T}) where T
     A = parent(V)
-    K,J = blocks(V)
+    K,J = _parent_blocks(V)
     Base.unsafe_convert(Ptr{T}, A.data) + sizeof(T)*(blockstart(A,K,J)-1)
 end
 
-strides(V::BlockBandedBlock) = (1,parent(V).block_sizes.block_strides[blocks(V)[2]])
+strides(V::BlockBandedBlock) = (1,parent(V).block_sizes.block_strides[_parent_blocks(V)[2]])
 
 @propagate_inbounds function getindex(V::BlockBandedBlock, k::Int, j::Int)
     @boundscheck checkbounds(V, k, j)
     A = parent(V)
-    K,J = blocks(V)
+    K,J = _parent_blocks(V)
     if -A.block_sizes.l[J] ≤ J-K ≤ A.block_sizes.u[J]
         b_start = blockstart(A,K,J)
         b_start == 0 && return zero(eltype(V))
@@ -440,7 +425,7 @@ end
 @propagate_inbounds function setindex!(V::BlockBandedBlock, v, k::Int, j::Int)
     @boundscheck checkbounds(V, k, j)
     A = parent(V)
-    K,J = blocks(V)
+    K,J = _parent_blocks(V)
     if -A.block_sizes.l[J] ≤ J-K ≤ A.block_sizes.u[J]
         b_start = blockstart(A,K,J)
         # TODO: What to do if b_start == 0 ?

src/BlockSkylineMatrix.jl

@@ -80,6 +80,62 @@ end
 
 
 _copyto!(_, ::AbstractBlockBandedLayout, dest::AbstractMatrix, src::AbstractMatrix) = blockbanded_copyto!(dest, src)
+function _copyto!(_, ::BlockLayout{<:BandedColumns}, dest::AbstractMatrix, src::AbstractMatrix)
+    if !blockisequal(axes(dest), axes(src))
+        copyto!(PseudoBlockArray(dest, axes(src)), src)
+        return dest
+    end
+
+    srcB = blocks(src)
+    srcD = bandeddata(srcB)
+
+    dl, du = colblockbandwidths(dest)
+    sl, su = bandwidths(srcB)
+    M,N = size(srcB)
+    # Source matrix must fit within bands of destination matrix
+    all(dl .≥ min(sl,M-1)) && all(du .≥ min(su,N-1)) || throw(BandError(dest))
+
+    for J = 1:N
+        for K = max(1,J-du[J]):min(J-su-1,M)
+            zero!(view(dest,Block(K),Block(J)))
+        end
+        for K = max(1,J-su):min(J+sl,M)
+            copyto!(view(dest,Block(K),Block(J)), srcD[K-J+su+1,J])
+        end
+        for K = max(1,J+sl+1):min(J+dl[J],M)
+            zero!(view(dest,Block(K),Block(J)))
+        end
+    end
+    dest
+end
+
+function _copyto!(_, ::BlockLayout{<:AbstractBandedLayout}, dest::AbstractMatrix, src::AbstractMatrix)
+    if !blockisequal(axes(dest), axes(src))
+        copyto!(PseudoBlockArray(dest, axes(src)), src)
+        return dest
+    end
+
+    srcB = blocks(src)
+
+    dl, du = colblockbandwidths(dest)
+    sl, su = bandwidths(srcB)
+    M,N = size(srcB)
+    # Source matrix must fit within bands of destination matrix
+    all(dl .≥ min(sl,M-1)) && all(du .≥ min(su,N-1)) || throw(BandError(dest))
+
+    for J = 1:N
+        for K = max(1,J-du[J]):min(J-su-1,M)
+            zero!(view(dest,Block(K),Block(J)))
+        end
+        for K = max(1,J-su):min(J+sl,M)
+            copyto!(view(dest,Block(K),Block(J)), inbands_getindex(srcB, K, J))
+        end
+        for K = max(1,J+sl+1):min(J+dl[J],M)
+            zero!(view(dest,Block(K),Block(J)))
+        end
+    end
+    dest
+end
 
 
 function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(identity)})

src/broadcast.jl

@@ -73,18 +73,4 @@ for op in (:-, :+)
             mortar(Tridiagonal(A.blocks.dl, broadcast($op, Ref(λ), A.blocks.d), A.blocks.du))
         end
     end
-end
-
-function replace_in_print_matrix(A::BlockDiagonal, i::Integer, j::Integer, s::AbstractString)
-    bi = findblockindex.(axes(A), (i,j))
-    I,J = block.(bi)
-    i,j = blockindex.(bi)
-    Int(J-I) == 0 ? s : Base.replace_with_centered_mark(s)
-end
-
-function replace_in_print_matrix(A::BlockTridiagonal, i::Integer, j::Integer, s::AbstractString)
-    bi = findblockindex.(axes(A), (i,j))
-    I,J = block.(bi)
-    i,j = blockindex.(bi)
-    -1 ≤ Int(J-I) ≤ 1 ? s : Base.replace_with_centered_mark(s)
 end

src/interfaceimpl.jl

@@ -97,6 +97,7 @@ isblockbanded(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{BlockRange1}, <:Blo
     isblockbanded(parent(V))
 
 sub_materialize(::AbstractBlockBandedLayout, V, _) = BlockBandedMatrix(V)
+sub_materialize(::BlockLayout{<:AbstractBandedLayout}, V, _) = BlockBandedMatrix(V)
 
 strides(V::SubBlockSkylineMatrix{<:Any,LL,UU,<:Union{BlockRange1,Block1},Block1}) where {LL,UU} =
     (1,parent(V).block_sizes.block_strides[Int(parentindices(V)[2].block)])