# Listable built-in functions - Why are they so fast on numerical vectors?

Consider the following small example, where I calculate Exp on a large vector by employing the listability which many mathematical built-in functions have:

data = RandomReal[1, 10^7];
RepeatedTiming[Exp[data];]
(* {0.0115, Null} *)


The list is packed but even unpacking it still has a good performance (it's most likely re-packed)

data2 = DeveloperFromPackedArray[data];
RepeatedTiming[Exp[data2];]
(* {0.0705, Null} *)


As soon as we append a non-numerical element or a number of a different type, the performance drops drastically

data3 = Append[data, 1];
RepeatedTiming[Exp[data3];]
(* {3.11, Null} *)


Now, I'm wondering if I missed something in the documentation or if there is a place, where it is explained how exactly these calls are made so performant. Specifically, I want to know

1. Does such a call use vectorization, is it parallelized or does it just call fast, but sequential code?
2. How does this compare to Compile with Parallelization turned on?
• This is just a guess, but Mathematica uses Intel's MKL, and MKL contains vectorized, parallelized, optimized functions for Exp and friends [software.intel.com/en-us/… - so my guess would be Mathematica uses MKL functions for this. I'd also guess that the MKL vector math functions beat the average C compiler - otherwise, why would C programmers use the (proprietary) MKL, for what is basically a one-liner in C. – Niki Estner Dec 3 '18 at 10:20
• @NikiEstner Yes, this is also my guess and I'm sure I already told this to some people here without knowing it for sure. Now I wonder if it is stated somewhere in the docs. One thing is sure, there are functions like LogGamma that don't perform so well which would support the MKL idea if these functions are not covered there. – halirutan Dec 3 '18 at 10:33
• Exp[data2] takes about as long as Exp[DeveloperToPackedArray@data2]. data3 cannot be repacked. What I would suggest is that Exp[data2] somehow takes advantage of the MKL without actually performing the ToPackedArray operation, since the result comes back unpacked. – Michael E2 Dec 3 '18 at 12:21
• From: support.wolfram.com/kb/39353 The kernel uses highly optimized multithreaded libraries such as Intel MKL and IPP, which are tuned for optimal performance and take advantage of advanced CPU features when available. This is important for vectorized machine arithmetic and numerical linear algebra (BLAS, …) routines, which are fundamental building blocks for many computational tasks." I wonder if further details can be found. Also: library.wolfram.com/infocenter/Conferences/6428/… – Moo Dec 3 '18 at 15:40
• FYI running a process monitor I can see that Exp[data] is calling the Intel IPP function ippsExp_64f_A53 with (on my machine) vectors of length 125000 (a total of 80 calls spread over 4 threads). – Simon Woods Dec 3 '18 at 20:55

I'm collecting what we have found out and I'm concentrating for the moment on the call Exp[data] where data is a vector of reals in machine precision. For this, the following two statements seem to apply in Mathematica 11.3 with a modern CPU:

1. Single instruction, multiple data (SIMD) instructions are used provided by the MMX technology and Streaming SIMD Extensions
2. Mathematica employs the Intel Integrated Performance Primitives library under the hood. For the Exp case, the function ippsExp_64fc_A53 is called.

Using a slightly more complex example, like

Exp[x]*Log[x + 1.0]


shows that all parts of the expression are relaid to the IPP libraries

• libippvm.so (ippsLn_64f_A53)
• libippvm.so (ippsExp_64f_A53)
• libippvm.so (ippsMul_64f_A53)

Evidently, not all functions are available in the IPP and then it gets complicated and there is no general rule. For instance

data = RandomReal[1, 10^7];
Do[LogGamma[data], {200}] // AbsoluteTiming
Do[Log[Gamma[data]], {200}] // AbsoluteTiming
(* {10.2816, Null} *)
(* {9.09973, Null} *)


It appears that the Gamma function is implemented in libWolframEngine and LogGamma is a combination of Wolfram code for Gamma and the log function from libm on my machine. The second call however uses IPP for the logarithm and appears to be a bit faster here.

As soon as we set a specific precision, different implementations are used that partly employ the The GNU Multiple Precision Arithmetic Library. As expected, the runtime drops significantly several orders of magnitude (100x - 1000x).

My second question was, how this relates to compiled code. I want to differentiate two cases

1. A parallelized compiled function that gets one number and is made listable.
2. A compiled function that receives a tensor of numbers and calls Exp[list] in its body

In the first case, it appears that no vectorization is used at all. The parallelism comes from spawning multiple threads. The second case is more interesting. My tests showed that even in the compiled code, IPP library functions are called. For simple cases, however, the runtime of the compiled code cannot compete with simply using high-level calls.

Here is one that I could find in Wolfram Support:

"The kernel uses highly optimized multithreaded libraries such as Intel MKL and IPP, which are tuned for optimal performance and take advantage of advanced CPU features when available. This is important for vectorized machine arithmetic and numerical linear algebra (BLAS, …) routines, which are fundamental building blocks for many computational tasks."

I wonder if further details can be found in the Wolfram documentation.

• You might mention that the WTC PPT is a 2006 presentation by the vendor of a defunct product. – Michael E2 Dec 3 '18 at 16:46
• @MichaelE2: Good point - in fact, that is enough to remove the link! Thanks – Moo Dec 3 '18 at 17:05