Merge branch 'master' into generalize_blockedunitrange

Author: mtfishman

.github/workflows/downstream.yml

@@ -31,6 +31,7 @@ jobs:
         os: [ubuntu-latest]
         package:
           - {repo: BlockBandedMatrices.jl, group: JuliaLinearAlgebra}
+          - {repo: ApproxFunBase.jl, group: JuliaApproximation}
           # - {repo: InfiniteLinearAlgebra.jl, group: JuliaLinearAlgebra}
 
     steps:

Project.toml

@@ -1,6 +1,6 @@
 name = "BlockArrays"
 uuid = "8e7c35d0-a365-5155-bbbb-fb81a777f24e"
-version = "0.16.39"
+version = "0.16.40"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -10,7 +10,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 [compat]
 Aqua = "0.8"
 ArrayLayouts = "1.0.8"
-Documenter = "0.27"
+Documenter = "1"
 FillArrays = "1"
 LinearAlgebra = "1.6"
 OffsetArrays = "1"

docs/src/man/blockarrays.md

@@ -2,32 +2,32 @@
 
 ```@meta
 DocTestSetup = quote
-    using BlockArrays
-    using Random
-    Random.seed!(1234)
+    using BlockArrays, SparseArrays
 end
 ```
 
 ## Creating `BlockArray`s from an array
 
 An `AbstractArray` can be repacked into a `BlockArray` with `BlockArray(array, block_sizes...)`.  The block sizes are each an `AbstractVector{Int}` which determines the size of the blocks in that dimension (so the sum of `block_sizes` in every dimension must match the size of `array` in that dimension).
 
-```julia
-julia> BlockArray(rand(4, 4), [2,2], [1,1,2])
-2×3-blocked 4×4 BlockMatrix{Float64}:
- 0.70393   │  0.568703  │  0.0137366  0.953038
- 0.24957   │  0.145924  │  0.884324   0.134155
- ──────────┼────────────┼─────────────────────
- 0.408133  │  0.707723  │  0.467458   0.326718
- 0.844314  │  0.794279  │  0.0421491  0.683791
-
-julia> block_array_sparse = BlockArray(sprand(4, 5, 0.7), [1,3], [2,3])
-2×2-blocked 4×5 BlockMatrix{Float64, Matrix{SparseMatrixCSC{Float64, Int64}}, Tuple{BlockedUnitRange{Int64, Vector{Int64}}, BlockedUnitRange{Int64, Vector{Int64}}}}:
- 0.0341601  0.374187  │  0.0118196  0.299058  0.0
- ---------------------┼-------------------------------
- 0.0945445  0.931115  │  0.0460428  0.0       0.0
- 0.314926   0.438939  │  0.496169   0.0       0.0
- 0.12781    0.246862  │  0.732      0.449182  0.875096
+```jldoctest
+julia> BlockArray(Array(reshape(1:16, 4, 4)), [2,2], [1,1,2])
+2×3-blocked 4×4 BlockMatrix{Int64}:
+ 1  │  5  │   9  13
+ 2  │  6  │  10  14
+ ───┼─────┼────────
+ 3  │  7  │  11  15
+ 4  │  8  │  12  16
+
+julia> S = spzeros(4,5); S[1,2] = S[4,3] = 1;
+
+julia> block_array_sparse = BlockArray(S, [1,3], [2,3])
+2×2-blocked 4×5 BlockMatrix{Float64, Matrix{SparseMatrixCSC{Float64, Int64}}, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}}:
+  ⋅   1.0  │   ⋅    ⋅    ⋅ 
+ ──────────┼───────────────
+  ⋅    ⋅   │   ⋅    ⋅    ⋅ 
+  ⋅    ⋅   │   ⋅    ⋅    ⋅ 
+  ⋅    ⋅   │  1.0   ⋅    ⋅ 
 ```
 
 
@@ -65,11 +65,11 @@ julia> BlockArray{Float32}(undef_blocks, [1,2], [3,2])
 
 The `block_type` should be an array type.  It specifies the internal block type, which defaults to an `Array` of the according dimension.  We can also use a `SparseVector` or any other user defined array type:
 
-```julia
+```jldoctest
 julia> BlockArray(undef_blocks, SparseVector{Float64, Int}, [1,2])
 2-blocked 3-element BlockVector{Float64, Vector{SparseVector{Float64, Int64}}, Tuple{BlockedUnitRange{Int64, Vector{Int64}}}}:
  #undef
- ------
+ ──────
  #undef
  #undef
 ```
