materialize as a block array whenever axes are blocked (#133)

* materialize as a block array whenever axes are blocked

* exp,log,sqrt overloads for PseudoBlockMatrix

* Update blocklinalg.jl

* add tests, support Colon

* Increase coverage

Author: dlfivefifty

src/blockaxis.jl

@@ -228,3 +228,15 @@ blocklengths(a::AbstractUnitRange) = blocklasts(a) .- blockfirsts(a) .+ 1
 
 Base.summary(a::BlockedUnitRange) = _block_summary(a)
 Base.summary(io::IO, a::BlockedUnitRange) =  _block_summary(io, a)
+
+
+###
+# Slice{<:BlockedUnitRange}
+###
+
+Base.axes(S::Base.Slice{<:BlockedUnitRange}) = (S.indices,)
+Base.unsafe_indices(S::Base.Slice{<:BlockedUnitRange}) = (S.indices,)
+Base.axes1(S::Base.Slice{<:BlockedUnitRange}) = S.indices
+blockaxes(S::Base.Slice) = blockaxes(S.indices)
+getindex(S::Base.Slice, b::Block{1}) = S.indices[b]
+getindex(S::Base.Slice, b::BlockRange{1}) = S.indices[b]

src/blocklinalg.jl

@@ -61,7 +61,13 @@ sublayout(BL::BlockLayout{MLAY,BLAY}, ::Type{<:Tuple{<:BlockSlice{BlockRange{1,T
 
 # materialize views, used for `getindex`
 sub_materialize(::AbstractBlockLayout, V, _) = BlockArray(V)
-sub_materialize(::AbstractStridedLayout, V, ::Tuple{<:BlockedUnitRange}) = PseudoBlockArray(V)
+sub_materialize(::AbstractBlockLayout, V, ::Tuple{<:BlockedUnitRange}) = BlockArray(V)
+sub_materialize(::AbstractBlockLayout, V, ::Tuple{<:BlockedUnitRange,<:BlockedUnitRange}) = BlockArray(V)
+# if it's not a block layout, best to use PseudoBlockArray
+sub_materialize(_, V, ::Tuple{<:BlockedUnitRange}) = PseudoBlockArray(V)
+sub_materialize(_, V, ::Tuple{<:BlockedUnitRange,<:BlockedUnitRange}) = PseudoBlockArray(V)
+sub_materialize(_, V, ::Tuple{<:AbstractUnitRange,<:BlockedUnitRange}) = PseudoBlockArray(V)
+sub_materialize(_, V, ::Tuple{<:BlockedUnitRange,<:AbstractUnitRange}) = PseudoBlockArray(V)
 
 conjlayout(::Type{T}, ::BlockLayout{MLAY,BLAY}) where {T<:Complex,MLAY,BLAY} = BlockLayout{MLAY,typeof(conjlayout(T,BLAY()))}()
 conjlayout(::Type{T}, ::BlockLayout{MLAY,BLAY}) where {T<:Real,MLAY,BLAY} = BlockLayout{MLAY,BLAY}()
@@ -279,4 +285,10 @@ for UNIT in ('U', 'N')
 end
 
 # For now, use PseudoBlockArray
-_inv(::AbstractBlockLayout, axes, A) = BlockArray(inv(PseudoBlockArray(A)))
+_inv(::AbstractBlockLayout, axes, A) = BlockArray(inv(PseudoBlockArray(A)))
+
+for op in (:exp, :log, :sqrt)
+    @eval begin
+       $op(A::PseudoBlockMatrix) = PseudoBlockMatrix($op(A.blocks), axes(A))
+    end
+end

test/test_blockindices.jl

@@ -230,6 +230,8 @@ end
         @test axes(b) == Base.unsafe_indices(b) == (b,)
         @test Base.dataids(b) == Base.dataids(blocklasts(b))
         @test_throws ArgumentError BlockedUnitRange(b)
+
+        @test summary(b) == "3-blocked 6-element BlockedUnitRange{Array{$Int,1}}"
     end
 
     @testset "OneTo interface" begin
@@ -261,6 +263,15 @@ end
         @test checkindex(Bool, b, Block(3))
         @test !checkindex(Bool, b, Block(4))
     end
+
+    @testset "Slice" begin
+        b = blockedrange([1,2,3])
+        S = Base.Slice(b)
+        @test blockaxes(S) == blockaxes(b)
+        @test S[Block(2)] == 2:3
+        @test S[Block.(1:2)] == 1:3
+        @test axes(S) == axes(b)
+    end
 end
 
 @testset "BlockSlice" begin

test/test_blockviews.jl

@@ -221,4 +221,41 @@ using BlockArrays, Test, Base64
 
         @test BlockArrays.hasmatchingblocks(view(B,Block.(3:3),Block.(2:2)))
     end
+
+    @testset "sub_materialize cases" begin
+        a = BlockArray(randn(6), 1:3)
+        b = PseudoBlockArray(randn(6), 1:3)
+        @test a[Block.(1:2)] isa BlockArray
+        @test b[Block.(1:2)] isa PseudoBlockArray
+        @test a[1:3] isa Array
+        @test b[1:3] isa Array
+        A = BlockArray(randn(6,6), 1:3, fill(3,2))
+        B = PseudoBlockArray(randn(6,6), 1:3, fill(3,2))
+        @test A[Block.(1:2),Block.(1:2)] isa BlockArray
+        @test B[Block.(1:2),Block.(1:2)] isa PseudoBlockArray
+        @test A[Block.(1:2),1:3] isa PseudoBlockArray
+        @test B[Block.(1:2),1:3] isa PseudoBlockArray
+        @test A[1:3,Block.(1:2)] isa PseudoBlockArray
+        @test B[1:3,Block.(1:2)] isa PseudoBlockArray
+        @test A[Block.(1:2),:] isa PseudoBlockArray
+        @test B[Block.(1:2),:] isa PseudoBlockArray
+        @test blockisequal(axes(A,2),axes(A[Block.(1:2),:],2))
+        @test blockisequal(axes(B,2),axes(B[Block.(1:2),:],2))
+        @test A[:,Block.(1:2)] isa PseudoBlockArray
+        @test B[:,Block.(1:2)] isa PseudoBlockArray
+        @test blockisequal(axes(A,1),axes(A[:,Block.(1:2)],1))
+        @test blockisequal(axes(B,1),axes(B[:,Block.(1:2)],1))
+        @test A[:,:] isa PseudoBlockArray
+        @test B[:,:] isa PseudoBlockArray
+        @test blockisequal(axes(A),axes(A[:,:]))
+        @test blockisequal(axes(B),axes(B[:,:]))
+        @test A[1:3,1:3] isa Array
+        @test B[1:3,1:3] isa Array
+        A = BlockArray(randn(6,6,6), 1:3, fill(3,2),1:3)
+        B = PseudoBlockArray(randn(6,6,6), 1:3, fill(3,2),1:3)
+        @test A[Block.(1:2),Block.(1:2),Block.(1:2)] isa BlockArray
+        @test B[Block.(1:2),Block.(1:2),Block.(1:2)] isa PseudoBlockArray
+        @test A[1:3,Block.(1:2),1:3] isa BlockArray
+        @test B[1:3,Block.(1:2),1:3] isa PseudoBlockArray
+    end
 end