Skip to content

Special case blocks for non-blocked arrays #474

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
35 changes: 27 additions & 8 deletions src/blocks.jl
Original file line number Diff line number Diff line change
Expand Up @@ -49,6 +49,7 @@
blocks(a::BlockArray) = a.blocks
blocks(A::Adjoint) = adjoint(blocks(parent(A)))
blocks(A::Transpose) = transpose(blocks(parent(A)))
blocks(A::StridedArray) = BlockView(A)

# convert a tuple of BlockRange to a tuple of `AbstractUnitRange{Int}`
_block2int(B::Block{1}) = Int(B):Int(B)
Expand All @@ -63,7 +64,7 @@
S, # eltype(eltype(BlocksView(...)))
N, # ndims
T<:AbstractArray{S,N}, # eltype(BlocksView(...)), i.e., block type
B<:AbstractArray{S,N}, # array to be wrapped
B<:AbstractArray{S,N}, # array to be wrapped
} <: AbstractArray{T,N}
array::B
end
Expand All @@ -78,13 +79,6 @@
Base.size(a::BlocksView) = blocksize(a.array)
Base.axes(a::BlocksView) = map(br -> Int.(br), blockaxes(a.array))

#=
This is broken for now. See: https://github.com/JuliaArrays/BlockArrays.jl/issues/120
# IndexLinear implementations
@propagate_inbounds getindex(a::BlocksView, i::Int) = view(a.array, Block(i))
@propagate_inbounds setindex!(a::BlocksView, b, i::Int) = copyto!(a[i], b)
=#

# IndexCartesian implementations
@propagate_inbounds getindex(a::BlocksView{T,N}, i::Vararg{Int,N}) where {T,N} =
view(a.array, Block.(i)...)
Expand All @@ -93,6 +87,31 @@
a
end

# Like `BlocksView` but specialized for a single block
# in order to avoid unnecessary wrappers when accessing the block.
# Note that it does not check the array being wrapped actually
# only has a single block, and will interpret it as if it just has one block.
# By default, this is what gets constructed when calling `blocks(::StridedArray)`.
struct BlockView{
S, # eltype(eltype(BlockView(...)))
N, # ndims
T<:AbstractArray{S,N}, # array to be wrapped
} <: AbstractArray{T,N}
array::T
end

Base.size(a::BlockView) = map(one, size(a.array))

# IndexCartesian implementations
@propagate_inbounds function getindex(a::BlockView{T,N}, i::Vararg{Int,N}) where {T,N}
@boundscheck checkbounds(a, i...)
a.array
end
@propagate_inbounds function setindex!(a::BlockView{T,N}, b, i::Vararg{Int,N}) where {T,N}
copyto!(a[i...], b)
a

Check warning on line 112 in src/blocks.jl

View check run for this annotation

Codecov / codecov/patch

src/blocks.jl#L112

Added line #L112 was not covered by tests
end

# AbstractArray version of `Iterators.product`.
# https://en.wikipedia.org/wiki/Cartesian_product
# https://github.com/lazyLibraries/ProductArrays.jl
Expand Down
24 changes: 20 additions & 4 deletions test/test_blocks.jl
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
module TestBlocks

using Test, BlockArrays
import BlockArrays: eachblockaxes1
using Test, BlockArrays, FillArrays, LinearAlgebra
import BlockArrays: BlockView, eachblockaxes1

@testset "blocks" begin
@testset "blocks(::BlockVector)" begin
Expand Down Expand Up @@ -90,17 +90,33 @@ import BlockArrays: eachblockaxes1
@testset "blocks(::Vector)" begin
v = rand(3)
@test size(blocks(v)) == (1,)
@test blocks(v)[1] ≡ v
@test blocks(v) isa BlockView{Float64,1,Vector{Float64}}
@test blocks(v')[1, 1] ≡ v'
@test blocks(v') isa Adjoint{Adjoint{Float64,Vector{Float64}},<:BlockView{Float64,1,Vector{Float64}}}
@test blocks(view(v, 1:2))[1, 1] ≡ view(v, 1:2)
@test blocks(view(v, 1:2)) isa BlockView{Float64,1,<:SubArray{Float64,1,Vector{Float64}}}
blocks(v)[1] = zeros(3)
@test iszero(v[1])
blocks(v)[1][1] = 123
@test v[1] == 123
@test parent(blocks(v)[1]) === v
end

@testset "blocks(::Matrix)" begin
m = rand(2, 4)
@test size(blocks(m)) == (1, 1)
@test blocks(m)[1] ≡ m
@test blocks(m)[1, 1] ≡ m
@test blocks(m) isa BlockView{Float64,2,Matrix{Float64}}
@test blocks(m')[1, 1] ≡ m'
@test blocks(m') isa Adjoint{Adjoint{Float64,Matrix{Float64}},BlockView{Float64,2,Matrix{Float64}}}
@test blocks(view(m, 1:2, 1:2))[1, 1] ≡ view(m, 1:2, 1:2)
@test blocks(view(m, 1:2, 1:2)) isa BlockView{Float64,2,<:SubArray{Float64,2,Matrix{Float64}}}
@test blocks(m)[1, 1] ≡ m
blocks(m)[1, 1] = zeros(2, 4)
@test iszero(m[1, 1])
blocks(m)[1, 1][1, 1] = 123
@test m[1, 1] == 123
@test parent(blocks(m)[1, 1]) === m
end

@testset "blocks(::Adjoint|Transpose)" begin
Expand Down
Loading