# What would you ask Mathematica to do on a big system?

Our hardware supplier has given us the opportunity to test Mathematica on a big Windows system. The specs will be 1 to 4 TB of working memory and 12 sockets with 18 cores each. I got 1 full week of computer time to do my testing. The question is what should I ask Mathematica to do because it is more a test of Mathematica and not a test of the hardware.

I could load a complete database in memory using the association structure and then test search and consolidate questions.

Another option is to run an algorithm that benefits from parallel computations.

The goal is to write a test plan on what to test and how to test it.

I welcome your suggestions on what to evaluate in Mathematica to measure how useful many processors and almost unlimited working memory is?

Kind regards,

Pieter van Bijnen

note: it is one system (no grid involved)

• I would run MY heaviest load and compare the timings with my current equip. Then I would come here to brag. Jan 25 '16 at 19:19
• Just feeling inquisitive. Will you mind posting a little snippet of the hardware spec from your OS. It must be a hell of a windows machine. Jan 25 '16 at 19:47
• You have more than 200 cores in a shared memory system? I do wonder how the parallel tools perform ... I was never sure it could scale well to such a high number of parallel threads. Jan 25 '16 at 19:52
• My half terabyte 32 core personal workstation just lost its Mathematicarection.....
– ciao
Jan 25 '16 at 20:52
• You should run BenchmarkReport[] and make the rest of us feel insignificant :-) Jan 25 '16 at 22:26

I've always wondered about the scalability of MathLink (now officially "Wolfram Symbolic Transfer Protocol"). This is the protocol used by Mathematica to communicate between the front end and the kernel, and the basis of the Parallel package. It has quite low bandwidth and high latency relative to, for example, MPI libraries. I also wonder how many MathLink connections the library can handle before it falls apart completely. This seems an ideal opportunity to test that, and so to explore the limits of Mathematica's scaling. You will need a great number of kernel licences to be able to do this, so you will probably have to ask WRI if they will grant you a temporary licence with unlimited kernels to be able to do these tests. I think it is in their interest to know the answer, so hopefully they will.

In the Parallel package, connections are made from a master kernel to each slave kernel (one per slave). But, for message passing implementations, it would be necessary to have one connection between each pair of kernels that could potentially wish to communicate with each other. I wouldn't be surprised if you can start 216 slave kernels (although, I also wouldn't be too shocked if this is already too many). But can you open another $216^2/2 = 23328$ MathLink connections? And if so, how quickly can you transfer data from each kernel to all the others?

I wrote some code that sets up all of the connections in this answer. The notebook is here.

Assuming you can open the connections, you might like to try some all-to-all messaging. In MPI, typically benchmarks are done for message sizes increasing as powers of two. Having run the code in the other answer so that the connections are established and the required definitions are in place, you can write

alltoall[n_Integer] := ParallelEvaluate[
LinkWrite[Rank[$AllRanks], ConstantArray[$MaxMachineNumber, n]];
LinkFlush@Rank[$AllRanks]; While[! And @@ LinkReadyQ@Rank[$AllRanks], Pause[0.001]];
LinkRead@Rank[$AllRanks];, Rank[$AllRanks]
];

Table[{2^n, First@AbsoluteTiming@alltoall[2^n]}, {n, 0, 12}]


to get a set of timings.

## Bug?

Unfortunately, on my computer (64-bit Windows), the LinkWrite call hangs when n goes above 12 (message size more than 32KB but less than 64KB), no matter how many kernels/connections are involved. This is curious, because LinkWrite is supposed to be non-blocking, and MathLink is meant to be thread-safe, as far as I understand it. The deadlock is definitely in MathLink rather than my code, because even with random delays inserted anywhere in the sequence, it always completes correctly for messages smaller than the critical size. What is even stranger is that this size depends on whether the links have LinkProtocol -> "SharedMemory" or LinkProtocol -> "TCPIP". In the latter case it sometimes works for n up to 15, but also sometimes hangs using 100% CPU (a livelock).

