getblock -> view (#139)

* ovlerload replace_in_print_matrix for PseudoBlockArray

* Update test/test_blockarrays.jl

Co-authored-by: Kristoffer Carlsson <[email protected]>

* test Static blocks

* Changes to support special views, non--allocating block axes

* make simplified, fast-tracked vector broadcast code

* blockbroadcast can depend on axes

* fix tests

* more tests pass, require Julia v1.5

* Start getblock -> view

* tests pass!

* Deprecations

* Add deprecations

* Update blocklinalg.jl

* support smarter view(view(A, Block.(1:5)), Block(5))

* general view of view of AbstractBlockArray

* More BlockStyle overloads

* Tests pass

* make view reindexing non-allocating

* Update Project.toml

* increase coverage, update docs

* Update test_blockbroadcast.jl

* increase coverage

* Remove single testset to avoid timeout on Travis

* increase coverage

* increase broadcast test coverage

* BlockIndexRange <: AbstractArray

Co-authored-by: Kristoffer Carlsson <[email protected]>

Author: dlfivefifty

.travis.yml

@@ -4,7 +4,6 @@ os:
   - linux
   - osx
 julia:
-  - 1.1
   - 1.5
   - nightly
 matrix:
@@ -23,7 +22,7 @@ after_success:
 jobs:
   include:
     - stage: "Documentation"
-      julia: 1.4
+      julia: 1.5
       os: linux
       script:
         - julia --project=docs/ --color=yes -e '

Project.toml

@@ -1,25 +1,25 @@
 name = "BlockArrays"
 uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
-version = "0.13.0"
+version = "0.14.0"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
-Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [compat]
-ArrayLayouts = "0.4.12"
-Compat = "2.2, 3"
-FillArrays = "0.10.1"
-julia = "1.1"
+ArrayLayouts = "0.5"
+FillArrays = "0.11"
+julia = "1.5"
 
 [extras]
 Base64 = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
+LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
 OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Base64", "FillArrays", "OffsetArrays", "SparseArrays", "Test"]
+test = ["Base64", "FillArrays", "LazyArrays", "OffsetArrays", "SparseArrays", "StaticArrays", "Test"]

docs/src/lib/public.md

@@ -38,9 +38,6 @@ blocklasts
 blocklengths
 blocks
 eachblock
-getblock
-getblock!
-setblock!
 blockcheckbounds
 ```
 

docs/src/man/abstractblockarrayinterface.md

@@ -40,9 +40,7 @@ It is possible to override additional functions to improve speed, however.
 | Methods to implement    | Brief description |
 | :---------------------- | :---------------- |
 | **Optional methods**    |           
-| `getblock(A, i...)`     | `X[Block(i...)]`, blocked indexing  |
-| `setblock!(A, v, i...)` | `X[Block(i...)] = v`, blocked index assignment |
-| `getblock!(x, A, i)`    | `X[i]`, blocked index assignment with in place storage in `x` |
+| `BlockArrays.viewblock(A, i::Block)`     | Specialised non-allocating `X[Block(i...)]`, blocked indexing  |
 
 For a more thorough description of the methods see the public interface documentation.
 

docs/src/man/blockarrays.md

@@ -87,8 +87,7 @@ julia> BlockArray(undef_blocks, SparseVector{Float64, Int}, [1,2])
 
 ## [Setting and getting blocks and values](@id setting_and_getting)
 
-A block can be set by `setblock!(block_array, v, i...)` where `v` is the array to set and `i` is the block index.
-An alternative syntax for this is `block_array[Block(i...)] = v` or
+A block can be set by  `block_array[Block(i...)] = v` or
 `block_array[Block.(i)...]`.
 
 ```jldoctest block_array
@@ -99,14 +98,14 @@ julia> block_array = BlockArray{Float64}(undef_blocks, [1,2], [2,2])
  #undef  #undef  │  #undef  #undef
  #undef  #undef  │  #undef  #undef
 
-julia> setblock!(block_array, rand(2,2), 2, 1)
+julia> block_array[Block(2,1)] = rand(2,2)
 2×2-blocked 3×4 BlockArray{Float64,2}:
  #undef      #undef      │  #undef  #undef
  ────────────────────────┼────────────────
    0.590845    0.566237  │  #undef  #undef
    0.766797    0.460085  │  #undef  #undef
 
-julia> block_array[Block(1, 1)] = [1 2];
+julia> block_array[Block(1),Block(1)] = [1 2];
 
 julia> block_array
 2×2-blocked 3×4 BlockArray{Float64,2}:
@@ -118,10 +117,15 @@ julia> block_array
 
 Note that this will "take ownership" of the passed in array, that is, no copy is made.
 
-A block can be retrieved with `getblock(block_array, i...)` or `block_array[Block(i...)]`:
+A block can be retrieved with `view(block_array, Block(i...))`, 
+or if a copy is desired, `block_array[Block(i...)]`:
 
 ```jldoctest block_array
-julia> block_array[Block(1, 1)]
+julia> view(block_array, Block(1, 1))
+1×2 Array{Float64,2}:
+ 1.0  2.0
+
+julia> block_array[Block(1, 1)] # makes a copy
 1×2 Array{Float64,2}:
  1.0  2.0
 
@@ -130,8 +134,6 @@ julia> block_array[Block(1), Block(1)]  # equivalent to above
  1.0  2.0
 ```
 
-Similarly to `setblock!` this does not copy the returned array.
-
 For setting and getting a single scalar element, the usual `setindex!` and `getindex` are available.
 
 ```jl
@@ -141,19 +143,25 @@ julia> block_array[1, 2]
 
 ## Views of blocks
 
-We can also view and modify views of blocks of `BlockArray` using the `view` syntax:
+To view and modify blocks of `BlockArray` use the `view` syntax.
 ```jldoctest
 julia> A = BlockArray(ones(6), 1:3);
 
 julia> view(A, Block(2))
-2-element view(::BlockArray{Float64,1,Array{Array{Float64,1},1},Tuple{BlockedUnitRange{Array{Int64,1}}}}, BlockSlice(Block(2),2:3)) with eltype Float64:
+2-element Array{Float64,1}:
  1.0
  1.0
 
 julia> view(A, Block(2)) .= [3,4]; A[Block(2)]
 2-element Array{Float64,1}:
  3.0
  4.0
+
+julia> view(A, Block.(1:2))
+3-element view(::BlockArray{Float64,1,Array{Array{Float64,1},1},Tuple{BlockedUnitRange{ArrayLayouts.RangeCumsum{Int64,UnitRange{Int64}}}}}, BlockSlice(Block{1,Int64}[Block(1), Block(2)],1:1:3)) with eltype Float64 with indices 1:1:3:
+ 1.0
+ 3.0
+ 4.0
 ```
 
 

docs/src/man/pseudoblockarrays.md

@@ -15,7 +15,7 @@ it can just return the wrapped array.
 
 When iteratively solving a set of equations with a gradient method the Jacobian typically has a block structure. It can be convenient
 to use a `PseudoBlockArray` to build up the Jacobian block by block and then pass the resulting matrix to
-a direct solver using `full`.
+a direct solver using `Matrix`.
 
 ## Creating PseudoBlockArrays
 
@@ -50,15 +50,16 @@ We can also any other user defined array type that supports `similar`.
 ## Setting and getting blocks and values
 
 Setting and getting blocks uses the same API as `BlockArrays`. The difference here is that setting a block will update the block in place and getting a block
-will extract a copy of the block and return it. For `PseudoBlockArrays` there is a mutating block getter called `getblock!` which updates a passed in array to avoid a copy:
+will extract a copy of the block and return it. Note to update a passed in array without allocating
+one can use views:
 
 ```jldoctest A
 julia> A = zeros(2,2)
 2×2 Array{Float64,2}:
  0.0  0.0
  0.0  0.0
 
-julia> getblock!(A, pseudo, 2, 1);
+julia> copyto!(A, view(pseudo, Block(2, 1)));
 
 julia> A
 2×2 Array{Float64,2}:
@@ -83,7 +84,7 @@ We can also view and modify views of blocks of `PseudoBlockArray` using the `vie
 julia> A = PseudoBlockArray(ones(6), 1:3);
 
 julia> view(A, Block(2))
-2-element view(::PseudoBlockArray{Float64,1,Array{Float64,1},Tuple{BlockedUnitRange{Array{Int64,1}}}}, BlockSlice(Block(2),2:3)) with eltype Float64:
+2-element view(::Array{Float64,1}, 2:3) with eltype Float64:
  1.0
  1.0
 
@@ -93,4 +94,4 @@ julia> view(A, Block(2)) .= [3,4]; A[Block(2)]
  4.0
 ```
 Note that, in memory, each block is in a BLAS-Level 3 compatible format, so
-that, in the future, algebra with blocks will be highly efficient.
+that algebra with blocks is highly efficient.

src/BlockArrays.jl

@@ -28,7 +28,7 @@ import Base: @propagate_inbounds, Array, to_indices, to_index,
             RangeIndex, Int, Integer, Number,
             +, -, *, /, \, min, max, isless, in, copy, copyto!, axes, @deprecate,
             BroadcastStyle, checkbounds, throw_boundserror, 
-            ones, zeros, intersect
+            ones, zeros, intersect, Slice
 using Base: ReshapedArray, dataids
 
 
@@ -41,10 +41,6 @@ import ArrayLayouts: _fill_lmul!, MatMulVecAdd, MatMulMatAdd, MatLmulVec, MatLdi
                         triangularlayout, triangulardata, _inv, _copyto!, axes_print_matrix_row,
                         colsupport, rowsupport, sub_materialize
 
-if !@isdefined(only)
-    using Compat: only
-end
-
 include("blockindices.jl")                        
 include("blockaxis.jl")
 include("abstractblockarray.jl")
@@ -60,4 +56,9 @@ include("show.jl")
 include("blockreduce.jl")
 include("blockdeque.jl")
 
+@deprecate getblock(A::AbstractBlockArray{T,N}, I::Vararg{Integer, N}) where {T,N} view(A, Block(I))
+@deprecate getblock!(X, A::AbstractBlockArray{T,N}, I::Vararg{Integer, N}) where {T,N} copyto!(X, view(A, Block(I)))
+@deprecate setblock!(A::AbstractBlockArray{T,N}, v, I::Vararg{Integer, N}) where {T,N} (A[Block(I...)] = v)
+
+
 end # module

src/abstractblockarray.jl

@@ -25,7 +25,7 @@ block2string(b, s) = string(join(map(string,b), '×'), "-blocked ", Base.dims2st
 _block_summary(a) = string(block2string(blocksize(a), size(a)), " ", typeof(a))
 Base.summary(a::AbstractBlockArray) = _block_summary(a)
 _show_typeof(io, a) = show(io, typeof(a))
-function _block_summary(io, a)    
+function _block_summary(io, a)
     print(io, block2string(blocksize(a), size(a)))
     print(io, ' ')
     _show_typeof(io, a)
@@ -37,98 +37,15 @@ Base.similar(a::AbstractBlockArray{T}) where {T}                             = s
 Base.similar(a::AbstractBlockArray, ::Type{T}) where {T}                     = similar(a, T, axes(a))
 Base.similar(a::AbstractBlockArray{T}, dims::Tuple) where {T}                = similar(a, T, dims)
 
+# If all we know is size, just return an Array which conforms to BlockArray interface
+Base.similar(::Type{<:AbstractBlockArray{T,N}}, dims::Dims) where {T,N} = similar(Array{T,N}, dims)
+
 Base.IndexStyle(::Type{<:AbstractBlockArray}) = IndexCartesian()
 
 # need to overload axes to return BlockAxis
 @inline size(block_array::AbstractBlockArray) = map(length, axes(block_array))
 @inline axes(block_array::AbstractBlockArray) = throw(error("axes for ", typeof(block_array), " is not implemented"))
 
-"""
-    getblock(A, inds...)
-
-Returns the block at blockindex `inds...`. An alternative syntax is `A[Block(inds...)].
-Throws a `BlockBoundsError` if this block is out of bounds.
-
-```jldoctest; setup = quote using BlockArrays end
-julia> v = Array(reshape(1:6, (2, 3)))
-2×3 Array{Int64,2}:
- 1  3  5
- 2  4  6
-
-julia> A = BlockArray(v, [1,1], [2,1])
-2×2-blocked 2×3 BlockArray{Int64,2}:
- 1  3  │  5
- ──────┼───
- 2  4  │  6
-
-julia> getblock(A, 2, 1)
-1×2 Array{Int64,2}:
- 2  4
-
-julia> A[Block(1, 2)]
-1×1 Array{Int64,2}:
- 5
-```
-"""
-function getblock(A::AbstractBlockArray{T,N}, ::Vararg{Integer, N}) where {T,N}
-    throw(error("getblock for ", typeof(A), " is not implemented"))
-end
-
-
-"""
-    getblock!(X, A, inds...)
-
-Stores the block at blockindex `inds` in `X` and returns it. Throws a `BlockBoundsError` if the
-attempted assigned block is out of bounds.
-
-```jldoctest; setup = quote using BlockArrays end
-julia> A = PseudoBlockArray(ones(2, 3), [1, 1], [2, 1])
-2×2-blocked 2×3 PseudoBlockArray{Float64,2}:
- 1.0  1.0  │  1.0
- ──────────┼─────
- 1.0  1.0  │  1.0
-
-julia> x = zeros(1, 2);
-
-julia> getblock!(x, A, 2, 1);
-
-julia> x
-1×2 Array{Float64,2}:
- 1.0  1.0
-```
-"""
-getblock!(X, A::AbstractBlockArray{T,N}, ::Vararg{Integer, N}) where {T,N} = throw(error("getblock! for ", typeof(A), " is not implemented"))
-
-@inline getblock!(X, A::AbstractBlockArray{T,N}, block::Block{N}) where {T,N}             = getblock!(X, A, block.n...)
-@inline getblock!(X, A::AbstractBlockVector, block::Block{1})                       = getblock!(X, A, block.n[1])
-@inline getblock!(X, A::AbstractBlockArray{T, N}, block::Vararg{Block{1}, N}) where {T,N} = getblock!(X, A, (Block(block).n)...)
-
-"""
-    setblock!(A, v, inds...)
-
-Stores the block `v` in the block at block index `inds` in `A`. An alternative syntax is `A[Block(inds...)] = v`.
-Throws a `BlockBoundsError` if this block is out of bounds.
-
-```jldoctest; setup = quote using BlockArrays end
-julia> A = PseudoBlockArray(zeros(2, 3), [1, 1], [2, 1]);
-
-julia> setblock!(A, [1 2], 1, 1);
-
-julia> A[Block(2, 1)] = [3 4];
-
-julia> A
-2×2-blocked 2×3 PseudoBlockArray{Float64,2}:
- 1.0  2.0  │  0.0
- ──────────┼─────
- 3.0  4.0  │  0.0
-```
-"""
-setblock!(A::AbstractBlockArray{T,N}, v, ::Vararg{Integer, N}) where {T,N} = throw(error("setblock! for ", typeof(A), " is not implemented"))
-
-@inline setblock!(A::AbstractBlockArray{T, N}, v, block::Block{N}) where {T,N}      = setblock!(A, v, block.n...)
-@inline setblock!(A::AbstractBlockVector, v, block::Block{1})                       = setblock!(A, v, block.n[1])
-@inline setblock!(A::AbstractBlockArray{T, N}, v, block::Vararg{Block{1}, N}) where {T,N} = setblock!(A, v, (Block(block).n)...)
-
 
 """
     BlockBoundsError([A], [inds...])
@@ -193,15 +110,27 @@ end
 end
 
 blockcheckbounds(A::AbstractArray{T, N}, i::Block{N}) where {T,N} = blockcheckbounds(A, i.n...)
+blockcheckbounds(A::AbstractArray{T, N}, i::Vararg{Block{1},N}) where {T,N} = blockcheckbounds(A, Int.(i)...)
+blockcheckbounds(A::AbstractVector{T}, i::Block{1}) where {T,N} = blockcheckbounds(A, Int(i))
 
 # Convert to @generated...
-@propagate_inbounds Base.getindex( block_arr::AbstractBlockArray{T, N}, block::Block{N}) where {T,N}       =  copy(getblock(block_arr, block.n...))
-@propagate_inbounds Base.setindex!(block_arr::AbstractBlockArray{T, N}, v, block::Block{N}) where {T,N}    =  setblock!(block_arr, v, block.n...)
-@propagate_inbounds Base.getindex( block_arr::AbstractBlockVector, block::Block{1})                  =  copy(getblock(block_arr, block.n[1]))
-@propagate_inbounds Base.setindex!(block_arr::AbstractBlockVector, v, block::Block{1})               =  setblock!(block_arr, v, block.n[1])
-@inline Base.getindex(block_arr::AbstractBlockArray{T,N}, block::Vararg{Block{1}, N}) where {T,N}     =  copy(getblock(block_arr, (Block(block).n)...))
-@inline Base.setindex!(block_arr::AbstractBlockArray{T,N}, v, block::Vararg{Block{1}, N}) where {T,N} =  setblock!(block_arr, v, (Block(block).n)...)
+@propagate_inbounds Base.getindex( block_arr::AbstractBlockArray{T,N}, block::Block{N}) where {T,N} = copy(view(block_arr, block))
+@propagate_inbounds Base.getindex( block_arr::AbstractBlockVector{T}, block::Block{1}) where {T} = copy(view(block_arr, block))
+@propagate_inbounds Base.setindex!(block_arr::AbstractBlockArray{T,N}, v, block::Block{N}) where {T,N} =
+    setindex!(block_arr, v, Block.(block.n)...)
+@inline @propagate_inbounds Base.getindex(block_arr::AbstractBlockArray{T,N}, block::Vararg{Block{1}, N}) where {T,N} =  block_arr[Block(block)]
+@inline @propagate_inbounds function Base.setindex!(block_arr::AbstractBlockArray{T,N}, v, block::Vararg{Block{1}, N}) where {T,N}
+    blockcheckbounds(block_arr, block...)
+    dest = view(block_arr, block...)
+    size(dest) == size(v) || throw(DimensionMismatch(string("tried to assign $(size(v)) array to $(size(dest)) block")))
+    copyto!(dest, v)
+    block_arr
+end
 
+viewblock(block_arr, block) = Base.invoke(view, Tuple{AbstractArray, Any}, block_arr, block)
+@inline Base.view(block_arr::AbstractBlockArray{<:Any,N}, block::Block{N}) where N = viewblock(block_arr, block)
+@inline Base.view(block_arr::AbstractBlockVector, block::Block{1}) = viewblock(block_arr, block)
+@inline @propagate_inbounds Base.view(block_arr::AbstractBlockArray, block::Block{1}...) = view(block_arr, Block(block))
 """
     eachblock(A::AbstractBlockArray)