Yep, the Jacobian is sparse. I'm using Julia a bit for work and play, but don't own a GPU to test out their GPU support quality.

Does anyone (I assume the poster I originally replied to) know which situations are even worth it to is solve the Jacobian on a GPU?

One of the questions I have about the Julia diffeq module is how many simultaneous equations can it take? Is it thousands?

DifferentialEquations.jl can take as many simultaneous equations as you have memory for. You can even use an MPIArray when you run out. When we solve PDEs, sometimes I'm testing with a state space of tens of millions of ODEs of course, which comes from a Method of Lines discretization. Right now we are doing a lot of tests in the thousands of ODEs range since we have a reader for chemical reaction networks which we are using to optimize on test models. This is all with directly writing the Julia functions.

If you question is, how far can ModelingToolkit be pushed? That's unknown. We made sure it's quite an efficient representation, but we haven't stress tested it to see the memory consumption at a million operations yet. Since it's able to generate thousands of operations in symbolic Jacobian inversions, it easy goes to thousands and tens of thousands of equations. How much larger? Who knows. What we can say is that its efficiency is more similar to SymEngine's symbolic types, so if SymEngine handles the millions of equations range (on a computer with enough RAM, symbolic representations may use up quite a bit here!) than we should be able to.

Wow...again this is all really impressive as is you taking the time to reply.

I'm looking forward to starting to play with these new libraries. Is there a more focused place to talk about Julia applied to numerical work (linear algebra and nonlinear equation solvers).

Looks like you can transfer sparse matrices to the GPU and use cuSOLVER to solve least-squares [1] but there isn't a nice interface for it yet (https://github.com/JuliaGPU/CuArrays.jl/blob/master/test/spa...).