@@ -78,10 +78,9 @@ julia> BlockArray(undef_blocks, SparseVector{Float64, Int}, [1,2])
 
     Note that accessing an undefined block will throw an "access to undefined reference"-error!  If you create an array with undefined blocks, you _have_ to [initialize it block-wise](@ref setting_and_getting)); whole-array functions like `fill!` will not work:
 
-    ```julia
+    ```jldoctest
     julia> fill!(BlockArray{Float32}(undef_blocks, [1,2], [3,2]), 0)
     ERROR: UndefRefError: access to undefined reference
-    …
     ```
 
 ## [Setting and getting blocks and values](@id setting_and_getting)
@@ -164,32 +163,35 @@ julia> view(A, Block.(1:2))
 
 An array can be repacked into a `BlockArray` with `BlockArray(array, block_sizes...)`:
 
-```jl
-julia> block_array_sparse = BlockArray(sprand(4, 5, 0.7), [1,3], [2,3])
-2×2-blocked 4×5 BlockArray{Float64, 2, Matrix{SparseMatrixCSC{Float64, Int64}}, Tuple{BlockedUnitRange{Int64, Vector{Int64}}, BlockedUnitRange{Int64, Vector{Int64}}}}:
- 0.0341601  0.374187  │  0.0118196  0.299058  0.0
- ---------------------┼-------------------------------
- 0.0945445  0.931115  │  0.0460428  0.0       0.0
- 0.314926   0.438939  │  0.496169   0.0       0.0
- 0.12781    0.246862  │  0.732      0.449182  0.875096
+```jldoctest repack
+julia> S = spzeros(4,5); S[1,2] = S[4,3] = 1;
+
+julia> block_array_sparse = BlockArray(S, [1,3], [2,3])
+2×2-blocked 4×5 BlockMatrix{Float64, Matrix{SparseMatrixCSC{Float64, Int64}}, Tuple{BlockedUnitRange{Vector{Int64}}, BlockedUnitRange{Vector{Int64}}}}:
+  ⋅   1.0  │   ⋅    ⋅    ⋅ 
+ ──────────┼───────────────
+  ⋅    ⋅   │   ⋅    ⋅    ⋅ 
+  ⋅    ⋅   │   ⋅    ⋅    ⋅ 
+  ⋅    ⋅   │  1.0   ⋅    ⋅ 
 ```
 
-To get back the underlying array use `Array`:
+To get back the underlying sparse array, use `sparse`:
 
-```jl
+```jldoctest repack
+julia> sparse(block_array_sparse)
+4×5 SparseMatrixCSC{Float64, Int64} with 2 stored entries:
+  ⋅   1.0   ⋅    ⋅    ⋅ 
+  ⋅    ⋅    ⋅    ⋅    ⋅ 
+  ⋅    ⋅    ⋅    ⋅    ⋅ 
+  ⋅    ⋅   1.0   ⋅    ⋅ 
+```
+
+To get a dense array, use `Array`:
+```jldoctest repack
 julia> Array(block_array_sparse)
