spzeros and speye are really slow

[jebej]

I searched and couldn't find any open issue on this, which surprised me.

Not too sure what the problem is with the code in Base, but if I write my own spzeros, it is much faster (the same thing happens with speye):

julia> using BenchmarkTools

julia> @benchmark spzeros(5,5)
BenchmarkTools.Trial:
  memory estimate:  1.11 KiB
  allocs estimate:  18
  --------------
  minimum time:     9.330 μs (0.00% GC)
  median time:      9.953 μs (0.00% GC)
  mean time:        10.215 μs (0.00% GC)
  maximum time:     60.647 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     1

julia> function my_spzeros(m,n)
           ((m < 0) || (n < 0)) && throw(ArgumentError("invalid Array dimensions"))
           rowval = Vector{Int}(0)
           colptr = ones(Int,n+1)
           nzval  = Vector{Float64}(0)
           return SparseMatrixCSC(m,n,colptr,rowval,nzval)
       end
my_spzeros (generic function with 1 method)

julia> @benchmark my_spzeros(5,5)
BenchmarkTools.Trial:
  memory estimate:  336 bytes
  allocs estimate:  4
  --------------
  minimum time:     88.951 ns (0.00% GC)
  median time:      97.654 ns (0.00% GC)
  mean time:        117.422 ns (13.87% GC)
  maximum time:     1.565 μs (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     965

[nalimilan]

Good catch. This is because spzeros(5, 5) ends up calling spzeros(Float64, Int, 5, 5), and the two first arguments to this method are not specialized on. Changing the signature of the method from

myspzeros(Tv::Type, Ti::Type, m::Integer, n::Integer)

to:

myspzeros(::Type{Tv}, ::Type{Ti}, m::Integer, n::Integer) where {Ti, Tv}

fixes the problem.

Would you make a pull request to make this change?

[jebej]

It seems to be an issue with having the type arguments being function calls.

I have actually noticed this problem in my own code. In this function, if I just call ptrace in the type constructor performance drops by a factor of 20.

[jebej]

Hm, that does make sense. What problem am I seeing with my code then?

[jebej]

Just to be extra clear what I am doing is replacing

function ptrace(ρ::Density, out)
    res = ptrace(full(ρ),out,dims(ρ))
    return Density(res[1],res[2])
end

with

function ptrace(ρ::Density, out)
    return Density(ptrace(full(ρ),out,dims(ρ))...)
end

[KristofferC]

Splatting penalty. Not related to the sparse issue, see: #13359.

[jebej]

I see. Is that penalty also expected if I was doing:

function ptrace(ρ::Density, out)
    res = ptrace(full(ρ),out,dims(ρ))
    return Density(res...)
end

This last form is slower than writing out res[1],res[2], but by 20ns 6ns, not microseconds.
Edit: it's pretty much the same after trying a few times.

[KristofferC]

I am surprised there is a difference between the two.

[jebej]

With the example in the docstrings, I get:

Density(res[1],res[2])

julia> @benchmark ptrace($ρ,2)
BenchmarkTools.Trial:
  memory estimate:  592 bytes
  allocs estimate:  7
  --------------
  minimum time:     283.362 ns (0.00% GC)
  median time:      300.149 ns (0.00% GC)
  mean time:        358.231 ns (14.33% GC)
  maximum time:     7.924 μs (94.38% GC)
  --------------
  samples:          10000
  evals/sample:     315

Density(ptrace(full(ρ),out,dims(ρ))...)

julia> @benchmark ptrace($ρ,2)
BenchmarkTools.Trial:
  memory estimate:  608 bytes
  allocs estimate:  8
  --------------
  minimum time:     6.531 μs (0.00% GC)
  median time:      6.842 μs (0.00% GC)
  mean time:        7.013 μs (0.69% GC)
  maximum time:     501.158 μs (96.20% GC)
  --------------
  samples:          10000
  evals/sample:     5

Density(res...)

julia> @benchmark ptrace($ρ,2)
BenchmarkTools.Trial:
  memory estimate:  592 bytes
  allocs estimate:  7
  --------------
  minimum time:     283.185 ns (0.00% GC)
  median time:      306.042 ns (0.00% GC)
  mean time:        365.921 ns (14.18% GC)
  maximum time:     9.788 μs (95.12% GC)
  --------------
  samples:          10000
  evals/sample:     313

[jebej]

So if I understand correctly, the issue only occurs when splatting directly on a function call?