Skip to content

Commit 657b480

Browse files
dlfivefiftydlfivefifty
authored
Fixes for new BlockArrays (#55)
* Fixes for new BlockArrays * Improve sublayout * Add tests * Update Project.toml
1 parent 4d3e366 commit 657b480

File tree

8 files changed

+161
-66
lines changed

8 files changed

+161
-66
lines changed

‎Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "BlockBandedMatrices"
22
uuid = "ffab5731-97b5-5995-9138-79e8c1846df0"
3-
version = "0.7"
3+
version = "0.7.1"
44

55
[deps]
66
BandedMatrices = "aae01518-5342-5314-be14-df237901396f"

‎src/AbstractBlockBandedMatrix.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -117,7 +117,7 @@ const BlockSlice1 = BlockSlice{Block{1,Int}}
117117
######################################
118118

119119
@inline colstart(A::AbstractBlockBandedMatrix, i::Integer) =
120-
first(axes(A,1)[blockcolstart(A,findblock(axes(A,2),i))])
120+
isempty(axes(A,1)) ? 1 : first(axes(A,1)[blockcolstart(A,findblock(axes(A,2),i))])
121121

122122
@inline function colstop(A::AbstractBlockBandedMatrix, i::Integer)
123123
CS = blockcolstop(A,findblock(axes(A,2),i))

‎src/BandedBlockBandedMatrix.jl

Lines changed: 116 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -85,6 +85,9 @@ function convert(::Type{<:BandedBlockBandedMatrix}, B::BandedMatrix)
8585
end
8686
end
8787

88+
convert(::Type{BandedBlockBandedMatrix{T,BLOCKS,RAXIS}}, A::BandedBlockBandedMatrix) where {T,BLOCKS,RAXIS} =
89+
_BandedBlockBandedMatrix(convert(BLOCKS, A.data), convert(RAXIS, A.raxis), (A.l, A.u), (A.λ, A.μ))
90+
8891
function BandedBlockBandedMatrix{T,B,R}(Z::Zeros, axes::NTuple{2,AbstractUnitRange{Int}},
8992
lu::NTuple{2,Int}, λμ::NTuple{2,Int}) where {T,B,R<:AbstractUnitRange{Int}}
9093
if size(Z) ≠ map(length,axes)
@@ -191,10 +194,7 @@ similar(A::BandedBlockBandedMatrix, ::Type{T}, axes::NTuple{2,AbstractUnitRange{
191194

192195

193196
similar(A::BandedBlockBandedMatrix{T}, axes::NTuple{2,AbstractUnitRange{Int}}) where T =
194-
BandedBlockBandedMatrix{T}(undef, axes, blockbandwidths(A), subblockbandwidths(A))
195-
196-
similar(A::BandedBlockBandedMatrix{T}, axes::NTuple{2,OneTo{Int}}) where T =
197-
Matrix{T}(undef, map(length,axes))
197+
similar(Matrix{T}, map(length,axes)...)
198198

199199
axes(A::BandedBlockBandedMatrix) = (A.raxis, axes(A.data,2))
200200

@@ -398,21 +398,91 @@ end
398398
# with BLASBandedMatrix.
399399
##################
400400

401+
const SubBandedBlockBandedMatrix{T,R1,R2} =
402+
SubArray{T,2,<:BandedBlockBandedMatrix{T},<:Tuple{<:BlockSlice{R1},<:BlockSlice{R2}}}
403+
404+
405+
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{Block1}}}) = BandedColumnMajor()
406+
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) = BandedColumnMajor()
407+
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) = BandedColumnMajor()
408+
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) = BandedColumnMajor()
409+
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockRange1}}}) = BandedBlockBandedColumnMajor()
410+
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{Block1},BlockSlice{BlockRange1}}}) = BandedBlockBandedColumnMajor()
411+
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{Block1}}}) = BandedBlockBandedColumnMajor()
412+
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{BlockSlice{BlockRange1},BlockSlice{BlockIndexRange1}}}) = BandedBlockBandedColumnMajor()
413+
414+
isbanded(A::SubArray{<:Any,2,<:BandedBlockBandedMatrix}) = MemoryLayout(typeof(A)) == BandedColumnMajor()
415+
isbandedblockbanded(A::SubArray{<:Any,2,<:BandedBlockBandedMatrix}) = MemoryLayout(typeof(A)) == BandedBlockBandedColumnMajor()
416+
417+
418+
subblockbandwidths(V::SubBandedBlockBandedMatrix) = subblockbandwidths(parent(V))
419+
420+
function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,BlockRange1,Block1})
421+
A = parent(V)
422+
423+
KR = parentindices(V)[1].block.indices[1]
424+
J = parentindices(V)[2].block
425+
shift = Int(KR[1])-Int(J)
426+
blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
427+
end
428+
429+
function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,Block1,BlockRange1})
430+
A = parent(V)
431+
432+
K = parentindices(V)[1].block
433+
JR = parentindices(V)[2].block.indices[1]
434+
shift = Int(K)-Int(JR[1])
435+
436+
blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
437+
end
438+
439+
function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,BlockRange1,BlockRange1})
440+
A = parent(V)
441+
442+
KR = parentindices(V)[1].block.indices[1]
443+
JR = parentindices(V)[2].block.indices[1]
444+
shift = Int(KR[1])-Int(JR[1])
445+
446+
blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
447+
end
448+
449+
401450
const BandedBlockBandedBlock{T, BLOCKS, RAXIS} = SubArray{T,2,BandedBlockBandedMatrix{T, BLOCKS, RAXIS},<:Tuple{<:BlockSlice1,<:BlockSlice1},false}
402451

