Skip to content

Commit cc5eb23

Browse files
authored
Special case blocks for non-blocked arrays (#474)
Closes #403.
1 parent d0e0e27 commit cc5eb23

File tree

2 files changed

+47
-12
lines changed

2 files changed

+47
-12
lines changed

src/blocks.jl

Lines changed: 27 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -49,6 +49,7 @@ blocks(a::AbstractArray) = BlocksView(a)
4949
blocks(a::BlockArray) = a.blocks
5050
blocks(A::Adjoint) = adjoint(blocks(parent(A)))
5151
blocks(A::Transpose) = transpose(blocks(parent(A)))
52+
blocks(A::StridedArray) = BlockView(A)
5253

5354
# convert a tuple of BlockRange to a tuple of `AbstractUnitRange{Int}`
5455
_block2int(B::Block{1}) = Int(B):Int(B)
@@ -63,7 +64,7 @@ struct BlocksView{
6364
S, # eltype(eltype(BlocksView(...)))
6465
N, # ndims
6566
T<:AbstractArray{S,N}, # eltype(BlocksView(...)), i.e., block type
66-
B<:AbstractArray{S,N}, # array to be wrapped
67+
B<:AbstractArray{S,N}, # array to be wrapped
6768
} <: AbstractArray{T,N}
6869
array::B
6970
end
@@ -78,13 +79,6 @@ Base.IteratorEltype(::Type{<:BlocksView}) = Base.EltypeUnknown()
7879
Base.size(a::BlocksView) = blocksize(a.array)
7980
Base.axes(a::BlocksView) = map(br -> Int.(br), blockaxes(a.array))
8081

81-
#=
82-
This is broken for now. See: https://github.com/JuliaArrays/BlockArrays.jl/issues/120
83-
# IndexLinear implementations
84-
@propagate_inbounds getindex(a::BlocksView, i::Int) = view(a.array, Block(i))
85-
@propagate_inbounds setindex!(a::BlocksView, b, i::Int) = copyto!(a[i], b)
86-
=#
87-
8882
# IndexCartesian implementations
8983
@propagate_inbounds getindex(a::BlocksView{T,N}, i::Vararg{Int,N}) where {T,N} =
9084
view(a.array, Block.(i)...)
@@ -93,6 +87,31 @@ This is broken for now. See: https://github.com/JuliaArrays/BlockArrays.jl/issue
9387
a
9488
end
9589

90+
# Like `BlocksView` but specialized for a single block
91+
# in order to avoid unnecessary wrappers when accessing the block.
92+
# Note that it does not check the array being wrapped actually
93+
# only has a single block, and will interpret it as if it just has one block.
94+
# By default, this is what gets constructed when calling `blocks(::StridedArray)`.
95+
struct BlockView{
96+
S, # eltype(eltype(BlockView(...)))
97+
N, # ndims
98+
T<:AbstractArray{S,N}, # array to be wrapped
99+
} <: AbstractArray{T,N}
100+
array::T
101+
end
102+
103+
Base.size(a::BlockView) = map(one, size(a.array))
104+
105+
# IndexCartesian implementations
106+
@propagate_inbounds function getindex(a::BlockView{T,N}, i::Vararg{Int,N}) where {T,N}
107+
@boundscheck checkbounds(a, i...)
108+
a.array
109+
end
110+
@propagate_inbounds function setindex!(a::BlockView{T,N}, b, i::Vararg{Int,N}) where {T,N}
111+
copyto!(a[i...], b)
112+
a
113+
end
114+
96115
# AbstractArray version of `Iterators.product`.
97116
# https://en.wikipedia.org/wiki/Cartesian_product
98117
# https://github.com/lazyLibraries/ProductArrays.jl

test/test_blocks.jl

Lines changed: 20 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,7 +1,7 @@
11
module TestBlocks
22

3-
using Test, BlockArrays
4-
import BlockArrays: eachblockaxes1
3+
using Test, BlockArrays, FillArrays, LinearAlgebra
4+
import BlockArrays: BlockView, eachblockaxes1
55

66
@testset "blocks" begin
77
@testset "blocks(::BlockVector)" begin
@@ -90,17 +90,33 @@ import BlockArrays: eachblockaxes1
9090
@testset "blocks(::Vector)" begin
9191
v = rand(3)
9292
@test size(blocks(v)) == (1,)
93+
@test blocks(v)[1] v
94+
@test blocks(v) isa BlockView{Float64,1,Vector{Float64}}
95+
@test blocks(v')[1, 1] v'
96+
@test blocks(v') isa Adjoint{Adjoint{Float64,Vector{Float64}},<:BlockView{Float64,1,Vector{Float64}}}
97+
@test blocks(view(v, 1:2))[1, 1] view(v, 1:2)
98+
@test blocks(view(v, 1:2)) isa BlockView{Float64,1,<:SubArray{Float64,1,Vector{Float64}}}
99+
blocks(v)[1] = zeros(3)
100+
@test iszero(v[1])
93101
blocks(v)[1][1] = 123
94102
@test v[1] == 123
95-
@test parent(blocks(v)[1]) === v
96103
end
97104

98105
@testset "blocks(::Matrix)" begin
99106
m = rand(2, 4)
100107
@test size(blocks(m)) == (1, 1)
108+
@test blocks(m)[1] m
109+
@test blocks(m)[1, 1] m
110+
@test blocks(m) isa BlockView{Float64,2,Matrix{Float64}}
111+
@test blocks(m')[1, 1] m'
112+
@test blocks(m') isa Adjoint{Adjoint{Float64,Matrix{Float64}},BlockView{Float64,2,Matrix{Float64}}}
113+
@test blocks(view(m, 1:2, 1:2))[1, 1] view(m, 1:2, 1:2)
114+
@test blocks(view(m, 1:2, 1:2)) isa BlockView{Float64,2,<:SubArray{Float64,2,Matrix{Float64}}}
115+
@test blocks(m)[1, 1] m
116+
blocks(m)[1, 1] = zeros(2, 4)
117+
@test iszero(m[1, 1])
101118
blocks(m)[1, 1][1, 1] = 123
102119
@test m[1, 1] == 123
103-
@test parent(blocks(m)[1, 1]) === m
104120
end
105121

106122
@testset "blocks(::Adjoint|Transpose)" begin

0 commit comments

Comments
 (0)