NDSolve in Mathematica won't use all the cores avaiable

Question

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
Gamma = csvread('Gamma.csv');
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.

Answer 1

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.

Answer 2

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.