Skip to content
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

Fix multiple problems with views-of-views #10505

Merged
merged 2 commits into from
Mar 22, 2015
Merged
Changes from 1 commit
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
Next Next commit
Fix multiple problems with views-of-views
Due to what was probably a copy/paste error, the tests
were grabbing a global variable A (rather than an argument)
that happened to be a plain array. So we weren't really testing
views-of-views.

Unsurprisingly, fixing this problem uncovered bugs.
Many "simple" things worked, but sufficiently complicated
constructs like
   A = rand(20, 20, 20)
   B = slice(A, 1:5, 6, 6:9)
   C = sub(B, :)
(specifically, those that use linear indexing in creation
of a view-of-view) certainly had problems.

Drat. Well, hopefully this has been rare up until now;
it seems rather likely that folks would have gotten BoundsErrors if
anyone had actually been using this. But we start returning views from
getindex it will become common, so better that it was caught now.

In debugging and fixing the problems, I took the opportunity to
do some cleanup to make this more maintainable, specifically:
adopting a consistent naming convention, improving many variable
names, and adding comments.
timholy committed Mar 13, 2015
commit 9a86698f7c120914f84203664a7996d0ba34e147
64 changes: 37 additions & 27 deletions base/multidimensional.jl
Original file line number Diff line number Diff line change
@@ -345,61 +345,71 @@ end
# the corresponding cartesian index into the parent, and then uses
# dims to convert back to a linear index into the parent array.
#
# However, a common case is linindex::UnitRange.
# Since div is slow and in(j::Int, linindex::UnitRange) is fast,
# However, a common case is linindex::Range.
# Since div is slow and in(j::Int, linindex::Range) is fast,
# it can be much faster to generate all possibilities and
# then test whether the corresponding linear index is in linindex.
# One exception occurs when only a small subset of the total
# is desired, in which case we fall back to the div-based algorithm.
stagedfunction merge_indexes(V, indexes::NTuple, dims::Dims, linindex::UnitRange{Int})
N = length(indexes)
#stagedfunction merge_indexes{T<:Integer}(V, parentindexes::NTuple, parentdims::Dims, linindex::Union(Colon,Range{T}), lindim)
stagedfunction merge_indexes_in{TT}(V, parentindexes::TT, parentdims::Dims, linindex, lindim)
N = length(parentindexes) # number of parent axes we're merging
N > 0 || throw(ArgumentError("cannot merge empty indexes"))
lengthexpr = linindex == Colon ? (:(prod(size(V)[lindim:end]))) : (:(length(linindex)))
L = symbol(string("Istride_", N+1)) # length of V's trailing dimensions
quote
n = length(linindex)
Base.Cartesian.@nexprs $N d->(I_d = indexes[d])
L = 1
dimoffset = ndims(V.parent) - length(dims)
Base.Cartesian.@nexprs $N d->(L *= dimsize(V.parent, d+dimoffset, I_d))
if n < 0.1L # this has not been carefully tuned
return merge_indexes_div(V, indexes, dims, linindex)
n = $lengthexpr
Base.Cartesian.@nexprs $N d->(I_d = parentindexes[d])
pdimoffset = ndims(V.parent) - length(parentdims)
Istride_1 = 1 # parentindexes strides
Base.Cartesian.@nexprs $N d->(Istride_{d+1} = Istride_d*dimsize(V.parent, d+pdimoffset, I_d))
Istridet = Base.Cartesian.@ntuple $(N+1) d->Istride_d
if n < 0.1*$L # this has not been carefully tuned
return merge_indexes_div(V, parentindexes, parentdims, linindex, lindim)
end
Pstride_1 = 1 # parent strides
Base.Cartesian.@nexprs $(N-1) d->(Pstride_{d+1} = Pstride_d*dims[d])
Istride_1 = 1 # indexes strides
Base.Cartesian.@nexprs $(N-1) d->(Istride_{d+1} = Istride_d*dimsize(V, d+dimoffset, I_d))
Base.Cartesian.@nexprs $N d->(counter_d = 1) # counter_0 is a linear index into indexes
Base.Cartesian.@nexprs $(N-1) d->(Pstride_{d+1} = Pstride_d*parentdims[d])
Base.Cartesian.@nexprs $N d->(counter_d = 1) # counter_0 is a linear index into parentindexes
Base.Cartesian.@nexprs $N d->(offset_d = 1) # offset_0 is a linear index into parent
k = 0
index = Array(Int, n)
Base.Cartesian.@nloops $N i d->(1:dimsize(V, d+dimoffset, I_d)) d->(offset_{d-1} = offset_d + (I_d[i_d]-1)*Pstride_d; counter_{d-1} = counter_d + (i_d-1)*Istride_d) begin
Base.Cartesian.@nloops $N i d->(1:dimsize(V.parent, d+pdimoffset, I_d)) d->(offset_{d-1} = offset_d + (I_d[i_d]-1)*Pstride_d; counter_{d-1} = counter_d + (i_d-1)*Istride_d) begin
if in(counter_0, linindex)
index[k+=1] = offset_0
end
end
index
end
end
merge_indexes(V, indexes::NTuple, dims::Dims, linindex) = merge_indexes_div(V, indexes, dims, linindex)

