Sure. Digging through my old emails, here's the particular example (a very stiff Van der Pol equation) I was looking at, taken from Intel's own examples, which for MATLAB would be:

  tic;
  [t, x] = ode15s(
    @(t, x) [x(2); ((1 - x(1)^2) * x(2) - x(1)) * 10^6],
    [0, 160], [2; 0], odeset('RelTol', 1e-5));
  toc;
  x(size(x, 1), :)

MATLAB takes 7.2 seconds to output the result [1.872, -0.748].

Intel's dodesol_rkm9mkn() gave me [1.877, -0.743], which is a more accurate result, in less than the blink of an eye: in 0.035 seconds.

That's more than 200x faster, by mere virtue of using a better algorithm. Apparently the algorithm is "Merson's method", which I'm completely unfamiliar with and which I'd never heard of until then.

I went ahead and ran some stuff on my laptop.

    using DifferentialEquations

    van = @ode_def_noinvhes VanDerPol2 begin
      dy = μ*((1-x^2)*y - x)
      dx = 1*y
    end μ=>1.

    prob = ODEProblem(VanDerPol2(μ=1e6),[0;2.],(0.0,160.0))
  
  
  solve(prob,CVODE_BDF(),reltol=1/10^5,save_everystep=false)
   solve(prob,Rodas5(),reltol=1/10^5,save_everystep=false)


    @time sol1 = solve(prob,CVODE_BDF(),reltol=1/10^5,save_everystep=false)
    @time sol2 = solve(prob,Rodas5(),reltol=1/10^5,save_everystep=false)
    sol3 = solve(prob,CVODE_BDF(),abstol=1/10^10,reltol=1/10^10)

On my laptop I'm getting ~0.1 seconds for both, with Sundials being a little less accurate and OrdinaryDiffEq.jl being a little more accurate against the 10^-10 reference solution than the reported Intel. But "spot checking" like this makes it hard to just timings without issues of more/less error involved, which is why more sophisticated benchmarks use work-precision diagrams. So the timings are in the same ballpark as Intel and I'd assume my (not so fast) laptop doesn't match what you benchmarked with, the point being that MATLAB is just slow here instead of Intel being unbelievably exceptional (though admirable). Still, maybe Intel does pull off some extra efficiency here, but then of course no method is universally good so I wouldn't expect that everywhere. In the end, it's probably a well optimized (probably in the DiffEq.jl/Sundials/FATODE range) single method without event handling or other extra features, but importantly it is abandoned. Probably would've been interesting if development kept up, but without it I think I'm just not going to try to get it installed.

Thanks for pointing me to it though. It's definitely an interesting nugget of ODE history that I didn't know existed.

> Probably would've been interesting if development kept up, but without it I think I'm just not going to try to get it installed.

Out of curiosity, why is this an issue? Wouldn't you have just downloaded the library and run it regardless of whether someone is actively developing it or not? I don't really understand this argument (and I see many people making it, not just you).

Well vanilla ODE solvers are a dime a dozen. Sufficiently optimized stiff ODE solvers tend to get around the same efficiency, with some methods doing better on some problems, and others doing better on other problems. So it's all problem dependent and this is just one implicit RK method, while things like ARKODE already have quite a few, FATODE has quite a few, and DifferentialEquations.jl has quite a few, etc, meaning I wouldn't expect it to be too different.

What would make it exciting to play with if it added something. Free parallelism or uses the Xeon Phi automatically, complex number support, event handling (many examples where this is necessary. One: PDE solvers need this for adding projection steps). But since development has halted, we're pretty much guaranteed that none of this will ever come.

The fact that it's not licensed for redistribution or even any regular use is a pretty good reason to not futz with it.

https://software.intel.com/en-us/articles/what-if-pre-releas...

> Intel grants to you a non-exclusive, non-assignable copyright license to make only the minimum number of copies of the Materials reasonably necessary for your internal testing and development of your products.

> You may not distribute any portion of Materials, whether in source or binary form, to any third party, except as specified in this Agreement.

