DifferentialEquations.jl supports GPUs. The question is what is the right way to utilize DiffEq's GPU support here. That's actually a much deeper question than

DifferentialEquations.jl supports GPUs. The question is what is the right way to utilize DiffEq's GPU support here. That's actually a much deeper question than you might at first think. The reason is because sparse factorizations are not an inherently parallel operation, so CUSOLVER is actually quite bad at doing them. It even seems that CUSOLVER does some CPU operations in there... ouch. But what can we do with ModelingToolkit to start pushing more into this domain? Well, here's a few things: 1) We plan to have new compilation targets, kind of like the SArray one that was shown, which target parallelism. For example, for large Jacobians, we want a target so you can say "generate the code that fills the Jacobian optimally with 16 threads", and then what it will do is split up the Jacobian into 16 parts with approximately equally many operations and thread between those. Since it's all scalar operations then this should be something that gets as efficient as you'd want to write by hand! Then a similar compilation target can be made for the GPU, where it would need to know the number of GPU threads you'd be using and it would split apart that Jacobian filling and write a gigantic `if threadid() == ...` setup, and use CUDAnative.jl to compile this to a .ptx kernel so that all of those threads are parallelized. Both of these trivially extend to sparse matrices doing what I showed in the writing. So filling sparse matrices is something we can and plan to have nice compilation targets for. 2) We can also have compilation targets for `J*v` operations, so that way it can accelerate Newton-Krylov solvers, which is probably the better way of using sparse matrices on the GPU since it would avoid sparse factorizations. 3) We can also implement a generic sparse factorization in Julia, to then build analytical solutions for sparse factorizations, run Rewrite.jl simplification passes on these, and generate the parallelized kernels for them. This would then be something that does sparse factorizations on the GPU with parallelized code but does not need to use CUSOLVER, and so that would be quite ideal. The issue is that it would require writing generic sparse factorization tools, since Julia currently relies on SuiteSparse for this with no generic fallback. This is not an easy task, but I think Yingbo will be up for it in a few years since we have a plan to not use standard sparse linear algebra routines in Julia's differential equation suites because there is a way to speed up operations specifically in this domain by using a new factorization, and so if we go this direction our first task will be to learn how the current sparse LU is implemented.