alternate getindex_general implementation

getindex_general.jl

function getindex_general{Tv,Ti}(A::SparseMatrixCSC{Tv,Ti}, I::AbstractVector, J::AbstractVector)
    (m, n) = size(A)
    nI = length(I)
    nJ = length(J)

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

    spI = sortperm(I)
    @inbounds sI = I[spI]
    spI = sortperm(spI)
    W = Array(Tv, nI)
    nzv = zv = zero(Tv)

    nR = 0
    for j = 1:nJ
        @inbounds col = J[j]
        ptrA::Int = colptrA[col]
        stopA::Int = colptrA[col+1]-1
        (ptrA > stopA) && continue

        for i in 1:nI
            @inbounds row = sI[i]
            @inbounds ptrA = (rowvalA[ptrA] < row) ? searchsortedfirst(rowvalA, row, ptrA, stopA, Forward) : ptrA
            (ptrA > stopA) && break
            @inbounds (rowvalA[ptrA] == row) && (nR += 1)
        end
    end

    nzvalR = Array(Tv, nR)
    rowvalR = Array(Ti, nR)
    colptrR = Array(Ti, nJ+1)
    ptrR::Int = 1
    rowR::Int = 1
    colptrR[1] = 1

    @inbounds for j = 1:nJ
        col = J[j]
        ptrA::Int = colptrA[col]
        stopA::Int = colptrA[col+1]-1

        if ptrA <= stopA
            rowR = 1
            for i in 1:nI
                row = sI[i]
                ptrA = (rowvalA[ptrA] < row) ? searchsortedfirst(rowvalA, row, ptrA, stopA, Forward) : ptrA
                if ptrA <= stopA
                    W[i] = (rowvalA[ptrA] == row) ? nzvalA[ptrA] : zv
                else
                    W[i:nI] = zv
                    break
                end
            end
            for i in 1:nI
                nzv = W[spI[i]]
                if nzv != zv
                    rowvalR[ptrR] = rowR
                    nzvalR[ptrR] = nzv
                    ptrR += 1
                end
                rowR += 1
            end
        end

        colptrR[j+1] = ptrR
    end

    SparseMatrixCSC{Tv,Ti}(nI, nJ, colptrR, rowvalR, nzvalR)
end