I don't think this really backs up your statement though. I am not at all surprised by this result. In the post I mentioned that MATLAB wasn't really designed for efficiency, and in this repo I show how even with explicit Runge-Kutta methods (i.e. the exact same methods, and these methods are simple enough that there isn't too much variation in implementation details sans getting it all to SIMD/FMA/etc.) you can get orders of magnitude better than MATLAB:

https://github.com/JuliaDiffEq/MATLABDiffEq.jl

So since stiff solvers have a lot more moving parts and thus tend to have even larger differences, 200x is easily explained by just language differences. One way to test this is to instead test the Intel solver against Sundials. I wouldn't be surprised if Sundials actually was faster in this case.

But there's another factor here. I explained in the blog post that multistep methods like BDF tend to do well on asymtopically large problems since they minimize the amount of user-function evaluations, but they don't tend to do well on smaller systems of ODEs. This is exactly that case: you're comparing an implicit RK (Intel) to MATLAB's BDF method, and see that the RK is more accurate and is fast. But this effect of the method choice on small (user-scale) problems is exactly why I emphasized it as a downside when solver suites only include multistep BDF methods as their stiff solvers.

That isn't to say there isn't interesting here. Intel's nonstiff method is a 4th order method that can add extra 1st order stages to boost the stability region, similar to a Runge-Kutta Chevyshev method. I honestly have no idea what this stiff solver is so I bookmarked it to read when I get to a campus where the paper is free:

http://www.sciencedirect.com/science/article/pii/00415553899...

But in general, RKC methods weren't widely adopted because usually using a full non-stiff or stiff solver ends of being better, and picking wildly different methods like this which are not widely known usually isn't a good sign in a very "crowded" field like ODE solvers where generally things are pretty figured out and it's all about the details.

> I don't think this really backs up your statement though.

I'm not sure which statement you're reading, but everything you said so far is in agreement with my statement:

>> Intel has an ODE solver that seems unbelievably good at solving certain kinds of differential equations. I'm surprised/disappointed it wasn't even mentioned. I would've loved to see an evaluation of it against MATLAB (or other solvers) more than any other pair in the list.

Note that I have specifically stated that this method seems better on (at least) certain kinds of equations, and that, therefore, I would have loved to see a comparison of it with everything else.

In particular, note that I never claimed the method Intel uses is better, at least (especially?) not in general. If I knew that to be the case, I wouldn't have even had a reason to hope for a comparison! The entire point of my comment was that I would be interested in an evaluation of its performance, and I'm glad to see that you feel similarly. :)

Gotta say, "long-abandoned implementation of some unclear method that was considerably faster than Matlab's known-to-be-slow implementation on some problems" does not sound particularly interesting. I suspect that Chris is only interested at all because he is spectacularly thorough (as evinced his comparisons with other ODE solvers and his code) and the method might not yet be included in DifferentialEquations.jl. There's little reason to believe that Intel was onto something "magical" here; it seems that they just had a good implementation of a decent method.

Finally got on campus to be able to read the paper in detail. It's an interesting method but it doesn't really seem to have a clear "best use case". It is a Rosenbrock-W (or ROW) type of method, which means that it's a Rosenbrock method where the order of accuracy is not dependent on having an "exact Jacobian". It then uses this fact to allow for "freezing", i.e. storing the same Jacobian across multiple steps. The reason this is nice is because then you can store the factorization and re-use it. (m,k) is with m being the number of linear solves, and k being the number of RHS function evaluations. The paper creates a bunch, and what Intel used was a (5,2) method.

However, the (5,2) method, in order to have an embedded method for adaptivity, has its coefficients (for the Rosenbrock tableau) almost completely determined. What this means is that the paper is not able to "optimize" the method for low error. Normally what's done is that a method is written in terms of tableaus, and the numbers for the tableau are chosen such that the principle terms in the error, the coefficients of the leading error terms, are minimized. Sufficiently good RK and Rosenbrock methods can get like 10^-8, meaning that for large enough dt they almost act like they are an order higher. Of course, the large amount of assumptions to get the given method means that this couldn't be done.