403452

404-
isbanded(::BandedBlockBandedBlock) = true
405-
MemoryLayout(::BandedBlockBandedBlock) = BandedColumnMajor()
406453
BroadcastStyle(::Type{<: BandedBlockBandedBlock}) = BandedStyle()
407454

408455

409-
function inblockbands(V::BandedBlockBandedBlock)
456+
function inblockbands(V::SubArray{<:Any,2,<:AbstractMatrix,<:Tuple{<:BlockSlice1,<:BlockSlice1},false})
410457
A = parent(V)
411458
K_sl, J_sl = parentindices(V)
412459
K, J = K_sl.block, J_sl.block
413-
-A.l ≤ Int(J-K) ≤ A.u
460+
l,u = blockbandwidths(A)
461+
-l ≤ Int(J-K) ≤ u
414462
end
415463

464+
function parentblock(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) where T
465+
A = parent(V)
466+
K_sl, J_sl = parentindices(V)
467+
view(A, K_sl.block.block, J_sl.block.block)
468+
end
469+
470+
function parentblock(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) where T
471+
A = parent(V)
472+
K_sl, J_sl = parentindices(V)
473+
view(A, K_sl.block, J_sl.block.block)
474+
end
475+
476+
function parentblock(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) where T
477+
A = parent(V)
478+
K_sl, J_sl = parentindices(V)
479+
view(A, K_sl.block.block, J_sl.block)
480+
end
481+
482+
# gives the columns of parent(V).data that encode the block
483+
parentblocks2Int(V::BandedBlockBandedBlock)::Tuple{Int,Int} = Int(first(parentindices(V)).block),
484+
Int(last(parentindices(V)).block)
485+
416486

417487
######################################
418488
# BandedMatrix interface for Blocks #
@@ -422,11 +492,23 @@ end
422492
(-720,-720)
423493
end
424494

495+
function bandwidths(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) where T
496+
B = parentblock(V)
497+
K_sl, J_sl = parentindices(V)
498+
bandwidths(B) .+ (-1,1) .* bandshift(K_sl.block.indices[1],J_sl.block.indices[1])
499+
end
425500