# HACK: dispatch seemingly wasn't working properly
function merge_indexes(V, parentindexes::NTuple, parentdims::Dims, linindex, lindim)
if isa(linindex, Colon) || isa(linindex, Range)
return merge_indexes_in(V, parentindexes, parentdims, linindex, lindim)
end
merge_indexes_div(V, parentindexes, parentdims, linindex, lindim)
end

# This could be written as a regular function, but performance
# will be better using Cartesian macros to avoid the heap and
# an extra loop.
stagedfunction merge_indexes_div(V, indexes::NTuple, dims::Dims, linindex)
N = length(indexes)
stagedfunction merge_indexes_div{TT}(V, parentindexes::TT, parentdims::Dims, linindex, lindim)
N = length(parentindexes)
N > 0 || throw(ArgumentError("cannot merge empty indexes"))
Istride_N = symbol("Istride_$N")
lengthexpr = :(length(linindex))
quote
Base.Cartesian.@nexprs $N d->(I_d = indexes[d])
Base.Cartesian.@nexprs $N d->(I_d = parentindexes[d])
Pstride_1 = 1 # parent strides
Base.Cartesian.@nexprs $(N-1) d->(Pstride_{d+1} = Pstride_d*dims[d])
Istride_1 = 1 # indexes strides
dimoffset = ndims(V.parent) - length(dims)
Base.Cartesian.@nexprs $(N-1) d->(Istride_{d+1} = Istride_d*dimsize(V.parent, d+dimoffset, I_d))
n = length(linindex)
L = $(Istride_N) * dimsize(V.parent, $N+dimoffset, indexes[end])
Base.Cartesian.@nexprs $(N-1) d->(Pstride_{d+1} = Pstride_d*parentdims[d])
Istride_1 = 1 # parentindexes strides
pdimoffset = ndims(V.parent) - length(parentdims)
Base.Cartesian.@nexprs $(N-1) d->(Istride_{d+1} = Istride_d*dimsize(V.parent, d+pdimoffset, I_d))
n = $lengthexpr
L = $(Istride_N) * dimsize(V.parent, $N+pdimoffset, parentindexes[end])
index = Array(Int, n)
for i = 1:n
k = linindex[i] # k is the indexes-centered linear index
k = linindex[i] # k is the parentindexes-centered linear index
1 <= k <= L || throw(BoundsError())
k -= 1
j = 0 # j will be the new parent-centered linear index
245 changes: 150 additions & 95 deletions base/subarray.jl
Original file line number Diff line number Diff line change
@@ -3,7 +3,7 @@ typealias ViewIndex Union(Int, NonSliceIndex)
typealias RangeIndex Union(Int, Range{Int}, UnitRange{Int}, Colon)

# LD is the last dimension up through which this object has efficient
# linear indexing. If LD==N, then the object itself has efficient
# linear indexing. If LD==length(I), then the object itself has efficient
# linear indexing.
immutable SubArray{T,N,P<:AbstractArray,I<:(ViewIndex...),LD} <: AbstractArray{T,N}
parent::P
@@ -40,61 +40,79 @@ function _slice(A, I)
slice_unsafe(A, I)
end

