sparse_getindex_range_of_cols.jl

# This performance regression seems to be fixed now.  Maybe one of the many commits of @vtjnash?
#
# sandbox/performance_regressions>> julia sg.jl # 664cab522c0
# elapsed time: 0.268316893 seconds (14911976 bytes allocated)
# elapsed time: 0.267490744 seconds (14837704 bytes allocated)
# elapsed time: 0.275272316 seconds (14793688 bytes allocated, 3.10% gc time)
# sandbox/performance_regressions>> juliamf sg.jl # around 66cfa99913
# elapsed time: 0.413370651 seconds (14868712 bytes allocated)
# elapsed time: 0.41186328 seconds (14828776 bytes allocated)
# elapsed time: 0.427139988 seconds (14916648 bytes allocated, 1.98% gc time)
# sandbox/performance_regressions>> julia0.2.1 sg.jl 
# elapsed time: 0.296622719 seconds (15920432 bytes allocated)
# elapsed time: 0.303510414 seconds (15881312 bytes allocated)
# elapsed time: 0.306049353 seconds (16140896 bytes allocated)

small = 10^3
med = 10^4
large = 10^5
huge = 10^6

n = small
sp = sprand(n, n, 1e-1)

function getindex_I_sorted2{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, I::Vector, J::AbstractVector)

    (m, n) = size(A)
    nI = length(I)
    nJ = length(J)

    colptrA = A.colptr; rowvalA = A.rowval; nzvalA = A.nzval

    I_ref = falses(m)
    I_ref[I] = true

    I_repeat = zeros(Int, m)
    for i=1:nI; I_repeat[I[i]] += 1; end

    colptrS = Array(Ti, nJ+1)
    colptrS[1] = 1
    nnzS = 0

    # Form the structure of the result and compute space
    for j = 1:nJ
        col = J[j]

        for k = colptrA[col]:colptrA[col+1]-1
            rowA = rowvalA[k]
            
            if I_ref[rowA]
                for r = 1:I_repeat[rowA]
                    nnzS += 1
                end
            end

        end
        colptrS[j+1] = nnzS+1
    end

    # Populate the values in the result
    rowvalS = Array(Ti, nnzS)
    nzvalS  = Array(Tv, nnzS)
    ptrS    = 1

    fI = zeros(Ti, m)
    for k=1:nI
        Ik = I[k]
        if fI[Ik] == 0; fI[Ik] = k; end
    end

    for j = 1:nJ
        col = J[j]

        for k = colptrA[col]:colptrA[col+1]-1
            rowA = rowvalA[k]
            
            if I_ref[rowA]
                for r = 1:I_repeat[rowA]
                    rowvalS[ptrS] = fI[rowA] + r - 1
                    nzvalS[ptrS] = nzvalA[k]
                    ptrS += 1
                end
            end

        end
    end

    return SparseMatrixCSC(nI, nJ, colptrS, rowvalS, nzvalS)
end

function fcols(sp)
    acc = zero(eltype(sp))
    for i=rand(1:500, 500)
        inds = [20:500+i]
        acc += sum(getindex_I_sorted2(sp, [i], inds))
    end
end

a =  fcols(sp)
gc()
@time fcols(sp)
@time fcols(sp)
@time fcols(sp)