501+
function bandwidths(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) where T
502+
B = parentblock(V)
503+
K_sl, J_sl = parentindices(V)
504+
bandwidths(B) .+ (-1,1) .* bandshift(Base.OneTo(1),J_sl.block.indices[1])
505+
end
426506

427-
# gives the columns of parent(V).data that encode the block
428-
blocks(V::BandedBlockBandedBlock)::Tuple{Int,Int} = Int(first(parentindices(V)).block),
429-
Int(last(parentindices(V)).block)
507+
function bandwidths(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) where T
508+
B = parentblock(V)
509+
K_sl, J_sl = parentindices(V)
510+
bandwidths(B) .+ (-1,1) .* bandshift(K_sl.block.indices[1],Base.OneTo(1))
511+
end
430512

431513

432514
function bandeddata(V::BandedBlockBandedBlock{T}) where T
@@ -438,6 +520,27 @@ function bandeddata(V::BandedBlockBandedBlock{T}) where T
438520
view(A.data, u + K - J + 1, J)
439521
end
440522

523+
function bandeddata(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{BlockIndexRange1}}}) where T
524+
A = parent(V)
525+
K_sl, J_sl = parentindices(V)
526+
view(bandeddata(parentblock(V)), :, J_sl.block.indices[1])
527+
end
528+
529+
function bandeddata(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{Block1},BlockSlice{BlockIndexRange1}}}) where T
530+
A = parent(V)
531+
K_sl, J_sl = parentindices(V)
532+
view(bandeddata(parentblock(V)), :, J_sl.block.indices[1])
533+
end
534+
535+
function bandeddata(V::SubArray{T,2,<:AbstractMatrix,<:Tuple{BlockSlice{BlockIndexRange1},BlockSlice{Block1}}}) where T
536+
A = parent(V)
537+
K_sl, J_sl = parentindices(V)
538+
bandeddata(parentblock(V))
539+
end
540+
541+
542+
543+
441544

442545
@inline function inbands_getindex(V::BandedBlockBandedBlock, k::Int, j::Int)
443546
A = parent(V)
@@ -453,7 +556,7 @@ end
453556
@propagate_inbounds function getindex(V::BandedBlockBandedBlock, k::Int, j::Int)
454557
@boundscheck checkbounds(V, k, j)
455558
A = parent(V)
456-
K,J = blocks(V)
559+
K,J = parentblocks2Int(V)
457560
if -A.l ≤ J-K ≤ A.u && -A.λ ≤ j-k ≤ A.μ
458561
inbands_getindex(V, k, j)
459562
else
@@ -464,7 +567,7 @@ end
464567
@propagate_inbounds function setindex!(V::BandedBlockBandedBlock, v, k::Int, j::Int)
465568
@boundscheck checkbounds(V, k, j)
466569
A = parent(V)
467-
K,J = blocks(V)
570+
K,J = parentblocks2Int(V)
468571
if -A.l ≤ J-K ≤ A.u && -A.λ ≤ j-k ≤ A.μ
469572
inbands_setindex!(V, v, k, j)
470573
elseif iszero(v) # allow setindex for 0 datya

‎src/BlockBandedMatrices.jl

Lines changed: 5 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -38,13 +38,17 @@ import BandedMatrices: isbanded, bandwidths, bandwidth, banded_getindex, colrang
3838
BandedColumnMajor,
3939
BandedSubBandedMatrix, bandeddata, tribandeddata,
4040
_BandedMatrix, colstart, colstop, rowstart, rowstop,
41-
BandedStyle, _fill_lmul!,
41+
BandedStyle, _fill_lmul!, bandshift,
4242
_banded_colval, _banded_rowval, _banded_nzval # for sparse
4343

4444
export BandedBlockBandedMatrix, BlockBandedMatrix, BlockSkylineMatrix, blockbandwidth, blockbandwidths,
4545
subblockbandwidth, subblockbandwidths, Ones, Zeros, Fill, Block, BlockTridiagonal, isblockbanded
4646

4747

