Wow, that's really cool. I feel like DiffEq is one of those libraries that seems to leverage every one of Julia's cool features that I can think of.
Making differential equations 'just work' on a GPU is super cool. I feel like if this were almost any other language, it'd be impossible to get this level of code-reuse where a user can just bring their own differential equation that may not even know about GPUs and then have the differential equation library solve it on the GPU by just giving a GPUarray as the initial condition.
Yeah, it is amazing to see how one can really rethink a core mathematical library with new tools and ideas. My mind continues to be blown every time I hear Chris talk about Measurements.jl and DifferentialEquations.jl just work, and you get sensitivity analysis for free (and neither package had planned for the other).
Do you mean they are simpler to use? Because those are just models of parallel computation. A little like saying Turing machines are simpler than python.
I'm contemplating writing some fun electromagnetic simulator, and so far had my eyes on Rust and/or Futhark with some prototypes in Python even Matlab, but now I may need to consider Julia too...
In my own recent experience, Julia is awesome but documentation could be sometimes clearer (eg. missing examples with each method, missing cross references). Also, when searching for details on some features, it's hard to figure which version someone is talking about, it seems some versions (0.6 for example) were there for ages and most forums refer to it but some stuff is now deprecated or replaced.
Otherwise, I would probably use it for everything now given the performance and the very good FFI story.
DiffEq is a lot more than a 1000 lines of code. By now it's about 75 git repo and some of those repos are themselves >10,000 lines of code. And there is a single file in the test suites that is 10,000 lines defining mathematical constants for 200 Runge-Kutta tableaus. So it's way way more, and actually quite large by now. However, we still get a lot of feature bang for our lines of code buck.
Noob question, having played around with CUDA and NVCC a bit: in languages such as Julia or Python, why does every algorithm have to be specifically adapted to be able to run on GPUs? Couldn't algorithms be written from higher level building blocks? Eg. implement map, reduce and scanl in CUDA/OpenCL, and have a default parameter like backend='cpu' than can be set to 'gpu'.
This is precisely correct. You just have to make sure in the thousands of lines of code that you only use higher level building blocks and make no assumptions greater than what a GPU-based array allows. GPU-based arrays don't even want to allow you to do scalar indexing for example, since every scalar indexing operation is a data transfer back to the CPU. So, with those rules in mind, getting the stiff ODE solvers to automatically work with GPU-based arrays amounted to:
1. Making sure every operation in the solvers used a higher level building block (@..) that automatically built fused broadcast kernels, and CuArrays.jl then supplies the broadcast overloads. That makes all of the non-stiff methods work.
2. Change the default choice of Jacobian to "dense array of size NxN" to "zero'd outer product of input array defined via broadcast". This translates to `J = false .* z .* z'` as a trick to make a Jacobian match whatever array type the array type library thinks is best, so J will be a GPUArray while DiffEq knows nothing about GPUs (also works for things like MPI, but GPUs is easier for people to reason about these days). So a one line change to the stiff ODE solvers made them all GPU compatible! But...
3. There's a small bug, or at least missing feature, in the GPU libraries that their QR doesn't implement backsolve. Technically the user can work around that because we allow someone to pass a linear solver routine with the `linsolve` argument, but it's easier to just handle it for everyone in the defaults. So using Requires.jl, if you have CuArrays installed, then we add a linsolve routine to the defaults for you, defined by https://github.com/JuliaDiffEq/DiffEqBase.jl/blob/master/src... . Technically this can get upstreamed to CuArrays.jl and get deleted from DiffEqBase, but for now it'll live there so that everything is easier on users.
So yes, make it high level, avoid assuming things like the Jacobian live on the CPU, and add one missing method to the GPU array library, and now when the user passes `u0 = CuArray(...)` for the initial condition, the algorithm compiles a version which does all state operations on the GPU (and all logic on the CPU, so it's a pretty good mix!).
Julia is doing exactly that (even better than setting a flag, it is based on subtypes). Algorithms are usually written to work with `AbstractArray`. The usual RAM/CPU arrays (i.e. the type `Array`) is a subtype. The CUDA array which is stored on the GPU is also a subtype. Most general purpose code does not care which type you use.
Definitely look at the whole JuliaGPU ecosystem, and the CUDAnative.jl package. There are two big questions. First, how do you get code generation for targets such as GPUs from a dynamic high level language such as Julia, and that is what CUDAnative.jl achieves. The second question is what abstractions should Julia present to programmers so that increasingly larger codebases can automatically leverage GPUs. Packages such as CuArrays.jl are early answers, but there is much work to be done here.
Incredible work. Thank you everyone for the great release. This massively improves our usecase (sparse matrix equations) and I am looking forward to playing with it.
Not yet, but we are looking into it. I'm hoping to get some funding for this, since it's a no-research software dev project though, which would be great for a just-out-of-undergrad student looking to transition to software dev, hence need money to pay such a student.
Which GPUs? I'm guessing by the "CuArrays" identifier that this is CUDA-only, it'd be nice if they said NVIDIA GPUs instead of getting my hopes up. :'- (
Well the implementation is generic to the array type, so if someone creates an array type which builds broadcast kernels on with OpenCL then it would work with OpenCL. Sadly, CLArrays.jl never got updated to Julia v1.0, so until someone does that it's just a dream. However, DifferentialEquations.jl will be waiting ready for that to happen :).
I'm honestly surprised that amd doesn't get more aggressive about developing towards Julia as a deep learning target, it would be way easier for their devs to optimize for performance parity than say tensorflows spaghetti code. Iirc, the raw flops performance of AMD gpus is even or better than Nvidia but their unoptimized libraries kick it in the shins (I could be misremembering thus)
There are folks working on the AMDgpu target for Julia, and I see that some early work has already happened in terms of building LLVM and such. I suspect the AMD toolchains are not as strong - but with the right contributors, it should be possible to leverage the work done on the Nvidia GPUs for the AMD GPUs and even other accelerators.
I bought an NVIDIA GPU based laptop just to be able to program in Julia. It's a price to pay, but worth it to be able to program efficiently in my favourite language.
That's my bad. Somehow an earlier draft of the post was used. The post was updated with working code. I guess that's what happens when you lose internet and push to a Github repo from your phone :P.
Are you on the latest versions? See what running a package update does. If that doesn't work, file an issue with your package Manifest so we can track this down. Thanks.
Making differential equations 'just work' on a GPU is super cool. I feel like if this were almost any other language, it'd be impossible to get this level of code-reuse where a user can just bring their own differential equation that may not even know about GPUs and then have the differential equation library solve it on the GPU by just giving a GPUarray as the initial condition.
________________________________
Also, since the blog post doesn't seem to link to it, here's a link to the actual DifferentialEquations.jl repo: https://github.com/JuliaDiffEq/DifferentialEquations.jl