# When does NDSolve parallelize ODE system solving?

I've long believed that NDSolve cannot make use of multiple cores to solve ODE system, but things seem to be different at least since v12. Consider the following toy example:

\$Version
(* "12.0.0 for Microsoft Windows (64-bit) (May 19, 2019)" *)
(* Tested on a 8-core machine. *)

eq = 2 y''[x] == y'[x] - 3 y[x] - 4;
ic = {y[0] == 5, y'[0] == 7};

With[{n = 4096}, sys = Table[{eq, ic} /. y -> y@i, {i, n}];
NDSolve[sys, y /@ Range@n, {x, 0, 10^3}]]


Seems that when n >= 4096, NDSolve automatically parallelizes. I'm pretty sure this isn't the case in v9.

So my question is:

2. In what situation does the parallelization happen? More tests show the threshold depends on the specific system being solved. In certain situations (e.g. this) parallelization seems to never happen.

3. Can we make NDSolve parallel for n < 4096? Is the parallelization controlled by any option?

• I have seen NDSolve parallelization on several occasions, and not just with enormous numbers of equations. However, I did not pay attention to the characteristics of the equations to be solved. Good question (+1). – bbgodfrey Oct 30 '19 at 22:01

The number 8192 is an autoparallelization threshold (see SystemOptions["ParallelOptions"]). There are 2*4096 = 8192 variables, y[k] and y'[k], in the system, perhaps a coincidence. The threshold can be lowered with SetSystemOptions, and one can see that when it is run, autoparallelization kicks in.

With[{ops = SystemOptions["ParallelOptions"]},
InternalWithLocalSettings[
SetSystemOptions[
"ParallelOptions" ->
"VectorParallelLengthThresholds" -> {8192, 8192, 8192, 8192,
8192}/2
],
With[{n = 4096/2}, sys = Table[{eq, ic} /. y -> y@i, {i, n}];
NDSolve[sys, y /@ Range@n, {x, 0, 10^3}]],
SetSystemOptions[ops]
]]


This is a system arithmetic feature (vectorization), not an NDSolve feature per se.

• +1, perhaps 'autocompilation' is not quite the right term – ilian Nov 4 '19 at 13:28
• @ilian Thanks for pointing it out. – Michael E2 Nov 4 '19 at 13:37
• Just a side note: in v11.3 there's no VectorParallelLengthThresholds option. There's a "MachineFunctionParallelThresholds" option whose option value is {196608, 98304, 98304, 65536, 49152}, but lowering the option value doesn't seem to trigger parallelization. – xzczd Nov 11 '19 at 7:01
• @xzczd Thank you for the bounty. Your comment explains why I couldn't find this chat message. Comparing the ParallelOptions" suggests internal parallelization has changed in V12. In V11.3 with n = 4096, I get up to 200% kernel activity (which I interpret as 2 of 4 cores being used), but I cannot find a way to control a threshold for it (neither "VectorArithmeticThresholds" nor "MachineFunctionParallelThresholds"` made a difference). – Michael E2 Nov 11 '19 at 11:56