48+
const Block1 = Block{1,Int}
49+
const BlockRange1 = BlockRange{1,Tuple{UnitRange{Int}}}
50+
const BlockIndexRange1 = BlockIndexRange{1,Tuple{UnitRange{Int}}}
51+
4852
include("AbstractBlockBandedMatrix.jl")
4953
include("broadcast.jl")
5054
include("BlockSkylineMatrix.jl")

‎src/broadcast.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -212,7 +212,7 @@ function blockbanded_axpy!(a, X::AbstractMatrix, Y::AbstractMatrix)
212212
end
213213

214214
function _combine_blockaxes(A, B)
215-
blockisequal(axes(A), axes(B)) || throw(DimensionMismatch())
215+
blockisequal(axes(A), axes(B)) || throw(DimensionMismatch("Block sizes do not agree"))
216216
axes(A)
217217
end
218218

‎src/linalg.jl

Lines changed: 0 additions & 49 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,7 @@
11
# BlockBandedMatrix with block range indexes is also block-banded
2-
const Block1 = Block{1,Int}
3-
const BlockRange1 = BlockRange{1,Tuple{UnitRange{Int}}}
4-
const BlockIndexRange1 = BlockIndexRange{1,Tuple{UnitRange{Int}}}
52
const SubBlockSkylineMatrix{T,LL,UU,R1,R2} =
63
SubArray{T,2,BlockSkylineMatrix{T,LL,UU},<:Tuple{<:BlockSlice{R1},<:BlockSlice{R2}}}
74

8-
const SubBandedBlockBandedMatrix{T,R1,R2} =
9-
SubArray{T,2,<:BandedBlockBandedMatrix{T},<:Tuple{<:BlockSlice{R1},<:BlockSlice{R2}}}
10-
11-
125

136
getindex(A::BandedBlockBandedMatrix, KR::BlockRange1, JR::BlockRange1) = BandedBlockBandedMatrix(view(A, KR, JR))
147
getindex(A::BandedBlockBandedMatrix, KR::BlockRange1, J::Block1) = BandedBlockBandedMatrix(view(A, KR, J))
@@ -180,48 +173,6 @@ sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:
180173
sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{BlockIndexRange1}}}) = ColumnMajor()
181174
sublayout(::BlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockIndexRange1}, <:BlockSlice{BlockIndexRange1}}}) = ColumnMajor()
182175

183-
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{Block1}, <:BlockSlice{Block1}}}) = BandedColumnMajor()
184-
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{BlockRange1}}}) = BandedBlockBandedColumnMajor()
185-
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{Block1}, <:BlockSlice{BlockRange1}}}) = BandedBlockBandedColumnMajor()
186-
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{Block1}}}) = BandedBlockBandedColumnMajor()
187-
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockRange1}, <:BlockSlice{BlockIndexRange1}}}) = BandedBlockBandedColumnMajor()
188-
sublayout(::BandedBlockBandedColumnMajor, ::Type{<:Tuple{<:BlockSlice{BlockIndexRange1}, <:BlockSlice{BlockIndexRange1}}}) = BandedBlockBandedColumnMajor()
189-
190-
isbanded(A::SubArray{<:Any,2,<:BandedBlockBandedMatrix}) = MemoryLayout(typeof(A)) == BandedColumnMajor()
191-
isbandedblockbanded(A::SubArray{<:Any,2,<:BandedBlockBandedMatrix}) = MemoryLayout(typeof(A)) == BandedBlockBandedColumnMajor()
192-
193-
194-
subblockbandwidths(V::SubBandedBlockBandedMatrix) = subblockbandwidths(parent(V))
195-
196-
function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,BlockRange1,Block1})
197-
A = parent(V)
198-
199-
KR = parentindices(V)[1].block.indices[1]
200-
J = parentindices(V)[2].block
201-
shift = Int(KR[1])-Int(J)
202-
blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
203-
end
204-
205-
function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,Block1,BlockRange1})
206-
A = parent(V)
207-
208-
K = parentindices(V)[1].block
209-
JR = parentindices(V)[2].block.indices[1]
210-
shift = Int(K)-Int(JR[1])
211-
212-
blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
213-
end
214-
215-
function blockbandwidths(V::SubBandedBlockBandedMatrix{<:Any,BlockRange1,BlockRange1})
216-
A = parent(V)
217-
218-
KR = parentindices(V)[1].block.indices[1]
219-
JR = parentindices(V)[2].block.indices[1]
220-
shift = Int(KR[1])-Int(JR[1])
221-
222-
blockbandwidth(A,1) - shift, blockbandwidth(A,2) + shift
223-
end
224-
225176

