# NDSolve in Mathematica won't use all the cores avaiable

When I solve a system of differential equations in MATLAB, the task manager shows that all the CPU cores are in use. This is not true when I solve the same system in Mathematica. I have six cores. MATLAB solve the system in 3.4 seconds and Mathematica solves it in 20.4 seconds -> 6*3.4 = 20.4.

Why doesn't Mathematica use all the cores? and How can I force it to use all of them?

Mathematica code:

SetDirectory[NotebookDirectory[]];
gamma = Import["gamma.csv"];
Y[t_] := ParallelTable[Subscript[y, i, i][t], {i, 1, 816}, {j, 1, 1}];
part1 = gamma.Y[t];
RHS = Table[part1[[k, 1]] - Total[gamma[[All, k]]]*(Y[t][[k, 1]]), {k, 1, 816}, {j, 1, 1}];
ini = {IdentityMatrix[816][[815]]}\[Transpose];
sol = NDSolve[{Y'[t] == RHS, Y[0] == ini}, Flatten[Y[t]], {t, 0, 800},
Method -> {"EquationSimplification" -> "Solve"}]; // AbsoluteTiming


MATLAB code:

The function inside ode:

function dy = test(t,y, Gamma)

dy = zeros(816,1);
part1 = Gamma*y;

for iter1 = 1:816
dy(iter1,1) = part1(iter1,1) - sum(Gamma(:,iter1)).*y(iter1,1);
end


and the code running the ode:

clc
clear all
kronDel = @(j, k) j==k ;
ini = zeros(1,816);
for iter1 = 1:816
ini(1,iter1) = kronDel(815,iter1);
end

tic
[t,y] = ode23tb(@(t,y)test(t,y,Gamma),[0 800],ini);
toc


You should download gamma matrix from here.

## Edit

I think this must be related to the fact that MATLAB uses vectorization and hence can use all the CPU cores but Mathematica fails to do the same thing. In this question,related to integration, MATLAB does integration faster because it uses vectorization and as a result uses all the cores however, Mathematica won't use all the cores in that question.

• I think I have given an answer in another question of you which asks the exact same question: you can't expect that NDSolve makes use of your parallel kernels, period. If matlab has a solver which does that and solves the system you want to solve faster and to good accuracy, then you might be better off to use matlab for that task. On the other hand I'm not sure whether you are comparing apples and oranges as the two programs most probably use very different methods and probably won't give results with comparable accuracy. Mar 27, 2015 at 10:04
• I read somewhere that Matlab cannot parallelize solving differential equations. But it uses vectorization in the function it want to solve. Because I forgot the source I did not include this.
– MOON
Mar 27, 2015 at 10:37
• So you think that the vectorized operations within test are responsible for the usage of several cores and the speedup? This could of course be, but honestly I could well imagine that it probably is just ode23tb which is fast and the parallelization doesn't actually contribute that much to matlab being faster than Mathematica for this. After all ode23tb is documented to be used "if using crude error tolerances to solve stiff systems." Trading speed for accuracy is of course the easiest way to solve fast -- and if the result is good enough there is nothing wrong with that... Mar 27, 2015 at 10:58

# Update

As of v12, it's possible to parallelize when using NDSolve to solve the problem, check this post for more information. (Sadly the gamma matrix in the question is deleted so I can't test. )

As mentioned by Albert Retey in the comment above, you can't expect that NDSolve makes use of your parallel kernels, period. However, since your equation set is just a system of 1st order linear ODEs, you can turn to MatrixExp, which seems to parallelize automatically:

coe = gamma - DiagonalMatrix@Total@gamma;
init = ConstantArray[0., 816];
init[[-2]] = 1.;
solu = With[{t = 1.}, MatrixExp[coe #, init] & /@ Range[0, t, t/24]]; //AbsoluteTiming

ListPlot3D[solu[[All, ;; 800]], PlotRange -> All]


The above code is about 1 order faster compared to your approach based on NDSolve, but notice I've intentionally set t to 1., To be honest, when t becomes larger (for example, t = 800., which is chosen by you), the code becomes slower than your approach, at least on my dual-core old laptop. I guess it's because for certain parameters MatrixExp has chosen a slower method (in order to guarantee precision?), but since the document doesn't say anything about the available Method of MatrixExp (it does accept Method option! ), I'd like to stop here. Anyway your question is just about parallelism.

Try checking the parallel kernel settings by clicking on the menu bar: Evaluation > Parallel Kernel Configuration

Click the tab Parallel in the window that pops up.

Uncheck Automatic as it may have fewer kernels than you want, subject to the limit imposed by your license.

Then click Manual setting, and set the number of kernels desired.

• Did you see Mathematica use all the cores? I did what you said, it didn't work for me. I even launched all the kernels before the start of the computation.
– MOON
Mar 26, 2015 at 18:56