Bad performance on large Rodas5 problem

[jonniedie]

ComponentArrays are up to 10x slower than LabelledArrays for this problem! ComponentArrays broadcast to ComponentArrays and LabelledArrays broadcast to plain Arrays, so there is that disadvantage for these types of problems, but the disparity shouldn't be nearly this high. All of the microbenchmarks show that ComponentArrays should be almost the same as normal Arrays. From profiling, it appears that there is are plain flat indices that are making it through to the specialized @generated _getindex, which shouldn't be happening. I suspect that's the main problem

using ComponentArrays
using DifferentialEquations
using LabelledArrays

function f1(du,u,p,t)
   du.x .= -1 .* u.x .* u.y .* p[1]
   du.y .= -1 .* u.y .* p[2]
end

const n = 1000

p = [0.1,0.1]

lu_0 = @LArray  fill(1000.0,2*n) (x = (1:n),y = (n+1:2*n))
cu_0 =ComponentArray(x=fill(1000.0, n), y=fill(1000.0, n))

lprob1 = ODEProblem(f1,lu_0,(0,100.0),p)
cprob1 = ODEProblem(f1,cu_0,(0,100.0),p)

lsol = solve(lprob1, Rodas5());
lsol = solve(lprob1, Rodas5(autodiff=false));
csol = solve(cprob1, Rodas5());
csol = solve(cprob1, Rodas5(autodiff=false));

[jonniedie]

So one big issue is that Colon wasn't counted as a flat index for fall through indexing. But we're still about 6x slower. It looks like the 'LabelledArray's version is using a LU factorization because of a if typeof(A) <: Matrix branch whereas the ComponentArrays version falls back to a QR factorization.

[jonniedie]

DiffEqBase > src > linear_nonlinear.jl > line 63

    if typeof(A) <: Matrix
      blasvendor = BLAS.vendor()
      # if the user doesn't use OpenBLAS, we assume that is a better BLAS
      # implementation like MKL
      #
      # RecursiveFactorization seems to be consistantly winning below 100
      # https://discourse.julialang.org/t/ann-recursivefactorization-jl/39213
      if ArrayInterface.can_setindex(x) && (size(A,1) <= 100 || ((blasvendor === :openblas || blasvendor === :openblas64) && size(A,1) <= 500))
        p.A = RecursiveFactorization.lu!(A)
      else
        p.A = lu!(A)
      end
    elseif typeof(A) <: Tridiagonal
      p.A = lu!(A)
    elseif typeof(A) <: Union{SymTridiagonal}
      p.A = ldlt!(A)
    elseif typeof(A) <: Union{Symmetric,Hermitian}
      p.A = bunchkaufman!(A)
    elseif typeof(A) <: SparseMatrixCSC
      p.A = lu(A)
    elseif ArrayInterface.isstructured(A)
      p.A = factorize(A)
    elseif !(typeof(A) <: AbstractDiffEqOperator)
      # Most likely QR is the one that is overloaded
      # Works on things like CuArrays
      p.A = qr(A)
    end

[jonniedie]

🚀🚀🚀

julia> @time sol = solve(lprob1, Rodas5());
  2.630749 seconds (1.30 M allocations: 134.332 MiB)

julia> @time sol = solve(cprob1, Rodas5());
  2.419173 seconds (1.30 M allocations: 141.137 MiB)

julia> @time sol = solve(lprob1, Rodas5(autodiff=false));
  1.824586 seconds (403.66 k allocations: 82.395 MiB)

julia> @time sol = solve(cprob1, Rodas5(autodiff=false));
  1.786986 seconds (405.70 k allocations: 89.200 MiB)

[jonniedie]

Looks like this is back to being slow. Not sure why.

julia> @time sol = solve(lprob1, Rodas5());
  2.711504 seconds (923.22 k allocations: 127.437 MiB, 0.76% gc time)

julia> @time sol = solve(cprob1, Rodas5());
 10.947427 seconds (1.22 M allocations: 146.974 MiB, 0.19% gc time)

julia> @time sol = solve(lprob1, Rodas5(autodiff=false));
  2.313289 seconds (2.10 k allocations: 63.985 MiB, 0.88% gc time)

julia> @time sol = solve(cprob1, Rodas5(autodiff=false));
 10.346782 seconds (302.46 k allocations: 82.064 MiB)

ComponentArrays don't go brrr.

[bgroenks96]

This was fixed by adding the diffeqbc dispatch, right?