Generating sparse matrices slow for combined operators

[bclyons12]

Calling sparse on a DiffEqOperator is much slower and more memory intensive than building the same SparseMatrixCSC from three vectors. For example

N = 10001
Rᵢ = 20.
r = collect(range(0.,Rᵢ,length=N))
Δr = r[2]

function by_hand()

    R = zeros(3N-2)
    C = zeros(3N-2)
    V = zeros(3N-2)
    for i in 2:(N-1)

        # diagonal - 1
        j = 3*(i-1)
        R[j], C[j], V[j] = i, i-1, 1.0/Δr^2

        # diagonal
        j = 3*(i-1) + 1
        R[j], C[j], V[j] = i, i, -2.0/Δr^2

        # diagonal - 1
        j = 3*(i-1) + 2
        R[j], C[j], V[j] = i, i+1, 1.0/Δr^2

    end
    # Zero derivative at boundaries (Neumann)
    j= 1
    R[j], C[j], V[j] = 1, 1, 1.
    j = 2
    R[j], C[j], V[j] = 1, 2, -1.

    j = 3N - 2
    R[j], C[j], V[j] = N, N, 1.
    j = 3N - 3
    R[j], C[j], V[j] = N, N-1, -1.

    A = sparse(R,C,V)
end

function by_DEO()
    Δ = CenteredDifference(2,2,Δr,N)
    bc = Neumann0BC(Δr,1)
    A,b = sparse(Δ*bc)
    return A
end

println("Create by hand")
@btime Abh = by_hand();
println("Create by DiffEqOperators")
@btime Abd = by_DEO();

Output:

Create by hand
  6.089 ms (185603 allocations: 5.27 MiB)
Create by DiffEqOperators
  5.218 s (257528 allocations: 6.01 GiB)

I would have thought sparse on a DerivativeOperator or GhostDerivativeOperator would have done something like my by_hand function. @ChrisRackauckas thought that it may be falling back to an iteration on the indices method in our Discourse.

[ChrisRackauckas]

It's not DiffEqOperators, it's specifically the type for the combined operator is missing a sparse matrix constructor.

[bclyons12]

@ChrisRackauckas I think the sparse matrix constructor for DerivativeOperator isn't working efficiently either. If I just call sparse(Δ) in my by_DEO function, it comes out at approximately the same timing.