Merge branch 'master' into offset1offsetarrays

Author: chriselrod

src/ArrayInterface.jl

@@ -25,6 +25,8 @@ else
     end
 end
 
+const CanonicalInt = Union{Int,StaticInt}
+
 if VERSION ≥ v"1.6.0-DEV.1581"
     _is_reshaped(::Type{ReinterpretArray{T,N,S,A,true}}) where {T,N,S,A} = true
     _is_reshaped(::Type{ReinterpretArray{T,N,S,A,false}}) where {T,N,S,A} = false
@@ -48,6 +50,7 @@ const LoTri{T,M} = Union{LowerTriangular{T,M},UnitLowerTriangular{T,M}}
 @inline static_step(x) = Static.maybe_static(known_step, step, x)
 
 include("ndindex.jl")
+include("array_index.jl")
 
 """
     parent_type(::Type{T})
@@ -264,240 +267,6 @@ function findstructralnz(x::Diagonal)
     (1:n, 1:n)
 end
 
-abstract type MatrixIndex end
-
-struct BidiagonalIndex <: MatrixIndex
-    count::Int
-    isup::Bool
-end
-
-struct TridiagonalIndex <: MatrixIndex
-    count::Int# count==nsize+nsize-1+nsize-1
-    nsize::Int
-    isrow::Bool
-end
-
-struct BandedMatrixIndex <: MatrixIndex
-    count::Int
-    rowsize::Int
-    colsize::Int
-    bandinds::Array{Int,1}
-    bandsizes::Array{Int,1}
-    isrow::Bool
-end
-
-function _bandsize(bandind, rowsize, colsize)
-    -(rowsize - 1) <= bandind <= colsize - 1 || throw(ErrorException("Invalid Bandind"))
-    if (bandind * (colsize - rowsize) > 0) & (abs(bandind) <= abs(colsize - rowsize))
-        return min(rowsize, colsize)
-    elseif bandind * (colsize - rowsize) <= 0
-        return min(rowsize, colsize) - abs(bandind)
-    else
-        return min(rowsize, colsize) - abs(bandind) + abs(colsize - rowsize)
-    end
-end
-
-function BandedMatrixIndex(rowsize, colsize, lowerbandwidth, upperbandwidth, isrow)
-    upperbandwidth > -lowerbandwidth || throw(ErrorException("Invalid Bandwidths"))
-    bandinds = upperbandwidth:-1:-lowerbandwidth
-    bandsizes = [_bandsize(band, rowsize, colsize) for band in bandinds]
-    BandedMatrixIndex(sum(bandsizes), rowsize, colsize, bandinds, bandsizes, isrow)
-end
-
-struct BlockBandedMatrixIndex <: MatrixIndex
-    count::Int
-    refinds::Array{Int,1}
-    refcoords::Array{Int,1}# storing col or row inds at ref points
-    isrow::Bool
-end
-
-function BlockBandedMatrixIndex(nrowblock, ncolblock, rowsizes, colsizes, l, u)
-    blockrowind = BandedMatrixIndex(nrowblock, ncolblock, l, u, true)
-    blockcolind = BandedMatrixIndex(nrowblock, ncolblock, l, u, false)
-    sortedinds = sort(
-        [(blockrowind[i], blockcolind[i]) for i = 1:length(blockrowind)],
-        by = x -> x[1],
-    )
-    sort!(sortedinds, by = x -> x[2], alg = InsertionSort)# stable sort keeps the second index in order
-    refinds = Array{Int,1}()
-    refrowcoords = Array{Int,1}()
-    refcolcoords = Array{Int,1}()
-    rowheights = cumsum(pushfirst!(copy(rowsizes), 1))
-    blockheight = 0
-    blockrow = 1
-    blockcol = 1
-    currenti = 1
-    lastrowind = sortedinds[1][1] - 1
-    lastcolind = sortedinds[1][2]
-    for ind in sortedinds
-        rowind, colind = ind
-        if colind == lastcolind
-            if rowind > lastrowind
-                blockheight += rowsizes[rowind]
-            end
-        else
-            for j = blockcol:blockcol+colsizes[lastcolind]-1
-                push!(refinds, currenti)
-                push!(refrowcoords, blockrow)
-                push!(refcolcoords, j)
-                currenti += blockheight
-            end
-            blockcol += colsizes[lastcolind]
-            blockrow = rowheights[rowind]
-            blockheight = rowsizes[rowind]
-        end
-        lastcolind = colind
-        lastrowind = rowind
-    end
-    for j = blockcol:blockcol+colsizes[lastcolind]-1
-        push!(refinds, currenti)
-        push!(refrowcoords, blockrow)
-        push!(refcolcoords, j)
-        currenti += blockheight
-    end
-    push!(refinds, currenti)# guard
-    push!(refrowcoords, -1)
-    push!(refcolcoords, -1)
-    rowindobj = BlockBandedMatrixIndex(currenti - 1, refinds, refrowcoords, true)
-    colindobj = BlockBandedMatrixIndex(currenti - 1, refinds, refcolcoords, false)
-    rowindobj, colindobj
-end
-
-struct BandedBlockBandedMatrixIndex <: MatrixIndex
-    count::Int
-    refinds::Array{Int,1}
-    refcoords::Array{Int,1}# storing col or row inds at ref points
-    reflocalinds::Array{BandedMatrixIndex,1}
-    isrow::Bool
-end
-
-function BandedBlockBandedMatrixIndex(
-    nrowblock,
-    ncolblock,
-    rowsizes,
-    colsizes,
-    l,
-    u,
-    lambda,
-    mu,
-)
-    blockrowind = BandedMatrixIndex(nrowblock, ncolblock, l, u, true)
-    blockcolind = BandedMatrixIndex(nrowblock, ncolblock, l, u, false)
-    sortedinds = sort(
-        [(blockrowind[i], blockcolind[i]) for i = 1:length(blockrowind)],
-        by = x -> x[1],
-    )
-    sort!(sortedinds, by = x -> x[2], alg = InsertionSort)# stable sort keeps the second index in order
-    rowheights = cumsum(pushfirst!(copy(rowsizes), 1))
-    colwidths = cumsum(pushfirst!(copy(colsizes), 1))
-    currenti = 1
-    refinds = Array{Int,1}()
-    refrowcoords = Array{Int,1}()
-    refcolcoords = Array{Int,1}()
-    reflocalrowinds = Array{BandedMatrixIndex,1}()
-    reflocalcolinds = Array{BandedMatrixIndex,1}()
-    for ind in sortedinds
-        rowind, colind = ind
-        localrowind =
-            BandedMatrixIndex(rowsizes[rowind], colsizes[colind], lambda, mu, true)
-        localcolind =
-            BandedMatrixIndex(rowsizes[rowind], colsizes[colind], lambda, mu, false)
-        push!(refinds, currenti)
-        push!(refrowcoords, rowheights[rowind])
-        push!(refcolcoords, colwidths[colind])
-        push!(reflocalrowinds, localrowind)
-        push!(reflocalcolinds, localcolind)
-        currenti += localrowind.count
-    end
-    push!(refinds, currenti)
-    push!(refrowcoords, -1)
-    push!(refcolcoords, -1)
-    rowindobj = BandedBlockBandedMatrixIndex(
-        currenti - 1,
-        refinds,
-        refrowcoords,
-        reflocalrowinds,
-        true,
-    )
-    colindobj = BandedBlockBandedMatrixIndex(
-        currenti - 1,
-        refinds,
-        refcolcoords,
-        reflocalcolinds,
-        false,
-    )
-    rowindobj, colindobj
-end
-
-Base.firstindex(ind::MatrixIndex) = 1
-Base.lastindex(ind::MatrixIndex) = ind.count
-Base.length(ind::MatrixIndex) = ind.count
-function Base.getindex(ind::BidiagonalIndex, i::Int)
-    1 <= i <= ind.count || throw(BoundsError(ind, i))
-    if ind.isup
-        ii = i + 1
-    else
-        ii = i + 1 + 1
-    end
-    convert(Int, floor(ii / 2))
-end
-
-function Base.getindex(ind::TridiagonalIndex, i::Int)
-    1 <= i <= ind.count || throw(BoundsError(ind, i))
-    offsetu = ind.isrow ? 0 : 1
-    offsetl = ind.isrow ? 1 : 0
-    if 1 <= i <= ind.nsize
-        return i
-    elseif ind.nsize < i <= ind.nsize + ind.nsize - 1
-        return i - ind.nsize + offsetu
-    else
-        return i - (ind.nsize + ind.nsize - 1) + offsetl
-    end
-end
-
-function Base.getindex(ind::BandedMatrixIndex, i::Int)
-    1 <= i <= ind.count || throw(BoundsError(ind, i))
-    _i = i
-    p = 1
-    while _i - ind.bandsizes[p] > 0
-        _i -= ind.bandsizes[p]
-        p += 1
-    end
-    bandind = ind.bandinds[p]
-    startfromone = !xor(ind.isrow, (bandind > 0))
-    if startfromone
-        return _i
-    else
-        return _i + abs(bandind)
-    end
-end
-
-function Base.getindex(ind::BlockBandedMatrixIndex, i::Int)
-    1 <= i <= ind.count || throw(BoundsError(ind, i))
-    p = 1
-    while i - ind.refinds[p] >= 0
-        p += 1
-    end
-    p -= 1
-    _i = i - ind.refinds[p]
-    if ind.isrow
-        return ind.refcoords[p] + _i
-    else
-        return ind.refcoords[p]
-    end
-end
-
-function Base.getindex(ind::BandedBlockBandedMatrixIndex, i::Int)
-    1 <= i <= ind.count || throw(BoundsError(ind, i))
-    p = 1
-    while i - ind.refinds[p] >= 0
-        p += 1
-    end
-    p -= 1
-    _i = i - ind.refinds[p] + 1
-    ind.reflocalinds[p][_i] + ind.refcoords[p] - 1
-end
-
 function findstructralnz(x::Bidiagonal)
     n = Base.size(x, 1)
     isup = x.uplo == 'U' ? true : false