-4×5 SparseMatrixCSC{Float64,Int64} with 13 stored entries:
-  [1, 1]  =  0.30006
-  [2, 1]  =  0.451742
-  [3, 1]  =  0.243174
-  [4, 1]  =  0.156468
-  [1, 2]  =  0.94057
-  [3, 2]  =  0.544175
-  [4, 2]  =  0.598345
-  [3, 3]  =  0.737486
-  [4, 3]  =  0.929512
-  [1, 4]  =  0.539601
-  [3, 4]  =  0.757658
-  [4, 4]  =  0.44709
-  [2, 5]  =  0.514679
+4×5 Matrix{Float64}:
+ 0.0  1.0  0.0  0.0  0.0
+ 0.0  0.0  0.0  0.0  0.0
+ 0.0  0.0  0.0  0.0  0.0
+ 0.0  0.0  1.0  0.0  0.0
 ```

src/blockaxis.jl

@@ -74,21 +74,21 @@ _blocklengths2blocklasts(blocks::Fill) = cumsum(blocks)
 _blocklengths2blocklasts(blocks::Ones) = cumsum(blocks)
 @inline blockedrange(blocks::Union{Tuple,AbstractVector}) = _BlockedUnitRange(_blocklengths2blocklasts(blocks))
 
-@inline blockfirsts(a::BlockedUnitRange) = [a.first; @views(a.lasts[1:end-1]) .+ 1]
+@inline blockfirsts(a::BlockedUnitRange) = [first(a); @views(blocklasts(a)[1:end-1]) .+ 1]
 # optimize common cases
-@inline function blockfirsts(a::BlockedUnitRange{<:Integer, <:Union{Vector, RangeCumsum{<:Any, <:UnitRange}}})
-    v = Vector{eltype(a)}(undef, length(a.lasts))
-    v[1] = a.first
-    v[2:end] .= @views(a.lasts[oneto(end-1)]) .+ 1
+@inline function blockfirsts(a::BlockedUnitRange{<:Union{Vector, RangeCumsum{<:Any, <:UnitRange}}})
+    v = Vector{eltype(a)}(undef, length(blocklasts(a)))
+    v[1] = first(a)
+    v[2:end] .= @views(blocklasts(a)[oneto(end-1)]) .+ 1
     return v
 end
 @inline blocklasts(a::BlockedUnitRange) = a.lasts
 
 _diff(a::AbstractVector) = diff(a)
 _diff(a::Tuple) = diff(collect(a))
-@inline blocklengths(a::BlockedUnitRange) = isempty(a.lasts) ? [_diff(a.lasts);] : [first(a.lasts)-a.first+one(eltype(a)); _diff(a.lasts)]
+@inline blocklengths(a::BlockedUnitRange) = isempty(blocklasts(a)) ? [_diff(blocklasts(a));] : [first(blocklasts(a))-first(a)+one(eltype(a)); _diff(blocklasts(a))]
 
-length(a::BlockedUnitRange) = isempty(a.lasts) ? zero(eltype(a)) : Integer(last(a.lasts)-a.first+one(eltype(a)))
+length(a::BlockedUnitRange) = isempty(blocklasts(a)) ? zero(eltype(a)) : Integer(last(blocklasts(a))-first(a)+one(eltype(a)))
 
 """
     blockisequal(a::AbstractUnitRange{Int}, b::AbstractUnitRange{Int})
@@ -142,7 +142,7 @@ Base.unitrange(b::BlockedUnitRange) = first(b):last(b)
 Base.promote_rule(::Type{BlockedUnitRange{T, CS}}, ::Type{Base.OneTo{Int}}) where {T, CS} = UnitRange{T}
 
 """
-    blockaxes(A)
+    blockaxes(A::AbstractArray)
 
 Return the tuple of valid block indices for array `A`.
 
@@ -169,19 +169,19 @@ julia> blockaxes(B)
 (BlockRange(Base.OneTo(2)), BlockRange(Base.OneTo(3)))
 ```
 """
-blockaxes(b::BlockedUnitRange) = _blockaxes(b.lasts)
+blockaxes(b::BlockedUnitRange) = _blockaxes(blocklasts(b))
 _blockaxes(b::AbstractVector) = (Block.(axes(b,1)),)
 _blockaxes(b::Tuple) = (Block.(Base.OneTo(length(b))),)
 blockaxes(b) = blockaxes.(axes(b), 1)
 
 """
-    blockaxes(A, d)
+    blockaxes(A::AbstractArray, d::Int)
 
 Return the valid range of block indices for array `A` along dimension `d`.
 
 # Examples
 ```jldoctest
-julia> A = BlockArray([1,2,3],[2,1])
+julia> A = BlockArray([1,2,3], [2,1])
 2-blocked 3-element BlockVector{Int64}:
  1
  2
@@ -202,8 +202,8 @@ julia> blockaxes(A,1) |> collect
 end
 
 """
-    blocksize(A)
-    blocksize(A,i)
+    blocksize(A::AbstractArray)
+    blocksize(A::AbstractArray, i::Int)
 
 Return the tuple of the number of blocks along each
 dimension. See also size and blocksizes.
@@ -229,8 +229,8 @@ blocksize(A,i) = length(blockaxes(A,i))
 @inline blocklength(t) = prod(blocksize(t))
 
 """
-    blocksizes(A)
-    blocksizes(A,i)
+    blocksizes(A::AbstractArray)
+    blocksizes(A::AbstractArray, i::Int)
 
 Return the tuple of the sizes of blocks along each
 dimension. See also size and blocksize.
@@ -352,8 +352,8 @@ Return the first index of each block of `a`.
 
 # Examples
 ```jldoctest