But does the algorithm have a place? This is the issue... it doesn't necessarily have a good fit anywhere. Standard Rosenbrock methods require a new Jacobian each timestep, but that's fine if factorization isn't sufficiently expensive when is true for smaller ODE systems. Also if the problem is "very nonlinear", i.e. not semi-stiff but solidly stiff, Jacobian re-usage doesn't help anyways. So standard coefficient-optimized Rosenbrock methods do really well when the equations are very stiff and the number of ODEs is small enough. (coefficient-optimized) Modern (E)SDIRK methods naturally allow for Jacobian re-usage but are otherwise similar to the Rosenbrock methods, making them do well in the small number of ODEs but semi-stiff case. The issue with this method is that it tries to sit in the middle here but it's not really optimized for either case.

But the issue really goes back to the 2 in the (5,2). Early on in ODE literature, people believed that reducing the number of steps was always a good idea. You can do 5th order RK methods in 7 steps (thus 7 f evaluations, where u'=f with f being the user given ODE function), so why would you ever do 8? However, in the early 90's when we finally got our hands on PCs, people found out that the assumption to keep stages minimized doesn't empirically do well. Instead, what you can do is have more stages, but then use that to design the method such that the local truncation error is orders of magnitude lower. The result is that you can take much larger steps in order to solve the equation but get the same error, making everything quicker. To compensate for the fact that you are now taking large steps you have to have "dense output", i.e. a continuous extension of the method which can use the internal stage data to turn the solution into a continuous function sol(t) which doesn't require any evaluations of f. While this is more cumbersome in to implement, these "higher stage with dense output" methods, the engineering is very much worth it and if you take a look at pretty much any production-quality method these days you'll see they've ended up with these Rosenbrock, (E)SDIRK, or FIRK methods... except for multistep methods.

Multistep methods though use history to achieve high order with just a single function evaluation in its implicit function per step. This is thus what's best when f is expensive, and is another reason why stage reduction isn't necessarily helpful because then you're ending up in an area where other methods are already pretty optimal. The classic is the BDF method. It can do the Jacobian re-usage as well. While these aren't L-stable above order 2, these can be adaptive order and do really well in mildly stiff situations.

So that said, there really isn't a good niche for these types of methods. I have other low-stage Rosenbrock tableaus implemented and it comes out pretty much as I say here: more error for the same amount of time as others. Rosenbrock-W (ROW) methods are not important anymore because with autodifferentiation we don't need to do numerical differentiation and thus we do get exact Jacobians easily, which negates the purpose of the ROW method (the reason why they were invented because of numerical errors with numerical differentiation causing order loss in Rosenbrock methods. The tradeoff is a lot more assumptions and a lot less coefficient optimization).

So the method Intel chose, extrapolating from other low-stage Rosenbrock methods and its other properties, would do best in a cases where the equation is not too nonlinear (you want to make use of the Jacobian freezing, otherwise there's no point since the standard Rosenbrock methods would simply get less error for the same steps), but also it's very stiff (otherwise you're in the regime where BDF dominates), but also not too stiff where you'd want to use the IMEX properties of the modern (E)SDIRK methods, and you'd also want f sufficiently costly that the higher-stage (E)SDIRK methods don't do better.

This puts it into a very very small niche, which is probably why this method was lost to the literature and not really adopted. But then again, since it does sit right in the middle of a lot of things, if you're to pick one method to implement, this is probably a pretty good choice. So while Intel probably did a very good implementation of it, I wouldn't see a reason to go looking for it somewhere else. If this method did well, then either a Rosenbrock, (E)SDIRK, or BDF method will probably do better depending on what regime it more falls into.

Wow! Thank you for the comprehensive review. Leaves me with nothing else to say but thanks :)

Do the initials stand for Runge-Kutta-Merson, order 9? I wonder about the other 3 (mkn)

I believe those are the initials, yes.

The "9" refers to the maximum number of "stages" which is not the same thing as the order of the solver.

The (m,k) refer to the corresponding parameters of the method (I think).

The "n" at the end is for "numerical" Jacobian (as opposed to "a" for user-provided/"analytical").

There should be a PDF in the download that explains these better than I possibly can (though it's still somewhat terse).

I had a look before asking but could not find anything before the download. Thanks, it piqued my curiosity (used to write & theoretically analyse several numerical integrators back in my teaching days)