Are you using Polyester.jl? Large numbers of threads are not optimized with Base threads usage due to GC interactions + the hierarchical threading adds overhead vs "unsafe" thread techniques which don't support the worksharing. Polyester is thus required to get very low overhead threading matching performance of non-worksharing scenarios.
The caveat to mention here is that this is with the standard IO interface. If you use Core.println instead of println you get something substantially smaller. A limitation of v1.12 is that the handling of IO is not that great, and this is something that IIUC is being addressed in the v1.13 updates.
No it only uses the same LLVM compiler passes and you enable certain optimizations locally via macros if you want to allow reordering in a given expression.
v1.13 already has a good number of improvements to trim in the cases which are less analyzable, but even then I think the best thing for the package ecosystem is to start to become really good for trim, JETLS, etc. through statically checking things.
We use it for binary deployments via building simulations into the FMI standard (i.e. building FMUs) https://fmi-standard.org/. We still need to get libraries updated for getting the smallest possible trim when using the more computationally difficult implicit methods, but at least the stack is matured for simple methods like explicit RK and Rosenbrock-type solvers already. For folks interested in --trim on SciML, the PR to watch is https://github.com/SciML/NonlinearSolve.jl/pull/665 which is currently held up by https://github.com/SciML/SciMLBase.jl/pull/1074 which is the last remaining nugget to get the vast majority of the SciML solver set trimming well. So hopefully very soon for anyone interested in this part of the world. Note that does not include inverse problems in the binary as part of small trim.
But all of this is more about maturing the ecosystem to be more amenable to static compilation and analysis. The whole SciML stack had an initiative starting at the beginning of this summer to add JET.jl testing for type inference everywhere and enforcing this to pass as part of standard unit tests, and using AllocCheck.jl for static allocation-free testing of inner solver loops. With this we have been growing the surface of tools that have static and real-time guarantees. Not done yet, some had to be marked as `@test_broken` for having some branch that can allocate if condition number hits a numerical fallback and such, but generally it's getting locked down. Think of it as "prototype in Julia, deploy in Rust", except instead of re-writing into Rust we're just locking down the behavior with incrementally enforcing the package internals to satisfy more and more static guarantees.
Julia v1.12, the unreleased version which is currently in release candidate stage (and has had a longer RC stage than expected but should be done at least by the end of the year) has the ability to fully ahead of time compile and generate binaries, like you would expect from a language like C. It also trims these binaries so that they are a reasonable size given the elements of the runtime that are used. Thus for good type-stable code you get callable binaries without JIT overhead, and this leads to a much better CLI experience. It will take a bit of time so that there's more package tooling and support built around this feature, but this is the path a lot of the ecosystem CLI tooling is now going and that will be a pretty dramatic shift in the usability for these kinds of use cases which Julia had traditionally ignored.
This has nothing to do with floating point error accumulation or numerical stability in the floating point sense. You can do this with arbitrary sized floating point values and you will still get the same non-convergence result.
Interesting to see this here! This example is one of the ones mentioned in the appendix of https://arxiv.org/abs/2406.09699, specifically "4.1.2.4 When AD is algorithmically correct but numerically wrong".
If people want a tl;dr, the main idea is you can construct an ODE where the forward pass is trivial, i.e. the ODE solution going forwards is exact, but its derivative is "hard". An easy way to do this is to make it so you have for example `x' = x - y, y' = y - x`, with initial conditions x(0)=y(0). If you start with things being the same value, then the solution to the ODE is constant since `x' = y' = 0`. But the derivative of the ODE solution with respect to its initial condition is very non-zero: a small change away from equality and boom the solution explodes. You write out the expression for dy(t)/dy(0) and what you get is a non-trivial ODE that has to be solved.
What happens in this case though is that automatic differentiation "runs the same code as the primal case". I.e., automatic differentiation has the property that it walks through your code and differentiates the steps of your code, and so the control flow always matches the control flow of the non-differentiated code, slapping the derivative parts on each step. But here the primal case is trivial, so the ODE solver goes "this is easy, let's make dt as big as possible and step through this easily". But this means that `dt` is not error controlled in the derivative (the derivative, being a non-trivial ODE, needs to have a smaller dt and take small steps in order to get an accurate answer), so the derivative then gets error due to this large dt (or small number of steps). Via this construction you can make automatic differentiation give as much error as you want just by tweaking the parameters around.
Thus by this construction, automatic differentiation has no bound to the error it can give, and no this isn't floating point errors this is strictly building a function that does not converge to the correct derivative. It just goes to show that automatic differentiation of a function that computes a nice error controlled answer does not necessarily give a derivative that also has any sense of error control, that is a property that has to be proved and is not true in general. This of course is a bit of a contrived example to show the point that the error is unbounded, but then it points to real issues that can show up in user code (in fact, this example was found because a user opened an issue with a related model).
Then one thing that's noted in here too is that the Julia differential equation solvers hook into the AD system to explicitly "not do forward-mode AD correctly", incorporating the derivative terms into the time stepping adaptivity calculation, so that it is error controlled. The property that you get is that for these solvers you get more steps to the ODE when running it in AD mode than outside of AD mode, and that is a requirement if you want to ensure that the user's tolerance is respected. But that's explicitly "not correct" in terms of what forward-mode AD is, so calling forward-mode AD on the solver doesn't quite do forward mode AD of the solver's code "correctly" in order to give a more precise solution. That of course is a choice, you could instead choose to follow standard AD rules in such a code. The trade-off is between accuracy and complexity.
B-S is just a (bad) explicit Runge-Kutta method. You can see this formally since if you take any order of the B-S method, you can write out its computation as a Runge-Kutta tableau. If you do this, you'll see it ends up being an explict RK method that has much higher f evaluations to achieve the leading truncation error coefficients and order that the leading RK methods get. For example, 9th order B-S uses IIRC about 2-3 times as many f evaluations as Vern9, while Vern9 gets like 3 orders of magnitude less error (on average w.r.t. LTE measurements). If you're curious about details of optimizing RK methods, check out https://www.youtube.com/watch?v=s_t6dIKjUUc.
That isn't to say B-S methods aren't worth anything. The interesting thing is that they are a way to generate arbitrary order RK methods (even if no specific order is efficient, it gives a general scheme) and these RK methods have a lot of parallelism that can be exploited. We really worked through this a few years ago building parallelized implementations. The results are here: https://ieeexplore.ieee.org/abstract/document/9926357/ (free version https://arxiv.org/abs/2207.08135) and are part of the DifferentialEquations.jl library. The results are a bit mixed though, which is why they aren't used as the defaults. For non-stiff equations, even with parallelism it's hard to ever get them to match the best explicit RK methods, they just aren't efficient enough.
For stiff equations, if you are in the regime where BLAS does not multithread yet (basically <256 ODEs, so 256x256 LU factorizations but this is BLAS/CPU-dependent), then you can use parallelization in the form of the implicit extrapolation methods in order to factorize different matrices at the same time. Because you're doing more work than something like an optimized Rosenbrock method, you don't get strictly better, but using this to get a high order method that exploits parallelism where other methods end up being serial, you can get a good 2x-4x if you tune the initial order and everything well. This is really nice because 4 ODEs - 200 ODEs, where this ends up being the fastest method for stiff ODEs that we know of right now according to the SciMLBenchmarks https://docs.sciml.ai/SciMLBenchmarksOutput/stable/, and this ends up being a sweet spot where "most" user code online tends to be. However, since it does require tuning a bit to your system, setting threads properly, having your laptop plugged in, using more memory, and possibly upping the initial order manually, it's not the most user-friendly method. Therefore it doesn't make a great default since you're putting a lot more dependencies into the user's mind for a 2x-4x... but it's a great option to have for the person really trying to optimize. This also shows different directions that would be more fruitful and beat this algorithm pretty handedly though... but that research is in progress.
The whole point of the article is that "methods are considered more robust and practical if they have larger stability regions" isn't (or at least shouldn't be) necessarily true :-/.
Yeah, that's a separate post https://scicomp.stackexchange.com/questions/29149/what-does-.... I wanted to keep this post as simple as possible. If you could show cases where explicit Runge-Kutta methods outperform (by some metrics) implicit Runge-Kutta methods, then it leads to a whole understanding of "what matters is what you're trying to measure". And then yes, symplectic integrators are for long time integrations (explicit RK methods will how lower error for shorter time, so its specifically longer time integrations on symplectic manifolds, though there are implicit RK methods which are symplectic but ... tradeoff tradeoff)
There were packages I wanted to add features to last week that I've used daily for years, and when I went to go check them out, I noticed they hadn't been touched since the 1.0 release which is now over a year ago. So I think the stability of 1.0 has worked out well since, in contrast to pre-1.0, these days core packages don't even need maitanance and will still continue to work perfectly well for years.
There is no CUDA API in the standard language. The tools for automatically generating CUDA code from Julia code is the CUDAnative.jl package, which at this point is pretty stable and has fairly widespread use.