I am not sure if this also happens on other platforms, but I have reported it to WRI under [CASE:3521773]. A self-contained notebook to reproduce it with just two kernels (not the arbitrary number proposed in this post) is available here.

## Not a bug; just poor documentation.

Despite many claims in the various sources of documentation about MathLink that LinkWrite is non-blocking, I have received a response today from WRI stating that this behavior is actually expected. The solution, as I figured out during a conversation with halirutan, is to know that LinkWrite will block until LinkRead is posted at the other end of the link, if the message is big enough. This seems not to be documented anywhere, as far as I could discern, and neither halirutan nor Szabolcs had heard of it either.

But the good thing is, when LinkRead and LinkWrite are used simultaneously, an arbitrarily large(?) message can be sent. With this knowledge it is at least possible to design a working all-to-all collective communication. I will do so and update this answer in due course, so that the 17 upvotes will not have been in vain.

• I'll try and test this on Ubuntu 14.04, 10.3.1 when I get the chance (not Windows I'm afraid) Jan 26 '16 at 9:11

I am sure you can easily install also Linux on it and then you could contact Vladyslav Shtabovenko, the current maintainer of FeynCalc (https://github.com/vsht) and ask him about hard problems in High Energy Physics he would like to benchmark on such a King-Kong machine. Either him or somebody else could also provide you with more complicated examples of FIRE (http://science.sander.su/FIRE-Examples.htm).

• Thanks Rolf, I will investigate the links! Jan 25 '16 at 21:24

I would choose a problem that exploits the unique power of Mathematica, in particular the natural functions involving graph theory, symbolic math, graphics, and the high compute power you have available. So I would choose some image recognition and clustering problem such as: Take some large number of images ($\sim\!\!10^8$) and perform deep learning on them but form clusters and find relationships that can be described by graphs. You can break the problem into parallel tasks, then merge appropriately.

• David, thanks for the suggestions. Do you have some example code I could use? Jan 25 '16 at 20:14
• Sorry, I work at a corporation and cannot give out code without lawyers getting involved. Jan 25 '16 at 20:53

4 TB memory? N-body of course;

dir = 2;

T = 10;
k = 1;
l = 10; n = l^3;
v = 6;

q[i_] := (-1)^i
m[i_] := RandomReal[5]
Rem[o_] := l (o - IntegerPart[o])

eqns = Table[
{D[Subscript[r, i][t], {t, 2}] m[i] == Sum[
Normalize[Subscript[r, i][t] - Subscript[r, j][t]] k q[
i] q[j]/(Subscript[r, i][t] -
Subscript[r, j][t]).(Subscript[r, i][t] -
Subscript[r, j][t]), {j, Delete[Table[k, {k, 1, n}], i]}
],
Subscript[r, i][0] == {Rem[IntegerPart[(i - 1)/l]/l],
Rem[(i - 1)/l], IntegerPart[(i - 1)/l^2]},
Subscript[r, i]'[0] ==
v Normalize[
Subscript[r, If[i == n - dir, n, Mod[i + dir, n]]][0] -
Subscript[r, i][0]]},
{i, 1, n}
];
sol = NDSolve[eqns, Table[Subscript[r, i][t], {i, 1, n}], {t, 0, T},
AccuracyGoal -> 5,
Method -> {"EquationSimplification" -> "Residual"}]; // Timing
Show[
Table[
ParametricPlot3D[Subscript[r, i][t] /. sol, {t, 0, T},
PlotStyle -> Hue[1 - i/(4 n)]], {i, 1, n}
],
ListPointPlot3D[
Flatten[Table[Subscript[r, i][t] /. sol /. t -> 0, {i, 1, n}], 1],
PlotStyle -> {Gray, PointSize[0.02]}],
PlotRange -> All, Background -> Black, BoxRatios -> Automatic,
Boxed -> False, Axes -> False
]
Beep[]


Try to up l, play with the definitions of forces and strain your machine. Here is the result with l=4:

What I am most interested in is described in Oleksandr's answer. Here's something else I would also try:

There are several built-in functions that take advantage of parallelization without any special settings (and without using the Parallel tools framework). I would like to know how well these scale to a high number of cores.

Examples:

Matrix operations, such as addition and matrix multiplication typically take advantage of multiple cores. I would try how Mathematica performs on tasks such as eigenvalue computations for dense and sparse matrices for very large matrix sizes. Large matrices take a lot of storage space, so you'll also make use of the memory. For large matrices it'll be realistic to compute only the few largest eigenvalues.

I noticed that NIntegrate can also use multiple cores. I would try how much speedup can be achieved on slow-to-compute integrals (e.g. multidimensional ones).

The point of this test would be to see how much speedup more cores can bring. Thus you would need to adjust the number of threads Mathematica uses for a computation and see how many times faster operations complete on 200 cores than 10 cores (i.e. how much less than 20 times).

The following system options are interesting:

In[130]:= SystemOptions["ParallelOptions"]

Out[130]= {"ParallelOptions" -> {"AbortPause" -> 2.,
"BusyWait" -> 0.01, "MachineFunctionParallelThreshold0" -> 32768,
"MachineFunctionParallelThreshold1" -> 32768,
"MachineFunctionParallelThreshold2" -> 16384,
"MachineFunctionParallelThreshold3" -> 16384,
"ParallelThreadNumber" -> 4, "RecoveryMode" -> "Retry",
"RelaunchFailedKernels" -> True}}


I believe that both "MKLThreadNumber" and "ParallelThreadNumber" are relevant. See SetSystemOptions on how to change them. I am not sure about their exact effect, but if you decide to do the test we can figure it out together.

Example benchmark:

We'll use this function to set the number of cores:

setThreadCount[
SetSystemOptions[


Eigenvalues for a dense symmetric matrix.

In[151]:= a = RandomReal[1, {5000, 5000}];
a = a + Transpose[a];


I have a laptop with a 4-core hyperthreading CPU. It appears as an 8-core one to the OS, but it only has 4 physical cores.

In[163]:= Table[
{t, First@AbsoluteTiming[Eigenvalues[a, 20];]},
{t, 1, 8}
]

Out[163]= {{1, 12.2941}, {2, 8.18082}, {3, 6.32726}, {4, 6.22842}, {5,
6.42646}, {6, 6.43579}, {7, 6.53777}, {8, 6.60184}}


Benchmarking is hard. Part of what we see here is due to the fact that as the CPU heats up it automatically throttles its clock rate. During the last few calculations it got hot so it has to slow down. If I repeat the 4-core test after waiting for it to cool down, I see better performance:

In[164]:= setThreadCount[4]

In[165]:= First@AbsoluteTiming[Eigenvalues[a, 20];]
Out[165]= 5.7838


AFAIK this won't be a problem with servers like the one you are testing.

But do pay attention to using only AbsoluteTiming and not Timing to measure parallel calculations. Timing measures the time used by each core separately then adds them up ...

• Thank you for your suggestions. Good ones. Especially to take into account the lower freq as the cpu's heat up. I'll get in contact. Feb 2 '16 at 11:56
• @pvanbijnen I don't know much about how modern CPUs work. I know that laptop CPUs slow down as they heat up, but I think (not sure) that this does not apply to server CPUs. Feb 2 '16 at 12:03

To take advantage of that kind of memory, you really want to do some parallel processing. Mathematica's parallel processing focuses on data-parallelism, or more simply, embarrassingly parallel problems. So you might try out various Monte Carlo simulations.

I don't recommend trying to reproduce full MPI functionality with LinkCreate et al. In the case mentioned, it sounds like all the links would be made to be SharedMemory. SharedMemory links handle indexing information through files in /tmp/MathLink (on Linux and OSX at least) and I suspect there might be some overhead there though this was improved greatly in 10.3.1. Ironically, you won't run into this issue with TCPIP connections though those are slower for more obvious reasons. Also, depending on the physical memory organization, the shared memory links between processes on different sockets may be significantly slower than the shared memory links between processes on the same socket, so obfuscating the architecture in an MPI-like fashion here may lead to oddities in performance.

As a last point, Mathematica's generic parallelization relies on the underlying operating system to appropriately give processes to cores on processors. On most modern laptops and desktops, this is straight forward. You have a chip with 8 logical (not including hyperthreading) cores and you launch eight subprocesses. The operating system sees those processes as active and distributes them to the 8 cores and they can run in parallel. I've never managed these high-socket count server systems before so I am not familiar with their process management. It might be more complex.

EDIT: Sorry to put this into my post as I don't have the reputation yet to leave a comment on someone else's post. Pertaining to Szalbolcs post below, if I am remembering correctly, MKLThreadNumber is the number of threads launched by the MKL library and ParallelThreadNumber is an override for when you are using subkernels.

• Could you elaborate a bit more on why Mathematica speed is determined by the OS it runs on. Would the OS Win2012R2 or a UNIX version greatly influence the parallel processing. Is process distribution different? I have no knowledge of this dependency. Feb 2 '16 at 12:05
• Your remarks on MathLink are interesting. However, I can't help thinking that doing some embarrassingly parallel problem is not a meaningful test, either of Mathematica or of a large single-image system. Fast communication among processors is the only reason one would use such a system rather than a cluster, and one pays through the nose for this capability. If I were interested in stressing the hardware, I would not do it using Mathematica. That's why I interpreted this question as asking about how to test Mathematica on this system, rather than the system itself. Feb 4 '16 at 4:56
• Sorry, pvanbijnen, I should have been more clear. Version 10.3.1 of Mathematica should improve speed significantly for systems with a lot of links, not version 10.3.1 of an OS. Also, the only reason I mention OSX and Linux is to point out that they have a /tmp/ folder that gets used by MathLink and I don't know off the top of my head where mathlink stores temp files on Windows. Feb 4 '16 at 15:34
• Oleksandr, the trouble with trying to implement something MPI-like with mathlink is that mathlink is only a communication protocol with a c-api attached. MPI does a lot more to accomplish high-performance communication, such as providing process-management, topological analysis of the hardware, cached data-communications, and many other tools. This is why MPI can scale quite well. Since mathlink does not include all of that, at best it will scale like sockets. Feb 4 '16 at 16:13
• Well, yes and no. Undoubtedly MPI libraries do take advantage of libraries like hwloc and limic2 to accomplish what one cannot easily do otherwise. But good performance is possible even without them, and process management is an MPI-2/3 addition, while many programs still only require MPI-1 functionality. The fundamental difference between MPI and MathLink, as I see it, is that MPI libraries do their utmost to minimize latency, contention, and per-communicator resource consumption. By the way, I believe that on Windows memory-mapped I/O is done through the pagefile, not separate files. Feb 7 '16 at 20:49

The more I've thought about this question, the more my answer (above) has changed. Now my answer is this: Assuming there is a high-bandwidth connection to the Wolfram server, choose a problem that relies on Mathematica's superior handling of curated data. Create an enormous problem that relies on curated financial data, geographic data, biological data, historical data, etc. Some problem that relies on search through combinatorically difficult conditions. (Other languages will never be able to deal with such a problem as well as Mathematica.)

A silly example, off the top of my head:

• Find two major cities whose mutual distance in meters is equal to their difference in population, and which traveling from the larger city to the smaller city is shortest by heading eastward, but where their average temperature difference in July 1950 is the same as in June 2015.

Silly, yes... but you get the idea.

Imagine trying to solve such a problem in Matlab or C++, or ...

Pieter, here is one suggestion: - try a FEM calculation on a huge and complicatd structure, eg. stress around clusters of cracks :)

You could just try to find the next largest prime number

http://www.iflscience.com/editors-blog/largest-ever-prime-number-found-gimps