stagedfunction slice_unsafe{T,NP,IndTypes}(A::AbstractArray{T,NP}, I::IndTypes)
# The most complicated part of this is matching the axes between the
# input index tuples (denoted by J), the index tuples that get stored
# in the view (denoted by I), and the overall dimensionality of the
# view.
# The complexities increase when you create a view-of-a-view, because
# then there is also the index tuple of the parent view (denoted IV)
# to consider.
#
# Examples:
# S1 = sub(A::Matrix, 2, 3:5) ndims(S1) == length(I) == length(J) == 2
# S2 = slice(A::Matrix, 2, 3:5) ndims(S2) == 1, length(I) == length(J) == 2
# S3 = sub(A::Matrix, 4:17) ndims(S3) == length(I) == length(J) == 1
# S4 = sub(S2, 1:2) ndims(S4) == length(J) == 1, length(I) == 2
# S3 addresses the trailing dimensions of the parent by linear indexing.
# For S4, J[1] corresponds to I[2], because of the slice along
# dimension 1 in S2

stagedfunction slice_unsafe{T,NP,IndTypes}(A::AbstractArray{T,NP}, J::IndTypes)
N = 0
sizeexprs = Array(Any, 0)
for k = 1:length(I)
i = I[k]
if !(i <: Real)
for Jindex = 1:length(J)
j = J[Jindex]
if !(j <: Real)
N += 1
push!(sizeexprs, dimsizeexpr(I[k], k, length(I), :A, :I))
push!(sizeexprs, dimsizeexpr(j, Jindex, length(J), :A, :J))
end
end
dims = :(tuple($(sizeexprs...)))
LD = subarray_linearindexing_dim(A, I)
strideexpr = stride1expr(A, I, :A, :I, LD)
exfirst = first_index_expr(:A, :I, length(I))
LD = subarray_linearindexing_dim(A, J)
strideexpr = stride1expr(A, J, :A, :J, LD)
exfirst = first_index_expr(:A, :J, length(J))
quote
$exfirst
SubArray{$T,$N,$A,$I,$LD}(A, I, $dims, f, $strideexpr)
SubArray{$T,$N,$A,$J,$LD}(A, J, $dims, f, $strideexpr)
end
end

# Conventional style (drop trailing singleton dimensions, keep any
# other singletons)
# other singletons by converting them to ranges, e.g., 3:3)
sub(A::AbstractArray, I::ViewIndex...) = _sub(A, I)
sub(A::AbstractArray, I::(ViewIndex...)) = _sub(A, I)
function _sub(A, I)
checkbounds(A, I...)
sub_unsafe(A, I)
end

stagedfunction sub_unsafe{T,NP,IndTypes}(A::AbstractArray{T,NP}, I::IndTypes)
stagedfunction sub_unsafe{T,NP,IndTypes}(A::AbstractArray{T,NP}, J::IndTypes)
sizeexprs = Array(Any, 0)
Itypes = Array(Any, 0)
Iexprs = Array(Any, 0)
N = length(I)
while N > 0 && I[N] <: Real
N = length(J)
while N > 0 && J[N] <: Real
N -= 1
end
for k = 1:length(I)
if k <= N
push!(sizeexprs, dimsizeexpr(I[k], k, length(I), :A, :I))
for Jindex = 1:length(J)
j = J[Jindex]
if Jindex <= N
push!(sizeexprs, dimsizeexpr(j, Jindex, length(J), :A, :J))
end
if k < N && I[k] <: Real
if Jindex < N && j <: Real
push!(Itypes, UnitRange{Int})
push!(Iexprs, :(Int(I[$k]):Int(I[$k])))
push!(Iexprs, :(Int(J[$Jindex]):Int(J[$Jindex])))
else
push!(Itypes, I[k])
push!(Iexprs, :(I[$k]))
push!(Itypes, j)
push!(Iexprs, :(J[$Jindex]))
end
end
dims = :(tuple($(sizeexprs...)))
Iext = :(tuple($(Iexprs...)))
It = tuple(Itypes...)
LD = subarray_linearindexing_dim(A, I)
strideexpr = stride1expr(A, I, :A, :I, LD)
exfirst = first_index_expr(:A, :I, length(It))
LD = subarray_linearindexing_dim(A, J)
strideexpr = stride1expr(A, J, :A, :J, LD)
exfirst = first_index_expr(:A, :J, length(It))
quote
$exfirst
SubArray{$T,$N,$A,$It,$LD}(A, $Iext, $dims, f, $strideexpr)
@@ -103,68 +121,87 @@ end

