Skip to content

Commit 79eac5e

Browse files
dlfivefiftydlfivefifty
authored
Support non-UnitRange based ranges (#94)
* Support non-UnitRange based ranges * BlockArrays new view * blockbandwidths and subblockbandwidths for Zeros * Add ==, sublayout * tests pass * v0.10 * Update Project.toml * simplify sublayout * add simple implementation for tests
1 parent 223cb15 commit 79eac5e

13 files changed

+245
-245
lines changed

‎Project.toml‎

Lines changed: 7 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,14 @@
11
name = "BlockBandedMatrices"
22
uuid = "ffab5731-97b5-5995-9138-79e8c1846df0"
3-
version = "0.9.6"
3+
version = "0.10.0"
44

55
[deps]
66
ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
77
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
88
BlockArrays = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
99
Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
1010
FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
11+
LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
1112
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
1213
MatrixFactorizations = "a3b82374-2e81-5b9e-98ce-41277c0e4c87"
1314
Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
@@ -16,11 +17,11 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
1617
Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
1718

1819
[compat]
19-
ArrayLayouts = "0.4.1"
20-
BandedMatrices = "0.15.16"
21-
BlockArrays = "0.13"
22-
FillArrays = "0.9.2, 0.10"
23-
MatrixFactorizations = "0.5.2, 0.6, 0.7, 0.8"
20+
ArrayLayouts = "0.5"
21+
BandedMatrices = "0.16"
22+
BlockArrays = "0.14"
23+
FillArrays = "0.11"
24+
MatrixFactorizations = "0.7.1, 0.8"
2425
julia = "1.5"
2526

2627
[extras]

‎src/AbstractBlockBandedMatrix.jl‎

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -172,3 +172,15 @@ function fill!(A::AbstractBlockBandedMatrix, val::Any)
172172
fill!(A.data, val)
173173
A
174174
end
175+
176+
177+
function ==(A::AbstractBlockBandedMatrix, B::AbstractBlockBandedMatrix)
178+
axes(A) == axes(B) || return false
179+
blockisequal(axes(A), axes(B)) || return Base.invoke(==, Tuple{AbstractArray,AbstractArray}, A, B)
180+
l,u = max.(blockbandwidths(A), blockbandwidths(B))
181+
N = blocksize(A,1)
182+
for J = blockaxes(A,2), K = max(Block(1),J-u):min(J+l,Block(N))
183+
view(A, K, J) == view(B, K, J) || return false
184+
end
185+
return true
186+
end

‎src/BandedBlockBandedMatrix.jl‎

Lines changed: 61 additions & 98 deletions
Large diffs are not rendered by default.

‎src/BlockBandedMatrices.jl‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -29,7 +29,7 @@ import ArrayLayouts: BlasMatLmulVec, MatLmulVec, MatLmulMat,
2929
_qr, _factorize, _copyto!, zero!, layout_replace_in_print_matrix
3030

3131
import BlockArrays: blocksize, blockcheckbounds, BlockedUnitRange, blockisequal, DefaultBlockAxis,
32-
Block, BlockSlice, getblock, unblock, setblock!, block, blockindex,
32+
Block, BlockSlice, unblock, block, blockindex,
3333
_blocklengths2blocklasts, BlockIndexRange, sizes_from_blocks, BlockSlice1,
3434
blockcolsupport, blockrowsupport, blockcolstart, blockcolstop, blockrowstart, blockrowstop,
3535
AbstractBlockLayout, BlockLayout, blocks, hasmatchingblocks, BlockStyle
@@ -38,7 +38,7 @@ import BandedMatrices: isbanded, bandwidths, bandwidth, banded_getindex, colrang
3838
inbands_setindex!, inbands_getindex, banded_setindex!,
3939
banded_generic_axpy!,
4040
BlasFloat, banded_dense_axpy!, MemoryLayout,
41-
BandedLayout, BandedColumnMajor, BandedColumns,
41+
BandedLayout, BandedColumnMajor, BandedColumns, bandedcolumns,
4242
BandedSubBandedMatrix, bandeddata,
4343
_BandedMatrix, colstart, colstop, rowstart, rowstop,
4444
BandedStyle, _fill_lmul!, bandshift,
@@ -49,8 +49,8 @@ export BandedBlockBandedMatrix, BlockBandedMatrix, BlockSkylineMatrix, blockband
4949

5050

5151
const Block1 = Block{1,Int}
52-
const BlockRange1 = BlockRange{1,Tuple{UnitRange{Int}}}
53-
const BlockIndexRange1 = BlockIndexRange{1,Tuple{UnitRange{Int}}}
52+
const BlockRange1{R<:AbstractUnitRange{Int}} = BlockRange{1,Tuple{R}}
53+
const BlockIndexRange1{R<:AbstractUnitRange{Int}} = BlockIndexRange{1,Tuple{R}}
5454

5555
blockcolrange(A...) = blockcolsupport(A...)
5656
blockrowrange(A...) = blockrowsupport(A...)

‎src/BlockSkylineMatrix.jl‎

Lines changed: 1 addition & 59 deletions
Original file line numberDiff line numberDiff line change
@@ -300,23 +300,6 @@ BroadcastStyle(::Type{<:BlockBandedMatrix}) = BlockBandedStyle()
300300
zeroblock(A::BlockSkylineMatrix, K::Int, J::Int) =
301301
Matrix(Zeros{eltype(A)}(length.(getindex.(axes(A),(Block(K),Block(J))))))
302302

303-
@inline function getblock(A::BlockSkylineMatrix, K::Int, J::Int)
304-
@boundscheck blockcheckbounds(A, K, J)
305-
if -A.block_sizes.l[J] ≤ J - K ≤ A.block_sizes.u[J]
306-
convert(Matrix, view(A, Block(K, J)))
307-
else
308-
zeroblock(A, K, J)
309-
end
310-
end
311-
312-
# @inline function Base.getindex(block_arr::BlockArray{T,N}, blockindex::BlockIndex{N}) where {T,N}
313-
# @boundscheck checkbounds(block_arr.blocks, blockindex.I...)
314-
# @inbounds block = block_arr.blocks[blockindex.I...]
315-
# @boundscheck checkbounds(block, blockindex.α...)
316-
# @inbounds v = block[blockindex.α...]
317-
# return v
318-
# end
319-
320303

321304
###########################
322305
# AbstractArray Interface #
@@ -366,48 +349,7 @@ end
366349
return A
367350
end
368351

369-
# @propagate_inbounds function Base.setindex!(block_array::BlockArray{T, N}, v, block_index::BlockIndex{N}) where {T,N}
370-
# getblock(block_array, block_index.I...)[block_index.α...] = v
371-
# end
372352

373-
########
374-
# Misc #
375-
########
376-
377-
# @generated function Base.Array(block_array::BlockArray{T, N, R}) where {T,N,R}
378-
# # TODO: This will fail for empty block array
379-
# return quote
380-
# block_sizes = block_array.block_sizes
381-
# arr = similar(block_array.blocks[1], size(block_array)...)
382-
# @nloops $N i i->(1:nblocks(block_sizes, i)) begin
383-
# block_index = @ntuple $N i
384-
# indices = globalrange(block_sizes, block_index)
385-
# arr[indices...] = getblock(block_array, block_index...)
386-
# end
387-
#
388-
# return arr
389-
# end
390-
# end
391-
#
392-
# @generated function Base.copyto!(block_array::BlockArray{T, N, R}, arr::R) where {T,N,R <: AbstractArray}
393-
# return quote
394-
# block_sizes = block_array.block_sizes
395-
#
396-
# @nloops $N i i->(1:nblocks(block_sizes, i)) begin
397-
# block_index = @ntuple $N i
398-
# indices = globalrange(block_sizes, block_index)
399-
# copyto!(getblock(block_array, block_index...), arr[indices...])
400-
# end
401-
#
402-
# return block_array
403-
# end
404-
# end
405-
#
406-
# function Base.fill!(block_array::BlockArray, v)
407-
# for block in block_array.blocks
408-
# fill!(block, v)
409-
# end
410-
# end
411353

412354

413355
##################
@@ -417,7 +359,7 @@ end
417359
# with StridedMatrix.
418360
##################
419361

420-
const BlockBandedBlock{T} = SubArray{T,2,<:BlockSkylineMatrix,<:Tuple{<:BlockSlice1,<:BlockSlice1},false}
362+
const BlockBandedBlock{T} = SubArray{T,2,<:BlockSkylineMatrix,<:Tuple{BlockSlice1,BlockSlice1},false}
421363

422364

423365

‎src/broadcast.jl‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -214,6 +214,15 @@ function _copyto!(_, ::BlockLayout{<:AbstractBandedLayout}, dest::AbstractMatrix
214214
dest
215215
end
216216

217+
function _copyto!(::BandedBlockBandedColumns, ::BandedBlockBandedColumns, dest::AbstractMatrix, src::AbstractMatrix)
218+
if blockbandwidths(dest) == blockbandwidths(src) && subblockbandwidths(dest) == subblockbandwidths(src)
219+
copyto!(bandedblockbandeddata(dest), bandedblockbandeddata(src))
220+
dest
221+
else # TODO: if subblockbandwidths match we can still do something better
222+
_copyto!(BandedBlockBandedLayout(), BandedBlockBandedLayout(), dest, src)
223+
end
224+
end
225+
217226

218227
function copyto!(dest::AbstractArray, bc::Broadcasted{<:AbstractBlockBandedStyle, <:Any, typeof(identity)})
219228
(A,) = bc.args

‎src/interfaceimpl.jl‎

Lines changed: 12 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -66,14 +66,16 @@ end
6666
blockbandwidths(A::BlockArray) = bandwidths(A.blocks)
6767
isblockbanded(A::BlockArray) = isbanded(A.blocks)
6868

69-
@inline function getblock(block_arr::BlockBidiagonal{T,VT}, K::Int, J::Int) where {T,VT<:AbstractMatrix}
69+
@inline function BlockArrays.viewblock(block_arr::BlockBidiagonal{T,VT}, KJ::Block{2}) where {T,VT<:AbstractMatrix}
70+
K,J = KJ.n
7071
@boundscheck blockcheckbounds(block_arr, K, J)
7172
l,u = blockbandwidths(block_arr)
7273
-l ≤ (J-K) ≤ u || return convert(VT, Zeros{T}(length.(getindex.(axes(block_arr),(Block(K),Block(J))))...))
7374
block_arr.blocks[K,J]
7475
end
7576

76-
@inline function getblock(block_arr::BlockTridiagonal{T,VT}, K::Int, J::Int) where {T,VT<:AbstractMatrix}
77+
@inline function BlockArrays.viewblock(block_arr::BlockTridiagonal{T,VT}, KJ::Block{2}) where {T,VT<:AbstractMatrix}
78+
K,J = KJ.n
7779
@boundscheck blockcheckbounds(block_arr, K, J)
7880
abs(J-K) ≥ 2 && return convert(VT, Zeros{T}(length.(getindex.(axes(block_arr),(Block(K),Block(J))))...))
7981
block_arr.blocks[K,J]
@@ -114,3 +116,11 @@ for op in (:-, :+)
114116
end
115117
end
116118
end
119+
120+
121+
###
122+
# Zeros
123+
###
124+
125+
blockbandwidths(::Zeros) = (-1,-1)
126+
subblockbandwidths(::Zeros) = (-1,-1)

‎src/linalg.jl‎

Lines changed: 10 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
# BlockBandedMatrix with block range indexes is also block-banded
22
const SubBlockSkylineMatrix{T,LL,UU,R1,R2} =
3-
SubArray{T,2,BlockSkylineMatrix{T,LL,UU},<:Tuple{<:BlockSlice{R1},<:BlockSlice{R2}}}
3+
SubArray{T,2,BlockSkylineMatrix{T,LL,UU},<:Tuple{BlockSlice{<:R1},BlockSlice{<:R2}}}
44

55

66

@@ -71,7 +71,7 @@ similar(M::MulAdd{<:AbstractBlockBandedLayout,<:DiagonalLayout}, ::Type{T}) wher
7171

7272

7373

74-
function blockbandwidths(V::SubBlockSkylineMatrix{<:Any,LL,UU,BlockRange1,BlockRange1}) where {LL,UU}
74+
function blockbandwidths(V::SubBlockSkylineMatrix{<:Any,LL,UU,<:BlockRange1,<:BlockRange1}) where {LL,UU}
7575
A = parent(V)
7676

7777
KR = parentindices(V)[1].block.indices[1]
@@ -87,13 +87,14 @@ end
8787
# BlockIndexRange subblocks
8888
####
8989

90+
sublayout(::AbstractBlockBandedLayout, ::Type{<:Tuple{<:BlockSlice{<:BlockRange1}, <:BlockSlice{<:BlockRange1}}}) = BlockBandedLayout()
9091

91-
sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{BlockRange1}}}) = BlockBandedColumnMajor()
92-
sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{Block1}}}) = ColumnMajor()
93-
sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{BlockIndexRange1}}}) = ColumnMajor()
94-
sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockIndexRange1}, <:BlockSlice{BlockIndexRange1}}}) = ColumnMajor()
92+
sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{<:BlockRange1}, <:BlockSlice{<:BlockRange1}}}) = BlockBandedColumnMajor()
93+
sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{<:BlockRange1}, <:BlockSlice{Block1}}}) = ColumnMajor()
94+
sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{<:BlockRange1}, <:BlockSlice{<:BlockIndexRange1}}}) = ColumnMajor()
95+
sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{<:BlockIndexRange1}, <:BlockSlice{<:BlockIndexRange1}}}) = ColumnMajor()
9596