226177
strides(V::SubBlockSkylineMatrix{<:Any,LL,UU,<:Union{BlockRange1,Block1},Block1}) where {LL,UU} =
227178
(1,parent(V).block_sizes.block_strides[Int(parentindices(V)[2].block)])

‎test/test_bandedblockbanded.jl

Lines changed: 24 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -444,6 +444,30 @@ import BlockBandedMatrices: _BandedBlockBandedMatrix, blockcolrange, blockrowran
444444
A = BandedBlockBandedMatrix{Float64}(undef, 1:5,1:5, (-1,1), (-1,1))
445445
@test size(sparse(A)) == size(A) == (15,15)
446446
end
447+
448+
@testset "non-standard blocks" begin
449+
A = BandedBlockBandedMatrix{Float64}(undef, Int[], 1:5,(-1,1), (-1,1))
450+
@test BlockBandedMatrices.colstart(A,1) == 1
451+
A = BandedBlockBandedMatrix{Float64}(undef, 1:2, 1:5,(-1,1), (-1,1))
452+
A.data .= randn.()
453+
V = view(A, Block(2,3))
454+
@test MemoryLayout(typeof(V)) == BandedMatrices.BandedColumnMajor()
455+
@test isbanded(V)
456+
@test bandwidths(V) == (-1,1)
457+
@test BandedMatrix(V) == A[2:3,4:6]
458+
V2 = view(V, :, 2:3)
459+
@test MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumnMajor()
460+
@test bandwidths(V2) == (0,0)
461+
@test BandedMatrix(V2) == A[2:3,5:6]
462+
V2 = view(V, 2:2, :)
463+
@test MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumnMajor()
464+
@test bandwidths(V2) == (-2,2)
465+
@test BandedMatrix(V2) == A[3:3,4:6]
466+
V2 = view(V, 2:2, 2:3)
467+
@test MemoryLayout(typeof(V2)) == BandedMatrices.BandedColumnMajor()
468+
@test bandwidths(V2) == (-1,1)
469+
@test BandedMatrix(V2) == A[3:3,5:6]
470+
end
447471
end
448472

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

‎test/test_linalg.jl

Lines changed: 13 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -136,6 +136,19 @@ end
136136

137137
@time BLAS.axpy!(1.0, A, B)
138138
@test B ≈ AB
139+
140+
@testset "degenerate" begin
141+
A = BandedBlockBandedMatrix{Float64}(undef, 1:5, 1:5,(-1,1), (-1,1))
142+
B = BandedBlockBandedMatrix{Float64}(undef, 1:5, 1:5,(1,-1), (1,-1))
143+
A.data .= randn.()
144+
B.data .= randn.()
145+
@test blockbandwidths(A^2) == subblockbandwidths(A^2) == (-2,2)
146+
@test A^2 == Matrix(A)^2
147+
@test blockbandwidths(B^2) == subblockbandwidths(B^2) == (2,-2)
148+
@test B^2 == Matrix(B)^2
149+
@test blockbandwidths(A*B) == subblockbandwidths(A*B) == (0,0)
150+
@test A*B == Matrix(A)*Matrix(B)
151+
end
139152
end
140153

141154
@testset "Rectangular block *" begin

0 commit comments

Comments
 (0)