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.

Sign up for GitHub

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

Improve sparse vcat #10206

Merged
merged 2 commits into from
Feb 16, 2015
Merged

Improve sparse vcat #10206

Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
3e1d2f2
Fix the type instability causing slowdown and extra memory
ViralBShah Feb 15, 2015
11b9872
Separate out the type unstable part into stuffcol!.
ViralBShah Feb 16, 2015
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
49 changes: 31 additions & 18 deletions base/sparse/sparsematrix.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1884,48 +1884,61 @@ function setindex!{Tv,Ti,T<:Real}(A::SparseMatrixCSC{Tv,Ti}, x, I::AbstractVecto
end



# Sparse concatenation

function vcat(X::SparseMatrixCSC...)
Tv = promote_type(map(x->eltype(x.nzval), X)...)
Ti = promote_type(map(x->eltype(x.rowval), X)...)

vcat(map(x->convert(SparseMatrixCSC{Tv,Ti}, x), X)...)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Rather than effectively converting the input arrays twice, another approach might be to break out the "stuff a single matrix into its slot" algorithm into a subfunction. That way the performance-intensive part will be type-stable.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I didn't fully understand. Input arrays that are of type {Tv, Ti} will anyways be no-op for convert. Worst case, most of them get converted once and then passed to vcat, which stuffs it in. How do we end up converting twice?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I shouldn't have said "convert twice," I should have said "copy twice."

To be clearer, my suggestion is:

function vcat(X::SparseMatrixCSC...)
    Tv = promote_type(map(x->eltype(x.nzval), X)...)
    Ti = promote_type(map(x->eltype(x.rowval), X)...)
    <compute sizes and allocate outputs here>
    moffset = 0
    for x in X
        stuff!(colptr, rowval, nzval, moffset, x)
        moffset += size(x,1)
    end
    < pack as sparse matrix >
end

That way, even though x might be type-unstable, you don't really care, because it will dispatch to a type-specialized implementation of stuff!, and it's that function that has all the compute-intensive operations.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am not really familiar with the sparse code, but I agree that @timholy is probably right. It should not be necessary to convert a complete array. By allocating the output array with the promoted type and then assigning into it, the convert should be taken care during the assignment at the level of single elements.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Converting is a pretty uncommon case, which is why I thought that it was worth the simplicity to just convert. Anyways, I tried implementing stuffcol! (code below in comments), but that is still slower than my current implementation, and takes more memory.

end

function vcat{Tv,Ti<:Integer}(X::SparseMatrixCSC{Tv,Ti}...)
num = length(X)
mX = [ size(x, 1) for x in X ]
nX = [ size(x, 2) for x in X ]
m = sum(mX)
n = nX[1]

for i = 2 : num
if nX[i] != n; throw(DimensionMismatch("")); end
if nX[i] != n
throw(DimensionMismatch("All inputs to vcat should have the same number of columns"))
end
end
m = sum(mX)

Tv = promote_type(map(x->eltype(x.nzval), X)...)
Ti = promote_type(map(x->eltype(x.rowval), X)...)

colptr = Array(Ti, n + 1)
nnzX = [ nnz(x) for x in X ]
nnz_res = sum(nnzX)
colptr = Array(Ti, n + 1)
rowval = Array(Ti, nnz_res)
nzval = Array(Tv, nnz_res)
nzval = Array(Tv, nnz_res)

colptr[1] = 1
@inbounds for c = 1 : n
for c = 1:n
mX_sofar = 0
rr1 = colptr[c]
ptr_res = colptr[c]
for i = 1 : num
XI = X[i]
rX1 = XI.colptr[c]
rX2 = XI.colptr[c + 1] - 1
rr2 = rr1 + (rX2 - rX1)
Xi = X[i]
colptrXi = Xi.colptr
rowvalXi = Xi.rowval
nzvalXi = Xi.nzval

col_length = (colptrXi[c + 1] - 1) - colptrXi[c]
ptrXi = colptrXi[c]
for k=ptr_res:(ptr_res + col_length)
@inbounds rowval[k] = rowvalXi[ptrXi] + mX_sofar
@inbounds nzval[k] = nzvalXi[ptrXi]
ptrXi += 1
end

rowval[rr1 : rr2] = XI.rowval[rX1 : rX2] .+ mX_sofar
nzval[rr1 : rr2] = XI.nzval[rX1 : rX2]
ptr_res += col_length + 1
mX_sofar += mX[i]
rr1 = rr2 + 1
end
colptr[c + 1] = rr1
colptr[c + 1] = ptr_res
end
SparseMatrixCSC(m, n, colptr, rowval, nzval)
end


function hcat(X::SparseMatrixCSC...)
num = length(X)
mX = [ size(x, 1) for x in X ]
Expand Down
2 changes: 2 additions & 0 deletions test/sparse/sparse.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,8 @@ do33 = ones(3)

# check vert concatenation
@test all([se33; se33] == sparse([1, 4, 2, 5, 3, 6], [1, 1, 2, 2, 3, 3], ones(6)))
se33_32bit = convert(SparseMatrixCSC{Float32,Int32}, se33)
@test all([se33; se33_32bit] == sparse([1, 4, 2, 5, 3, 6], [1, 1, 2, 2, 3, 3], ones(6)))

# check h+v concatenation
se44 = speye(4)
Expand Down