96-
isblockbanded(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{BlockRange1}}}) =
97+
isblockbanded(V::SubArray{<:Any,2,<:Any,<:Tuple{<:BlockSlice{<:BlockRange1}, <:BlockSlice{<:BlockRange1}}}) =
9798
isblockbanded(parent(V))
9899

99100
sub_materialize(::AbstractBlockBandedLayout, V, _) = BlockBandedMatrix(V)
@@ -113,10 +114,10 @@ function unsafe_convert(::Type{Ptr{T}}, V::SubBlockSkylineMatrix{T,LL,UU,<:Union
113114
p = unsafe_convert(Ptr{T}, view(A, first(KR), J))
114115
end
115116

116-
strides(V::SubBlockSkylineMatrix{<:Any,LL,UU,BlockRange1,BlockIndexRange1}) where {LL,UU} =
117+
strides(V::SubBlockSkylineMatrix{<:Any,LL,UU,<:BlockRange1,<:BlockIndexRange1}) where {LL,UU} =
117118
(1,parent(V).block_sizes.block_strides[Int(Block(parentindices(V)[2]))])
118119

119-
function unsafe_convert(::Type{Ptr{T}}, V::SubBlockSkylineMatrix{T,LL,UU,BlockRange1,BlockIndexRange1}) where {T,LL,UU}
120+
function unsafe_convert(::Type{Ptr{T}}, V::SubBlockSkylineMatrix{T,LL,UU,<:BlockRange1,<:BlockIndexRange1}) where {T,LL,UU}
120121
A = parent(V)
121122
JR = parentindices(V)[2]
122123
K = first(parentindices(V)[1].block)

‎test/test_bandedblockbanded.jl‎

Lines changed: 22 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,7 @@
11
using BlockArrays, BandedMatrices, BlockBandedMatrices, FillArrays, SparseArrays, Test, ArrayLayouts , LinearAlgebra
22
import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolsupport, blockrowsupport, colsupport, rowsupport,
33
isbandedblockbanded, bandeddata, BandedBlockBandedColumns
4+
import ArrayLayouts: RangeCumsum
45

56
@testset "BandedBlockBandedMatrix" begin
67
@testset "constructors" begin
@@ -181,7 +182,7 @@ import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolsupport, blockrows
181182
A = _BandedBlockBandedMatrix(data, rows,cols, (l,u), (λ,μ))
182183

183184
@test A[Block(1), Block(1)] isa BandedMatrix
184-
@test A[Block(1), Block(1)] == A[Block(1,1)] == BlockArrays.getblock(A, 1, 1) == BandedMatrix(view(A, Block(1,1)))
185+
@test A[Block(1), Block(1)] == A[Block(1,1)] == view(A, Block(1, 1)) == BandedMatrix(view(A, Block(1,1)))
185186
@test A[1,1] == view(A,Block(1),Block(1))[1,1] == view(A,Block(1,1))[1,1] == A[Block(1,1)][1,1] == A[Block(1),Block(1)][1,1] == 5
186187
@test A[2,1] == view(A,Block(2),Block(1))[1,1] == view(A,Block(2,1))[1,1] == 8
187188
@test A[3,1] == view(A,Block(2),Block(1))[2,1] == 9
@@ -454,20 +455,20 @@ import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolsupport, blockrows
454455
A = BandedBlockBandedMatrix{Float64}(undef, 1:2, 1:5,(-1,1), (-1,1))
455456
A.data .= randn.()
456457
V = view(A, Block(2,3))
457-
@test MemoryLayout(typeof(V)) == BandedMatrices.BandedColumnMajor()
458+
@test MemoryLayout(typeof(V)) == BandedMatrices.BandedColumns{ColumnMajor}()
458459
@test isbanded(V)
459460
@test bandwidths(V) == (-1,1)
460461
@test BandedMatrix(V) == A[2:3,4:6]
461462
V2 = view(V, :, 2:3)
462-
@test MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumnMajor()
463+
@test MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumns{ColumnMajor}()
463464
@test bandwidths(V2) == (0,0)
464465
@test BandedMatrix(V2) == A[2:3,5:6]
465466
V2 = view(V, 2:2, :)
466-
@test MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumnMajor()
467+
@test MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumns{ColumnMajor}()
467468
@test bandwidths(V2) == (-2,2)
468469
@test BandedMatrix(V2) == A[3:3,4:6]
469470
V2 = view(V, 2:2, 2:3)
470-
@test MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumnMajor()
471+
@test MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumns{ColumnMajor}()
471472
@test bandwidths(V2) == (-1,1)
472473
@test BandedMatrix(V2) == A[3:3,5:6]
473474
end
@@ -487,6 +488,22 @@ import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolsupport, blockrows
487488
A = _BandedBlockBandedMatrix(PseudoBlockVector([1,2,3],[1,2])', blockedrange([1,2]), (-1,1), (-1,1))
488489
@test MemoryLayout(A) isa BandedBlockBandedColumns{RowMajor}
489490
end
491+
492+
@testset "1:N blocks" begin
493+
N = 10
494+
A = BandedBlockBandedMatrix{Float64}(undef, 1:N,1:N, (1,1), (1,1))
495+
@test axes(A) isa NTuple{2,BlockedUnitRange{<:RangeCumsum}}
496+
end
497+
498+
@testset "change bandwidths" begin
499+
l , u = 1,1
500+
λ , μ = 1,1
501+
N = M = 4
502+
cols = rows = 1:N
503+
data = reshape(collect(1:(λ+μ+1)*(l+u+1)*sum(cols)), ((λ+μ+1)*(l+u+1), sum(cols)))
504+
A = _BandedBlockBandedMatrix(data, rows,cols, (l,u), (λ,μ))
505+
@test A == BandedBlockBandedMatrix(A, (2,1), (2,1)) == BandedBlockBandedMatrix{Float64}(A, (2,1), (2,1))
506+
end
490507
end
491508

492509
if false # turned off since tests have check-bounds=yes

‎test/test_blockbanded.jl‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -47,7 +47,7 @@ import Base.Broadcast: materialize!
4747

4848

4949
@test A[Block(1), Block(1)] isa Matrix
50-
@test A[Block(1), Block(1)] == A[Block(1,1)] == BlockArrays.getblock(A, 1, 1) == Matrix(view(A, Block(1,1)))
50+
@test A[Block(1), Block(1)] == A[Block(1,1)] == view(A, Block(1, 1)) == Matrix(view(A, Block(1,1)))
5151
@test A[1,1] == view(A,Block(1),Block(1))[1,1] == view(A,Block(1,1))[1,1] == A[Block(1,1)][1,1] == A[Block(1),Block(1)][1,1] == 1
5252
@test A[2,1] == view(A,Block(2),Block(1))[1,1] == view(A,Block(2,1))[1,1] == 2
5353
@test A[3,1] == view(A,Block(2),Block(1))[2,1] == 3

0 commit comments

Comments
 (0)