Grid concatenation of Sparse Matrices in julia.

gridcat.jl

function cat2d(X::Array{SparseMatrixCSC,2})
    iX = ones(Int, size(X)...)

    mX, nX = size(X)

    m = sum([size(X[i,1], 1) for i in 1:size(X, 1)])
    n = sum([size(X[1,i], 2) for i in 1:size(X, 2)])

    Tv = eltype(X[1].nzval)
    Ti = eltype(X[1].rowval)
    for i = 2:length(X)
        Tv = promote_type(Tv, eltype(X[i].nzval))
        Ti = promote_type(Ti, eltype(X[i].rowval))
    end

    # todo: Check for dimensional equality in rows and columns.
    nzl = sum(map(nnz, X[1:end]))
    nzval = Array(Tv, nzl)
    rowval = Array(Ti, nzl)
    colptr = Array(Ti, n+1)

    ptr = 1
    cp_i = 1

    colptr[1] = 1

    for j in 1:mX
        for col in 1:size(X[1,j], 2)
            column_count = zero(Tv)
            for i in 1:nX
                x = X[i,j]
                ix = iX[i,j]
                xcolptr = x.colptr
                xrowval = x.rowval
                xnzval = x.nzval

                col_length = xcolptr[col + 1] - xcolptr[col]

                for k = ptr : (ptr + col_length - 1)
                    @inbounds rowval[k] =  xrowval[ix] + column_count
                    @inbounds nzval[k] = xnzval[ix]
                    ix += 1
                end

                ptr += col_length
                column_count += size(x, 1)
                iX[i,j] = ix
            end
            cp_i += 1
            colptr[cp_i] = ptr
        end
    end
    SparseMatrixCSC{Tv, Ti}(m, n, colptr, rowval, nzval)
end

ma = Array{SparseMatrixCSC, 2}(2, 2)
ma[1:4] = [round(Int, sprand(1000, 1000, 0.1)*1) for i in 1:4]

println("2D concatenation")
for i in 1:10
    @time cat2d(ma);
    gc();
end
println("Hvcat concatenation")
for i in 1:10
    @time hvcat((2,2), ma[1,1], ma[1,2], ma[2,1], ma[2,2]);
    gc();
end

@time spw = cat2d(ma);
spr = hvcat((2,2), ma[1,1], ma[1,2], ma[2,1], ma[2,2]);

println(spw == spr)