-julia> b = blockedrange(1:3)
-3-blocked 6-element BlockedUnitRange{Int64, ArrayLayouts.RangeCumsum{Int64, UnitRange{Int64}}}:
+julia> b = blockedrange([1,2,3])
+3-blocked 6-element BlockedUnitRange{Int64, Vector{Int64}}:
  1
  ─
  2
@@ -378,8 +378,8 @@ Return the last index of each block of `a`.
 
 # Examples
 ```jldoctest
-julia> b = blockedrange(1:3)
-3-blocked 6-element BlockedUnitRange{Int64, ArrayLayouts.RangeCumsum{Int64, UnitRange{Int64}}}:
+julia> b = blockedrange([1,2,3])
+3-blocked 6-element BlockedUnitRange{Int64, Vector{Int64}}:
  1
  ─
  2
@@ -390,7 +390,7 @@ julia> b = blockedrange(1:3)
  6
 
 julia> blocklasts(b)
-3-element ArrayLayouts.RangeCumsum{Int64, UnitRange{Int64}}:
+3-element Vector{Int64}:
  1
  3
  6
@@ -404,8 +404,8 @@ Return the length of each block of `a`.
 
 # Examples
 ```jldoctest
-julia> b = blockedrange(1:3)
-3-blocked 6-element BlockedUnitRange{Int64, ArrayLayouts.RangeCumsum{Int64, UnitRange{Int64}}}:
+julia> b = blockedrange([1,2,3])
+3-blocked 6-element BlockedUnitRange{Int64, Vector{Int64}}:
  1
  ─
  2
@@ -453,24 +453,24 @@ Base.BroadcastStyle(::Type{BlockedUnitRange{T, R}}) where {T, R} = Base.Broadcas
 
 _blocklengths2blocklasts(blocks::AbstractRange) = RangeCumsum(blocks)
 function blockfirsts(a::BlockedUnitRange{<:Integer, <:Base.OneTo})
-    a.first == 1 || error("Offset axes not supported")
-    Base.OneTo{eltype(a)}(length(a.lasts))
+    first(a) == 1 || error("Offset axes not supported")
+    Base.OneTo{eltype(a)}(length(blocklasts(a)))
 end
 function blocklengths(a::BlockedUnitRange{<:Integer, <:Base.OneTo})
-    a.first == 1 || error("Offset axes not supported")
-    Ones{eltype(a)}(length(a.lasts))
+    first(a) == 1 || error("Offset axes not supported")
+    Ones{eltype(a)}(length(blocklasts(a)))
 end
 function blockfirsts(a::BlockedUnitRange{<:Integer, <:AbstractRange})
-    st = step(a.lasts)
-    a.first == 1 || error("Offset axes not supported")
-    @assert first(a.lasts)-a.first+1 == st
-    range(one(eltype(a)); step=st, length=eltype(a)(length(a.lasts)))
+    st = step(blocklasts(a))
+    first(a) == 1 || error("Offset axes not supported")
+    @assert first(blocklasts(a))-first(a)+1 == st
+    range(one(eltype(a)); step=st, length=eltype(a)(length(blocklasts(a))))
 end
 function blocklengths(a::BlockedUnitRange{<:Integer, <:AbstractRange})
-    st = step(a.lasts)
-    a.first == 1 || error("Offset axes not supported")
-    @assert first(a.lasts)-a.first+1 == st
-    Fill(st,length(a.lasts))
+    st = step(blocklasts(a))
+    first(a) == 1 || error("Offset axes not supported")
+    @assert first(blocklasts(a))-first(a)+1 == st
+    Fill(st,length(blocklasts(a)))
 end
 
 

src/blockbroadcast.jl

@@ -67,7 +67,7 @@ julia> import BlockArrays: SubBlockIterator, BlockIndexRange
 
 julia> A = BlockArray(1:6, 1:3);
 
-julia> subblock_lasts = axes(A, 1).lasts;
+julia> subblock_lasts = blocklasts(axes(A, 1));
 
 julia> @assert subblock_lasts == [1, 3, 6];