Skip to content

Commit 86572bb

Browse files
dlfivefiftydlfivefifty
authored
Add zeros/ones for BlockedUnitRange, Avoid calling axes for subarrays (#124)
* Add zeros/ones for BlockedUnitRange, Avoid calling axes for subarrays * Update test_blockviews.jl * Add tests * Remove unused code * v0.12.9
1 parent 0255982 commit 86572bb

File tree

8 files changed

+152
-19
lines changed

8 files changed

+152
-19
lines changed

‎.gitignore

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -7,3 +7,4 @@ docs/site/
77
benchmark/*.md
88
src/.DS_Store
99
/Manifest.toml
10+
.DS_Store

‎Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "BlockArrays"
22
uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
3-
version = "0.12.8"
3+
version = "0.12.9"
44

55
[deps]
66
ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"

‎src/BlockArrays.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,8 @@ import Base: @propagate_inbounds, Array, to_indices, to_index,
2727
@_inline_meta, _maybetail, tail, @_propagate_inbounds_meta, reindex,
2828
RangeIndex, Int, Integer, Number,
2929
+, -, *, /, \, min, max, isless, in, copy, copyto!, axes, @deprecate,
30-
BroadcastStyle, checkbounds, throw_boundserror
30+
BroadcastStyle, checkbounds, throw_boundserror,
31+
ones, zeros
3132
using Base: ReshapedArray, dataids
3233

3334

‎src/blockarrayinterface.jl

Lines changed: 16 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,20 @@
11

22
getindex(a::Number, ::Block{0}) = a
33

4+
function _sym_axes(A)
5+
ax = axes(parent(A),2)
6+
(ax, ax)
7+
end
8+
9+
# Symmetric and Triangular should inherit blocks from parent
10+
axes(A::HermOrSym{<:Any,<:AbstractBlockMatrix}) = _sym_axes(A)
11+
axes(A::HermOrSym{<:Any,<:SubArray{<:Any,2,<:AbstractBlockMatrix}}) = _sym_axes(A)
412
axes(A::AbstractTriangular{<:Any,<:AbstractBlockMatrix}) = axes(parent(A))
5-
axes(A::HermOrSym{<:Any,<:AbstractBlockMatrix}) = axes(parent(A))
13+
axes(A::AbstractTriangular{<:Any,<:SubArray{<:Any,2,<:AbstractBlockMatrix}}) = axes(parent(A))
14+
15+
blocksize(A::AbstractTriangular) = blocksize(parent(A))
16+
blocksize(A::AbstractTriangular, i::Int) = blocksize(parent(A), i)
17+
blockaxes(A::AbstractTriangular) = blockaxes(parent(A))
18+
19+
hasmatchingblocks(A::AbstractTriangular) = hasmatchingblocks(parent(A))
20+
hasmatchingblocks(A::HermOrSym) = true

‎src/pseudo_blockarray.jl

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -328,4 +328,12 @@ Base.unsafe_convert(::Type{Ptr{T}}, A::PseudoBlockArray) where T = Base.unsafe_c
328328
###
329329

330330
colsupport(A::PseudoBlockArray, j) = colsupport(A.blocks, j)
331-
rowsupport(A::PseudoBlockArray, j) = rowsupport(A.blocks, j)
331+
rowsupport(A::PseudoBlockArray, j) = rowsupport(A.blocks, j)
332+
333+
###
334+
# zeros/ones
335+
###
336+
337+
for op in (:zeros, :ones)
338+
@eval $op(::Type{T}, axs::Tuple{BlockedUnitRange,Vararg{Any}}) where T = PseudoBlockArray($op(T, map(length,axs)...), axs)
339+
end

‎src/views.jl

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -177,3 +177,24 @@ unsafe_convert(::Type{Ptr{T}}, V::SubArray{T,N,PseudoBlockArray{T,N,AT},<:Tuple{
177177
unsafe_convert(Ptr{T}, V.parent) + (Base.first_index(V)-1)*sizeof(T)
178178

179179

180+
# The default blocksize(V) is slow for views as it calls axes(V), which
181+
# allocates. Here we work around this.
182+
183+
_sub_blocksize() = ()
184+
_sub_blocksize(ind::BlockSlice{<:BlockRange{1}}, inds...) = tuple(length(ind.block),_sub_blocksize(inds...)...)
185+
_sub_blocksize(ind::BlockSlice{<:Block{1}}, inds...) = tuple(1,_sub_blocksize(inds...)...)
186+
_sub_blocksize(ind::AbstractVector, inds...) = tuple(blocksize(ind,1),_sub_blocksize(inds...)...)
187+
_sub_blocksize(ind::Integer, inds...) = _sub_blocksize(inds...)
188+
blocksize(V::SubArray) = _sub_blocksize(parentindices(V)...)
189+
blocksize(V::SubArray, i::Int) = _sub_blocksize(parentindices(V)[i])[1]
190+
191+
function hasmatchingblocks(V::SubArray{<:Any,2,<:Any,<:NTuple{2,BlockSlice{<:BlockRange{1}}}})
192+
a,b = axes(parent(V))
193+
kr,jr = parentindices(V)
194+
KR,JR = (kr.block),(jr.block)
195+
length(KR) == length(JR) || return false
196+
for (K,J) in zip(KR,JR)
197+
length(a[K]) == length(b[J]) || return false
198+
end
199+
true
200+
end

‎test/test_blockarrayinterface.jl

Lines changed: 55 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -30,21 +30,61 @@ end
3030
end
3131

3232
@testset "Triangular/Symmetric/Hermitian block arrays" begin
33-
A = PseudoBlockArray{ComplexF64}(undef, (1:4), (1:4))
34-
A .= reshape(1:length(A), size(A))
35-
36-
@test blocksize(UpperTriangular(A)) == blocksize(Symmetric(A)) == blocksize(A)
37-
@test UpperTriangular(A)[Block(2,2)] == UpperTriangular(A[2:3,2:3])
38-
@test UpperTriangular(A)[Block(2,3)] == A[2:3,4:6]
39-
@test UpperTriangular(A)[Block(3,2)] == zeros(3,2)
40-
@test Symmetric(A)[Block(2,2)] == Symmetric(A[2:3,2:3])
41-
@test Symmetric(A)[Block(2,3)] == A[2:3,4:6]
42-
@test Symmetric(A)[Block(3,2)] == transpose(A[2:3,4:6])
43-
@test Hermitian(A)[Block(2,2)] == Hermitian(A[2:3,2:3])
44-
@test Hermitian(A)[Block(2,3)] == A[2:3,4:6]
45-
@test Hermitian(A)[Block(3,2)] == A[2:3,4:6]'
46-
if VERSION ≥ v"1.2"
47-
@test stringmime("text/plain", UpperTriangular(A)) == "10×10 UpperTriangular{Complex{Float64},PseudoBlockArray{Complex{Float64},2,Array{Complex{Float64},2},Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}}}} with indices 1:1:10×1:1:10:\n 1.0+0.0im │ 11.0+0.0im 21.0+0.0im │ 31.0+0.0im 41.0+0.0im 51.0+0.0im │ 61.0+0.0im 71.0+0.0im 81.0+0.0im 91.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n ⋅ │ 12.0+0.0im 22.0+0.0im │ 32.0+0.0im 42.0+0.0im 52.0+0.0im │ 62.0+0.0im 72.0+0.0im 82.0+0.0im 92.0+0.0im\n ⋅ │ ⋅ 23.0+0.0im │ 33.0+0.0im 43.0+0.0im 53.0+0.0im │ 63.0+0.0im 73.0+0.0im 83.0+0.0im 93.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n ⋅ │ ⋅ ⋅ │ 34.0+0.0im 44.0+0.0im 54.0+0.0im │ 64.0+0.0im 74.0+0.0im 84.0+0.0im 94.0+0.0im\n ⋅ │ ⋅ ⋅ │ ⋅ 45.0+0.0im 55.0+0.0im │ 65.0+0.0im 75.0+0.0im 85.0+0.0im 95.0+0.0im\n ⋅ │ ⋅ ⋅ │ ⋅ ⋅ 56.0+0.0im │ 66.0+0.0im 76.0+0.0im 86.0+0.0im 96.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n ⋅ │ ⋅ ⋅ │ ⋅ ⋅ ⋅ │ 67.0+0.0im 77.0+0.0im 87.0+0.0im 97.0+0.0im\n ⋅ │ ⋅ ⋅ │ ⋅ ⋅ ⋅ │ ⋅ 78.0+0.0im 88.0+0.0im 98.0+0.0im\n ⋅ │ ⋅ ⋅ │ ⋅ ⋅ ⋅ │ ⋅ ⋅ 89.0+0.0im 99.0+0.0im\n ⋅ │ ⋅ ⋅ │ ⋅ ⋅ ⋅ │ ⋅ ⋅ ⋅ 100.0+0.0im"
33+
@testset "square blocks" begin
34+
A = PseudoBlockArray{ComplexF64}(undef, 1:4, 1:4)
35+
A .= reshape(1:length(A), size(A))
36+
U = UpperTriangular(A)
37+
S = Symmetric(A)
38+
H = Hermitian(A)
39+
40+
@test blocksize(U) == blocksize(S) == blocksize(A)
41+
@test blockaxes(U) == blockaxes(S) == blockaxes(A)
42+
43+
@test U[Block(2,2)] == UpperTriangular(A[2:3,2:3])
44+
@test U[Block(2,3)] == A[2:3,4:6]
45+
@test U[Block(3,2)] == zeros(3,2)
46+
@test S[Block(2,2)] == Symmetric(A[2:3,2:3])
47+
@test S[Block(2,3)] == A[2:3,4:6]
48+
@test S[Block(3,2)] == transpose(A[2:3,4:6])
49+
@test H[Block(2,2)] == Hermitian(A[2:3,2:3])
50+
@test H[Block(2,3)] == A[2:3,4:6]
51+
@test H[Block(3,2)] == A[2:3,4:6]'
52+
53+
if VERSION ≥ v"1.2"
54+
@test stringmime("text/plain", UpperTriangular(A)) == "10×10 UpperTriangular{Complex{Float64},PseudoBlockArray{Complex{Float64},2,Array{Complex{Float64},2},Tuple{BlockedUnitRange{Array{Int64,1}},BlockedUnitRange{Array{Int64,1}}}}} with indices 1:1:10×1:1:10:\n 1.0+0.0im │ 11.0+0.0im 21.0+0.0im │ 31.0+0.0im 41.0+0.0im 51.0+0.0im │ 61.0+0.0im 71.0+0.0im 81.0+0.0im 91.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n ⋅ │ 12.0+0.0im 22.0+0.0im │ 32.0+0.0im 42.0+0.0im 52.0+0.0im │ 62.0+0.0im 72.0+0.0im 82.0+0.0im 92.0+0.0im\n ⋅ │ ⋅ 23.0+0.0im │ 33.0+0.0im 43.0+0.0im 53.0+0.0im │ 63.0+0.0im 73.0+0.0im 83.0+0.0im 93.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n ⋅ │ ⋅ ⋅ │ 34.0+0.0im 44.0+0.0im 54.0+0.0im │ 64.0+0.0im 74.0+0.0im 84.0+0.0im 94.0+0.0im\n ⋅ │ ⋅ ⋅ │ ⋅ 45.0+0.0im 55.0+0.0im │ 65.0+0.0im 75.0+0.0im 85.0+0.0im 95.0+0.0im\n ⋅ │ ⋅ ⋅ │ ⋅ ⋅ 56.0+0.0im │ 66.0+0.0im 76.0+0.0im 86.0+0.0im 96.0+0.0im\n ───────────┼──────────────────────────┼──────────────────────────────────────┼─────────────────────────────────────────────────\n ⋅ │ ⋅ ⋅ │ ⋅ ⋅ ⋅ │ 67.0+0.0im 77.0+0.0im 87.0+0.0im 97.0+0.0im\n ⋅ │ ⋅ ⋅ │ ⋅ ⋅ ⋅ │ ⋅ 78.0+0.0im 88.0+0.0im 98.0+0.0im\n ⋅ │ ⋅ ⋅ │ ⋅ ⋅ ⋅ │ ⋅ ⋅ 89.0+0.0im 99.0+0.0im\n ⋅ │ ⋅ ⋅ │ ⋅ ⋅ ⋅ │ ⋅ ⋅ ⋅ 100.0+0.0im"
55+
end
56+
57+
V = view(A, Block.(1:2), Block.(1:2))
58+
@test blockisequal(axes(Symmetric(V)), axes(view(A,Block.(1:2), Block.(1:2))))
59+
end
60+
61+
@testset "rect blocks" begin
62+
A = PseudoBlockArray{ComplexF64}(undef, 1:3, fill(3,2))
63+
A .= randn(6,6) .+ im * randn(6,6)
64+
U = UpperTriangular(A)
65+
S = Symmetric(A)
66+
H = Hermitian(A)
67+
68+
@test blockisequal(axes(S), (axes(A,2),axes(A,2)))
69+
@test blockisequal(axes(U), axes(A))
70+
@test blocksize(S) == (blocksize(S,1),blocksize(S,2)) == (2,2)
71+
@test blocksize(U) == (blocksize(U,1),blocksize(U,2)) == blocksize(A)
72+
@test blockaxes(S) == (blockaxes(S,1),blockaxes(S,2)) == (Block.(1:2), Block.(1:2))
73+
@test blockaxes(U) == (blockaxes(U,1),blockaxes(U,2)) == blockaxes(A)
74+
75+
@test U[Block(2,2)] == A[Block(2,2)]
76+
@test U[Block(3,2)] == triu(A[Block(3,2)])
77+
@test U[Block(3,1)] == zeros(3,3)
78+
@test S[Block(2,2)] == Symmetric(A[4:6,4:6])
79+
@test S[Block(1,2)] == A[1:3,4:6]
80+
@test S[Block(2,1)] == transpose(A[1:3,4:6])
81+
@test H[Block(2,2)] == Hermitian(A[4:6,4:6])
82+
@test H[Block(1,2)] == A[1:3,4:6]
83+
@test H[Block(2,1)] == A[1:3,4:6]'
84+
85+
@test !BlockArrays.hasmatchingblocks(U)
86+
@test BlockArrays.hasmatchingblocks(S)
87+
@test BlockArrays.hasmatchingblocks(H)
4888
end
4989
end
5090

‎test/test_blockviews.jl

Lines changed: 47 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -174,4 +174,51 @@ using BlockArrays, Test, Base64
174174
@test A[Block.(2:3)] isa PseudoBlockArray
175175
@test A[Block.(2:3)] == A[2:end]
176176
end
177+
178+
@testset "non-allocation blocksize" begin
179+
A = BlockArray(randn(5050), 1:100)
180+
@test blocksize(A) == (100,)
181+
@test @allocated(blocksize(A)) ≤ 40
182+
V = view(A, Block(3))
183+
@test blocksize(V) == (1,)
184+
@test @allocated(blocksize(V)) ≤ 40
185+
V = view(A, Block.(1:30))
186+
@test blocksize(V) == (30,)
187+
@test @allocated(blocksize(V)) ≤ 40
188+
V = view(A, 3:43)
189+
@test blocksize(V) == (1,)
190+
V = view(A, 5)
191+
@test blocksize(V) == ()
192+
193+
A = BlockArray(randn(5050,21), 1:100, 1:6)
194+
@test blocksize(A) == (100,6)
195+
@test @allocated(blocksize(A)) ≤ 40
196+
V = view(A, Block(3,2))
197+
@test blocksize(V) == (1,1)
198+
@test @allocated(blocksize(V)) ≤ 40
199+
V = view(A, Block.(1:30), Block(3))
200+
@test blocksize(V) == (30,1)
201+
@test @allocated(blocksize(V)) ≤ 40
202+
V = view(A, Block.(1:30), Block.(1:3))
203+
@test blocksize(V) == (30,3)
204+
@test @allocated(blocksize(V)) ≤ 40
205+
V = view(A, 3:43,1:3)
206+
@test blocksize(V) == (1,1)
207+
V = view(A, 5, 1:3)
208+
@test blocksize(V) == (1,)
209+
end
210+
211+
@testset "hasmatchingblocks" begin
212+
A = BlockArray{Int}(undef, 1:20, 1:20)
213+
B = BlockArray{Int}(undef, 1:3, fill(3,2))
214+
V = view(A,Block.(1:10),Block.(1:10))
215+
216+
@test BlockArrays.hasmatchingblocks(A)
217+
@test BlockArrays.hasmatchingblocks(V)
218+
@test @allocated(BlockArrays.hasmatchingblocks(V)) == 0
219+
@test !BlockArrays.hasmatchingblocks(view(A,Block.(1:2),1:3))
220+
@test !BlockArrays.hasmatchingblocks(view(A,Block.(1:2),Block.(2:3)))
221+
222+
@test BlockArrays.hasmatchingblocks(view(B,Block.(3:3),Block.(2:2)))
223+
end
177224
end

0 commit comments

Comments
 (0)