# Constructing from another SubArray
# This "pops" the old SubArray and creates a more compact one
stagedfunction slice_unsafe{T,NV,PV,IV,PLD,IndTypes}(V::SubArray{T,NV,PV,IV,PLD}, I::IndTypes)
stagedfunction slice_unsafe{T,NV,PV,IV,PLD,IndTypes}(V::SubArray{T,NV,PV,IV,PLD}, J::IndTypes)
N = 0
sizeexprs = Array(Any, 0)
indexexprs = Array(Any, 0)
Itypes = Array(Any, 0)
k = 0 # index into I
LD, die_next_vector, Ilast, isLDdone = 0, false, Void, false # for linear indexing inference
for j = 1:length(IV)
if IV[j] <: Real
push!(indexexprs, :(V.indexes[$j]))
push!(Itypes, IV[j])
# The next two Ints, if nonzero, record information about the place
# in the index tuple at which trailing dimensions got packed into a
# single Vector{Int}. For stride1 computation, we need to keep track
# of whether the index that triggered this had uniform stride.
# Iindex_lin is the spot in the resulting index tuple
# Jindex_lin is the corresponding spot in the input index tuple
Iindex_lin = Jindex_lin = 0
# Linear indexing inference makes use of the following variables:
# LD: the last dimension up through which linear indexing is efficient
# isLDdone: true if we've quit incrementing LD
# die_next_vector: if true, stop incrementing LD on the next
# "extended" input index
# jprev: holds the previous input index type
LD, die_next_vector, jprev, isLDdone = 0, false, Void, false # for linear indexing inference
Jindex = 0
for IVindex = 1:length(IV)
iv = IV[IVindex]
if iv <: Real
push!(indexexprs, :(V.indexes[$IVindex]))
push!(Itypes, iv)
if !isLDdone
LD += 1
end
else
k += 1
Jindex += 1
j = J[Jindex]
if Jindex < length(J) || Jindex == NV || IVindex == length(IV)
if !(j <: Real)
N += 1
push!(sizeexprs, dimsizeexpr(j, Jindex, length(J), :V, :J))
end
push!(indexexprs, :(reindex(V.indexes[$IVindex], J[$Jindex])))
push!(Itypes, rangetype(iv, j))
else
# We have a linear index that spans more than one
# dimension of the parent
N += 1
push!(sizeexprs, dimsizeexpr(j, Jindex, length(J), :V, :J))
push!(indexexprs, :(merge_indexes(V, V.indexes[$IVindex:end], size(V.parent)[$IVindex:end], J[$Jindex], $Jindex)))
push!(Itypes, Array{Int, 1})
Iindex_lin = length(Itypes)
Jindex_lin = Jindex
break
end
if !isLDdone
if LD < PLD
LD += 1
Ilast, LD, die_next_vector, isdone = nextLD(Ilast, I[k], LD, die_next_vector)
jprev, LD, die_next_vector, isdone = nextLD(jprev, j, LD, die_next_vector)
isLDdone |= isdone
else
if I[k] <: Real
if j <: Real
LD += 1
else
isLDdone = true
end
end
end
if k < length(I) || k == NV || j == length(IV)
if !(I[k] <: Real)
N += 1
push!(sizeexprs, dimsizeexpr(I[k], k, length(I), :V, :I))
end
push!(indexexprs, :(reindex(V.indexes[$j], I[$k])))
push!(Itypes, rangetype(IV[j], I[k]))
else
# We have a linear index that spans more than one dimension of the parent
N += 1
push!(sizeexprs, dimsizeexpr(I[k], k, length(I), :V, :I))
push!(indexexprs, :(merge_indexes(V, V.indexes[$j:end], size(V.parent)[$j:end], I[$k])))
push!(Itypes, Array{Int, 1})
break
end
end
end
for i = k+1:length(I)
if !(I[i] <: Real)
for Jind = Jindex+1:length(J)
j = J[Jind]
if !(j <: Real)
N += 1
push!(sizeexprs, dimsizeexpr(I[i], i, length(I), :V, :I))
push!(sizeexprs, dimsizeexpr(j, Jind, length(J), :V, :J))
isLDdone = true
elseif !isLDdone
LD += 1
end
push!(indexexprs, :(I[$i]))
push!(Itypes, I[i])
push!(indexexprs, :(J[$Jind]))
push!(Itypes, j)
end
Inew = :(tuple($(indexexprs...)))
dims = :(tuple($(sizeexprs...)))
It = tuple(Itypes...)
LD = max(LD, subarray_linearindexing_dim(PV, It))
strideexpr = stride1expr(PV, It, :(V.parent), :Inew, LD)
strideexpr = stride1expr(PV, It, :(V.parent), :Inew, LD, :J, Iindex_lin, Jindex_lin)
exfirst = first_index_expr(:(V.parent), :Inew, length(It))
quote
Inew = $Inew
@@ -173,52 +210,59 @@ stagedfunction slice_unsafe{T,NV,PV,IV,PLD,IndTypes}(V::SubArray{T,NV,PV,IV,PLD}
end
end

stagedfunction sub_unsafe{T,NV,PV,IV,PLD,IndTypes}(V::SubArray{T,NV,PV,IV,PLD}, I::IndTypes)
N = length(I)
while N > 0 && I[N] <: Real
stagedfunction sub_unsafe{T,NV,PV,IV,PLD,IndTypes}(V::SubArray{T,NV,PV,IV,PLD}, J::IndTypes)
N = length(J)
while N > 0 && J[N] <: Real
N -= 1
end
sizeexprs = Array(Any, 0)
indexexprs = Array(Any, 0)
Itypes = Array(Any, 0)
ItypesLD = Array(Any, 0)
preexprs = Array(Any, 0)
k = 0
LD, die_next_vector, Ilast, isLDdone = 0, false, Void, false # for linear indexing inference
for j = 1:length(IV)
if IV[j] <: Real
push!(indexexprs, :(V.indexes[$j]))
push!(Itypes, IV[j])
LD, die_next_vector, jprev, isLDdone = 0, false, Void, false
Jindex = 0
for IVindex = 1:length(IV)
iv = IV[IVindex]
if iv <: Real
push!(indexexprs, :(V.indexes[$IVindex]))
push!(Itypes, iv)
push!(ItypesLD, iv)
if !isLDdone
LD += 1
end
else
k += 1
if k <= N
push!(sizeexprs, dimsizeexpr(I[k], k, length(I), :V, :I))
Jindex += 1
j = J[Jindex]
if Jindex <= N
push!(sizeexprs, dimsizeexpr(j, Jindex, length(J), :V, :J))
end
if k < N && I[k] <: Real
if Jindex < N && j <: Real
# convert scalar to a range
sym = gensym()
push!(preexprs, :($sym = reindex(V.indexes[$j], Int(I[$k]))))
push!(preexprs, :($sym = reindex(V.indexes[$IVindex], Int(J[$Jindex]))))
push!(indexexprs, :($sym:$sym))
push!(Itypes, UnitRange{Int})
elseif k < length(I) || k == NV || j == length(IV)
push!(ItypesLD, j)
elseif Jindex < length(J) || Jindex == NV || IVindex == length(IV)
# simple indexing
push!(indexexprs, :(reindex(V.indexes[$j], I[$k])))
push!(Itypes, rangetype(IV[j], I[k]))
push!(indexexprs, :(reindex(V.indexes[$IVindex], J[$Jindex])))
push!(Itypes, rangetype(iv, j))
push!(ItypesLD, Itypes[end])
else
# We have a linear index that spans more than one dimension of the parent
push!(indexexprs, :(merge_indexes(V, V.indexes[$j:end], size(V.parent)[$j:end], I[$k])))
push!(indexexprs, :(merge_indexes(V, V.indexes[$IVindex:end], size(V.parent)[$IVindex:end], J[$Jindex], $Jindex)))
push!(Itypes, Array{Int, 1})
push!(ItypesLD, Itypes[end])
break
end
if !isLDdone
if LD < PLD
LD += 1
Ilast, LD, die_next_vector, isdone = nextLD(Ilast, I[k], LD, die_next_vector)
jprev, LD, die_next_vector, isdone = nextLD(jprev, j, LD, die_next_vector)
isLDdone |= isdone
else
if I[k] <: Real
if j <: Real
LD += 1
else
isLDdone = true
@@ -227,18 +271,21 @@ stagedfunction sub_unsafe{T,NV,PV,IV,PLD,IndTypes}(V::SubArray{T,NV,PV,IV,PLD},
end
end
end
for i = k+1:length(I)
if i <= N
push!(sizeexprs, dimsizeexpr(I[i], i, length(I), :V, :I))
for Jind = Jindex+1:length(J)
j = J[Jind]
if Jind <= N
push!(sizeexprs, dimsizeexpr(j, Jind, length(J), :V, :J))
end
push!(indexexprs, :(I[$i]))
push!(Itypes, I[i])
push!(indexexprs, :(J[$Jind]))
push!(Itypes, j)
push!(ItypesLD, Itypes[end])
end
Inew = :(tuple($(indexexprs...)))
dims = :(tuple($(sizeexprs...)))
It = tuple(Itypes...)
ItLD = tuple(ItypesLD...)
LD = max(LD, subarray_linearindexing_dim(PV, It))
strideexpr = stride1expr(PV, It, :(V.parent), :Inew, LD)
strideexpr = stride1expr(PV, ItLD, :(V.parent), :Inew, LD)
preex = isempty(preexprs) ? nothing : Expr(:block, preexprs...)
exfirst = first_index_expr(:(V.parent), :Inew, length(It))
quote
@@ -294,6 +341,7 @@ getindex(::Colon, i) = i

step(::Colon) = 1
first(::Colon) = 1
in(::Int, ::Colon) = true

## Strides
stagedfunction strides{T,N,P,I}(V::SubArray{T,N,P,I})
@@ -315,7 +363,7 @@ end

stride(V::SubArray, d::Integer) = d <= ndims(V) ? strides(V)[d] : strides(V)[end] * size(V)[end]

function stride1expr(Atype::Type, Itypes::Tuple, Aexpr, Inewsym, LD)
function stride1expr(Atype::Type, Itypes::Tuple, Aexpr, Isym, LD, Jsym=Isym, Iindex_lin=0, Jindex_lin=0)
if LD == 0
return 0
end
@@ -325,13 +373,20 @@ function stride1expr(Atype::Type, Itypes::Tuple, Aexpr, Inewsym, LD)
if I <: Real
ex = :($ex * size($Aexpr, $d))
else
ex = :($ex * step($Inewsym[$d]))
if d == Iindex_lin
ex = :($ex * step_sa($Jsym[$Jindex_lin]))
else
ex = :($ex * step($Isym[$d]))
end
break
end
end
ex
end

step_sa(arg) = step(arg)
step_sa(::Integer) = 1

# This might be conservative, but better safe than sorry
function iscontiguous{T,N,P,I,LD}(::Type{SubArray{T,N,P,I,LD}})
LD == length(I) || return false
@@ -373,41 +428,41 @@ function first_index_expr(Asym, Isym::Symbol, n::Int)
end

# Detecting whether one can support fast linear indexing
function nextLD(Ilast, I, LD, die_next_vector)
function nextLD(jprev, j, LD, die_next_vector)
isdone = false
if I <: Real
if Ilast != Void && !(Ilast <: Real)
if j <: Real
if jprev != Void && !(jprev <: Real)
die_next_vector = true
end
elseif die_next_vector || I <: Vector
elseif die_next_vector || j <: Vector
LD -= 1
isdone = true
elseif I == Colon
elseif I <: UnitRange
elseif j == Colon
elseif j <: UnitRange
die_next_vector = true
elseif I <: Range
if !(Ilast == Void || Ilast <: Real)
elseif j <: Range
if !(jprev == Void || jprev <: Real)
LD -= 1
isdone = true
end
die_next_vector = true
else
error("This shouldn't happen (linear indexing inference)")
end
Ilast = I
return Ilast, LD, die_next_vector, isdone
jprev = j
return jprev, LD, die_next_vector, isdone
end

function subarray_linearindexing_dim{A<:AbstractArray}(::Type{A}, It::Tuple)
isa(Base.linearindexing(A), Base.LinearSlow) && return 0
isempty(It) && return 0
Ilast = Void
jprev = Void
LD = 0
die_next_vector = false
while LD < length(It)
LD += 1
I = It[LD]
Ilast, LD, die_next_vector, isdone = nextLD(Ilast, I, LD, die_next_vector)
jprev, LD, die_next_vector, isdone = nextLD(jprev, I, LD, die_next_vector)
if isdone
break
end
61 changes: 41 additions & 20 deletions test/subarray.jl
Original file line number Diff line number Diff line change
@@ -167,8 +167,8 @@ function test_linear(A, B)
end

# "mixed" means 2 indexes even for N-dimensional arrays
test_mixed{T}(A::AbstractArray{T,1}, B::Array) = nothing
test_mixed{T}(A::AbstractArray{T,2}, B::Array) = nothing
test_mixed{T}(::AbstractArray{T,1}, ::Array) = nothing
test_mixed{T}(::AbstractArray{T,2}, ::Array) = nothing
test_mixed(A, B::Array) = _test_mixed(A, reshape(B, size(A)))
function _test_mixed(A, B)
L = length(A)
@@ -194,10 +194,11 @@ function err_li(I::Tuple, ld::Int, ldc::Int)
error("Linear indexing inference mismatch")
end

function err_li(S::SubArray, ldc::Int)
function err_li(S::SubArray, ld::Int, szC)
println(summary(S))
@show S.indexes
@show ldc
@show ld
@show szC
error("Linear indexing inference mismatch")
end

@@ -228,18 +229,35 @@ function runtests(A::SubArray, I...)
AA = copy_to_array(A)
# Direct test of linear indexing inference
C = Agen_nodrop(AA, I...)
ld = single_stride_dim(C)
Cld = ld = single_stride_dim(C)
Cdim = AIindex = 0
while Cdim <= Cld && AIindex < length(A.indexes)
AIindex += 1
if isa(A.indexes[AIindex], Real)
ld += 1
else
Cdim += 1
end
end
# sub
S = sub(A, I...)
local S
try
S = sub(A, I...)
catch err
@show typeof(A)
@show A.indexes
@show I
rethrow(err)
end
ldc = getLD(S)
ldc <= ld || err_li(S, ld)
ldc <= ld || err_li(S, ld, size(C))
test_linear(S, C)
test_cartesian(S, C)
test_mixed(S, C)
# slice
S = slice(A, I...)
ldc = getLD(S)
ldc <= ld || err_li(S, ld)
ldc <= ld || err_li(S, ld, size(C))
test_linear(S, C)
test_cartesian(S, C)
test_mixed(S, C)
@@ -248,28 +266,28 @@ end
# indexN is a cartesian index, indexNN is a linear index for 2 dimensions, and indexNNN is a linear index for 3 dimensions
function runviews{T}(SB::AbstractArray{T,3}, indexN, indexNN, indexNNN)
for i3 in indexN, i2 in indexN, i1 in indexN
runtests(A, i1, i2, i3)
runtests(SB, i1, i2, i3)
end
for i2 in indexNN, i1 in indexN
runtests(A, i1, i2)
runtests(SB, i1, i2)
end
for i1 in indexNNN
runtests(A, i1)
runtests(SB, i1)
end
end

function runviews{T}(SB::AbstractArray{T,2}, indexN, indexNN, indexNNN)
for i2 in indexN, i1 in indexN
runtests(A, i1, i2)
runtests(SB, i1, i2)
end
for i1 in indexNN
runtests(A, i1)
runtests(SB, i1)
end
end

function runviews{T}(SB::AbstractArray{T,1}, indexN, indexNN, indexNNN)
for i1 in indexN
runtests(A, i1)
runtests(SB, i1)
end
end

@@ -279,21 +297,24 @@ runviews{T}(SB::AbstractArray{T,0}, indexN, indexNN, indexNNN) = nothing

### Views from Arrays ###

A = reshape(1:5*7*11, 11, 7, 5)
index5 = (2, :, 2:5, 1:2:5, [4,1,5]) # all work with at least size 5
index25 = (8, :, 2:11, 12:3:22, [4,1,5,9])
index125 = (113, :, 85:121, 2:15:92, [99,14,103])
runviews(A, index5, index25, index125)

let A = reshape(1:5*7*11, 11, 7, 5)
runviews(A, index5, index25, index125)
end

### Views from views ###

B = reshape(1:13^3, 13, 13, 13)
# "outer" indexes create snips that have at least size 5 along each dimension, with the exception of Int-slicing
oindex = (:, 6, 3:7, 13:-2:1, [8,4,6,12,5,7])

for o3 in oindex, o2 in oindex, o1 in oindex
sliceB = slice(B, o1, o2, o3)
runviews(sliceB, index5, index25, index125)
let B = reshape(1:13^3, 13, 13, 13)
for o3 in oindex, o2 in oindex, o1 in oindex
sliceB = slice(B, o1, o2, o3)
runviews(sliceB, index5, index25, index125)
end
end